Pseudobulk GSEA with limma contrasts and per-patient ssGSEA

The statistically rigorous alternative to RunGSEA(). Instead of treating individual cells as replicates, this function aggregates counts to the donor/sample level (pseudobulk), fits a limma-voom linear model, and ranks genes by the moderated t-statistic before running fgsea. This respects the true unit of replication and avoids the inflated degrees of freedom that arise from cell-level comparisons.

Optionally also computes per-patient pathway activity scores (ssGSEA via GSVA) and tests them across groups, producing boxplots and heatmaps alongside the GSEA results.

Usage

RunGSEA_pseudobulk(
  seurat_object = NULL,
  group.by = NULL,
  sample.by = NULL,
  contrast.by = NULL,
  split.by = NULL,
  subset.by = NULL,
  run_label = NULL,
  label.by = NULL,
  pseudobulk = NULL,
  covariates = NULL,
  contrasts = 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,
  identity_column = NULL,
  sample_column = NULL,
  group_column = NULL,
  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 = "human",
  output_dir,
  assay = NULL,
  min.cells = 10L,
  min.samples.per.group = 3L,
  run.ssgsea = TRUE,
  ssgsea.use.padj = TRUE,
  gsea.use.padj = TRUE,
  top_n_heatmap = 30L,
  top_n_ssgsea = 30L,
  nes.cutoff = 1,
  group_colors = NULL,
  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,
  resume = FALSE,
  resume_folder = NULL,
  timestamp = TRUE,
  save.rds = TRUE,
  redo_plots = FALSE,
  caffeinate = FALSE
)

Arguments

seurat_object

A Seurat object. Either seurat_object or a precomputed pseudobulk list must be supplied.

group.by

Character or NULL. Metadata column holding cell-type labels. One independent GSEA is run per level (e.g. "Assignment"). Pass NULL when the object already contains only one cell type (e.g. a Microglia-only object); all cells are treated as a single group called "All" and a single GSEA is run.

sample.by

Character. Metadata column identifying the biological replicate (donor/patient ID). This is the statistical unit; each unique value becomes one pseudobulk sample in the limma-voom model (e.g. "DonorID").

contrast.by

Character. Metadata column holding the comparison group label (e.g. "Disease" with levels "AD" and "Control"). Cells with NA in this column are excluded before aggregation.

split.by

Character or NULL. Optional metadata column used to create interaction groups before fitting the model. When supplied, the contrast groups become <contrast.by>.<split.by> combinations (e.g. split.by = "Sex" turns "AD" and "Control" into "AD.Female", "AD.Male", "Control.Female", "Control.Male"), allowing sex-stratified or condition-stratified contrasts. Default NULL.

subset.by

Character or NULL. Optional metadata column used to run the full pipeline separately for each level. Each level gets its own subdirectory under output_dir. Use this to produce independent analyses per broad cell class, brain region, or any other grouping variable without calling the function multiple times. Default NULL.

run_label

Character or NULL. Short label appended to output filenames, plot titles, AND (when timestamp = TRUE) the run's subfolder name, to distinguish this run from others saved to the same output_dir (e.g. run_label = "with_age_covariate"). Set automatically when subset.by is used; otherwise leave NULL unless you need a custom tag.

pseudobulk

Optional precomputed pseudobulk. A named list with elements $counts (genes x samples matrix) and $samples (sample-level metadata data frame). Supply the value returned by a previous RunGSEA_pseudobulk() call to skip the aggregation step on reruns. Default NULL (aggregation is performed from seurat_object).

covariates

Character vector of sample-level metadata columns to include as covariates in the limma-voom model (e.g. c("Age", "PMI")). Numeric covariates are median-imputed for NA; factor/character covariates receive an "unknown" level. Default NULL.

contrasts

Named list of limma contrast expressions. Every element must be a single string built from model coefficient names of the form make.names(paste0("group", level)). A small helper makes this concise:


  cf <- function(l) make.names(paste0("group", l))
  contrasts <- list(
    AD_vs_Ctrl   = paste(cf("AD"),     "-", cf("Control")),
    Interaction  = paste0("(", cf("AD.F"), "-", cf("Ctrl.F"), ")",
                          " - ",
                          "(", cf("AD.M"), "-", cf("Ctrl.M"), ")")
  )

Positive NES = enrichment in the first term. When NULL a standard set is generated automatically from the levels of contrast.by. When redo_plots = TRUE, a user-supplied list takes priority over the list saved in the RDS. Default NULL.

gene_sets

Named list of character vectors. Each element is one gene set (e.g. list(MySet = c("GENE1", "GENE2", ...))). When supplied, pathway_sets and deg_df are ignored and msigdbr is not required. Default NULL.

deg_df

Data frame of differential expression results. Top DEGs per group are extracted and used as gene sets. When supplied, pathway_sets is ignored and msigdbr is not required. Use the deg_* parameters below to map your data frame's column names. Default NULL.

deg_gene_column

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

deg_group_column

Character (length 1 or 2). Column(s) in deg_df defining the group each row belongs to. Two columns create a two-level label used for heatmap row splitting. Default "group".

deg_fc_column

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

deg_padj_column

Character or NULL. Column in deg_df containing adjusted p-values for pre-filtering. Set to NULL to skip filtering. Default "padj".

deg_padj_cutoff

Numeric. Adjusted p-value threshold applied to deg_df before selecting top genes. Default 0.05.

deg_top_n

Integer. Number of top genes per group (ranked by deg_fc_column, descending) to include in each gene set. Default 20.

pathway_sets

Named list of MSigDB collections to test. Each element must have a category field and an optional subcategory field (see msigdbr::msigdbr_collections()). The element name becomes the label used in output filenames. Default: Hallmark, KEGG, Reactome, WikiPathways. Ignored when gene_sets or deg_df is supplied.

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 name passed to msigdbr::msigdbr(). Use "human" or "mouse" (shorthand accepted) or any full name returned by msigdbr::msigdbr_species(). Default "human".

output_dir

Character. Root directory for all CSV, PDF, and RDS output. Per-contrast and per-database subdirectories are created automatically.

assay

Character or NULL. Seurat assay containing raw counts to aggregate. NULL uses the default assay. For BPCells sketch workflows, always pass "RNA" explicitly since pseudobulk aggregation requires raw counts from all cells. Default NULL.

min.cells

Integer. Minimum number of cells a sample must contribute to a given cell type to be included in the pseudobulk for that cell type. Samples below this threshold are excluded to avoid noisy pseudo-replicates. Default 10.

min.samples.per.group

Integer. Minimum number of samples required per group level within a cell type. Cell types where any group has fewer than this many samples are skipped entirely. Default 3.

run.ssgsea

Logical. Also compute per-patient pathway activity via single-sample GSEA (ssGSEA, using GSVA::gsva())? Produces score heatmaps and contrast boxplots per cell type per database. Requires the GSVA Bioconductor package. Default TRUE.

ssgsea.use.padj

Logical. When TRUE (default), pathway significance in ssGSEA is judged by BH-adjusted p-value (padj < 0.05). Set to FALSE to use the raw p-value instead – useful when many contrasts make BH correction too aggressive and real signals drop out. Affects pathway selection for heatmaps, boxplot panel selection, and significance bracket thresholds. Default TRUE.

gsea.use.padj

Logical. Same choice as ssgsea.use.padj, but for the main fgsea NES heatmaps: which column (padj vs. raw pval) drives the significance stars (* / ** / ***) drawn in each pathway x cell-type (or pathway x contrast, for the summary heatmap) cell. Default TRUE.

top_n_heatmap

Integer. Number of pathways shown in the NES heatmap per contrast/database combination, selected by maximum |NES| across cell types. Default 30.

top_n_ssgsea

Integer. Number of pathways shown in the ssGSEA score heatmap per cell type, selected by significance or variance. Default 30.

nes.cutoff

Numeric. Minimum absolute NES required for a pathway to appear in the summary NES heatmap (pathways x contrasts). The candidate row set is first built as the union of the top top_n_heatmap pathways from each individual per-contrast heatmap; nes.cutoff then removes rows whose maximum |NES| across all contrasts falls below this threshold. Default 1.0. Set to 0 to skip the cutoff and show all union rows.

group_colors

Named character vector mapping contrast.by levels to colors. Used consistently across all boxplots and heatmap annotations. NULL auto-generates colors from the Nour palette. Default NULL.

heatmap_params

Named list of additional arguments forwarded to ComplexHeatmap::Heatmap(). Any default set internally can be overridden here (e.g. list(clustering_distance_rows = "pearson")).

heatmap_colors

A circlize::colorRamp2 color function for all NES and ssGSEA score heatmaps. NULL (default) uses a 5-stop blue-white-red diverging palette centered at 0. Breaks are fixed at c(-2, -1, 0, 1, 2) for ssGSEA row-z-score heatmaps and auto-scaled to the data for NES heatmaps. Supply any colorRamp2 object to override all heatmaps at once. Examples:

Classic RdBu: circlize::colorRamp2(c(-3, 0, 3), c("#2166ac", "white", "#b2182b"))
Purple-green (PRGn): circlize::colorRamp2(c(-3, 0, 3), RColorBrewer::brewer.pal(3, "PRGn"))
Viridis plasma (for one-sided scores): circlize::colorRamp2(seq(0, 3, length.out = 9), viridis::plasma(9))
Tighter scale for subtle signals: circlize::colorRamp2(c(-1.5, -0.75, 0, 0.75, 1.5), c("#007dd1", "#b3d9f5", "white", "#f5c08a", "#ab3000"))

width, height

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

resume

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

FALSE (default): Always start a fresh, new timestamped subfolder (see timestamp).
TRUE: Scan output_dir for previous RunGSEA_pseudobulk_* subfolders whose call-time parameters (group.by, sample.by, contrast.by, covariates, gene sets, etc. - NOT purely cosmetic ones like heatmap_colors) match this call exactly. Zero matches falls back to a fresh subfolder; one match resumes it (skipping recomputation the same way pseudobulk reuse already does); 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, cell types whose GSEA CSV files already exist are additionally skipped (as before) so a crashed run can continue.

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. Useful in scripts where you already know which previous run to continue.

timestamp

Logical. When TRUE (default), every fresh call gets its own subfolder named RunGSEA_pseudobulk[_run_label]_<YYYYMMDD-HHMMSS> under output_dir, so re-running with different settings (e.g. with vs. without covariates) into the same output_dir never overwrites a previous run's files or methods JSON. Set FALSE to use output_dir directly (legacy behavior) - in that mode you are responsible for giving different configurations different output_dir values yourself.

save.rds

Logical. Save the pseudobulk counts, sample metadata, GSEA results, and ssGSEA results as RunGSEA_pseudobulk_results.rds in output_dir. Required for redo_plots = TRUE. Default TRUE.

redo_plots

Logical. Reload a previously saved RDS and regenerate all plots without rerunning the full pipeline. Default FALSE. When TRUE:

RunGSEA_pseudobulk_results.rds is loaded from output_dir.
NES heatmaps and lollipop plots are regenerated from the saved GSEA data frame.
ssGSEA score CSVs on disk are reread and used to regenerate score heatmaps and boxplots.
seurat_object, group.by, sample.by, and contrast.by are not needed and can be omitted.

Useful after changing group_colors, heatmap_params, or any plot parameter without rerunning the hours-long pipeline.

caffeinate

Logical. Prevent the system from sleeping during the run (macOS only, uses the built-in caffeinate utility). Recommended for overnight runs. Default FALSE.

Value

Invisibly returns a named list:

pseudobulk: Genes x samples aggregated count matrix.
samples: Sample-level metadata data frame used in the model.
gsea: Long-format data frame of all fgsea results with columns pathway, NES, padj, cell_type, contrast, and db.
ssgsea: ssGSEA contrast results (one row per pathway/cell type/ contrast), or NULL if run.ssgsea = FALSE.
contrasts: The named contrast list actually used in the model (useful when contrasts were auto-generated).

Three required columns

Three metadata columns drive the analysis:

group.by: Cell type column. GSEA is run independently for each level (e.g. "Assignment"). Cells are aggregated within each level.
sample.by: Donor/sample ID column. This is the statistical unit. Each unique value becomes one pseudobulk replicate in the limma model (e.g. "DonorID", "SampleID").
contrast.by: Comparison group column. Levels define the groups being contrasted (e.g. "Disease" with levels "AD" and "Control"). Cells with NA in this column are dropped before aggregation.

Gene set input modes

Supply gene sets using exactly one of the following (priority order):

gene_sets - a named list of character vectors, one vector per gene set.
deg_df - a differential expression data frame; top DEGs per group are extracted and used as gene sets (see deg_* parameters).
pathway_sets (default) - one or more MSigDB collections fetched via msigdbr. Does not require gene_sets or deg_df.