代谢通路富集分析（KEGG/MetPA）¶

# === 使用Python进行代谢通路富集分析 ===
import pandas as pd
from scipy.stats import fisher_exact
import numpy as np

# 1. 准备数据
# 差异代谢物列表（已映射到KEGG compound ID）
diff_metabolites = ["C00031", "C00036", "C00042", "C00074", "C00149"]  # 差异代谢物KEGG ID
# 背景代谢物列表（所有检测到的代谢物）
all_metabolites = ["C00031", "C00036", "C00042", "C00074", "C00149",
                   "C00022", "C00024", "C00158", "C00311", "C00417",
                   "C00026", "C00122", "C00091", "C00044", "C00035"]

# 2. 定义KEGG通路及其包含的代谢物
pathways = {
    "TCA cycle": ["C00036", "C00042", "C00074", "C00149", "C00158", "C00311", "C00026", "C00122", "C00091"],
    "Glycolysis": ["C00031", "C00022", "C00024", "C00074", "C00036"],
    "Amino acid metabolism": ["C00044", "C00035", "C00417"]
}

# 3. 对每条通路做Fisher精确检验
results = []
n_total = len(all_metabolites)  # 背景代谢物总数
n_diff = len(diff_metabolites)   # 差异代谢物总数

for pathway_name, pathway_mets in pathways.items():
    # 构建2x2列联表
    pathway_in_bg = len(set(pathway_mets) & set(all_metabolites))  # 通路中属于背景的
    diff_in_pathway = len(set(diff_metabolites) & set(pathway_mets))  # 差异且在通路中
    diff_not_pathway = n_diff - diff_in_pathway  # 差异但不在通路中
    not_diff_in_pathway = pathway_in_bg - diff_in_pathway  # 不差异但在通路中
    not_diff_not_pathway = n_total - n_diff - not_diff_in_pathway  # 不差异也不在通路中

    table = [[diff_in_pathway, diff_not_pathway],
             [not_diff_in_pathway, not_diff_not_pathway]]
    odds_ratio, p_value = fisher_exact(table, alternative="greater")  # Fisher精确检验

    results.append({
        "Pathway": pathway_name,
        "Hits": diff_in_pathway,
        "Total": pathway_in_bg,
        "P_value": p_value
    })

# 4. 输出结果
df_results = pd.DataFrame(results).sort_values("P_value")  # 按p值排序
print(df_results.to_string(index=False))

# === 使用MetaboAnalystR做通路分析 ===
library(MetaboAnalystR)  # MetaboAnalyst的R版本

# 1. 初始化分析对象
mSet <- InitDataObjects("conc", "pathora", FALSE)  # 浓度数据，通路ORA分析

# 2. 设置代谢物列表
cmpd_list <- c("Glucose", "Pyruvate", "Citrate", "Succinate", "Fumarate")
mSet <- Setup.MapData(mSet, cmpd_list)  # 设置代谢物
mSet <- CrossReferencing(mSet, "name")  # 名称匹配到KEGG ID

# 3. 运行通路分析
mSet <- SetKEGG.PathLib(mSet, "hsa")  # 选择人类KEGG通路
mSet <- SetMetabolomeFilter(mSet, FALSE)  # 不过滤
mSet <- CalculateOraScore(mSet, "rbc", "hyperg")  # ORA + 超几何检验