bio-core-pathways

name	bio-core-pathways
description	Run GO/KEGG pathway enrichment via ORA (hypergeometric test) or GSEA, with BH-FDR and KEGG REST API queries. Use for "GO enrichment", "GSEA", "ORA vs GSEA", or "pathway enrichment" requests.
tool_type	python
primary_tool	scipy

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

Ontology	Abbreviation	Meaning
Molecular Function	MF	Biochemical activity of the gene product
Biological Process	BP	Pathway or larger biological process
Cellular Component	CC	Where in the cell the product is active

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

Tier	Codes	Quality
Experimental	EXP, IDA, IPI, IMP, IGI, IEP	High
Computational (curated)	ISS, ISO, ISA, ISM, IBA	Medium
Traceable/Non-traceable	TAS, NAS	Low
Electronic (IEA)	IEA	Lowest — filter out

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, : k, : K,
            : k / expected  expected >   (),
            : p_value, : (overlap),
        })

    results.sort(key= r: r[])
    _add_bh_fdr(results)
     results


 ():
    
    m = (results)
     i, r  (results):
        r[] = (r[] * m / (i + ), )
    running_min = 
     r  (results):  
        running_min = (running_min, r[])
        r[] = 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 = (running_sum), (running_sum)
    es = max_pos  (max_pos) >= (min_neg)  min_neg
     {: es, : running_sum, : 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

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

Scenario	Tool
Gene list, binary significant/not	ORA (hypergeometric / Fisher exact)
Ranked gene list (all tested genes)	GSEA (fgsea, clusterProfiler, GSEApy)
Human pathways with reactions	Reactome (ReactomePA)
Community-curated, open-access	WikiPathways
Organism-specific metabolism/signaling	KEGG
All GO terms + FDR in one call	goatools, clusterProfiler, g:Profiler

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.