Single-cell RNA-seq gives you thousands of gene-level measurements per cell, but biological interpretation usually happens at the level of pathways, programs, and processes. Pathway analysis helps compress noisy, high-dimensional expression data into interpretable scores, and decoupler is a good fit for this because it works directly with AnnData and integrates cleanly with Scanpy workflows.
Gene set enrichment asks whether a predefined set of genes, such as a pathway or hallmark signature, is coordinately shifted relative to the background. In decoupler‘s GSEA implementation, genes are ranked by a continuous statistic and the enrichment score is computed from a running-sum statistic over that ranked list.
Summary
In this tutorial you will:
- load a ready-to-use scRNA-seq
AnnDataobject - retrieve Hallmark pathway gene sets
- run per-cell GSEA with
decoupler - extract the score matrix from
adata.obsm - visualize pathway activity on UMAP
- rank pathways by cell type or cluster
Expected outcome: a Scanpy-compatible AnnData object containing per-cell pathway scores in adata.obsm["score_gsea"], plus plots that show which pathways are enriched in different cell populations. Hallmark provides 50 curated gene sets, which makes it a practical first resource for this kind of analysis.
Requirements
- Basic Python and pandas familiarity
- Basic knowledge of Scanpy and
AnnData - A recent Python environment
- Packages:
scanpy,decoupler - No GPU required for this example
- Jupyter Notebook or JupyterLab recommended
Step 1: Install and import the libraries
pip install scanpy decoupler
import scanpy as sc
import decoupler as dc
import pandas as pd
sc.set_figure_params(figsize=(4, 4), frameon=False)
decoupler exposes enrichment methods through the dc.mt.* API, including dc.mt.gsea, and is designed to work with AnnData objects used throughout Scanpy.
Step 2: Load a scRNA-seq dataset
To keep the tutorial focused on pathway analysis, we will start from a ready-to-use PBMC dataset that already contains log-normalized expression values and cell annotations.
adata = dc.ds.pbmc3k()
adata
adata.obs.head()
sc.pl.umap(adata, color="celltype")
This dataset is a PBMC 3k single-cell RNA-seq AnnData object with log-normalized transcript counts and metadata such as celltype and leiden, which makes it convenient for pathway scoring and Scanpy visualization.
If you already have your own Scanpy object, replace dc.ds.pbmc3k() with your adata. The key requirement is that adata.X should contain normalized, log-transformed expression values rather than raw counts.
Step 3: What gene set enrichment means in practice
For scRNA-seq, gene set enrichment turns a cell-by-gene matrix into a cell-by-pathway matrix. Instead of asking whether GeneA or GeneB is highly expressed, you ask whether a whole pathway, such as inflammatory signaling or hypoxia, is coordinately active in a cell population.
For this tutorial, we will use GSEA. It works well when your prior knowledge is an unweighted set of genes and you want a pathway score per cell.
Step 4: Download pathway gene sets
We will use the Hallmark collection.
hallmark = dc.op.hallmark(organism="human")
hallmark.head()
Check how many unique pathways are available:
hallmark["source"].nunique()
Prune the resource so that only pathways with enough overlapping genes remain:
hallmark = dc.pp.prune(adata.var_names, hallmark, tmin=5)
hallmark["source"].nunique()
dc.op.hallmark() returns Hallmark gene sets for a chosen organism, and the resource contains 50 curated sets. dc.pp.prune() removes sources that do not share at least tmin genes with your expression matrix, which prevents empty or weakly supported pathway scores.
Step 5: Run per-cell GSEA
Now compute pathway enrichment scores for each cell.
dc.mt.gsea(
data=adata,
net=hallmark,
tmin=5,
times=1,
seed=42,
verbose=True,
)
This creates per-cell GSEA scores and stores them in adata.obsm["score_gsea"].
If you also want permutation-based significance, rerun GSEA with more than one permutation:
dc.mt.gsea(
data=adata,
net=hallmark,
tmin=5,
times=100,
seed=42,
verbose=True,
)
[k for k in adata.obsm.keys() if "gsea" in k]
dc.mt.gsea() ranks features using a continuous statistic and computes an enrichment score from the running sum over that ranking. When times > 1, decoupler also performs empirical testing and adjusted p-value correction.
Step 6: Extract the score matrix
decoupler stores results in obsm, so the next step is to pull them into a standalone AnnData object for plotting.
gsea = dc.pp.get_obsm(adata, key="score_gsea")
gsea
gsea.to_df().iloc[:5, :5]
dc.pp.get_obsm() converts the selected obsm entry into a new AnnData, which lets you reuse standard Scanpy plotting and analysis functions on the pathway score matrix.
Step 7: Visualize pathway activity on UMAP
A quick way to explore results is to plot the most variable pathways across cells.
top_var = (
gsea.to_df()
.var()
.sort_values(ascending=False)
.head(4)
.index
.tolist()
)
top_var
sc.pl.umap(gsea, color=top_var, cmap="RdBu_r", ncols=2)
This produces pathway-level maps instead of gene-level maps, which is often easier to interpret when comparing immune subpopulations.
Step 8: Find pathways that define each cell type
You can rank pathways by group to see which signatures are most characteristic of each annotated cell type.
ranked = dc.tl.rankby_group(
adata=gsea,
groupby="celltype",
reference="rest",
method="wilcoxon",
)
ranked.head()
Keep the strongest positive pathways per cell type:
top_pathways = (
ranked[ranked["stat"] > 0]
.sort_values(["group", "stat"], ascending=[True, False])
.groupby("group")
.head(3)
)
top_pathways[["group", "name", "stat", "padj"]].head(15)
Plot the selected pathways as a matrix plot:
plot_terms = top_pathways["name"].drop_duplicates().tolist()
sc.tl.dendrogram(gsea, groupby="celltype")
sc.pl.matrixplot(
gsea,
var_names=plot_terms,
groupby="celltype",
standard_scale="var",
cmap="RdBu_r",
dendrogram=True,
)
dc.tl.rankby_group() compares pathway scores across groups and returns a table with statistics and adjusted p-values, which makes it useful for turning a large score matrix into a compact list of interpretable pathway markers.
Step 9: Adapt the workflow to your own data
For your own dataset, the workflow is almost identical:
- Start from a Scanpy
AnnData - Make sure expression is normalized and log-transformed
- Choose a gene set resource
- Run
dc.mt.gsea(...) - Extract
score_gsea - Plot and compare scores by cluster, cell type, condition, or pseudotime
Minimal template:
# adata = your Scanpy AnnData object
gene_sets = dc.op.hallmark(organism="human")
gene_sets = dc.pp.prune(adata.var_names, gene_sets, tmin=5)
dc.mt.gsea(data=adata, net=gene_sets, tmin=5, times=1)
scores = dc.pp.get_obsm(adata, key="score_gsea")
Recap
You now have a complete pattern for pathway analysis with Scanpy and decoupler: load an AnnData object, fetch a pathway resource, run GSEA, extract the score matrix, and visualize pathway activity per cell and per group. For unweighted gene sets, GSEA is a strong first choice; for weighted resources, methods such as ULM or MLM are often a better fit because GSEA does not model edge weights.
Further Reading
- decoupler single-cell enrichment tutorial
- decoupler GSEA API
- decoupler Hallmark gene sets API
- Scanpy PBMC3k dataset docs
FAQ
1. Can I run this on raw counts?
Not directly. decoupler methods assume appropriately normalized input for observational count data, and the docs explicitly recommend normalization such as library-size normalization followed by log1p. Raw integer counts are not recommended for GSEA input. (decoupler.readthedocs.io)
2. Why are some pathways missing from my result?
Usually this means there was not enough overlap between your gene identifiers and the gene set resource. Make sure species and gene symbols match, and use dc.pp.prune() so only pathways with enough shared targets are retained.
3. Should I score each cell or use pseudobulk?
Per-cell scoring is great for exploratory analysis and visualization. But when your experiment has multiple samples or conditions, the recommended approach is often pseudobulk analysis so biological replication is handled at the sample level rather than treating cells as independent replicates. (decoupler.readthedocs.io)
