评估标签丰度差异的适度统计检验

摘要

动机：数字基因表达（DGE）技术通过计算序列标签来测量基因表达。它们是在基因组尺度上测量基因表达的敏感技术，无需事先了解基因组序列。随着DNA测序成本的降低，预计DGE数据集的数量将急剧增加。

对于使用二项式、泊松、负二项式或伪似然（PL）模型进行计数的重复DGE数据，已经提出了各种差异表达测试，但当重复次数非常少时，这些都不可用。

结果：我们开发了使用负二项分布的测试来模拟相对于泊松的过度分散，并使用条件加权似然来调节跨基因的过度分散水平。我们的策略不仅适用于最少数量的库，而且在有更多库可用的情况下，它也比以前的策略更强大。该方法同样适用于其他计数技术，例如蛋白质组光谱计数。

可利用性：可以从以下位置访问R包http://bioinf.wehi.edu.au/resources网站/

联系人： smyth@wehi.edu.au

补充信息： http://bioinf.wehi.edu.au/resources网站/

1简介

1.1 DGE技术

DGE技术通过生成序列标签来测量基因表达。目前已有多种DGE技术，包括基因表达的系列分析（SAGE；Velculescu等。,1995)，大规模并行签名测序（MPSS；Brenner等。,2000)、合成测序（SBS、Margulies等。,2005)基因表达的polony多重分析（PMAGE，Kim等。,2007). DGE的负担能力似乎有望在未来几年实现突破。同样旨在生产1000美元基因组的技术也可能应用于表达谱研究，因为直到最近，成本限制的步骤一直是测序。此外，有许多有前途的基于序列分析的方法可用于基因表达谱分析或在不久的将来（Shaffer，2007).

现在有几种不同的技术（王，2007)，但本质上DGE系统可以量化基因表达的快照，而不需要知道基因序列或设计探针序列，这是微阵列所需要的。在SAGE的情况下，信使RNA（mRNA）从感兴趣的样本（例如癌症肿瘤）中提取，反转录成cDNA，由酶片段化成小片段标签从成绩单中的固定位置。这些标签的长度为10–20 bp，具体取决于协议。然后对标签进行测序，要么将它们连接起来并对其中的一段进行测序，要么并行测序。对于每个标签，在图书馆计数越大，表达越高。在存在mRNA数据库或基因组的地方，标签可以映射到基因组中的特定mRNA或位置。

SAGE最初用于测定胰腺（Velculescu）中表达的转录物等。,1995). 自那时以来，SAGE及其变体已成功应用于许多领域，包括建立人类癌症基因表达数据库（Lal等。,1999)发现癌症的预后因素（Aung等。,2006)以及创建小鼠组织表达图谱（Siddiqui等。,2005).

像微阵列一样，许多基于序列分析的技术除了转录谱分析之外还有其他应用，我们希望这里开发的方法在其他地方也有应用。例如，量化微小RNA（称为miRAGE）（康明斯等。,2006)，拷贝数分析（陈等。,2002)，全基因组DNA甲基化分析（Hu等。,2005)染色质占有率系列分析（SACO）（Impey等。,2004).

我们将注意力集中在推断两组库之间差异表达的问题上（例如癌症与正常），假设复制最小（至少一个类有多个样本）。由于测序的成本，传统的SAGE既费力又昂贵。即使在高通量测序的最新发展中，通常大多数“不动产”都用于对更多标签进行测序，而不是对更多库（样本或复制品）进行测序。因此，很少有大量的库可以进行比较。因此，统计分析方法必须在小样本中保持稳定。

在本文中，我们引用DGE数据。然而，我们的方法应该同样适用于此类其他计数数据，例如质谱数据Lu中的肽计数等。(2007). 我们已经优化了对两类比较问题的计算，但扩展到许多类或计算协变量是非常简单的。

许多全基因组统计推断方法共享所有基因的信息，包括特别的方式（Tusher等。,2001)或通过层次模型（Smyth，2004)，已被证明比标准方法更敏感。据我们所知，这是第一次探索用于标签计数数据差异表达分析的适度测试统计。我们方法的新颖之处在于，我们共享所有标签的信息，以稳定小样本中的分散估计。

2准备工作

2.1多个DGE库之间的差异表达

早期用于多个库之间差异表达的方法涉及将每个类中的库合并，并使用标准的两样本比例差异测试或Fisher精确测试。如前所述（Baggerly等。,2003,2004、卢等。,2005)但这种池并没有充分处理类内的可变性，因此提出了更灵活的模型。后来的方法（Ryu等。,2002)计算出的两个样本t吨-统计比例，从而考虑到图书馆之间的差异。然而，t吨-当数据真正非正常时，对于非常小的样本的统计可能会有问题。

标签统计模型的自然选择计数可以是泊松或二项式。实际上，泊松分布或二项式分布的均值-方差关系可能无法提供足够的灵活性。通常会观察到比模型允许的更多变化，即过度分散。因此，最近的方法探索了β二项式（Baggerly等。,2003)[以及更广泛地说，过度分散的物流（巴格利等。,2004)]和过分散对数线性（即伽马-泊松或负二项）模型（Lu等。,2005). （Lu）的模拟研究等。)建议负二项（NB）假设即使在非NB抽样情况下也是可靠的，因此应该为实际数据提供更灵活的框架。因此，我们将我们的模型与Lu的模型进行了比较等。(2005).

2.2统计框架：负二项模型

为了便于记谱，我们首先考虑单个标记。让Y（Y）_伊吉表示类的观察计数我和库j个用于特定标记。在这里j个= 1, …,n个_我现在，我们假设只进行两组比较，以便我= 1, 2. 我们分析的一个特点是，我们只需要其中一个n个₁或n个₂大于1。严格来说，以前的方法（Baggerly等。,2003,2004; 卢等。,2005)可以在此设置下操作。然而，在两个库与一个库的极端情况下，一次一段的推理需要从三个观测值中估计三个参数，这是一个相当徒劳的练习。

设λ_我是类RNA中该标签的真实相对丰度我。然后μ_伊吉=米_伊吉λ_我哪里米_伊吉是样本的库大小j个为了评估相对丰度的差异，零假设H（H）₀: λ₁=λ₂根据双面选择进行测试，并对每个标签重复此操作。

2.3色散估计

在库大小不相等的实际情况下，计数分布不一致，条件参数也不准确。罗宾逊和史密斯(2007)使用分位数调整将观察到的计数向上或向下调整，具体取决于相应的库大小是否低于或高于几何平均值（分位数调整的条件最大似然称为qCML）。这将创建大致相同的分布伪数据可以插入方程式中的(1)，对所有标签进行求和，并针对φ进行最大化，从而得出一个通用的估计值。即使只有100个标签，qCML估计值在一组常用估计值中在广泛的条件下也是偏差最小的（Robinson和Smyth，2007).

2.4统计测试

为了测试两种条件之间表达式的差异，我们在下面比较了两种统计测试。默认情况下，我们使用Lu中使用的Wald测试等。(2005). Wald测试简单地划分根据其估计的标准误差。其次，我们使用我们之前开发的精确测试（Robinson和Smyth，2007).

简单地说，精确测试的工作原理如下。用于将标签计数调整为用于估计的公共库大小的相同分位数调整用于构造精确测试。使用此伪数据，我们再次使用独立且同分布的NB随机变量之和也是NB这一事实。通过对总伪和（一个NB随机变量）进行条件处理，我们可以计算观察计数的概率为或大于我们观察到的计数，从而得出精确的P（P）-值。

3通过加权似然法进行中度弥散估计

3.1加权条件似然框架

共同分散假设，如罗宾逊和斯密思(2007)与标记估计相比，提供了显著的稳定性，尤其是在非常小的样本中。然而，通常情况下并非每个标签都有相同的离散度，这表明可以通过更复杂、更不激烈的稳定化技术来改进推断。对于微阵列数据，经验贝叶斯（EB）层次模型通过共享所有基因的结构（Smyth，2004). 这种战略是适应性的。如果方差没有很大差异，EB模型基本上是一个集合估计。然而，如果差异非常大，EB模型缩小的幅度较小。对于我们的NB模型，由于NB不属于指数族且φ不存在共轭先验，EB解受到阻碍。布拉德洛等。(2002)建议使用多项式近似，以避免随机马尔可夫链推理方法的计算开销。我们建议不要在所有标签上强制实施通用分散，而是挤压每个标记分散（表示为φ_克，用一个额外的下标表示标记）指向公共值（φ）。我们雇佣加权似然并选择似然权重以近似EB解。

公共似然在WL中的作用与φ的先验相同_克将在贝叶斯层次模型中使用，α为优先精度。如果α=0英寸(2)，然后我们得到标记式qCML估计。在另一个极端，如果α选择得足够大，则任何单个对数似然的贡献都会被共同的可能性所抵消，结果是共同的离散。在这两个极端之间有一个估计方案，其中标记估计介于个别估计和普通估计之间。

3.2选择α作为近似EB规则

我们希望选择一个适当的α，使估计具有自适应性。如果有证据表明分散度没有太大差异，则α的选择应足够高，以鼓励所有标签强烈收缩，接近一般估计值。然而，如果有证据表明分散度可变，则应选择收缩量较小的α。

我们无法直接使用选择α的规则，因为qCML估计量远离正态分布，没有已知方差，实际上可以在样本空间的边界上取零或无穷大的值，概率为正。为了避免这些困难，我们利用分数统计（对数似然导数）比最大似然估计更快地收敛到正态性这一事实。我们还注意到，估算公式(4)可以用可能性得分来表示S公司_克(φ) = ∂我_克（φ） /dφ和预期信息我_克(φ) =E类(J型_克),J型_克= − ∂²我_克(φ)/∂ φ²，φ的函数_克。这使我们可以将估算算法陈述如下。该算法与(三)和(4)但它具有更广泛的应用，因为它只使用在估计τ₀.

1. 找到常见的色散估计量最大化我_C.

2. 评估和对于每个标签。

3. 估计τ₀通过求解

   

   如果那么τ₀= 0.
4. 设置

   

5. 获得加权似然估计量通过最大化WL（φ_克).

预期信息我_克很难直接计算，但可以使用观测信息很好地近似J型_克.对于φ的任何给定值_克,我_克应该与总数几乎成正比z（z）₁+z（z）₂因此，我们计算截距为零的线性回归J型_克总伪计数（参见图1)，并使用拟合值表示我_克.

该算法实际上可以应用于φ的任何变换。我们发现在δ=φ/（φ+1）尺度上实现该算法很方便，因为δ取严格有界值。

3.3方法的解释

上述算法具有很好的统计解释。如果分散度真的相等（所有φ_克= φ₀)，然后以便将被估计为接近零，因此α将较大。然而，如果分散度确实不同，那么E类(S公司_克)将为非零，并且将大于我_克平均而言，强制大于零且权重较小是对常见可能性的赋值。分散度越不相同，越大将进行估算，并完成较少的收缩。事实上在零假设下是一个不依赖渐近正态性的精确结果。这确保了我们的算法即使在库数量很少的情况下也具有良好的定性行为。

4个结果

4.1压缩改进了负二项模型中离散度的估计

对于标签计数的NB模型，φ的估计是一个关键步骤，它会影响差异表达的确定重要性。请注意，φ对统计测试的直接影响不如微阵列中的方差，因为我们数据中的方差也是均值的函数。卢等。(2005)使用PL模型估算φ，并分别估算每个标签。在这里，我们计算了所有标签上的通用分散度，并在一种新的近似EB策略中缩小了向它的标签分散度。

我们首先表明，我们的近似EB方法改进了以均方误差（MSE）表示的分散度的总体估计。我们比较了四种估计策略：标记式qCML、使用近似EB规则的WL、通用qCML（Robinson和Smyth，2007)和标签式PL（Lu等。,2005)，超过三个真的情况。对于这里的所有比较，我们将库大小固定为50000，平均值为10（λ=0.0002），从NB中抽取1000个标签，并重复模拟50次。计算每个模拟的总MSE。MSE是根据标度，因为存在无限tag-wise qCML估计的非零概率。首先考虑的是中等数量的图书馆(n个=4）和固定色散，如所示图2答：这种情况显然有利于共同的分散估计，而共同的分散估算具有最低的MSE。但是，请注意，近似EB策略在这里也做得很好，为通用模型提供了足够的权重，以便将标签分散几乎完全压缩到通用值。其余两种情况涉及随机分散，这里是从伽马分布中提取的。我们选择伽马参数进行模拟，以匹配Zhang数据集（Zhang）上的经验分散估计分布（使用近似EB估计）等。,1997)，并使用了介质(n个=4）和大量库(n个= 10) (图2B和C）。我们可以看到，近似EB解决方案为估计离散度提供了显著优势，并适应了这种情况，显示了分数和信息中包含的证据的相关性。

重量很适合这种情况。在固定分散度的情况下，权重足够大，可以将估计值缩小到几乎相同的值。面对更分散的真实分散，重量如预期般减少。在随机分散和较大样本的情况下，权重再次降低，以自动调整，以获得标记估计中包含的更多信息，表明近似EB系统似乎正在实现预期目标，因此在实践中提供了合适的规则。

4.2方法比较：模拟数据

离散估计的改进可以提高我们从非DE中分离差异表达（DE）标签的能力等。(2005)并考虑进行更广泛、更现实的研究。

Lu中的模拟等。(2005)考虑在两种条件下采样10000个标签，λ固定₁，库大小在30000到9000之间均匀抽样，比较另一个库的1-5个库中的5个库。对于5000个标签，植入的差异为λ₂=b条.λ₁使用，其余标记没有差异（λ₂= λ₁). 当从NB取样时，他们选择0.17、0.42和0.95的固定分散度。我们重复了他们对λ的性能分析₁=0.0002和b条=4，与图2Lu的等。(2005). 在图3在本文中，我们使用Wald检验统计量比较了接收者操作特征（ROC）曲线，这两种统计量都用于PL估计和收缩估计。也就是说，我们拟合了Lu的广义线性模型等。(2005)用我们的φ估计来证明改进的估计是有益的。在所有情况下，改进的离散度估计器提高了分离真正的DE和非DE标记的能力。在这里，精确测试和Wald测试之间几乎没有区别。

接下来，我们以一些简单但重要的方式扩展了他们的模拟研究。首先，我们通过非固定分散度使问题更加现实。再次，我们将随机分散度设置为伽马近似经验估计分布（形状=0.85，尺度=0.5）。我们使用Zhang数据集λ估计的经验分布而不是固定的λ，并随机将伽马抽样分散度分配给λs。我们对模拟做出的一个微妙变化是，植入差异的乘数，b条，并不总是增加真正的意思是，因为计数越大，估计问题就越容易。相反，对于具有真正差异的采样标签，我们使用和按照真实比例。最后，我们考虑一个更现实的10%差异表达标签，并比较小(n个₁=n个₂=2）和中等(n个₁=n个₂=5）图书馆数量。

我们更喜欢错误发现（FD）图，因为它们突出了排名靠前的标签的性能，而不是ROC曲线。图4显示了四种情况下的FD图：中小真差(b条=4，8）和中小型图书馆。标签根据其测试统计（Wald）或P（P）-值（近似精确测试）。由于有1000个真正的DE标签，因此选择了前1000个(x个-轴），并绘制错误检测的数量(年-轴）。显然，首选较少的FD。我们发现，无论使用何种统计检验，使用近似EB估计进行离散度估计都会导致所有情况下的平均FD更少。如果我们要验证前200个标签，例如“Exact”的比率。EB到Wald。中所述四种情况的PL’FD分别为80/95、42/59、31/35和3/7图4从而显示出持续且实际有意义的改进。

当然，这只表明统计数据的顺序符合要求。由于精确测试不依赖于其分布假设的渐近性，因此最好能够实现一组假阳性率，允许设定合理的截止值，可能是在调整多次测试之后。在一项关于小样本NB假阳性率的研究中，Robinson和Smyth(2007)证明了精确检验最能达到标称假阳性率，事实上，Wald检验在所有渐近检验中的假阳性率最高。置换试验（Tusher等。,2001)使用如此少的库不太可能创建合理的null分布。

4.3 SAGE数据的应用

我们将该方法应用于来自Zhang的SAGE数据等。(1997)，因为有真正的生物复制品。我们只是比较了两个正常结肠文库和两个结肠肿瘤文库。图5A显示了分散度和比例（λ）估计值，假设库集合之间没有差异。在这里，我们可以看到，在低丰度下，存在少量分散性较大的标记。大多数情况下，两个复制计数中的一个为零，另一个为非零，在某些情况下，标签计数在一个条件下都为零，而在另一个条件中只有一个非零。在这些情况下，一次性qCML估计是φ=∞（δ=1），因此向公共估计压缩是必要的。注意，随着丰度的增加，分散度呈下降趋势并不奇怪，因为观测值的方差是平均值的函数。

将精确测试应用于正常结肠样本与结肠肿瘤的比较，我们发现肿瘤中有49个基因上调，115个基因下调，错误发现率为5%（使用Benjamini–Hochberg校正）。图5B显示了微阵列数据的“MA”图的类似图，其中x个-轴表示丰富年-轴表示两种条件之间的变化幅度。

5结论

NB数据的离散度估计对于评估平均值变化的重要性至关重要。对于复制量最小的标签计数数据，我们引入了一个加权条件似然估计量，该估计量将单个标签分散向共同分散。这个过程可以被认为是使用数据相关的先验值并找到最大后验估计值，或者简单地说是加权似然。我们根据近似经验贝叶斯规则选择收缩量。EB规则在实践中运行良好，根据样本大小和条件对数似然的一阶和二阶导数（分数和信息）中包含的证据，适应了离散度的真实相似性。估计弥散度的精度提高了实验条件之间测试的功率。对于测试两种实验条件之间的差异，精确测试的性能与Wald测试一样好或更好。精确测试还有一个额外的优点，即它可以达到接近其标称错误率的水平。

我们的加权似然收缩算法具有非常普遍的应用，只需要在公共参数估计下评估对数似然函数及其前两个导数。这种方法在许多其他基因组规模的估计和推断问题中可能被证明是有用的。

致谢

作者感谢Terry Speed的宝贵讨论。这项研究得到了NHMRC项目拨款406657的支持。

利益冲突：未声明。

参考文献

作者注释