一种准确估计差异表达基因识别中错误发现率的通用方法

摘要

总结：基因组数据、转录组数据、蛋白质组数据和单核苷酸多态性数据等“组学”数据一直在快速增长。经济数据是大规模、高通量的数据。这些数据挑战了传统的统计方法，需要进行多次测试。已开发出Bonferroni程序、Benjamini–Hochberg（BH）程序和Westfall–Young程序等多项测试程序，其中一些控制家族错误率，另一些控制错误发现率（FDR）。这些程序在某些情况下有效，不能应用于所有类型的大规模数据。为了解决组学数据分析中这一具有统计挑战性的问题，我们提出了一种生成一组多重测试程序的通用方法。该方法基于BH定理。通过选择C值，可以实现特定的多重测试程序。例如，通过设置C=1.22，我们的方法生成BH过程。当C<1.22时，我们的方法生成弱控制FDR的过程，当C>1.22时程序强控制FDR。C=G（基因或测试数量）和C=0的分别是Bonferroni程序和单一测试程序。这是这个家庭中的两个极端程序。为了让人们在实践中选择合适的多次测试程序，我们开发了一种算法，通过该算法可以正确可靠地估计FDR。仿真结果表明，我们的方法在各种场景下都能很好地精确估计FDR，并且我们用三个实际数据集说明了我们的方法的应用。

可用性和实施：我们的程序在Matlab中实现，可根据要求提供。

联系人:hxu@gru.edu

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

1简介

“组学”技术的进步导致大规模数据的巨大发展，如微阵列数据、转录组数据、代谢组数据和蛋白质组数据。这些“组学”数据使我们能够全面了解复杂的生物过程、遗传和环境因素之间的相互作用以及复杂疾病的病理机制，如糖尿病、心脏病、高血压和各种癌症，通过识别差异表达的基因，对功能基因、通路或网络进行聚类或分类。这种大规模数据或高通量数据挑战了传统的统计检验(埃夫隆等。, 2001;牛顿等。, 2001;平移等。, 2003;棕褐色等。, 2006;塔瑟等。, 2001)因为传统的统计测试都是单项测试。在单假设检验中，阈值α对于确定检验是否显著有意义，因为α是人群中事件发生的概率阈值。然而，阈值α=0.05对于组学数据的多重统计分析仍然无效，因为例如，在鉴定10000个基因时，在显著性水平为0.05时，至少有500个偶然发现(尼科尔斯和哈亚斯卡，2003年;塔瑟等。, 2001). 因此，需要针对多次测试调整阈值α。目前，已经开发了几种用于调整阈值α或调整的统计程序P（P）-多个测试中的值。其中，Bonferroni程序，Holm程序(霍尔姆，1979年)Hochberg程序(霍克伯格，1988年)，Westfall和Young程序(Westfal and Young，1993年)Sidak程序(Nichols和Hayasaka，2003年)Benjamini–Hochberg（BH）程序(本杰米尼和霍奇伯格，1995年)和Benjamini–Liu程序(Benjamini和Liu，1999年)是典型的多重测试方法。在这些程序中，BH程序在实践中应用最广泛。但我们发现，真假发现率（FDR）并不总是与BH程序的估计FDR保持一致，所有现有的多重测试程序实际上都是多重测试程序家族的特例，范围从单一测试程序到Bonferroni程序；换句话说，这些程序给出了不同的拒绝空间(图1). 在实践中，我们的模拟还表明，BH程序通常会大大高估或低估真实的FDR。为了在大规模统计分析中准确估计FDR，我们提出了一种基于BH定理的方法，在Bonferroni程序的一系列单一测试程序中生成一组多重测试程序，并提出了一个模拟算法来选择所需的程序。

2方法

2.1假设检验的拒绝区

在这里，我们重点讨论了单一测试程序[A（s）]、Bonferroni[A（B）]和BH多重测试程序[B（BH）]的拒绝范围。如中所示图1，用于G公司待测试假设中，单一测试程序的拒绝范围最大：哪里α是一个概率阈值，而Bonferroni多次试验程序的拒绝面积最小：这是两个极端的拒绝区域。在单一测试程序中，假设会被偶然拒绝，而Bonferroni多重测试程序有机会拒绝G公司假设，这被称为家庭错误率。BH多次试验程序的拒收面积是对角线下的面积，是单次试验程序的一半，哪里我是我-假设。但不同于单一测试程序为所有假设和Bonferroni多重测试程序设置，其中/G公司对于所有假设，BH多重测试程序是一个程序，从G公司到1，调整后的值位于对角线上，其拒绝区域为/2.因此，伯克希尔哈撒韦多重测试程序介于单一测试程序和邦费罗尼多重测试程序之间。因此，从理论上讲，可以给出具有拒绝区域的任何多次试验程序到.

2.2 Benjamini-Hochberg定理

对于米无效假设，测试它们有相应的P（P）-值.让服从命令P（P）-值及其相应的有序空假设然后，我们有一个BH程序：

让k个成为最大的我对于其中，然后全部将被拒绝。 

 

2.3功能

方程式（9）给出了从单一测试程序到Bonferroni程序的多个程序系列。C值是所谓的控制值，因为它对; 换句话说，不同的C值对假阳性结果的控制程度不同。例如，当C=0时，则和为所有人P（P）-值，这是单一测试程序；什么时候，那么和或C=G>1000，和为所有人P（P）-值，这只是Bonferroni程序；如果C=1.22，我们有，这只是BH程序(图2). 因此，在C类> 1.22,方程式（9）生成一组比BH程序具有更强的FDR控制的多重测试程序，而C类<1.22，它生成了另一组多重测试程序，与BH程序相比，FDR控制较弱(图2). 因此，该方法不仅可以实现现有的多次试验程序，而且可以通过选择适当的C值来产生新的程序(图2).

2.4 FDR控制的C值选择

为了控制FDR，我们需要选择一个C值来生成多重测试过程。理论上，可以使用置换方法(雷纳等。, 2003)确定多次试验程序。但正如我们将在第3节中看到的，排列并不是一个好方法，因为它往往会大大低估罗斯福。因此，我们在此提出了一种替代方法，用于为实际数据集选择理想的C值。为了便于讨论，我们将算法集中在传统的两个条件（或两类）上t吨-测试方法。

我们的算法包括以下九个步骤，其中第2步是计算每个基因的均值、方差和距离，第3步和第4步是使用参数引导程序模拟数据，第5步和第6步是使用一组多个程序估计FDR： 

        

绘制由M多项测试程序确定的实际和估计比率集，并拟合具有交集的线性回归模型和回归系数.如果0.1和，则估计的比率集的结果几乎与M多次测试过程中的实际比率集相同，并且真实的FDR集也将与估计的FDR集合一起给出；C值将通过选择校准的FDR来选择。否则，通过重置参数重复步骤4–7一、b和直到0.1和.

2.5统计方法的性能

使用baySeq包进行经验贝叶斯分析；在edgeR软件包中实现了Fisher精确检验和广义线性模型（GLM）。程序包baySeq和edgeR从Bioconductor下载。这两个组件使用BH程序进行调整P（P）-值（以估计FDR）。贝塔t吨-在Matlab中进行了测试。我们还使用BH程序来调整P（P）-值。使用从Bioconductor下载的qvalue包执行q值方法。此软件包具有GUI界面。我们加载了P（P）-两个样本的值t吨-测试qvalue.gui接口并执行q值计算，获得q值。我们选择q值≤0.05的DE基因。

我们的多重过程方法和C值选择算法在Matlab中实现。

3结果

3.1仿真评估

我们使用塔瑟等。的（2001年）微阵列数据来说明该算法。我们从网站上下载了数据http://www-stat.stanford.edu/～tibs/SAM/数据包括7129个基因和8个人类淋巴母细胞系样本，其中4个样本用电离辐射处理，其余样本用非电离辐射处理作为对照。从实际的微阵列数据中，获得了两种条件下每个基因表达值的平均值和方差，并使用一个正常的伪随机生成器、两个平均值中的任何一个、每个基因的两个标准差中的任意一个和四个重复来创建两个具有相等方差的7129个基因的噪声表达数据集。然后，通过线性和随机分配条件影响生成模拟数据关于两种条件对两个噪声数据集中30%基因的差异表达U型是（0,1]中的一个统一变量。请注意，我们的模拟程序在分配其表达值时在这30%的基因之间产生了强烈的依赖性。我们应用了一个经典的两个条件t吨-对该模拟数据集进行测试，并将M=140的C值从0.01设置为1.4，以创建140个用于多次测试的程序，并获得140个结果比率。通过应用上述算法并调整a、b和值，我们发现当一= 0.001,b条=5和=0.054，140个调查结果的估计比率和实际比率非常接近。因此，结果在15次重复模拟中平均，并显示在图3.图3A表明，估计值与实际发现比率的绘图线几乎是一条带交点的对角线=−0.009≈0和回归系数=1.0576≈1，因此估计的比率结果集与实际结果几乎相同。发件人图3B、 还可以看到，估计的FDR与真实的FDR的图点几乎在对角线上，这表明估计的FDR在每个FDR点处与其真实值大致相同。

接下来，我们用负二项分布模拟转录组数据。我们通过设置13000个mRNA亚型、两种条件（控制和疾病状态）、每种条件下的六个复制文库、概率来进行模拟=U型哪里U型是随机均匀变量，尺寸=1000U型对于每个亚型。负二项伪随机生成器和这些设置参数用于创建两个方差相等的噪声转录组数据集，然后对于条件（或治疗），mRNA亚型差异转录的影响值随机线性分配给两个噪声数据集中10%的mRNA亚型别，其中u个=（0,1]是统一变量。因此，获得了10%差异转录亚型、六次重复和两种条件的模拟数据。有一批现有方法可用于寻找两种条件下差异转录的mRNA亚型：经验贝叶斯方法(Hardcastle和Kelly，2010年)，Fisher精确测试(Robinson和Smyth，2008年)、GLM(麦卡锡等。, 2012)和Betat吨-测试(巴格利等。, 2003). 这些方法都使用BH程序估计FDR。类似于图3B、 我们在这里使用模拟转录组数据集来评估基于四种统计方法的BH程序。我们在这些方法的结果中绘制了真实的FDR，并使用BH程序从cuttop=～0到～0.21绘制了估计的FDR。图4结果表明，在经验Bayesian、GLM和Fisher精确检验中，尤其是在Beta中，BH估计的FDR曲线t吨-测试结果远低于其理论线（真实FDR与真实FDR），表明在这些方法中，BH程序严重低估了FDR。

为了在各种微阵列实验中充分显示我们的多重测试方法的统计特性，我们还生成了27个场景中的模拟数据集。我们的模拟基于塔瑟等。的（2001）微阵列数据。我们设置了三个条件（或治疗）效应水平：A=100、200和300；差异表达基因的比例分为三个等级：P=10、20和30%，样本量分为三级：n=4、6和18。条件效应随机分配到基因的P。这种分配在基因之间产生了强相关性表达。我们的多步骤方法和FDR估计算法应用于这些模拟数据集。由于有大量模拟数据集，我们无法显示调查结果的估计和实际比率以及估计和真实FDR的概况，如所示图3A和B，针对每个场景；相反，我们列出了通过传统t检验方法获得的结果数量，以及通过选择一个程序获得的真实和估计的FDR数量，该程序的估计FDR小于但最接近于0.05的截止值。为了演示我们的算法，我们还应用了置换(雷纳等。, 2003)和q值(Storey和Tibshirani，2003年)这27个模拟数据集的处理方法。结果总结如下补充表S1.

可以在中看到补充表S1所选C值在0.9–1.32范围内变化，包括C=1.22。如中所示图1和2，C=1.22生成BH程序。因此，BH程序只是我们的多个程序之一。在每个给定场景中，我们的方法在每个C值处估计的FDR接近其真实值，但在场景3和场景17中，真实FDR值略大于其估计值；所有其他FDR都被稍微高估了，这意味着我们对FDR的估计是保守的。补充表S1还给出了通过执行q值方法获得的结果(Storey和Tibshirani，2003年). 我们的方法和q值方法在大多数模拟数据集中都有相似数量的发现。然而，q值方法在九种情况下低估了FDR（33.3%）。因此，由于FDR估计的保守性，我们的方法明显优于q值方法。通过使用置换算法，FDR在所有27种情况下都被严重低估(补充表S1)，表明置换不是估计FDR的理想方法。

3.2应用

应用多重过程方法的第一个示例是Tusher公司等。的（2001）微阵列数据。如上所示，数据由7129个基因组成。塔瑟等。(2001)试图在四个经电离辐射处理的人类淋巴母细胞系样品和四个未经电离辐射的样品之间找到差异表达的基因。我们将140个C值从0.1设置为1.39，以创建140个多次测试程序。由于C值>0.92没有产生有意义的结果（没有发现），我们使用经典的两个条件从实际数据中获得了83个发现的实际比率t吨-试验方法。调整后，我们发现一=0.02，b=2.9和=0.3，并重复我们的算法过程15次，生成一组83个估计的发现比率和一组83个估计的FDR。我们在图5回归线几乎是一条有交点的对角线= 0.00080和回归系数= 0.9451≈1. 因此，结果的估计比率与每个C值的实际比率几乎相同。因此，中给出的模拟FDR图5A2可用于估计83个C值的真实FDR。图5A2表明，当C值=0.92时，最小FDR为0.36，远大于截止值0.05。因此，在FDR截止值为0.05时，在经过电离辐射处理的人类淋巴母细胞系和对照细胞系之间，未发现该数据集中的基因存在差异表达。该结果与塔瑟等。(2001)其中，SAM中估计的最小FDR为0.12。

第二个例子是白血病微阵列数据下载自http://www.broadinstitute.org/cgi-bin/cancer/datasets.cgi该数据集由12582个基因组成，通过寡核苷酸阵列检测了它们在24个急性淋巴细胞白血病（ALL）细胞系、20个髓系淋巴细胞白血病细胞系和28个急性髓系白血病（AML）细胞株中的表达值。在这里，我们选择了两种白血病细胞类型ALL和AML进行差异表达测试。这是两个样本大小较大且不相等的条件的示例。对于这个数据集，我们输入了140个从0.1到1.39的C值，并从经典的不成对双条件中生成了140个实际比率t吨-测试。调整后，我们获得一= 0.5,b条=1.5和=0.05，并应用我们的算法通过140个程序估计一组FDR，重复15次。结果如所示图5B.如中所示图5B1，调查结果的估计与实际比率的绘图点几乎落在具有交叉点的对角线上和回归系数这表明，在每个C级水平上，估计的发现比率几乎与实际比率相同。因此，真实的FDR集合可以由估计的FDR集合给出。图5B2显示，当C=1.31时，FDR=0.049723接近0.05，在ALL和AML细胞系之间鉴定出582个基因差异表达。

第三个例子是拟南芥微阵列数据。探索是否拟南芥基因对由转移DNA（T-DNA）编码的癌基因或由农杆菌属进入植物细胞，李等。(2009)接种受伤幼体后3h和6d进行微阵列实验拟南芥含有C58株和GV3101株的植物，GV310株是C58株的同源株，仅缺乏转铁蛋白-DNA，但具有病毒D2、病毒E2、病毒E3和病毒F等蛋白毒力因子(Vergunst公司等。, 2000,2003). 受伤但未感染的秸秆被用作对照。例如，我们刚刚从GEO（GEO登录：GSE14106）网站下载了接种后6天的数据(网址：http://www.ncbi.nlm.nih.gov/geo/). 该数据由22810个基因组成，分别来自受感染的C58和对照秸秆，每个秸秆有三个重复。我们执行了经典的两个条件t吨-在此数据集上测试方法，并应用我们的多过程方法进行多次测试。我们将140个C值从0.01设置为1.4，以创建140个程序，并从两种条件中获得140个实际发现比率t吨-测试。调整后，我们发现一= 0.08,b条=15和= 0.035. 使用这些参数，我们通过140个程序对15个重复获得了140个估计的发现和FDR比率。图5C1显示，观察到的估计与实际发现比率的绘图线与预测线完全重叠。因为十字路口和回归系数，估计结果与实际结果比率的观察线是一条对角线，因此，在每个C值处，估计结果比率与实际比率是一一相同的。因此，真正的FDR集可以由估计的FDR集合给出。图5C2显示，在C=1.2，FDR=0.0497的条件下，在年轻肿瘤和对照茎之间鉴定出1427个基因的差异表达(表1). 然而，使用折叠变换方法，李等。(2009)刚刚鉴定出196个差异表达基因P（P）< 0.01. 两者都有很大的不同。为了验证我们的发现，我们显示了1427个DE基因的原始数据的热图图6.图6结果表明，在下调过程中，所有498个基因在年轻冠瘿瘤中的表达值明显低于对照（参考）茎秆（分别以深红色和黄色显示），而在上调过程中，全部929个基因在冠瘿肿瘤中的表达也明显高于对照茎秆。也就是说，在1427个基因中，肿瘤和对照明显由基因表达分开。在这三个重复中，1427个基因有差异表达。

4讨论

如中所示表1，我们的方法有点类似于SAM(塔瑟等。, 2001)和RAM(棕褐色等。, 2006)因为一组C值（过程）等同于一组阈值()用于声明SAM和RAM中观察到的排名统计数据和估计的排名空统计数据之间的重要测试。但是，如中所示图3和中补充表S1我们的方法通过约束发现的估计比率与每个给定C值的实际比率相同或近似，从而确保估计的FDR接近其真实值，而SAM和RAM没有任何约束来确保FDR的估计是可靠的。27个模拟数据集表明，我们的方法对于准确估计不同样本大小、不同条件效应和不同比例基因在不同研究生物系统中差异表达的FDR非常有效。三个实际数据集也表明，即使在样本大小较大且不相等或样本大小很小但阵列上检测到大量基因的情况下，我们的方法也表现良好。

在微阵列和转录组学实验中，基因的数量极为庞大（特征空间的维数达到数十甚至数十万），而生物条件的数量极低，导致基因在表达上高度相关。这意味着这种全基因组数据是高度依赖的数据。我们的模拟和实际结果表明，基因表达之间的依赖性并不影响FDR的估计，即使多重程序基于假设的独立性（BH定理）。这是因为对发现比率的估计也是基于高度相关的基因表达。

理论上，我们的递增或递减多重过程不仅适用于任何统计方法，而且我们的C值选择算法（过程）也可以推广到任何其他统计测试，或者通过改变分布适用于任何数据。例如，如果实际数据集是计数数据，例如转录组数据或基因表达序列分析（SAGE）数据，则可以通过在泊松、伽马、二项式或负二项式分布上选择适当的参数来进行我们算法中的模拟。但在实践中，我们的方法不适用于统计分布不明确的情况。

此外，与现有的大规模统计方法如SAM相比(塔瑟等。, 2001)、RAM(棕褐色等。, 2006)，ODP公司(故事等。, 2007)、Cyber T(Baldi and Long，2001年)，利马(Smyth，2004年;斯迈思等。, 2005)和q值(Storey和Tibshirani，2003年)，我们的方法需要更多的时间来确定参数值，以便在选定的一组C值中，发现的估计比率与实际比率相同或非常接近，或者说，接近0并且至1。实际上，,一和b条可以通过注意以下内容轻松确定值，如中所示方程式10，较大的-值将允许调整更大的d值和更小的一个-值会使这些d值变小，然后变大-M端附近的值会变小；根据公式12，b值越大，两个条件之间所有微分表达式的差异越大，结果会越多-靠近1端的值将被放大。我们的经验是，首先我们看一下-值集和-使用设置的值=一=b条=1并确定b值，以便在1端-值接近-值。然后我们设置了一个小一-值和调整因此-值和-M端的值小于公差值，例如0.01。此外，可以先修复较小的-值，然后调整一-值。然而，约束条件是接近0并且至1。

确认

作者感谢两位匿名评论员提供了许多建设性建议。

基金：H.X.得到了佐治亚摄政大学的校内资助。

利益冲突：未声明。

参考文献

作者注释