/sc-enrichment

name: sc-enrichment
description: Load when running bulk-style pathway enrichment (ORA / GSEA / GSEA-R / GSVA-R) on a per-group
  ranked DE / marker list against a gene-set library. Skip when computing per-cell pathway scores in-place
  (use sc-pathway-scoring); de-novo gene-program discovery (use sc-gene-programs).
trigger: sc enrichment, single-cell enrichment, GO enrichment, KEGG enrichment, GSEA, ORA, pathway enrichment
tags:
- singlecell
- scrna
- enrichment
- gsea
- ora
- gsva
- decoupler
- pathway-enrichment

sc-enrichment

When to use

Test named gene sets against per-group marker/DE rankings with ORA or GSEA. GSVA scores mean expression per group in R. For scores per cell use `sc-pathway-scoring`; for programmes discovered without gene sets use `sc-gene-programs`.

Steps

Start a step with `python skills/_sdk/notebook/run.py new enrichment`. Load `enrichment = load_skill("sc-enrichment")`, then call `rank_groups`, `load_gene_sets` (GMT, JSON or Enrichr library), and `ora` or `gsea`. Pass the complete tested gene universe as `background` for ORA; its API default is all genes in the supplied ranking. Write returned tables and Figures with `write_output`. See `examples/example_step.py`.

The default engine is Python. Explicit `engine="auto"` prefers R when its packages are available; `engine="r"` requires clusterProfiler. The legacy `gsea_r` method uses GO_BP/KEGG/Reactome annotation databases, not custom GMT. `gsva` accepts a gene-set mapping, or None for the GO_BP/KEGG R database route.

Inputs & Outputs

Rankings need gene/names and a score or effect, with an optional group column. `rank_groups` reads log-normalised `X`, never `raw`. For the processed PBMC demo, use `.raw.to_adata()` in new steps; the CLI retains its old scaled-X demo ranking for compatibility. The CLI accepts H5AD or an upstream output directory containing `processed.h5ad` and marker/DE tables.

ORA/GSEA CLI writes `processed.h5ad`, `report.md`, `result.json`, `reproducibility/`, and five tables: `enrichment_results.csv`, `enrichment_significant.csv`, `group_summary.csv`, `ranking_input.csv`, `top_terms.csv`. Python GSEA adds `tables/gsea_running_scores.csv` when running curves can be built. `figure_data/` mirrors plot tables. Figures depend on available terms. GSVA uses `tables/gsva_r_scores.csv` and `figures/gsva_r_heatmap.png` instead. See `references/output_contract.md` for R-specific files.

Gotchas

• `result.json.summary.resolved_engine` records the engine actually used.

The default is now `python`; ask for `--engine auto` to retain auto-selection.

• `tables/enrichment_results.csv` records Python GSEA's local fallback in

`engine`, with the reason in `result.json.summary.warnings`. This is not interchangeable with gseapy's permutation implementation.

• `tables/enrichment_significant.csv` retains the legacy fixed FDR 0.05

filter. `--fdr-threshold` controls `group_summary.csv` and the report summary.

• `_api.py:37`: `rank_groups` does not infer a suitable matrix. Using scaled `X` can yield

invalid log-fold changes; use a log-normalised snapshot instead.

• `_api.py:103`: empty gene sets fail validation. Failed R annotation mapping raises an

error; the R scripts no longer manufacture pathways to fill an empty result.

• `_api.py:205`: GSVA's R database route supports GO_BP and KEGG, not Enrichr's Hallmark

alias. For Hallmark, explicitly load gene sets and pass the mapping to `gsva`.

Key CLI

python skills/singlecell/scrna/sc-enrichment/sc_enrichment.py --demo --output /tmp/sc_enrich_demo
python skills/singlecell/scrna/sc-enrichment/sc_enrichment.py --input clustered.h5ad --output results/ora --gene-set-db hallmark --groupby cell_type
python skills/singlecell/scrna/sc-enrichment/sc_enrichment.py --input clustered.h5ad --output results/gsea --method gsea --gene-sets pathways.gmt --gsea-seed 123
python skills/singlecell/scrna/sc-enrichment/sc_enrichment.py --input clustered.h5ad --output results/gsva --method gsva_r --groupby cell_type --gene-sets pathways.gmt

API

<!-- api:begin generated from _api.py; regenerate with run.py api <skill dir> --write -->

`rank_groups(adata, *, groupby: str, method: str='wilcoxon') -> pd.DataFrame`

Rank each group's genes against the rest using log-normalised ``X``.

``method`` is wilcoxon (default), t-test or logreg. This uses ``X``, not ``raw``; convert a scaled object with ``adata.raw.to_adata()`` first when raw holds log-normalised values. Scanpy ranking statistics remain in ``uns``. Return a table with group, gene and available scores/p-values.

`load_gene_sets(source, *, species: str='human', universe=None) -> dict[str, list[str]]`

Load GMT/JSON or fetch an Enrichr library; optionally match genes to a universe.

``source`` accepts a local path or hallmark/kegg/reactome/go_bp aliases and full Enrichr library names. Remote sources need network access and gseapy; local files do not. ``species`` is human (default) or mouse. Return a term-to-gene-list mapping; no output files are written.

`demo_gene_sets(*, species: str='human') -> dict[str, list[str]]`

Return the six named PBMC demo signatures, human by default or mouse symbols.

`marker_gene_sets(markers: pd.DataFrame, *, groups=None, top_n: int | None=100, universe=None) -> dict[str, list[str]]`

Convert marker groups to gene sets, sorted by adjusted p-value or score.

``markers`` needs group and names (or gene). ``groups=None`` selects all; ``top_n=None`` keeps all genes. ``universe`` optionally restricts and canonicalises gene symbols. Raise ValueError when no sets remain.

`ora(ranking: pd.DataFrame, gene_sets, *, background=None, engine: str='python', padj_cutoff: float=0.05, log2fc_cutoff: float=0.25, max_genes: int=200, source: str='custom', library_mode: str='local', n_top: int=18) -> pd.DataFrame`

Test over-representation with a local hypergeometric test or clusterProfiler.

``ranking`` accepts marker/DE columns and optional group (default