用python编写生词本 用python做生信

基于Python实现通路富集模型

• 前言
• 一：超几何分布介绍
• 二：富集原理
• 三：代码计算
• 四：总结

前言

本文章主要涉及基因功能富集分析的原理解释，统计检验以及最终基于Python代码的整体逻辑实现。富集分析应该算生信里是最常用的分析方法之一了，很多做生信的都是基于R或者Spss等软件，所以这次想用Python来回顾每一步处理任务。

一：超几何分布介绍

超几何分布时一种离散型概率分布，也许中学就学过最经典的例子：假设一个袋子有10个球，其中红色球6个，白色球4个，那么我一次从袋子中抓取4个，请问这4个球中至少3个是红球的概率是多少？

解：至少3个红球，那就是3个或者4个，设随机变量 X 为4个中含有红球的个数，即求解 $P_{res} = P(X=3)+P(X=4)$
通过超几何概率公式： $P(X=k)=\frac{C_M^{k}C_{N-M}^{n-k}}{C_N^n}$ ，这里 $N$ 等价于10个球， $n$ 等价于抓取的4个球， $M$ 等价于6个红色球，那么最后的 $k$

• 计算结果

我们先手动计算此结果，然后用python计算，由于此公式比较简单，直接调用相关函数即可。

手动结果 $P_{res} = P(X=3)+P(X=4)=\frac{80}{210}+\frac{15}{210}=0.45$

我们用代码函数计算

from scipy import stats

p1 = stats.hypergeom.pmf(k=3, M=10, n=6, N=4)
p2 = stats.hypergeom.pmf(k=4, M=10, n=6, N=4)
print(round(p1 + p2,2))###0.45

• 期望与方差

对超几何分布 $X\sim H(N,M,n)$ ，随机变量 $X$

$E(X)=nM/N$

方差由定义推导比较麻烦，我们直接给出

$D(X)=E(X^2)-E(X)^2=n\frac{M}{N}(1-\frac{M}{N})\frac{N-n}{N-1}$

二：富集原理

富集分析原理就是超几何分布的使用。假设人类背景基因3W个，已知共500个基因参与某一通路，由已知我们有一个自己的含有200个基因的基因群，且这200个基因都来自背景基因，恰巧这200个基因里面有6个参与此通路，那么这6个基因到底是不是偶然还是显著富集在此通路上？于是我们要做相应的显著性检验。

• p 值检验

原假设 $H_0$: 此一组基因并不显著富集在此通路上

$p$ 值计算公式： $p=1-\sum_{i=0}^{m-1}{\frac{C_{M}^{i}C_{N-M}^{n-i}}{C_{N}^{n}}}$ ，若 $p<0.5$，则为统计意义上的显

著，即拒绝原假设。

统计上，一般原假设的建立都是要被推翻的。

很显然上面的 $p$

• 温和的BH（Benjaminiand Hochberg）法

我们知道尽管计算出来的 $p$值往往都很小，一次“不好的出现”也许只有0.001，但如果重复1万次，那么再低的概率也许都能发生，按照平均期望最起码也能出现10次 。所以基于这种可能，在做检验时，我们也怕重复过多而出现不显著结果过多出现。

$BH$方法一定程度上能解决这种事件过多的出现。
$BH$ 法需要将总计 $M$ 次检验的结果按由小到大进行排序， $k$ 为其中一次检验结果的 $p$ 值所对应的排名。找到符合原始阈值 $\alpha$ 的最大的 $k$ 值，满足 $p(k)<\alpha*k/M$，认为排名从1到 $k$ 的所有检验存在显著差异，并计算对应的 $q$ 值公式为 $q=p*(M/k)$

举个例子就能很清楚地说明此修正过程：

图1

比如 $k=1$ ，则 $\alpha*k/M=0.05*1/5=0.025$ ， $p(1)=0.001<0.025$ 是显著的，同理 $k=2$ 时$\alpha*k/M=0.05*2/5=0.02>P_2=0.023$ 不再显著，则$G1$、$G4$、$G5$、$G2$都不是显著表达基因，那么我们就把 $k$ 找到了，再根据公式能算出对应的 $q$

三：代码计算

我们依旧使用上面图中的数据，看代码计算结果：

图2

从图2可以发现，计算过程完全正确。但是我们要注意一点，就是并不是任何 $p$ 值结果都能进行$BH$检验的，从不等式 $p(k)<\alpha*k/M$ 来看，此不等式并不是一直成立的（很有可能 $p$ 值过于大导致不等式不存在），那么我们只能用原始 $p$

def BH_test(alpha,pvals,gn):
    l = len(pvals)
    k = [i for i in range(1,l+1)]
    k = np.array(k)
    sort_id = np.argsort(np.array(pvals))
    pvals = sorted(np.array(pvals))##升序排序
    gn = np.array(gn)[sort_id]##基因名字按照p值从小到大排序
    P_k = alpha * k / len(pvals)##计算P(k)向量值
    q =  pvals * (l / k)##计算qvals
    comp_res = pvals < P_k
    idx = [i for i,x in enumerate(list(comp_res)) if x==1]
    df = pd.DataFrame([gn,pvals,P_k,comp_res,q],index=['order_name','pvals','P(k)','compres','qvals']).T
    try:
        loc = idx[-1] + 1###取最后的一个索引
        print(df)
        print('\033[1;33m对应位置坐标最大k值：\033[0m',loc)
        print('\033[1;34m最终筛选出基因跟疾病有关的疾病列表为(取一般置信度=%s水平下)：%s\033[0m'%(alpha,gn[:loc]))
        return df
    except Exception as ee:###异常处理
        print('\033[1;31m{0:*^80}\033[1;0m'.format('BH认证失败，可能是所有原始P_value过大，如下结果展示'))
        df = pd.DataFrame([pvals,P_k],index=['P_value','P(k)']).T
        print(df)
        return df

pls = [0.023,0.12,0.5,0.15,0.082]
namels = ['g1','g2','g3','g4','g5']
BH_test(0.05,pls,namels)

图3

四：总结

富集分析的原理其实很简单，在这里只是根据原理做相关总结实验，并仅仅提供些具体数值的计算代码。但实际运用中可能要考虑具体场景问题做相应的数据处理。