跳到主要内容
访问密钥 NCBI主页 MyNCBI主页 主要内容 主导航
BMC生物信息学。2008; 9: 117.
2008年2月25日在线发布。 doi(操作界面):10.1186/1471-2105-9-117
预防性维修识别码:PMC2335282型
PMID:18298817

两样本比较微阵列实验中样本量估计的混合模型方法

摘要

背景

选择合适的样本大小是设计微阵列实验的一个重要步骤,最近提出了估计样本大小以控制错误发现率(FDR)的方法。其中许多方法需要了解差异表达基因之间的效应大小分布。如果可以确定此分布,则可以计算准确的样本量要求。

结果

我们提出了一种混合模型方法来估计两样本比较研究数据中的效应大小分布。具体来说,我们提出了一种新的闭合形式算法,用于估计差异表达基因的测试统计分布中的非中心参数。然后,我们展示了如何使用我们的模型来估计控制FDR的样本大小以及其他统计指标,如平均功率或假不发现率。通过与现有样本量估计方法的比较,评估了方法的性能,发现该方法非常好。

结论

提出了一种新的方法来估计两样本比较微阵列研究的适当样本大小。与现有方法相比,该方法表现得非常好。

背景

微阵列最常用的实验装置之一是双样本比较研究,即比较两种不同实验条件下样本的表达水平的研究。在重复两样本比较的情况下,可以使用统计检验来评估测量的差异表达的重要性。这样做的一个自然测试统计量是t统计量(参见示例[1]),这将是我们这里的重点。在两个样本比较的背景下,引入“效应大小”的概念也很方便。本文中的效应大小是指:基因平均表达水平的两个条件之间的差异除以表达水平测量的共同标准偏差。

在一个普通的微阵列实验中,同时测量数千个基因。对每个基因进行统计测试会导致多重假设测试问题,因此需要一种策略来控制测试中的假阳性数量。一种成功的方法是控制错误发现率(FDR)[2]或FDR变化,如阳性假发现率(pFDR)[].

为了从实验中获得想要的结果,重要的是要使用适当的样本大小,即生物复制次数。例如,可以根据指定的FDR和平均功率以及选择的样本大小来设定目标,以便实现目标[4].

在过去的几年中,提出了许多方法来帮助估计所需的样本量。一些早期方法[5,6]依靠模拟来观察样本大小对FDR的影响。后来的工作在样本量和FDR之间建立了明确的关系。最新方法的一个共同特点是,它们需要了解要运行的实验中效果大小的分布。在缺乏这种分布的情况下,有两种选择。第一种选择是简单地指定要使用的分布。这种选择可能与人们感兴趣的差异表达的特定模式相对应,也可能基于对效应大小如何分布的事先知识。许多可用的方法讨论了特定分布的样本大小估计[7-11]。第二种选择是从试验数据集估计所需的分布。费雷拉和兹温德曼[12]讨论一种这样的方法。假设测试统计量的概率密度函数是对称的,并且属于一个位置族,它们使用反褶积估计器获得想要的分布。应该注意的是,对于微阵列实验中经常使用的样本量,t统计量的非中心密度函数偏离了这些假设。. [13]英镑和程[14]讨论两种不同的方法。这两种方法都认识到差异调节基因的检验统计量是非中心分布的,目的是估计相应的非中心参数。从非中心参数可以找到效应大小。和我们一样,考虑t-统计量,并通过将三组分混合模型拟合到试验数据的观测统计量来估计非中心性参数。Pounds和Cheng考虑了F-统计量,并为每个观察值估计了一个非中心参数。然后,他们根据贝叶斯q值解释重新调整估计值[15]。最后一个需要提及的方法是Pawitan. [16]它使用基于似然的准则将混合模型拟合到观测到的t-统计量。Pawitan在论文中未将该方法作为样本量估计方法进行探讨.,但它可以用于此用途。

在本文中,我们介绍了一种混合模型方法,用于估计试验数据中观察到的t-统计量的非中心参数的潜在分布,从而也估算了影响大小。使用的混合成分的数量不受限制,我们提出了一种新颖的、封闭形式的模型参数估计算法。然后,我们演示了如何将此模型用于样本量估计。通过检查FDR和平均功率之间的关系,以及FDR和假不发现率(FNR)之间的关系。我们可以选择成对控制这些测量值的样本大小。为了验证我们的模型和样本量估计,我们在模拟数据上测试了它的性能。我们包括Hu方法的估计.、Pounds和Cheng以及Pawitan.

结果和讨论

符号、假设和测试统计

贯穿本文t吨ν(λ)表示t分布随机变量的概率密度函数(pdf)ν自由度和非中心参数λ.一个中心t pdf,t吨ν(0),也可以写入t吨ν.A型t吨ν(λ)评估时间:x个已写入t吨ν(x个;λ).

假设可以在初步研究中进行基因表达测量。对于特定基因,我们表示n个1条件1和n个2根据条件2X(X)1, (= 1, ...,n个1)和X(X)2j个, (j个= 1, ...,n个2). 让(μ1,σ12)和(μ2,σ22)是每个人的期望和方差X(X)1X(X)2j个分别为。为了简单起见,我们在本文中重点讨论以下情况σ12=σ22=σ2。正如微阵列数据分析中常见的那样X(X)1s和X(X)2j个假设s是正态分布的随机变量。

测量的表达水平通常在假设正态性之前进行转换。

在这种情况下,经常用来检测差分表达式的统计量是t-统计量。根据实验设置,可以使用两种版本的t-统计量。在第一个设置中,对每个条件分别进行测量。推断基于双样本统计

T型1=n个1n个2(n个1+n个2)1(X(X)¯1X(X)¯2)(n个1+n个22)1[(n个11)S公司12+(n个21)S公司22],
(1)

哪里X(X)¯k=n个k1=1n个kX(X)kS公司k2=(n个k1)1=1n个k(X(X)kX(X)¯k)2.在零假设下H(H)0:μ1=μ2,T型1有一个t吨n个1+n个22pdf格式。然而,如果,H(H)0不是真的,有区别μ1-μ2=ξ,统计数据的pdf是t吨n个1+n个22(ξσ1n个1n个2(n个1+n个2)1此设置包括比较单色阵列和双色阵列参考设计的测量值。在第二种实验中,测量值是成对的。例如,在以下情况下,n个直接比较两种条件的双色幻灯片,然后n个1=n个2=n个,使用的统计数据是

T型2=d日¯S公司d日2/n个,
(2)

哪里d日¯=n个1=1n个(X(X)1X(X)2)S公司d日2=(n个1)1=1n个(d日d日¯)2现在,在H(H)0:μ1=μ2,T型2具有作为t吨n个-1个pdf格式。如果μ1-μ2=ξ然而T型2是一个t吨n个- 1(ξ σ-1个n个).

在这两个实验装置中,我们注意到真正不受调控基因的t-统计量的pdf是t吨ν(0). 对于真正受调控的基因,pdf是t吨ν(δ),使用δ≠0反映基因的差异表达水平。我们还注意到δ与基因的效应大小成正比,ξ/σ. Theδs可以被认为是一些潜在随机变量Δ的实现,分布为小时(δ). 因此,根据我们的假设,观察到的t型芯应建模为t吨ν(δ)-分布,使用小时(δ)作为其混合分布。这个小时(δ)未直接观察到,必须进行估计。

在下文中,假设实验中计算的t-统计量是独立的。这个假设,以及这个假设σ12=σ22=σ2,在微阵列设置中可能不成立。在测试部分,我们将检查这些假设不满足的情况,以了解结果是如何受到影响的。

估算效果大小的算法

j个,j个= 1, ...,表示观察到的t统计量实验基因,具有Y(Y)j个s作为相应的随机变量。(f)()成为他们的pdf。我们的混合模型通常可以表述为

(f)(;小时) = ∫t吨ν(;δ)dh公司(δ),

哪里小时(δ)是任何概率测度,离散或连续。估计小时(δ)我们想找到一个概率测度,使可能性最大化,L(左)(小时),其中

L(左)(小时)==1(f)(;小时).

它已经显示出来了[17]为了解决这个最大化问题L(左)(小时)是有界的,只考虑离散概率测度就足够了或更少的支持点。受此激励,我们选择小时(δ)作为离散概率测度,并旨在拟合形式的混合模型

(f)()==0πt吨ν(;δ)=π0t吨ν(;0)+=1πt吨ν(;δ).
(3)

这个小时(δ)因此是一个分布,其中Pr(δ=δ)=π(= 0, ...,)、和=0π=1.第二种形式(f)()在(3)中是由于知道δ=0表示未调控基因。

我们现在的目标是找到一个具有固定分量数的模型(3)的参数+ 1. 很明显,找到这些参数可以用缺失数据问题来表示,这表明需要使用EM算法[18]。尽管在早期的工作中已经讨论过这种方法[13,16]到目前为止,还没有一种解决该问题的闭合形式EM算法。构造该算法的主要困难是缺乏非中心t-pdf的闭合形式表达式。在本节的其余部分中,我们将展示如何获得所需的算法。

与EM一样,随机分量标签向量Z轴1, ...,Z轴定义了Y(Y)1, ...,Y(Y)。这些有Z轴ij公司= (Z轴j个)根据是否等于一或零Y(Y)j个属于第个组件。一个Z轴j个根据

公共关系(Z轴j个=z(z)j个)=π0z(z)0j个π1z(z)1j个πz(z)j个,

其中z(z)j个是随机的已实现值Z轴j个.

我们首先认识到这样一个事实,即的非中心t-pdf本身就是一种混合物。变量的累积分布根据t吨ν(δ)是(请参见[19])

F类ν(;δ)=12ν2Γ(ν2)0v(v)ν1e(电子)v(v)2212πv(v)νe(电子)12(δ)2d日d日v(v).

差异化F类ν()关于,并替换v(v)=νu个,收益率

F类ν(;δ)=0u个2πe(电子)12(δu个)2u个2ν2ν21Γ(ν2)(νu个)ν1e(电子)νu个22d日u个.
(4)

这种形式的非中心t pdf可以被识别为正常的混合N个(δu个,1u个2)分布,按比例χν随机变量的混合分布U型基于(4)中的特征,我们可以引入一组新的缺失数据uij公司, (= 0, ...,;j个= 1, ...,)实现了U型ij公司,并定义为(Y(Y)j个|u个ij公司,z(z)ij公司=1)在a之后N个(δu个j个,1u个j个2)分配。以这种形式重述模型,作为混合物的混合物,是找到闭合形式算法的关键步骤。

这个j个s由z(z)ij公司s和u个ij公司s构成完整的数据集。可以编写完整的数据日志

日志L(左)c(c)(π,δ)==0j个=1z(z)j个日志(π(f)c(c)(j个|u个j个,z(z)j个=1)c(c)(u个j个)),
(5)

哪里π= (π1, ...,π),δ= (δ1。。。,δ)和(f)c(c)(|u个,z(z))和c(c)(u个)上述是正常的和有刻度的吗χν分布。EM算法的E步要求L(左)c(c)在(5)中,以数据为条件。结合(4)和(5),我们发现我们需要期望

E类(Z轴ij公司|j个)和E类(U型ij公司|j个,z(z)ij公司= 1).

计算第一个期望值很简单(参见示例[20])被发现是

E类(Z轴j个|j个)=πt吨ν(j个;δ)=0πt吨ν(j个;δ).
(6)

计算第二个期望值比较困难,但通过使用贝叶斯定理,我们发现可以将其表述为(为了清楚起见,现在省略了指数)

E类(U型j个|j个,z(z)j个=1)=0u个(f)c(c)(u个|,z(z))d日u个
(7)

=0u个(f)c(c)(|u个,z(z))c(c)(u个)0(f)c(c)(|u个,z(z))c(c)(u个)d日u个d日u个
(8)

哪里(f)c(c)(u个|,z(z))是的pdf吗(U型ij公司|j个,z(z)ij公司= 1). 请注意c(c)(u个|z(z))=c(c)(u个)自U型ij公司Z轴ij公司都是独立的。

现在必须计算(8)中的积分。为此,我们注意到被积函数的分母本身是一个积分,它不依赖于积分变量u个实际上,(8)是两个积分的比值。在两个积分中替换变量后,(w个=u个2+ν)(w个=u个2+ν),我们发现这个比率可以改写为小时-功能为

E类(U型j个|j个,z(z)j个=1)=ν+1j个2+νH(H)小时v(v)+1(j个δj个2+v(v))H(H)小时v(v)(j个δj个2+v(v)),
(9)

哪里H(H)小时k(x个)=0w个kk!e(电子)12(w个+x个)2数据仓库对于整数k≥ 0. 的属性小时k-功能在中进行了讨论[21]。一个特别好的性质是它满足递归关系

(k+ 1)小时k+ 1(x个) = -x小时k(x个) +小时k- 1(x个).

易于计算H(H)小时1(x个)=e(电子)12x个2,H(H)小时0(x个)=x个e(电子)12u个2d日u个,我们有一种方便的计算方法(9)。

M步骤要求最大化L(左)c(c)关于πδ这是通过简单的微分和产量最大化来实现的

π^=1j个=1z(z)j个
(10)

δ^=j个=1z(z)j个j个u个j个j个=1z(z)j个.
(11)

等式(6)和(9)–(11)构成了所需闭合形式EM算法的主干,以拟合类似于(3)的混合模型+1个组件。根据方案更新参数估计

.)单位j个(k)=E类π^(k)δ^(k)(U型j个|j个,z(z)j个=1),z(z)j个(k)=E类π^(k)δ^(k)(Z轴j个|j个).)π^(k+1)=1j个=1z(z)j个(k),δ^(k)=j个=1z(z)j个(k)j个u个j个(k)j个=1z(z)j个(k).

在收敛时,估计πδ用作小时(δ)带有+1点质量。如上所述,我们修复了δ0= 0.

一个备受关注的问题是,当许多假设测试同时进行时,估计真零假设的比例(例如[,22,23]). 在微阵列环境中,这相当于估计所有考虑的基因中真正不受调控的基因的比例。参考(3),我们看到这个数量作为π0.利用π0-估计,我们建议使用已知的保守值π0-估计以指导一些模型参数的选择。下文将对此进行讨论。在我们的实现中,我们使用在[24],但可能会输入不同的估计值。

评估混合物模型中适当数量的组分是一个尚未完全解决的难题。经常采用的策略是使用Akaike信息标准(AIC)等措施[25]或贝叶斯信息标准(BIC)[26]。我们从模拟中发现,在我们的设置中,这些标准似乎低估了所需的组件数量(有关一些模拟结果,请参阅测试部分)。这可能是由于在微阵列数据中经常发现大量未调控基因。由于调控基因相对较少,根据这些标准,拟合额外成分的可能性增加不足以证明所使用的额外参数是合理的。在我们的实现中,我们使用=对数2((1 -π^0)),其中π^0是上述估算。这种选择是由经验驱动的,但在我们的数值研究中已经证明了它是足够的。它还反映了这样一个事实,即一个单一的成分应该为未调控的基因提供充分的解释,而其余的成分组件解释了受监管的组件。一个不同的可由样本量方法的用户指定。

拟合混合模型时可能出现的一个复杂情况是{δ}=1,可以分配给δ=0,或非常接近,从而影响中心构件的配合。为了避免这种情况,我们在附近定义了一个小邻里δ0=0,其中所有其余组件被排除在外。一个δ,≠0,即在拟合模型时尝试穿过此邻域,但仅在邻域边界处停止。通过找到最小值来确定边界δ˜因此可以区分t吨ν(0)来自t吨ν(δ˜)一个,基于尺寸样本π^0和(1-π^0)米/克分别为。后一个样本量假设受调控基因在其组分中均匀分布。样本作为两个分布的0.05和0.95分位数内的均匀间隔点。用于检查两个样本是否来自同一分布的标准是两样本Kolmogorov-Smirnov检验,显著性水平为0.05。这个标准背后的理由是,对于与调控基因相关的组分,我们只想让那些具有合理确定性的组分与中心组分区分开来。再一次π^0使用的是上面讨论的估计值。

与拟合混合模型相关的另一个困难是,最佳拟合不是唯一的,这可能会导致收敛问题。困难是由于置换混合物成分不会改变似然函数这一事实。在我们的实施中,我们没有任何约束来解决这个问题。然而,我们在一些测试运行中跟踪了我们的混合物成分的更新,并没有发现出现此问题。

总之,我们的方法提供了估计,{δ^}=0{π^}=0数据中的非中心参数和一组权重。这些数量一起构成了对分布的估计小时(δ). 如测试统计a部分所示δ与效果大小成正比,ξ/σ。因此,不对数据中平均值或方差的数值大小进行估计或假设。我们只估计了一组均值-方差比。

FDR控制中样本量估计算法

实验设计中的一个重要问题是估计实验成功所需的样本量。我们现在概述了如何选择控制FDR的样本大小以及其他度量,以及如何将上述模型用于此目的。

表11总结了假设检验。除之外的所有表元素是随机变量。从表表11常用的阳性假发现率(pFDR)[]可以定义为E(R(右)0/M(M)R(右)|M(M)R(右)> 0). 在[]也证明了,对于给定的显著性区域,在假设测试独立性的情况下,我们有

表1

的结果假设检验

H(H)0已接受H(H)0被拒绝总计
H(H)0真的一个0R(右)00
H(H)0一个1R(右)11
总计M(M)一个M(M)R(右)
pFDR=Pr(H0真|H0拒绝),
(12)

就p值而言第页和选定的p值截止α,(12)可以重写为

pFDR公司(α)=π0α公共关系(第页<α).
(13)

方程(13)在样本量估计中很重要,因为它提供了pFDR、显著性区域和样本量之间的关系。样本大小将决定Pr(第页<α)因为t分布的形状取决于n个1,n个2(本节末尾提供了解释)。Hu首先讨论了(13)在样本量估计中的应用. [13].

使用(13)和拟合混合模型,可以估计达到指定值的样本量α和pFDR。剩下的问题是选择合适的α和pFDR。pFDR是一种易于理解的度量,其大小可以由样本大小估计方法的用户直接设置。如何挑选α另一方面,情况并不那么清楚。一个解决方案是将(13)重申为αpFDR和其他统计指标。.提出一种这样的关系。他们建议选择α通过指定被拒绝的假设的预期数量,E(M(M)R(右)). 他们的想法是取代Pr(第页<α)=E(M(M)R(右))/在(13)中获得

α=pFDR公司π0E类(M(M)R(右)).
(14)

用文字表示,而不是指定α,可以指定E(M(M)R(右)). 这种获取方式α然而,它有一个缺点。它几乎没有向用户提供有关实验识别受调控基因能力的直接信息。在我们看来,选择一种更具信息性的方式α将允许用户指定数量,例如平均功率或假不发现率(FNR)。我们现在讨论如何实现这一点。

权力被定义为拒绝错误假设的概率。在微阵列多重假设设置中,平均功率控制正确识别的受调控基因的比例。设置α通过一种直观而有吸引力的衡量标准,平均功率将因此而有所帮助。之间的关系α平均功率可以通过将(13)中右侧的分母改写为

优先级(第页<α)=优先级(第页<α|小时0true)π0+优先级(第页<α|H(H)0假)(1-π0).

识别Pr(第页<α|H(H)0false)作为平均功率,(13)可以倒置求出

α=pFDR公司1pFDR公司1π0π0平均功率.
(15)

结合(15)和(13)可以找到达到指定pFDR和平均功率的样本大小,而无需指定α.

另一个有趣的控制指标是虚假不被发现率(FNR),即假设中虚假否定的预期比例被拒绝。换言之,FNR控制被错误接受为不受调控的受调控基因的比例。我们使用中讨论的FNR版本[15]称为pFNR=E(一个1/M(M)一个|M(M)一个>0),在与pFDR相同的假设下,可以表示为

pFNR=Pr(H0假|H0接受)。

根据pFDR和α产量

α=1pFDR公司π01pFNR公司pFDR公司.
(16)

同样,指定pFNR将对应于特定的选择α.pFNR设置方法α使用起来可能很有趣。然而,应该注意的是,在微阵列设置中,这种测量有时很难应用。原因是潜在的M(M)一个由于非调节基因的比例很高。A大型M(M)一个使pFNR数值较小,并且可能很难设置合理的大小。

选择pFDR和α(根据平均功率或pFNR),我们需要找到解决(13)问题的样本量。为此,我们表示Pr(第页<α)通过我们的混合模型(3)as

公共关系(第页<α)=>|0|(f)(;小时)d日,
(17)

哪里0是对应于p值截止值的t核α用于检验零假设。(f)()是以下各项的加权和t吨ν(δ)s、 需要取的积分只是(非中心)t分布的分位数。样本大小影响(17),因为t吨ν(δ)表中列出了要集成的t吨1(n个)(ξσ12(n个)),使用1(n个)和2(n个)依赖于n个1n个2(参考(1)和(2)的pdf)。如果效果大小{ξσ1}=0和重量{π}=0我们知道我们可以为所有人进行评估1(n个),2(n个).

然而,这些参数可以从拟合模型中找到。为了获得效果大小,我们可以将其相等δ=ξσ12(n个˜),=(0,...,),其中ñ是拟合中使用的样本大小。直接估算重量。表达了Pr(第页<α)就样本量而言,我们的目标是找到一个能够解决(13)问题的方法。这个问题可以重新表述为查找函数的根

S公司(n个)=pFDR公司=0(π>|0|t吨1(n个)(ξσ12(n个))d日)π0α.

为了找到的根S公司(n个)我们实现了一种二分法。

测试

为了评估我们的方法,我们实现了EM算法和上面描述的样本大小估计方法。然后我们对模拟数据集进行测试。由于上述原因,我们在样本量估计中重点控制了平均功率和pFDR。

当使用模拟数据集时,可以计算实现pFDR和平均功率给定组合所需的真实样本大小。为了评估我们方法的性能,我们将估计值与真实值进行了比较。为了与现有方法进行比较,我们还包括了Hu的方法所做的估计. [13]、磅和成[14]和帕伊坦. [16].

使用现有方法

在他们的论文中Hu讨论三种不同的混合模型。在我们的比较中,我们使用了他们的截断正态模型,因为这似乎是最受欢迎的模型。要使用此模型生成样本量估计值,需要输入想要的pFDR和E(M(M)R(右)). 由于我们希望控制pFDR和平均功率,我们计算了E(M(M)R(右))对应于pFDR的每个选择,平均功率为E(M(M)R(右))=平均功率·(1 -π0)/(1–pFDR)。所有测试都使用源代码中的默认参数运行。在胡锦涛的实施中但是,有三个参数(diffthres、iternum、gridsize)缺少默认值。作者提供了合理的选择(0.01,10100)。

对于Pounds and Cheng方法,需要输入称为预期错误发现率(aFDR)的数量和预期平均功率。对于这里给出的估计,我们使用了相应的pFDR和平均功率组合。由于Pounds和Cheng使用F-统计量,因此在我们的测试中计算的t-统计量也进行了相应的转换。(如果T跟随t吨ν分布,然后T2作为分发F类1,ν). 在实施过程中,Pounds和Cheng设定了样本大小估计值的上限nmax,并用nmax本身替换所有高于nmax的估计值。在我们的测试中,默认值nmax=50被替换为nmax=1000。原因是,在测试中,样本量估计可能超过50,而且确实如此。使用低nmax将影响与其他方法的比较。除了nmax之外,所有测试都使用源代码中的默认参数运行。

在[16]帕维坦人讨论一种将混合模型拟合到t-统计量的方法。拟合模型用于估计未调节基因的比例,π0。提到了使用他们的方法估计样本量,但没有进一步探索或测试。拟合模型所需的唯一输入是混合组分的数量。在我们的测试中,这个数字是根据Pawitan的建议使用AIC确定的在他们的论文中。在执行Pawitan.靠近中心部件的非中心部件的分配不受限制。使用未经调整的模型拟合进行的初步测试表明,这大大恶化了样本量估计。因此,在我们的测试中,我们通过在给定阈值内折叠所有非中心组件来调整其拟合模型(|δ|<0.75)进入中央。我们的模型调整与π0-帕维坦实施中使用的估算程序。对于测试运行,我们使用我们的程序根据此方法拟合的模型生成样本量估计值。

测试程序和结果

对于此处显示的测试结果,我们使用=10000个基因,以及n个1=n个2=每组5次测量。对于未调节基因的比例,我们检查了π0=0.7和π0= 0.9.

在第一组测试中,我们考虑了正态分布测量和等方差情况下的样本量估计。在这种情况下,我们模拟了有关联结构和无关联结构的数据。非中心参数的真实分布1= (1 -π0)受调控基因的生成方式如下:从N个(1,0.52)和aN个(-1,0.52)分配。两个样本的大小都相同1/2,他们一起组成了{ξk}k=11用于数据集中的调控基因。测量差异设置为σ2=σ12=σ22=0.52。然后将受调控基因的非中心参数计算为{ξkσ1n个1n个2(n个1+n个2)1}k=11(请参阅(1)的讨论)。根据我们的微阵列数据分析经验,上述选择似乎适用于对数转换微阵列测量。由于真实的非中心参数都是已知的,因此可以计算出实现特定pFDR和平均功率所需的真实样本大小。

相关性是使用块大小为50的块相关结构引入的。这种结构背后的原因在[27]。所有的基因,不管是受调控的还是不受调控的,都被随机分配到它们的区块。然后,通过从统一的U型(-1,1)分布到对称矩阵的非对角元素中。然后,使用[28],我们迭代找到与随机生成的矩阵最接近的单位对角线半正定矩阵。

使用上述方法,我们模拟了数据。对于每个测试用例,我们生成了50个数据集,并基于这些数据集进行了样本大小估计。在最初的测试中,我们想评估我们选择使用更多的混合物成分,,而不是AIC或BIC标准所建议的。为此,使用我们的算法将两个模型拟合到每个模拟数据集。一个模型有我们选择的混合组分数量,另一个模型具有AIC指示的数量。然后对这两个模型进行了样本量估计。与AIC而不是BIC进行比较的原因是,在这种情况下,BIC所支持的组件比AIC更少。在这个初始运行中,只使用了不相关的数据。样本量估计结果如表所示表2。2AIC和我们的方法选择的平均成分数量分别为6.1和11.8π0=0.7、5.0和10.1π0= 0.9. 根据我们的发现,我们得出结论,使用比AIC建议的更多的组件可能会有好处,并且我们在剩余的测试中使用了更多的组件。然后,我们将注意力转向比较不同的样本量估算方法。使用所有四种方法生成不相关和相关数据,并生成样本大小估计值。结果见表的上半部分表3。一般来说,我们的估计值似乎接近其真实值。当基因之间没有相关性时,结果稍好。正如预期的那样,随着功率的增加,精度降低,标准偏差增加。这可能与描述非中心参数在无调节点附近的分布的困难问题有关,即接近δ= 0. 胡的估计与我们的估计值相比,似乎在很大程度上偏离了真实值,并且更加保守,但标准差更低。在高功率估计中,与真实值的偏差特别大。Pounds和Cheng的估计值似乎偏离了真实值,比我们的估计值更保守,标准偏差更高。Pounds和Cheng的估计值的保守性也可以从他们自己的数值测试中看出,其中估计的实际功率超过了期望功率。帕维坦的估计看起来比胡强.和Pounds和Cheng,但似乎仍然比我们的估计值离真实值更远,并且有更高的标准偏差。对于高功率,使用此方法有低估所需样本大小的趋势。

表2

评估混合物成分的数量。

π0pFDR公司权力真的JMB(标准偏差)AIC(标准偏差)
0.70.050.666 (0.3)7 (0.5)
0.70.050.788 (0.5)8 (0.8)
0.70.050.81110 (1.1)15 (4.0)
0.70.050.92422 (2.5)52 (22.5)
0.70.010.699 (0.7)9(0.5)
0.70.010.71111 (0.8)12 (1.0)
0.70.010.81615 (1.8)25 (8.4)
0.70.010.93535 (3.6)79 (34.4)
0.90.050.698 (1.1)11 (3.8)
0.90.050.71111 (2.4)15(5.6)
0.90.050.81615 (4.1)21 (8.5)
0.90.050.93524 (7.9)33 (13.7)
0.90.010.61111 (1.9)15 (6.3)
0.90.010.71414 (3.2)21 (8.8)
0.90.010.82120 (6.3)29 (12.7)
0.90.010.94532 (11.1)44 (21.8)

模拟数据集的真实和估计的每组样本大小π0=0.7和π0=0.9,对于不同的pFDR和平均功率截止值。报告的样本量估计值是50个此类估计值的平均值,四舍五入为最接近的整数。标准偏差(sd)基于相应的50个数据集。对于每个数据集,本文介绍的估计方法都有两种不同的选择,混合物成分的数量。JMB列(来自作者姓名)使用如本文所述。AIC列列出了使用AIC标准进行选择的结果.

表3

评估不同方法的样本量估计值。

无相关性具有相关性
π0pFDR公司权力真的JMB(标准偏差)HZW(标准差)个人电脑(sd)PMMP(标准偏差)真的JMB(标准偏差)HZW(标准偏差)个人电脑(sd)PMMP(标准偏差)
0.70.050.666 (0.4)11 (0.0)11 (1.6)6(0.4)66 (0.5)11 (0.0)11 (1.4)6 (0.5)
0.70.050.788 (0.6)18 (0.2)14 (2.8)7 (1.0)88 (1.0)18 (0.0)14 (2.3)7 (0.9)
0.70.050.81111(1.5)39 (0.7)20 (4.8)9 (2.5)1111 (2.0)38 (0.5)19 (4.1)9 (2.0)
0.70.050.92324 (3.4)146 (2.4)30 (9.5)13 (6.5)2323 (4.5)145 (1.6)29 (7.9)15(6.8)
0.70.010.699 (0.5)17 (0.5)16 (2.7)9 (0.9)99 (0.7)17 (0.5)16 (2.4)8 (0.8)
0.70.010.71111 (1.1)28(0.5)21 (4.6)10 (1.7)1111 (1.7)28 (0.5)21 (3.8)10 (1.7)
0.70.010.81616 (2.4)60 (1.1)28 (7.9)13 (4.5)1616 (3.3)60 (0.8)27 (6.4)13(3.4)
0.70.010.93437 (6.0)231 (4.2)42 (14.5)18 (10.0)3236 (7.5)229 (2.8)40 (11.9)22 (11.3)
0.90.050.699 (2.1)16 (0.5)24 (7.8)9 (3.9)89 (2.6)16 (0.5)23(6.6)9 (3.1)
0.90.050.71111 (3.5)27 (0.8)31 (11.4)11 (7.1)1012 (4.1)27 (0.8)30 (9.7)11 (6.6)
0.90.050.81616(5.2)59 (1.7)41 (16.7)15 (11.1)1416 (5.9)58 (1.7)41 (14.4)15 (11.8)
0.90.050.93426 (7.6)227(6.6)60 (26.8)21 (17.9)2925 (8.5)225(6.7)59 (23.1)22 (18.8)
0.90.010.61112 (3.3)22 (0.7)33 (11.8)12 (6.0)1112 (4.1)22 (0.6)32 (9.9)12 (4.9)
0.90.010.71415 (5.3)38 (1.2)43 (16.9)15 (10.4)1316(6.0)37 (1.2)42 (14.4)15 (9.9)
0.90.010.82121 (7.4)82 (2.6)56 (24.3)20 (15.6)1921 (8.4)81 (2.7)55(20.9)21 (16.8)
0.90.010.94635 (10.0)318 (10.0)79 (37.2)27 (24.3)3833 (11.3)316 (11.3)78 (32.0)29 (25.1)
0.70.050.666 (0.2)6 (0.0)10 (1.0)6(1.1)66 (0.3)6 (0.0)10 (1.1)5 (0.7)
0.70.050.778 (0.7)8 (0.3)12 (1.7)8 (2.4)78 (0.7)8 (0.1)12 (1.8)7 (1.6)
0.70.050.8911 (1.4)16 (0.4)16(2.9)10 (4.9)911 (1.5)16 (0.2)16 (3.2)10 (3.8)
0.70.050.91423 (4.1)56 (1.2)24 (5.5)18 (11.3)1524(5.3)56 (1.0)25 (6.1)14 (7.8)
0.70.010.688 (0.6)8 (0.0)11 (1.7)8 (2.2)88 (0.7)8 (0.0)14 (1.8)8 (1.0)
0.70.010.71011 (1.1)12 0.1)14 (2.8)11 (4.3)1011 (1.3)12 (0.1)18 (2.8)10 (3.1)
0.70.010.81316 (2.5)23 (0.5)24 (4.6)15 (8.4)1516 (2.5)23 (0.3)24 (5.0)13 (6.4)
0.70.010.92636(7.2)84 (1.8)35 (8.4)27 (17.8)3037 (8.8)83 (1.3)34 (9.4)20 (12.0)
0.90.050.6810 (2.0)7 (0.4)21(8.2)8 (1.6)810 (2.3)8 (0.4)24 (8.5)9 (3.6)
0.90.050.7913 (3.3)11 (0.7)27 (12.1)10 (3.9)913 (3.4)12 (0.7)32 (12.8)11(5.7)
0.90.050.81218 (5.5)21 (1.2)37 (19.0)13 (8.0)1319 (5.0)23 (1.5)44 (19.7)14 (8.7)
0.90.050.92431 (8.5)76 (4.6)54 (30.2)18 (14.5)2531 (7.9)85 5.4)65(32.6)20 (14.5)
0.90.010.61113 (3.1)9 (0.5)29 (12.6)11 (2.4)1113 (3.6)10 (0.6)34 (13.2)12 (5.6)
0.90.010.71317(5.0)16 (0.6)38 (18.5)14 (6.2)1419 (5.0)17 (0.8)45 (19.6)15 (8.2)
0.90.010.81825 (8.1)28 (1.5)50 (27.2)17 (11.8)1926 (7.1)31(1.9)60 (29.2)19 (12.3)
0.90.010.95143 (11.4)103 (6.2)72 (42.8)23 (19.8)5343 (10.9)115 (7.9)88 (46.3)26 (19.7)

具有π0=0.7和π0=0.9,对于不同的pFDR和平均功率截止值。报告的样本量估计值是50个此类估计值的平均值,四舍五入为最接近的整数。标准偏差(sd)基于相应的50个数据集。使用本文讨论的方法进行的估算在表中称为JMB(来自作者姓名),而使用胡讨论的方法所做的估算[13],Pounds and Cheng[14]和Pawitan[16]分别称为HZW、PC和PMMP。

还运行了具有正态分布测量值的测试,其中{ξk}k=11从伽马分布中得出,根据下面讨论的模型,方差有所不同。结果与上述结果相似(未显示)。

在第二组测试中,我们希望模拟具有真正微阵列实验特征的模型中的数据。我们还想了解如果正态性和方差相等的假设不成立,样本量估计会受到什么样的影响。为了实现这一点,我们基于包含在limma软件包中的涡流数据集进行模拟[29],以及Rocke和Durbin讨论的基因表达水平测量模型[30]。Roke和Durbin的模型表明

w个=α+μe(微秒)η+,
(18)

哪里w个是强度测量值,μ是表达式级别,α是背景级别η错误项分布为N个(0,σε2)和N个(0,ση2)分别是。使用他们在论文中讨论的估计方法,我们估计了涡流数据集(18)的参数。估计参数(α^,σ^ε,σ^η)突变型为(394.12,150.84,0.18),野生型为(612.99,291.40,0.19)。为了生成一组代表相同数据的对数比率,我们进行了limma用户指南中概述的显著性分析。对于FDR调节的p值,使用0.10的截止水平,我们获得了一组可能被调节的基因的280个对数比。在我们的测试中,受调控基因的对数比率是从这组样本中进行替换的。真正的表达水平是通过从突变数据集中的基因的背景校正平均强度中取样而产生的。为了模拟两种条件下的微阵列数据,使用了以下程序:对一组对数比率和一种条件下真实表达水平进行采样。然后计算其他条件的真实表达水平。使用上述模型及其各自的参数集,对两种条件下的测量进行了模拟,然后进行对数转换。为了在这个设置中引入相关性,我们添加了一个随机效果,γ到每个相关基因块的对数转换测量。这个γ是从N个(0,0.5ση2)分配。再次假设区块大小为50,并将基因随机分配到每个区块。在这种情况下,通过反复从给定模型中提取数据集,并在t-统计量的截止值的精细网格上计算其平均功率和FDR,来估计真实的样本量要求。由于调控基因已知,因此可以进行直接计算。

在生成如上所述的模型后,我们再次模拟了有关联结构和无关联结构的数据。对于每个测试设置,我们采样了50个数据集,并使用所有四种方法从数据中估计了样本量。结果汇总在表的下半部分表3。.对于胡的方法.和Pounds and Cheng的趋势与第一组测试中的趋势相同。我们的方法似乎稍微高估了所需的样本量,而Pawitan方法现在有趣的是,提供了最接近真实值的估计值。帕维坦估计值的标准偏差仍然比我们的高一些。

注意,自从Hu的实现.和Pounds and Cheng仅支持基于双样本t统计(1)的样本大小估计,列出的所有测试都基于此统计。我们的实现支持这两种类型,并且还运行了测试来检查单样本t统计量的情况(2)。结果与双样本t统计量的结果相似(未显示)。

结论

在这篇文章中,我们讨论了一种混合模型方法,用于估计两样本比较微阵列研究中受调控基因之间的非中心参数分布,从而估计效应大小。该模型可以拟合从试验数据中计算出的用于研究的t-统计量。我们还说明了如何使用该模型来估计控制pFDR的样本大小以及其他统计指标,如平均功率或pFNR。在微阵列设置中,当使用FDR和FNR而不是pFDR和pFNR时,我们的结果通常也是近似有效的。这是由于,参考表表1,1,经常会有Pr(M(M)R(右))≈1和Pr(M(M)一个)≈1。在任何大规模微阵列实验的规划中,所提出的样本大小估计方法都是有用的。

我们通过进行一系列数值研究来检验样本量估计的准确性。结论是,对于所用误差度量的适度截止值,我们的估计是相当准确的,并且方差很低。对于严格的截止值,我们看到较大的方差和较低的精度。总的来说,我们的方法似乎比胡的可用样本量估计方法提供了更好的结果.和Pounds and Cheng。我们还评估了Pawitan的方法用于拟合混合模型及其在样本量估计中的使用。使用该方法的结果很好,我们的测试表明优化Pawitan方法用于样本量估计可能会很有趣。

估计中发现的严格截止值的准确性降低,可能是因为很难描述受调控基因在接近无调控点时的分布,也就是接近无调控的点δ= 0. 重要的是,这些基因的特征描述是准确的,但它不会影响未调控基因的估计分布。我们在本文中的解决方案是在δ=0,其中未安装其他部件。区分接近0的两个分布的更好方法可能是进一步研究的主题。

在我们的测试中,我们还检查了基因之间的相关性如何影响样本大小估计。我们发现,这些估计仅受到轻微影响。尽管如此,我们认为,结合基因之间相关性的估计方法是未来研究的一个重要课题。

对于微阵列设置,使用调节t统计量,如[31]通常比普通的t统计量更合适。为了使我们的方法直接适用于这类统计,我们需要知道非调节基因的调节t-统计遵循中心t-分布,并且需要分布的自由度。在[31]这个分布结果是成立的,中心t分布的自由度增加了。人们还需要知道,在一定程度的调控下,受调控基因的调节t-统计量遵循非中心t-分布,并且需要分布的自由度和非中心参数。对于第二个结果,我们没有发现任何结果。如果第二个结果被证明成立,那么我们的方法也可以应用于慢化t统计量。即使没有第二个结果,用适度的t-统计量测试我们算法的性能也可能很有趣。

在我们的工作中,我们将重点放在平均功率上,作为与pFDR一起使用的控制措施。应该注意的是,平均功率为0.9意味着能够正确识别90%的受调控基因。在许多实验中,实现这一点可能并不有趣。一个例子是,研究的目的仅仅是确定最能区分样本与其他样本的一小组标记基因。其他例子是,在比较处理过的和未处理过的组织样本时,只对受影响最严重的调节途径感兴趣的实验。在这两个例子中,远低于0.9的功率就足够了。在这种情况下,我们的估计特别有用,因为对于pFDR较低的中等功率切断,它们似乎都是准确的,并且方差较小。

尽管我们在本文中的主要目标是拟合一个用于样本量估计的模型,但我们要强调的是,拟合模型本身确实提供了一些有趣的信息。一个例子是对π0,未调控基因的比例,我们经常发现它比现有方法制作的更好。

未来工作的其他主题包括加速算法。众所周知,简单EM算法的收敛速度相当慢,将实施改进方法。另一个可能需要进一步研究的问题是选择模型组件的数量。如前所述,使用比传统使用的信息标准建议的更多的成分通常会大大提高样本量估计。在此设置中,可能需要使用不同的方法来选择组件的数量。

本文讨论的方法适用于使用t检验的任何两样本比较多假设检验情况,而不仅仅适用于微阵列设置中的问题。

可利用性

作者可根据要求提供本文中讨论的R环境方法的实现。

作者的贡献

HM和TSJ开发了统计算法,AMB提供了生物学见解。TSJ实施了这些方法,进行了测试并起草了手稿。HM和AMB修订了手稿。所有作者阅读并批准了最终手稿。

致谢

这项工作得到了挪威研究委员会(AMB,TSJ)生物技术和功能基因组学(FUGE)项目NFR 143250/140和NFR 151991/S10的资助。挪威科技大学(HM)的跨学科项目“BIOEMIT–功能基因组学的预测和修改:结合生物信息学、生物伦理学、生物医学和生物技术研究”也提供了资金支持。作者感谢Mette Langaas阅读手稿并提出改进建议。我们还要感谢两位匿名评论员的仔细阅读和有益的建议。

工具书类

  • Callow MJ,Dudoit S,Gong EL,Speed TP,Rubin EM。微阵列表达谱鉴定高密度脂蛋白缺乏小鼠中表达改变的基因。基因组研究。2000年;10:2022–2029. [PMC免费文章][公共医学][谷歌学者]
  • Benjamini Y,Hochberg Y。控制错误发现率:一种实用而强大的多重测试方法。J R Stat Soc Ser B公司。1995;57:289–300. [谷歌学者]
  • Storey JD公司。错误发现率的直接方法。J R Stat Soc Ser B公司。2002年;64:479–498. [谷歌学者]
  • Jörstad TS,Langaas M,Bones AM。了解样本大小:什么决定了实验所需的微阵列数量?植物科学趋势。2007;12:46–50.[公共医学][谷歌学者]
  • Gadbury GL、Page GP、Edwards J、Kayo T、Prolla TA、Weindruch R、Permana PA、Mountz JD、Allison DB。高维生物学中的功率和样本量估计。统计方法医学研究。2004;13:325–338. [谷歌学者]
  • Müller P,Parmigiani G,Robert C,Rousseau J.多重测试的最佳样本大小:基因表达微阵列案例。美国统计协会。2004;99:990–1001. [谷歌学者]
  • Jung SH.微阵列数据分析中FDR控制的样本量。生物信息学。2005;21:3097–3104.[公共医学][谷歌学者]
  • Li SS,Bigler J,Lampe JW,Potter JD,Feng Z.微阵列的FDR控制测试程序和样本量测定。统计医学。2005;24:2267–2280.[公共医学][谷歌学者]
  • Pawitan Y、Michiels S、Koscielny S、Gusnanto A、Ploner A。微阵列研究的错误发现率、灵敏度和样本量。生物信息学。2005;21:3017–3024.[公共医学][谷歌学者]
  • Tibshirani R.评估微阵列实验中样本大小的简单方法。BMC生物信息学。2006:7. [PMC免费文章][公共医学][谷歌学者]
  • Liu P,Hwang JTG。快速计算样本量,同时控制错误发现率,并应用于微阵列分析。生物信息学。2007;23:739–746.[公共医学][谷歌学者]
  • Ferreira JA,Zwinderman AH。使用Benjamini-Hochberg方法进行近似功率和样本量计算。国际生物统计杂志。2007;2:第8条。 [谷歌学者]
  • 胡杰,邹峰,赖特。微阵列实验中基于FDR的实用样本量计算。生物信息学。2005;21:3264–3272.[公共医学][谷歌学者]
  • Pounds S,Cheng C.错误发现率的样本量测定。生物信息学。2005;21:4263–4271.[公共医学][谷歌学者]
  • Storey JD公司。正误发现率:贝叶斯解释和q值。Ann统计。2003;31:2013–2035. [谷歌学者]
  • Pawitan Y、Murthy KRK、Michiels S、Ploner A.基因芯片研究中错误发现率估计的偏差。生物信息学。2005;21:3865–3872.[公共医学][谷歌学者]
  • 林赛·BG。混合可能性的几何:一般理论。Ann统计。1983;11:86–94. [谷歌学者]
  • Dempster AP、Laird NM、Rubin DB。不完整数据通过相对长度单位算法。J R Stat Soc Ser B公司。1977;39:1–38. [谷歌学者]
  • Johnson NL、Kotz S、Balakrishnan N。连续单变量分布。第二。第2卷。约翰·威利父子公司;1995[谷歌学者]
  • McLachlan G,Peel D。有限混合模型。约翰·威利父子公司;2000[谷歌学者]
  • 杰弗里斯·H,杰弗里斯理学学士。数学物理方法。第三。剑桥大学出版社;1972[谷歌学者]
  • Schweder T,Spjötvoll E.绘制p值以同时评估许多测试。生物特征。1982;69:493–502. [谷歌学者]
  • Allison DB、Gadbury GL、Heo M、Fernandez JR、Lee CK、Prolla TA、Weindruch R.微阵列基因表达数据分析的混合模型方法。计算统计数据。2002年;39:1–20. [谷歌学者]
  • Langaas M,Lindqvist BH,Ferkingstad E.估计真零假设的比例,并应用于DNA微阵列数据。J R Stat Soc Ser B公司。2005;67:555–572. [谷歌学者]
  • Akaike H.信息理论和最大似然原理的扩展。收件人:Petrov BN,Csaki F,编辑。第二届信息理论国际研讨会。1973年,第267-281页。[谷歌学者]
  • Schwarz G.估算模型的维数。Ann统计。1978年;6:461–464. [谷歌学者]
  • Storey JD公司。Ge、Dudoit和Speed受邀就“DNA微阵列数据分析的基于重采样的多重测试”发表评论。测试。2003;12:1–77. [谷歌学者]
  • 新泽西州海姆。计算最近的相关矩阵——这是金融学的一个问题。IMA J数字分析。2002年;22:329–343. [谷歌学者]
  • Smyth GK公司。Limma:微阵列数据的线性模型。作者:绅士R、凯里V、杜多伊特S、伊里扎里R、胡贝尔W,编辑。使用R和生物导体的生物信息学和计算生物学解决方案。纽约州施普林格;2005年,第397–420页。[谷歌学者]
  • Rocke DM,Durbin B.基因表达阵列测量误差模型。计算机生物学杂志。2001;8:557–569.[公共医学][谷歌学者]
  • Smyth GK公司。微阵列实验中评估差异表达的线性模型和经验贝叶斯方法。统计应用基因分子生物学。2004;:第3条。[公共医学][谷歌学者]

文章来自BMC生物信息学由以下人员提供BMC公司