引言

今天给大家带来一个关于“撒币“的问题，说起”撒币“，最经典就是二项分布，也就是你拿起一枚硬币抛向空中，当它转体n个360°，最后华丽的落在地上，出现正面或者反面的概率。具体的一些计算公式各位可自行找度娘索要。

在这里插入图片描述

除“撒币“问题，在数学中也有一个很好玩的实验，叫“取球实验”，也是就是在盒子里有黑球黑球黑球黑球白球白球白球…，当我们取n个球时，黑球数目为m时的概率，这个实验换个高级点的词汇，就是超几何分布。超几何分布和二项分布同为离散概率分布，超几何分布描述的是不放回的取样，而二项分布的话是放回取样。

超几何分布

在这里，我们为什么要提这个超几何分布呢，这就要从我们的转录组差异基因的下游分析说起。在差异基因下游分析中，GO富集分析和KEGG分析是最常见的，其Pvalue计算都是基于超几何分布。

首先，我们来看超几何分布的计算公式：

超几何分布计算公式

在这里呢，我们可以将N，M，n和k转化为的差异基因的参数，这也便于我们理解怎么计算某个GO term的Pvalue。

• N：基因组中与该GO term同属于同一层面（BP、CC或者MF）的基因数目
• M：基因组中含有该GO term的基因数目
• n：差异基因中与该GO term同属于同一层面（BP、CC或者MF）的基因数目
• k：n中含有该GO term的基因数目

Python计算超几何分布

目前python中Numpy的random包中提供了产生超几何分布结果的函数：
numpy.random.hyermetric(ngood, nbad, nsample, size=None)
具体可见这里.

• ngood：做出好的选择的数量，也就是上面的含有该GO term的基因数目，即M
• nbad：做出坏的选择的数目，也就是上面的不含有该GO term的基因数目，即N-M
• nsample：取样的数目，也就是差异基因的数目，即n
• size：采样的组数，即试验组数组数

上面的hypermetric()函数返回一组size大小的数组,数组中的每个数是在一组采样中合格产品的数量,即上面的k。
如：

s=np.random.hypergeometric(10,20,5,size=1000)
print(s)
p1=sum(s>=4)/1000
print(p1)

返回结果如下：

代码返回结果

我们可以看到s返回的一个列表，包含了1000组试验中抽到好样品的数目。其次，当我们选取好样品数>=4时，概率为0.036。

根据这个原理的话，我们也可以自定义一个函数，用于计算超几何分布的概率。

def p_cal(N,M,n,k,x):
    goodList=[]
    knum=0
    for i in range(0,x):            #循环x次，与hypergeometric的size参数相同
        good = 0
        for j in range(0,n):
            random_num=random.randint(1,N)  #随机生成1-N的随机数，当它<=M时，认为是好的样品
            if random_num<=M:
                good+=1
        goodList.append(good)
        if good>=k:
            knum+=1
    #return(goodList)  ##也可以返回类似上面s的列表
    return(knum/x)

这里的话我们就是采用了与超几何分布的相同的参数，我们利用相同的数字来测试：

s=np.random.hypergeometric(10,20,5,size=1000)
p1=sum(s>=4)/1000
print(p1)
p2=p_cal(30,10,5,4,1000)
print(p2)
p3=p_cal(30,10,5,5,1000)
print(p3)

输出结果：

Pvalue

我们可以看到概率相差不大，这表明可以达到超几何分布概率的计算。当知道怎么算pvalue，我们也离富集分析更进一步啦！

---

[图片上传失败...(image-9ad0db-1588929912338)]