|
American Journal of Botany 99(2): 248–256. 2012. A COMPARISON OF STATISTICAL METHODS FOR DETECTING DIFFERENTIALLY EXPRESSED GENES FROM RNA-SEQ DATA
作为一个阶段性总结，该文出现得很及时。请看摘要： “RNA-seq技术正在快速革新基因组学研究，而RNA-seq数据的统计方法正在不断发展。适时地回顾和比较最近提出的统计方法可以提供一个有用的指南，以便选择合适的方法进行数据分析。人们对检测基因的差异表达的能力情有独钟。这里我们通过基于不同的分布模型或真实数据进行的一系列模拟，比较了四种近期提出的统计方法，edgeR，DESeq，baySeq和一个两步Poisson模型(TSPM)的方法。我们按照基因的显著性排序和假阳性率控制比较了这些方法检测差异表达基因的能力。所有进行比较的方法都用可免费获取的软件来实现。我们还讨论了这些软件当前可获取版本的可用性和功能。”
首先提纲挈领地回顾历史，综述研究背景，突出该文章的研究目的和重要性。 1、RNA-seq与微阵列比较 “相比于基于杂交的微阵列技术，RNA-seq有若干优势，包括更大的表达水平范围，更多的信息来检测等位基因特异的表达，新启动子，新亚型，更低噪音，更高通量。因此，RNA-seq正准备在未来几年取代微阵列技术成为研究基因表达的主要平台。” 2、RNA-seq实验与分析技术概要 “典型的RNA-seq实验中，一个RNA样品被转换成一个cDNA片段文库，然后在高通量商用平台上测序，这样的平台诸如Illumina的Genome Analyzer，Helicos BioSciences的HeliScope，Applied Biosystems的SOLiD，Pacific Biosciences的SMRT，以及Roche的454Life Sciences测序系统。原始数据由大量的DNA片段序列（称为reads）组成，这些数据要经历一系列的分析步骤。Oshlack等(2010)提供了一个极好的分析流程的评论，流程包括映射reads，汇总每个基因的reads计数，归一化和检测差异表达基因。其中的表1提供了每一步分析的软件列表。通常，RNA-seq研究产生的reads要基于映射到目标基因组或de novo组装的转录组的情况分配给基因或其他分类单元。有一些RNA-seq数据的基因表达水平定量方法仍在研究中。可变剪切转录本和亚型表达的复杂性使得其成为一个活跃的研究领域。因为亚型检测不是本文的重点，有兴趣的读者可参考Hiller等(2009)和Salzman等(2011)及其估计RNA-seq数据中的亚型丰度的参考文献。‘基因’是我们贯穿本文剩余部分都采用的一个广义术语，可以指一个基因模型的单外显子或外显子的子集合或所有外显子。基因表达用映射到一个基因的reads数来度量。因此，RNA-seq产生了一种基因表达的离散度量，这与微阵列技术中可被视为连续变量的荧光强度度量不同。因此，用于分析微阵列数据的统计方法不能直接应用，而迫切需要发展合适的统计方法来处理海量的RNA-seq数据。” 3、聚焦于差异表达 “检测跨处理/条件的差异表达基因是一个关键步骤，而且有时是RNA-seq数据统计分析的主要目标。差异表达基因的确定有助于我们阐明基因功能，当细胞响应不同的处理和条件时。此外，检测差异表达基因是聚类基因表达谱或检验基因集富集性的事先步骤。由于RNA-seq历史尚短并且在不断发展，目前还没有可用的标准方法基于这些数据来检测差异表达基因。很多统计工作者在为此而努力。一些文章已经发表，更多的可能还在研究中。” 4、文章研究内容与结果预告 “本文中我们首先回顾了目前可用的检测差异表达基因的方法，包括edgeR，DESeq，baySeq和一个基于两步Poisson模型的方法。我们提供了如何下载对应的包或代码以在R软件中应用这些方法的信息。然后我们在模拟真实数据的各种设置下通过模拟研究，比较了它们在基因的显著性排序上的表现。我们还检查了不同过程的假阳性率控制，对这样一个基因组数据分析中的高维检验问题来说是一个必要的步骤。” “我们的结果表明baySeq在低假阳性条件下具有最高的真阳性率。还发现相比之下，TSPM在样本大小为2时不能像其他方法一样有效执行。在这些方法的FDR控制中，我们发现真实的FDR事实上比期望值大很多。我们对这些方法的优劣的研究是透明的，这可能对科学家分析将来从RNA-seq研究中得到的数据是有用的。”
第二部分是基于RNA-seq数据检测差异表达基因的若干方法的回顾 “要检测基因的差异表达，统计假设检验已经准备好了。基于正态分布已经发展了很多统计方法用于归一化微阵列的基因表达度量。例如，广泛使用的修正t检验就是R包limma基于正态假设实现的。如前所述，RNA-seq技术产生了基因表达的离散度量，因此基于正态分布开发的统计方法不能直接应用。对数变换可能会使高度扭曲的离散RNA-seq数据更接近正态分布。不过，这就不得不添加一个任意的小数给一些样本中的那些零计数基因以完成对数转换。即便如此，转换后的数据仍可能不太吻合正态分布。研究者不再聚焦于寻求转换以使现存微阵列分析方法可以被应用，而是基于可直接对基因计数建模的离散分布开发出了一些方法。” 有三个离散概率分布被提出来对RNA-seq研究的计数数据进行建模：二项分布、Poisson分布、负二项分布。数学上可以证明：如果reads数充分大（对RNA-seq来说是真的），而且一条read映射到一个给定基因的概率足够小，那么二项分布可以用Poisson分布很好滴逼近。在早期使用单来源RNA的RNA-seq研究报道中，对大部分基因来说技术重复之间计数的分布用Poisson分布拟合得很好。但是Poisson分布的一个性质是方差等于均值。当有生物学重复时，RNA-seq数据会表现出比Poisson分布期望的更高的变异性，对相当多的基因来说方差可能超过均值。这种现象叫做过离散。对过离散数据，基于Poisson的分析容易因取样误差低估而产生高假阳性率。这里可用拟似然方法，因为它引入了一个关于方差的尺度因子以允许它不同于均值。假设一个负二项模型而不是Poisson模型是处理过离散数据的另一方法，因为负二项分布可以设定方差大于均值。因为在研究生物学上有意义的结果时生物学重复是至关重要的，我们希望所有试验都能设计成包含生物学重复。因此，我们只回顾目前能处理过离散的可用方法。 应该提到的是，为Serial analysis of gene expression（SAGE）数据而开发的方法可应用于RNA-seq数据分析。例如，edgeR方法最初就是为SAGE数据分析开发的，现在用于RNA-seq数据分析。对用于分析SAGE数据的方法的综合分析超出了本文范围，有兴趣的读者可参考Lu等(2005)，Robinson & Smyth（2008）和Baggerly等(2004). 基于Poisson分布 最近，Auer和Doerge(2011)提出了一种基于一个两步Poisson模型的方法。它们的原理是一些基因可能有过离散而另一些可能没有。因此，该方法的第一步是检验每个基因的过离散。如果检验表明有过离散，则在第二阶段应用一个准Poisson似然方法来检验差异表达。否则，在第二阶段应用基于Poisson模型的检验。他们对两个基因列表分别控制FDR，因为评估差异表达显著性的两种不同方法被应用于不同的基因。 Srivastava & Chen(2010)有一篇很有意思的文章，用了一个广义Poisson分布对位置水平的reads计数进行建模。GPseq实现的这个方法考虑了进行差异表达分析时潜在的位置偏倚，这与我们回顾的所有其他方法不同。我们可以获得的真实数据集均未提供位置水平的计数，而GPseq不能处理基因水平的计数。因此，我们的分析不包括GPseq。
基于负二项分布 有三种R包（edgeR，DESeq，baySeq）实现的方法是基于负二项模型的。edgeR所用方法是最先提出的，原本是为SAGE数据而开发的，SAGE数据可视为小规模的RNA-seq数据。负二项分布被用作Poisson分布的自然推广，只要加上一个散度参数而比Poisson分布允许额外的变异。因为重复很昂贵，这就导致了RNA-seq研究的很小样本量，所以散度参数的估计是一个很有挑战性的问题。Robinson和Smyth(2008)提出对所有基因使用一个公共散度来达到对散度参数的更好估计。如果散度参数(度量了相比于均值额外的变异)对于所有基因相同的假设成立，则公共散度参数可以非常精确地估计出来，因为很多数据可用于该估计。但是，所有基因有公共散度这一点在实践中很少是一个合适的假设。可能更好的策略是允许不同的基因有不同的散度参数，而这些散度参数的估计可以用一些合适的统计方法借助基因间的信息来改进。一些策略用在微阵列数据分析中，发展出很多检验借助基因间信息来更好地估计方差或平均表达与方差。这些检验在与不借助这些信息的检验进行比较时，被证明有很好的性能。遵循类似的策略，一种修正了的检验被提出来用于RNA-seq数据，可以用edgeR包来实现。 Anders和Huber(2010)也尝试了借助基因间信息来更好地估计散度参数。他们假设了一个方差和平均表达水平之间的局部线性关系。该假设允许使用具有相似表达水平的混合数据来估计方差（或等价地过离散参数）。这种方法在R/Bioconductor包DESeq中实现。edgeR和DESeq都提供基于精确检验或精确检验的逼近的检验p-value。 Hardcastle和Kelly (2010)提出的方法也假设数据服从负二项分布，但在显著性估计上与另两种方法不同。他们遵循一个经验Bayes方法，将基因按照后验概率估计进行排序，该模型对每个基因定义了差异表达。关于负二项分布参数的先验概率是利用数据找到的。相互之间表现相似的样本应该对潜在的基因参数具有相同的先验分布，而表现不同的样本应该具有不同的先验分布。该方法在R/Bioconductor包中已实现。
归一化 对于使用RNA-seq数据来比较样本间的表达，进行归一化来调整不同的测序深度和潜在的其他重复间的技术性影响，在上面的四种方法中都需要归一化。一个例子是用每个样本的总reads数和基因长度来归一化reads计数。该方法用RPKM来定量转录本水平。但是在进行样本间相同基因的差异表达分析而不在基因之间比较时，相对于基因长度的归一化就不重要了（也就是说，这种偏倚对不同样本中会以同样的方式影响相同的基因）。如果不考虑基因长度，目前可用的归一化方法可以用一些相对于平均表达水平的尺度因子来进行。最简单最常用的归一化因子是文库的总reads数，这是考虑到样本测序越深则每个基因分配到的reads越多。但是，通常reads总数主要是一小部分大量表达的基因贡献的。如果这小部分基因是差异表达的，则使用reads总数就会极大地影响检测差异表达的结果。Bullard等(2010)比较了几种归一化方法并发现在每个lane内使用非零计数分布的75th分位数作为归一化因子是一种更稳健的选择，相比于标准的总计数归一化来说，总体表现是所研究的这些方法中最优的。R包DESeq通过scaled counts的中位数来估计归一化因子，这是一种与75th分位数归一化类似的思想。R/Bioconductor包edgeR使用了对数表达比率的加权截断均值（trimmed mean of M values，TMM），这是另一种稳健的归一化方法。根据我们的经验，75th分位数和TMM方法的性能相似。
第三部分是模拟结果 模拟数据的好处是我们知道产生数据的真实潜在机制；因此，我们能评估结果，例如，一个宣称的阳性事实上是正确的（真阳性）。 为baySeq和TSPM，我们遵循Bullard等(2010)、估计归一化因子为计数的第三个四分位数。而DESeq和edgeR自有方法估计归一化因子。 我们基于两个标准比较和评价不同统计方法的分析结果。首先，看基因的显著性排序。选取不同的阈值得到TPR、FPR对，作出receiver operating characteristic (ROC) curve。其次，比较FDR。其中baySeq没有FDR，不参与此项比较。 文章中共有四次模拟，每次模拟都超过100个数据集，每个数据及不少于1万个gene。四次模拟各有侧重，分别是一半Poisson混合一半overdispersed Poisson模拟数据、基于玉米研究真实Illumina数据生成NB参数再产生模拟数据、同样基于玉米真实数据用经验分布得到NB参数、基于幼成淋巴细胞系真实数据模拟差异表达数据。 总体上都表明baySeq表现最好。
第四部分是数据分析 edgeR和DESeq吻合度较高（86%、91%），而与TSPM一致性较差（70?G未被DESeq和edgeR检测到）。
第五部分是讨论 模拟数据的分析能给事件中选择合适的方法分析RNA-seq数据提供有益指导。ROC曲线结果表明baySeq优于edgeR、DESeq和TSPM，而edgeR和DESeq性能相似，且与baySeq接近。TSPM方法可变性大，对小样本产生的结果不好，但是样本量增加时性能会改善很多。 FDR控制结果常常比真实的FDR要高。需要研究为什么没控制好，开发更好地FDR控制方法。建议实践者采用更严格的FDR控制方法来避免太多的假发现。 对于处理不同的实验设计，所用R包的灵活性各异。都能比较完全随机化设计的两组比较。而有一些提供允许更复杂的实验设计，还有的允许不同的估计模式。baySeq可以分析涉及多个处理组的实验设计。edgeR可用于两组或多组，至少一个组有重复。但是差异分析只能成对比较。edgeR有两种方法估计散度参数——普通的和tag层次的。除了DESeq可以没有重复外，所有方法都要求至少一个重复。但是不推荐无重复的实验设计。 鉴于这些选项，方法的选取还取决于实验设计。 |
|
|
