/sc-pathway-scoring

name: sc-pathway-scoring
description: Load when computing per-cell pathway / gene-set scores on a normalised scRNA AnnData via
  AUCell (R or Python) or Scanpy score_genes. Skip when running condition-vs-control bulk-style enrichment
  on top of a DE table (use sc-enrichment); de-novo gene-program discovery (use sc-gene-programs).
trigger: pathway score, pathway scoring, gene set score, module score, pathway activity, signature score
tags:
- singlecell
- scrna
- pathway-scoring
- aucell
- scanpy-score-genes
- gene-sets
- decoupler

sc-pathway-scoring

When to use

The user has a normalised scRNA AnnData and a gene-set library (GMT file or one of the built-in DB aliases: `hallmark`, `kegg`, `reactome`, `go_bp`, ...) and wants per-cell scores quantifying how active each gene set is. Three methods:

• `aucell_r` (default) — R-backed AUCell via `decoupler-py`-style

bridge. Best statistical foundation; requires R env.

• `aucell_py` — Python AUCell (`--aucell-py-auc-threshold`). Pure

Python.

• `score_genes_py` — Scanpy `tl.score_genes` per gene set. Lightest

and fastest.

Output: `tables/enrichment_scores.csv` (cells × gene_sets), plus group-mean / group-high-fraction tables when `--groupby` is provided.

For *bulk-style* condition-vs-control GSEA / ORA on a DE table use `sc-enrichment`. For de-novo gene-program discovery use `sc-gene-programs`.

Use in an analysis step

Use `load_skill` from the notebook SDK; write returned objects with `write_output`. This runnable example is also in `examples/example_step.py`. The CLI remains available for standalone reports and galleries.

# Score PBMC lineage gene sets with the Python AUCell implementation.
# Reads pbmc3k_processed and uses its log-normalized raw snapshot.
# Calls sc-pathway-scoring: score_gene_sets, group_scores, score_distribution_figure.

from skills._sdk.notebook import load_demo, load_skill, write_output

pathways = load_skill('sc-pathway-scoring')
adata = load_demo('pbmc3k_processed').raw.to_adata()
gene_sets = {
    'B_cell': ['MS4A1', 'CD79A', 'CD79B', 'CD74', 'HLA-DRA'],
    'T_cell': ['CD3D', 'CD3E', 'CD3G', 'TRAC', 'IL7R'],
    'Myeloid': ['LYZ', 'S100A8', 'S100A9', 'FCN1', 'CTSS'],
}

scores = pathways.score_gene_sets(adata, gene_sets, method='aucell_py')
write_output(scores, 'tables/pathway_scores.csv')
write_output(pathways.group_scores(adata, scores, groupby='louvain'), 'tables/group_scores.csv')
write_output(pathways.score_distribution_figure(scores), 'figures/score_distribution.png')

assert scores.index.equals(adata.obs_names)
assert set(scores.columns) == set(gene_sets)
assert scores.max().max() > 0
assert scores.min().min() >= 0

Inputs & Outputs

Input is an AnnData and gene sets keyed by name. The CLI reads H5AD plus a GMT file or an Enrichr library. Python scoring returns a cell-by-set DataFrame; `attach_scores` returns a copy with `obs["enrich__..."]` columns.

The CLI writes `processed.h5ad`, `report.md`, `result.json`, and `tables/enrichment_scores.csv`, `gene_set_overlap.csv`, `top_pathways.csv`. Grouped runs also write `group_mean_scores.csv` and `group_high_fraction.csv`. Plot source tables live under `figure_data/`, not `tables/`. R exchange matrices and AUCell CSVs are temporary.

Flow

1. Load AnnData (`--input`) or build a demo. 2. Load gene sets: parse `--gene-sets` GMT, OR fetch via `--gene-set-db <alias>` and write a resolved GMT. 3. Validate: at least one gene-set member overlaps the input features (`feature_label_source` chosen from `var_names` / `var["gene_symbol"]` / etc.). 4. Run preflight; resolve `--groupby` (auto-pick from `leiden` / `louvain` / `cell_type` if unset). 5. Dispatch to method:

• `aucell_r`: shell out to bundled R script via `RScriptRunner`.
• `aucell_py`: AUCell-Python with `--aucell-py-auc-threshold`.
• `score_genes_py`: Scanpy `tl.score_genes` per gene set.

6. Compute group-aware aggregates if `--groupby` is set. 7. Save tables, figures, `processed.h5ad`, `report.md`, `result.json`.

Gotchas

• `run_info(scores)["skipped_gene_sets"]` lists sets with no matching features. Inspect `tables/gene_set_overlap.csv` before interpreting scores; all-unmatched Python input raises.
• `score_gene_sets(..., method="score_genes_py")` requires normalized X. AUCell ranks expression and can also use counts.
• `load_gene_sets` downloads named libraries through gseapy. Local GMT/JSON avoids network access. This skill's mouse `kegg` alias is `KEGG_2021_Mouse`, intentionally separate from sc-enrichment's aliases.
• `aucell_r` needs AUCell and GSEABase. Its temporary result must contain a `Cell` column; otherwise the API raises instead of misaligning rows.
• API scoring uses seed 42 by default. The CLI preserves its historical `score_genes_py` seed 0; its AUCell seed defaults to 42. Set `random_state=0` for API/CLI score_genes comparisons.
• `processed.h5ad` adds an `enrich__` obs column per gene set. `attach_scores` copies its input; `score_gene_sets` leaves it unchanged.
• CLI `--demo` slices the first 60 feature names into four arbitrary sets. They exercise the pipeline, not biological pathways. The step below uses named PBMC lineage genes.

Key CLI

# Demo (built-in gene sets)
python skills/singlecell/scrna/sc-pathway-scoring/sc_pathway_scoring.py --demo --output /tmp/sc_pw_demo

# AUCell-R with MSigDB Hallmark, grouped by cell type
python skills/singlecell/scrna/sc-pathway-scoring/sc_pathway_scoring.py \
  --input clustered.h5ad --output results/ \
  --gene-set-db hallmark --groupby cell_type

# AUCell-Python (no R needed) with custom GMT
python skills/singlecell/scrna/sc-pathway-scoring/sc_pathway_scoring.py \
  --input clustered.h5ad --output results/ \
  --method aucell_py --gene-sets pathways.gmt --groupby leiden

# Scanpy score_genes for fast prototyping
python skills/singlecell/scrna/sc-pathway-scoring/sc_pathway_scoring.py \
  --input clustered.h5ad --output results/ \
  --method score_genes_py --gene-sets pathways.gmt

##