这个时候已经不是表达矩阵的事情了，要从新从fastq测序数据开始。
对测序后的fastq数据进行转录本的组装。基于组装后的转录本，通过数据库注释去掉编码蛋白质的mRNA以及数据库中收集的已知的lncRNA，对剩余的转录本进行生物信息学分析，最终鉴定出全新的lncRNA，作为后续研究的起点。
第一步 ：重构转录本 --stringtie
STEP4: 得到表达矩阵的流程用比对软件hisat2将reads比对到参考基因组得到bam文件，如果要鉴定新的转录本，需要重新组装转录本，可以用的软件有cufflinks，stringtie，这里用stringtie。
REF=/pnas/fangxd_group/renyx/macaque/00ref
assemble_out=/pnas/fangxd_group/renyx/macaque/07assemble_out
align_out=/pnas/fangxd_group/renyx/macaque/03align_out/hisat2_mapping
stringtie -p 4 -G $REF/Macaca_mulatta.Mmul_8.0.1.91.gtf -o $assemble_out/OC_1yrF.stringtie.gtf -l $align_out/OC_1yrF SRR4042230_sorted.bam
stringtie -p 4 -G $REF/Macaca_mulatta.Mmul_8.0.1.91.gtf -o $assemble_out/OC_1yrM.stringtie.gtf -l $align_out/OC_1yrM SRR4042231_sorted.bam
第二步：预测新的转录本 --Cuffcompare
cuffcompare是cufflinks其中的一个软件，Cuffcompare提供了一种有效的分类和注释方法，即将重建转录组与现有基因注释进行比较，以获取重建转录组的分类，并用类别代码(class code)加以标示。
cuffcompare 用法及参数说明：
cuffcompare [-r <reference_mrna.gtf>] [-R] [-T] [-V] [-s <seq_path>] [-o <outprefix>] [-p <cprefix>] {-i <input_gtf_list> | <input1.gtf> [<input2.gtf> .. <inputN.gtf>]}
-r 参考基因组的注释文件
-o 输出文件的前缀
-i 输入文件
cuffcompare -r $REF/Macaca_mulatta.Mmul_8.0.1.91.gtf -o $assemble_out/cufcompF $assemble_out/OC_1yrF.stringtie.gtf
cuffcompare -r $REF/Macaca_mulatta.Mmul_8.0.1.91.gtf -o $assemble_out/cufcompM $assemble_out/OC_1yrM.stringtie.gtf
输出文件包括6个：
-
cufcompF.combined.gtf
*.combind.gtf结果包含很多信息，如exon的位置信息，gene_id和transcript_id(stringtie内部给的ID), gene_name(ensemble_ID), class_code等。
-
cufcompF.loci
此文件中包含了stringtie给的gene_ID (XLOC_000001)与ensemble gene_id 和transcript_id，及exon序列位置的信息等。
-
cufcompF.OC_1yrF.stringtie.gtf.refmap
这个文件包含四列信息，第一列ref_gene_id是gene symbol ，无symbol的给出的是ensemble的gene id; 第二列ref_id是指ensemble的transcript id; 第三列class_code 是“=”和“c”;第四列是cuff_id_list。这个文件指组装后与参考基因组几乎完全匹配的转录本。
-
cufcompF.OC_1yrF.stringtie.gtf.tmap
这个文件很重要，包括很多有用的信息，如FPKM,coverage,length,gene id ,class code等，可用于lncRNA初步筛选。
-
cufcompF.stats
一些基本的统计信息，可以看到novel exons和novel introns 的比例。
-
cufcompF.tracking
这个文件的信息都包含在*gtf.tmap.
第三步：筛选coverage,length,FPKM
对于单个转录本的组装结果，按一下要求筛选转录本：
1）FPKM>=0.5
2）coverage >1
3）Length > 200
awk '{if($7>=0.5 && $10 > 1 && $11 >200) print $0}' cufcomp.OC_1yrF.stringtie.gtf.tmap > filter.OC_1yrF
awk '{if($7>=0.5 && $10 > 1 && $11 >200) print $0}' cufcomp.OC_1yrM.stringtie.gtf.tmap > filter.OC_1yrM
第四步 ：class code分类
class_code分类的具体含义: “=”代码表示此预测转录本与注释基因的所有内含子完全吻合，但它们在第一外显子(first exon)的起始端或最后外显子(last exon)的末端可能有差别。然而，这并不影响将“=”类重建转录本判定为已注释转录本。又如，转录本标有“j”类别代码，表明此转录本至少有一个内含子与已注释基因的内含子相同，而其他位置可能不同，据此可推断此类转录本可能是注释基因的一个新异构体(novel isoform)。另外“i,o,u,x”的分类符合lncRNA的特征，可用于lncRNA的识别过程。因此，“i,j,o,u,x”这5类转录本表示可能是新的转录本，符合lncRNA的要求，保留作为后续分析。
1 = Complete match of intron chain
2 c Contained
3 j Potentially novel isoform (fragment): at least one splice junction is shared with a reference transcript
4 e Single exon transfrag overlapping a reference exon and at least 10 bp of a reference intron, indicating a possible pre-mRNA fragment.
5 i A transfrag falling entirely within a reference intron
6 o Generic exonic overlap with a reference transcript
7 p Possible polymerase run-on fragment (within 2Kbases of a reference transcript)
8 r Repeat. Currently determined by looking at the soft-masked reference sequence and applied to transcripts where at least 50% of the bases are lower case
9 u Unknown, intergenic transcript
10 x Exonic overlap with reference on the opposite strand
11 s An intron of the transfrag overlaps a reference intron on the opposite strand (likely due to read mapping errors)
12 . (.tracking file only, indicates multiple classifications)
awk '{if ($3=="u" || $3=="x" || $3=="i" || $3=="j" || $3=="o"){print $0}}' filter.OC_1yrF > class.OC_1yrF
awk '{if ($3=="u" || $3=="x" || $3=="i" || $3=="j" || $3=="o"){print $0}}' filter.OC_1yrM > class.OC_1yrM
参考资料：
基于RNA-Seq的lncRNA预测流程介绍
cuffcompare介绍
转录组的组装STINGTIE和CUFFLINKS
网友评论