转录组的组装Stingtie和Cufflinks

首先这两款软件都是用于基于参考基因组的转录组组装，当然也可用于转录本的定量。前者于2016年的 protocol上发表的转录组流程HISAT, StringTie and Ballgown后被广泛使用，后者则是老牌的RNA分析软件了。在算法上来说Stringtie使用的是流神经网络算法，Cufflinks则是吝啬算法；从组装效果上来看Stringtie在灵敏度和准确度上表现较好，能够拼接出更完整、更准确的基因；从定量上来说，两者相差不大，但是cufflinks在一些特殊情况下会有异常的表达量；从运行速度上来说，Stringtie远远快了cufflinks了。。。

以上均为听说。。。

不仅以上的所述，最近的一篇Gaining comprehensive biological insight into the
transcriptome by performing a broad-spectrum
RNA-seq analysis转录组各个分析流程对比的文章中也提到HISAT2+StringTie的搭配使用效果是在基于参考基因组转录本定量中最好的一个。所以我果断选择HISAT2+Stingtie来代替TopHat+Cufflinks来用于转录组的组装

这里主要是记下StringTie的使用，HISAT2则先略过了；当然还有一个Cuffcompare不得不提

下载及安装

直接下载二进制软件

wget http://ccb.jhu.edu/software/stringtie/dl/stringtie-1.3.3b.Linux_x86_64.tar.gz
tar zxvf stringtie-1.3.3b.Linux_x86_64.tar.gz

Stringtie的使用

1. 将bam文件通过stringTie进行组装，以一个公共数据小鼠为例

   stringtie -p 20 -G ~/reference/genome/mm10/gencode.vM13.annotation.gtf -o control1.gtf control1_sorted.bam

   这里的参数比较简单：

   -p 线程数  
   -G 基因组注释文件  
   -o 输出的gtf文件

   还有一些其他参数，比如：

   -m 组装预测的最小的转录本长度  
   -B 用于下游Ballgown做差异分析  
   -A 用于输出Gene abundances文件
   -e 表示只对参考基因组注释文件中的转录组进行定量

   还有需要注意的是，如果输入的bam文件是来自于HISAT2比对的话，需要在HISAT2使用时加上–dta，这有利于stringtie的组装

2. 输出gtf格式文件，内容如下（可看官网的上的说明）：1.seqname ：Denotes the chromosome, contig, or scaffold for this transcript

   2.source ：The source of the GTF file, default stringtie  
   3.feature ：Feature type; e.g., exon, transcript, mRNA, 5'UTR)  
   4.start ：Start position of the feature (exon, transcript, etc), using a 1-based index  
   5.end ：End position of the feature, using a 1-based index  
   6.score ：A confidence score for the assembled transcript. Currently this field is not used, and StringTie reports a constant value of 1000 if the transcript has a connection to a read alignment bundle  
   7.strand ： If the transcript resides on the forward strand, '+'. If the transcript resides on the reverse strand, '-'  
   8.frame ：Frame or phase of CDS features. StringTie does not use this field and simply records a "."  
   9.attributes ：

   • gene_id: A unique identifier for a single gene and its child transcript and exons based on the alignments’ file name.
   • transcript_id: A unique identifier for a single transcript and its child exons based on the alignments’ file name.
   • exon_number: A unique identifier for a single exon, starting from 1, within a given transcript.
   • reference_id: The transcript_id in the reference annotation (optional) that the instance matched.
   • ref_gene_id: The gene_id in the reference annotation (optional) that the instance matched.
   • ref_gene_name: The gene_name in the reference annotation (optional) that the instance matched.
   • cov: The average per-base coverage for the transcript or exon.
   • FPKM: Fragments per kilobase of transcript per million read pairs. This is the number of pairs of reads aligning to this feature, normalized by the total number of fragments sequenced (in millions) and the length of the transcript (in kilobases).
   • TPM: Transcripts per million. This is the number of transcripts from this particular gene normalized first by gene length, and then by sequencing depth (in millions) in the sample. A detailed explanation and a comparison of TPM and FPKM can be found here, and TPM was defined by B. Li and C. Dewey here
3. 用Transcript merge mode合并所有样本组装的转录本，类似于cufflinks套件中的cuffmerge

   stringtie --merge -p 20 -o stringtie_merged.gtf mergelist.txt

   mergelist包含了所有样本的组装后的gtf文件

   还可以通过-m，-c，-F，-T以及-f等参数对组装后的转录本在merge时进行过滤

   如果需要加入参考基因组的注释文件，可以用-G参数加入

4. 后续可以接差异分析在stringtie组装时使用了-B/-b，-e参数，然后stringtie merge时使用-G参数，接着用Ballgown进行differential expression分析。并且在官网内还提到可以使用一个python脚本来从stringtie结果中提取read count，用于DESeq2和edgeR包做差异分析。。。但是没试过。。

Cuffcompare

这个软件是cufflinks套件中的一个，为什么要使用这个软件呢，是因为它有stringtie没有的功能，但是我觉得蛮实用的，就是用于预测新转录本或者对基因结构进行优化等（当然只能用于参考而已）。Cuffcompare是将组装后的转录本与参考基因组的转录本进行比较，从而对比对结果进行分类

1. 使用cuffcompare将组装转录本与参考基因组的进行比较

   cuffcompare -r ~/reference/genome/mm10/gencode.vM13.annotation.gtf -s ~/reference/genome/mm10/GRCm38.p5.genome.fa stringtie_merged.gtf

   -r 参考基因组的注释文件
   -s 参考基因组序列

2. 输出结果文件
   • cuffcmp.loci ：locus id的位置信息，个人理解为组装后的gene的local信息
   • cuffcmp.stats ：统计结果展示，看看就行了
   • cuffcmp.combined.gtf ：如果输入多个gtf组装文件的话，这个文件是将多个文件中的转录本并集输入到这个文件中
   • cuffcmp.stringtie_merged_new.gtf.refmap ：如果组装的gtf文件中的转录本either fully or partially match参考基因组转录本，则输入该文件中，也就是说组装后与参考基因组几乎完全匹配的转录本
   • cuffcmp.stringtie_merged_new.gtf.tmap :这个文件跟cuffcmp.stringtie_merged_new.gtf.refmap形式相似，但是输入的则是组装后与参考基因组较为近似的转录本
   • cuffcmp.tracking ：这个是比较重要的一个文件，里面含有cuffcompare软件对组装转录本分类定义的描述信息

   2.1 cuffcmp.tracking文件

   • 第1列 ：Cufflinks transfrag id，A unique internal id for the transfrag
   • 第2列 ：Cufflinks locus id，A unique internal id for the locus
   • 第3列 ：Reference gene id（如果没有则为”-“）
   • 第4列 ：Reference transcript id（如果没有则为”-“）
   • 第5列 ：class code，也就是组装转录本的分类信息

     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)

     理解cuffcmp.tracking文件中的分类信息，然后参考一些文献中的预测新转录本的阈值条件，即可筛选出潜在的新转录本（作为参考），我暂时是这么理解的