lncRNA组装流程的软件介绍之Stringtie

咱们《生信技能树》的B站有一个lncRNA数据分析实战，缺乏配套笔记，所以我们安排了100个lncRNA组装案例文献分享，以及这个流程会用到的100个软件的实战笔记教程！

下面是100个lncRNA组装流程的软件的笔记教程

StringTie 是用于 RNA-seq 的转录本组装和定量软件，StringTie 可以看做是cufflinks软件的升级版本，其功能和Cufflinks是一样的，包括下面两个主要功能：转录本组装和定量；相比Cuffinks, 其运行速度更快。

该软件的官网：https://ccb.jhu.edu/software/stringtie/index.shtml。

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中的每一个外显子需要有相似的覆盖度，如果忽略这个参数可能会产生一些保守但是错误的转录本，其中含有大量剪切位点的基因组装起来尤其困难。

一、软件安装

使用conda安装

conda install stringtie


二、输入文件介绍

输入的文件必须是一个根据基因组位置排好序的BAM文件，可以是Tophat或Hisat2的输出文件，在通过samtools排序

注意：一定要使用-dta选项来运行HISAT2，否则结果将会受到影响。

作为选项，可以向StringTie提供GTF / GFF3格式的参考注释基因组文件。在这种情况下，StringTie更喜欢使用注释文件中的这些“已知”基因，对于那些被表达的基因，它将计算coverage，TPM和FPKM值。它还会产生额外的转录本，而注释文件中并没有这些转录本。请注意，如果不使用选项-e，那么参考转录本就需要被reads 完全覆盖，以便包含在StringTie的输出中。在这种情况下，其他通过StringTie从数据中组装的转录本，且不在注释文件中的转录本也会输出。

三、stringtie的用法

安装完成以后，可以使用stringtie -h来查看软件的帮助文档。

1. 软件用法：

2. 常用参数：

-o [<path/>]<out.gtf> # 设置StringTie组装转录本的输出GTF文件的路径和文件名。此处可指定完整路径，在这种情况下，将根据需要创建目录。默认情况下，StringTie将GTF写入标准输出。

-p <int>  # 指定组装转录本的线程数（CPU）。默认值是1

-G <ref_ann.gff>  # 使用参考注释基因文件指导组装过程，格式GTF/GFF3。输出文件中既包含已知表达的转录本，也包含新的转录本。选项-B，-b，-e，-C需要此选项（详情如下）

--rf  # 链特异性建库方式：fr-firststrand(最常用的是dUTP测序方式，其他有NSR，NNSR).

--fr  # 链特异性建库方式：fr-secondstrand(如 Ligation,Standard SOLiD).

-l <label>  # 将<label>设置为输出转录本名称的前缀。默认：STRG

-f <0.0-1.0>  #将预测转录本的最低isoform的丰度设定为在给定基因座处组装的丰度最高的转录本的一部分。较低丰度的转录物通常是经加工的转录本的不完全剪接前体的artifacts。默认值为0.1。

-m <int>  # 设置预测的转录本所允许的最小长度.默认值为200

-A <gene_abund.tab> # 输出基因丰度的文件（制表符分隔格式）

-C <cov_refs.gtf> # 输出所有转录本对应的reads覆盖度的文件，此处的转录本是指参考注释基因文件中提供的转录本。(需要参数 -G).

-a <int>  # Junctions that don't have spliced reads that align across them with at least this amount of bases on both sides are filtered out. Default: 10

-j <float>  # 连接点的覆盖度，即设置至少有这么多的spliced reads 比对到连接点(align across a junction)。 这个数字可以是分数, 因为有些reads可以比对到多个地方。 当一个read 比对到 n 个地方是，则此处连接点的覆盖度为1/n 。默认值为1。

-t  # 该参数禁止修剪组装的转录本的末端。默认情况下，StringTie会根据组装的转录本的覆盖率的突然下降来调整预测的转录本的开始和/或停止坐标。

-c <float>  # 设置预测转录本所允许的最小read 覆盖度。 当一个转录本的覆盖度低于阈值，则输出文件中不含该转录本。默认值为 2.5

-g <int>  #设置ga最小值。Reads that are mapped closer than this distance are merged together in the same processing bundle. Default: 50 (bp)

-B  # 应用该选项，则会输出Ballgown输入表文件（* .ctab），其中包含用-G选项给出的参考转录本的覆盖率数据。（有关这些文件的说明，请参阅Ballgown文档。）
    如果选项-o 给出输出转录文件的完整路径，则* .ctab文件与输出GTF文件在相同的目录下。

-b <path>   # 指定 *.ctab 文件的输出路径, 而非由-o选项指定的目录。
    # 注意: 建议在使用-B/-b选项中同时使用-e选项，除非StringTie GTF输出文件中仍需要新的转录本。

-e  # 限制reads比对的处理，仅估计和输出与用-G选项给出的参考转录本匹配的组装转录本。使用该选项，则会跳过处理与参考转录本不匹配的组装转录本，这将大大的提升了处理速度。

--merge #转录本合并模式。在合并模式下，StringTie将所有样品的GTF/GFF文件列表作为输入，并将这些转录本合并/组装成非冗余的转录本集合。这种模式被用于新的差异分析流程中，用以生成一个跨多个RNA-Seq样品的全局的、统一的转录本。如果提供了-G选项（参考注释基因组文件），则StringTie将从输入的GTF文件中将参考转录本组装到transfrags中。(个人理解：transfrags可能指的是拼接成更大的转录本片段，tanscript fragments)。在此模式下可以使用以下附加选项：
-G <guide_gff> #参考注释基因组文件(GTF/GFF3)
-o <out_gtf>   #指定输出合并的GTF文件的路径和名称 (默认值：标准输出)
-m <min_len>   #合并文件中，指定允许最小输入转录本的长度 (默认值: 50)
-c <min_cov>   #合并文件中，指定允许最低输入转录本的覆盖度(默认值: 0)
-F <min_fpkm>  #合并文件中，指定允许最低输入转录本的FPKM值 (默认值: 0)
-T <min_tpm>   #合并文件中，指定允许最低输入转录本的TPM值 (默认值: 0)
-f <min_iso>   #minimum isoform fraction (默认值: 0.01)
-i             #合并后，保留含retained introns的转录本 (默认值: 除非有强有力的证据，否则不予保留)
-l <label>     #输出转录本的名称前缀 (默认值: MSTRG)


四、软件运行命令

1、当对新型转录本有需求时

1）对每个样本单独进行转录本预测

注意不要设置-e参数！

使用-e就不会有新转录本

批量处理

ls /home/data/lihe/lncRNA_project/04.mapping/*.bam > 1
ls /home/data/lihe/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config

cat > 05.stringtie.sh

config=$1
number1=$2
number2=$3
gtf=/home/data/lihe/reference/human/gtf/gencode.v37.annotation.gtf
cat $1 | while read id
do
    if((i%$number1==$number2))
    then
    arr=(${id})
    input=${arr[1]}
    sample=${arr[0]}
    stringtie -B -p 4 -G $gtf -o ./${sample}.stringtie.gtf -A ./${sample}.tab -l ${sample} $input
    fi    ## end for number1
    i=$((i+1))
done

for i  in {0..4}
do 
(nohup bash 05.stringtie.sh  config 5 $i 1>log.$i.txt 2>&1 & )
done


参数说明

-G $gtf # 提供的参考gtf文件，指导组装 
-o ./${sample}.stringtie.gtf # 输出文件名
-l ${sample} # 输出的转录本前缀名
$input # 输入bam文件
-A ./${sample}.tab # 对输出的gtf统计基因表达量，并以一个tab分割的文件输出，这里需要提交输出的文件名{sample}


2）merge

ls ../01.stringtie_gtf/*.gtf > mergelist.txt
gtf=/home/data/lihe/reference/human/gtf/gencode.v36.annotation.gtf
nohup stringtie --merge -p 8 -G $gtf -o ./stringtie_merged.gtf ./mergelist.txt > merge.log 2>&1 &

3)  利用merge得到的gtf重新对各个样本做定量

注意这里一定要设置-e参数！

ls /home/data/lihe/lncRNA_project/04.mapping/*.bam > 1
ls /home/data/lihe/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config

cat > 05.expression_stringtie.sh

config=$1
number1=$2
number2=$3
cat $1 | while read id
do
    if((i%$number1==$number2))
    then
    arr=(${id})
    input=${arr[1]}
    sample=${arr[0]}
    stringtie -e -B -p 4 -G ~/lncRNA_project/05.stringtie/02.merge_gtf/stringtie_merged.gtf -A ./${sample}_gene_abundances.tsv -o ./${sample}/${sample}.gtf  $input

fi    ## end for number1
    i=$((i+1))
done

for i  in {0..4}
do 
(nohup bash 05.expression_stringtie  config 5 $i 1>log.$i.txt 2>&1 & )
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


五、输出文件

主要输出文件有：

1、 GTF文件： 记录组装的转录本信息

2、 Tab文件： 记录基因丰度信息

3、 GTF文件：完全覆盖与参考注释基因组文件所匹配的转录本信息

4、 *.ctab文件：用于下游Ballgown软件做差异表达分析的输入文件

5、 GTF文件：在合并模式下，生成一个合并的GTF文件


1.GTF文件：记录组装的转录本信息

• seqname: 染色体，contig, 或 scaffold

• source: GTF文件的源文件。

• feature: 特征类型；如：exon, transcript, mRNA, 5'UTR。

• start: 开始位置，使用基于1的索引

• end: 结束位置，使用基于1的索引

• score: 组装的转录本的可信度分数。目前这个字段没有被使用，并且如果转录本 与a read alignment bundle

  有连接，则StringTie输出常数值1000。

• strand: 正向链：'+'；反向链：'-'.

• frame: CDS特征的 Frame or phase 。StringTie不使用该字段，只记录一个“.”。

• 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.

2.Tab文件：记录基因丰度信息

如果StringTie使用-A <gene_abund.tab>选项运行，则返回包含基因丰度的文件。

• Column 1 / Gene ID: The gene identifier comes from the reference annotation provided with the -G option. If no reference is provided this field is replaced with the name prefix for output transcripts (-l).
• Column 2 / Gene Name: This field contains the gene name in the reference annotation provided with the -G option. If no reference is provided this field is populated with '-'.
• Column 3 / Reference: Name of the reference sequence that was used in the alignment of the reads. Equivalent to the 3rd column in the .SAM alignment.
• Column 4 / Strand: '+' denotes that the gene is on the forward strand, '-' for the reverse strand.
• Column 5 / Start: Start position of the gene (1-based index).
• Column 6 / End: End position of the gene (1-based index).
• Column 7 / Coverage: Per-base coverage of the gene.
• Column 8 / FPKM: normalized expression level in FPKM units (see previous section).
• Column 9 / TPM: normalized expression level in RPM units (see previous section).

3.GTF文件：完全覆盖与参考注释基因组文件所匹配的转录本信息

如果StringTie与 -C <cov_refs.gtf> 选项一起运行（需要选项-G


4.*.ctab文件：用于下游Ballgown软件做差异表达分析的输入文件**

如果StringTie与-B选项一起运行，它将返回Ballgown输入文件，包含以下文件：
(1) e2t.ctab 
(2) e_data.ctab
(3) i2t.ctab
(4) i_data.ctab
(5) t_data.ctab。


5.GTF文件：在合并模式下，生成一个合并的GTF文件

如果StringTie使用--merge选项运行，它将多个GTF / GFF文件作为输入，并将这些转录本合并和组装成非冗余转录本集合。

文末友情推荐

与十万人一起学生信，你值得拥有下面的学习班：

• 数据挖掘（GEO,TCGA,单细胞）2021第4期
• 生信爆款入门-2021第4期