转录组测序—数据分析与解读

转录组测序分析流程

流程概览

转录组测序的分析流程大致可以分成三类，包括基因组比对（Genome mapping）、转录组比对（Transcriptome mapping）、转录组组装（Reference-free assembly），见下图。其中第三种主要是用于分析没有参考基因组和基因注释的物种，应用场合较少且不适合新手入门。对于人、小鼠、大鼠等模式物种，通常用前两种方法进行分析。

1. Data Cleaning

用FastQ文件来储存测序数据，下图为文件的格式：

FastQ文件格式

每一条序列（read）包含四行，第一行是read的ID，第二行是序列，第四行是序列中每个碱基的测序质量。质量评估用FastQC，测序数据处理是：Trim Galore或Trimmomatics 可以自行选择。

2. 比对

2.1什么是比对

由于二代测序的reads长度通常介于50~300个碱基，因此即便使用双端测序，也基本不可能覆盖完整的mRNA转录本，因此想直接用FastQ文件从头分析测到了哪些转录本需要非常复杂的分析和计算。好在通常情况下，公共数据库已经提供了测序样品的基因组和转录本的序列。因此我们只需要知道，每一条reads来自哪一条转录本就可以了，这个将reads与参考（Reference）基因组/转录组的序列进行比较和匹配的过程，我们通常称之为“比对”（文献中提到的read alignment和mapping通常说的都是这个）。

基因组比对就是把reads比对到完整的基因组序列上，而转录组比对则是把reads比对到所有已知的转录本序列上。建议使用基因组比对的方法进行分析，理由如下：

① 转录组比对需要准确的已知转录本的序列，对于来自未知转录本（比如一些未被数据库收录的lncRNA）或序列不准确的reads无法正确比对；

② 与上一条类似，转录组比对不能对转录本的可变剪接进行分析，数据库中未收录的剪接位点会被直接丢弃；

③ 由于同一个基因存在不同的转录本，因此很多reads可以同时完美比对到多个转录本，reads的比对评分会偏低，可能被后续计算表达量的软件舍弃，影响后续分析

④ 由于与DNA测序使用的参考序列不同，因此不利于RNA和DNA数据的整合分析。

此外，值得注意的是，RNA测序并不能直接使用DNA测序常用的BWA、Bowtie等比对软件，这是由于真核生物内含子的存在，导致测到的reads并不与基因组序列完全一致（如下图所示），因此需要使用Tophat/HISAT/STAR等专门为RNA测序设计的软件进行比对。

比对结果会展示为BAM/SAM文件，其中BAM格式是SAM格式的二进制版本（请理解为压缩后的版本，用Samtools可以打开），打开之后长这样：

BAM文件

BAM文件中每行代表一条reads的比对信息，其中第一列是read的ID，第二列为FLAG（包括是否双端比对，比对位点是否唯一等信息），第三列为比对的染色体，第四列为比对的起始位置，第六列为CIGAR值，代表比对的具体方式（例60M2D80M代表60个碱基完美匹配+2个碱基缺失+80个碱基完美匹配）等等。

表达量计算

基因表达量最直接的分析手段就是计算比对到每个基因的reads有多少条，在转录组测序中，我们通常称这个数字为count。前文说过，使用基因组比对的方法，我们获得了每条read比对到基因组上的位置，因此，我们需要把基因组的位置对应到基因。这个提供基因位置的注释文件通常以GTF或GFF3格式呈现，GTF格式示例如下：

GTF文件格式

第一列是染色体编号，第三列是本行的特征（feature），如gene、transcript、exon、CDS等（实际上大多数情况下，计算表达量只要带exon的行就够了），第四列和第五列是基因组起始和终止位置，第七列是正负链，第九列是注释信息（可以包括类似基因ID、转录本ID、基因名等信息）。

GTF的作用：有GTF文件后，就可以利用注释信息计算每个基因/转录本/外显子比对了多少reads，从而获取counts值。值得一提的是，为提高比对的准确性，HISAT2和STAR等软件在比对的过程中就已经结合GTF文件中提供的转录本剪接信息进行了优化。

然而，由于每个测序样品的起始RNA量不同，文库量不同，测序数据量不同，不适合直接作为表达量用于样品间比较。因此，我们需要对counts值进行校正，或均一化。最早流行的是RPKM（Reads Per Kilobase per Million mapped reads）或在此基础上衍生的FPKM（Fragments Per Kilobase per Million mapped reads）。此外，近几年TPM（Transcripts Per Million）也越来越常用了，第三类常见的校正方法是使用TMM（trimmed mean of M values）这类方法通常是差异表达计算软件生成的（edgeR、DESeq等），并且符合负二相分布/泊松分布，因此非常适合直接用于后续分析。

我个人比较推崇DESeq等软件校正出来的Normalized Counts，一方面是评价比较高，一方面用来算差异表达的p值也更方便。但Normalized Counts最大的问题在于，每次加入新样本时都需要重新校正，对于多批次完成的大样本量研究、或数据库的维护非常麻烦。因此，作为TCGA数据首选的校正方式，TPM同样是一个不错的选择。此外，人工合成的Spike-in同样可以作为外参校正表达量。

#常用软件#

Counts值计算常用HTSeq和featureCounts，此外部分软件自带counts值计算，如RSEM、Salmon等。

TPM和RPKM用RSEM都能算，或者其实直接写个代码手算都可以。

TMM之类的校正有不少R包可以用，我一般用DESeq（DESeq1和DESeq2没区别）来计算，edgeR也可以。

Heatmap同样可以用R包画，pheatmap应该是里面最简单的，通常三到四行代码就能画一张最简单的带有聚类结果的heatmap。如果不擅长R的话，可以用Cluster和Treeview的组合。

差异表达分析

通常认为基因的counts值服从某种统计学的分布模型（负二相分布、泊松分布等），然后利用该分布模型分析每个基因的差异表达情况，并计算组间差异。计算差异后，软件会输出每个基因对应的Fold-change以及p值。

在大多数转录组测序文章中，通常使用q值或FDR来衡量表达差异的显著水平。这是使用多重检验校正控制错误率后的结果。举个例子，当我们使用p<0.05为基准时，那么每20个基因，就可能有一个单纯因为概率被认为p<0.05，而当我们对上万个基因分析表达差异时，必须排除这种因为概率导致的假阳性，因此需要使用Bonferroni校正等方法控制false discovery rate（FDR），并得到更加严格的q值。（当然，无论是p值还是q值，都只是基于统计学的计算，虽然p<0.05是一个约定俗成的阈值，但q值完全没有必要限制使用0.05的阈值。另外，我们样品的差异大小与样本本身性质相关，比如病人组织间的表达差异通常较大，而用比较温和的药物刺激细胞系的前后表达差异通常较小，因此在做后续筛选时建议根据实际情况设定阈值。）

得到Fold-change以及p值/q值后，就可以绘制火山图（Volcano Plot）和其它散点图，火山图中通常横轴代表Fold change，纵轴代表p值。具体参见：差异表达分析图标结果释义。

最后，在转录组测序中，唯p值/q值论是不可取的，一方面p/q值的可信区间的范围是随着|Fold-change|的增大而缩小的，因此只有|Fold-change|足够大时，才能确保差异分析的假阳性率足够小。况且在|Fold-change|太小的情况下，即使p/q值显著，这种差异也未必具有生物学意义。此外在筛选差异表达基因时，每个基因的表达水平也需要考虑，因为丰度越低的基因，受测序误差的影响也越显著（同样是1个count的误差，对于表达量为1和100的基因，影响显然是不同的）；虽然有些软件会低丰度的基因进行校正，但是低丰度基因的p值计算在大多数软件中都存在较大误差，况且一个本底表达量比较低的基因，生物学功能可能也不会太强（用qPCR也很难检测）。因此，建议大家手动去除丰度过低的基因，我个人习惯counts>50（12G数据量），也可以使用RPKM>1等条件。

#常用软件#

差异表达分析（都是R包）：

DESeq/DESeq2，大概是使用率最高的R包之一，但不支持非整数的counts值；

edgeR：这个有点老了，但还是有人在用，不太推荐；

EBSeq：RSEM使用的分析软件，可以使用非整数counts值进行分析；

此外还有sSeq、BaySeq等等使用率较高的R包，就不一一介绍了。

Cufflinks自带计算差异表达的Cuffdiff，可以直接使用FKPM值进行计算，但个人建议如果不是万不得已，还是不要使用Cuffdiff进行分析。

火山图和其它散点图可以用R包ggplot做（其实用Excel都能做），因为我自己不常做这个图，用R做也不是很麻烦，所以也没有专门找过可以直接出结果的软件。

转载：聊聊转录组测序 聊聊转录组测序