Step1：链特异性比对并组装转录本
用HISAT2对每个样本进行链特异性比对，随后利用StringTie对每个样本的转录本进行单独组装
mkdir cleandata hisat2 stringtie count target
#利用HISAT2进行链特异性比对，--rna-strandness RF
#以sample1为例
hisat2 --dta -x ~/genome/hisat2_index/genome.hisat2 -p 4 --rna-strandness RF -1 ./cleandata/lnc.rna.sample1_R1.fq.gz -2 ./cleandata/lnc.rna.sample1_R2.fq.gz -S ./hisat2/lnc.rna.sample1.sam
samtools view -@ 4 -bS -q 30 ./hisat2/lnc.rna.sample1.sam | samtools sort -@ 4 - -o ./hisat2/lnc.rna.sample1.bam
rm ./hisat2/lnc.rna.sample1.sam
samtools flagstat -@ 4 ./hisat2/lnc.rna.sample1.bam > ./hisat2/lnc.rna.sample1.stat.txt
#利用StringTie对每个样本单独组装转录本
stringtie -p 4 --rf -G ~/genome/genome.gff3 -o ./stringtie/lnc.rna.sample1.gtf -l lnc.rna.sample1 ./hisat2/lnc.rna.sample1.bam
Step2：过滤lncRNA
将所有样本组装的转录本合并，随后根据lncRNA的特征进行过滤
#利用StringTie将所有样本的转录本合并
##其中gtf.list.txt包含所有要合并的gtf文件，一行一个
ls ./stringtie/*gtf > gtf.list.txt
stringtie --merge -p 4 --rf -G ~/genome/genome.gff3 -o ./stringtie/all.merged.gtf -l lncrna gtf.list.txt
#利用gffcompare将组装好的转录本与原始基因注释文件进行比较
gffcompare -r ~/genome/genome.gff3 -o ./stringtie/gffcmp ./stringtie/all.merged.gtf
#根据gffcompare的结果过滤转录本，根据lncRNA的特性，只保留calss code为i,u,x,o并且长度大于200bp的转录本
#i: 完全包含在参考的内含子内
#u: 位于基因间区
#x: 与参考外显子存在overlap但与参考转录本方向相反
#o: 与参考外显子存在overlap且与参考转录本方向一致
awk '$3=="i" || $3=="u" || $3=="x" || $3=="o" && $10>200 {print "\"" $5 "\""}' ./stringtie/gffcmp.all.merged.gtf.tmap > ./stringtie/filter.iuxo.l200.transcript.id.txt
LC_ALL=C fgrep -f ./stringtie/filter.iuxo.l200.transcript.id.txt ./stringtie/gffcmp.annotated.gtf > ./stringtie/filter.iuxo.l200.gtf
#根据gtf文件提取转录本的序列
gffread ./stringtie/filter.iuxo.l200.gtf -g ~/genome/genome.fa -w ./stringtie/filter.iuxo.l200.transcript.fa
#预测转录本的编码能力
##CPC2，依赖python3
python ~/software/CPC2/CPC2_standalone-1.0.1/bin/CPC2.py -i ./stringtie/filter.iuxo.l200.transcript.fa -o ./stringtie/CPC2.out --ORF
#CNCI，依赖python2.7
#-m pl指定模型为植物
python2.7 ~/software/CNCI/CNCI-master/CNCI.py -f ./stringtie/filter.iuxo.l200.transcript.fa -o ./stringtie/CNCI.out -p 4 -m pl
##LGC，依赖python2.7
python2.7 ~/software/LGC/LGC-master/scr/lgc-1.0.0.py ./stringtie/filter.iuxo.l200.transcript.fa ./stringtie/LGC.out.txt
#过滤掉具有蛋白编码能力的转录本
awk '$9=="noncoding"{print "\"" $1 "\""}' ./stringtie/CPC2.out.txt > ./stringtie/tmp.CPC2.txt
awk '$5=="Non-coding"{print "\"" $1 "\""}' ./stringtie/LGC.out.txt > ./stringtie/tmp.LGC.txt
awk '$3=="noncoding"{print "\"" $1 "\""}' ./stringtie/CNCI.out/CNCI.index > ./stringtie/tmp.CNCI.txt
cat ./stringtie/tmp.CPC2.txt ./stringtie/tmp.LGC.txt | sort | uniq -d > ./stringtie/tmp1
cat ./stringtie/tmp.CPC2.txt ./stringtie/tmp.CNCI.txt | sort | uniq -d > ./stringtie/tmp2
cat ./stringtie/tmp1 ./stringtie/tmp2 | sort | uniq -d > ./stringtie/filter.iuxo.l200.ncd.transcript.txt
LC_ALL=C fgrep -f ./stringtie/filter.iuxo.l200.ncd.transcript.txt ./stringtie/filter.iuxo.l200.gtf > ./stringtie/filter.iuxo.l200.ncd.gtf
gffread ./stringtie/filter.iuxo.l200.ncd.gtf -g ~/genome/genome.fa -w ./stringtie/filter.iuxo.l200.ncd.transcript.fa
#利用pfam数据库进行过滤
##用Transeq将转录本翻译为6种可能的蛋白序列
transeq ./stringtie/filter.iuxo.l200.ncd.transcript.fa ./stringtie/filter.iuxo.l200.ncd.pep6.fa -frame=6
##设置E-value < 1e-5为阈值
pfam_scan.pl -cpu 4 -fasta ./stringtie/filter.iuxo.l200.ncd.pep6.fa -dir ~/database/pfam-A/ -o ./stringtie/pfam.out
##根据pfam结果过滤存在蛋白结构域的transcripts
sed "/#/d" ./stringtie/pfam.out | sed '/^$/d' | awk '$13 < 1e-5 {print $1}' | sed "s/_/\t/g" | awk '{print "\"" $1 "\""}' > ./stringtie/pfam.coding.txt
LC_ALL=C fgrep -v -f ./stringtie/pfam.coding.txt ./stringtie/filter.iuxo.l200.ncd.gtf > ./stringtie/final.lncRNA.gtf
Step3：转录本定量
将lncRNA与基因组注释整合，随后用Stringtie进行转录本水平FPKM定量
#进行转录本水平定量(以sample1为例)
cat ./stringtie/final.lncRNA.gtf ~/genome/genome.gtf > count/final.all.gtf
##转录本的FPKM值在sample1目录下的t_data.ctab文件中
stringtie -p 4 -b ./count/sample1 --rf -G count/final.all.gtf -e -o ./count/tmp.sample1.gtf -l sample1 ./hisat2/lnc.rna.sample1.bam
从中提取count值用于后续差异分析
#提取count值
##其中lncrna.sample.txt包含了所有样本的信息，两列，第一列为样本名，第二列为样本的gtf名
prepDE.py -i ./count/lncrna.sample.txt -g ./count/all.gene.count.csv -t ./count/all.transcript.count.csv
Step4：lncRNA靶基因预测
分别预测lncRNA的cis和trans的靶基因
#lncRNA靶基因预测（只针对FPKM>1的lncRNA和mRNA）
##cis-lncRNA靶基因，一般根据lncRNA和mRNA的位置信息进行判定，比如上下游100k
awk '{print $1 "\t" $2-100000 "\t" $3+100000 "\t" $4}' ~/genome/all.gene.bed | awk '$2<0 {print $1 "\t" 0 "\t" $3 "\t" $4} $2>=0 {print $0}' > ./target/all.gene.+-100k.bed
awk '$3=="transcript" {print $1 "\t" $4 "\t" $5 "\t" $10 }' ./stringtie/final.lncRNA.gtf | sed "s/\"//g" | sed "s/\;//g" > ./stringtie/final.lncRNA.bed
bedtools intersect -wo -a ./stringtie/final.lncRNA.bed -b ./target/all.gene.+-100k.bed | cut -f 4,8 > ./target/lncrna.cis.target.txt
##trans-lncRNA靶基因
##利用LucTar预测lncRNA经过trans作用的靶基因
#-d设置ndG（the normalized deltaG）阈值，小于该阈值则认为存在互作，大于该阈值则认为不存在互作。建议设为：-0.08(lowest);-0.10(low);-0.13(medium);-0.15(high);-0.20(highest)
#-s F/T，设置是否输出配对的序列（批量运算时不建议输出）
##保证fa文件的基因ID行只有基因ID，没有空格后面的，否则会报错“the format of file is error, Please check! The name and the sequence do not put in a row!”
gffread final.lncRNA.gtf -g ~/genome/genome.fa -w final.lncRNA.transcript.fa
perl LncTar.pl -p 1 -l final.lncRNA.transcript.fa -m mRNA.transcript.fa -d ndg -s F -o lncRNA.trans.mRNA.LucTar.txt
