轉錄組的組裝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的使用

將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的組裝

輸出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
用Transcript merge mode合並所有樣本組裝的轉錄本，類似於cufflinks套件中的cuffmerge
```
stringtie --merge -p 20 -o stringtie_merged.gtf mergelist.txt
```
mergelist包含了所有樣本的組裝后的gtf文件

還可以通過-m，-c，-F，-T以及-f等參數對組裝后的轉錄本在merge時進行過濾

如果需要加入參考基因組的注釋文件，可以用-G參數加入
后續可以接差異分析

在stringtie組裝時使用了-B/-b，-e參數，然后stringtie merge時使用-G參數，接着用Ballgown進行differential expression分析。並且在官網內還提到可以使用一個python腳本來從stringtie結果中提取read count，用於DESeq2和edgeR包做差異分析。。。但是沒試過。。

Cuffcompare

這個軟件是cufflinks套件中的一個，為什么要使用這個軟件呢，是因為它有stringtie沒有的功能，但是我覺得蠻實用的，就是用於預測新轉錄本或者對基因結構進行優化等（當然只能用於參考而已）。Cuffcompare是將組裝后的轉錄本與參考基因組的轉錄本進行比較，從而對比對結果進行分類

使用cuffcompare將組裝轉錄本與參考基因組的進行比較

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

-r 參考基因組的注釋文件
-s 參考基因組序列

輸出結果文件
- 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文件中的分類信息，然后參考一些文獻中的預測新轉錄本的閾值條件，即可篩選出潛在的新轉錄本（作為參考），我暫時是這么理解的

補充一點，才發現其實還有一個軟件gffcompare也能做cuffcompare的工作，stringtie開發者也推薦使用

參考文章：
http://blog.csdn.net/hill_night/article/details/44829965
http://www.bio-info-trainee.com/2073.html

本文出自於http://www.bioinfo-scrounger.com轉載請注明出處