如何在R中使用超几何分布计算p-value

问题,假设有两组ChIP-seq的数据A和B,我想知道它们之间是否相关。A组一共有50个数据,B组一共有60个数据,它们之间重叠的数据有30个,总计的binding位点如果是200个,那么如果它们相关的p-value如何计算呢?这个问题其实就是转变成抽到重叠数比30或者比30大的可能性是多大。所以我们需要计算抽到重叠数为29的可能性,然后用1减去它就可以了。
公式:p.value = phyper(p1.and.p2 -1, p2, totalTest-p2, p1, lower.tail = FALSE,log.p = FALSE)

> pval <- 1-phyper(30-1, 60, 200-60, 50)
> pval
[1] 2.555588e-07
> 1-phyper(30-1, 50, 200-50, 60)
[1] 2.555588e-07
> phyper(29, 60, 200-60, 50, lower.tail=FALSE)
[1] 2.555588e-07
> sum(sapply(30:50, dhyper, m=60, n=200-60, k=50))
[1] 2.555588e-07

解释:phyper函数返回一个概率,这个概率与lower.tail参数有关。当lower.tail的值为真时,phyper返回一个当抽取到的白球数小于等于q值的概率,当lower.tail值为假时,phyper返回一个当抽取到的白球数大于q值的概率。phyper的返回值是一个累积的概率值。

lower.tail logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].

函数dhyper返回当抽取到的白球数等q值的概率,所以计算从30至50的dhyper,然后加起来就一定等于phyper取lower.tail为假时的结果。

2 thoughts on “如何在R中使用超几何分布计算p-value

  1. 这里说的“A组一共有50个数据,B组一共有60个数据,它们之间重叠的数据有30个,总计的binding位点如果是200个”。对于这句描述我有些不解:

    1)50个数据,指的是什么?是指50个peak吗?

    2)200个binding site 指的是A和B的总和吗?如果A组和B组中的数据指的是peak,那不是只有50+60=110吗?

    3)最终的分析,是想说明两个factor的co-binding具有显著性吗?或者说两套同一factor的ChIP-seq 的 结果重复性很好吗?

    水平有限,请赐教。谢谢!

    1. 谢谢您的留言。依次回答您的问题如下:
      1)是指50个peaks
      2)200不是A和B的总和。是所有可能的binding位点。
      3)是指两组数据A和B的binding sites相关显著,它们可能是co-binding,也可能是重复实验。这里的结果只表明,A和Bbind到了同一(或者相邻)位置的可能性显著。

发表评论

电子邮件地址不会被公开。 必填项已用*标注