Run GSEA pathway analysis on a Seurat object

Performs the following steps for each combination of label.by (outer loop) and split.by (inner loop) levels:

Subset the object to the current level.
Run Wilcoxon rank-sum DE with presto::wilcoxauc().
For each pathway database in pathway_sets, run fgsea::fgsea() using AUC scores as ranking statistics.
Generate a lollipop NES plot (top/bottom 10 pathways per cluster).
Build NES heatmaps showing top top_n pathways per cluster.

Usage

RunGSEA(
  seurat_object,
  group.by = "Assignment",
  split.by = "GlobalAssignment",
  label.by = NULL,
  gene_sets = NULL,
  deg_df = NULL,
  deg_gene_column = "feature",
  deg_group_column = "group",
  deg_fc_column = "logFC",
  deg_padj_column = "padj",
  deg_padj_cutoff = 0.05,
  deg_top_n = 20L,
  pathway_sets = list(Hallmark = list(category = "H"), KEGG = list(category = "C2",
    subcategory = "CP:KEGG"), Reactome = list(category = "C2", subcategory =
    "CP:REACTOME"), WP = list(category = "C2", subcategory = "CP:WIKIPATHWAYS")),
  search_terms = NULL,
  species = "Homo sapiens",
  assay = "RNA",
  min_cells = 20L,
  padj_thresh = 0.05,
  logfc_thresh = 0.1,
  top_n = 5L,
  nes.cutoff = 1,
  output_dir = NULL,
  resume = FALSE,
  resume_folder = NULL,
  timestamp = TRUE,
  heatmap_params = list(row_names_side = "left", show_row_dend = FALSE,
    row_names_max_width = grid::unit(15, "cm")),
  heatmap_colors = NULL,
  width = NULL,
  height = NULL,
  caffeinate = FALSE
)

Arguments

seurat_object

A Seurat object.

group.by

Character. Metadata column whose levels are compared within each subset (used as DE groups). Default "Assignment".

split.by

Character or NULL. Metadata column whose levels each get their own independent GSEA (e.g., run GSEA separately per broad cell type such as "GlobalAssignment"). NULL runs a single GSEA on the whole object (or on each label.by subset), comparing all group.by levels without any further cell-type subsetting. Default "GlobalAssignment".

label.by

Character or NULL. Outer metadata column for further subsetting (e.g., condition). NULL treats the whole object as one group.

gene_sets

Named list of character vectors, where each name is a gene set label and each vector contains gene symbols. When supplied, MSigDB is not queried and pathway_sets is ignored. Takes priority over deg_df. Default NULL.

deg_df

Data frame of differential expression results. When supplied (and gene_sets is NULL), gene sets are built by taking the top deg_top_n up- and down-regulated genes per group from this table. Columns are identified by deg_gene_column, deg_group_column, deg_fc_column, and deg_padj_column. Default NULL.

deg_gene_column

Character. Column in deg_df containing gene symbols. Default "feature".

deg_group_column

Character. Column in deg_df containing the group label used to split genes into per-group sets. Default "group".

deg_fc_column

Character. Column in deg_df containing log-fold change values used to rank genes within each group. Default "logFC".

deg_padj_column

Character. Column in deg_df containing adjusted p-values used to filter genes before ranking. Default "padj".

deg_padj_cutoff

Numeric. Adjusted p-value threshold applied to deg_df rows before building gene sets. Default 0.05.

deg_top_n

Integer. Maximum number of genes per group taken from deg_df after filtering and ranking. Default 20.

pathway_sets

Named list of MSigDB gene-set databases to test. Each element is a list with a category field and an optional subcategory field; the element name becomes the short label used in output filenames and the return list. The default covers four curated collections:

list(
  Hallmark = list(category = "H"),
  KEGG     = list(category = "C2", subcategory = "CP:KEGG"),
  Reactome = list(category = "C2", subcategory = "CP:REACTOME"),
  WP       = list(category = "C2", subcategory = "CP:WIKIPATHWAYS")
)

Other useful collections (pass as additional list elements):

list(category = "C2", subcategory = "CP:BIOCARTA") - BioCarta
list(category = "C2", subcategory = "CP:PID") - NCI Pathway Interaction Database
list(category = "C5", subcategory = "GO:BP") - Gene Ontology Biological Process
list(category = "C5", subcategory = "GO:MF") - GO Molecular Function
list(category = "C5", subcategory = "GO:CC") - GO Cellular Component
list(category = "C7", subcategory = "IMMUNESIGDB") - ImmuneSigDB (immune gene sets; recommended for immune datasets)
list(category = "C8") - Cell-type signature gene sets

See msigdbr::msigdbr_collections() for the full catalog.

search_terms

Optional keyword filter for gene-set names (case- insensitive, regex), applied after fetching each database in pathway_sets (or to the names of a custom gene_sets list). Two forms: a character vector is an OR search where any set matching at least one term is kept (e.g. c("APOPTOSIS", "CELL_DEATH")); a list of character vectors is a union of AND-groups where all terms in an element must appear in the same name (e.g. list(c("T_CELL", "ACTIVATION"), "APOPTOSIS")). Ignored when deg_df is supplied. NULL (default) tests all gene sets.

species

Character. Species passed to msigdbr::msigdbr(). Must match the species names recognized by that function. Common values:

"Homo sapiens" - human (default)
"Mus musculus" - mouse
"Rattus norvegicus" - rat
"Danio rerio" - zebrafish

Run msigdbr::msigdbr_species() for the complete list.

assay

Character. Seurat assay for DE via presto::wilcoxauc(). Default "RNA". For BPCells sketch workflows use "sketch".

min_cells

Integer. Skip subsets with fewer cells. Default 20. A separate guard also skips any (label, split) combination that has fewer than 2 levels of group.by after subsetting - this prevents presto::wilcoxauc from crashing when e.g. one label.by level contains only one group.by category (e.g., a Cognitive.Status group with a single Sex). A message is emitted for each skipped combination.

padj_thresh

Numeric. Adjusted p-value threshold for DE and fgsea filtering. Default 0.05.

logfc_thresh

Numeric. Minimum log-fold-change for DE filtering. Default 0.1.

top_n

Integer. Number of top and bottom pathways per cluster shown in summary heatmaps. Default 5.

nes.cutoff

Numeric. After selecting the union of top/bottom top_n pathways per cluster, rows whose maximum |NES| across all clusters falls below this threshold are removed from the summary heatmap. Set to 0 to skip the filter. Default 1.0.

output_dir

Character or NULL. If provided, CSVs and PDFs are written here in a structured subdirectory tree. If NULL, results are only returned as an R list.

resume

Logical or "last". Controls both WHICH subfolder this call uses and whether cached work is reused:

FALSE (default): Always start a fresh, new timestamped subfolder (see timestamp).
TRUE: Scan output_dir for previous RunGSEA_* subfolders whose call-time parameters (group.by, split.by, species, gene sets, etc.) match this call exactly. Zero matches falls back to a fresh subfolder; one match resumes it; several matches prompt interactively (utils::menu()) when the session is interactive, or stop() with instructions otherwise.
"last": Same scan, but silently resumes the most recent match instead of asking.

Within a resumed subfolder, work already completed is additionally skipped:

DE cache: Wilcoxon DE results for each (label, split) pair are saved to output_dir/label/de_cache_<split>_<label>.rds and reloaded on subsequent runs instead of re-running presto::wilcoxauc.
DB-level skip: If the full NES heatmap PDF for a (label, split, db_name) triple already exists, that database is skipped entirely (its slot in the return list is set to NULL).

resume_folder

Character or NULL. Explicitly resume this exact subfolder (a bare folder name under output_dir, or a full path), skipping the scan/match described under resume.

timestamp

Logical. When TRUE (default), every fresh call gets its own subfolder named RunGSEA_<YYYYMMDD-HHMMSS> under output_dir, so re-running with different settings into the same output_dir never overwrites a previous run's files or methods JSON. Set FALSE to use output_dir directly (legacy behavior).

heatmap_params

Named list of additional arguments forwarded to ComplexHeatmap::Heatmap() for both the full and summary NES heatmaps. Any default set by scSidekick can be overridden here - e.g. clustering distance, font size, or row ordering.

heatmap_colors

A circlize::colorRamp2 color function for the heatmap fill scale. NULL (default) uses a 5-stop blue-white-red diverging palette where the limits are auto-scaled to the data range (capped at +/-4):


  circlize::colorRamp2(
    c(-lim, -lim/2, 0, lim/2, lim),
    c("#007dd1", "#b3d9f5", "white", "#f5c08a", "#ab3000")
  )

Supply any colorRamp2 object to override. The breaks you set in the custom function are used as-is (no auto-scaling). Examples:

Classic RdBu (3-stop): circlize::colorRamp2(c(-3, 0, 3), c("#2166ac", "white", "#b2182b"))
Viridis plasma (sequential, for one-sided scores): circlize::colorRamp2(seq(0, 3, length.out = 9), viridis::plasma(9))
Purple-green (PRGn, for a cooler look): circlize::colorRamp2(c(-3, 0, 3), RColorBrewer::brewer.pal(3, "PRGn"))
Same palette, tighter breaks (good for subtle signals): circlize::colorRamp2(c(-1.5, -0.75, 0, 0.75, 1.5), c("#007dd1", "#b3d9f5", "white", "#f5c08a", "#ab3000"))

Note: for NES heatmaps, diverging palettes (centered at 0) are semantically correct because 0 means "no enrichment".

width, height

Numeric or NULL. Override the auto-calculated PDF dimensions (in inches) applied to every NES heatmap this run saves. NULL (default) auto-sizes each heatmap from its own row/column count and label lengths.

caffeinate

Logical. If TRUE, prevents the Mac from sleeping during the run via caffeinate. Useful for long multi-database runs overnight. Default FALSE.

Value

A nested named list structured as results[[label]][[split]][[db_name]] with elements:

de_table: Full DE results from presto::wilcoxauc.
gsea_table: Full fgsea NES table (all pathways, all clusters).
nes_matrix: Wide NES matrix (pathways × clusters).
lollipop_plots: Named list of ggplot2 lollipop plots per cluster.

Details

No setwd() calls are made. All output paths are constructed from output_dir.