快速计算样本量，同时控制错误发现率，并应用于微阵列分析

摘要

动机：样本量计算在实验设计中非常重要，在微阵列或蛋白质组实验中更是如此，因为只能重复几次。在涉及这些实验的多重测试问题中，控制假发现率（FDR）或正发现率（pFDR）比控制I类错误（如家族错误率（FWER））更有效、更合理。在控制FDR时，通过控制I型误差估计样本量的传统方法不再适用。

结果：我们提出的方法适用于控制FDR。样本量计算简单明了，需要最少的计算，如两个样本所示t吨-测试和F类-测试。基于对所得样本量的模拟，证明通过q值程序可以实现功率。

可利用性：如有要求，可提供实现所述方法的Matlab代码。

联系人： pliu@iastate.edu

补充信息：补充数据可在生物信息学在线。

1 简介

微阵列和蛋白质组实验在许多生物学科中变得流行和重要，例如神经科学（Mandel等.,2003)药物基因组学、遗传病和癌症诊断（Heller，2002). 这些实验在材料（样品、试剂、设备等）和实验室人力方面都相当昂贵。许多微阵列实验仅使用少量重复（2-8）（Yang和Speed，2003). 在许多情况下，样本量不足以实现可靠的统计推断，导致资源浪费。因此，科学家经常提出以下问题。样本量应该有多大？

为了回答这个问题，我们将计算控制一些错误率并达到所需功率的样本大小。在计算单个测试的样本量时，要控制的错误率通常是I型错误率，即通过拒绝真零假设得出假阳性结论的概率。然而，我们同时测试了大量的假设，每个假设都与一个基因有关。因此，多重测试通常应用于微阵列数据的分析。在这种情况下，有几种错误率需要控制，例如系列错误率（FWER）或错误发现率（FDR）。假设有米基因芯片上的基因和每个基因的差异表达意义进行测试。测试结果总结如下表1，其中，例如，V（V）是误报数R（右）是被拒绝的数量米测试（Benjamini和Hochberg，1995).

在许多基因组数据的情况下，如微阵列，Storey和Tibshirani认为(2003)更合理、更强大地控制FDR或pFDR，而不是FWER。然而，传统上，样本量是以特定的I类错误率计算的，不能直接应用于FDR控制。

几篇文章已经解决了微阵列实验中的样本量计算问题（Hwang等.,2002; Lee和Whitmore，2002; Warnes和Liu，2006.Lee和Whitmore(2002)在控制假阳性数时，使用方差分析模型计算样本量表(E类[V（V）]). 黄星京等(2002)提出了一种首先识别差异表达基因，然后在Fisher判别分析缩小的空间上计算功率和样本大小的方法。Warnes和Liu(2006)提出了一种基于累积图的方法来可视化功率和样本量之间的权衡。一些文章解决了不同设计中的样本量计算问题（Dobbin和Simon，2005)或特定设置，如分类（华等.,2005). 上述方法控制类型I错误，而不是FDR。

最近，有几篇文章研究了在控制FDR的同时计算样本量的必要性，并提出了实现这一目标的方法。杨等(2003)应用几个不等式得到与FDR控制水平相对应的I型错误率。由于应用的不平等，样本量可能被高估。帕维坦人等(2005)研究了几个工作特性曲线，以可视化FDR、灵敏度和样本量之间的关系。尽管他们的方法在计算样本量时很有用，但没有提供简单的直接算法。荣格（Jung）(2005)导出了FDR与I型错误率之间的关系式。然后，FDR由适当级别的I类错误率控制。英镑和程(2005)提出了一种迭代搜索样本大小的算法，在该样本大小下可以实现FDR的期望功率和受控水平。由于FDR控制程序在包括微阵列分析在内的许多问题的多项测试中越来越流行，因此在设计实验时能够计算控制FDR所需的样本量非常重要。

在这里，我们提出了一种在控制FDR的同时计算多次测试的样本量的程序。首先，对于非差异表达基因的比例和要控制的FDR水平的任何估计，我们为每个样本大小找到一个排斥区域。然后计算每个样本大小的所选拒绝区域的功率。根据期望功率，最终确定样本大小。

荣格的方法(2005)我们在完成初稿后就知道了这一点，与其他方法相比，它与我们提出的方法更为相关。Jung和我们的方法都基于相同的模型假设，导致相同的FDR表达式。然后通过研究FDR表达式与一个量的关系来控制它，这个量是Jung的I型错误率和我们的临界值（拒绝区域）Z轴-测试和t吨-测试。将我们的方法应用于Jung的设置时，也会产生相同的结果。然而，我们的方法比Jung的方法更具图形性。这使得功率和样本大小之间的权衡得以可视化，并在修改用户定义的数量（如功率）时提供快速答案。

尽管有相似之处，本文将该方法进一步扩展到几个不同的方向，我们发现我们的方法非常令人满意。首先，我们将我们的方法应用于F类-广泛用于微阵列数据分析的测试（崔等.,2005). 其次，当表达水平的均值和方差在基因之间发生变化时，我们仔细研究了我们的方法，这是微阵列的一个重要而实用的设置。第三，我们还通过仿真表明q个-Storey提出的控制FDR的价值程序等(2004)使用我们建议的样本量达到了令人满意的目标功率。这积极地回答了这样一个问题，即是否有任何统计程序可以实现该方法所声称的目标功率。最后，我们还将我们的方法与Yang进行了比较等(2003)英镑和程(2005)与其他文章相比，它提供了更多定义良好的算法。我们的仿真表明，我们提出的方法是优越的。

文章结构如下。第2节用两个示例描述了我们提出的方法t吨-测试和F类-测试。在第3节，我们报告了模拟研究的结果，将基于该方法的功率与q值程序的实际结果进行了比较。第4节总结了我们的结果。

Matlab中的建议方法代码可用于实现该方法。

2 方法

在本节中，我们首先说明我们的想法，然后展示如何将所提出的方法应用于两种微阵列实验设计。

2.1 建议的方法

2.2 所提方法的应用

微阵列实验通常是为了发现不同处理之间差异表达的基因。如Smyth中所述，微阵列的扫描强度数据通常经过质量控制、转换和标准化等(2003)和Quackenbush(2002). 我们假设数据在应用统计测试之前首先经过这些步骤。在实验之前，我们没有观察结果来检查分布。似乎有理由作出一个方便的假设，即预处理数据的分布是正常的，因此是两个样本t吨-测试和F类-测试是适用的。计算样本量的其他拟议方法也做出了相同的假设（Dobbin和Simon，2005; 荣格，2005; 华等.,2005; 黄星京等.,2002).

2.2.1 用t检验进行两样本比较

（7）的右侧在c（c）因此c（c）存在时是唯一的。同样的注释也适用于这两个方程式(14)和(17)在后面的章节中。有关证明，请参见附录C。在回答裁判的问题时，我们发现最小（超过c（c）)FDR水平为正，发生在c（c）→ ∞. 这很有趣，因为I型错误没有这样的正下限。然而，随着样本量的增加，最小FDR很快收敛到零。参见附录中的图S1。

图1当FDR控制在5%时，绘制功率与样本大小的关系图。例如，我们想确定当π₀= 90 %. 假设需要2倍的变化（相应地，Δ=log₂（2） =1），σ=0.5，则Δ/σ=2。在中使用中间曲线图1a、 80%的期望功率需要每组9个样本。

我们已经包括了π₀相对较小（50%）图1.当π₀因为微阵列的标准化方法通常依赖于大π的假设，所以微阵列数据应该谨慎地进行标准化₀即少量差异表达的基因。在这种情况下，我们建议使用家政基因来执行标准化。如果使用适当的σ估计（基于适当的标准化值），我们的方法仍然适用。

我们将σ取为0.2，这是Warnes和Liu中U133微阵列数据集标准偏差的中位数(2006)即健康志愿者的人类平滑肌细胞的基因表达水平。（一位裁判向我们提到，人体样本和U133A阵列的中值σ通常在0.7左右。在这种情况下，我们将σ设置为0.7。）等(2005)，0.2大约是颗粒细胞肿瘤微阵列数据的剩余标准偏差的第90个百分位。（这里第90个百分位是一个保守的选择，因为如果我们使用的百分比小于90%，所需的样本量就会更小。）如果仍然是2倍的变化（Δ=log₂（2） =1）被认为是真实的影响大小，则Δ/σ=5。从中间曲线图1b、 对应π₀=0.9，可以确定需要4个样本才能获得至少80%的功率。

2.2.2 与F检验的多样本比较

对于比较几种治疗方法的微阵列实验，应用了不同的设计方案（Yang和Speed，2003). 假设没有任何复制，设计需要秒幻灯片。我们称之为秒为这个设计幻灯片一套。例如，我们想比较三种独立治疗方法之间的基因表达，例如三种基因型小鼠的肝脏（Horton等.,2003). 如果我们应用如中所示的回路设计图2，双色微阵列实验需要一组三张幻灯片。无论复制品是不同的生物样品还是不同的技术重复，只要在计算中使用适当的参数（均值和方差），我们的方法都是适用的。我们建议在实验中使用不同的生物样品，因为这将提供更一般的结论。问题是有多少组幻灯片足以获得足够的功率和受控的FDR。

2.2.3 Δ不相等的情况_克和σ_克

到目前为止，我们一直认为所有基因都有相同的参数集。在这种情况下，所有基因的平均功率将与单个基因的功率相同。实际上，每个基因可能有一组不同的参数。如果我们以双样本比较为例，基因特异性参数包括σ_克、标准偏差和Δ_克治疗组和对照组的平均表达之间的真正差异。

对于以下方面的集成σ_克，我们可以将自适应Lobatto求积应用于数值积分，这样可以进行稳定的计算，以获得c（c）。此数值积分的计算可立即提供答案。一旦我们得到答案c（c）对于每个样本量，我们相应地计算幂，并基于幂找到所需的样本量。

三 模拟

本文中提出的计算样本量有多现实？更具体地说，如果所需功率为80%，FDR = 5%，我们的方法得出两个样本比较的样本量为9t吨-测试，是否有一个统计测试可以用9张幻灯片实现所有的操作特性？为了找出答案，我们用计算出的样本量模拟数据，并用FDR控制程序进行多次测试。然后我们检查：

• 多次测试是否实际产生了计算样本量所需的功率，以及

• 观察到的FDR是否与我们想要控制的水平相当。

如果我们能找到一个统计程序，在计算的样本量下达到所需的FDR和功率，那么我们的程序就被证明是实用的。事实确实如此。

有几个程序可以控制FDR，例如Storey和Tibshirani提出的q值程序(2003)和Storey等(2004)以及Benjamini和Hochberg提出的程序(1995，2000). 这些程序都由罗斯福保守控制（Storey等.,2004). 为了进行模拟研究，我们应用Storey中概述的q值程序等(2004)控制FDR。手稿的早期版本在Storey和Tibshirani中应用了该程序(2003)结果与这里的报告相似。

当观测值（基因）相互独立时，我们首先测试所提出的方法。在微阵列环境中，我们假设共有5000个基因，并且治疗组和对照组的样本大小相等(n个₁=n个₂=n个). 基因特异性方差，⁠，是根据逆伽马分布模拟的。与Wright和Simon中的相同(2003)，我们选择了1/σ²～Γ（3,1），因为这种分布很好地逼近了我们所分析的几个微阵列数据集。对于对照组，基因表达值由⁠对于治疗组，我们设置Δ_克=0（对于非差异表达基因），并模拟Δ_克从对于差异表达的基因，基因表达值由⁠.

模拟涉及几个参数，π₀（非差异表达基因的比例），σ_Δ（效应大小的标准偏差），对于相关情况，相关系数ρ。为了评估样本量计算方法的准确性，我们采用析因设计进行模拟，每个因子（参数）的水平（值）总结如下表2对于48个参数设置中的每一个，对于多次测试，FDR控制在5%。

对于独立案例的每个参数设置，我们计算每个样本大小的预期功率，并生成功率曲线，如所述第2节。我们还模拟了200组数据并执行t吨-用q值程序测试每个数据集（Storey等.,2004)控制FDR。对200个模拟数据集的观测功率进行平均，并记录观测到的错误发现比例。与仿真结果相比，基于我们计算的预期功率曲线与所有参数设置的仿真结果几乎没有区别。示例如所示图3a.因此，我们提出的方法提供了对样本大小的准确估计。观察到的FDR也接近控制水平（5%），如图3b、 证明Storey程序的有效性等(2004).

由于许多基因可能作为群体发挥作用，因此很可能在基因表达数据中存在依赖性。为了检查所提出的方法在违反独立性假设时的性能，还根据依赖性结构模拟了基因表达水平（Ibrahim等.,2002). 然后应用与上述相同的测试程序，并将所得功率曲线与我们的计算进行比较。

当Δ_克和所有基因都是一样的，仿真结果表明，我们的方法可以为独立基因和依赖数据提供准确的样本大小估计，与上述仿真结果类似。

有几篇文章讨论了在控制FDR的同时计算样本量的问题。在这些文章中，杨等(2003)英镑和程(2005)提供了明确定义的算法。我们在两个样本的背景下将我们的方法与这些方法进行了比较t吨-固定Δ试验_克和⁠.表3结果表明，基于本文方法计算的样本量与基于仿真结果所需的实际样本量相吻合。杨的方法得到了与我们相似的答案，只是在某些情况下，它有点保守。Pounds和Cheng算法的答案在一种情况下过于宽松（当Δ/σ= 1） 与其他两种方法相比，与正确答案的偏差更大。

4 讨论

微阵列实验中包含的阵列数量直接影响数据分析的能力。制定选择样本大小的指南至关重要。由于与这些数据集相关的巨大维度，在许多情况下，控制FWER是非常保守的（Storey和Tibshirani，2003). 相反，由Benjamini和Hochberg提出的FDR(1995)和Storey(2002)似乎是一种更适合控制的错误率，并已广泛应用于微阵列分析。因此，获得一种能够控制FDR并保证一定功率的样本大小的方法非常重要。

该方法很容易应用，如所述第2节对于t吨-和F类-测试。只要有一个明确的形式来计算单个测试的I型误差和功率，所提出的方法可以推广到其他测试。本文提出的方法允许在设计实验时以最小的工作量计算准确的样本量。

致谢

作者感谢两位评论员和格雷戈里·沃恩斯博士的深刻评论和建议。我们还感谢王冲博士指出数值积分的洛巴托求积。

参考文献

作者注释