kegg-pathway-analysis

KEGG 通路富集分析指南

概述

KEGG（京都基因与基因组百科全书）通路富集分析用于识别在差异表达基因中统计学上显著过度代表的生物通路。本指南介绍两种主要的富集方法（ORA 和 GSEA）、关键工作流决策（例如按变化方向拆分基因）、R clusterProfiler 与 Python gseapy 之间的工具选择，以及应对以不可靠著称的 KEGG REST API 的策略。指南还讨论了导致通路计数错误或分析停滞的常见故障模式。

KEGG 通路分析中最常见的三种错误是：（1）将上调和下调基因合并到一次富集分析中，从而掩盖真实的通路信号；（2）KEGG REST API 超时且没有备用策略，导致分析失败；（3）尝试查询仅用于美化显示的通路名称时，迟迟不报告结果，而此类查询可能永远无法完成。本指南为每种情况提供了具体解决方案。

关键概念

ORA 与 GSEA

过度代表分析（ORA）和基因集富集分析（GSEA）是两种主要的通路富集方法，它们在输入和统计方法上都各不相同。

ORA 接收经过预筛选的基因列表（例如 padj < 0.05 且 |log2FC| > 1.5 的基因），并检验相对于背景基因总体，KEGG 通路成员在该列表中是否过度代表。ORA 使用超几何检验（Fisher 精确检验）。这种方法简单直接，但会丢弃变化幅度信息，并且高度依赖所选的显著性阈值。

GSEA 接收所有基因的排序列表（通常按 log2 倍数变化或带符号的显著性统计量排序），不设置任何阈值。它通过遍历排序列表计算累积富集分数，并识别成员集中在排序顶部或底部的通路。GSEA 能捕捉 ORA 可能遗漏的细微协同变化。

实际应用中，通常会从通过 enrichKEGG()（clusterProfiler）或 gp.enrichr()（gseapy）进行 ORA 开始。如果希望避免任意设定阈值，或者效应量较小，则更适合通过 gseKEGG() 或 gp.prerank() 进行 GSEA。

富集分析中的变化方向

进行 ORA 时，基因的变化方向——即基因是上调还是下调——至关重要。单个通路可能包含变化方向相反的基因。如果将上调和下调基因合并到一个列表中，它们相反的信号会相互抵消，削弱富集信号并掩盖真正富集的通路。分别对上调和下调基因集运行富集分析，可以得到更准确且更易解释的结果。对于 ORA，这种拆分是必须的。GSEA 通过带符号的排序天然处理变化方向，不过仍需按方向解读前沿基因。

KEGG 生物代码

KEGG 使用三字母（或四字母）生物代码来标识物种特异性通路数据库。使用错误代码时，结果会悄然为空。常见代码：

生物	代码
人类	hsa
小鼠	mmu
大鼠	rno
斑马鱼	dre
果蝇	dme
秀丽隐杆线虫	cel
大肠杆菌 K-12	eco
铜绿假单胞菌 PA14	pau
铜绿假单胞菌 PAO1	pae
酿酒酵母	sce
拟南芥	ath

基因 ID 格式也因生物而异：真核生物通常需要 Entrez 基因 ID，而细菌则使用基因座标签。ID 类型不匹配是一种不易察觉的失败模式。

KEGG API 可靠性

KEGG REST API（rest.kegg.jp）受到速率限制，经常响应缓慢，并且容易超时。clusterProfiler::enrichKEGG() 和直接向 KEGG 发送 HTTP 请求都可能不可预测地失败。为 API 故障做好准备并非可选项，而是任何基于 KEGG 的工作流中的必要环节。可采取的策略包括预先获取并缓存通路数据、使用 gseapy 附带的离线基因集数据库，以及实现带超时设置的重试逻辑。

前置访谈

在编写任何分析代码之前，先与用户确认以下事项。

decisions:
  - id: D1
    param: organismCode
    kind: required
    source: upstream
    ask: "Which species should pathways be looked up for?"
    default: "carried from the dataset"

  - id: D2
    param: identifierMapping
    kind: required
    source: data
    depends_on: [D1]
    ask: "Which identifier type do the genes carry, and how should unmapped ones be handled?"
    default: null

  - id: D3
    param: testingApproach
    kind: required
    source: upstream
    ask: "Test a significant-gene list against pathways (over-representation), or the whole ranked list without a cut (enrichment)?"
    default: "over-representation when the input is a gene list"

  - id: D4
    param: backgroundSet
    kind: required
    source: upstream
    depends_on: [D3]
    ask: "Which genes form the background - everything measured in this experiment, or the whole annotated genome?"
    default: "the genes tested in this experiment"
    skip_if: "ranked-list enrichment, which uses the ranking rather than a background"

  - id: D5
    param: pathwayScope
    kind: optional
    source: user
    ask: "Should all pathway categories be tested, or only metabolism, signalling, or disease maps?"
    default: "all categories"

  - id: D6
    param: significanceThreshold
    kind: required
    source: user
    ask: "How strong must a pathway's evidence be, after correcting for the number of pathways tested?"
    default: "adjusted p-value 0.05"

D4 是最常被忽略、也最能改变答案的决策。 如果实验只测量了表达基因，却将整个基因组作为背景，那么每个组织特异性通路看起来都会显著富集——富集是相对于那些根本无法测量的基因，而非相对于实验本身。

决策框架

Question: What enrichment analysis do you need?
|
+-- Have a pre-filtered DEG list (with cutoffs applied)?
|   +-- Yes --> ORA
|   |   +-- Using R? --> clusterProfiler::enrichKEGG()
|   |   +-- Using Python? --> gseapy.enrichr()
|   |   +-- KEGG API failing? --> gseapy with offline gene sets
|   +-- No, want cutoff-free analysis --> GSEA
|       +-- Using R? --> clusterProfiler::gseKEGG()
|       +-- Using Python? --> gseapy.prerank()
|
+-- Need to split by direction?
|   +-- ORA --> YES, always split up/down (mandatory)
|   +-- GSEA --> No split needed (direction encoded in ranking)
|
+-- KEGG API unreliable?
    +-- Try cached/pre-fetched data first
    +-- Fall back to gseapy offline databases
    +-- Use retry logic with short timeouts

场景	推荐方法	理由
使用 R 进行标准 ORA	clusterProfiler::enrichKEGG()，按方向拆分	使用最广泛，与 Bioconductor 生态系统集成
使用 Python 进行标准 ORA	使用 KEGG_2021_Human 的 gseapy.enrichr()	离线基因集无需依赖 API
无阈值富集	通过 gseKEGG() 或 gp.prerank() 进行 GSEA	能捕捉细微的协同变化，无需任意阈值
KEGG API 不可用	切换到 gseapy 离线数据库	gseapy 在本地捆绑了 KEGG 基因集
比较不同条件	对每个条件分别进行上调和下调富集分析	支持跨条件进行考虑方向的集合运算
非模式生物	核实生物体代码，并使用 KEGGREST 检查可用性	错误代码会静默返回空结果

最佳实践

1. 始终按基因变化方向拆分 ORA。 分别对上调和下调基因运行 enrichKEGG() 或 gp.enrichr()。将它们合并会扩大基因列表、稀释富集信号，并产生错误的通路计数。报告两个方向显著通路的并集。

2. 明确指定背景基因集合。 将背景设为所有经过测试的基因（即差异表达分析中的完整基因集合），而非仅包含显著基因。省略背景集合时，默认会使用 KEGG 数据库中的所有基因，这会夸大研究充分的通路的显著性。

3. 运行富集分析前，预先获取并缓存 KEGG 数据。 在分析开始时下载通路与基因的映射关系并保存到本地。这样可以避免 KEGG API 无响应时分析中途失败，并使分析具有可复现性。

4. 解析通路名称之前，先报告数值结果。 计算出显著通路 ID 的数量或列表后，立即输出该结果。通过额外的 KEGG API 调用将 ID 解析为易读名称只是表面处理，调用可能超时，导致主要结果丢失。

5. 始终一致地进行多重检验校正。 使用校正后的 p 值（通常采用 BH 方法，p.adjust < 0.05），而非原始 p 值。clusterProfiler 和 gseapy 默认都会进行校正，但务必确认使用的是校正后的值作为阈值。

6. 核实基因 ID 格式与生物体相匹配。 真核生物的 KEGG 通路需要 Entrez 基因 ID；细菌则需要位点标签。格式不匹配会静默返回零条富集通路。如果输入使用基因符号，请使用 clusterProfiler 中的 bitr() 或等效方法转换 ID。

7. clusterProfiler 失败时，使用 gseapy 作为备用方案。 当 enrichKEGG() 因 KEGG API 问题失败时，gseapy 的 enrichr() 函数搭配捆绑的离线基因集（例如 KEGG_2021_Human），无需任何网络依赖即可提供等效的 ORA 结果。

常见陷阱

1. 将上调和下调基因合并到同一次富集分析中。 受调控方向相反的基因会相互抵消，导致显著通路数量少于实际数量。合并列表得出的结果不可靠。反面示例如下：

   # WRONG: combining up and down genes into one list
   all_sig_genes <- rownames(subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5))
   ekegg <- enrichKEGG(gene = all_sig_genes, ...)  # Will miss pathways
   

   • 避免方法：进行 ORA 前，务必按调控方向拆分 DEG。分别运行两次富集分析（一次针对上调基因，一次针对下调基因），然后取显著通路的并集。

2. 使用了错误的 KEGG 生物体代码。 对于无效或不匹配的生物体代码，KEGG 会悄无声息地返回空结果。细菌物种有多个菌株专属代码时，这种情况尤其常见（例如，PAO1 对应 pae，PA14 对应 pau）。

   • 避免方法：运行富集分析前，先从 KEGG 生物体列表中确认代码。对于细菌，请核实具体菌株的代码。

3. 基因 ID 类型不匹配。 当 KEGG 需要 Entrez ID（细菌则需要位点标签）时，如果提供的是基因符号，结果会悄无声息地显示富集通路数量为零，且不会报错。

   • 避免方法：检查你的生物体所需的 ID 类型。使用 clusterProfiler::bitr() 或同等工具，在富集分析前将基因符号转换为 Entrez ID。

4. 未处理 KEGG API 超时。 KEGG REST API 经常超时，导致 enrichKEGG() 在分析过程中失败。如果没有错误处理，整个分析结果都会丢失。

   • 避免方法：使用带有短超时时间（30 秒）的重试逻辑封装 KEGG API 调用。在开始时预先获取通路数据。准备好 gseapy 作为备用方案。

5. 为了解析通路名称而延迟返回答案。 结果计算完成后调用 keggGet() 将通路 ID 转换为易读名称，可能会超时，导致数值结果也一并丢失。反面示例如下：

   # BAD: answer delayed by name lookup that may hang
   result_ids <- setdiff(pathways_condA, pathways_condB)
   names <- keggGet(result_ids)  # Can timeout -- answer never emitted
   count <- length(result_ids)
   

   • 避免方法：立即报告通路 ID 和数量。只有在主要结果得到妥善保存后，才在 tryCatch() 或 try/except 块中解析名称。

6. 未指定背景基因集。 未指定 universe 参数时，默认使用该生物体的整个 KEGG 基因数据库，这会夸大含有注释完善的管家基因的通路的统计显著性。

   • 避免方法：始终向 enrichKEGG() 传入 universe = rownames(res)（差异表达分析中所有经过检验的基因）。

7. 筛选时使用原始 p 值，而非校正后的 p 值。 未经多重检验校正便报告 p < 0.05 的通路，会大幅增加假阳性。

   • 避免方法：始终根据 p.adjust < 0.05 进行筛选。确认筛选的列是校正后的值，而非原始 p 值。

工作流

1. 步骤 1：准备基因列表

   • 从差异表达结果中筛选显著差异表达基因（DEG）（例如，padj < 0.05、|log2FC| > 1.5）
   • 将基因拆分为上调和下调基因集
   • 将基因 ID 转换为 KEGG 所需的格式（Entrez ID 或基因座标签）

   library(clusterProfiler)
   
   # Filter significant DEGs
   sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5)
   
   # MANDATORY: Split by direction BEFORE running enrichment
   up_genes <- rownames(subset(sig_genes, log2FoldChange > 0))
   dn_genes <- rownames(subset(sig_genes, log2FoldChange < 0))
   
   cat("Up-regulated genes:", length(up_genes), "\n")
   cat("Down-regulated genes:", length(dn_genes), "\n")
   

2. 步骤 2：预先获取 KEGG 数据（推荐）

   • 在运行富集分析前，将通路-基因映射缓存到本地
   • 决策点：如果可以访问 KEGG API，则继续进行实时查询；否则，切换到离线基因集（步骤 2b）。

   R（KEGGREST 包）：

   library(KEGGREST)
   pathway_list <- tryCatch(
     keggList("pathway", organism_code),
     error = function(e) NULL
   )
   

   Python（使用 requests 和缓存）：

   import requests
   import json
   import os
   
   cache_file = f"kegg_{organism_code}_pathways.json"
   if os.path.exists(cache_file):
       with open(cache_file) as f:
           pathway_map = json.load(f)
   else:
       resp = requests.get(
           f"https://rest.kegg.jp/list/pathway/{organism_code}", timeout=30
       )
       if resp.ok:
           pathway_map = dict(
               line.split("\t") for line in resp.text.strip().split("\n")
           )
           with open(cache_file, 'w') as f:
               json.dump(pathway_map, f)
   

3. 步骤 3：针对每个方向分别运行富集分析

   • 运行两次 ORA：一次分析上调基因，一次分析下调基因
   • 始终指定背景基因集（所有经过检验的基因）

   R（clusterProfiler）：

   ekegg_up <- enrichKEGG(gene = up_genes, organism = organism_code,
                          universe = rownames(res), pvalueCutoff = 0.05)
   ekegg_dn <- enrichKEGG(gene = dn_genes, organism = organism_code,
                          universe = rownames(res), pvalueCutoff = 0.05)
   
   up_pathways <- subset(as.data.frame(ekegg_up), p.adjust < 0.05)$ID
   dn_pathways <- subset(as.data.frame(ekegg_dn), p.adjust < 0.05)$ID
   
   all_sig_pathways <- union(up_pathways, dn_pathways)
   cat("Significant pathways (up):", length(up_pathways), "\n")
   cat("Significant pathways (down):", length(dn_pathways), "\n")
   cat("Total unique significant pathways:", length(all_sig_pathways), "\n")
   

   Python（gseapy）：

   import gseapy as gp
   
   up_genes = sig_genes[sig_genes['log2FoldChange'] > 0].index.tolist()
   dn_genes = sig_genes[sig_genes['log2FoldChange'] < 0].index.tolist()
   
   enr_up = gp.enrichr(gene_list=up_genes, gene_sets='KEGG_2021_Human',
                        organism='human', outdir=None)
   enr_dn = gp.enrichr(gene_list=dn_genes, gene_sets='KEGG_2021_Human',
                        organism='human', outdir=None)
   

4. 步骤 4：立即报告结果

   • 在尝试解析通路名称前，先输出通路 ID 和数量
   • 可以采用非阻塞方式解析通路名称

   result_ids <- setdiff(pathways_condA, pathways_condB)
   count <- length(result_ids)
   cat("Answer:", count, "pathways unique to condition A\n")
   cat("Pathway IDs:", paste(result_ids, collapse = ", "), "\n")
   
   # Optional: resolve names (non-blocking)
   names <- tryCatch(keggGet(result_ids), error = function(e) NULL)
   

5. 步骤 5：跨条件比较（如适用）

   • 针对每个条件分别运行上调/下调富集分析（共调用 4 次 enrichKEGG）
   • 在相同方向内比较通路集
   • 报告各方向特有通路的并集

   # Condition A
   up_pathways_A <- subset(as.data.frame(ekegg_up_A), p.adjust < 0.05)$ID
   dn_pathways_A <- subset(as.data.frame(ekegg_dn_A), p.adjust < 0.05)$ID
   
   # Condition B
   up_pathways_B <- subset(as.data.frame(ekegg_up_B), p.adjust < 0.05)$ID
   dn_pathways_B <- subset(as.data.frame(ekegg_dn_B), p.adjust < 0.05)$ID
   
   # Find pathways unique to condition A in each direction
   unique_up <- setdiff(up_pathways_A, up_pathways_B)
   unique_dn <- setdiff(dn_pathways_A, dn_pathways_B)
   
   # Union
   unique_to_A <- union(unique_up, unique_dn)
   cat("Pathways in A but not B:", length(unique_to_A), "\n")
   cat("  From up-regulated:", length(unique_up), "-",
       paste(unique_up, collapse = ", "), "\n")
   cat("  From down-regulated:", length(unique_dn), "-",
       paste(unique_dn, collapse = ", "), "\n")
   

6. 步骤 6：通过重试逻辑处理 API 故障

   • 如果 KEGG API 调用失败，请先设置超时后重试，再回退到离线数据

   fetch_kegg_with_retry <- function(organism_code, max_retries = 3,
                                      timeout_sec = 30) {
     for (i in seq_len(max_retries)) {
       result <- tryCatch({
         R.utils::withTimeout(
           keggList("pathway", organism_code),
           timeout = timeout_sec
         )
       }, error = function(e) NULL)
       if (!is.null(result)) return(result)
       Sys.sleep(2)
     }
     warning("KEGG API unreachable after retries. Proceeding without pathway names.")
     return(NULL)
   }
   

延伸阅读

• KEGG：京都基因与基因组百科全书 -- KEGG 官方数据库，提供通路图、物种列表和 REST API 文档
• clusterProfiler：用于比较基因簇间生物学主题的 R 包 -- Guangchuang Yu 介绍 enrichKEGG 和 gseKEGG 函数的原始论文
• gseapy 文档 -- 用于基因集富集分析的 Python 包，支持离线 KEGG 基因集
• KEGG REST API 参考 -- 以编程方式访问 KEGG 通路数据的技术文档

相关技能

• gseapy-gene-enrichment -- 基于 Python 的基因集富集分析；当 clusterProfiler 的 KEGG API 调用失败时可用作后备方案，也可作为基于 Python 的工作流的主要工具
• deseq2-differential-expression / pydeseq2-differential-expression -- 上游差异表达分析，用于生成作为 KEGG 通路富集分析输入的 DEG 列表