【生物信息学】基因富集分析enrichment

conda create -n DL python=3.7 

conda activate DL

pip install torch==1.8.1+cu102 torchvision==0.9.1+cu102 torchaudio==0.8.1 -f https://download.pytorch.org/whl/torch_stable.html

conda install matplotlib

conda install scikit-learn

conda install pandas

conda install seaborn

conda install networkx

conda install statsmodels

pip install pyHSICLasso

from scipy.stats import hypergeom
import numpy as np
import pandas as pd
from statsmodels.stats.multitest import fdrcorrection
import matplotlib.pyplot as plt

[M, n, N] = [8000, 100, 300]

x = np.arange(0, n + 1)
pmf_deg = []
for i in x:
    pmf_deg.append(hypergeom.pmf(i, 8000, 100, 300))
p_A = np.sum(pmf_deg[30:300])

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, pmf_deg, 'bo')
ax.vlines(x, 0, pmf_deg, lw=1)
ax.set_xlabel('# of genes in pathway A of chosen deg')
ax.set_ylabel('hypergeom PMF')
plt.show()

all_pathway_canonica_GO_bp_unify = pd.read_csv("all_pathway_canonica_GO_bp_unify.csv", index_col=0, header=0)
deg_gene = pd.read_csv('deg_gene.csv', index_col=0, header=0)
gene_symbol = pd.read_table("gene_symbol_list.txt", header=None)all_pathway_canonica_GO_bp_unify = pd.read_csv("all_pathway_canonica_GO_bp_unify.csv", index_col=0, header=0)
deg_gene = pd.read_csv('deg_gene.csv', index_col=0, header=0)
gene_symbol = pd.read_table("gene_symbol_list.txt", header=None)
gene_symbol.columns = ['gene']
gene_symbol['label'] = 0
gene_symbol.index = gene_symbol['gene']
gene_symbol.loc[deg_gene['0'], 'label'] = 1

def Enrichment(gene_symbol, all_pathway_canonica_GO_bp_unify):
    data = gene_symbol
    up = data[data['label'] == 1].index.tolist()
    all_pathway = all_pathway_canonica_GO_bp_unify
    all_symbol = data.index
    all_list = np.unique(all_pathway['V1'])
    pathway_enrich_up = pd.DataFrame()
    pvalue = []
    for i in range(len(all_list)):
        # 获取每一个通路的基因集合
        gene_pathway = all_pathway[all_pathway['V1'] == all_list[i]]['V2']
        # 通路基因与差异表达基因交集
        up_no = len(set(gene_pathway).intersection(set(up)))
        # 背景基因与通路基因的交集
        hit_no = len(set(gene_pathway).intersection(set(all_symbol)))

        # 做统计模型时，统计不是分布概率，而是累积概率，所以计算的p_value 值计算公式为P(X＞x-1)或者 1- P（X≤x）
        # sf(k, M, n, N) = 1-cdf
        # cdf ，hypergeom.cdf表示：总共有M件产品，n件次品，从M件中随机挑出N件，这N件中最多包含n件中的k件的概率
        # M  is the total number of objects
        # n is total number of Type I objects
        # Scipy中的stats.hypergeom.sf(x, m+n, m, k)相当于R中的phyper(x, m, n, k, lower.tail=FALSE)(P X>x)
        pval_up = hypergeom.sf(up_no - 1, len(all_symbol), hit_no, len(up))
        pvalue.append(pval_up)
        res = pd.DataFrame((all_list[i], pval_up, up_no, hit_no, set(gene_pathway).intersection(set(up)))).transpose()
        pathway_enrich_up = pd.concat([pathway_enrich_up, res], axis=0)
    pathway_enrich_up.columns = ['pathway', 'pval_up', 'up_no', 'hit_no', 'gene_list']


    # 在富集分析中，每一个结果都进行了很多次的差异比较（一个通路一次），这种多重比较下的假阳性会急剧升高（这个假阳性的比例：FDR，false discovery rate，
    # 其含义是拒绝零假设的事件中错误拒绝事件的所占比例），因此要校正p值。
    # 举个栗子，如果检验1000次，设定的阈值为0.05（5%），那么无论得到多少个显著富集通路，这些显著富集通路出现假阳性的概率保持在2%之内，这就叫FDR＜0.02
    # BH: Benjamini & Hochberg法,简称BH法。将n个p值按照从小到大排序，k为p值的顺序值，找到一个最大的k使得p*n/k ＜α ，认为1,2，...k个通路是显著富集的
    # 最大的p值保持不变，对于排名n-1的，p_adj（n-1） = min( p * n / k, p_adj（n）)。其影响是除了最大的p值保持不变外，所有的其他p值都会被增加，变得没那么显著
    fdr = fdrcorrection(pvalue, alpha=0.05, is_sorted=False)
    p_adjust = fdr[1]
    pathway_enrich_up['fdr_p_adjust'] = list(p_adjust)
    final_pathway_enrich_up = pathway_enrich_up.sort_values(by='pval_up', ascending=True)
    return (final_pathway_enrich_up)

res = Enrichment(gene_symbol, all_pathway_canonica_GO_bp_unify)

res.to_csv('enrichment.csv')
print(res)