Kegg Pathway Analysis

Common Pitfalls

1. Combining up-regulated and down-regulated genes into a single enrichment run. Pathways with genes regulated in opposite directions cancel out, producing fewer significant pathways than the true count. Results from combined lists are unreliable. Example of the anti-pattern:

   # 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

   • How to avoid: Always split DEGs by direction before ORA. Run enrichment twice (once for up, once for down) and take the union of significant pathways.

2. Using the wrong KEGG organism code. KEGG silently returns empty results for invalid or mismatched organism codes. This is especially common for bacterial species with multiple strain-specific codes (e.g., pae for PAO1 vs pau for PA14).

   • How to avoid: Confirm the organism code from the KEGG organism list before running enrichment. For bacteria, verify the specific strain code.

3. Gene ID type mismatch. Providing gene symbols when KEGG expects Entrez IDs (or locus tags for bacteria) silently yields zero enriched pathways with no error message.

   • How to avoid: Check the expected ID type for your organism. Use clusterProfiler::bitr() or equivalent to convert gene symbols to Entrez IDs before enrichment.

4. Not handling KEGG API timeouts. The KEGG REST API frequently times out, causing enrichKEGG() to fail mid-analysis. Without error handling, the entire analysis is lost.

   • How to avoid: Wrap KEGG API calls in retry logic with short timeouts (30 seconds). Pre-fetch pathway data at the start. Have gseapy as a fallback.

5. Delaying the answer to resolve pathway names. Calling keggGet() to convert pathway IDs to human-readable names after computing results can timeout, losing the numeric answer entirely. Example of the anti-pattern:

   # 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)

   • How to avoid: Report pathway IDs and counts immediately. Resolve names only after the primary result is secured, inside a tryCatch() or try/except block.

6. Omitting the background universe. Not specifying the universe parameter defaults to the entire KEGG gene database for that organism, inflating statistical significance for pathways containing well-annotated housekeeping genes.

   • How to avoid: Always pass universe = rownames(res) (all tested genes from the DE analysis) to enrichKEGG().

7. Using raw p-values instead of adjusted p-values for filtering. Reporting pathways with p < 0.05 without multiple testing correction dramatically increases false positives.

   • How to avoid: Always filter on p.adjust < 0.05. Verify that the column you are filtering is the adjusted value, not the raw p-value.

Workflow

1. Step 1: Prepare gene lists

   • Filter significant DEGs from differential expression results (e.g., padj < 0.05, |log2FC| > 1.5)
   • Split into up-regulated and down-regulated gene sets
   • Convert gene IDs to the format expected by KEGG (Entrez IDs or locus tags)

   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. Step 2: Pre-fetch KEGG data (recommended)

   • Cache pathway-gene mappings locally before running enrichment
   • Decision point: If KEGG API is accessible, proceed with live queries. If not, switch to offline gene sets (Step 2b).

   R (KEGGREST package):

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

   Python (requests with caching):

   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. Step 3: Run enrichment separately for each direction

   • Run ORA twice: once for up-regulated genes, once for down-regulated genes
   • Always specify the universe (all tested genes)

   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. Step 4: Report results immediately

   • Emit pathway IDs and counts before attempting name resolution
   • Optionally resolve pathway names in a non-blocking manner

   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. Step 5: Cross-condition comparison (if applicable)

   • Run separate up/down enrichment for each condition (4 enrichKEGG calls total)
   • Compare pathway sets within the same direction
   • Report the union of direction-specific unique pathways

   # 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. Step 6: Handle API failures with retry logic

   • If KEGG API calls fail, retry with timeout before falling back to offline data

   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)
   }

Further Reading

• KEGG: Kyoto Encyclopedia of Genes and Genomes -- Official KEGG database with pathway maps, organism lists, and REST API documentation
• clusterProfiler: an R package for comparing biological themes among gene clusters -- Original publication by Guangchuang Yu describing the enrichKEGG and gseKEGG functions
• gseapy documentation -- Python package for gene set enrichment analysis with offline KEGG gene set support
• KEGG REST API reference -- Technical documentation for programmatic access to KEGG pathway data

Related Skills

• gseapy-gene-enrichment -- Python-based gene set enrichment analysis; use as a fallback when clusterProfiler KEGG API calls fail, or as the primary tool for Python-based workflows
• deseq2-differential-expression / pydeseq2-differential-expression -- Upstream differential expression analysis that produces the DEG lists used as input to KEGG pathway enrichment