← All skills

pathway-enrichment

analysis

Pathway and gene-set enrichment on a gene list or ranked table — over-representation (Enrichr, g:Profiler), preranked GSEA, and ssGSEA/GSVA against GO, KEGG, Reactome, and MSigDB, with background choice and FDR handled correctly.

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

Use this skill for full, defensible enrichment workflows. For a single quick lookup, hitting the Enrichr or Reactome API directly is lighter weight.

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 pathway score per sample/cell ssGSEA / GSVA gp.ssgsea(...), gp.gsva(...)
You need a custom background or 500+ 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 gprofiler-official
# 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).

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")   # 200+ 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

# Enrichr libraries expect HGNC gene SYMBOLS (human: UPPERCASE). Map IDs first if needed.
genes = [g.strip() for g in open("deg_symbols.txt") if g.strip()]

enr = gp.enrichr(
    gene_list=genes,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2023",
               "KEGG_2021_Human", "Reactome_2022"],
    organism="human",
    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", "Overlap", "Adjusted P-value", "Combined Score", "Genes"]].head(20))

Preranked GSEA from DESeq2 results

import gseapy as gp
import pandas as pd

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)
rnk.index = rnk.index.str.upper()
rnk = rnk[~rnk.index.duplicated(keep="first")]

pre = gp.prerank(
    rnk=rnk,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2023"],
    min_size=15, max_size=500,        # drop tiny/huge sets (noisy or generic)
    permutation_num=1000, seed=123,   # seed = reproducible p-values
    threads=4, 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))

If you have no stat column, build the rank from sign(log2FoldChange) * -log10(pvalue).

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

Enrichr/MSigDB libraries are keyed by gene symbols (human UPPERCASE, mouse Title-case). If you have Ensembl/Entrez IDs, convert first. 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 inflates significance. Enrichr uses a fixed background; when background matters, use g:Profiler with domain_scope='custom' + your background, or gp.enrich() with an explicit background. Rationale in references/interpretation.md.

Step 5 — Run the analysis

Use the Quick Start patterns. For GSEA always set a seed and report permutation_num.

Step 6 — Filter on adjusted p-values

Use Adjusted P-value (ORA, Benjamini–Hochberg) or FDR q-val (GSEA), not raw p-values. Typical cutoff 0.05; also check the overlap/gene count so a "hit" isn't 1 gene out of a 2000-gene set.

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 many near-duplicate terms. Collapse with an enrichment map (term–term similarity), leading-edge overlap, or parent terms, and report representative terms. Interpretation framework and a publication-table format are in references/interpretation.md.

Boilerplate That Is Easy To Get Wrong

Symbol cleanup and rank construction cause most of the silent failures. Do both explicitly rather than trusting the input file.

import pandas as pd

def clean_symbols(genes, organism="human"):
    """Dedup, drop NA/blank, and match the casing the libraries use."""
    s = pd.Series(list(genes), dtype="string").dropna().str.strip()
    s = s[s.ne("")]
    s = s.str.upper() if organism == "human" else s.str.capitalize()
    return s.drop_duplicates().tolist()

def rank_from_deseq2(path):
    """Preranked GSEA input from a DESeq2 table, `stat` preferred over LFC."""
    res = pd.read_csv(path, index_col=0)
    if "stat" in res:
        rnk = res["stat"].dropna()
    else:
        import numpy as np
        d = res.dropna(subset=["pvalue", "log2FoldChange"])
        rnk = np.sign(d.log2FoldChange) * -np.log10(d.pvalue.clip(lower=1e-300))
    rnk.index = rnk.index.str.upper()
    return rnk[~rnk.index.duplicated(keep="first")].sort_values(ascending=False)

FDR is computed within a library, so filter per library rather than across the concatenated table:

sig = (res.groupby("Gene_set", group_keys=False)
          .apply(lambda g: g[g["Adjusted P-value"] < 0.05])
          .sort_values("Adjusted P-value"))

Common Pitfalls

These cause most wrong or irreproducible results:

  1. Gene-ID / organism mismatch — symbols vs Ensembl, human vs mouse casing. Map IDs and set organism correctly, or matches silently drop to ~zero.
  2. Wrong background (ORA) — using the whole genome instead of the tested/expressed gene set inflates p-values. Set a custom background when it matters.
  3. Thresholding before GSEA — GSEA needs the full ranked list; only ORA uses a cut list.
  4. Ranking GSEA by log2FoldChange alone — unstable for low-count genes; prefer stat or sign(LFC) * -log10(p).
  5. Multiple-testing across libraries — FDR is computed within a library; running many libraries multiplies tests. Report per-library FDR and stay conservative.
  6. Redundant GO terms — don't report 40 variants of the same term; collapse and show representatives.
  7. Significance ≠ relevance — check the overlap count and gene-set size; tiny sets reach significance trivially.
  8. List too short/long for ORA — <10 genes is underpowered; >2000 loses specificity (consider GSEA instead).
  9. No reproducibility metadata — Enrichr/GO libraries are versioned and drift over time. Record library names+date and set a GSEA seed.

Where The Genes Come From

Enrichment is always a second step. The usual upstream sources and what to carry forward from each:

  • Bulk RNA-seq differential expression — the results table. Carry the stat column for preranked GSEA, and the significant subset for ORA.
  • Single-cell marker detection — per-cluster marker genes with their scores.
  • CRISPR or drug screens — the hit list, plus the set of genes actually targeted by the library, which is the correct ORA background.
  • Proteomics — the identified protein set mapped to gene symbols first.

Whatever the source, record which genes were testable, not just which were hit — Step 4 needs it.

Reference Files

Read the relevant file when you need depth:

  • references/gseapy.md — full gseapy API: enrichr, offline enrich, prerank, gsea, ssgsea, gsva, Msigdb, Biomart, get_library_name/read_gmt, every plot, result-column meanings, GMT/offline usage, and troubleshooting (rate limits, empty results).
  • references/databases-and-gene-sets.md — GO, KEGG, Reactome, WikiPathways, MSigDB collections, Enrichr library naming, g:Profiler sources, organism handling, gene-ID conversion, library selection by question, and pointers to Reactome/STRING APIs and decoupler activity inference.
  • references/interpretation.md — ORA vs GSEA statistics, background-universe choice, multiple-testing methods (BH vs g:SCS vs Bonferroni), leading-edge genes, redundancy reduction, effect vs significance, a publication-table template, and reproducibility checklist.

Resources