转录组数据分析笔记（5）——stringtie转录本组装和定量

本篇笔记主要记录如何用Stringtie做转录本的组装和定量，以及如何制作样本列表矩阵，为后面DESeq2分析做铺垫。

stringtie转录本组装和定量

1 转录本组装

Stringtie是一个基因和转录本组装定量的软件，stringtie的输入文件有两个，一个是经过排序的bam文件，排序可以用前面说到的samtools sort命令完成，还有一个是参考基因组的注释文件（gff或gtf格式）。

在使用Stringtie进行基因或者转录本组装定量的过程中，有一个非常重要的参数 - e，我之前跑了一遍流程没有加参数-e，结果组装的结果非常差，还有大量的未注释的基因。我请教了度娘，网上的教程攻略也都是抄来抄去的没解决什么问题，官网只有这么一句解释：

-e this option directs StringTie to operate in expression estimation mode; this limits the processing of read alignments to estimating the coverage of the transcripts given with the -G option (hence this option requires -G).

对于加了参数-e之后如何做的比对和组装处理还是不明了，不知道表达评估模式的原理是什么，只能自己做个大概的总结（不知正确与否）：

如果我们研究的样本没有很好的注释信息，研究的人少，现有的注释信息都不完善，那么我们就需要重建转录本进行注释，这个时候就不需要加参数-e。
如果样品的注释信息非常完整，比如拟南芥这种模式生物，我们不需要重建新的转录本进行注释，只对现有的参考基因组注释文件就足够了，那就要用-e参数，不需要预测新的转录本。

-e参数还有个比较重要的地方，只有用了-e参数后，才可以运行prepDE.py3脚本得到read count矩阵（也就是进行定量），这个脚本后面会说。

我们首先创建一个shell脚本进行转录本组装：

#!/bin/bash
for i in `seq 194 209`
do
    stringtie -p 4 -e \
    -G /media/sf_/data/ref/Arabidopis_thaliana.gtf \    # 参考基因组注释文件
    -o /media/sf_/data/fastq/gtf/ERR1698"$i".gtf \      # 自定义输出文件
    /media/sf_/data/fastq/bam/ERR1698"$i".bam   # 输入的bam文件
done

保存，运行，我们可以得到.gtf格式文件，less一下查看里面的内容：

我们这里因为加了参数-e，不会有新的基因和转录本，可以看到每个read比对上的基因的信息。（不加参数-e会组装新基因和转录本，默认采用STRG加数字编号进行区分）。每行数据会给出coverage，FPKM和TPM三个信息，后两者都可以用来定量。FPKM和TPM都是对read counts数目进行的标准化，如果是单端测序数据可以用RPKM进行标准化，不进行数据标准化的比较是没有意义的。

2 合并转录本(重构转录本才需要)

这一步要注意下，如果需要重构转录本才需要合并所有的转录本的组装结果，得到一个非冗余的转录本合集，也就是获得跨多个RNA-seq样品的全局的转录本。这里需要分两步：

$ ls *.gtf > mergelist.txt # 将所有组装的转录本文件名合并到一个文件

$ stringtie --merge -p 4 -G /media/sf_/data/ref/Arabidopis_thaliana.gtf -o merge.gtf ./mergelist.txt #这一步是用--merge指令将所有转录本合并输出到merge.gtf文件中

我们最后得到的merge.gtf就是全局的转录本。这里只是记录一下这步操作，我们只关注参考基因组的注释结果就不需要merge。

3 获得定量表达矩阵

DESeq2要求输入的定量结果为raw count形式，raw count是根据mapping到基因或转录本的reads数计算得到，而stringtie只提供了转录本水平的表达量，定量方式包括TPM和FPKM值两种。为了进行raw count定量，stringtie官方提供了prepDE.py脚本（两个版本，我选择的python 3版本，在我base环境下不会冲突），可以计算出raw count的表达量。

下载这个python脚本，如果你用的是windows浏览器，在官网找到脚本直接右键复制链接，用wget直接下到linux系统里，千万不要在windows上直接复制粘贴代码过去。因为windows的换行符和linux的不一样，两个系统间直接粘贴代码会出现错行和莫名其妙的缩进导致程序报错（可以用cat -A看两个系统换行符的区别，血的教训，排查了老半天才发现）推荐用prepDE.py3，不用再切python 2 的环境了。

官方给出的prepDE.py脚本有两种运行方式（如下图所示），一种是建立Ballgown能识别的目录结构，一种是建立sample_lst文件并指定所有样品数据的路径。两种方法都可行，Ballgown现在用的比较少，比较主流的还是Stringtie+DESeq2的分析方法。演示一下如何创建sample_lst和解释一下这个文件要求的格式。

3.1 sample_lst文件准备

简单来说，sample_lst.txt要求第一列为样品编号，第二列为对应编号的样品gtf文件所在路径，中间用制表符tab隔开，如下图（命名不一定要完全一样，注意格式，后面要导入prepDE脚本，能找到就行）：

这个文件准备工作比较简单，不再赘述

3.2 运行prepDE.py3

将prepDE.py3脚本放在上面gtf文件的目录下，运行以下命令：

$ python prepDE.py3 -i sample_lst.txt -g gene_count_matrix.csv -t transcript_count_matrix.csv

解释一下：

参数含义

-i # 输入文件，就是前面做的sample_lst.txt
-g # 自定义基因组表达矩阵名字，默认也是gene_count_matrix.csv
-t # 自定义转录本表达矩阵名字，默认也是transcript_count_matrix.csv

得到的这两个文件就是基因和转录水平的raw count表达量矩阵，我们都可以用于后面的DESeq2分析。

4. 制作样本列表矩阵

这里需要和前面为了运行prepDE.py脚本而制作的sample_lst文件区分开，要做下一步DESeq2差异基因分析，我们需要自己手动创建一个DESeq2能识别的样本列表矩阵，包含两列信息：一列是样本名称，一列是样本分组。样本分组信息我们可以直接从下载样本数据的地方（EBI官网）得到，只需要自己改一下格式。