RNA-seq数据差异表达分析方法的比较

RNA-seq数据差异表达分析方法的比较

BMC Bioinformatics 2013, 14:91 doi:10.1186/1471-2105-14-91

A comparison of methods for differential expression analysis of RNA-seq data

Charlotte Soneson (Charlotte.Soneson@isb-sib.ch) Mauro Delorenzi (Mauro.Delorenzi@unil.ch)

摘要说的背景是：

“发现条件间差异表达的基因是理解表型变异的分子基础的一个有机部分。过去几十年中，DNA微阵列被广泛用于定量不同基因的mRNA丰度，更近期的RNA-seq作为一个强有力的竞争者冒了出来。随着测序成本持续下降，可以想象使用RNA-seq做差异表达分析会迅速增加。为了探索可能性和解决这种相对新型的数据提出的挑战，大量软件包特别为RNA-seq数据的差异表达分析开发出来了。”

而本文的结果是：

“我们广泛比较了RNA-seq数据的差异表达分析的7种方法。所有方法都可以在R框架下免费获得，并以一个计数矩阵作为输入，计数即每个样品中映射到每个感兴趣的基因组特征上的reads数目。我们基于模拟数据和实际RNA-seq数据评价了这些方法。”

结论就是：

“极小样本量仍是RNA-seq实验的普遍状况，对所有评价方法造成了困难；而任何在这样的条件下获得的结果都应该谨慎解释。对于更大的样本量，组合稳定方差变换和limma方法来进行差异表达分析会在很多不同的条件下表现良好，正如非参数的SAMseq方法一样。”

到2013年还说这种话，这些结论实在有点鸡肋啊~ 貌似为SAMseq摇旗呐喊来的……不过：

比较了11种软件包，这还是前所未有的：DESeq、edgeR、NBPSeq、TSPM、baySeq、EBSeq、NOISeq、SAMseq、ShrinkSeq这9种可直接处理计数数据，另两种分别是voom(+limma)和vst(+limma)，转换数据后用limma做差异表达分析。

正如很多文章已经提到的那些，RNA-seq比起微阵列有三大优点：

1、更大的动态范围

2、更低的背景噪音

3、能检测和定量先前未知的转录本及亚型

RNA-seq也有一些难题：

1、样本内不均匀性：基因组区域之间核苷酸组成的变异性导致沿基因组的read覆盖深度并不均匀；

2、同样表达水平下，长基因比短基因有更多的reads；

3、对于条件之间的表达差异，分别对各个基因进行差异表达分析，而忽略了样本内的偏倚（它们被假设对所有样本有类似的影响）

4、样本间不均匀性：测序深度或文库大小

5、少数高表达基因抑制了其他基因的read计数比例，可能导致低表达基因的差异表达假阳性

相应的解决办法是：

1、对上述4,5两点，估计样本特异的归一化因子，用于重新调整观测计数。这些归一化方法是为了使样本间的非差异表达基因的归一化计数是相似的。本研究中使用的是DESeq 包中的TMM方法。归一化因子和TMM两种方法的性能相似，也是仅有的两个在文献9的评价中对所有度量都提供了满意结果的方法。

2、对于2、4两点，基因长度和文库大小，有的软件采用了RPKM或相关的FPKM方法。只有非参数方法才适用于RPKM值。

致命的假设：大部分基因在样本间的表达是相等的。于是差异表达基因分成上调、下调及其之间的或多或少相等的三部分。

差异表达已经提出的模型中，Poisson分布和负二项分布最常见，还有β分布也被提出来。Poisson分布很适合技术变异的分析；而更高的生物学变异需要合并过离散的分布，如负二项分布。

RNA-seq数据的差异表达分析仍处于它的婴儿期，需要不断提出新的方法来。目前没有一致认可的表现最好的方法，只发表了很少的现有方法的广泛比较。文献19中，依照区分真实差异/非差异表达基因的能力比较了四种参数方法。还有研究评价了测序深度的影响，并与样本量进行了比较，并发现后者的影响相当大。

本文中比较了为不同条件下RNA-seq数据的差异表达分析开发的11种方法。其中9种直

接对计数数据进行建模，而另两个先对计数进行变换再应用微阵列数据的差异表达分析的传统方法。研究限于R框架下实现的可应用于计数矩阵的可用方法。进一步我们聚焦于发

现两条件之间的差异表达基因，因为这是最常见的应用，虽然大多数方法也允许更复杂的试验设计。

对NB和Poisson分布模拟的数据和分别加了一些例外点的数据共四种数据集，研究了在不同实验条件下方法的下列方面：

1、排序真实DEGs在nonDEGs之前的能力；

2、在给定水平控制I型错误和假发现率的能力；

3、计算时间。

对于真实RNA-seq数据，比较了DEGs集，各自的数目和重叠的数目。还研究了不同方法获得的基因排序的一致性。

六种方法有名义p-value（edgeR、DESeq、NBPSeq、TSPM、voom+limma、vst+limma），我们定义分数值为1- pvalue。对于SAMseq，定义平均Wilcoxon统计量的绝对值为排序分值，而对于baySeq、EBSeq、ShrinkSeq使用估计的差异表达后验概率，或者等价地，

1-BFDR，其中BFDR表示估计的Bayesian FDR。对于NOISeq，使用统计量q_NOISeq。所有这些分数都是双侧的，即不被差异表达的方向影响。给定上述分数的阈值，我们就选择出阈值以上的为DEGs，其余的为non-DEGs。

SAMseq使用了重抽样策略使文库大小一致，因此隐含假设了所有归一化因子是相等的，在整体上调的模拟研究中表现最好。上下调基因都有的时候，所有方法的AUC（Area under the ROC curve (AUC)，Receiver Operating Characteristic）性能类似。TSPM和EBSeq

在所有方法中表现出最强烈的样本量依赖性，其次是SAMseq和baySeq。对于最小样本量（每条件下2样本），最佳结果是DESeq、edgeR、NBPSeq、voom+limma和vst+limma。

当所有DEGs上调时，baySeq结果的变异性很高；而DEGs向不同的方向调整时这种变异就会减小。

将真实DEGs排的很靠前的方法是基于变换的voom+limma和vst+limma方法和ShrinkSeq，但是TSPM和NOISeq还会把一些真实的nonDEGs排的很靠前。SAMseq也表现不错，但是有一些真实DEGs和nonDEGs被返回靠前的相同的值。

……实在是相当繁琐，直接看结论吧……

Conclusions

文本评价和比较了11种RNA-seq数据的差异表达分析方法。主要结果如下：

Table 2 Summary of the main observations

DESeq - Conservative with default settings. Becomes more conservative when outliers are introduced.

- Generally low TPR.

- Poor FDR control with 2 samples/condition, good FDR control for larger sample sizes, also with

outliers.

- Medium computational time requirement, increases slightly with sample size.

edgeR - Slightly liberal for small sample sizes with default settings. Becomes more liberal when outliers are introduced.

- Generally high TPR.

- Poor FDR control in many cases, worse with outliers.

- Medium computational time requirement, largely independent of sample size.

NBPSeq - Liberal for all sample sizes. Becomes more liberal when outliers are introduced. - Medium TPR.

- Poor FDR control, worse with outliers. Often truly non-DE genes are among those with smallest p-

values.

- Medium computational time requirement, increases slightly with sample size.

TSPM - Overall highly sample-size dependent performance.

- Liberal for small sample sizes, largely unaffected by outliers.

- Very poor FDR control for small sample sizes, improves rapidly with increasing sample size.

Largely unaffected by outliers.

- When all genes are overdispersed, many truly non-DE genes are among the ones with smallest p-

values. Remedied when the counts for some genes are Poisson distributed.

- Medium computational time requirement, largely independent of sample size.

voom / vst

- Good type I error control, becomes more conservative when outliers are introduced.

- Low power for small sample sizes. Medium TPR for larger sample sizes.

- Good FDR control except for simulation study B04000. Largely unaffected by introduction of outliers.

- Computationally fast.

baySeq - Highly variable results when all DE genes are regulated in the same direction. Less variability when the DE genes are regulated in different directions.

- Low TPR. Largely unaffected by outliers.

- Poor FDR control with 2 samples/condition, good for larger sample sizes in the absence of outliers. Poor FDR control in the presence of outliers.

- Computationally slow, but allows parallelization.

EBSeq - TPR relatively independent of sample size and presence of outliers.

- Poor FDR control in most situations, relatively unaffected by outliers.

- Medium computational time requirement, increases slightly with sample size.

NOISeq - Not clear how to set the threshold for qNOISeq to correspond to a given FDR threshold.

- Performs well, in terms of false discovery curves, when the dispersion is different between the

conditions (see supplementary material).

- Computational time requirement highly dependent on sample size.

SAMseq - Low power for small sample sizes. High TPR for large enough sample sizes. - Performs well also for simulation study B04000.

- Largely unaffected by introduction of outliers.

- Computational time requirement highly dependent on sample size.

ShrinkSeq - Often poor FDR control, but allows the user to use also a fold change threshold in the inference procedure.

- High TPR.

- Computationally slow, but allows parallelization.

没有哪种单独的方法对所有情形都是最优的，特定情形下方法的选择取决于实验条件。本文评价的这些方法中，基于稳定方差的变换与limma组合的方法在很多情况下都表现不错，而且不受例外点影响、计算很快，但是要求每条件下至少3个样本来提供充分的检定力。而且在两条件下散度不同时表现更糟糕。非参数方法SAMseq在大样本量时是性能最优的方法，需要至少每条件下4-5个样本提供充分的检定力。对于高表达基因，SAMseq的统计

显著性所需的倍数变化比很多其他方法要低，这可能潜在地折中了一些统计显著的DEGs 的生物学显著性。对ShrinkSeq也是一样，不过它有一个选项在推断过程中强加一个倍数变化要求。

小样本导致一些方法的误报率远超FDR阈值。对于参数方法，这可能是因为均值和方差估计不精确。TSPM受样本量影响最大，可能因为它使用了渐进估计。尽管发展指向大样本量，而且barcoding和multiplexing创造了固定成本分析更多样本的机会，但是目前为止RNA-seq实验仍然太贵而不允许广泛的重复。本研究所传达的结果强烈建议小样本差异表达基因应该谨慎解释，真实FDR可能超出所选FDR阈值数倍。

DESeq、edgeR和NBPSeq基于类似的原理，因此基因排序的精确度很类似。但是相同阈值选取出的DEGs有很大不同，这是因为它们估计散度参数的方法不同。在缺省设置和合理的大样本量下，DESeq通常过于保守而edgeR和NBPSeq通常过于慷慨而得出大量假DEGs。分析表明参数选择影响很大，而且缺省推荐参数事实上选择的很好通常能得到最佳结果。

EBSeq、baySeq、ShrinkSeq使用了不同的推断方法来估计每个基因差异表达的后验概率。baySeq一些条件下表现不错，但是高度可变，特别是所有基因都上调或都下调时。大样本量条件下有异常值时，EBSeq比baySeq的假阳性低，小样本量时baySeq比EBSeq的假阳性低。