From 30ce58fe8fa81c6f24b4a960f898367ec3d07791 Mon Sep 17 00:00:00 2001 From: A Wokaty Date: Tue, 28 Apr 2026 08:53:13 -0400 Subject: [PATCH 1/7] bump x.y.z version to even y prior to creation of RELEASE_3_23 branch --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6ca152c..fef0d89 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: escape Title: Easy single cell analysis platform for enrichment -Version: 2.7.3 +Version: 2.8.0 Authors@R: c( person(given = "Nick", family = "Borcherding", role = c("aut", "cre"), email = "ncborch@gmail.com"), person(given = "Jared", family = "Andrews", role = c("aut"), email = "jared.andrews07@gmail.com"), From 15bf933fefc6c815eab2bca66f457edc1541a3e6 Mon Sep 17 00:00:00 2001 From: A Wokaty Date: Tue, 28 Apr 2026 08:53:13 -0400 Subject: [PATCH 2/7] bump x.y.z version to odd y following creation of RELEASE_3_23 branch --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index fef0d89..10c590e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: escape Title: Easy single cell analysis platform for enrichment -Version: 2.8.0 +Version: 2.9.0 Authors@R: c( person(given = "Nick", family = "Borcherding", role = c("aut", "cre"), email = "ncborch@gmail.com"), person(given = "Jared", family = "Andrews", role = c("aut"), email = "jared.andrews07@gmail.com"), From 01c44a9048047dcee87c019e3291489d14a1ce3c Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 04:48:49 -0500 Subject: [PATCH 3/7] Spatal Support Per issue #180 --- R/utils.R | 175 ++++++++++++++++++++++++++------ man/performNormalization.Rd | 24 ++++- tests/testthat/helper-spatial.R | 52 ++++++++++ tests/testthat/test-spatial.R | 134 ++++++++++++++++++++++++ 4 files changed, 348 insertions(+), 37 deletions(-) create mode 100644 tests/testthat/helper-spatial.R create mode 100644 tests/testthat/test-spatial.R diff --git a/R/utils.R b/R/utils.R index 3a7636a..63e0316 100644 --- a/R/utils.R +++ b/R/utils.R @@ -10,6 +10,20 @@ stop("Expecting a Seurat or SummarizedExperiment object") } +# `base::%||%` only exists in R >= 4.4.0 and escape declares R (>= 4.1). +`%||%` <- function(x, y) if (is.null(x)) y else x + +# Single error message for every "you asked for an assay that isn't there" case, +# so a missing name never reaches assay() as NULL. +.stop_missing_assay <- function(requested, available, + what = "assay", where = "the object") { + stop("Could not find ", what, " '", requested, "' in ", where, ". ", + if (length(available)) + paste0("Available: ", paste(available, collapse = ", "), ".") + else "No assays are present.", + call. = FALSE) +} + # ----------------------------------------------------------------------------- # ORDERING UTILITY # ----------------------------------------------------------------------------- @@ -200,6 +214,68 @@ split(vec, ceiling(seq_along(vec) / chunk.size)) } +# ----------------------------------------------------------------------------- +# INPUT ASSAY RESOLUTION +# ----------------------------------------------------------------------------- +# Translate the user-facing `input.assay` into the layer (Seurat) or main assay +# (SummarizedExperiment) name that .cntEval() should read. +# +# "auto" counts for the native backend, logcounts for plaid +# "counts" Seurat layer "counts" / SCE assay "counts" +# "logcounts" Seurat layer "data" / SCE assay "logcounts" +# taken literally +.resolve_input_assay <- function(obj, input.assay = "auto", + backend = "native") { + if (!is.character(input.assay) || length(input.assay) != 1L) + stop("`input.assay` must be a single string.", call. = FALSE) + + ia <- if (identical(input.assay, "auto")) { + if (identical(backend, "plaid")) "logcounts" else "counts" + } else { + input.assay + } + + type <- if (.is_seurat(obj)) { + switch(ia, counts = "counts", logcounts = "data", data = "data", ia) + } else { + ia + } + + .assert_input_assay(obj, type, ia) + type +} + +.assert_input_assay <- function(obj, type, requested) { + if (.is_seurat(obj)) { + if (!requireNamespace("SeuratObject", quietly = TRUE)) return(invisible(TRUE)) + lay <- tryCatch(SeuratObject::Layers(obj, assay = "RNA"), + error = function(e) character()) + # Seurat v5 split objects expose "data.1", "data.2", ... rather than "data" + ok <- length(lay) == 0L || + any(lay == type | startsWith(lay, paste0(type, "."))) + if (!ok) + stop("input.assay = \"", requested, "\" needs layer '", type, + "' in the 'RNA' assay, which has: ", paste(lay, collapse = ", "), + ".\n Run Seurat::NormalizeData() (and Seurat::JoinLayers() if the ", + "object is split), or set input.assay = \"counts\".", call. = FALSE) + + } else if (.is_sce(obj)) { + if (!requireNamespace("SummarizedExperiment", quietly = TRUE)) + return(invisible(TRUE)) + an <- SummarizedExperiment::assayNames(obj) + if (!type %in% an) + stop("input.assay = \"", requested, "\" needs assay '", type, + "', which is not present. Available: ", paste(an, collapse = ", "), + ".\n Run scuttle::logNormCounts(), or set input.assay = \"counts\".", + call. = FALSE) + + } else if (identical(type, "logcounts")) { + message("`input.data` is a bare matrix; escape assumes it is already ", + "log-normalized (input.assay = \"logcounts\").") + } + invisible(TRUE) +} + # ----------------------------------------------------------------------------- # EXPRESSION MATRIX EXTRACTOR # ----------------------------------------------------------------------------- @@ -223,12 +299,21 @@ } else if (.is_sce(obj)) { if (requireNamespace("SummarizedExperiment", quietly = TRUE) && requireNamespace("SingleCellExperiment", quietly = TRUE)) { - pos <- if (assay == "RNA") "counts" else assay - cnts <- if (assay == "RNA") { - SummarizedExperiment::assay(obj, pos) + ## `assay = "RNA"` means "the main expression assay"; `type` names it. + ## Every in-package caller passes type = "counts", so this is inert + ## unless a caller deliberately asks for logcounts. + avail <- SummarizedExperiment::assayNames(obj) + if (!type %in% avail) + .stop_missing_assay(type, avail, "assay", + "the SummarizedExperiment") + SummarizedExperiment::assay(obj, type) } else { - SummarizedExperiment::assay(SingleCellExperiment::altExp(obj, pos)) + avail <- SingleCellExperiment::altExpNames(obj) + if (!assay %in% avail) + .stop_missing_assay(assay, avail, "altExp", + "the SingleCellExperiment") + SummarizedExperiment::assay(SingleCellExperiment::altExp(obj, assay)) } } else { stop("SummarizedExperiment and SingleCellExperiment packages are required but not installed.") @@ -281,8 +366,15 @@ } .pull.Enrich <- function(sc, name) { + if (is.null(name) || !is.character(name) || length(name) != 1L) + stop("`assay` must be a single enrichment assay name.", call. = FALSE) + if (.is_seurat(sc)) { if (requireNamespace("Matrix", quietly = TRUE)) { + avail <- SeuratObject::Assays(sc) + if (!name %in% avail) + .stop_missing_assay(name, avail, "enrichment assay", + "the Seurat object") so_version <- utils::packageVersion("SeuratObject") if (so_version >= "5.0.0") { Matrix::t(SeuratObject::GetAssayData(sc, assay = name, layer = "data")) @@ -292,15 +384,22 @@ } else { stop("Matrix package is required to transpose Seurat assay data.") } - + } else if (.is_sce(sc)) { if (requireNamespace("SummarizedExperiment", quietly = TRUE) && requireNamespace("SingleCellExperiment", quietly = TRUE)) { - Matrix::t(SummarizedExperiment::assay(SingleCellExperiment::altExp(sc)[[name]])) + avail <- SingleCellExperiment::altExpNames(sc) + if (!name %in% avail) + .stop_missing_assay(name, avail, "enrichment altExp", + "the SingleCellExperiment") + # altExp(sc, name), not altExp(sc)[[name]] - the latter indexes colData of + # the *first* altExp and silently returns NULL. + Matrix::t(SummarizedExperiment::assay( + SingleCellExperiment::altExp(sc, name))) } else { stop("SummarizedExperiment and SingleCellExperiment packages are required to pull enrichment from SCE object.") } - + } else { stop("Unsupported object type for pulling enrichment.") } @@ -317,6 +416,43 @@ gene.sets } +# Align a gene-set list to the columns of an enrichment matrix. +# +# Seurat coerces underscores in feature names to hyphens when an assay is built, +# so scores pulled back off a Seurat object have mangled column names while +# scores held as a plain matrix (or in an SCE altExp) keep the originals. Match +# the literal names first and fall back to the mangled form only for the columns +# that are still unmatched - mangling unconditionally silently drops every +# HALLMARK_/GO_/REACTOME_ set for non-Seurat input. +# +# Returns the sets re-keyed and re-ordered to `cols`, one per column. +.match_sets_to_cols <- function(gene.sets, cols) { + nm <- names(gene.sets) + if (is.null(nm)) + stop("`gene.sets` must be a named list.", call. = FALSE) + if (is.null(cols)) + stop("Enrichment matrix has no column names; cannot match gene sets.", + call. = FALSE) + + idx <- match(cols, nm) + miss <- is.na(idx) + if (any(miss)) + idx[miss] <- match(cols[miss], gsub("_", "-", nm, fixed = TRUE)) + + if (all(is.na(idx))) + stop("None of the supplied gene sets match enrichment columns.", + call. = FALSE) + if (anyNA(idx)) + stop("No gene set supplied for enrichment column(s): ", + paste(cols[is.na(idx)], collapse = ", "), + ". Supply the same `gene.sets` used to compute the scores.", + call. = FALSE) + + out <- gene.sets[idx] + names(out) <- cols + out +} + .grabMeta <- function(sc) { if (.is_seurat(sc)) { if (!requireNamespace("SeuratObject", quietly = TRUE)) { @@ -433,31 +569,6 @@ ) } -#─ Split a matrix into equal-sized column chunks ------------------------------ -.split_cols <- function(mat, chunk) { - if (ncol(mat) <= chunk) return(list(mat)) - idx <- split(seq_len(ncol(mat)), ceiling(seq_len(ncol(mat)) / chunk)) - lapply(idx, function(i) mat[, i, drop = FALSE]) -} - -.match_summary_fun <- function(fun) { - if (is.function(fun)) return(fun) - - if (!is.character(fun) || length(fun) != 1L) - stop("'summary.fun' must be a single character or a function") - - kw <- tolower(fun) - fn <- switch(kw, - mean = base::mean, - median = stats::median, - max = base::max, - sum = base::sum, - geometric = function(x) exp(mean(log(x + 1e-6))), - stop("Unsupported summary keyword: ", fun)) - attr(fn, "keyword") <- kw # tag for fast matrixStats branch - fn -} - .computeRunningES <- function(gene.order, hits, weight = NULL) { N <- length(gene.order) hit <- gene.order %in% hits diff --git a/man/performNormalization.Rd b/man/performNormalization.Rd index f1f08b5..d69122f 100644 --- a/man/performNormalization.Rd +++ b/man/performNormalization.Rd @@ -17,15 +17,19 @@ performNormalization( \arguments{ \item{input.data}{A raw-counts matrix (genes x cells), a \link[SeuratObject]{Seurat} object, or a -\link[SingleCellExperiment]{SingleCellExperiment}. Gene identifiers must +\link[SingleCellExperiment]{SingleCellExperiment} (including a +\link[SpatialExperiment]{SpatialExperiment}). Gene identifiers must match those in \code{gene.sets}.} \item{enrichment.data}{Matrix. Output of \code{\link{escape.matrix}} or -\code{NULL} if enrichment scores are already stored in \code{input.data}.} +\code{NULL} if enrichment scores are already stored in \code{input.data}. +When supplied it takes precedence over anything stored in +\code{input.data}.} \item{assay}{Character. Name of the assay holding enrichment scores when \code{input.data} is a single-cell object. Default is \code{"escape"}. -Ignored otherwise.} +Ignored when \code{input.data} is a matrix. Set to \code{NULL} to return +the normalized matrix rather than attaching it to the object.} \item{gene.sets}{A named list of character vectors, the result of \code{\link{getGeneSets}}, or the built-in data object @@ -44,8 +48,9 @@ Larger values reduce overhead but increase memory usage. Default is \code{NULL} (process all cells at once).} } \value{ -If `input.data` is an object, the same object with a new assay - "_normalized". Otherwise a matrix of normalized scores. +If `input.data` is an object and `assay` is not `NULL`, the same + object with a new assay "_normalized". Otherwise a matrix of + normalized scores. } \description{ Scales each enrichment value by the \strong{number of genes from the set @@ -53,6 +58,15 @@ that are expressed} in that cell (non-zero counts). Optionally shifts results into a positive range and/or applies a natural-log transform for compatibility with log-based differential tests. } +\section{Which expression values are used}{ + +The per-cell scale factor is the number of genes from each set with a +\strong{non-zero raw count}, so this function always reads the \code{counts} +assay regardless of any \code{input.assay} used when the scores were +computed. Detection is identical in count and log space, so this is +deliberate rather than an oversight. +} + \examples{ gs <- list(Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"), Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A")) diff --git a/tests/testthat/helper-spatial.R b/tests/testthat/helper-spatial.R new file mode 100644 index 0000000..1da5473 --- /dev/null +++ b/tests/testthat/helper-spatial.R @@ -0,0 +1,52 @@ +# Toy SpatialExperiment fixture. +# +# Built on the fly rather than shipped as testdata: the object is tiny, fully +# deterministic, and keeps SpatialExperiment out of the hard test dependencies. +# Everything is fully qualified so sourcing this helper never fails when +# SpatialExperiment is absent - call skip_if_not_installed("SpatialExperiment") +# inside the test_that() block instead. + +make_toy_spe <- function(n.genes = 40, n.cells = 60, n.samples = 2, + seed = 42) { + set.seed(seed) + cnts <- Matrix::rsparsematrix( + n.genes, n.cells, density = 0.4, + rand.x = function(n) stats::rpois(n, 5) + 1 + ) + dimnames(cnts) <- list(paste0("gene", seq_len(n.genes)), + paste0("cell", seq_len(n.cells))) + + # deterministic log values so tests never depend on scuttle/scater + logc <- log1p(cnts) + + SpatialExperiment::SpatialExperiment( + assays = list(counts = cnts, logcounts = logc), + colData = S4Vectors::DataFrame( + group = rep(c("a", "b"), length.out = n.cells), + row.names = colnames(cnts) + ), + spatialCoords = matrix( + seq_len(2 * n.cells), ncol = 2, + dimnames = list(colnames(cnts), c("x", "y")) + ), + sample_id = rep(paste0("s", seq_len(n.samples)), length.out = n.cells) + ) +} + +# Underscored name on purpose: HALLMARK_/GO_/REACTOME_ sets are the case that +# used to be silently dropped during normalization for non-Seurat input. +toy_spe_sets <- function() { + list(HALLMARK_SET_A = paste0("gene", 1:8), + SetB = paste0("gene", 9:16)) +} + +# Matching SingleCellExperiment, for the regression guard that the same bug was +# never spatial-specific. +make_toy_sce <- function(...) { + spe <- make_toy_spe(...) + SingleCellExperiment::SingleCellExperiment( + assays = list(counts = SummarizedExperiment::assay(spe, "counts"), + logcounts = SummarizedExperiment::assay(spe, "logcounts")), + colData = SummarizedExperiment::colData(spe)[, "group", drop = FALSE] + ) +} diff --git a/tests/testthat/test-spatial.R b/tests/testthat/test-spatial.R new file mode 100644 index 0000000..f7a9fc6 --- /dev/null +++ b/tests/testthat/test-spatial.R @@ -0,0 +1,134 @@ +# test script for SpatialExperiment / SingleCellExperiment input - +# testcases are NOT comprehensive! +# +# Regression coverage for the bug reported in #180: escape.matrix(normalize = +# TRUE) threw an S4 dispatch error ("assay" on a NULL) for every +# SummarizedExperiment-derived input, not just SpatialExperiment. + +gs <- toy_spe_sets() + +test_that("escape.matrix(normalize = TRUE) works on a SpatialExperiment", { + skip_if_not_installed("SpatialExperiment") + spe <- make_toy_spe() + + res <- escape.matrix(spe, gene.sets = gs, method = "UCell", + normalize = TRUE, min.size = NULL) + + expect_true(is.matrix(res)) + expect_equal(dim(res), c(ncol(spe), length(gs))) + expect_equal(colnames(res), names(gs)) + expect_equal(rownames(res), colnames(spe)) + expect_true(all(is.finite(res))) +}) + +test_that("escape.matrix(normalize = TRUE) works on a plain SingleCellExperiment", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + + res <- escape.matrix(sce, gene.sets = gs, method = "UCell", + normalize = TRUE, min.size = NULL) + + expect_equal(dim(res), c(ncol(sce), length(gs))) + expect_equal(colnames(res), names(gs)) + expect_true(all(is.finite(res))) +}) + +test_that("normalized scores agree across SPE, SCE and raw matrix input", { + skip_if_not_installed("SpatialExperiment") + spe <- make_toy_spe() + sce <- make_toy_sce() + mat <- as.matrix(SummarizedExperiment::assay(spe, "counts")) + + f <- function(x) escape.matrix(x, gene.sets = gs, method = "UCell", + normalize = TRUE, min.size = NULL) + + expect_equal(f(spe), f(sce), tolerance = 1e-10) + expect_equal(f(spe), f(mat), tolerance = 1e-10) +}) + +test_that("runEscape() preserves spatial metadata", { + skip_if_not_installed("SpatialExperiment") + spe <- make_toy_spe() + out <- runEscape(spe, gene.sets = gs, method = "UCell", min.size = NULL) + + expect_s4_class(out, "SpatialExperiment") + expect_true("escape" %in% SingleCellExperiment::altExpNames(out)) + expect_equal(SpatialExperiment::spatialCoords(out), + SpatialExperiment::spatialCoords(spe)) + expect_equal(out$sample_id, spe$sample_id) + expect_equal(nrow(SpatialExperiment::imgData(out)), + nrow(SpatialExperiment::imgData(spe))) + + # colnames must stay aligned between the altExp and its parent + alt <- SingleCellExperiment::altExp(out, "escape") + expect_equal(colnames(alt), colnames(spe)) + expect_equal(rownames(alt), names(gs)) +}) + +test_that("multi-sample SpatialExperiment completes with chunking", { + skip_if_not_installed("SpatialExperiment") + spe <- make_toy_spe(n.cells = 60, n.samples = 3) + expect_gt(length(unique(spe$sample_id)), 1L) + + res <- escape.matrix(spe, gene.sets = gs, method = "UCell", + normalize = TRUE, min.size = NULL, groups = 25) + expect_equal(dim(res), c(ncol(spe), length(gs))) +}) + +test_that("performNormalization() round-trips through an SCE altExp", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + obj <- runEscape(sce, gene.sets = gs, method = "UCell", min.size = NULL) + + out <- performNormalization(obj, assay = "escape", gene.sets = gs) + + expect_true("escape_normalized" %in% SingleCellExperiment::altExpNames(out)) + expect_equal( + dim(SingleCellExperiment::altExp(out, "escape_normalized")), + c(length(gs), ncol(sce)) + ) +}) + +test_that("performPCA() works on a SingleCellExperiment", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + obj <- runEscape(sce, gene.sets = gs, method = "UCell", min.size = NULL) + + out <- performPCA(obj, assay = "escape") + expect_s4_class(out, "SingleCellExperiment") +}) + +test_that("input.assay selects the expression matrix and errors informatively", { + skip_if_not_installed("SpatialExperiment") + spe <- make_toy_spe() + + res <- escape.matrix(spe, gene.sets = gs, method = "UCell", + min.size = NULL, input.assay = "logcounts") + expect_equal(dim(res), c(ncol(spe), length(gs))) + expect_equal(colnames(res), names(gs)) + + expect_error( + escape.matrix(spe, gene.sets = gs, method = "UCell", + min.size = NULL, input.assay = "nope"), + "needs assay 'nope'" + ) +}) + +test_that(".resolve_input_assay maps names per object class", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + + expect_equal(escape:::.resolve_input_assay(sce, "auto", "native"), "counts") + expect_equal(escape:::.resolve_input_assay(sce, "auto", "plaid"), "logcounts") + expect_equal(escape:::.resolve_input_assay(SeuratObject::pbmc_small, + "logcounts", "native"), "data") + expect_equal(escape:::.resolve_input_assay(SeuratObject::pbmc_small, + "auto", "native"), "counts") + + # an SCE without logcounts must say how to get them + bare <- SingleCellExperiment::SingleCellExperiment( + assays = list(counts = SummarizedExperiment::assay(sce, "counts")) + ) + expect_error(escape:::.resolve_input_assay(bare, "auto", "plaid"), + "logNormCounts") +}) From 5a65356091a79bb36d2853c659eef74ff0aa0a27 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 04:49:14 -0500 Subject: [PATCH 4/7] Add plaid per issue #167 --- R/plaid.R | 212 ++++++++++++++++ R/runEscape.R | 225 ++++++++++++----- tests/testthat/test-performNormalization.R | 94 +++++++ tests/testthat/test-plaid.R | 278 +++++++++++++++++++++ tests/testthat/test-runEscape.R | 33 ++- 5 files changed, 770 insertions(+), 72 deletions(-) create mode 100644 R/plaid.R create mode 100644 tests/testthat/test-plaid.R diff --git a/R/plaid.R b/R/plaid.R new file mode 100644 index 0000000..957b2e3 --- /dev/null +++ b/R/plaid.R @@ -0,0 +1,212 @@ +# ----------------------------------------------------------------------------- +# PLAID BACKEND +# +# plaid (Zito et al., Bioinformatics 2025, btaf621) reimplements several +# single-sample enrichment scores on top of a single sparse crossprod. escape +# exposes it as an opt-in backend rather than a replacement: the replaid.* +# functions are fast *approximations* of the methods they are named after, not +# drop-in numerical equivalents. +# ----------------------------------------------------------------------------- + +# canonical key -> label used in messages and errors +.METHOD_LABELS <- c(SSGSEA = "ssGSEA", + GSVA = "GSVA", + UCELL = "UCell", + AUCELL = "AUCell", + PLAID = "PLAID", + SINGSCORE = "singscore", + SCSE = "scSE") + +# canonical key -> exported plaid function name +.PLAID_FUNS <- c(SSGSEA = "replaid.ssgsea", + GSVA = "replaid.gsva", + UCELL = "replaid.ucell", + AUCELL = "replaid.aucell", + SINGSCORE = "replaid.sing", + SCSE = "replaid.scse", + PLAID = "plaid") + +# methods that exist only through plaid +.PLAID_ONLY <- c("PLAID", "SINGSCORE", "SCSE") + +# methods escape can compute itself +.NATIVE_METHODS <- c("SSGSEA", "GSVA", "UCELL", "AUCELL") + +# one-time-per-session message bookkeeping +.escape_env <- new.env(parent = emptyenv()) + +# ----------------------------------------------------------------------------- +# METHOD / BACKEND RESOLUTION +# ----------------------------------------------------------------------------- +.resolve_backend <- function(method, backend = c("native", "plaid")) { + backend <- match.arg(backend) + + if (!is.character(method) || length(method) != 1L) + stop("`method` must be a single string.", call. = FALSE) + + key <- toupper(method) + if (!key %in% names(.METHOD_LABELS)) + stop("Unknown `method`: '", method, "'. One of: ", + paste(.METHOD_LABELS, collapse = ", "), ".", call. = FALSE) + + ## PLAID / singscore / scSE have no native implementation, so `backend` is + ## not a meaningful choice for them - route to plaid and say so. + if (key %in% .PLAID_ONLY && backend != "plaid") { + message("method = '", .METHOD_LABELS[[key]], "' is provided by the plaid ", + "backend; using backend = \"plaid\".") + backend <- "plaid" + } + + list(method = key, backend = backend) +} + +# ----------------------------------------------------------------------------- +# DEPENDENCY GUARD +# ----------------------------------------------------------------------------- +.require_plaid <- function() { + if (!requireNamespace("plaid", quietly = TRUE)) + stop("plaid not installed. Install it with ", + "BiocManager::install(\"plaid\") to use backend = \"plaid\".", + call. = FALSE) + invisible(TRUE) +} + +.plaid_fun <- function(key) { + .require_plaid() + utils::getFromNamespace(.PLAID_FUNS[[key]], "plaid") +} + +# Warn once per session that plaid scores are approximations. +.plaid_fidelity_note <- function(key) { + if (isTRUE(.escape_env$plaid_noted)) return(invisible(NULL)) + .escape_env$plaid_noted <- TRUE + message("Scoring with plaid::", .PLAID_FUNS[[key]], "(). The replaid.* ", + "functions are fast reimplementations and are not guaranteed to ", + "reproduce the native escape scores exactly - see ", + "?escape.matrix for the per-method fidelity notes.") + invisible(NULL) +} + +# ----------------------------------------------------------------------------- +# DISPATCH +# ----------------------------------------------------------------------------- +# expr : genes x cells (sparse is kept sparse) +# gene_sets : named list of character vectors, already min.size-filtered +# fn : injection point - supply a stub to unit test the dispatch logic +# without plaid installed +# returns : cells x gene-sets +.plaid_dispatch <- function(expr, gene_sets, method, + min.size = 5, + backend.args = list(), + groups = NULL, + fn = NULL) { + key <- toupper(method) + if (!key %in% names(.PLAID_FUNS)) + stop("Unknown `method`: '", method, "'.", call. = FALSE) + + if (is.null(fn)) { + fn <- .plaid_fun(key) + .plaid_fidelity_note(key) + } + fmls <- names(formals(fn)) + + if (!length(gene_sets)) + stop("No gene sets to score.", call. = FALSE) + if (is.null(names(gene_sets))) + stop("`gene.sets` must be a named list.", call. = FALSE) + + ## ---- defaults escape owns ------------------------------------------------- + ## plaid defaults to max.genes = 500, which silently drops most HALLMARK, GO + ## and REACTOME sets. Set a real cap above anything scoreable rather than + ## relying on a sentinel whose meaning lives inside plaid. + cap <- max(c(lengths(gene_sets), nrow(expr), 1L)) + args <- list() + if ("min.genes" %in% fmls) + args$min.genes <- max(1L, as.integer(min.size %||% 1L)) + if ("max.genes" %in% fmls) + args$max.genes <- as.integer(cap) + ## `chunk` is documented only on plaid(); do not guess it through replaid dots + if (key == "PLAID" && "chunk" %in% fmls && !is.null(groups)) + args$chunk <- as.integer(groups) + ## plaid's `assay=` is deliberately never set - escape always hands plaid a + ## bare matrix, so plaid never has to reach into an object. + + ## ---- user overrides ------------------------------------------------------- + if (length(backend.args)) { + if (is.null(names(backend.args)) || any(!nzchar(names(backend.args)))) + stop("`backend.args` must be a fully named list.", call. = FALSE) + bad <- setdiff(names(backend.args), fmls) + if (length(bad)) + stop("Unknown `backend.args` for method '", .METHOD_LABELS[[key]], + "' (plaid::", .PLAID_FUNS[[key]], "): ", + paste(bad, collapse = ", "), ".\n Accepted: ", + paste(setdiff(fmls, c("X", "matG", "...")), collapse = ", "), + call. = FALSE) + if ("normalize" %in% names(backend.args)) + message("`backend.args$normalize` is plaid's median normalization of ", + "the scores. escape's own `normalize =` argument is post-hoc ", + "drop-out scaling (performNormalization()). They are independent.") + args <- utils::modifyList(args, backend.args, keep.null = TRUE) + } + + ## ---- call ----------------------------------------------------------------- + out <- do.call(fn, c(list(X = expr, matG = gene_sets), args)) + out <- as.matrix(out) + + ## ---- validate + orient (plaid returns gene sets x cells) ------------------ + ## an empty result loses its (zero-length) rownames, so check it before the + ## dimnames guard or the user gets a misleading message + if (!nrow(out)) + stop("plaid scored none of the ", length(gene_sets), " gene set(s). ", + "Check that identifiers in `gene.sets` match rownames of the ", + "expression matrix, and inspect backend.args$min.genes / ", + "backend.args$max.genes.", call. = FALSE) + if (is.null(rownames(out)) || is.null(colnames(out))) + stop("plaid returned a matrix without dimnames; cannot align results.", + call. = FALSE) + if (!setequal(colnames(out), colnames(expr))) + stop("plaid returned ", ncol(out), " cell(s) but ", ncol(expr), + " were supplied.", call. = FALSE) + out <- out[, colnames(expr), drop = FALSE] + + dropped <- setdiff(names(gene_sets), rownames(out)) + if (length(dropped) == length(gene_sets)) + stop("plaid scored none of the ", length(gene_sets), " gene set(s). ", + "Check that identifiers in `gene.sets` match rownames of the ", + "expression matrix.", call. = FALSE) + if (length(dropped)) + warning("plaid dropped ", length(dropped), " of ", length(gene_sets), + " gene set(s): ", paste(utils::head(dropped, 5L), collapse = ", "), + if (length(dropped) > 5L) ", ..." else "", + ". Inspect backend.args$min.genes / backend.args$max.genes.", + call. = FALSE) + + out <- out[intersect(names(gene_sets), rownames(out)), , drop = FALSE] + t(out) # cells x gene-sets +} + +# Record which engine produced a score matrix. .adding.Enrich() drops +# attributes, so runEscape() also stashes this on the object itself. +.stamp_backend <- function(mat, method, backend) { + attr(mat, "escape.backend") <- list( + backend = backend, + method = unname(.METHOD_LABELS[[toupper(method)]]), + plaid.version = if (identical(backend, "plaid") && + requireNamespace("plaid", quietly = TRUE)) + as.character(utils::packageVersion("plaid")) else NA_character_ + ) + mat +} + +# Persist provenance on the object, since assay containers drop attributes. +.record_backend <- function(sc, name, prov) { + if (is.null(prov)) return(sc) + if (.is_seurat(sc)) { + if (requireNamespace("SeuratObject", quietly = TRUE)) + SeuratObject::Misc(sc, slot = paste0(name, "_backend")) <- prov + } else if (.is_sce(sc)) { + if (requireNamespace("S4Vectors", quietly = TRUE)) + S4Vectors::metadata(sc)[[paste0(name, "_backend")]] <- prov + } + sc +} diff --git a/R/runEscape.R b/R/runEscape.R index 334ffe8..4e144b0 100644 --- a/R/runEscape.R +++ b/R/runEscape.R @@ -14,21 +14,69 @@ #' \item{\code{"ssGSEA"}}{Single-sample GSEA.} #' \item{\code{"UCell"}}{Rank-based UCell scoring.} #' \item{\code{"AUCell"}}{Area-under-the-curve ranking score.} +#' \item{\code{"PLAID"}}{Average log-intensity of set members (plaid only).} +#' \item{\code{"singscore"}}{Rank-based singscore (plaid only).} +#' \item{\code{"scSE"}}{Single-cell signature explorer score (plaid only).} #' } #' +#' @section Backends: +#' The first four methods run on escape's own engines by default +#' (\code{backend = "native"}). Setting \code{backend = "plaid"} routes them to +#' \pkg{plaid}'s \code{replaid.*} family instead, which is substantially faster +#' and lighter on memory for large objects. \code{"PLAID"}, \code{"singscore"} +#' and \code{"scSE"} have no native implementation and always use \pkg{plaid}. +#' +#' \strong{The plaid backend approximates rather than reproduces, and the gap is +#' wider than the plaid documentation suggests.} On a simulated 2000-gene, +#' 120-cell matrix, pooled Pearson correlation between the native and plaid +#' scores for the same method was roughly 0.85 (\code{ssGSEA}), 0.83 +#' (\code{UCell}, \code{AUCell}) and 0.73 (\code{GSVA}); on the 230-gene +#' \code{pbmc_small} it was lower still. The scores are also on different +#' scales, not just noisier. Notably \code{replaid.ssgsea} is documented as +#' exact at \code{alpha = 0}, but compared directly against +#' \code{GSVA::gsva()} on identical input it correlated at 0.85, not 1. The +#' input assay is not the cause - these are rank-based scores, and counts +#' versus logcounts correlate at exactly 1. +#' +#' Treat \code{backend = "plaid"} as a fast screen, not as a drop-in +#' replacement. Do not mix backends within an analysis, and do not compare +#' plaid scores against previously published \pkg{escape} results. Tuning that +#' may narrow the gap for individual methods: +#' \describe{ +#' \item{\code{ssGSEA}}{\code{backend.args = list(alpha = 0)}.} +#' \item{\code{GSVA}}{the empirical CDF row transform is approximated by a +#' z-transform (\code{rowtf = "z"}); pass +#' \code{backend.args = list(rowtf = "ecdf")} for the slower exact form.} +#' \item{\code{UCell}}{\code{backend.args = list(rmax = ...)} shifts the score +#' scale but did not change rank agreement in testing.} +#' \item{\code{scSE}}{plaid documents a match to the original with +#' \code{backend.args = list(removeLog2 = TRUE, scoreMean = FALSE)}.} +#' } +#' +#' Users of the plaid backend should cite Zito \emph{et al.}, \emph{Bioinformatics} +#' 2025, 41(12):btaf621 in addition to the original method paper and +#' \pkg{escape}. +#' #' @param input.data A raw-counts matrix (genes x cells), a #' \link[SeuratObject]{Seurat} object, or a -#' \link[SingleCellExperiment]{SingleCellExperiment}. Gene identifiers must +#' \link[SingleCellExperiment]{SingleCellExperiment} (including a +#' \link[SpatialExperiment]{SpatialExperiment}). Gene identifiers must #' match those in \code{gene.sets}. #' @param gene.sets A named list of character vectors, the result of #' \code{\link{getGeneSets}}, or the built-in data object #' \code{\link{escape.gene.sets}}. List names become column names in the #' result. #' @param method Character. Scoring algorithm (case-insensitive). One of -#' \code{"GSVA"}, \code{"ssGSEA"}, \code{"UCell"}, or \code{"AUCell"}. -#' Default is \code{"ssGSEA"}. +#' \code{"GSVA"}, \code{"ssGSEA"}, \code{"UCell"}, \code{"AUCell"}, +#' \code{"PLAID"}, \code{"singscore"}, or \code{"scSE"}. The last three are +#' available only through \code{backend = "plaid"} and select it +#' automatically. Default is \code{"ssGSEA"}. #' @param groups Integer. Number of cells per processing chunk. Larger values #' reduce overhead but increase memory usage. Default is \code{1000}. +#' Meaning depends on the backend: chunk size for the \pkg{BiocParallel} loop +#' when \code{backend = "native"}, forwarded to \code{plaid::plaid(chunk=)} +#' for \code{method = "PLAID"}, and ignored for the \code{replaid.*} paths +#' (use \code{backend.args$chunk} there). #' @param min.size Integer or \code{NULL}. Minimum number of genes from a set #' that must be detected in the expression matrix for that set to be scored. #' Default is \code{5}. Use \code{NULL} to disable filtering. @@ -46,9 +94,27 @@ #' which the \code{min.expr.cells} rule is applied. Default is \code{NULL}. #' @param BPPARAM A \pkg{BiocParallel} parameter object describing the #' parallel backend. Default is \code{NULL} (serial execution). -#' @param ... Extra arguments passed verbatim to the chosen back-end scoring +#' @param ... Extra arguments passed verbatim to the chosen native scoring #' function (\code{gsva()}, \code{ScoreSignatures_UCell()}, or -#' \code{AUCell_calcAUC()}). +#' \code{AUCell_calcAUC()}). Not forwarded when \code{backend = "plaid"} - +#' use \code{backend.args} there. +#' @param backend Character. Scoring engine, \code{"native"} (default) or +#' \code{"plaid"}. Ignored for methods that only exist in \pkg{plaid}. +#' @param backend.args Named list of method-specific tuning arguments passed to +#' the underlying \pkg{plaid} function, e.g. \code{alpha}, \code{tau}, +#' \code{rowtf}, \code{aucMaxRank}, \code{rmax}, \code{nsmooth}, +#' \code{stats}, \code{chunk}, \code{removeLog2}, \code{scoreMean}. Names are +#' validated against the target function. Note that +#' \code{backend.args$normalize} is \pkg{plaid}'s median normalization of the +#' scores and is unrelated to escape's \code{normalize} argument. Default is +#' \code{list()}. +#' @param input.assay Character. Which expression matrix to score. +#' \code{"auto"} (default) reads raw counts for the native backend and +#' log-normalized values for \pkg{plaid}. \code{"counts"} and +#' \code{"logcounts"} map to the right layer for both \pkg{Seurat} +#' (\code{counts} / \code{data}) and \pkg{SummarizedExperiment}-derived +#' objects, including \link[SpatialExperiment]{SpatialExperiment}. Any other +#' string is taken literally. #' #' @return A numeric matrix with one row per cell and one column per gene set, #' ordered as in \code{gene.sets}. @@ -81,13 +147,23 @@ escape.matrix <- function(input.data, min.expr.cells = 0, min.filter.by = NULL, BPPARAM = NULL, - ...) { + ..., + backend = c("native", "plaid"), + backend.args = list(), + input.assay = "auto") { if(is.null(min.size)) min.size <- 0 - + + # ---- 0) resolve method / backend ------------------------------------------ + res <- .resolve_backend(method, backend) + method <- res$method + backend <- res$backend + # ---- 1) resolve gene-sets & counts ---------------------------------------- egc <- .GS.check(gene.sets) - cnts <- .cntEval(input.data, assay = "RNA", type = "counts") # dgCMatrix - + cnts <- .cntEval(input.data, assay = "RNA", + type = .resolve_input_assay(input.data, input.assay, + backend)) # dgCMatrix + if (is.null(min.filter.by)) { cnts <- .filter_genes(cnts, min.expr.cells) } else { @@ -110,44 +186,69 @@ escape.matrix <- function(input.data, stop("No gene-sets meet the size threshold (min.size = ", min.size, ")") } - # ---- 3) split cells into chunks ------------------------------------------- - chunks <- .split_cols(cnts, groups) - message("escape.matrix(): processing ", length(chunks), " chunk(s)...") - - # ---- 4) compute enrichment in parallel ------------------------------------ - res_list <- .plapply( - chunks, - function(mat) - .compute_enrichment(mat, egc, method, BPPARAM, ...), - BPPARAM = BPPARAM - ) - - # ---- 5) combine + orient (rows = cells) ----------------------------------- - all_sets <- names(egc) - res_mat <- do.call(cbind, lapply(res_list, function(m) { - m <- as.matrix(m) - m <- m[match(all_sets, rownames(m)), , drop = FALSE] - m - })) - res_mat <- t(res_mat) - colnames(res_mat) <- all_sets - + # ---- 3-5) score ----------------------------------------------------------- + if (backend == "plaid") { + ## plaid does its own chunked crossprod and its own forking. Running + ## escape's chunk loop on top would (a) make rank- and median-based scores + ## depend on `groups`, (b) fork inside a fork, and (c) throw away the sparse + ## speed-up that is the whole point of the backend. + if (...length()) + stop("Extra arguments in `...` are not forwarded when ", + "backend = \"plaid\". Pass method-specific tuning through ", + "`backend.args = list(...)`.", call. = FALSE) + if (!is.null(BPPARAM) && !inherits(BPPARAM, "SerialParam")) + message("backend = \"plaid\" manages its own parallelism; ", + "`BPPARAM` is ignored.") + if (normalize && !identical(backend.args$normalize, FALSE)) + message("normalize = TRUE stacks escape's drop-out scaling on top of ", + "plaid's median normalization. Set ", + "backend.args = list(normalize = FALSE) to use escape's only.") + + .require_plaid() + message("escape.matrix(): scoring ", ncol(cnts), + " cells with the plaid backend...") + res_mat <- .plaid_dispatch(cnts, egc, method, + min.size = min.size, + backend.args = backend.args, + groups = groups) + } else { + chunks <- .split_cols(cnts, groups) + message("escape.matrix(): processing ", length(chunks), " chunk(s)...") + + res_list <- .plapply( + chunks, + function(mat) + .compute_enrichment(mat, egc, method, BPPARAM, ...), + BPPARAM = BPPARAM + ) + + ## combine + orient (rows = cells) + all_sets <- names(egc) + res_mat <- do.call(cbind, lapply(res_list, function(m) { + m <- as.matrix(m) + m <- m[match(all_sets, rownames(m)), , drop = FALSE] + m + })) + res_mat <- t(res_mat) + colnames(res_mat) <- all_sets + } + # ---- 6) optional dropout scaling ------------------------------------------ if (normalize) { + ## assay = NULL keeps this on the matrix path for every input class - the + ## previous round trip through .adding.Enrich()/.pull.Enrich() is what broke + ## SingleCellExperiment and SpatialExperiment input. res_mat <- performNormalization( input.data = input.data, enrichment.data = res_mat, assay = NULL, - gene.sets = gene.sets, + gene.sets = egc, make.positive = make.positive, groups = groups ) - if (.is_seurat_or_sce(input.data)) { - res_mat <- .pull.Enrich(res_mat, "escape_normalized") - } } - - res_mat + + .stamp_backend(res_mat, method, backend) } #' Calculate Enrichment Scores Using Seurat or SingleCellExperiment Objects @@ -188,7 +289,8 @@ escape.matrix <- function(input.data, #' @export runEscape <- function(input.data, gene.sets, - method = c("ssGSEA", "GSVA", "UCell", "AUCell"), + method = c("ssGSEA", "GSVA", "UCell", "AUCell", + "PLAID", "singscore", "scSE"), groups = 1e3, min.size = 5, normalize = FALSE, @@ -197,14 +299,34 @@ runEscape <- function(input.data, min.expr.cells = 0, min.filter.by = NULL, BPPARAM = NULL, - ...) { - method <- match.arg(method) + ..., + backend = c("native", "plaid"), + backend.args = list(), + input.assay = "auto") { + method <- match.arg(method) + backend <- match.arg(backend) .checkSingleObject(input.data) - esc <- escape.matrix(input.data, gene.sets, method, groups, min.size, - normalize, make.positive, min.expr.cells, - min.filter.by, BPPARAM, ...) - + + ## named, not positional - inserting an argument above must never silently + ## shift what the callee receives + esc <- escape.matrix(input.data = input.data, + gene.sets = gene.sets, + method = method, + groups = groups, + min.size = min.size, + normalize = normalize, + make.positive = make.positive, + min.expr.cells = min.expr.cells, + min.filter.by = min.filter.by, + BPPARAM = BPPARAM, + ..., + backend = backend, + backend.args = backend.args, + input.assay = input.assay) + + prov <- attr(esc, "escape.backend") input.data <- .adding.Enrich(input.data, esc, new.assay.name) + input.data <- .record_backend(input.data, new.assay.name, prov) return(input.data) } @@ -233,18 +355,3 @@ runEscape <- function(input.data, return(colData(obj)[[col]]) stop("min.filter.by requires a Seurat or SingleCellExperiment object") } - -.filter_genes <- function(m, min.expr.cells) { - if (is.null(min.expr.cells) || identical(min.expr.cells, 0)) - return(m) # nothing to do - - ncells <- ncol(m) - - thr <- if (min.expr.cells < 1) - ceiling(min.expr.cells * ncells) # proportion → absolute - else - as.integer(min.expr.cells) - - keep <- Matrix::rowSums(m > 0) >= thr - m[keep, , drop = FALSE] -} diff --git a/tests/testthat/test-performNormalization.R b/tests/testthat/test-performNormalization.R index 7bbc1f8..db308ce 100644 --- a/tests/testthat/test-performNormalization.R +++ b/tests/testthat/test-performNormalization.R @@ -107,6 +107,100 @@ test_that("error handling works", { ), "None of the supplied gene sets match" ) + + # a gene set for only some of the columns must not silently misalign + expect_error( + performNormalization( + input.data = toy_counts, + enrichment.data = toy_enrich, + gene.sets = list(Set1 = c("g1", "g2")) + ), + "No gene set supplied for enrichment column" + ) +}) + +# -------------------------------------------------------------------------- +# Underscored set names (HALLMARK_*, GO_*, REACTOME_*) used to be mangled to +# hyphens unconditionally, which dropped every such set for non-Seurat input. +under_sets <- list(HALLMARK_SET_ONE = c("g1", "g2"), + Set2 = c("g2", "g3")) +under_enrich <- toy_enrich +colnames(under_enrich) <- names(under_sets) + +test_that("underscored gene-set names normalize on matrix input", { + norm <- performNormalization( + input.data = toy_counts, + enrichment.data = under_enrich, + gene.sets = under_sets + ) + expect_equal(dim(norm), dim(under_enrich)) + expect_equal(colnames(norm), names(under_sets)) + expect_true(all(is.finite(norm))) + + gs_counts_c1 <- c( + sum(toy_counts[c("g1", "g2"), "c1"] != 0), + sum(toy_counts[c("g2", "g3"), "c1"] != 0) + ) + manual <- log1p(under_enrich["c1", ] / gs_counts_c1 + 1e-6) + expect_equal(unname(norm["c1", ]), unname(manual)) +}) + +test_that("underscored names give identical scores across input classes", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + gs <- toy_spe_sets() # HALLMARK_SET_A + SetB + mat <- as.matrix(SummarizedExperiment::assay(sce, "counts")) + + f <- function(x) escape.matrix(x, gene.sets = gs, method = "UCell", + normalize = TRUE, min.size = NULL) + + from_sce <- f(sce) + from_mat <- f(mat) + expect_equal(from_sce, from_mat, tolerance = 1e-10) + expect_equal(colnames(from_sce), names(gs)) +}) + +# -------------------------------------------------------------------------- +test_that("supplied enrichment.data wins over scores held on the object", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + gs <- toy_spe_sets() + obj <- runEscape(sce, gene.sets = gs, method = "UCell", min.size = NULL) + + # a matrix that is deliberately nothing like the stored scores + fake <- matrix(1, nrow = ncol(sce), ncol = length(gs), + dimnames = list(colnames(sce), names(gs))) + + expect_warning( + out <- performNormalization(obj, enrichment.data = fake, + assay = "escape", gene.sets = gs), + "using `enrichment.data`" + ) + + from_fake <- performNormalization(SummarizedExperiment::assay(sce, "counts"), + enrichment.data = fake, gene.sets = gs) + expect_equal( + Matrix::t(SummarizedExperiment::assay( + SingleCellExperiment::altExp(out, "escape_normalized"))), + from_fake, + tolerance = 1e-10, ignore_attr = TRUE + ) +}) + +test_that("a nonexistent enrichment assay names what is available", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + gs <- toy_spe_sets() + obj <- runEscape(sce, gene.sets = gs, method = "UCell", min.size = NULL) + + expect_error( + performNormalization(obj, assay = "not_there", gene.sets = gs), + "Could not find enrichment assay 'not_there'" + ) + expect_error( + performNormalization(obj, assay = "not_there", gene.sets = gs), + "Available: escape" + ) }) diff --git a/tests/testthat/test-plaid.R b/tests/testthat/test-plaid.R new file mode 100644 index 0000000..332c579 --- /dev/null +++ b/tests/testthat/test-plaid.R @@ -0,0 +1,278 @@ +# test script for plaid.R - testcases are NOT comprehensive! +# +# The suite runs at testthat edition 2 (no Config/testthat/edition in +# DESCRIPTION), so local_mocked_bindings() is unavailable and assignInNamespace() +# would fail R CMD check. .plaid_dispatch() therefore takes an injectable `fn`, +# which lets the dispatch logic be tested on CI whether or not plaid is +# installed. Everything needing the real package is skipped explicitly. + +# genes chosen to exist in pbmc_small so UCell does not warn about imputation +mini_gs <- list(B = c("MS4A1", "CD79B", "CD79A", "HLA-DRA", "TCL1A"), + T = c("CD3E", "CD3D", "CD7", "CD8A", "IL7R")) + +# A stand-in for plaid's replaid.* family: honors min.genes/max.genes, returns +# gene sets x cells like the real thing, and records what it was called with. +fake_plaid <- function(X, matG, min.genes = 5, max.genes = 500, alpha = 0, ...) { + keep <- names(matG)[lengths(matG) >= min.genes & lengths(matG) <= max.genes] + m <- matrix(seq_len(length(keep) * ncol(X)), + nrow = length(keep), ncol = ncol(X), + dimnames = list(keep, colnames(X))) + attr(m, "call.args") <- list(min.genes = min.genes, max.genes = max.genes, + alpha = alpha) + m +} + +big_sets <- list(Big = paste0("G", 1:900), + Small = paste0("G", 1:10)) +big_X <- matrix(1, nrow = 1000, ncol = 4, + dimnames = list(paste0("G", 1:1000), paste0("c", 1:4))) + +# -------------------------------------------------------------------------- +# Dispatch logic - no plaid needed +# -------------------------------------------------------------------------- +test_that("the max.genes cap is lifted above the largest gene set", { + out <- escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", + min.size = 5, fn = fake_plaid) + # plaid's own default (500) would have silently dropped the 900-gene set + expect_setequal(colnames(out), names(big_sets)) + expect_equal(dim(out), c(4L, 2L)) +}) + +test_that("output is cells x gene sets with the input dimnames", { + out <- escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", + min.size = 5, fn = fake_plaid) + expect_equal(rownames(out), colnames(big_X)) + expect_equal(colnames(out), names(big_sets)) +}) + +test_that("min.size is translated to plaid's min.genes", { + captured <- NULL + spy <- function(X, matG, min.genes = 5, max.genes = 500, ...) { + captured <<- list(min.genes = min.genes, max.genes = max.genes) + fake_plaid(X, matG, min.genes = min.genes, max.genes = max.genes) + } + escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", min.size = 8, fn = spy) + expect_equal(captured$min.genes, 8L) + expect_gte(captured$max.genes, 900L) +}) + +test_that("backend.args reach the plaid function and typos are rejected", { + captured <- NULL + spy <- function(X, matG, min.genes = 5, max.genes = 500, alpha = 0, ...) { + captured <<- alpha + fake_plaid(X, matG, min.genes = min.genes, max.genes = max.genes) + } + escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", + backend.args = list(alpha = 0.25), fn = spy) + expect_equal(captured, 0.25) + + expect_error( + escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", + backend.args = list(alhpa = 1), fn = fake_plaid), + "Unknown `backend.args`" + ) +}) + +test_that("dropped gene sets warn rather than silently shrinking the output", { + # a user-supplied cap that excludes the 900-gene set + expect_warning( + out <- escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", + backend.args = list(max.genes = 100), + fn = fake_plaid), + "plaid dropped 1 of 2 gene set" + ) + expect_equal(colnames(out), "Small") + + # nothing scored at all is an error, not an empty matrix + expect_error( + escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", + backend.args = list(max.genes = 1), + fn = fake_plaid), + "scored none of the" + ) +}) + +test_that("a plaid result with the wrong number of cells is rejected", { + wrong <- function(X, matG, ...) { + m <- matrix(1, nrow = length(matG), ncol = 2, + dimnames = list(names(matG), c("c1", "c2"))) + m + } + expect_error( + escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", fn = wrong), + "returned 2 cell" + ) +}) + +# -------------------------------------------------------------------------- +# Method / backend resolution +# -------------------------------------------------------------------------- +test_that(".resolve_backend normalizes method names and backends", { + expect_equal(escape:::.resolve_backend("ssgsea", "plaid"), + list(method = "SSGSEA", backend = "plaid")) + expect_equal(escape:::.resolve_backend("AUCell", "native"), + list(method = "AUCELL", backend = "native")) + expect_error(escape:::.resolve_backend("nope", "native"), "Unknown `method`") +}) + +test_that("plaid-only methods route to plaid regardless of backend", { + for (m in c("PLAID", "singscore", "scSE")) { + expect_message(res <- escape:::.resolve_backend(m, "native"), + "provided by the plaid backend") + expect_equal(res$backend, "plaid") + } + # asking for plaid explicitly is silent + expect_silent(escape:::.resolve_backend("scSE", "plaid")) +}) + +test_that("escape.matrix rejects a plaid-only method under the native engine only via routing", { + expect_error(escape.matrix(SeuratObject::pbmc_small, mini_gs, method = "nope"), + "Unknown `method`") +}) + +# -------------------------------------------------------------------------- +# Guards +# -------------------------------------------------------------------------- +test_that("the plaid backend errors informatively when plaid is missing", { + skip_if(requireNamespace("plaid", quietly = TRUE)) + expect_error( + escape.matrix(SeuratObject::pbmc_small, mini_gs, min.size = 0, + backend = "plaid"), + "plaid not installed" + ) +}) + +test_that("dots are refused on the plaid path", { + expect_error( + escape.matrix(SeuratObject::pbmc_small, mini_gs, min.size = 0, + backend = "plaid", maxRank = 100), + "not forwarded when backend" + ) +}) + +test_that("backend.args must be a fully named list", { + expect_error( + escape:::.plaid_dispatch(big_X, big_sets, "ssGSEA", + backend.args = list(1), fn = fake_plaid), + "fully named list" + ) +}) + +# -------------------------------------------------------------------------- +# Provenance and native-path regression +# -------------------------------------------------------------------------- +test_that("scores are stamped with the engine that produced them", { + res <- escape.matrix(SeuratObject::pbmc_small, mini_gs, method = "UCell", + min.size = 0) + prov <- attr(res, "escape.backend") + expect_equal(prov$backend, "native") + expect_equal(prov$method, "UCell") +}) + +test_that("runEscape records provenance on the object", { + skip_if_not_installed("SingleCellExperiment") + sce <- make_toy_sce() + out <- runEscape(sce, gene.sets = toy_spe_sets(), method = "UCell", + min.size = NULL) + prov <- S4Vectors::metadata(out)[["escape_backend"]] + expect_equal(prov$backend, "native") + + obj <- runEscape(SeuratObject::pbmc_small, gene.sets = mini_gs, + method = "UCell", min.size = 0) + expect_equal(SeuratObject::Misc(obj, "escape_backend")$backend, "native") +}) + +test_that("the native path is unchanged by the new arguments", { + a <- escape.matrix(SeuratObject::pbmc_small, mini_gs, min.size = 0) + b <- escape.matrix(SeuratObject::pbmc_small, mini_gs, min.size = 0, + backend = "native", input.assay = "auto") + expect_equal(a, b) +}) + +# -------------------------------------------------------------------------- +# Integration - needs the real package +# -------------------------------------------------------------------------- +test_that("plaid backend reproduces escape's orientation and dimnames", { + skip_if_not_installed("plaid") + pbmc <- SeuratObject::pbmc_small + + for (m in c("ssGSEA", "GSVA", "UCell", "AUCell")) { + res <- escape.matrix(pbmc, mini_gs, method = m, min.size = 0, + backend = "plaid") + expect_equal(dim(res), c(ncol(pbmc), length(mini_gs)), + info = paste("method:", m)) + expect_equal(colnames(res), names(mini_gs), info = paste("method:", m)) + expect_equal(rownames(res), colnames(pbmc), info = paste("method:", m)) + expect_true(all(is.finite(res)), info = paste("method:", m)) + } +}) + +test_that("plaid-only methods return sane matrices", { + skip_if_not_installed("plaid") + pbmc <- SeuratObject::pbmc_small + + for (m in c("PLAID", "singscore", "scSE")) { + res <- escape.matrix(pbmc, mini_gs, method = m, min.size = 0) + expect_equal(dim(res), c(ncol(pbmc), length(mini_gs)), + info = paste("method:", m)) + expect_true(all(is.finite(res)), info = paste("method:", m)) + } +}) + +test_that("no gene sets are dropped at Hallmark scale under defaults", { + skip_if_not_installed("plaid") + pbmc <- SeuratObject::pbmc_small + feats <- rownames(pbmc) + # a set larger than plaid's default max.genes = 500 + wide <- list(Wide = rep(feats, length.out = 600), Narrow = feats[1:10]) + + # the drop warning must not fire - that is what "nothing dropped" means here + expect_no_warning( + res <- suppressMessages( + escape.matrix(pbmc, wide, method = "ssGSEA", min.size = 0, + backend = "plaid") + ) + ) + expect_equal(colnames(res), names(wide)) +}) + +test_that("plaid backend is broadly concordant with the native engines", { + skip_if_not_installed("plaid") + skip_if_not_installed("UCell") + pbmc <- SeuratObject::pbmc_small + + nat <- escape.matrix(pbmc, mini_gs, method = "UCell", min.size = 0) + pla <- escape.matrix(pbmc, mini_gs, method = "UCell", min.size = 0, + backend = "plaid") + + # Deliberately loose. plaid documents replaid.ucell as near-identical to + # UCell, but measured pooled Pearson is ~0.79 on pbmc_small (230 genes) and + # ~0.83 on a 2000-gene simulation - the scores are correlated but on a + # different scale. This guards against the backend breaking outright, not + # against drift, and must not be tightened into an equivalence claim the + # package cannot make. + expect_gt(stats::cor(as.vector(nat), as.vector(pla)), 0.6) +}) + +test_that("backend.args change the scores", { + skip_if_not_installed("plaid") + pbmc <- SeuratObject::pbmc_small + a <- escape.matrix(pbmc, mini_gs, method = "ssGSEA", min.size = 0, + backend = "plaid") + b <- escape.matrix(pbmc, mini_gs, method = "ssGSEA", min.size = 0, + backend = "plaid", backend.args = list(alpha = 0.25)) + expect_false(isTRUE(all.equal(unname(a), unname(b)))) +}) + +test_that("plaid scores feed the downstream workflow", { + skip_if_not_installed("plaid") + pbmc <- SeuratObject::pbmc_small + obj <- runEscape(pbmc, gene.sets = mini_gs, method = "UCell", min.size = 0, + backend = "plaid") + + expect_true("escape" %in% SeuratObject::Assays(obj)) + expect_equal(SeuratObject::Misc(obj, "escape_backend")$backend, "plaid") + expect_no_error(performNormalization(obj, assay = "escape", + gene.sets = mini_gs)) + expect_no_error(performPCA(obj, assay = "escape")) +}) diff --git a/tests/testthat/test-runEscape.R b/tests/testthat/test-runEscape.R index 3d7b794..f5791e7 100644 --- a/tests/testthat/test-runEscape.R +++ b/tests/testthat/test-runEscape.R @@ -8,15 +8,17 @@ mini_gs <- list( pbmc_small <- SeuratObject::pbmc_small get_score <- function(method = "ssGSEA", ...) { - escape.matrix(pbmc_small, - gene.sets = mini_gs, - method = method, - groups = 200, # small chunk for speed - min.size = 0, - normalize = FALSE, - make.positive = FALSE, - min.filter.by = NULL, - BPPARAM = BiocParallel::SerialParam()) + defaults <- list(input.data = pbmc_small, + gene.sets = mini_gs, + method = method, + groups = 200, # small chunk for speed + min.size = 0, + normalize = FALSE, + make.positive = FALSE, + min.filter.by = NULL, + BPPARAM = BiocParallel::SerialParam()) + # `...` must actually reach escape.matrix(), and must win over the defaults + do.call(escape.matrix, utils::modifyList(defaults, list(...))) } # ------------------------------------------------------------- interface ----- @@ -50,10 +52,15 @@ test_that("gene-sets failing min.size are dropped with message", { # --------------------------------------------------- min.expr.cells (global) - test_that("min.expr.cells filters genes globally", { sc0 <- get_score(min.expr.cells = 0) - sc5 <- get_score(min.expr.cells = 0.5) # keep genes in >= 50% of cells - expect_true(is.matrix(sc5) && is.matrix(sc0)) - # dimension equality (gene filter should not affect cell × set shape) - expect_equal(dim(sc0), dim(sc5)) + sc2 <- get_score(min.expr.cells = 0.2) # keep genes in >= 20% of cells + expect_true(is.matrix(sc0) && is.matrix(sc2)) + # the gene filter must not change the cell x gene-set shape ... + expect_equal(dim(sc0), dim(sc2)) + # ... but it must actually change the scores, or it is doing nothing + expect_false(isTRUE(all.equal(unname(sc0), unname(sc2)))) + + # pushing the threshold past every set member is an error, not silence + expect_error(get_score(min.expr.cells = 0.5), "could be matched") }) # --------------------------------------------------------- chunk invariance -- From 7e699d06955a7822208529a96c0550254e3f1c89 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 04:49:30 -0500 Subject: [PATCH 5/7] Update github actions --- .github/workflows/R-CMD-check.yaml | 12 +++++++++++- .github/workflows/test-coverage.yaml | 14 +++++++++++++- 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 9cd83b0..2b4b287 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -45,7 +45,17 @@ jobs: - name: Add GSVA repo run: Rscript -e 'remotes::install_github("rcastelo/GSVA@c9fb985eb555f06a4195260b902747fe0a1ade9f")' - + + # plaid is Bioconductor 3.23 (devel) only, so BiocManager will not resolve + # it on the R release runner. Fall back to source so the plaid backend + # tests execute rather than skip. Note the GitHub HEAD may lag Bioconductor. + - name: Install plaid + run: Rscript -e 'if (!requireNamespace("plaid", quietly = TRUE)) tryCatch(BiocManager::install("plaid", ask = FALSE, update = FALSE), error = function(e) NULL)' -e 'if (!requireNamespace("plaid", quietly = TRUE)) remotes::install_github("bigomics/plaid", upgrade = "never")' + + - name: Install SpatialExperiment + run: Rscript -e 'if (!requireNamespace("SpatialExperiment", quietly = TRUE)) BiocManager::install("SpatialExperiment", ask = FALSE, update = FALSE)' + + - uses: r-lib/actions/check-r-package@v2 with: upload-snapshots: true diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 35544f0..acfcaf7 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -23,9 +23,21 @@ jobs: - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: any::covr + extra-packages: any::covr, any::remotes, any::BiocManager needs: coverage + - name: Add GSVA repo + run: Rscript -e 'remotes::install_github("rcastelo/GSVA@c9fb985eb555f06a4195260b902747fe0a1ade9f")' + + # plaid is Bioconductor 3.23 (devel) only, so BiocManager will not resolve + # it on the R release runner. Fall back to source so the plaid backend + # tests execute rather than skip. Note the GitHub HEAD may lag Bioconductor. + - name: Install plaid + run: Rscript -e 'if (!requireNamespace("plaid", quietly = TRUE)) tryCatch(BiocManager::install("plaid", ask = FALSE, update = FALSE), error = function(e) NULL)' -e 'if (!requireNamespace("plaid", quietly = TRUE)) remotes::install_github("bigomics/plaid", upgrade = "never")' + + - name: Install SpatialExperiment + run: Rscript -e 'if (!requireNamespace("SpatialExperiment", quietly = TRUE)) BiocManager::install("SpatialExperiment", ask = FALSE, update = FALSE)' + # This step is corrected to use the right functions - name: Test coverage and generate report run: | From 8f9848a303a472242f6b398b7062973181ebf4a1 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 04:49:38 -0500 Subject: [PATCH 6/7] Update escape.Rmd --- vignettes/escape.Rmd | 144 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 132 insertions(+), 12 deletions(-) diff --git a/vignettes/escape.Rmd b/vignettes/escape.Rmd index 4350d5b..75026b3 100644 --- a/vignettes/escape.Rmd +++ b/vignettes/escape.Rmd @@ -36,6 +36,8 @@ The core workflow is: 3. (Optional) Normalize for drop-out (```performNormalization()```) 4. Explore with the built-in visualization gallery +Counts can arrive as a plain matrix, a Seurat object, a `SingleCellExperiment`, or a `SpatialExperiment` - spots are treated as cells throughout, so nothing about the workflow changes for spatial data. Scoring runs on escape's own engines by default, with an optional [plaid](https://bigomics.github.io/plaid/) backend for large objects. + # Installation ```{r eval=FALSE} @@ -47,6 +49,13 @@ if (!require("BiocManager", quietly = TRUE)) BiocManager::install("escape") ``` +Two optional packages unlock features used below. Neither is required for the default workflow: + +```{r eval=FALSE} +BiocManager::install("SpatialExperiment") # spatial objects +BiocManager::install("plaid") # the fast scoring backend +``` + Load escape alongside a single-cell container (Seurat or SingleCellExperiment) and a plotting backend: ```{r} @@ -142,7 +151,7 @@ In contrast to ssGSEA and GSVA, AUCell takes the gene rankings for each cell and Please see the following [citation](https://pubmed.ncbi.nlm.nih.gov/28991892/) for more information. -### UCell +### **UCell** UCell calculates a Mann-Whitney U statistic based on the gene rank list. **Importantly**, UCell has a cut-off for ranked genes ($$r_{max}$$) at 1500 - this is per design as drop-out in single-cell can alter enrichment results. This also substantially speeds the calculations up. @@ -155,22 +164,113 @@ $$ Please see the following [citation](https://pubmed.ncbi.nlm.nih.gov/34285779/) for more information. +## Choosing a backend + +Each method above runs on escape's own engine by default. Setting `backend = "plaid"` reroutes the calculation through [plaid](https://bigomics.github.io/plaid/), which reimplements these scores as a single sparse matrix operation and uses substantially less memory. On a 166,000-spot Xenium `SpatialExperiment` with three gene sets, UCell scoring took 28.6 s natively and 6.8 s through plaid, a 4.2x speedup. plaid also provides three scores escape has no native implementation for, which select the plaid backend on their own. + +| Method | Backends | Measured agreement with the native score | Input assay | +|---|---|---|---| +| `ssGSEA` | native, plaid | pooled Pearson ~0.85 | counts (native) / logcounts (plaid) | +| `GSVA` | native, plaid | pooled Pearson ~0.73 | counts (native) / logcounts (plaid) | +| `AUCell` | native, plaid | pooled Pearson ~0.83 | counts (native) / logcounts (plaid) | +| `UCell` | native, plaid | pooled Pearson ~0.83 | counts (native) / logcounts (plaid) | +| `PLAID` | plaid only | n/a - a distinct score | logcounts | +| `singscore` | plaid only | n/a - no native implementation | logcounts | +| `scSE` | plaid only | n/a - no native implementation | logcounts | + +Those correlations come from a simulated 2000-gene, 120-cell matrix with two planted gene sets. On narrower matrices they drop: `pbmc_small` (230 genes) gives roughly 0.79 for UCell and 0.50 for ssGSEA. + + +```{r eval = FALSE} +# native (default) +scores <- escape.matrix(pbmc_small, gene.sets = GS.hallmark, method = "UCell") + +# same method, plaid engine, log-normalized input picked automatically +fast <- escape.matrix(pbmc_small, gene.sets = GS.hallmark, method = "UCell", + backend = "plaid") + +# method-specific tuning goes through backend.args, never through ... +tuned <- escape.matrix(pbmc_small, gene.sets = GS.hallmark, method = "ssGSEA", + backend = "plaid", backend.args = list(alpha = 0.25)) +``` + +Two arguments are easy to confuse. `normalize` is escape's post-hoc drop-out scaling (see `performNormalization()` below). `backend.args$normalize` is plaid's median normalization of the raw scores, applied inside the calculation. They are independent, and turning both on stacks them. + +Install plaid with `BiocManager::install("plaid")`. If it is missing, the native path is unaffected and any `backend = "plaid"` call errors with install instructions. Anyone using a plaid backend should cite Zito *et al.*, *Bioinformatics* 2025, 41(12):btaf621 in addition to the original method paper and escape. + +## Choosing the input assay + +`input.assay` controls which expression matrix is scored. The default `"auto"` reads raw counts for the native backend and log-normalized values for plaid, which is what each expects. `"counts"` and `"logcounts"` resolve to the right place for both object types - the `counts` and `data` layers in Seurat, the `counts` and `logcounts` assays in a SummarizedExperiment - so the same call works across classes. Any other string is taken literally as a layer or assay name. + +Note that the four native methods are rank or CDF based, so a monotone transform of the input generally leaves their scores unchanged. `input.assay` matters most for the plaid paths, which average in log space. + +## Spatial data + +`SpatialExperiment` objects work anywhere a `SingleCellExperiment` does, spots being treated as cells. Nothing special is required - pass the object straight to `runEscape()`. + +To keep this vignette self-contained, we build a small `SpatialExperiment` from `pbmc_small` with made-up coordinates rather than downloading a real slide: + +```{r tidy=FALSE, eval = requireNamespace("SpatialExperiment", quietly = TRUE)} +library(SpatialExperiment) + +set.seed(42) +spe <- SpatialExperiment( + assays = list(counts = counts(sce.pbmc)), + colData = colData(sce.pbmc), + spatialCoords = matrix(runif(2 * ncol(sce.pbmc)), ncol = 2, + dimnames = list(colnames(sce.pbmc), c("x", "y"))) +) + +spe <- runEscape(spe, + gene.sets = escape.gene.sets, + method = "UCell", + min.size = 3, + new.assay.name = "escape") + +spe +``` + +`runEscape()` attaches the scores as an altExp named by **new.assay.name** and leaves the rest of the object alone - `spatialCoords()`, `imgData()`, and `sample_id` all survive unchanged: + +```{r tidy=FALSE, eval = requireNamespace("SpatialExperiment", quietly = TRUE)} +altExpNames(spe) +dim(spatialCoords(spe)) +``` + +The scores themselves carry no spatial information. Coordinates stay on the parent object, so to plot enrichment against position you pull both from `spe`: + +```{r tidy=FALSE, eval = requireNamespace("SpatialExperiment", quietly = TRUE)} +scores <- t(assay(altExp(spe, "escape"))) + +df <- data.frame(spatialCoords(spe), + Proinflammatory = scores[, "Proinflammatory"]) + +ggplot(df, aes(x = x, y = y, color = Proinflammatory)) + + geom_point(size = 2) + + scale_color_gradientn(colors = hcl.colors(7, "inferno")) + + coord_fixed() + + theme_classic() +``` + +Objects carrying several `sample_id` values need no special treatment. Note that all downstream escape functions - `performNormalization()`, `performPCA()`, and the visualization gallery below - read the altExp the same way they would for a `SingleCellExperiment`, so the rest of this vignette applies unchanged. + ## escape.matrix escape has 2 major functions - the first being ```escape.matrix()```, which serves as the backbone of enrichment calculations. Using count-level data supplied from a single-cell object or matrix, ```escape.matrix()``` will produce an enrichment score for the individual cells with the gene sets selected and output the values as a matrix. -**method** +**method** -* AUCell -* GSVA +* AUCell +* GSVA * ssGSEA * UCell +* PLAID, singscore, scSE (these require the plaid backend and select it on their own) -**groups** +**groups** -* The number of cells to calculate at once. +* The number of cells to calculate at once. -**min.size** +**min.size** * The minimum size of detectable genes in a gene set. Gene sets less than the **min.size** will be removed before the calculation. @@ -180,10 +280,22 @@ escape has 2 major functions - the first being ```escape.matrix()```, which serv **make.positive** -* During normalization, whether to shift the enrichment values to a positive range (**TRUE**) or not (**FALSE**). The default value is **FALSE**. +* During normalization, whether to shift the enrichment values to a positive range (**TRUE**) or not (**FALSE**). The default value is **FALSE**. *Cautionary note:* **make.positive** was added to allow for differential analysis downstream of enrichment as some methods may produce negative values. It preserves log-fold change, but ultimately modifies the enrichment values and should be used with caution. +**backend** + +* Which engine computes the scores, **"native"** (the default) or **"plaid"**. See [Choosing a backend](#choosing-a-backend) above. + +**backend.args** + +* A named list of method-specific tuning passed to the underlying plaid function, e.g. `list(alpha = 0.25)`. Names are checked against the target function, so a typo errors rather than being ignored. Ordinary `...` arguments go to the native engines and are refused on the plaid path. + +**input.assay** + +* Which expression matrix to score. **"auto"** (the default) reads counts for the native backend and logcounts for plaid. See [Choosing the input assay](#choosing-the-input-assay). + ```{r tidy = FALSE} enrichment.scores <- escape.matrix(pbmc_small, @@ -261,7 +373,15 @@ pbmc_small <- performNormalization(input.data = pbmc_small, scale.factor = pbmc_small$nFeature_RNA) ``` -```performNormalization()``` has an additional parameter **make.positive**. Across the individual gene sets, if negative normalized enrichment scores are seen, the minimum value is added to all values. For example if the normalized enrichment scores (after the above accounting for drop out) ranges from -50 to 50, **make.positive** will adjust the range to 0 to 100 (by adding 50). This allows for compatible log2-fold change downstream, but can alter the enrichment score interpretation. +`performNormalization()` has an additional parameter **make.positive**. Across the individual gene sets, if negative normalized enrichment scores are seen, the minimum value is added to all values. For example if the normalized enrichment scores (after the above accounting for drop out) ranges from -50 to 50, **make.positive** will adjust the range to 0 to 100 (by adding 50). This allows for compatible log2-fold change downstream, but can alter the enrichment score interpretation. + +Three details are worth knowing: + +* **It always reads raw counts.** The scale factor is the number of genes in each set with a non-zero count, so `performNormalization()` ignores whatever **input.assay** was used to compute the scores. Detection is identical in count and log space, so this is deliberate. +* **Supplied scores win.** If you pass **enrichment.data** explicitly it is used even when the named **assay** also exists on the object, and you get a warning that both were available. +* **Setting `assay = NULL` returns a matrix** rather than attaching a new assay to the object. This is how `escape.matrix(normalize = TRUE)` calls it internally. + +Gene-set names are matched to the enrichment columns literally first, then against the hyphenated form Seurat produces when it builds an assay. That means underscored library names such as `HALLMARK_HYPOXIA` normalize correctly whether the scores live in a Seurat assay, an altExp, or a bare matrix. **** @@ -407,7 +527,7 @@ splitEnrichment(pbmc_small, ## gseaEnrichment -```gseaEnrichment()``` reproduces the two-panel GSEA graphic from Subramanian et al. (2005): +`gseaEnrichment()` reproduces the two-panel GSEA graphic from Subramanian et al. (2005): * Panel A – the running enrichment score (RES) as you “walk” down the ranked list. * Panel B – a rug showing exact positions of each pathway gene. @@ -432,7 +552,7 @@ gseaEnrichment(pbmc_small, ## densityEnrichment -```densityEnrichment()``` is a method to visualize the mean rank position of the gene set features along the total feature space by group. Instead of the classic GSEA running-score, it overlays **kernel-density traces** of the *gene ranks* (1 = most highly expressed/ranked gene) for every group or cluster. High densities at the *left-hand* side mean the pathway is collectively **up-regulated**; peaks on the *right* imply down-regulation. +`densityEnrichment()` is a method to visualize the mean rank position of the gene set features along the total feature space by group. Instead of the classic GSEA running-score, it overlays **kernel-density traces** of the *gene ranks* (1 = most highly expressed/ranked gene) for every group or cluster. High densities at the *left-hand* side mean the pathway is collectively **up-regulated**; peaks on the *right* imply down-regulation. **Anatomy of the plot** @@ -492,7 +612,7 @@ pcaEnrichment(pbmc_small, y.axis = "PC2") ``` -```pcaEnrichment()``` can plot additional information on the principal component analysis. +`pcaEnrichment()` can plot additional information on the principal component analysis. **add.percent.contribution** will add the relative percent contribution of the x and y.axis to total variability observed in the PCA. From 61e2d56413968ed7e8eafd7e233b6e7bfccc6de4 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 04:49:53 -0500 Subject: [PATCH 7/7] Documentation update --- .Rbuildignore | 1 + DESCRIPTION | 14 ++++--- NEWS.md | 23 +++++++++++ R/performNormalization.R | 68 ++++++++++++++++++++++--------- README.md | 9 ++++ inst/WORDLIST | 22 ++++++++++ man/escape.matrix.Rd | 88 ++++++++++++++++++++++++++++++++++++---- man/runEscape.Rd | 47 +++++++++++++++++---- 8 files changed, 233 insertions(+), 39 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index b13432f..ddcb771 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,4 +1,5 @@ ^\.github$ +^\.claude$ ^www$ ^codecov\.yml$ ^.*\.Rproj$ diff --git a/DESCRIPTION b/DESCRIPTION index 10c590e..22070f5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,17 +1,17 @@ Package: escape -Title: Easy single cell analysis platform for enrichment -Version: 2.9.0 +Title: Easy single-cell analysis platform for enrichment +Version: 2.9.1 Authors@R: c( person(given = "Nick", family = "Borcherding", role = c("aut", "cre"), email = "ncborch@gmail.com"), person(given = "Jared", family = "Andrews", role = c("aut"), email = "jared.andrews07@gmail.com"), person(given = "Tobias", family = "Hoch", role = c("ctb"), email = "tobias@hoch.earth"), - person(given = "Alexei", family = "Martsinkovskiy", role = c("ctb"), email = "alexei.martsinkovskiy@gmail.com") + person(given = "Alexei", family = "Martsinkovskiy", role = c("ctb"), email = "alexei.martsinkovskiy@gmail.com"), + person(given = "Cathal", family = "King", role = c("ctb")) ) -Description: A bridging R package to facilitate gene set enrichment analysis (GSEA) in the context of single-cell RNA sequencing. Using raw count information, Seurat objects, or SingleCellExperiment format, users can perform and visualize ssGSEA, GSVA, AUCell, and UCell-based enrichment calculations across individual cells. Alternatively, escape supports use of rank-based GSEA, such as the use of differential gene expression via fgsea. +Description: A bridging R package to facilitate gene set enrichment analysis (GSEA) in the context of single-cell RNA sequencing. Using raw count information, Seurat objects, SingleCellExperiment or SpatialExperiment format, users can perform and visualize ssGSEA, GSVA, AUCell, and UCell-based enrichment calculations across individual cells. An optional plaid backend provides fast reimplementations of these methods plus PLAID, singscore, and scSE scoring. Alternatively, escape supports use of rank-based GSEA, such as the use of differential gene expression via fgsea. License: MIT + file LICENSE Encoding: UTF-8 LazyData: false -RoxygenNote: 7.3.3 biocViews: Software, SingleCell, Classification, Annotation, GeneSetEnrichment, Sequencing, GeneSignaling, Pathways Depends: R (>= 4.1) Imports: @@ -42,12 +42,15 @@ Suggests: knitr, msigdb, patchwork, + plaid, rmarkdown, rlang, + S4Vectors, scran, SeuratObject, Seurat, SingleCellExperiment, + SpatialExperiment, spelling, stringr, testthat (>= 3.0.0), @@ -55,3 +58,4 @@ Suggests: VignetteBuilder: knitr Language: en-US BugReports: https://github.com/BorchLab/escape/issues +Config/roxygen2/version: 8.0.0 diff --git a/NEWS.md b/NEWS.md index 696c3fa..6115957 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,26 @@ +# 2.9.1 + +## NEW FEATURES +* **`plaid` backend**: `escape.matrix()` and `runEscape()` gain `backend = c("native", "plaid")`. Setting `backend = "plaid"` routes `ssGSEA`, `GSVA`, `UCell`, and `AUCell` to the `replaid.*` family for large speed and memory gains. `plaid` is a `Suggests`, so the native path is unaffected if it is not installed. **These are approximations, and the gap is wider than the plaid documentation suggests** — measured pooled Pearson agreement with the native scores was roughly 0.85 (ssGSEA), 0.83 (UCell, AUCell) and 0.73 (GSVA) on a 2000-gene simulation, and lower on narrower matrices. See `?escape.matrix` and the vignette before substituting one backend for the other. +* **Three new methods**: `method = "PLAID"`, `"singscore"`, and `"scSE"`. These have no native implementation and select the `plaid` backend automatically. +* **`backend.args`**: a named list carrying method-specific tuning (`alpha`, `tau`, `rowtf`, `aucMaxRank`, `rmax`, `nsmooth`, `stats`, `chunk`, `removeLog2`, `scoreMean`). Names are validated against the target function, so typos error instead of being silently swallowed. Note `backend.args$normalize` is `plaid`'s median normalization of the scores and is unrelated to escape's `normalize`. +* **`input.assay`**: choose which expression matrix to score. `"auto"` (default) keeps raw counts for the native backend and picks log-normalized values for `plaid`. `"counts"` and `"logcounts"` map to the correct layer for both Seurat and `SummarizedExperiment`-derived objects. +* **`SpatialExperiment` support** is now explicit and tested, including `spatialCoords()`, `imgData()`, and `sample_id` preservation through `runEscape()`. +* Enrichment scores carry an `escape.backend` attribute recording the engine, method, and `plaid` version; `runEscape()` stores the same record in `metadata()` (SCE) or `Misc()` (Seurat). + +## BUG FIXES +* **`escape.matrix(normalize = TRUE)` failed on every `SummarizedExperiment`-derived input** with `unable to find an inherited method for function 'assay' for signature 'x = "NULL"'`. `.pull.Enrich()` used `altExp(sc)[[name]]`, which indexes the `colData` of the *first* altExp rather than selecting the named one, and returned `NULL`. Now uses `altExp(sc, name)`. This also fixes `performNormalization()` and `performPCA()` on `SingleCellExperiment` and `SpatialExperiment` objects. Reported by Cathal King, who also provided the Xenium test data. +* **Underscored gene-set names were silently dropped during normalization.** `performNormalization()` rewrote every `_` to `-` to match Seurat's feature-name coercion, which meant `HALLMARK_*`, `GO_*`, and `REACTOME_*` sets matched nothing for matrix, `SingleCellExperiment`, or `SpatialExperiment` input. Matching is now literal first with the Seurat-mangled form as a fallback, and per-set scale factors are aligned to the enrichment columns explicitly so a partial match can never divide the wrong column. +* Requesting an assay that does not exist now reports what was asked for and what is available, instead of failing with an S4 dispatch error. +* `performNormalization()` gives a supplied `enrichment.data` precedence over scores stored on the object, and warns when both are present. +* `%||%` was used but never defined; it only exists in base R >= 4.4 while the package declares R >= 4.1. +* `runEscape()` now calls `escape.matrix()` with named rather than positional arguments. + +## ENHANCEMENTS +* Removed duplicate internal definitions of `.split_cols()`, `.match_summary_fun()`, and `.filter_genes()`. +* `escape.matrix()` validates `method` up front rather than failing inside the scoring switch. +* Test suite: added a toy `SpatialExperiment` fixture built at run time (no new package data), regression coverage for all three bug fixes, and `plaid` dispatch tests that run whether or not `plaid` is installed. + # 2.7.3 ## BUG FIXES diff --git a/R/performNormalization.R b/R/performNormalization.R index 6b1e830..e7a2e45 100644 --- a/R/performNormalization.R +++ b/R/performNormalization.R @@ -7,13 +7,17 @@ #' #' @param input.data A raw-counts matrix (genes x cells), a #' \link[SeuratObject]{Seurat} object, or a -#' \link[SingleCellExperiment]{SingleCellExperiment}. Gene identifiers must +#' \link[SingleCellExperiment]{SingleCellExperiment} (including a +#' \link[SpatialExperiment]{SpatialExperiment}). Gene identifiers must #' match those in \code{gene.sets}. #' @param enrichment.data Matrix. Output of \code{\link{escape.matrix}} or #' \code{NULL} if enrichment scores are already stored in \code{input.data}. +#' When supplied it takes precedence over anything stored in +#' \code{input.data}. #' @param assay Character. Name of the assay holding enrichment scores when #' \code{input.data} is a single-cell object. Default is \code{"escape"}. -#' Ignored otherwise. +#' Ignored when \code{input.data} is a matrix. Set to \code{NULL} to return +#' the normalized matrix rather than attaching it to the object. #' @param gene.sets A named list of character vectors, the result of #' \code{\link{getGeneSets}}, or the built-in data object #' \code{\link{escape.gene.sets}}. List names must match column names in the @@ -39,8 +43,16 @@ #' assay = "escape", #' gene.sets = gs) #' -#' @return If `input.data` is an object, the same object with a new assay -#' "_normalized". Otherwise a matrix of normalized scores. +#' @section Which expression values are used: +#' The per-cell scale factor is the number of genes from each set with a +#' \strong{non-zero raw count}, so this function always reads the \code{counts} +#' assay regardless of any \code{input.assay} used when the scores were +#' computed. Detection is identical in count and log space, so this is +#' deliberate rather than an oversight. +#' +#' @return If `input.data` is an object and `assay` is not `NULL`, the same +#' object with a new assay "_normalized". Otherwise a matrix of +#' normalized scores. #' @export performNormalization <- function(input.data, @@ -62,28 +74,40 @@ performNormalization <- function(input.data, } } else if (.is_sce(input.data)) { if (requireNamespace("SingleCellExperiment", quietly = TRUE)) { - assay.present <- assay %in% names(SingleCellExperiment::altExps(input.data)) + assay.present <- assay %in% SingleCellExperiment::altExpNames(input.data) } else { warning("SingleCellExperiment package is required but not installed.") } } } - - enriched <- if (assay.present) .pull.Enrich(input.data, assay) else enrichment.data - if (is.null(enriched)) { + + ## Supplied scores always win - never re-pull something the caller handed in. + if (!is.null(enrichment.data)) { + if (assay.present) + warning("Both `enrichment.data` and assay '", assay, "' are available; ", + "using `enrichment.data`.", call. = FALSE) + enriched <- enrichment.data + } else if (assay.present) { + enriched <- .pull.Enrich(input.data, assay) + } else { + ## nothing to work with - say what was asked for and what exists + if (!is.null(assay) && .is_seurat_or_sce(input.data)) { + avail <- if (.is_seurat(input.data)) SeuratObject::Assays(input.data) + else SingleCellExperiment::altExpNames(input.data) + .stop_missing_assay(assay, avail, "enrichment assay", + "`input.data`, and `enrichment.data` was not supplied") + } stop("Could not obtain enrichment matrix, please set `assay` or supply `enrichment.data`.") } - + ## 2. Validate / derive scale factors ---------------------------------- if (!is.null(scale.factor) && length(scale.factor) != nrow(enriched)) stop("Length of 'scale.factor' must match number of cells.") - + if (is.null(scale.factor)) { - egc <- .GS.check(gene.sets) - names(egc) <- gsub("_", "-", names(egc), fixed = TRUE) - egc <- egc[names(egc) %in% colnames(enriched)] - if (!length(egc)) stop("None of the supplied gene sets match enrichment columns.") - + ## one gene set per enrichment column, in column order + egc <- .match_sets_to_cols(.GS.check(gene.sets), colnames(enriched)) + ## counts matrix (genes x cells) - drop after use to save RAM cnts <- .cntEval(input.data, assay = "RNA", type = "counts") message("Computing expressed-gene counts per cell...") @@ -92,7 +116,12 @@ performNormalization <- function(input.data, vec[vec == 0] <- 1L # avoid /0 vec })) + colnames(scale.mat) <- names(egc) rm(cnts) + + ## alignment is guaranteed by .match_sets_to_cols(); assert it anyway so a + ## future refactor cannot silently divide the wrong column + stopifnot(identical(colnames(scale.mat), colnames(enriched))) ## optionally split large matrices to spare memory chunksize <- if (is.null(groups)) nrow(enriched) else min(groups, nrow(enriched)) sf.split <- .split_rows(scale.mat, chunk.size = chunksize) @@ -120,9 +149,10 @@ performNormalization <- function(input.data, } ## 6. Return ------------------------------------------------------------ - if (.is_seurat_or_sce(input.data)) { - input.data <- .adding.Enrich(input.data, normalized, paste0(assay %||% "escape", "_normalized")) - } else { - normalized + ## `assay = NULL` means "hand the matrix back" - that is how escape.matrix() + ## calls this, and it avoids attaching an altExp only to pull it straight off. + if (.is_seurat_or_sce(input.data) && !is.null(assay)) { + return(.adding.Enrich(input.data, normalized, paste0(assay, "_normalized"))) } + normalized } diff --git a/README.md b/README.md index aeaa4f3..c5374a4 100644 --- a/README.md +++ b/README.md @@ -21,9 +21,18 @@ Single-cell sequencing (SCS) is a fundamental technology in investigating a dive * Single-sample Gene Set Enrichment Analysis (ssGSEA) - [citation](https://pubmed.ncbi.nlm.nih.gov/19847166/) * AUCell - [citation](https://pubmed.ncbi.nlm.nih.gov/28991892/) * UCell -[citation](https://pubmed.ncbi.nlm.nih.gov/34285779/) +* PLAID, singscore and scSE via the optional [plaid](https://bigomics.github.io/plaid/) backend ([source](https://github.com/bigomics/plaid), [citation](https://pubmed.ncbi.nlm.nih.gov/41223139/)) More information on each method is available in the *escape* manual for ```escape.matrix()``` and the citation links. If using these methods, users should cite the original works as well. +#### Backends + +Setting ```backend = "plaid"``` reroutes GSVA, ssGSEA, AUCell and UCell through *plaid*, which is faster and lighter on memory for large objects (4.2x on a 166,000-spot Xenium object). These are approximations rather than exact reproductions, and measured agreement with the native scores is lower than *plaid*'s own documentation suggests - see the vignette before substituting one backend for the other. The default backend remains ```"native"```. + +#### Object types + +Raw count matrices, [Seurat](https://satijalab.org/seurat/), [SingleCellExperiment](https://bioconductor.org/books/release/OSCA/book-contents.html#basics), and [SpatialExperiment](https://bioconductor.org/packages/SpatialExperiment/) objects are all supported. For spatial objects, ```runEscape()``` attaches scores as an altExp and leaves ```spatialCoords()```, ```imgData()``` and ```sample_id``` untouched. + ### Installation #### Install Via GitHub diff --git a/inst/WORDLIST b/inst/WORDLIST index 9602499..0c0fa33 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -5,7 +5,9 @@ BIOCARTA BPPARAM Bcells BioC +Bioc BiocParallel +Bioinformatics CGN CGP CMD @@ -31,6 +33,7 @@ NES NG Nebulosa OpenMP +PLAID PNAS Parallelization Precomputed @@ -41,7 +44,9 @@ SYM ScoreSignatures SeuratObject SingleCellExperiment +SpatialExperiment Subramanian +SummarizedExperiment Tcells TukeyHSD UCell @@ -49,10 +54,16 @@ Vishwakarma Visualisation Visualising Voigt +Xenium +Zito al +altExp args +backend +backends bioconductor biocparallel +btaf calcAUC centred cnet @@ -90,6 +101,7 @@ heatmapEnrichment hexbin https ident +imgData jk leadingEdge limma @@ -97,12 +109,14 @@ linewidth lm loadings logFC +logcounts masterPCAPlot microenvironment msigdb msigdbr multithread musculus +natively ncbi nih nlm @@ -114,27 +128,35 @@ pcaEnrichment performNormalization performPCA phenotypes +plaid plyr pubmed pval reclustering +reimplementations +reimplements rescaling ridgeEnrichment rlang +roxygen runEscape runPCA scRNA +scSE scater's scatterEnrichment singScore singscore +spatialCoords splitEnrichment ssGSEA standardises stringr stromal subcollection +summarise summarization +theming tibble tidyverse vectorised diff --git a/man/escape.matrix.Rd b/man/escape.matrix.Rd index 4c7dd08..f2c931a 100644 --- a/man/escape.matrix.Rd +++ b/man/escape.matrix.Rd @@ -15,13 +15,17 @@ escape.matrix( min.expr.cells = 0, min.filter.by = NULL, BPPARAM = NULL, - ... + ..., + backend = c("native", "plaid"), + backend.args = list(), + input.assay = "auto" ) } \arguments{ \item{input.data}{A raw-counts matrix (genes x cells), a \link[SeuratObject]{Seurat} object, or a -\link[SingleCellExperiment]{SingleCellExperiment}. Gene identifiers must +\link[SingleCellExperiment]{SingleCellExperiment} (including a +\link[SpatialExperiment]{SpatialExperiment}). Gene identifiers must match those in \code{gene.sets}.} \item{gene.sets}{A named list of character vectors, the result of @@ -30,11 +34,17 @@ match those in \code{gene.sets}.} result.} \item{method}{Character. Scoring algorithm (case-insensitive). One of -\code{"GSVA"}, \code{"ssGSEA"}, \code{"UCell"}, or \code{"AUCell"}. -Default is \code{"ssGSEA"}.} +\code{"GSVA"}, \code{"ssGSEA"}, \code{"UCell"}, \code{"AUCell"}, +\code{"PLAID"}, \code{"singscore"}, or \code{"scSE"}. The last three are +available only through \code{backend = "plaid"} and select it +automatically. Default is \code{"ssGSEA"}.} \item{groups}{Integer. Number of cells per processing chunk. Larger values -reduce overhead but increase memory usage. Default is \code{1000}.} +reduce overhead but increase memory usage. Default is \code{1000}. +Meaning depends on the backend: chunk size for the \pkg{BiocParallel} loop +when \code{backend = "native"}, forwarded to \code{plaid::plaid(chunk=)} +for \code{method = "PLAID"}, and ignored for the \code{replaid.*} paths +(use \code{backend.args$chunk} there).} \item{min.size}{Integer or \code{NULL}. Minimum number of genes from a set that must be detected in the expression matrix for that set to be scored. @@ -59,9 +69,30 @@ which the \code{min.expr.cells} rule is applied. Default is \code{NULL}.} \item{BPPARAM}{A \pkg{BiocParallel} parameter object describing the parallel backend. Default is \code{NULL} (serial execution).} -\item{...}{Extra arguments passed verbatim to the chosen back-end scoring +\item{...}{Extra arguments passed verbatim to the chosen native scoring function (\code{gsva()}, \code{ScoreSignatures_UCell()}, or -\code{AUCell_calcAUC()}).} +\code{AUCell_calcAUC()}). Not forwarded when \code{backend = "plaid"} - +use \code{backend.args} there.} + +\item{backend}{Character. Scoring engine, \code{"native"} (default) or +\code{"plaid"}. Ignored for methods that only exist in \pkg{plaid}.} + +\item{backend.args}{Named list of method-specific tuning arguments passed to +the underlying \pkg{plaid} function, e.g. \code{alpha}, \code{tau}, +\code{rowtf}, \code{aucMaxRank}, \code{rmax}, \code{nsmooth}, +\code{stats}, \code{chunk}, \code{removeLog2}, \code{scoreMean}. Names are +validated against the target function. Note that +\code{backend.args$normalize} is \pkg{plaid}'s median normalization of the +scores and is unrelated to escape's \code{normalize} argument. Default is +\code{list()}.} + +\item{input.assay}{Character. Which expression matrix to score. +\code{"auto"} (default) reads raw counts for the native backend and +log-normalized values for \pkg{plaid}. \code{"counts"} and +\code{"logcounts"} map to the right layer for both \pkg{Seurat} +(\code{counts} / \code{data}) and \pkg{SummarizedExperiment}-derived +objects, including \link[SpatialExperiment]{SpatialExperiment}. Any other +string is taken literally.} } \value{ A numeric matrix with one row per cell and one column per gene set, @@ -83,9 +114,52 @@ lazily, keeping them in the package's \strong{Suggests} field. \item{\code{"ssGSEA"}}{Single-sample GSEA.} \item{\code{"UCell"}}{Rank-based UCell scoring.} \item{\code{"AUCell"}}{Area-under-the-curve ranking score.} + \item{\code{"PLAID"}}{Average log-intensity of set members (plaid only).} + \item{\code{"singscore"}}{Rank-based singscore (plaid only).} + \item{\code{"scSE"}}{Single-cell signature explorer score (plaid only).} } } +\section{Backends}{ + +The first four methods run on escape's own engines by default +(\code{backend = "native"}). Setting \code{backend = "plaid"} routes them to +\pkg{plaid}'s \code{replaid.*} family instead, which is substantially faster +and lighter on memory for large objects. \code{"PLAID"}, \code{"singscore"} +and \code{"scSE"} have no native implementation and always use \pkg{plaid}. + +\strong{The plaid backend approximates rather than reproduces, and the gap is +wider than the plaid documentation suggests.} On a simulated 2000-gene, +120-cell matrix, pooled Pearson correlation between the native and plaid +scores for the same method was roughly 0.85 (\code{ssGSEA}), 0.83 +(\code{UCell}, \code{AUCell}) and 0.73 (\code{GSVA}); on the 230-gene +\code{pbmc_small} it was lower still. The scores are also on different +scales, not just noisier. Notably \code{replaid.ssgsea} is documented as +exact at \code{alpha = 0}, but compared directly against +\code{GSVA::gsva()} on identical input it correlated at 0.85, not 1. The +input assay is not the cause - these are rank-based scores, and counts +versus logcounts correlate at exactly 1. + +Treat \code{backend = "plaid"} as a fast screen, not as a drop-in +replacement. Do not mix backends within an analysis, and do not compare +plaid scores against previously published \pkg{escape} results. Tuning that +may narrow the gap for individual methods: +\describe{ + \item{\code{ssGSEA}}{\code{backend.args = list(alpha = 0)}.} + \item{\code{GSVA}}{the empirical CDF row transform is approximated by a + z-transform (\code{rowtf = "z"}); pass + \code{backend.args = list(rowtf = "ecdf")} for the slower exact form.} + \item{\code{UCell}}{\code{backend.args = list(rmax = ...)} shifts the score + scale but did not change rank agreement in testing.} + \item{\code{scSE}}{plaid documents a match to the original with + \code{backend.args = list(removeLog2 = TRUE, scoreMean = FALSE)}.} +} + +Users of the plaid backend should cite Zito \emph{et al.}, \emph{Bioinformatics} +2025, 41(12):btaf621 in addition to the original method paper and +\pkg{escape}. +} + \examples{ gs <- list(Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"), Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A")) diff --git a/man/runEscape.Rd b/man/runEscape.Rd index fabf090..69bfb60 100644 --- a/man/runEscape.Rd +++ b/man/runEscape.Rd @@ -7,7 +7,7 @@ runEscape( input.data, gene.sets, - method = c("ssGSEA", "GSVA", "UCell", "AUCell"), + method = c("ssGSEA", "GSVA", "UCell", "AUCell", "PLAID", "singscore", "scSE"), groups = 1000, min.size = 5, normalize = FALSE, @@ -16,13 +16,17 @@ runEscape( min.expr.cells = 0, min.filter.by = NULL, BPPARAM = NULL, - ... + ..., + backend = c("native", "plaid"), + backend.args = list(), + input.assay = "auto" ) } \arguments{ \item{input.data}{A raw-counts matrix (genes x cells), a \link[SeuratObject]{Seurat} object, or a -\link[SingleCellExperiment]{SingleCellExperiment}. Gene identifiers must +\link[SingleCellExperiment]{SingleCellExperiment} (including a +\link[SpatialExperiment]{SpatialExperiment}). Gene identifiers must match those in \code{gene.sets}.} \item{gene.sets}{A named list of character vectors, the result of @@ -31,11 +35,17 @@ match those in \code{gene.sets}.} result.} \item{method}{Character. Scoring algorithm (case-insensitive). One of -\code{"GSVA"}, \code{"ssGSEA"}, \code{"UCell"}, or \code{"AUCell"}. -Default is \code{"ssGSEA"}.} +\code{"GSVA"}, \code{"ssGSEA"}, \code{"UCell"}, \code{"AUCell"}, +\code{"PLAID"}, \code{"singscore"}, or \code{"scSE"}. The last three are +available only through \code{backend = "plaid"} and select it +automatically. Default is \code{"ssGSEA"}.} \item{groups}{Integer. Number of cells per processing chunk. Larger values -reduce overhead but increase memory usage. Default is \code{1000}.} +reduce overhead but increase memory usage. Default is \code{1000}. +Meaning depends on the backend: chunk size for the \pkg{BiocParallel} loop +when \code{backend = "native"}, forwarded to \code{plaid::plaid(chunk=)} +for \code{method = "PLAID"}, and ignored for the \code{replaid.*} paths +(use \code{backend.args$chunk} there).} \item{min.size}{Integer or \code{NULL}. Minimum number of genes from a set that must be detected in the expression matrix for that set to be scored. @@ -63,9 +73,30 @@ which the \code{min.expr.cells} rule is applied. Default is \code{NULL}.} \item{BPPARAM}{A \pkg{BiocParallel} parameter object describing the parallel backend. Default is \code{NULL} (serial execution).} -\item{...}{Extra arguments passed verbatim to the chosen back-end scoring +\item{...}{Extra arguments passed verbatim to the chosen native scoring function (\code{gsva()}, \code{ScoreSignatures_UCell()}, or -\code{AUCell_calcAUC()}).} +\code{AUCell_calcAUC()}). Not forwarded when \code{backend = "plaid"} - +use \code{backend.args} there.} + +\item{backend}{Character. Scoring engine, \code{"native"} (default) or +\code{"plaid"}. Ignored for methods that only exist in \pkg{plaid}.} + +\item{backend.args}{Named list of method-specific tuning arguments passed to +the underlying \pkg{plaid} function, e.g. \code{alpha}, \code{tau}, +\code{rowtf}, \code{aucMaxRank}, \code{rmax}, \code{nsmooth}, +\code{stats}, \code{chunk}, \code{removeLog2}, \code{scoreMean}. Names are +validated against the target function. Note that +\code{backend.args$normalize} is \pkg{plaid}'s median normalization of the +scores and is unrelated to escape's \code{normalize} argument. Default is +\code{list()}.} + +\item{input.assay}{Character. Which expression matrix to score. +\code{"auto"} (default) reads raw counts for the native backend and +log-normalized values for \pkg{plaid}. \code{"counts"} and +\code{"logcounts"} map to the right layer for both \pkg{Seurat} +(\code{counts} / \code{data}) and \pkg{SummarizedExperiment}-derived +objects, including \link[SpatialExperiment]{SpatialExperiment}. Any other +string is taken literally.} } \value{ The input single-cell object with an additional assay containing the