r - 使用R计算维恩图超几何p值

Question

嗨，我看到有人计算维恩图重叠 p 值，如下例所示。他们使用超几何分布和 R。当我在 R 中应用他们的函数时，我无法得到相同的结果。谁能帮我解决这个问题？

我在别人发表的文章中看到的样本:

从15220个基因中，A组为1850+195个基因，B组为195+596个基因，overlap为195个基因。他们的 p 值为 2e-26。

他们的方法是:给定总共N个基因，如果基因集A和B分别包含m和n个基因，其中k个是共同的，那么富集的p值计算如下:

p = Σ (m,i)(N-m,n-i)/(N,n)

对于 i 从 k 到 min(m,n)，其中“(m,i) "表示二项式。

我使用 R 的方式是:

sum(选择(596+195,195:(195+596))*选择(15220-596-195,(1850+195)-195:(195+596)))/选择(15220,1850 +195)。

我得到了 NaN。

或者使用:phyper(195,1850+195,15220-1850-195,596+195)，我得到了 1。

我也引用了链接 http://www.pangloss.com/wiki/VennSignificance但是当我计算

1 - phyper(448,1000,13800,2872) 在 R 中，我得到 0 而不是链接的 1.906314e-81。

我对 R 和统计完全陌生，很抱歉在这里发布了很多错误。

Answer 1

使用包 gmp，并将 choose 替换为 chooseZ，我们可以将您的 p 值实现为:

require(gmp)

enrich_pvalue <- function(N, A, B, k)
{
    m <- A + k
    n <- B + k
    i <- k:min(m,n)

    as.numeric( sum(chooseZ(m,i)*chooseZ(N-m,n-i))/chooseZ(N,n) )
}

结果:

> enrich_pvalue(15220, 1850, 596, 195)
[1] 1.91221e-18

使用您的 pangloss 链接中的示例(使用您的符号)，我们得到:

> enrich_pvalue(N=14800, A=1000-448, B=2872-448, k=448)
[1] 7.289388e-81