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)

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)

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.

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:

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

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:

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",
             # "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