404_多组学通路分析¶

import pandas as pd
import numpy as np
from scipy import stats

# ========== 1. 准备多组学差异分析结果 ==========
# 模拟转录组差异基因
rna_results = pd.DataFrame({
    "gene": [f"GENE_{i}" for i in range(200)],
    "log2FC": np.random.normal(0, 1.5, 200),  # 对数倍数变化
    "padj": np.random.uniform(0.001, 0.5, 200)  # 校正后p值
})

# 模拟蛋白组差异蛋白
protein_results = pd.DataFrame({
    "protein": [f"GENE_{i}" for i in range(0, 150)],  # 前150个与基因重叠
    "log2FC": np.random.normal(0, 1.0, 150),
    "padj": np.random.uniform(0.001, 0.3, 150)
})

# 筛选显著差异特征
sig_rna = rna_results[
    (rna_results["padj"] < 0.05) & (rna_results["log2FC"].abs() > 1)
]["gene"].tolist()

sig_protein = protein_results[
    (protein_results["padj"] < 0.05) & (protein_results["log2FC"].abs() > 0.5)
]["protein"].tolist()

print(f"显著差异基因数: {len(sig_rna)}")
print(f"显著差异蛋白数: {len(sig_protein)}")

# ========== 2. 模拟KEGG通路基因集 ==========
# 实际应用中使用 KEGGREST 或 msigdbr 包获取
pathways = {
    "hsa00010_Glycolysis": [f"GENE_{i}" for i in range(0, 30)],
    "hsa04151_PI3K_Akt": [f"GENE_{i}" for i in range(20, 70)],
    "hsa04010_MAPK": [f"GENE_{i}" for i in range(50, 110)],
    "hsa00190_Oxidative_phosphorylation": [f"GENE_{i}" for i in range(80, 130)],
    "hsa04110_Cell_cycle": [f"GENE_{i}" for i in range(100, 160)],
}

# ========== 3. 超几何检验（ORA） ==========
def ora_test(query_genes, pathway_genes, total_genes=20000):
    """
    超几何检验：检验query_genes在pathway_genes中的富集情况
    就像从牌堆中随机抽牌，统计抽到红心的概率
    """
    k = len(set(query_genes) & set(pathway_genes))  # 交集数量（抽到红心数）
    M = total_genes     # 总基因数（牌堆总数）
    n = len(pathway_genes)  # 通路基因数（红心总数）
    N = len(query_genes)    # 查询基因数（抽牌数）

    if k == 0:
        return 1.0, k  # 没有交集，不显著

    # 超几何检验：计算至少k个基因在通路中的概率
    p_value = stats.hypergeom.sf(k - 1, M, n, N)
    return p_value, k

# 对每条通路做ORA
print("\n===== RNA组学通路分析 =====")
rna_enrichment = []
for pathway_name, pathway_genes in pathways.items():
    p_val, overlap = ora_test(sig_rna, pathway_genes)
    rna_enrichment.append({
        "pathway": pathway_name,
        "overlap": overlap,
        "pathway_size": len(pathway_genes),
        "p_value": p_val,
        "odds_ratio": (overlap / len(sig_rna)) / (len(pathway_genes) / 20000)
    })

rna_enrich_df = pd.DataFrame(rna_enrichment).sort_values("p_value")

# 多重检验校正（Benjamini-Hochberg方法）
from scipy.stats import rankdata
p_values = rna_enrich_df["p_value"].values
n_tests = len(p_values)
ranks = rankdata(p_values)  # 从小到大排序获得名次
fdr = p_values * n_tests / ranks  # BH校正
fdr = np.minimum(fdr, 1.0)  # FDR不超过1
rna_enrich_df["FDR"] = fdr
print(rna_enrich_df[["pathway", "overlap", "p_value", "FDR"]].to_string(index=False))

# ========== 4. 多组学联合富集 ==========
print("\n===== 多组学联合通路分析 =====")
# 策略：要求通路在RNA和蛋白两个层次都富集
protein_enrichment = []
for pathway_name, pathway_genes in pathways.items():
    p_val, overlap = ora_test(sig_protein, pathway_genes)
    protein_enrichment.append({
        "pathway": pathway_name,
        "p_value_protein": p_val,
        "overlap_protein": overlap
    })

protein_enrich_df = pd.DataFrame(protein_enrichment)

# 合并两组学的富集结果
joint_enrichment = rna_enrich_df.merge(
    protein_enrich_df, on="pathway"
)

# 联合显著性：两个p值的乘积（Fisher合并法）
joint_enrichment["joint_p"] = (
    joint_enrichment["p_value"] * joint_enrichment["p_value_protein"]
)
joint_enrichment = joint_enrichment.sort_values("joint_p")

print("\n多组学联合富集结果（按联合p值排序）：")
print(joint_enrichment[["pathway", "overlap", "overlap_protein", "joint_p"]].to_string(index=False))

# ========== 5. Pathview可视化（R语言） ==========
r_pathview_code = """
library(pathview)  # 在KEGG通路图上叠加表达量数据

# 准备基因数据（Entrez ID为行名，log2FC为值）
gene.data <- setNames(rna_results$log2FC, rna_results$entrez_id)

# 准备代谢物数据（KEGG化合物ID为行名）
cpd.data <- setNames(metabolite_results$log2FC, metabolite_results$cpd_id)

# 绘制Glycolysis通路图，颜色代表上调（红）/下调（绿）
pathview(
    gene.data = gene.data,
    cpd.data = cpd.data,
    pathway.id = "hsa00010",  # KEGG通路ID
    species = "hsa",           # 人类
    out.suffix = "multiomics", # 输出文件后缀
    kegg.native = TRUE         # 使用KEGG原始通路图
)
"""
print("\nR-Pathview代码（需R环境运行）:")
print(r_pathview_code)