Stringtie进行转录本组装和定量

参考 http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual​​​​​​​

stringtie <aligned_reads.bam> [options]*

1. 分别对每个样本进行转录本预测

stringtie ${sample_id}.sorted.bam -o ${sample_id}.gtf -p ${threads} -G ${gtf} -B -l ${sample_id} 

# 可以批量处理，测序数据放在项目目录下不同的子文件夹下，以SRR数据为例,在项目文件夹下运行：
ls -d SRR*|while read id;do input=$id/$id'.sorted.bam'; output=$id/$id'.gtf' ;stringtie  -o $output -p ${threads} -G ${gtf} -B -l $id $input;echo $id ;done &

主要参数说明：

-G 参考基因组注释信息，GTF/GFF3格式，

-B enable output of Ballgown table files which will be created in the

    same directory as the output GTF (requires -G, -o recommended)

-l 生成的.gtf文件中transcript_id的前缀，默认为STRG，一般设置sample号。

注: .gtf注释数据要和比对时的索引基因组组要统一。hg38.refGene.gtf 不等同于Homo_sapiens.GRCh38.84.gtf，组装版本号也要一致。

2. 合并转录本

# 生成gtf.list文件
find ./SRR* -name *gtf > gtf.list

# merge
stringtie --merge  -p ${threads} -G ${gtf}  -o total_merged.gtf -l merge gtf.list &

注：StringTie将GTF/GFF文件列表作为输入，并将这些转录本合并/组合成一组非冗余的转录本。这种模式用于新的差异分析管道，以生成跨多个RNA序列样本的一组全局统一转录本（异构体）。

3. 根据全局统一转录本重新计算各个样本中转录本的表达量。

for id in `ls -d SRR*`;do stringtie ${id}/${id}.sorted.bam -e -A ${id}/${id}_gene_abund.tab -B -p 50 -G total_merged.gtf -o ${id}/${id}.merged.gtf;done &

主要参数说明：

 ${id}/${id}.sorted.bam 各个样本的排好序的bam文件

-e only estimate the abundance of given reference transcripts (requires -G)

-A 输出的转录本丰度text文件名

-B enable output of Ballgown table files which will be created in the

    same directory as the output GTF (requires -G, -o recommended)

-G 注意后面的文件是 total_merged.gtf

注：这次是根据全局统一转录本来计算转录本丰度，一定要加-e 参数

4. 预测指定转录本的丰度

for id in `ls -d SRR*`;do stringtie ${id}/${id}.sorted.bam -e -A ${id}/${id}_gene_abund.tab -B -p 50 -G ${gtf} -o ${id}/${id}.merged.gtf;done &

注：即使是人的转录本也没能充分注释，一般用上面的流程1-3，如果只关注参考注释中的转录本，可以加-e参数，一步完成。${gtf}为制定的转录本注释数据，一般为参考基因组的注释信息。

5. 比较gtf文件

Usage:

gffcompare [-r <reference_mrna.gtf> [-R]] [-T] [-V] [-s <seq_path>]

    [-o <outprefix>] [-p <cprefix>]

    {-i <input_gtf_list> | <input1.gtf> [<input2.gtf> .. <inputN.gtf>]}

## 比较软件得到的gtf文件和参考gtf的差异

# 单次比较
gffcompare -r -r ${gtf} -o ${id}/strtcmp ${id}/${id}.merged.gtf

# 批量比较
ls -d SRR*|while read id;do gffcompare -r ${gtf} -o ${id}/strtcmp ${id}/${id}.merged.gtf;done &