参考文献：stringtie enables improved reconstruction of a transcriptome from rNA-seq reads
帮助文档：http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual
另附参考文献：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
注意：第三篇文章中有完整的hisat2-stringtie-ballgown的代码，一定要看！
一、原理
摘要
1、二代测序产生了大量的短读段，对于转录组定量来说，通过将短读段组转成转录本是一个定量的方法
2、Stringtie应用了起源于最优化理论的网络流算法，与可选择的从头组装策略一起来将这些短读段组装成转录本
3、与目前其他的转录本组装软件相比，stringtie具有更精准的基因组装效果以及更好的基因表达估计，同时通过它获得的组装好的转录本的数目也比其它软件多。
背景
组装目前遇到的问题
1、RNA-seq产生了大量150bp左右的read
2、人类基因组中90%的多外显子蛋白编码基因和30%的ncRNA都具有可变剪切体
3、外显子可以在转录本间共享、由于旁系同源导致的模糊的read比对以及低表达的基因都会阻碍组装的进行
4、错误组装的转录本会进一步干扰isoform的表达量的估计
5、上述问题目前很多软件都有解决方案，如专注于转录本确定的Trinity、专注于表达定量的RSEM以及二者兼顾的Cufflink等
6、但有研究发现，上述软件即使确定了一个转录本的所有外显子，也很难把它们组装成正确的isoform的形式，同时多isoform的表达和新的剪切位点也会干扰这些软件的组装
7、目前解决转录组组装主要有两种策略1）reference/genome指导的组装，事先需要完成序列比对工作 ，如果使用双端测序的数据可以确定进一步提高组转成功率 2）没有参考的从头组装 这种组装可以帮助那些没有参考基因组的区域完成组装，但技术上实现更困难，因为存在多拷贝基因家族和表达水平上的变化以及可变剪切的影响，这种方法在精确性上不如上一种因此多用于没有参考基因组的物种的组装
Stringtie
1、Stringtie通过使用genome指导的组装的方法与从头组装的概念结合的新方法来改善转录组组装
2、Stringtie的输入不仅可以是经过比对的结果，也可以是Stringtie通过从头组装read所得到的contig，当这两种输入都用到的时候，我们下面称之为stringtie+SR
3、对于很多使用参考基因组辅助组装的方法，组装的的策略都是先对read进行cluter，然后建立一个graph model来推测每个基因所有可能的isoform，最终通过不同的graph的解析方法得到对转录本的组装结果
4、有名的cufflinks用的是overlap graph，该模型中nodes代表fragment，如果两个fragment存在overlap并存在兼容的剪切模式，则对应的node连接起来。其解析方法为一种保守的算法，可以产生能够解释所有read的最少的转录本，尽管这种方法很吸引人，但是它没有考虑到转录本的丰度并且对于某些isoform来说该方法没有办法组装！
5、stringtie采用了组装转录本和估计表达量同步进行的方法，这不同于cufflinks的先组装后定量的策略。
6、首先，stringtie将read聚成cluster，然后采用了splice graph，其中node代表外显子或外显子的一部分，path将graph中可能 的剪切现象都连起来，最终对每个转录本通过创建一个网络流的方法，利用最大流算法（maximum flow algorithm）估计转录本的表达量
7、最大流的问题是最优理论中的经典问题，但是目前还没有应用到转录本定量中。
8、与其它组装软件相比，stringtie具有很高的准确性和新型isoform的发现能力，其优势在于使用了网络流算法，同时stringtie也支持将read从头组装成更长的片段，这进一步提高了其组装的正确性
9、其另一个优势在于它的最优化策略，它平衡了每次组装中每条转录本的覆盖度，这样可以对组装算法产生一定的限制，因为在组装基因组时，覆盖度是很重要的一个参数因为它需要被用来限制算法，否则组装器可能将重复的片段错误地堆叠到一起，相似地转录组装也是如此，在isoform中的每一个外显子需要有相似的覆盖度，如果忽略这个参数可能会产生一些保守但是错误的转录本，其中含有大量剪切位点的基因组装起来尤其困难。
二、操作说明
Input:输入的文件必须是一个根据基因组位置排好序的BAM文件，可以是Tophat或Hisat2的输出文件，在通过samtools排序
命令行为：
stringtie <aligned_reads.bam> [options]*
常用参数说明
| 参数 | 描述 |
|---|---|
| -G <ref_ann.gff> | 使用注释好的gtf文件辅助组装，在-e未设置的条件下，输出中包括注释文件中的转录本和预测的新型转录本 |
| -o [<path/>]<out.gtf> | 输出文件的名字，最好是全路径，默认输出为标准输出 |
| -l <label> | 为输出的转录本设置前缀名，默认为STRG |
| -p <int> | 线程数，默认为1 |
| -A <gene_abund.tab> | 对输出的gtf统计基因表达量，并以一个tab分割的文件输出，这里需要提交输出的文件名 |
| -C <cov_refs.gtf> | 对输出的gtf中属于-G提交的参考gtf的转录本统一输出到该文件，这里需要提交一个文件名 |
| -B | 是否需要输出Ballgown可以识别的文件,在-b设置的情况下，使用-o的路径输出 |
| -b <path> | 对Ballgown输出的文件指定一个路径保存 |
| -e | 我认为是最需要注意的参数！！只统计可以匹配-G提交的参考gtf中的转录本，不再对新的转录本做预测，这可以加快程序的运行速度 |
| -m <int> | 对预测的转录本设置最小长度，默认为200 |
stringtie --merge [options] gtf.list :转录组merge模式，在该模式下，Stringtie可以利用输入的一个gtf list并将他们中的转录本进行非冗余的整理。可以在处理多个RNA-seq样本的时候，由于转录组存在时空特异性，可以将每个样本各自的转录组进行非冗余的整合，如果-G提供了参考gtf文件，可以将其一起整合到一个文件中,最终输出成一个完整的gtf文件
| 参数 | 说明 |
|---|---|
| -G <guide_gff> | 提供的参考gtf文件，指导整合 |
| -o <out_gtf> | 输出文件名，默认是输出到标准输出中 |
| -l <label></label> | 输出的转录本前缀名，默认是MSTRG |
三、实例
1、当对新型转录本有需求时
1）对每个样本单独进行转录本预测
注意不要设置-e参数！
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
#做一个gtf.list
ls -d SRR*|while read id ;do find ./ -path './'$id'*' -name *.gtf ;done >gtf.list
#对所有gtf merge
stringtie --merge -p 10 -G $gtf -o total_merged.gtf -l merge gtf.list
3)利用merge得到的gtf重新对各个样本做定量
注意这里一定要设置-e参数！
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、当不需要预测新型转录本时
注意，这里直接使用-e参数并且-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
网友评论