# Bio Core Pathways

> 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.

- Skill: `pavel-kravchenko/bio-core-pathways` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-pathways`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-pathways/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Integrations & APIs
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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):

| 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.

```python
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.

```python
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}
```

```r
# 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).

```python
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.

```python
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.