通路富集分析秒出图 - AI科研技能库

🎯

技能简介

差异基因分析完成后，下一步就是回答"这些基因参与什么生物学通路？"。传统通路富集分析流程繁琐：上传基因列表 → 等待 → 下载结果 → 手动画图， 一套流程下来，半天时间没了。

本教程教你使用三大通路数据库—— KEGG（京都基因与基因组百科全书）、Reactome（反应组通路数据库）、STRING（蛋白质互作网络），通过Scientific Skills MCP自动化完成从差异基因到发表级气泡图的全流程。 将分析时间从半天缩短至10分钟！

半天

传统富集分析

10分钟

AI辅助分析

30倍

效率提升

💡 使用场景

🧬

差异基因功能解读

RNA-seq或芯片数据分析后，快速了解差异基因参与的生物学通路

🔍

疾病机制研究

发现疾病相关的异常通路，为疾病机制研究提供线索

💊

药物作用机制

分析药物处理后差异基因的通路变化，揭示药物作用机制

🕸️

关键调控因子发现

通过蛋白质互作网络发现hub基因，识别关键调控因子

🛠️ 核心技能调用

通过Scientific Skills MCP，可以一键调用三大通路分析工具：

🔗 KEGG Database - 通路富集分析

KEGG是最常用的通路富集数据库，涵盖代谢通路、信号通路、疾病通路等。支持气泡图、条形图等多种可视化方式。

# 调用 KEGG Database 技能
# scientific-skills:kegg-database

# 示例：差异基因KEGG通路富集
kegg_enrichment(
    gene_list = ["TP53", "EGFR", "MYC", "KRAS", "PIK3CA", # ... 差异基因列表],
    organism = "hsa",  # 人类
    pvalue_cutoff = 0.05,
    fdr_cutoff = 0.05,
    top_n = 20  # 返回top 20通路
)

# 返回结果包含：
# - 通路ID和名称
# - 基因数量（输入基因/背景基因）
# - p值和FDR（Benjamini-Hochberg校正）
# - 富集因子
# - 涉及的基因列表

# 自动生成气泡图
plot_enrichment_bubble(
    kegg_results,
    x = "GeneRatio",  # x轴：基因比例
    y = "Pathway",     # y轴：通路名称
    size = "Count",     # 气泡大小：基因数量
    color = "p.adjust"  # 气泡颜色：FDR值
)

💡 KEGG气泡图解读：

• 气泡大小：该通路中差异基因的数量
• 气泡颜色：FDR值，红色越深表示显著性越高
• 横轴位置：基因比例，反映富集程度

⚡ Reactome Database - 反应组通路分析

Reactome是高质量的反应组通路数据库，包含细胞信号、代谢、转录调控等通路。支持层级结构可视化。

# 调用 Reactome Database 技能
# scientific-skills:reactome-database

# Reactome通路富集分析
reactome_enrichment(
    gene_list = ["TP53", "EGFR", "MYC", "KRAS", "PIK3CA"],
    species = "Homo sapiens",
    include_low_levels = False,  # 只返回高层级通路
    pvalue_cutoff = 0.05
)

# 返回结果包含：
# - 通路ID、名称、层级
# - p值和FDR
# - 富集统计信息
# - 通路层级结构

# 生成条形图
plot_enrichment_bar(
    reactome_results,
    top_n = 15,
    color_by = "p.adjust",
    show_counts = True
)

📌 KEGG vs Reactome：

• KEGG：覆盖更广，包含代谢通路、疾病通路
• Reactome：更精细的反应级描述，适合信号转导研究
• 推荐：两个数据库同时使用，交叉验证结果

🕸️ STRING Database - 蛋白质互作网络

STRING数据库包含已知和预测的蛋白质互作关系，可以构建PPI网络并识别hub基因。

# 调用 STRING Database 技能
# scientific-skills:string-database

# 构建蛋白质互作网络
string_network(
    gene_list = ["TP53", "EGFR", "MYC", "KRAS", "PIK3CA"],
    species = 9606,  # 人类
    confidence_threshold = 0.7,  # 互作置信度阈值
    add_nodes = 50  # 添加额外互作蛋白
)

# 返回结果包含：
# - 蛋白质互作网络（节点和边）
# - 互作分数（综合分数）
# - 互作证据类型（实验、数据库、共表达等）

# 识别hub基因
hub_genes = identify_hub_genes(
    string_network,
    method = "betweenness",  # 或"degree", "closeness"
    top_n = 10
)

# 可视化网络
plot_network(
    string_network,
    layout = "fruchterman_reingold",
    highlight_hubs = True,
    node_size_by = "degree",
    label_hubs = True
)

💡 Hub基因识别方法：

• Degree：连接数最多的节点
• Betweenness：最短路径经过次数最多的节点
• Closeness：到所有其他节点平均距离最短的节点
• Hub基因通常是关键调控因子，是进一步研究的候选靶点

📖 实战示例：癌症差异基因通路分析

以下是完整的差异基因通路富集分析流程，从KEGG/Reactome富集到STRING网络构建：

1 数据准备

准备差异基因列表（通常来自RNA-seq差异分析）：

# 示例差异基因列表（假设来自癌症vs正常组织对比）
differential_genes = [
    "TP53", "EGFR", "MYC", "KRAS", "PIK3CA",
    "PTEN", "AKT1", "MTOR", "CDK4", "CCND1",
    "RB1", "MDM2", "BCL2", "BAX", "CASP3",
    # ... 总共200个差异基因
]

# 基因ID转换（如需要）
convert_gene_ids(
    gene_list = differential_genes,
    from_type = "symbol",  # 基因符号
    to_type = "entrez",   # Entrez ID
    species = "human"
)

# 过滤背景基因（用于富集分析）
background_genes = get_background_genes(
    species = "human",
    gene_type = "protein_coding"
)

2 KEGG通路富集

使用KEGG数据库进行通路富集分析：

# Step 2: KEGG富集
kegg_results = kegg_enrichment(
    gene_list = differential_genes,
    background = background_genes,
    organism = "hsa",
    pvalue_cutoff = 0.05,
    fdr_method = "BH"
)

# 查看显著富集的通路
print(kegg_results[head(10)])

# 示例输出：
# ID          Description                  GeneRatio  BgRatio  pvalue  p.adjust
# hsa05200    Pathways in cancer           45/200    500/20000  1.2e-15  3.4e-12
# hsa04151    PI3K-Akt signaling pathway  28/200    350/20000  5.6e-12  8.9e-09
# hsa05206    MicroRNAs in cancer         32/200    400/20000  2.1e-10  2.7e-07
# hsa04068    FoxO signaling pathway      18/200    200/20000  4.5e-08  4.8e-05

# 生成气泡图
plot_enrichment_bubble(
    kegg_results,
    title = "KEGG通路富集分析",
    x = "GeneRatio",
    size = "Count",
    color = "p.adjust"
)

3 Reactome通路富集

使用Reactome数据库进行交叉验证：

# Step 3: Reactome富集
reactome_results = reactome_enrichment(
    gene_list = differential_genes,
    species = "Homo sapiens",
    pvalue_cutoff = 0.05
)

# 对比KEGG和Reactome结果
common_pathways = find_common_pathways(
    kegg_results,
    reactome_results
)

print(f"共同富集的通路数: {len(common_pathways)}")

# 示例共同通路：
# - PI3K-Akt signaling (KEGG) / Signaling by PI3K (Reactome)
# - Cell cycle (KEGG) / Cell Cycle (Reactome)
# - p53 signaling (KEGG) / p53 dependent DNA damage (Reactome)

# 生成对比条形图
plot_comparison_bar(
    kegg_results,
    reactome_results,
    top_n = 15
)

4 STRING蛋白互作网络

构建蛋白质互作网络并识别hub基因：

# Step 4: STRING网络构建
string_net = string_network(
    gene_list = differential_genes,
    species = 9606,
    confidence_threshold = 0.7,
    network_type = "full"  # 包含所有证据类型
)

# 计算网络拓扑参数
network_stats = calculate_network_stats(string_net)
print(network_stats)

# 示例输出：
# - 节点数: 185
# - 边数: 1247
# - 平均度: 13.5
# - 网络密度: 0.073
# - 聚类系数: 0.42

# 识别hub基因（使用多种指标）
hub_genes = identify_hub_genes(
    string_net,
    methods = ["degree", "betweenness", "closeness"],
    top_n = 10
)

# Top hub基因示例：
# Rank  Gene    Degree  Betweenness  Closeness
# 1     TP53    67      0.234        0.512
# 2     MYC     54      0.187        0.489
# 3     EGFR    48      0.156        0.467
# 4     AKT1    45      0.142        0.451
# 5     PIK3CA  42      0.128        0.438

# 可视化网络
plot_network(
    string_net,
    highlight_genes = hub_genes['Gene'].head(5),
    node_color = 'betweenness',
    node_size = 'degree',
    title = "差异基因蛋白互作网络"
)

5 整合可视化与报告

生成整合的多图报告：

# Step 5: 整合可视化
generate_enrichment_report(
    kegg_results = kegg_results,
    reactome_results = reactome_results,
    string_network = string_net,
    hub_genes = hub_genes,
    output_format = "pdf"
)

# 报告包含：
# 1. KEGG气泡图 (Top 20通路)
# 2. Reactome条形图 (Top 15通路)
# 3. 共同富集通路对比图
# 4. STRING蛋白互作网络图
# 5. Hub基因统计表
# 6. 通路-基因对应表

# 导出数据表格（用于论文补充材料）
export_tables(
    kegg_results,
    file = "kegg_enrichment.xlsx"
)
export_tables(
    hub_genes,
    file = "hub_genes.xlsx"
)