一、 基本介绍

StringTie是一个转录组组装定量软件。StringTie利用了最优化算法中研究比较多的网络流算法，并且兼顾从头组装方法，来将这些复杂的数据组装成为转录本。

有参考基因组的转录本拼接，首先将reads比对到参考基因组，然后根据定位的坐标信息和跨越内含子的spliced reads中反映的连接关系建立备选的拓扑图，最后设计相应的算法在拓扑图中选择合理的转录本形成最终的转录组数据集。

• Cufflinks根据每个片段在基因组上的回贴位置，建立重叠图（overlap graph），设计一个最小路覆盖模型，从图中找到能覆盖所有的点并且数目最少的路，作为最终预测的表达转录本集合。但是已有文献报道，最小路覆盖模型并不是最优模型，在很多情况下，最小路覆盖中的路并不代表表达的转录本，而某些表达的转录本并不出现在最小路覆盖集合中。

• StringTie改进了拼接的图结构，变成更为典型直观的剪接图（splicing graph），利用一个贪婪算法逐条拼接出转录本，并采用广义最大流模型对找到的每条转录本逐条估计表达量，之后根据表达量更新剪接图，迭代进行。与Cufflinks相比，StringTie速度更快，预测的转录本更多。

• 无参考基因组转录本的重构，如Trinity、OASES，都是利用了de Bruijn图。

二、 hisat2-stringtie流程

(1) 对于每个RNA-Seq样本，使用--dta选项，用hisat2将reads比对到参考基因组。使用samtools对每次HISAT2运行的SAM输出进行排序并转换为BAM。

(2) 对于每个RNA-Seq样本，运行StringTie将上一步获得的BAM组装起来。如果参考注释可用，建议使用-G选项运行StringTie。

(3) 使用--merge运行StringTie，以产生之前组装的RNA-Seq样本的一个非冗余转录本集合GTF。stringtie --merge模式接受所有以前的每个样本获得的组装转录本文件(GTF格式)的列表作为输入，以及一个参考注释文件(-G选项)。

(4) 对于每个RNA-Seq样本，使用-B/-b和-e选项运行StringTie，以估计转录丰度并生成Ballgown的read coverage tables。-e选项不是必需的，但建议在此运行中使用，以便对输入转录本产生更准确的丰度估计。在此步骤中运行的每个StringTie都将输入在步骤1、2中获得的相应样本的sorted read alignments (BAM file)，以及在步骤3中StringTie --merge生成的合并文本(GTF文件)作为-G选项输入。请注意，只有在这种情况下-G选项没有与参考注释一起使用，而是与在所有样本中观察到的全局、合并的转录本集合一起使用。此步骤相当于原始Ballgown pipeline中描述的Tablemaker步骤。

(5) 现在可以使用Ballgown来加载上一步中生成的coverage tables，并对差异表达执行各种统计分析，生成图等等。

(6) 如果对新的转录本不感兴趣，或者分析只针对一组已知的感兴趣的转录本，则可以采用另一种更快的差异表达分析工作流程。这个简化的流程只有3个步骤，因为它绕过了每个RNA-Seq样本的单独组装和转录本合并步骤。这个简化的工作流程试图直接估计和分析参考注释文件中给定的已知转录本的表达。

三、 使用方法

(1) 转录本组装

stringtie <in.bam ..> [-G <guide_gff>] [-l <prefix>] [-o <out.gtf>] [-p <cpus>] [-v] [-a <min_anchor_len>] [-m <min_len>] [-j <min_anchor_cov>] [-f <min_iso>] [-c <min_bundle_cov>] [-g <bdist>] [-u] [-L] [-e] [--viral] [-E <err_margin>] [--ptf <f_tab>] [-x <seqid,..>] [-A <gene_abund.out>] [-h] {-B|-b <dir_path>} [--mix] [--conservative] [--rf] [--fr]

stringtie <input.bam ..> [-G <guide_gff>] [-l <label>] [-o <out_gtf>] [-p <cpus>] [-e] [-A <gene_abund.out>] {-B | -b <dir_path>}

#组装
stringtie  -p 4  -G genome.gtf  -o sample_out.gtf  (-l sample)  sample.sorted.bam

<in.bam >：按基因组位置排序的bam文件。如hisat2 --dta得到的sam文件，经过samtools转换和排序后的bam文件。输入SAM文件中的每个剪接read比对（spliced read alignment，即至少跨一个junction的比对）都必须包含标记tag XS，以指明产生这个测序read来源的RNA的基因组链，由TopHat和HISAT2(--dta)生成的比对已经包含了这个标签。
<out.gtf>：是Stringtie的主要输出，一个包含组装好的转录本的GTF文件。GTF格式记录了转录本表达量信息，在attributes里面会给出3种表达量coverage、TPM、FPKM。
-p：线程数（默认1）。
-G：用于指导装配过程的参考注释（GTF/GFF）。给出参考基因组注释文件，StringTie将会更倾向于使用已知的转录本。对于表达的已知转录本，它将计算覆盖率、TPM和FPKM值。输出将同时包括给出的参考转录本以及新组装的转录本。-B、-b、-e、-C需要此选项。
-o：组装的转录本GTF的输出路径/文件名。会针对每个bam文件生成一个GTF文件。默认情况下，StringTie将GTF写入标准输出。
-l：输出转录本名称的前缀（默认STRG）。
-e：只对注释中给出的转录本进行定量（需要-G）。如果我们只关心已在GTF中注释出的转录本和基因，则可以使用-e参数指定不进行转录本的de novo组装，提升处理速度。（第二次定量时）建议使用。如果没有给出-e，参考转录本需要被reads完全覆盖才会包含在StringTie's output，并且新组装的、参考转录本中没有的转录本也会被输出。
-A：输出的基因丰度估计文件（gene abundance estimation output file）。输出的基因丰度信息文件（Coverage、FPKM、TPM值），以Tab分割，包括九列。

(2) 转录本整合

转录本合并模式，可以将所有样品的GTF/GFF文件列表作为输入，并将这些转录本合并/组装，生成一个统一的非冗余的转录本集合。输入的多个gtf文件合并成一个gtf，得到多个样品的全部转录本。这种模式被用于新的差异分析流程中，用以生成一个跨多个RNA-Seq样品的全局的、统一的转录本。（This mode is used in the new differential analysis pipeline to generate a global, unified set of transcripts (isoforms) across multiple RNA-Seq samples.）

该步骤为所有样本创建了一套统一的转录本，以便于下游在不同实验条件下计算所有转录本的差异表达水平。这一步是必须的，因为有些样本的转录本可能仅仅被部分reads覆盖，无法被第二步的StringTie组装出来。输出是一个合并的GTF文件，包含所有合并的基因模型，但是在覆盖率、FPKM和TPM上没有任何数值结果。然后，有了这个合并的GTF，StringTie可以重新估计丰度，使用-e选项对原始的比对文件集再次运行它。

注：StringTie 的merge模式能够合并不同的来源的结果，但在合并的同时会根据FPKM，TPM和转录本长度过滤，最终结果可以认为是在所有样本里面都是有所表达的基因，因此最终的数目会少一些。同时由于某些基因表达量低，单个样本里由于read数少无法覆盖基因，因此最终的预测结构还完整。因此，可以先将BAM合并后，然后用StringTie进行预测，如果为了输出结果的可靠性，还可以根据FPKM和TPM做过滤。

stringtie --merge  [Options]  { gtf_list | strg1.gtf ...}

#整合
stringtie --merge  -p 4  (-m 40  -F 0.4  -T 0.4)  -G genome.gtf  -o merge.gtf  mergelist.txt

-G <guide_gff>：reference annotation to include in the merging (GTF/GFF3)。如果提供了-G选项（参考注释），StringTie将把来自输入GTF文件的transfrags与参考转录组合在一起。
-o <out_gtf>：output file name for the merged transcripts GTF (default: stdout).
-l <label>：name prefix for output transcripts (default: MSTRG).
merge.gtf：从一系列的GTF文件中得到了一个融合的GTF文件。
mergelist.txt：包含上一步的sample_out.gtf列表。需要预先将要整合的GTF文件的文件名录入到mergelist.txt文件中，每行包含一个gtf路径。
-m <min_len >：要包含在merge中的最小输入transcript length (默认值50)
-c <min_cov >：要包含在merge中的最小输入transcript coverage (默认值0)
-F <min_fpkm >：要包含在merge中的最小输入transcript FPKM (默认值1.0)
-T <min_tpm >：要包含在merge中的最小输入transcript TPM (默认值1.0)
-f <最小iso >：最小isoform fraction (默认值0.01)
-g <gap_len >：要合并在一起的转录本之间的gap (默认值250)
-i：保留带有retained introns的merged transcripts；默认情况下，除非有强有力的证据，否则不会保留这些信息

(3) 转录本定量

#定量
stringtie  -p 4  -e  (-B)  (-A gene_abundance.out)  -G merge.gtf  -o sample.gtf  sample.sorted.bam

(4) prepDE.py

软件官网提供了一个Python脚本（prepd.py或prepDE.py3）来直接从StringTie生成的文件中提取reads计数信息（使用-e参数运行），以用于其他差异表达分析软件，如DESeq2和edgeR等。前期要进行到Ballgown前的一步，得到一个根据merge文件重新评估的GTF文件。prepDE.py需要python2，prepDE.py3用python2和3都行。

# counts
python prepDE.py  -i samplelist.txt  -g gene_count_matrix.csv  -t transcript_count_matrix.csv

-i：输入文件。samplelist.txt以制表符分隔的格式列出样本ID和GTF文件的路径，如sample1\t(即Tab)./sample1.gtf。
-g：输出基因raw count矩阵。
-t：输出转录本raw count矩阵。

四、 举个栗子

(1) 转录本组装

# stringtie_1.sh
for i in `ls ./04_bam/*.bam`
do 
    i=`basename $i`
    i=${i%.bam}
    echo $i >>log/stringtie_log.txt
    stringtie  ./04_bam/$i.bam  -p 4  -G ./annotation/genome/Drosophila_melanogaster.BDGP6.32.109.gtf  -o ./05_gtf_1/$i.out.gtf
done

(2) 整合上一步的转录本

ls 05_gtf_1/* > 05_gtf_2/mergelist.txt

stringtie --merge  -p 4  -m 40  -F 0.4  -T 0.4  -G ./annotation/genome/Drosophila_melanogaster.BDGP6.32.109.gtf  -o ./05_gtf_2/merge.gtf  ./05_gtf_2/mergelist.txt  &

(3) 转录本定量

# stringtie_3.sh
for i in `ls ./04_bam/*.bam`
do 
    i=`basename $i`
    i=${i%.bam}
    echo $i >>log/stringtie_log.txt
    stringtie  ./04_bam/$i.bam  -p 4  -e  -G ./05_gtf_2/merge.gtf  -o ./05_gtf_3/$i.gtf  -A ./05_gtf_3/$i.gene_abundance.txt
done

(4) 获得reads数量

ls 05_gtf_3/*.gtf > 05_gtf_4/samplelist.txt   #格式为sample1\t./sample1.gtf
cat -A 05_gtf_4/samplelist.txt

python /mnt/f/bioinfo_GEO/prepDE.py3  -i ./05_gtf_4/samplelist.txt  -g ./05_gtf_4/gene_count_matrix.csv  -t ./05_gtf_4/transcript_count_matrix.csv