# Bio Core Gene Ontology

> Decode GO evidence codes (IDA/IEA/TAS), propagate GO annotations via the true path rule, and run ORA (hypergeometric test + BH-FDR) in Python/R. Use when doing GO enrichment, evidence-code checks, or DAG annotation propagation.

- Skill: `pavel-kravchenko/bio-core-gene-ontology` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-gene-ontology`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-gene-ontology/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-core-gene-ontology

---


# Gene Ontology: Structure, Evidence Codes & Enrichment

## When to Use

- Explaining what a GO evidence code (IDA, IEA, TAS, ISS, ...) means and how reliable it is.
- Propagating a gene's direct GO annotations up the ontology DAG (true path rule) before counting.
- Running or debugging an over-representation analysis (ORA): hypergeometric test + BH-FDR on a DE gene list.
- Auditing an existing GO enrichment result for a wrong background set, missing FDR correction, or unpropagated annotations.
- Choosing between ORA (binary gene list) and GSEA (ranked list) for a given input.

## Version Compatibility

- Python: scipy ≥1.11, statsmodels ≥0.14, goatools ≥1.4, gseapy ≥1.1, Python ≥3.10
- R: clusterProfiler ≥4.10 (Bioconductor ≥3.18), org.Hs.eg.db ≥3.18, R ≥4.3

## Prerequisites

- `pip install scipy statsmodels goatools gseapy pandas` (goatools/gseapy for real OBO/GAF-backed enrichment; scipy is enough for the mechanics below)
- R: `BiocManager::install(c("clusterProfiler", "org.Hs.eg.db"))`
- Concepts: DAG traversal, hypergeometric distribution, multiple-testing correction
- Related skill for a full production pipeline against the real `go-basic.obo` + GAF: `bio-pathway-analysis-go-enrichment`

## Key Concepts

- **True Path Rule**: annotation to a specific term implies annotation to all its ancestors. Enrichment must propagate up the DAG or it undercounts genes in broad parent terms.
- **Background set**: use all genes expressed/tested in your experiment, not the whole genome — the wrong background inflates or deflates every p-value.
- **Multiple testing**: testing ~20,000 GO terms at α=0.05 gives ~1,000 false positives by chance alone. Always report BH-adjusted p-values (FDR), not raw p-values.
- **ORA vs GSEA**: ORA needs a binary gene list (e.g., DE genes above a threshold) and is simpler but throws away magnitude. GSEA uses the full ranked list (fold change, p-value) and is more powerful — prefer it when you have a quantitative ranking.

### Evidence codes (most to least reliable)

```text
EXP, IDA, IPI, IMP, IGI, IEP  — Experimental (gold standard, quality 4-5)
ISS, ISO, ISA, ISM, IBA        — Computational/sequence-based (quality 3)
TAS                            — Traceable Author Statement (quality 2)
NAS                            — Non-traceable Author Statement (quality 1)
IEA                            — Inferred from Electronic Annotation, automated, unreviewed (quality 1)
```

```python
EVIDENCE_QUALITY = {
    'EXP': 5, 'IDA': 5, 'IPI': 5, 'IMP': 5, 'IGI': 5, 'IEP': 4,
    'ISS': 3, 'ISO': 3, 'ISA': 3, 'ISM': 3, 'IBA': 3,
    'TAS': 2, 'NAS': 1, 'IEA': 1,
}

def filter_annotations_by_quality(gene_annotations, min_quality=3):
    """Keep only (GO_ID, evidence_code) pairs at or above min_quality.

    gene_annotations: dict of gene -> list of (go_id, evidence_code)
    Returns a dict with the same shape, genes with zero surviving
    annotations dropped entirely.
    """
    filtered = {}
    for gene, annots in gene_annotations.items():
        kept = [(go_id, ec) for go_id, ec in annots if EVIDENCE_QUALITY.get(ec, 0) >= min_quality]
        if kept:
            filtered[gene] = kept
    return filtered
```

## Propagating Annotations: the True Path Rule

**Goal:** turn each gene's direct GO annotations into the full set of terms it is implicitly annotated to (direct + all ancestors), so enrichment counts broad parent terms correctly.

**Approach:** walk each term's `parents` list with BFS, collecting every ancestor reached, then union that with the gene's direct terms.

```python
def get_ancestors(go_id, go_terms):
    """Return the set of all ancestor GO IDs of go_id (BFS over `parents`).

    go_terms: dict of GO_ID -> {'name': str, 'domain': str, 'parents': [GO_ID, ...]}
    Does not include go_id itself.
    """
    seen = set()
    queue = list(go_terms.get(go_id, {}).get('parents', []))
    while queue:
        parent = queue.pop()
        if parent in seen:
            continue
        seen.add(parent)
        queue.extend(go_terms.get(parent, {}).get('parents', []))
    return seen


def propagate_annotations(gene_annotations, go_terms):
    """Apply the true path rule: propagate each annotation to all ancestor terms.

    Returns dict gene -> set of GO IDs (direct + propagated).
    """
    propagated = {}
    for gene, annots in gene_annotations.items():
        all_terms = set()
        for go_id, _ec in annots:
            all_terms.add(go_id)
            all_terms |= get_ancestors(go_id, go_terms)
        propagated[gene] = all_terms
    return propagated
```

## GO Enrichment (ORA): Hypergeometric Test + BH-FDR

2x2 contingency table per term:

```text
                  | In GO term | Not in term | Total
In gene list      |     k      |    n - k    |   n
Not in gene list  |   K - k    |  N-K-n+k    | N - n
Total             |     K      |    N - K    |   N

N = background size; K = genes annotated to term in background
n = genes in your list; k = overlap
P(X >= k) = hypergeometric survival function -> stats.hypergeom.sf(k-1, N, K, n)
```

**Goal:** score every GO term for over-representation in a gene-of-interest list against a background, with FDR-corrected significance.

**Approach:** propagate annotations first (true path rule), build term→genes from the propagated set, run `hypergeom.sf` per term, then apply Benjamini-Hochberg to the sorted p-values.

```python
from scipy import stats


def go_enrichment(gene_list, gene_annotations, go_terms, background_size=20000):
    """Over-representation analysis with the hypergeometric test + BH-FDR.

    gene_list: genes of interest (e.g., DE genes)
    gene_annotations: gene -> list of (GO_ID, evidence_code)
    go_terms: GO_ID -> {'name', 'domain', 'parents'}
    background_size: total genes in the tested background (NOT whole genome)

    Returns a list of dicts sorted by p-value, each with 'fdr' added.
    """
    propagated = propagate_annotations(gene_annotations, go_terms)

    term_to_genes = {}
    for gene, terms in propagated.items():
        for t in terms:
            term_to_genes.setdefault(t, set()).add(gene)

    gene_set = set(g.upper() for g in gene_list)
    n, N = len(gene_set), background_size

    results = []
    for term_id, term_genes in term_to_genes.items():
        annotated_in_list = gene_set & {g.upper() for g in term_genes}
        k = len(annotated_in_list)
        if k == 0:
            continue
        # In production, K comes from the full annotation DB (goatools/GAF),
        # not scaled from this toy set.
        K = max(k, int(len(term_genes) / len(propagated) * N)) if propagated else k
        p_value = stats.hypergeom.sf(k - 1, N, K, n)
        expected = n * K / N
        info = go_terms.get(term_id, {})
        results.append({
            'go_id': term_id, 'name': info.get('name', 'unknown'), 'domain': info.get('domain', '?'),
            'k': k, 'K': K, 'p_value': p_value,
            'fold_enrichment': k / expected if expected > 0 else float('inf'),
            'genes': sorted(annotated_in_list),
        })

    results.sort(key=lambda r: r['p_value'])
    m = len(results)
    for i, r in enumerate(results):  # Benjamini-Hochberg
        r['fdr'] = min(r['p_value'] * m / (i + 1), 1.0)
    return results


def demo():
    """Runnable self-check: a small DAG + annotation set with a known signal."""
    go_terms = {
        'GO:0008150': {'name': 'biological_process', 'domain': 'BP', 'parents': []},
        'GO:0009987': {'name': 'cellular process', 'domain': 'BP', 'parents': ['GO:0008150']},
        'GO:0006915': {'name': 'apoptotic process', 'domain': 'BP', 'parents': ['GO:0009987']},
        'GO:0097193': {'name': 'intrinsic apoptotic signaling pathway', 'domain': 'BP', 'parents': ['GO:0006915']},
        'GO:0007049': {'name': 'cell cycle', 'domain': 'BP', 'parents': ['GO:0009987']},
        'GO:0000278': {'name': 'mitotic cell cycle', 'domain': 'BP', 'parents': ['GO:0007049']},
    }
    gene_annotations = {
        'TP53': [('GO:0006915', 'IDA'), ('GO:0000278', 'IMP')],
        'BAX':  [('GO:0097193', 'IDA')],
        'BCL2': [('GO:0097193', 'IMP')],
        'CDK2': [('GO:0000278', 'IDA')],
    }
    ancestors = get_ancestors('GO:0097193', go_terms)
    assert ancestors == {'GO:0006915', 'GO:0009987', 'GO:0008150'}, ancestors

    results = go_enrichment(['TP53', 'BAX', 'BCL2'], gene_annotations, go_terms, background_size=20000)
    top = results[0]
    assert top['go_id'] in ('GO:0006915', 'GO:0009987', 'GO:0008150')
    assert all(r['fdr'] >= r['p_value'] for r in results)
    print(f"OK: top term {top['go_id']} ({top['name']}), FDR={top['fdr']:.2e}")


if __name__ == '__main__':
    demo()
```

## Real-World Enrichment in R (clusterProfiler)

For production ORA against the live GO/KEGG databases (real background, real annotation DB), use `clusterProfiler::enrichGO` instead of hand-rolled hypergeometric code:

```r
library(clusterProfiler)
library(org.Hs.eg.db)

## de_genes: character vector of DE gene symbols; universe: all tested genes
de_genes <- c("TP53", "BAX", "CASP3", "CASP9", "BCL2", "CDK2")
universe  <- rownames(res)  # e.g., all genes tested in DESeq2 results

ego <- enrichGO(
  gene          = de_genes,
  universe      = universe,        # correct background, not the whole genome
  OrgDb         = org.Hs.eg.db,
  keyType       = "SYMBOL",
  ont           = "BP",            # "BP", "MF", "CC", or "ALL"
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  qvalueCutoff  = 0.2
)

ego <- clusterProfiler::simplify(ego, cutoff = 0.7)  # collapse redundant/near-duplicate terms
head(as.data.frame(ego))
```

## Pitfalls

- **Skipping propagation**: counting only direct annotations undercounts broad parent terms and misses the true path rule — always propagate before testing.
- **Wrong background**: using the whole genome/proteome as N when only a subset was tested (e.g., only expressed genes) biases every p-value.
- **No FDR correction**: reporting raw hypergeometric p-values across thousands of GO terms guarantees hundreds of false positives — always use BH (or `simplify()` in clusterProfiler to also cut redundant terms).
- **ID mismatch**: GO annotation databases key on Entrez/UniProt; your DE list may use symbols or Ensembl IDs — a silent mapping failure quietly shrinks the gene list before it ever gets tested.
- **Mixing evidence quality unfiltered**: IEA (automated, unreviewed) annotations dominate most GO databases numerically; blending them with experimental (IDA/IMP) evidence without tracking quality can manufacture enrichment from weak computational calls.

## See Also

- `bio-pathway-analysis-go-enrichment` — full production ORA pipeline against real `go-basic.obo` + GAF files
- `bio-pathway-analysis-gsea` — ranked-list gene set enrichment (GSEA) instead of binary ORA
- `bio-pathway-analysis-kegg-pathways` — pathway-level enrichment (KEGG) alongside GO
- `bio-expression-matrix-gene-id-mapping` — resolving symbol/Ensembl/Entrez ID mismatches before enrichment