如何使用 R 对基因列表重叠应用多重测试校正

机器算法验证 r 多重比较 遗传学 邦费罗尼 超几何分布
2022-03-23 03:04:36

我有 2 项研究着眼于患者对相同药物的反应。研究 1 发现 10,000 个基因在背景之上表达,其中 500 个基因存在差异表达,称为药物反应特征。研究 2 发现了 1,000 个代表药物反应特征的基因。两个签名之间的重叠是 100 个基因。

我想计算签名之间重叠的统计显着性。如果我理解正确,一种方法(基于此处的帖子:计算 RNA seq 和 ChLP 芯片数据集之间的基因列表重叠的概率和此处:使用 R 的 phyper 获取列表重叠的概率)是通过phyper()

> overlap  <- 100
> list1    <- 500
> totalPop <- 10000
> list2    <- 1000
> 
> 1-phyper(overlap-1, list1, totalPop-list1, list2)
[1] 4.103051e-12
  1. 这听起来合理吗?

  2. 如果我想应用 Bonferroni 校正,我需要将此 p 值乘以比较次数。在这种情况下,比较次数对应于多少?清单 2?或者,有什么方法可以快速进行不那么保守的校正(例如,Benjamini-Hochberg)?

2个回答

我对基因表达研究一无所知,但我确实对多重推理感兴趣,所以无论如何我都会冒险回答这部分问题。

就个人而言,我不会以这种方式解决问题。我会调整原始研究中的错误级别,计算新的重叠并将测试留在最后。如果差异表达基因的数量(以及您正在使用的任何其他结果)已经基于调整后的测试,我认为您不需要做任何事情。

如果您无法返回原始数据并且确实想要调整p值,您确实可以将其乘以测试次数,但我不明白为什么它应该与 list2 的大小有关。调整两项研究中进行的测试总数(即人口的两倍)会更有意义。不过,这将是残酷的。

要调整R 中的p值,您可以使用p.adjust(p),其中pp值的向量。

p.adjust(p, method="bonferroni") # Bonferroni method, simple multiplication
p.adjust(p, method="holm") # Holm-Bonferroni method, more powerful than Bonferroni
p.adjust(p, method="BH") # Benjamini-Hochberg

如帮助文件中所述,没有理由不使用 Holm-Bonferroni 而不是 Bonferroni,因为它在任何情况下也提供了对家庭错误率的强大控制,但更强大。Benjamini-Hochberg 控制错误发现率,这是一个不太严格的标准。


在下面的评论后编辑:

我对这个问题思考得越多,我就越认为在这种情况下对多重比较进行校正是不必要和不合适的。这就是假设“家庭”的概念出现的地方。您的最后一个测试与所有早期测试都不太可比,没有“利用机会”或挑选重要结果的风险,只有一个感兴趣的测试,并且对此使用普通错误级别是合法的。

即使您对之前进行的许多测试进行了积极的纠正,您仍然不会直接解决主要问题,即两个列表中的某些基因可能被错误地检测为差异表达。较早的测试结果仍然“站得住脚”,如果您想在控制家庭错误率的同时解释这些结果,您仍然需要纠正所有这些结果。

但是,如果所有基因的原假设都为真,那么任何显着的结果都将是假阳性,并且您不会期望在下一个样本中再次标记相同的基因。因此,两个列表之间的重叠只会偶然发生,而这正是基于超几何分布的测试正在测试的内容。因此,即使基因列表完全是垃圾,最后一次测试的结果也是安全的。直觉上,似乎介于两者之间的任何事情(真假假设的混合)也应该没问题。

也许在这个领域有更多经验的人可能会参与进来,但我认为只有在你想比较检测到的基因总数或找出哪些基因表达差异时才需要进行调整,即如果你想解释成千上万的个体在每项研究中进行的测试。

您无需为单次重叠检验校正 p 值。但是,假设您有兴趣确定药物是否会影响同一途径中的基因。您将如何确定哪条路径重叠最多?假设您有 500 个通路基因组。您运行超几何集合重叠测试 500 次并按 p 值对它们进行排序。由于您运行此测试 500 次(甚至更多取决于您拥有的数据量),因此您有可能只是偶然获得好分数(误报)。因此,您需要对此进行纠正并执行 pvalue 调整……Bonferroni(最保守的)或 Benjamini Hochberg。