pathway-enrichment – Data & Science skill for AI agents

Pathway Enrichment

Overview

Enrichment analysis answers "what biology is over-represented in my genes?" It is the standard last step after differential expression, a screen, or clustering. There are two core methods, and choosing correctly is the single most important decision:

• ORA (over-representation analysis) — take a *thresholded* gene list (e.g., padj < 0.05) and test which gene sets it overlaps more than chance, using Fisher's exact / hypergeometric tests. Tools: Enrichr, g:Profiler.
• GSEA (gene set enrichment analysis) — take the *whole ranked list* of genes (no threshold) and test whether each gene set is concentrated toward the top or bottom. Preranked GSEA uses a per-gene score (e.g., the DESeq2 stat). Better when effects are broad and subtle.

This skill orchestrates these analyses, the gene-set databases behind them, and the interpretation pitfalls that make results wrong or unpublishable.

When to Use This Skill

Use this skill when the user wants to:

• Find enriched GO terms / KEGG / Reactome / WikiPathways / MSigDB Hallmark sets in a gene list.
• Run GSEA / preranked GSEA on DESeq2, edgeR, limma, or Scanpy rank_genes_groups output.
• Score gene-set expression per sample/cell (ssGSEA, GSVA).
• Interpret, deduplicate, and visualize enrichment results, or build a publication table/figure.
• Decide between ORA and GSEA, pick gene-set libraries, choose a background, or fix gene-ID problems.

For quick one-off Enrichr lookups the gget skill (gget enrichr) is lighter weight; for raw pathway/interaction APIs (Reactome, KEGG, STRING) see the database-lookup skill. Use this skill for full, defensible enrichment workflows.

Choosing the Right Method

| Situation | Method | Tool / entry point |

|-----------|--------|--------------------|

| You have a discrete hit list (DE genes, screen hits, cluster markers) | ORA | gp.enrichr(...) or g:Profiler |

| You have a full ranked list (every tested gene + a score) | Preranked GSEA | gp.prerank(...) |

| You have an expression matrix + class labels | GSEA | gp.gsea(...) |

| You want a gene-set score per sample/cell | ssGSEA / GSVA | gp.ssgsea(...), gp.gsva(...) |

| You need a custom background or additional organisms | ORA with custom domain | g:Profiler (domain_scope='custom') |

| You want TF / signaling *activity* (PROGENy, DoRothEA) | activity inference | see references/databases-and-gene-sets.md (decoupler) |

When in doubt: a thresholded list → ORA; a ranked table with scores → GSEA. Never threshold a list and then feed it to GSEA — that discards the ranking GSEA depends on.

Setup

uv pip install gseapy==1.3.1 gprofiler-official==1.0.0 mygene==3.2.2 lxml
# gseapy pulls pandas, numpy, scipy, matplotlib. Network access is needed for
# Enrichr, g:Profiler, and MSigDB downloads. For fully offline ORA, use a local
# GMT file with gp.enrich() (see references/gseapy.md).

The examples target GSEApy 1.3.1 (released stable); the helper runs classic

permutation GSEA. Current MSigDB is 2026.1.Hs/2026.1.Mm. Local synthetic workflows

and small public mapping/catalog queries were executed; Enrichr submission

routes were verified from released source with mocked transport. BioMart returned

service-unavailable HTML during review; validate its output schema. See

[verified API contracts](references/verified-api.md) for evidence and boundaries.

Verify and list available gene-set libraries (names change over time — never hardcode blindly):

import gseapy as gp
names = gp.get_library_name(organism="human")   # discover current Enrichr libraries
print([n for n in names if "Reactome" in n or "KEGG" in n or "Hallmark" in n])

Quick Start

ORA on a hit list (gseapy + Enrichr)

import gseapy as gp

# Match the library organism and identifier namespace; map aliases explicitly.
# Keep original spelling: capitalization is not gene-ID or orthology mapping.
genes = [g.strip() for g in open("deg_symbols.txt") if g.strip()]
tested = [g.strip() for g in open("tested_symbols.txt") if g.strip()]
assert set(genes) <= set(tested)

enr = gp.enrichr(
    gene_list=genes,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2026",
               "KEGG_2026", "Reactome_Pathways_2024"],
    organism="human", background=tested,
    outdir=None,            # in-memory; set a path to also write tables/plots
)
res = enr.results
sig = res[res["Adjusted P-value"] < 0.05].sort_values("Adjusted P-value")
print(sig[["Gene_set", "Term", "Adjusted P-value", "Combined Score", "Genes"]].head(20))
# Speedrichr results omit Overlap; retain its actual response schema.

Preranked GSEA from DESeq2 results

import gseapy as gp
import pandas as pd
import numpy as np

res = pd.read_csv("deseq2_results.csv", index_col=0)   # index = gene symbols
# Rank by the test statistic (sign = direction, magnitude = evidence). This is
# more stable than ranking by log2FoldChange, which is noisy for low-count genes.
rnk = res["stat"].dropna().sort_values(ascending=False, kind="stable")
rnk.index = rnk.index.str.strip()
assert rnk.index.is_unique, "Resolve duplicate mappings before ranking"
assert np.isfinite(rnk.to_numpy()).all()

pre = gp.prerank(
    rnk=rnk,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2026"],
    organism="human", method="permutation",
    min_size=15, max_size=500,        # size AFTER intersection with ranked genes
    permutation_num=1000, seed=123,   # seed = reproducible p-values
    threads=4, ascending=None, outdir=None,
)
out = pre.res2d.sort_values("FDR q-val")
print(out[["Term", "ES", "NES", "NOM p-val", "FDR q-val", "Lead_genes"]].head(20))

Use a signed Wald statistic, not the unsigned DESeq2 likelihood-ratio statistic.

If unavailable, sign(log2FoldChange) * -log10(raw pvalue) is a fallback: validate

p-values in [0, 1], bound numerical zeros (the helper uses 1e-300), and report ties

and exclusions. Do not use adjusted p-values or select only significant genes.

Core Workflow

For a defensible analysis, work through these steps. The middle steps (ID type, background) are where results most often silently go wrong.

Step 1 — Pin down inputs and pick the method

Confirm: which genes, what organism, is there a per-gene score (→ GSEA) or just a list (→ ORA), and what comparison they represent (direction matters for interpretation).

Step 2 — Get gene IDs into the right namespace

Match the identifiers actually stored in the selected library: MSigDB offers symbol and Entrez GMTs. Human/mouse capitalization is a convention, not a conversion. Preserve original IDs, resolve one-to-many mappings deliberately, and map the ORA query and background identically. See references/databases-and-gene-sets.md for gp.Biomart, g:Profiler g:Convert, and mygene. A silent ID mismatch is the #1 cause of "nothing is significant".

Step 3 — Choose gene-set libraries to match the question

Hallmark (broad themes) → GO:BP (mechanism) → KEGG/Reactome/WikiPathways (curated pathways) → C7 (immune), etc. Don't run 50 libraries; pick 2–4 that fit the biology. Catalog and selection guidance: references/databases-and-gene-sets.md.

Step 4 — Set the background universe (ORA only)

The background must be the genes that *could* have been detected in your assay (e.g., all expressed/tested genes), not the whole genome. The wrong background biases significance. Every query gene must belong to that universe. GSEApy 1.3.1 uses Speedrichr for an explicit online background; local gp.enrich() + a pinned GMT gives a directly inspectable universe. g:Profiler also accepts domain_scope='custom' + background. Do not assume a background gene count or service default represents the assay. Rationale in references/interpretation.md.

Step 5 — Run the analysis

Use the Quick Start patterns or the bundled scripts/run_enrichment.py. For GSEA always set a seed and report permutation_num.

Step 6 — Filter on adjusted p-values

Use the correction returned by the selected method: Enrichr ORA reports BH-adjusted

p-values, g:Profiler defaults to g:SCS, and classic GSEA estimates FDR q-val from its

permutation distributions. GSEApy 1.3.1 method="multilevel" instead returns

BH-adjusted p-values across all tested terms and a log2err diagnostic. These are not interchangeable BH outputs. Report the

method and permutation type with the cutoff; GSEA's exploratory 0.25 convention

is for phenotype permutations, while its documentation recommends 0.05 for

gene-set permutations such as preranked analyses. Also inspect overlap and

gene-set size. See the [GSEA FAQ](https://docs.gsea-msigdb.org/GSEA/GSEA_FAQ/).

Step 7 — Visualize

Dotplots, bar plots, enrichment maps, and GSEA running-score plots are built into gseapy (gp.dotplot, gp.barplot, gp.enrichment_map, gp.gseaplot). See references/gseapy.md.

Step 8 — Reduce redundancy and interpret

GO especially returns m

… (full content is in the download package)