Skip to contents

Single-cell ssGSEA with significance testing and heatmap visualization

RunSCssGSEA.Rd

Runs per-cell single-sample Gene Set Enrichment Analysis (ssGSEA) on a Seurat object using GSVA (compatible with both GSVA v1 and v2 APIs), computes group-level significance (ANOVA or Kruskal-Wallis), saves enrichment scores and significance tables as CSVs, and produces a ComplexHeatmap PDF of mean z-scored enrichment per group.

Usage

RunSCssGSEA(
  seurat_object,
  group.by,
  split.by = NULL,
  subset_cells = FALSE,
  subset_by = NULL,
  subset_values = 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 = 20,
  gene_set_library = NULL,
  gene_set_subcategory = NULL,
  pathway_sets = NULL,
  search_terms = NULL,
  species = "Homo sapiens",
  method = "gsva",
  assay = "RNA",
  downsample = NULL,
  min.size = 5,
  ssgsea.norm = TRUE,
  ucell_max_rank = 1500L,
  cores = 1,
  chunk_size = 1000L,
  resume = FALSE,
  fit = "ANOVA",
  sig_group_by = NULL,
  p_cutoff = 0.05,
  use.padj = TRUE,
  show_only_significant = TRUE,
  pathways = NULL,
  group_colors = NULL,
  split_colors = NULL,
  heatmap_params = list(row_names_side = "left", show_row_dend = FALSE,
    row_names_max_width = grid::unit(15, "cm")),
  heatmap_colors = NULL,
  add_boxplots = TRUE,
  add_pvalues = TRUE,
  add_to_object = FALSE,
  output_dir,
  resume_folder = NULL,
  timestamp = TRUE,
  object_name = "Analysis",
  subset_name = "",
  caffeinate = FALSE
)

Arguments

seurat_object

A Seurat object.

group.by

Character. Metadata column used as the primary grouping variable for the heatmap columns and significance testing (e.g. "Assignment", "celltype").

split.by

Character or NULL. Optional secondary metadata column for splitting heatmap columns (e.g. "Health", "Condition"). When set, significance is also tested within each split level in addition to the pooled test.

subset_cells

Logical. If TRUE, subset the Seurat object to subset_values before scoring. Default FALSE.

subset_by

Character or NULL. Metadata column used for subsetting (required when subset_cells = TRUE). NULL treats subset_values as cell barcodes directly.

subset_values

Character vector. Values to retain in subset_by, or cell barcodes when subset_by = NULL.

gene_sets

Named list of gene vectors (Mode 1). NULL if using another mode. When search_terms is also supplied, it filters these custom names too (matching RunGSEA()/RunGSEA_pseudobulk() behavior) rather than being silently ignored.

deg_df

Data frame of DEG results (Mode 2). NULL otherwise. Compatible with standard outputs from presto::wilcoxauc(), edgeR, or any table with gene, group, and log-fold-change columns.

deg_gene_column

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

deg_group_column

Character scalar or length-2 character vector. Column(s) in deg_df defining the gene-set groups. One column → one gene set per unique value. Two columns → gene sets named "col1__NKSEP__col2"; heatmap rows split by col1 automatically. Default "group".

deg_fc_column

Character. Log-fold-change column for selecting top genes. Default "logFC".

deg_padj_column

Character or NULL. Adjusted p-value column for pre- filtering. NULL skips filtering. Default "padj".

deg_padj_cutoff

Numeric. Adjusted p-value threshold for DEG pre-filtering. Default 0.05.

deg_top_n

Integer. Top N genes per group (by deg_fc_column) to include in each gene set. Default 20.

gene_set_library

Character or NULL. MSigDB collection code for a single collection (Mode 3 or 4). Ignored when pathway_sets is supplied. Common options:

  • "H" - Hallmark (50 coherent biological processes)

  • "C2" - Curated: canonical pathways (combine with gene_set_subcategory = "CP:KEGG", "CP:REACTOME", "CP:WIKIPATHWAYS", "CP:BIOCARTA", "CP:PID")

  • "C5" - Gene Ontology (subcategory = "GO:BP", "GO:MF", "GO:CC")

  • "C6" - Oncogenic signatures

  • "C7" - Immunologic signatures (IMMUNESIGDB)

  • "C8" - Cell-type signatures

NULL (and pathway_sets also NULL) fetches all collections (slow; use with search_terms).

gene_set_subcategory

Character or NULL. MSigDB sub-collection code, paired with gene_set_library. Ignored when pathway_sets is supplied. Default NULL (entire collection).

pathway_sets

Named list or NULL. Test multiple MSigDB collections at once - the same argument and format RunGSEA() and RunGSEA_pseudobulk() use, e.g. list(Hallmark = list(category = "H"), KEGG = list(category = "C2", subcategory = "CP:KEGG")). All listed collections are fetched and pooled into a single set of gene sets (RunSCssGSEA always scores everything into one heatmap/boxplot set, unlike RunGSEA's one-heatmap-per-database output). When supplied, this takes priority over gene_set_library / gene_set_subcategory. NULL (default) uses those two instead.

search_terms

Search filter applied to MSigDB gene-set names (case- insensitive, regex). Accepts two forms:

Character vector - OR search. Any gene set whose name matches at least one term is kept: search_terms = c("APOPTOSIS", "CELL_DEATH").

List of character vectors - AND/OR search. Each list element is an AND-group: all terms in that element must appear in the same pathway name (order and separation do not matter). Results are the union (OR) across AND-groups. Examples:

  • list(c("T_CELL", "ACTIVATION")) - both words in same name

  • list(c("T_CELL", "ACTIVATION"), "APOPTOSIS") - previous OR APOPTOSIS

  • list(c("T_CELL","ACTIVATION"), c("B_CELL","DIFFERENTIATION")) - two independent AND-groups, union of results

NULL (default) keeps all gene sets in the fetched collection.

species

Character. Species for MSigDB fetching. Passed to msigdbr::msigdbr(). Default "Homo sapiens". Use "Mus musculus" for mouse datasets. See msigdbr::msigdbr_species() for all supported species.

method

Character. Scoring backend. "gsva" (default) uses the GSVA package's ssGSEA algorithm; "ucell" uses the UCell package's rank-based U statistic, which is faster and memory-efficient for large datasets but does not produce normalized enrichment scores.

assay

Character. Seurat assay to pull normalized ("data" layer) expression from. Default "RNA". BPCells/sketch-aware: pass "sketch" to score a Seurat v5 sketch-subsampled assay for large datasets instead of the full "RNA" assay.

downsample

Integer or NULL. Maximum cells per group.by identity to use for scoring. ssGSEA is \(O(n_{\text{cells}})\) - downsampling to 500-2000 cells per group dramatically reduces runtime with minimal loss of group-level signal. NULL uses all cells. Default NULL.

min.size

Integer. Minimum gene-set size after intersecting with the expression matrix. Gene sets below this threshold are silently dropped. Default 5.

ssgsea.norm

Logical. Normalize ssGSEA enrichment scores to the [0, 1] range? Default TRUE (matches escape::enrichIt default). Ignored when method = "ucell".

ucell_max_rank

Integer. UCell only: only the top ucell_max_rank genes per cell (by expression rank) are considered when computing the U statistic. Lower values reduce memory use. Default 1500L.

cores

Integer. Number of parallel cores. For method = "gsva", requires BiocParallel; falls back to serial if unavailable. For method = "ucell", sets the number of UCell threads. Default 1.

chunk_size

Integer. UCell only: number of cells processed per parallel chunk. Tune to balance memory and speed. Default 1000L.

resume

Logical. If TRUE and a combined-scores CSV checkpoint from a previous run exists in output_dir, load it and skip scoring. Useful for re-running plots after adjusting significance thresholds. Default FALSE.

fit

Character. Significance test applied per pathway: "ANOVA" (default, parametric F-test) or "Kruskal" (non-parametric Kruskal-Wallis). Both use BH adjustment for multiple testing.

sig_group_by

Character or NULL. Metadata column to use as the grouping factor for significance testing. NULL (default) uses group.by. Override when your display grouping (group.by) differs from your statistical grouping (e.g. display by fine cell type, test by broad disease category).

p_cutoff

Numeric. Significance threshold for selecting pathways to display in the heatmap - compared against BH-adjusted or raw p-value depending on use.padj. Default 0.05.

use.padj

Logical. When TRUE (default), p_cutoff is compared against the BH-adjusted p-value (p_adj). Set to FALSE to use the raw p-value instead – useful when many pathways make BH correction too aggressive and real signals drop out. Affects pathway selection for the heatmap and the significance wording in legend text.

show_only_significant

Logical. If TRUE, heatmap shows only pathways passing p_cutoff (see use.padj). Falls back to top 30 by statistic if no pathways pass the threshold. If FALSE, all pathways are shown. Ignored when pathways is supplied. Default TRUE.

pathways

Character vector or NULL. Force these specific pathway names into the heatmap/boxplots regardless of significance, bypassing the p_cutoff filter entirely (e.g. a pathway of interest that didn't clear the threshold). Names not found among the tested gene sets are dropped with a warning. Default NULL (use the significance-based selection).

group_colors

Named character vector of colors for group.by levels. NULL (default) auto-generates colors from the Nour palette.

split_colors

Named character vector of colors for split.by levels. NULL (default) auto-generates colors from the Nour palette. Ignored when split.by = NULL.

heatmap_params

Named list of additional arguments forwarded to ComplexHeatmap::Heatmap(). Overrides internal defaults for row label placement, dendrogram display, and row-name width.

heatmap_colors

A circlize::colorRamp2 color function for the heatmap fill scale. NULL (default) uses a 5-stop blue-white-red diverging palette at breaks c(-2, -1, 0, 1, 2) (row z-scores):


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

Supply any colorRamp2 object to override. The breaks you set are used as-is - for row z-score heatmaps, symmetric breaks around 0 are recommended. Examples:

  • Classic RdBu (3-stop): circlize::colorRamp2(c(-2, 0, 2), c("#2166ac", "white", "#b2182b"))

  • Purple-green (PRGn): circlize::colorRamp2(c(-2, 0, 2), RColorBrewer::brewer.pal(3, "PRGn"))

  • Viridis plasma (for raw non-z-scored scores): circlize::colorRamp2(seq(0, 1, length.out = 9), viridis::plasma(9))

  • Tighter scale for subtle signals: circlize::colorRamp2(c(-1, -0.5, 0, 0.5, 1), c("#007dd1", "#b3d9f5", "white", "#f5c08a", "#ab3000"))

add_boxplots

Logical. If TRUE (default), generate per-pathway boxplot pages showing per-cell score distributions across group.by groups (and split by split.by if set).

add_pvalues

Logical. If TRUE (default), add pairwise t-test p-value brackets to boxplots using ggpubr::stat_compare_means(). These are raw (unadjusted) p-values for each displayed pair - not BH-corrected, and unrelated to the use.padj/p_cutoff pathway-selection threshold above. Requires add_boxplots = TRUE.

add_to_object

Logical. If FALSE (default), returns the enrichment scores data frame invisibly and does not modify the Seurat object (scores are saved to CSV only - recommended for large datasets to avoid metadata bloat). If TRUE, adds one metadata column per gene set to seurat_object@meta.data and returns the modified Seurat object. When downsample is also set, cells not included in the downsampled scoring set receive NA.

output_dir

Character. Directory for output files (created if absent).

object_name

Character. First component of output file name prefix. Default "Analysis".

subset_name

Character. Second component of output file name prefix. Default "".

caffeinate

Logical. If TRUE, prevents the Mac from sleeping during the run using the caffeinate system command. Default FALSE.

Value

If add_to_object = FALSE (default): invisibly returns a data frame of enrichment scores (cells × gene sets) with additional columns for group.by and split.by identity.

If add_to_object = TRUE: invisibly returns the Seurat object with one metadata column per gene set.

Files saved to output_dir: written into a dedicated subfolder named SC ssGSEA <gene-set label> (e.g. SC ssGSEA STAT3, SC ssGSEA custom) - this subfolder name is how to recognize a figure was made by RunSCssGSEA() specifically (not RunGSEA_pseudobulk()'s own run.ssgsea = TRUE sub-step, which does not use this naming convention - see RunGSEA_pseudobulk()) when working backward from a saved file without its analysis_params.json.

  • <prefix> ssGSEA_scores.csv - per-cell enrichment scores + group labels

  • <prefix> ssGSEA_significance.csv - per-pathway significance (pooled across all cells)

  • <prefix> ssGSEA_significance_by_<split.by>.csv - per-pathway significance within each split.by level (only when split.by is set)

  • <prefix> ssGSEA heatmap.pdf - ComplexHeatmap of mean z-scored enrichment per group

Details

Unlike pseudobulk ssGSEA (see RunGSEA_pseudobulk()), this function scores every cell individually rather than aggregating by sample first. It is appropriate when you have no meaningful biological replicates (sample IDs), or when you want to use the single-cell resolution for visualization. For datasets with proper replicates, pseudobulk ssGSEA is statistically more rigorous.

Gene-set input (choose exactly one mode)

Mode 1 - pre-built named list

gene_sets = list(SetA = c("GENE1", "GENE2"), SetB = c("GENE3"))

Mode 2 - DEG data frame
Converts a marker/DEG table into one gene set per group. The deg_group_column argument can name one column (simple group, e.g. "celltype") or two columns (e.g. c("celltype", "cytokine")). When two columns are supplied, gene set names are built as "celltype__NKSEP__cytokine" internally and the heatmap rows are automatically split by the first column (celltype).

Mode 3 - MSigDB collection

gene_set_library = "H"                        # Hallmark
gene_set_library = "C2", gene_set_subcategory = "CP:KEGG"

Mode 4 - MSigDB term search (fetch all gene sets matching a pattern)

search_terms = c("CELL_DEATH", "APOPTOSIS")
# optionally restrict to a collection:
gene_set_library = "C2", search_terms = "CELL_DEATH"