參考文獻(xiàn)：stringtie enables improved reconstruction of a transcriptome from rNA-seq reads
幫助文檔：http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual
另附參考文獻(xiàn)：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
注意：第三篇文章中有完整的hisat2-stringtie-ballgown的代碼扫倡，一定要看嫁审！
一捧搞、原理
摘要
1扇商、二代測序產(chǎn)生了大量的短讀段，對于轉(zhuǎn)錄組定量來說，通過將短讀段組轉(zhuǎn)成轉(zhuǎn)錄本是一個(gè)定量的方法
2、Stringtie應(yīng)用了起源于最優(yōu)化理論的網(wǎng)絡(luò)流算法，與可選擇的從頭組裝策略一起來將這些短讀段組裝成轉(zhuǎn)錄本
3耍共、與目前其他的轉(zhuǎn)錄本組裝軟件相比烫饼，stringtie具有更精準(zhǔn)的基因組裝效果以及更好的基因表達(dá)估計(jì)，同時(shí)通過它獲得的組裝好的轉(zhuǎn)錄本的數(shù)目也比其它軟件多划提。
背景
組裝目前遇到的問題
1枫弟、RNA-seq產(chǎn)生了大量150bp左右的read
2、人類基因組中90%的多外顯子蛋白編碼基因和30%的ncRNA都具有可變剪切體
3鹏往、外顯子可以在轉(zhuǎn)錄本間共享淡诗、由于旁系同源導(dǎo)致的模糊的read比對以及低表達(dá)的基因都會(huì)阻礙組裝的進(jìn)行
4、錯(cuò)誤組裝的轉(zhuǎn)錄本會(huì)進(jìn)一步干擾isoform的表達(dá)量的估計(jì)
5伊履、上述問題目前很多軟件都有解決方案韩容，如專注于轉(zhuǎn)錄本確定的Trinity、專注于表達(dá)定量的RSEM以及二者兼顧的Cufflink等
6唐瀑、但有研究發(fā)現(xiàn)群凶，上述軟件即使確定了一個(gè)轉(zhuǎn)錄本的所有外顯子，也很難把它們組裝成正確的isoform的形式哄辣，同時(shí)多isoform的表達(dá)和新的剪切位點(diǎn)也會(huì)干擾這些軟件的組裝
7请梢、目前解決轉(zhuǎn)錄組組裝主要有兩種策略1）reference/genome指導(dǎo)的組裝，事先需要完成序列比對工作 力穗，如果使用雙端測序的數(shù)據(jù)可以確定進(jìn)一步提高組轉(zhuǎn)成功率 2）沒有參考的從頭組裝 這種組裝可以幫助那些沒有參考基因組的區(qū)域完成組裝毅弧，但技術(shù)上實(shí)現(xiàn)更困難，因?yàn)榇嬖诙嗫截惢蚣易搴捅磉_(dá)水平上的變化以及可變剪切的影響当窗，這種方法在精確性上不如上一種因此多用于沒有參考基因組的物種的組裝
Stringtie
1够坐、Stringtie通過使用genome指導(dǎo)的組裝的方法與從頭組裝的概念結(jié)合的新方法來改善轉(zhuǎn)錄組組裝
2、Stringtie的輸入不僅可以是經(jīng)過比對的結(jié)果崖面，也可以是Stringtie通過從頭組裝read所得到的contig元咙，當(dāng)這兩種輸入都用到的時(shí)候，我們下面稱之為stringtie+SR
3巫员、對于很多使用參考基因組輔助組裝的方法庶香，組裝的的策略都是先對read進(jìn)行cluter，然后建立一個(gè)graph model來推測每個(gè)基因所有可能的isoform简识，最終通過不同的graph的解析方法得到對轉(zhuǎn)錄本的組裝結(jié)果
4脉课、有名的cufflinks用的是overlap graph，該模型中nodes代表fragment财异，如果兩個(gè)fragment存在overlap并存在兼容的剪切模式，則對應(yīng)的node連接起來唱遭。其解析方法為一種保守的算法戳寸，可以產(chǎn)生能夠解釋所有read的最少的轉(zhuǎn)錄本，盡管這種方法很吸引人拷泽，但是它沒有考慮到轉(zhuǎn)錄本的豐度并且對于某些isoform來說該方法沒有辦法組裝疫鹊！
5袖瞻、stringtie采用了組裝轉(zhuǎn)錄本和估計(jì)表達(dá)量同步進(jìn)行的方法，這不同于cufflinks的先組裝后定量的策略拆吆。
6聋迎、首先，stringtie將read聚成cluster枣耀，然后采用了splice graph霉晕，其中node代表外顯子或外顯子的一部分，path將graph中可能 的剪切現(xiàn)象都連起來捞奕，最終對每個(gè)轉(zhuǎn)錄本通過創(chuàng)建一個(gè)網(wǎng)絡(luò)流的方法牺堰，利用最大流算法（maximum flow algorithm）估計(jì)轉(zhuǎn)錄本的表達(dá)量
7、最大流的問題是最優(yōu)理論中的經(jīng)典問題颅围，但是目前還沒有應(yīng)用到轉(zhuǎn)錄本定量中伟葫。
8、與其它組裝軟件相比院促，stringtie具有很高的準(zhǔn)確性和新型isoform的發(fā)現(xiàn)能力筏养，其優(yōu)勢在于使用了網(wǎng)絡(luò)流算法，同時(shí)stringtie也支持將read從頭組裝成更長的片段常拓，這進(jìn)一步提高了其組裝的正確性
9渐溶、其另一個(gè)優(yōu)勢在于它的最優(yōu)化策略，它平衡了每次組裝中每條轉(zhuǎn)錄本的覆蓋度墩邀，這樣可以對組裝算法產(chǎn)生一定的限制掌猛，因?yàn)樵诮M裝基因組時(shí)，覆蓋度是很重要的一個(gè)參數(shù)因?yàn)樗枰挥脕硐拗扑惴级茫駝t組裝器可能將重復(fù)的片段錯(cuò)誤地堆疊到一起荔茬，相似地轉(zhuǎn)錄組裝也是如此，在isoform中的每一個(gè)外顯子需要有相似的覆蓋度竹海，如果忽略這個(gè)參數(shù)可能會(huì)產(chǎn)生一些保守但是錯(cuò)誤的轉(zhuǎn)錄本慕蔚，其中含有大量剪切位點(diǎn)的基因組裝起來尤其困難。
二斋配、操作說明
Input:輸入的文件必須是一個(gè)根據(jù)基因組位置排好序的BAM文件孔飒，可以是Tophat或Hisat2的輸出文件，在通過samtools排序
命令行為：
stringtie <aligned_reads.bam> [options]*
常用參數(shù)說明
| 參數(shù) | 描述 |
|---|---|
| -G <ref_ann.gff> | 使用注釋好的gtf文件輔助組裝艰争，在-e未設(shè)置的條件下坏瞄，輸出中包括注釋文件中的轉(zhuǎn)錄本和預(yù)測的新型轉(zhuǎn)錄本 |
| -o [<path/>]<out.gtf> | 輸出文件的名字，最好是全路徑甩卓，默認(rèn)輸出為標(biāo)準(zhǔn)輸出 |
| -l <label> | 為輸出的轉(zhuǎn)錄本設(shè)置前綴名鸠匀，默認(rèn)為STRG |
| -p <int> | 線程數(shù)，默認(rèn)為1 |
| -A <gene_abund.tab> | 對輸出的gtf統(tǒng)計(jì)基因表達(dá)量逾柿，并以一個(gè)tab分割的文件輸出缀棍，這里需要提交輸出的文件名 |
| -C <cov_refs.gtf> | 對輸出的gtf中屬于-G提交的參考gtf的轉(zhuǎn)錄本統(tǒng)一輸出到該文件宅此，這里需要提交一個(gè)文件名 |
| -B | 是否需要輸出Ballgown可以識(shí)別的文件,在-b設(shè)置的情況下，使用-o的路徑輸出 |
| -b <path> | 對Ballgown輸出的文件指定一個(gè)路徑保存 |
| -e | 我認(rèn)為是最需要注意的參數(shù)Ｅ婪丁父腕！只統(tǒng)計(jì)可以匹配-G提交的參考gtf中的轉(zhuǎn)錄本，不再對新的轉(zhuǎn)錄本做預(yù)測青瀑，這可以加快程序的運(yùn)行速度 |
| -m <int> | 對預(yù)測的轉(zhuǎn)錄本設(shè)置最小長度璧亮，默認(rèn)為200 |
stringtie --merge [options] gtf.list :轉(zhuǎn)錄組merge模式，在該模式下狱窘，Stringtie可以利用輸入的一個(gè)gtf list并將他們中的轉(zhuǎn)錄本進(jìn)行非冗余的整理杜顺。可以在處理多個(gè)RNA-seq樣本的時(shí)候蘸炸，由于轉(zhuǎn)錄組存在時(shí)空特異性躬络，可以將每個(gè)樣本各自的轉(zhuǎn)錄組進(jìn)行非冗余的整合，如果-G提供了參考gtf文件搭儒，可以將其一起整合到一個(gè)文件中,最終輸出成一個(gè)完整的gtf文件
| 參數(shù) | 說明 |
|---|---|
| -G <guide_gff> | 提供的參考gtf文件穷当，指導(dǎo)整合 |
| -o <out_gtf> | 輸出文件名，默認(rèn)是輸出到標(biāo)準(zhǔn)輸出中 |
| -l <label></label> | 輸出的轉(zhuǎn)錄本前綴名淹禾，默認(rèn)是MSTRG |
三馁菜、實(shí)例
1、當(dāng)對新型轉(zhuǎn)錄本有需求時(shí)
1）對每個(gè)樣本單獨(dú)進(jìn)行轉(zhuǎn)錄本預(yù)測
注意不要設(shè)置-e參數(shù)铃岔！
ls -d SRR*|while read id;do input=$id/$id'.sorted.bam' ; output=$id/$id'.gtf' ;stringtie -o $output -p 10 -G $gtf -B -l $id $input;echo $id ;done
2）merge
#做一個(gè)gtf.list
ls -d SRR*|while read id ;do find ./ -path './'$id'*' -name *.gtf ;done >gtf.list
#對所有g(shù)tf merge
stringtie --merge -p 10 -G $gtf -o total_merged.gtf -l merge gtf.list
3)利用merge得到的gtf重新對各個(gè)樣本做定量
注意這里一定要設(shè)置-e參數(shù)汪疮！
ls -d SRR*|while read id;do stringtie -e -A $id/$id'_gene_abund.tab' -C $id/$id'_cov_refs.gtf' -B -p 10 -G total_merged.gtf -o $id/$id'.merged.gtf' $id/$id'.sorted.bam' ;done
2、當(dāng)不需要預(yù)測新型轉(zhuǎn)錄本時(shí)
注意毁习，這里直接使用-e參數(shù)并且-G傳遞參考的gtf
ls -d SRR*|while read id;do stringtie -e -A $id/$id'_gene_abund.tab' -C $id/$id'_cov_refs.gtf' -B -p 10 -G $gtf -o $id/$id'.merged.gtf' $id/$id'.sorted.bam' ;done
