Bio Core Pathways

Gene Ontology and Pathway Analysis

When to Use

• Testing whether a DE gene list (RNA-seq, proteomics) is over-represented in GO terms or KEGG/Reactome pathways.
• Running GSEA on a full ranked gene list (log fold-change or t-statistic) instead of an arbitrary significance cutoff.
• Querying the KEGG REST API for pathway membership, or mapping a gene to all pathways/GO terms it belongs to.
• Choosing between ORA, GSEA, KEGG, Reactome, or WikiPathways for a given input, and applying correct multiple-testing correction.

Version Compatibility

• Python ≥3.10, scipy ≥1.11, statsmodels ≥0.14, numpy ≥1.26.
• R: clusterProfiler ≥4.10, fgsea ≥1.28 (Bioconductor ≥3.18), org.Hs.eg.db for human gene ID mapping.
• KEGG REST API (rest.kegg.jp) has no version — always log the query date since pathway curation changes over time.

Prerequisites

• pip install scipy statsmodels numpy (Python) or BiocManager::install(c("clusterProfiler","fgsea","org.Hs.eg.db")) (R).
• A gene list (for ORA) or a full ranked gene table with a numeric score, e.g. from bio-differential-expression-deseq2-basics or bio-core-gene-ontology.
• Gene symbols normalized to one namespace (HGNC symbol or Entrez ID) before joining against KEGG/GO gene sets.

GO Structure

Three orthogonal ontologies, each a Directed Acyclic Graph (DAG):

Evidence code quality (filter for high-confidence work — keep IDA/IMP/IPI/IGI/IEP, drop IEA):

ORA: Hypergeometric Test with GO Propagation

Goal: test whether a fixed gene list is over-represented in a GO term or pathway, correcting for the true-path rule (a gene annotated to a term is implicitly annotated to every ancestor term).

Approach: propagate annotations up the DAG, build a term→genes map, then run a one-tailed hypergeometric test per term and apply BH-FDR across all terms tested.

from scipy import stats


def propagate_annotations(gene_annotations, go_terms, get_ancestors):
    """Apply the true-path rule: add every ancestor GO term to each gene's set.

    gene_annotations: dict gene -> list of (go_id, evidence_code)
    go_terms: dict go_id -> {'name', 'domain', 'parents'}
    get_ancestors: callable(go_id, go_terms) -> set of ancestor go_ids
    Returns dict gene -> set of GO IDs (direct + propagated).
    """
    propagated = {}
    for gene, annots in gene_annotations.items():
        all_terms = set()
        for go_id, _evidence in annots:
            all_terms.add(go_id)
            all_terms |= get_ancestors(go_id, go_terms)
        propagated[gene] = all_terms
    return propagated


def go_enrichment(gene_list, term_to_genes, background_size=20000):
    """Hypergeometric (ORA) enrichment of gene_list across GO terms/pathways.

    term_to_genes: dict term_id -> set of annotated genes (already propagated)
    background_size (N): number of genes actually tested in the assay,
        NOT the whole genome -- wrong background inflates/deflates results.
    Returns list of dicts sorted by p-value, with BH FDR added.
    """
    gene_set = {g.upper() for g in gene_list}
    n = len(gene_set)
    N = background_size

    results = []
    for term_id, term_genes in term_to_genes.items():
        term_genes_upper = {g.upper() for g in term_genes}
        overlap = gene_set & term_genes_upper
        k = len(overlap)
        if k == 0:
            continue
        K = len(term_genes_upper)
        # P(X >= k): probability of seeing k or more hits by chance
        p_value = stats.hypergeom.sf(k - 1, N, K, n)
        expected = n * K / N
        results.append({
            'term_id': term_id, 'k': k, 'K': K,
            'fold_enrichment': k / expected if expected > 0 else float('inf'),
            'p_value': p_value, 'genes': sorted(overlap),
        })

    results.sort(key=lambda r: r['p_value'])
    _add_bh_fdr(results)
    return results


def _add_bh_fdr(results):
    """Add Benjamini-Hochberg FDR to a p-value-sorted (ascending) results list, in place."""
    m = len(results)
    for i, r in enumerate(results):
        r['fdr'] = min(r['p_value'] * m / (i + 1), 1.0)
    running_min = 1.0
    for r in reversed(results):  # enforce monotonicity
        running_min = min(running_min, r['fdr'])
        r['fdr'] = running_min

GSEA: Ranked Gene List, No Threshold

Goal: score whether a gene set is enriched at the top or bottom of a fully ranked gene list (e.g. all genes by log fold-change), avoiding an arbitrary significance cutoff.

Approach: walk the ranked list, increment a running-sum statistic when a gene is in the set (weighted by |rank metric|) and decrement when it isn't; the enrichment score (ES) is the maximum deviation from zero. The genes contributing up to the ES are the "leading edge" — the core driver signal.

def gsea_enrichment_score(ranked_genes, gene_set):
    """Compute the GSEA running-sum enrichment score for one gene set.

    ranked_genes: list of (gene_name, rank_metric), sorted descending by metric
        (e.g. log fold-change or signed -log10(p)).
    gene_set: set of gene symbols in the pathway/term of interest.
    Returns dict with 'ES' (signed enrichment score), 'running_sum', 'hits'
        (indices of gene_set members in ranked_genes -- the leading edge is
        the hits up to argmax|running_sum|).
    """
    gene_set_upper = {g.upper() for g in gene_set}
    n_total = len(ranked_genes)
    n_hits = sum(1 for g, _ in ranked_genes if g.upper() in gene_set_upper)
    n_miss = n_total - n_hits
    if n_hits == 0 or n_miss == 0:
        return {'ES': 0.0, 'running_sum': [], 'hits': []}

    norm = sum(abs(m) for g, m in ranked_genes if g.upper() in gene_set_upper)

    running_sum, hits, current = [], [], 0.0
    for i, (gene, metric) in enumerate(ranked_genes):
        if gene.upper() in gene_set_upper:
            current += abs(metric) / norm if norm > 0 else 0.0
            hits.append(i)
        else:
            current -= 1.0 / n_miss
        running_sum.append(current)

    max_pos, min_neg = max(running_sum), min(running_sum)
    es = max_pos if abs(max_pos) >= abs(min_neg) else min_neg
    return {'ES': es, 'running_sum': running_sum, 'hits': hits}

# Same task in R with fgsea (production-grade, uses permutation p-values)
library(fgsea)

# ranked_genes: named numeric vector, names = gene symbols, values = logFC or stat,
# sorted is not required -- fgsea sorts internally.
# pathways: named list, name = pathway/GO id, value = character vector of gene symbols.
res <- fgsea(pathways = pathways, stats = ranked_genes,
             minSize = 10, maxSize = 500, eps = 0)
res <- res[order(res$padj), ]
# res$leadingEdge is the leading-edge gene list per pathway (list-column)

KEGG REST API

Goal: query live KEGG pathway data (no local database needed) for pathway lists, gene membership, or gene→pathway lookups.

Approach: hit rest.kegg.jp endpoints directly; each returns tab-separated text, so parse line-by-line and always handle network failure explicitly (the API is a shared public resource and can be rate-limited or briefly down).

import urllib.request
import urllib.error


def kegg_api_request(operation, *args, timeout=15):
    """GET a KEGG REST API endpoint: rest.kegg.jp/<operation>/<args...>.

    operation: one of list, get, find, link, conv.
    Returns the raw response text, or None on network failure.
    """
    url = "https://rest.kegg.jp/" + "/".join([operation, *args])
    try:
        with urllib.request.urlopen(url, timeout=timeout) as response:
            return response.read().decode("utf-8")
    except (urllib.error.URLError, urllib.error.HTTPError) as exc:
        print(f"KEGG API error: {exc}")
        return None


def kegg_pathway_genes(pathway_id):
    """Return KEGG gene IDs (e.g. 'hsa:7157') linked to a pathway, e.g. 'hsa04210'."""
    organism = pathway_id[:3]  # 'hsa' for human
    text = kegg_api_request("link", organism, f"pathway:{pathway_id}")
    if text is None:
        return []
    return [line.split("\t")[1].strip() for line in text.strip().split("\n") if "\t" in line]


# pathways = kegg_api_request("list", "pathway", "hsa")           # all human pathways
# genes_in_pathway = kegg_pathway_genes("hsa04210")                 # apoptosis
# pathway_for_gene = kegg_api_request("link", "pathway", "hsa:7157")  # TP53

Multiple Testing Correction

Bonferroni is overly conservative for GO/pathway terms because they are correlated (a gene's annotation propagates to many ancestors) — always prefer Benjamini-Hochberg FDR.

from statsmodels.stats.multitest import multipletests

# p_values: list/array of raw p-values, one per term tested
_, adj_pvals, _, _ = multipletests(p_values, method="fdr_bh")

Tool Decision Table

Pitfalls

• KEGG is dynamic: pathways are manually curated and updated; a gene list classified in 2020 may differ in 2024 — always log the KEGG release date/query date.
• ORA vs GSEA: ORA needs a binary gene list from an arbitrary threshold. GSEA uses the full ranking and is more powerful — no threshold, and it surfaces a "leading edge" subset as the core signal.
• True path rule: propagate GO annotations to all ancestor terms before testing, or broad terms ("cell death") will be undercounted relative to specific ones ("apoptosis").
• Background gene set matters: use the genes actually expressed/tested in your assay, not the whole genome, or enrichment is inflated/deflated.
• IEA evidence codes are low-quality: electronically inferred (IEA) annotations are unreviewed; filter to experimental (IDA/IMP/IPI/IGI/IEP) or curated computational (ISS/IBA) codes for high-confidence calls.
• Gene ID mismatches silently shrink results: GO/KEGG use Entrez or KEGG-prefixed IDs (hsa:7157); map symbols/Ensembl IDs first and check the mapping success rate.

See Also

• bio-core-gene-ontology — GO DAG structure, evidence codes, and ancestor traversal in more depth.
• bio-differential-expression-deseq2-basics — producing the ranked/DE gene list that feeds ORA or GSEA.
• bio-pathway-analysis-reactome-pathways — Reactome-specific reaction-level pathway queries.
• bio-workflows-expression-to-pathways — end-to-end pipeline from counts to pathway enrichment.