基于RNA-seq测序数据的植物lncRNA鉴定

Oligo-dT 法的测序数据无法捕获无poly尾的RNA，所获得的lncRNA不能完全代替全转测序

但是用来糊弄组会还是可以的

请引用

Plant-LncPipe: a computational pipeline providing significant improvement in plant lncRNA identification

还要使用工具 salmon、等，全程wsl ubuntu2204，Python 3.8.11，开始吧

1. Plant-LncRNA-pipline的教程讲的非常详细，首先安装

https://github.com/xuechantian/Plant-LncRNA-pipline

hisat2安装时候可能报错，存在一些依赖关系，建议照着官方走一遍。
StringTie 照着来就行

LncFinder-plant安装也有点麻烦，ubuntu可能需要一些依赖比如：

sudo apt-get install libfreetype6-dev libpng-dev libtiff5-dev libjpeg-dev
sudo apt-get install libfribidi-dev libharfbuzz-dev libpng-dev libjpeg-de

CPAT稍微麻烦点，可以直接使用pip安装，而不用创建python=2.7环境，后面会改脚本。我自动下载的是cpat 3.0.5
Diamond 用mambar或者conda安装
FEELnc安装也有点小坑，建议参照官方安装指南，主要是其中的linux 和MAC的区别

可能还需要安装的东西“

sudo cpan Parallel::ForkManager
sudo cpan Bio::DB::SeqFeature::Store

2.开始鉴定lncRNA

首先下载数据，仍然强烈建议使用iseq。下好了以后转换成fastq（这里采用fastaq.gz）

将下载时候的ssr号丢进sample.list，注意LF分隔符，这里假设你的文件都是正反向测序（strand-specific）

1.把你的基因组命名为genome.fasta，建立index

hisat2-build -p 24 genome.fasta genome.index

-p 24 表示线程是24，对于400mb的基因组执行了大约3min后

2.开始比对，


for i in `cat sample.list`; do hisat2 --new-summary --rna-strandness RF -p 24 -x genome.index -1 ${i}_1.fastq.gz -2 ${i}_2.fastq.gz -S ${i}.sam; done
for i in `cat sample.list`; do samtools sort -o ${i}.bam ${i}.sam; done
######或者##########
#samtools可能会和环境下一些包冲突，比如找不到libssl1.0.0，没事切个新的conda环境运行
#-l 9 表示最大压缩，但是速度会减慢
cat sample.list | parallel -j 16 'samtools sort -l 9 -m 4G --threads 24 -o {1}.bam {1}.sam'

合并上面的步骤


# 先确保 fastaq2bam.log 存在
touch fastaq2bam.log

# 使用 parallel 执行任务
cat sample.list | parallel -j 4 --joblog job.log --progress '
    # 检查 BAM 文件是否存在
    if [ ! -f {1}.bam ]; then
        # 如果 BAM 文件不存在，检查 SAM 文件是否存在
        if [ ! -f {1}.sam ]; then
            # 如果 SAM 文件也不存在，运行 hisat2
            hisat2 --new-summary --rna-strandness RF -p 24 -x genome.index -1 {1}_1.fastq.gz -2 {1}_2.fastq.gz -S {1}.sam
        fi
        # 排序 SAM 文件生成 BAM 文件
        samtools sort -l 9 -m 3G --threads 48 -o {1}.bam {1}.sam
        # 删除中间生成的 SAM 文件
        rm {1}.sam
    fi
    # 检查并记录 BAM 文件的完整性
    samtools quickcheck {1}.bam 2>&1 | tee -a fastaq2bam.log | grep -v "^$"
'

对于4.5g每个的双端测序数据，生成了大概45g的sam文件，压缩9后又回到大概是原来的大小
amd 7900x 对于1.5g的fq gz文件大概运行了20min

3.转录本重建

将基因组注释重命名为genome.gtf

parallel -j 4 'stringtie -p 24 --rf -G genome.gtf -o {}.gtf {}.bam' :::: sample.list
#这里假设文件格式是SRR*.gtf
stringtie --merge -p 24 -o candidate_transcript.gtf -G genome.gtf SRR*.gtf
gffread -w candidate_transcript.fasta -g genome.fasta candidate_transcript.gtf
grep '>' candidate_transcript.fasta | awk '{print $1}' | sed 's/>//g' | sort -u > candidate_transcript.txt

这里有个小坑是你的gff文件可能是有问题的，建议使用agat_convert_sp_gff2gtf.pl 修复。可能提示找不到transcript_id，可以通过sed命令替换下手动批量添加信息。

4.鉴定lncRNA，照着官方教程来

FEELnc_filter.pl -i candidate_transcript.gtf -a genome.gtf --monoex=-1 -s 200 -p 48 > candidate_lncRNA.gtf
cut -d ";" -f 2 candidate_lncRNA.gtf |sed 's/ transcript_id //g' | sed 's/"//g' | sort -u > candidate_lncRNA.txt

进R操作

这里把Plant-LncRNA-pipline中的Model文件夹和example文件夹拷贝到当前文件夹下

#R
library(LncFinder)
library(seqinr)
mRNA <- seqinr::read.fasta(file ="example_data/training_mRNA.fasta")
lncRNA <- seqinr::read.fasta(file ="example_data/training_lncRNA.fasta")

frequencies <- make_frequencies(cds.seq = mRNA, lncRNA.seq = lncRNA, SS.features = FALSE, cds.format = "DNA", lnc.format = "DNA", check.cds = TRUE, ignore.illegal = TRUE)  

plant = readRDS("./Model/Plant_model.rda")
Seqs <- seqinr::read.fasta(file ="candidate_transcript.fasta")
Plant_results <- LncFinder::lnc_finder(Seqs, SS.features = FALSE, format = "DNA", frequencies.file = frequencies, svm.model = plant, parallel.cores = 24)
write.table(Plant_results, file ="plant-lncFinder.txt", sep ="\t",row.names =TRUE, col.names =TRUE,quote =FALSE)

退出R

5.这里要注意下，cpat3.0以上，不用教程的内容

cpat -x ./Model/Plant_Hexamer.tsv -d ./Model/Plant.logit.RData -g candidate_transcript.fasta -o CPAT_plant.outputcpat -x ./Model/Plant_Hexamer.tsv -d ./Model/Plant.logit.RData -g candidate_transcript.fasta -o CPAT_plant.output


#下载慢的话挂代理或者windows下载了拖进来
wget https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gunzip uniprot_sprot.fasta.gz
diamond makedb --in uniprot_sprot.fasta -d uniprot_out
diamond blastx -d uniprot_out -q candidate_transcript.fasta -o uniprotoutput.txt

这里在当前文件夹创建一个脚本 insersection.change.sh（sudo nano insersection.change.sh）

#!Rscript
# Usage: Rscript myscript.R
##################################################################################################################
#Taking the intersection of lncRNA identification results, to obtain high confidence lncRNAs.
##################################################################################################################
    
args<-commandArgs(TRUE)

#################################################
## R package
myPaths <- .libPaths()
library(tidyverse)


#################################################
##### 
candidate_lncRNA <- read.table(args[1], header = FALSE, col.names = "genes") %>% pull(genes)


##### shorter than 200 bp and overlapping with known mRNAs.
length <- read_delim(args[2], delim = "\t", col_names = FALSE) %>% pull(X1)


##### CPAT-plant###change
CPAT <- read.table(args[3], header = FALSE, col.names = "mRNA") %>% pull(mRNA)




##### LncFinder-plant
LncFinder <- read_delim(args[4], delim = "\t") %>% filter(Coding.Potential == "NonCoding") %>% pull(Pred)


##### Swissprot
uniprot <- read_delim(args[5], delim = "\t", col_names = FALSE) %>% filter(X3>60, X11<1e-5) %>% pull(X1) %>% unique()
uniprot <- dplyr::setdiff(candidate_lncRNA, uniprot)


##### Take intersection
lncRNA = Reduce(intersect, list(length, CPAT, LncFinder, uniprot))
write.table(lncRNA, file = "final_lncRNA_results.txt",
            row.names = F,
            col.names = F,
            quote = F)

合并结果（之前的脚本是判断CPAT输出文件中<0.46，新版本改为了直接输出no_ORF文件，所以这里要改）


Rscript insersection.change.sh candidate_lncRNA.txt candidate_lncRNA.txt CPAT_plant.output.no_ORF.txt plant-lncFinder.txt uniprotoutput.txt

grep -Fwf final_lncRNA_results.txt candidate_transcript.gtf > lncRNA.gtf

FEELnc_classifier.pl -i lncRNA.gtf -a genome.gtf > lncRNA_classes.txt

6.其他的都一样了，输出结果到output文件夹

mkdir -p ./output/

FEELnc_classifier.pl -i lncRNA.gtf -a genome.gtf > ./output/lncRNA_classes.txt

awk -F '\t' '{if($1==1 && $6 == "antisense" && $10 == "exonic") {print $0}}' lncRNA_classes.txt > ./output/LncRNA_antisense_exonic.txt

awk -F '\t' '{if($1==1 && $6 == "sense" && $10 == "intronic") {print $0}}' lncRNA_classes.txt > ./output/LncRNA_intronic.txt

awk -F '\t' '{if($1==1 && $6 == "sense" && $7 == "intergenic" && $8 <= 2000 && $10 == "upstream") {print $0}}' lncRNA_classes.txt > ./output/LncRNA_upstream.txt

awk -F '\t' '{if($1==1 && $6 == "sense" && $7 == "intergenic" && $8 <= 2000 && $10 == "downstream") {print $0}}' lncRNA_classes.txt > ./output/LncRNA_downstream.txt

awk -F '\t' '{if($1==1 && $7 == "intergenic" && $8 > 2000) {print $0}}' lncRNA_classes.txt > ./output/LncRNA_intergenic.txt

awk -F '\t' '{if($1==1 && $6 == "antisense" && $7 == "intergenic" && $8 <= 2000 && $10 == "upstream") {print $0}}' lncRNA_classes.txt > ./output/LncRNA_Bidirectional.txt

cp lncRNA.gtf ./output/

3. 定量

使用salmon定量的文献：

Identification and comparative analysis of long non-coding RNAs in the brain of fire ant queens in two different reproductive states

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容