diff --git a/DESCRIPTION b/DESCRIPTION
index 85f3527..72072a2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: SpaceTrooper
Type: Package
Title: SpaceTrooper performs Quality Control analysis of Image-Based spatial
-Version: 1.1.8
+Version: 1.1.9
Authors@R:
c(person("Dario", "Righelli",
email="dario.righelli@gmail.com",
@@ -23,7 +23,7 @@ Description: SpaceTrooper performs Quality Control analysis using data driven
License: MIT + file LICENSE
Encoding: UTF-8
Depends:
- R (>= 4.4.0),
+ R (>= 4.6.0),
SpatialExperiment
Imports:
DropletUtils,
diff --git a/NAMESPACE b/NAMESPACE
index 5911763..a5e4a18 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,10 +1,11 @@
# Generated by roxygen2: do not edit by hand
-export(.applyQScoreModel)
export(.getActiveGeometryName)
export(.renameGeometry)
export(.setActiveGeometry)
export(addPolygonsToSPE)
+export(applyQCScoreModel)
+export(applyQScoreModel)
export(checkOutliers)
export(computeAreaFromPolygons)
export(computeAspectRatioFromPolygons)
@@ -12,7 +13,10 @@ export(computeCenterFromPolygons)
export(computeLambda)
export(computeMissingMetricsMerfish)
export(computeMissingMetricsXenium)
+export(computeOutliersQCScore)
export(computeOutliersQScore)
+export(computeQCScore)
+export(computeQCScoreFlags)
export(computeQScore)
export(computeQScoreFlags)
export(computeSpatialOutlier)
@@ -132,3 +136,4 @@ importFrom(stats,complete.cases)
importFrom(stats,model.matrix)
importFrom(stats,predict)
importFrom(stats,quantile)
+importFrom(stats,terms)
diff --git a/NEWS.md b/NEWS.md
index f7aff73..2e92b31 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,8 +1,38 @@
-# Changes in version 1.1.8
-
-* implementing Quality Score computation with custom formula
-* fixing naming of spacetrooper utilities vignette
-* adding citation file with biorxiv paper
+# Changes in version 1.1.9
+
+* Added the canonical `computeQScore()`, `computeQScoreFlags()`,
+ `computeOutliersQScore()`, and `applyQScoreModel()` APIs.
+* Retained `computeQCScore()`, `computeQCScoreFlags()`,
+ `computeOutliersQCScore()`, and `applyQCScoreModel()` as deprecated
+ compatibility APIs. They preserve historical arguments, `QC_score`,
+ `low_qcscore`, `low_threshold_qcscore`, and `QCScore_model*` metadata.
+* Restored historical plotting arguments (`size`, `alpha`, `alphaNumbers`,
+ `csize`, `calpha`, `mapNumbersCol`, and `mapAlphaNumbers`) as aliases for
+ the canonical plotting arguments. Canonical arguments take precedence when
+ both forms are supplied.
+* Restored `getModelFormula(formulaVars, verbose=FALSE)` compatibility while
+ accepting `metricList` as a named replacement.
+* Restored the public `computeLambda(trainDF, modelFormula)` signature,
+ complete-case filtering, and support for both `qcscore_train` and
+ `QScore_train`.
+* User-supplied `modelFormula` values are now fitted without adding, removing,
+ or rebuilding terms. Supported subsets, additive formulas, and selected
+ interactions are preserved.
+* Supported base QC metrics selected by a custom `modelFormula` now determine
+ the outliers, good and bad cells, and `QScore_train` used for model training;
+ omitted metrics no longer contribute to training-label construction.
+* Formula predictors are limited to `log2SignalDensity`, `Area_um`,
+ `log2AspectRatio`, and `log2Ctrl_total_ratio`. Unsupported predictors and
+ transformations now produce informative errors, and CosMx-only border terms
+ are checked against dataset technology.
+* Model transfer now has a canonical public API, preserves training matrix
+ column order, supports validated custom formulas, and retains the historical
+ deprecated interface.
+* Fixed the SpaceTrooper utilities vignette name and added the bioRxiv
+ citation.
+* Follow-up: intercept handling remains unchanged in this release and should be
+ reviewed consistently across training, lambda selection, prediction, and
+ model transfer.
# Changes in version 1.1.7
diff --git a/R/QC.R b/R/QC.R
index 7bd3e76..e63cde1 100644
--- a/R/QC.R
+++ b/R/QC.R
@@ -352,10 +352,10 @@ computeThresholdFlags <- function(spe, totalThreshold=0,
#' @param trainDF `data.frame`
#' A data frame for QS model training that must include:
#' Predictor columns: All columns referenced in the formula returned by `getModelFormula()`.
-#' `qscore_train` A binary (0/1) response vector to be modeled.
-#' @param modelFormula `character`
-#' A character string representing the model formula
-#' `~ log2SignalDensity + ...`, as returned by `getModelFormula()`.
+#' `QScore_train` or historical `qcscore_train`, containing the binary
+#' (0/1) response to be modeled.
+#' @param modelFormula A one-sided formula or character string describing the
+#' supported predictor terms.
#'
#' @return
#' `numeric`
@@ -363,25 +363,48 @@ computeThresholdFlags <- function(spe, totalThreshold=0,
#' `cv.glmnet` that minimizes the cross-validation error.
#'
#' @details
-#' Internally, the function:
-#' constructs the design matrix via \code{model.matrix()},
-#' runs k-fold cross-validation of ridge logistic regression using `cv.glmnet` with `alpha = 0`,
-#' extracts and returns `ridge_cv$lambda.min`.
+#' Complete cases for the formula variables and response are retained. The
+#' function constructs a model matrix and passes it to an internal
+#' matrix-based implementation that runs ridge logistic-regression
+#' cross-validation using `cv.glmnet` with `alpha = 0`.
#'
#' @examples
#' example(computeTrainDF)
-#' modform <- getModelFormula(metadata(spe)$formula_variables)
+#' modform <- getModelFormula(names(metadata(spe)$formula_variables))
#' best_lambda <- computeLambda(df_train, modform)
#' print(best_lambda)
#'
#'
#' @export
-computeLambda <- function(modelMatrix, trainDF) {
- ridge_cv <- cv.glmnet(modelMatrix, trainDF$QScore_train,
- family="binomial", alpha=0, lambda=NULL)
- bestLambda <- ridge_cv$lambda.min
- return(bestLambda)
+computeLambda <- function(trainDF, modelFormula) {
+ stopifnot(is.data.frame(trainDF))
+ model_formula <- stats::as.formula(modelFormula)
+ response <- .getQScoreResponse(trainDF)
+ train_ok <- .filterCompleteModelCases(
+ df=trainDF,
+ modelFormula=model_formula,
+ response=response,
+ context="training cells for lambda selection"
+ )
+ train_df <- trainDF[train_ok, , drop=FALSE]
+ response <- response[train_ok]
+ model_matrix <- stats::model.matrix(
+ model_formula,
+ data=train_df
+ )
+ .computeLambda(modelMatrix=model_matrix, response=response)
+}
+
+.computeLambda <- function(modelMatrix, response) {
+ ridge_cv <- glmnet::cv.glmnet(
+ modelMatrix,
+ response,
+ family="binomial",
+ alpha=0,
+ lambda=NULL
+ )
+ ridge_cv$lambda.min
}
#' computeQScore
@@ -414,9 +437,9 @@ computeLambda <- function(modelMatrix, trainDF) {
#' metrics in the formula based on their availability in the `colData` of the
#' `SpatialExperiment` object.
#'
-#' Inclusion of metrics in the formula depends also on the number of available
-#' outliers. If the number of outliers for each metric is less than 0.1\% out of the
-#' entire dataset, the metric will be excluded from the QS formula.
+#' For the default formula, inclusion of metrics also depends on the number of
+#' available outliers. If the number of outliers for a metric is less than
+#' 0.1\% of the dataset, that metric is excluded from the default QS formula.
#'
#' - Model fitting: ridge (L2) logistic regression is fitted (via `glmnet`) on
#' the balanced training set. The function uses `trainModel()` for fitting
@@ -428,19 +451,28 @@ computeLambda <- function(modelMatrix, trainDF) {
#' `computeQScore`. Otherwise, a fixed value of lambda previously computed with
#' `computeLambda` preceeded by `computeTrainDF` and `getModelFormula` can be set.
#'
-#' - Model formula details: the model formula is automatically generated as
-#' follows:
+#' - Model formula details: the most complete default formula is:
#' `~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2`.
-#' When user-provided, the formula must follow the same default syntax and
-#' removed (or added) terms should be written exactly as in the default formula,
-#' e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))` must have spaces
-#' around the `*` and `<` operators.
-#' In any case, metrics with insufficient outliers (less than 0.1\% of the dataset)
-#' will be excluded from the QS formula.
+#' A user-supplied formula is fitted as supplied: additive structures, selected
+#' interactions, `:` and `*` are preserved. The formula must include
+#' `log2SignalDensity` and may use any subset of the other supported predictors:
+#' `Area_um`,
+#' `log2AspectRatio`, and `log2Ctrl_total_ratio`. Arbitrary predictors and
+#' transformations are rejected. The supported CosMx-only border expression is
+#' `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`.
+#'
+#' When a custom `modelFormula` is supplied, the supported base QC metrics in
+#' that formula are also used to compute outliers, define good and bad training
+#' cells, and construct `QScore_train`. Interactions do not create additional
+#' training metrics. In the supported border expression, `dist_border` is an
+#' auxiliary variable and `log2AspectRatio` is the corresponding training QC
+#' metric. Supported metrics omitted from the custom formula do not contribute
+#' to training-label construction. When `modelFormula=NULL`, the established
+#' default metric-selection and training-label workflow is unchanged.
#'
#' The computed model output is stored in `metadata(spe)$QScore_model`.
#' For inspection and model coefficient transfer between datasets
-#' see also \code{\link{.applyQScoreModel}}. However, transferring models between
+#' see also \code{\link{applyQScoreModel}}. However, transferring models between
#' datasets is not recommended. As described in the paper, QS model training is
#' dataset-specific and does not generalize well across datasets. Moreover,
#' this step is computationally efficient. Please, refer to the paper for
@@ -449,15 +481,15 @@ computeLambda <- function(modelMatrix, trainDF) {
#' @param spe A `SpatialExperiment` object with spatial omics data.
#' @param verbose logical for having a verbose output. Default is FALSE.
#' @param bestLambda the best lambda typically computed using `computeLambda`.
-#' @param modelFormula a character string representing the formula to be used for
-#' training the model. If NULL, the formula is automatically generated
-#' based on the available metrics and their outliers in the dataset. See details
-#' for more information.
+#' @param modelFormula A one-sided formula or character string. If `NULL`, the
+#' default formula is generated from available supported metrics and their
+#' outliers. If supplied, its supported terms and interactions are fitted
+#' without rewriting. See Details.
#' @return The `SpatialExperiment` object with added Quality Score in `colData`.
#' @export
#' @importFrom dplyr case_when filter mutate distinct pull
#' @importFrom glmnet glmnet cv.glmnet
-#' @importFrom stats as.formula model.matrix quantile predict coef
+#' @importFrom stats as.formula model.matrix quantile predict coef terms
#' @examples
#' example(spatialPerCellQC)
#' set.seed(1998)
@@ -470,25 +502,38 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE
" cells with 0 counts were found. These cells will be removed."))
spe <- spe[,spe$total > 0]
}
- metricList <- c("log2SignalDensity", "Area_um",
- "log2AspectRatio",
- "log2Ctrl_total_ratio")
-
- if (!is.null(modelFormula)) {
- model_formula <- modelFormula
- metricList <- attr(terms(as.formula(modelFormula)), "term.labels")
- metricList <- metricList[!grepl(":", metricList, fixed=TRUE)]
- if("I(abs(log2AspectRatio) * as.numeric(dist_border < 50))" %in% metricList) {
- metricList <- gsub("I\\(abs\\((log2AspectRatio)\\) \\* as\\.numeric\\((dist_border) < 50\\)\\)",
- "log2AspectRatio", metricList)
- }
+ formula_info <- NULL
+ if (is.null(modelFormula)) {
+ metric_list <- intersect(
+ .qscoreSupportedPredictors(),
+ names(colData(spe))
+ )
+ } else {
+ formula_info <- .validateQScoreFormula(
+ modelFormula=modelFormula,
+ dataNames=names(colData(spe)),
+ technology=metadata(spe)$technology
+ )
+ metric_list <- formula_info$training_metrics
+ }
+
+ if (!"log2SignalDensity" %in% metric_list) {
+ stop(
+ "'log2SignalDensity' is required to construct Quality Score ",
+ "training labels."
+ )
}
- stopifnot("Not all required metrics in the colData.\nPlease run spatialPerCellQC first." = all(metricList %in% names(colData(spe))))
- ctx <- .prepQCContext(spe, metricList, verbose)
+ ctx <- .prepQCContext(spe, metric_list, verbose)
df <- ctx$df; out_var <- ctx$out_var; tech <- ctx$tech
- model_formula <- getModelFormula(names(out_var))
+ if (is.null(modelFormula)) {
+ model_formula <- getModelFormula(names(out_var))
+ model_formula_object <- stats::as.formula(model_formula)
+ } else {
+ model_formula <- formula_info$text
+ model_formula_object <- formula_info$formula
+ }
if (verbose) {
message("Using final model formula:")
@@ -499,21 +544,27 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE
train_ok <- .filterCompleteModelCases(
df=train_df,
- modelFormula=model_formula,
+ modelFormula=model_formula_object,
response=train_df$QScore_train,
context="training cells"
)
train_df <- train_df[train_ok, , drop=FALSE]
- model_matrix <- model.matrix(as.formula(model_formula), data=train_df)
+ model_matrix <- stats::model.matrix(
+ model_formula_object,
+ data=train_df
+ )
model <- trainModel(model_matrix, train_df)
if(is.null(bestLambda)) {
- bestLambda <- computeLambda(model_matrix, train_df)
+ bestLambda <- .computeLambda(
+ modelMatrix=model_matrix,
+ response=train_df$QScore_train
+ )
}
- coefs <- coef(model, s=bestLambda)[coef(model, s=bestLambda)[,1]!= 0,,drop=FALSE]
+ coefs <- coef(model, s=bestLambda)
if (verbose) {
message("Model coefficients for every term used in the formula:")
@@ -521,11 +572,15 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE
collapse=" "))
}
- full_ok <- .filterCompleteModelCases(df=df, modelFormula=model_formula,
+ full_ok <- .filterCompleteModelCases(
+ df=df,
+ modelFormula=model_formula_object,
context="cells")
- full_matrix <- model.matrix(as.formula(model_formula),
- data=df[full_ok, , drop=FALSE])
+ full_matrix <- stats::model.matrix(
+ model_formula_object,
+ data=df[full_ok, , drop=FALSE]
+ )
## NAs may come from model variables used in `model_formula`, e.g. cells with
## missing `log2AspectRatio` when aspect ratio cannot be computed from polygons.
@@ -563,8 +618,8 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE
#' using \pkg{glmnet}, given a design matrix and a training data frame.
#'
#' @param trainDF `data.frame`
-#' A data frame containing at least the response column
-#' `qscore_train`, coded as 0/1.
+#' A data frame containing `QScore_train` or historical `qcscore_train`,
+#' coded as 0/1. If both are present, they must agree.
#' @param modelMatrix a matrix describing the model variables, tipically created
#' with `getModelFormula` and `model.matrix` functions.
#'
@@ -582,7 +637,8 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE
#' coef(fit, s = 0.01)
trainModel <- function(modelMatrix, trainDF)
{
- model <- glmnet(x=modelMatrix, y=trainDF$QScore_train,
+ response <- .getQScoreResponse(trainDF)
+ model <- glmnet(x=modelMatrix, y=response,
family="binomial", lambda=NULL, alpha=0)
return(model)
}
@@ -616,7 +672,8 @@ trainModel <- function(modelMatrix, trainDF)
#'
#' @return
#' A \code{data.frame} with one row per cell, including:
-#' \code{QScore_train} (0/1) indicating “bad” vs “good”,
+#' \code{QScore_train} and its compatibility alias
+#' \code{qcscore_train} (0/1) indicating “bad” vs “good”,
#' relevant \code{colData} columns used for modeling.
#' Deduplicates and down-samples “good” cells to match the number of “bad” cells.
#'
@@ -745,6 +802,9 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) {
train_df <- rbind(train_bad, train_good)
train_df <- train_df |>
dplyr::distinct(cell_id, .keep_all=TRUE)
+ ## `qcscore_train` remains available because `computeTrainDF()` is a
+ ## published API. Canonical internals use `QScore_train`.
+ train_df$qcscore_train <- train_df$QScore_train
return(train_df)
}
@@ -752,31 +812,183 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) {
#' @name getModelFormula
#' @rdname getModelFormula
#' @description
-#' Returns the right‐hand side of a model formula string based on formula
-#' variables found in the `metadata` of a `SpatialExperiment` object.
-#' @param formulaVars A named character vector mapping variable names
-#' (e.g. `"log2SignalDensity"`, `"Area_um"`, etc.) to their corresponding
-#' outlier label columns, typically from
-#' `metadata(spe)$formula_variables`.
+#' Returns the right‐hand side of a model formula string based on a vector of
+#' metric names.
+#' @param formulaVars A named character vector mapping metric names to outlier
+#' columns, as historically stored in
+#' `metadata(spe)$formula_variables`. An unnamed character vector of metric
+#' names is also accepted.
+#' @param verbose Logical. If `TRUE`, prints the generated formula.
+#' @param ... May contain the named replacement `metricList`. If both
+#' `formulaVars` and `metricList` are supplied, `metricList` takes precedence
+#' with a warning.
#' @return `character`
#' A one‐sided formula as a string (e.g. "~ log2SignalDensity + ...").
#' @export
#' @examples
#' example(checkOutliers)
#' getModelFormula(metadata(spe)$formula_variables)
-getModelFormula <- function(metricList)
+getModelFormula <- function(formulaVars, verbose=FALSE, ...)
{
- out_var <- metricList
+ has_formula_vars <- !missing(formulaVars)
+ dots <- list(...)
+ unknown_dots <- setdiff(names(dots), "metricList")
+ if (length(unknown_dots) > 0L) {
+ stop(
+ "Unused argument(s): ",
+ paste(unknown_dots, collapse=", "),
+ "."
+ )
+ }
+ has_metric_list <- "metricList" %in% names(dots)
+ if (!has_formula_vars && !has_metric_list) {
+ stop("'formulaVars' or 'metricList' must be supplied.")
+ }
+ if (has_formula_vars && has_metric_list) {
+ warning(
+ "Both 'formulaVars' and 'metricList' were supplied; ",
+ "using 'metricList'."
+ )
+ }
+ out_var <- if (has_metric_list) dots$metricList else formulaVars
+ if (!is.character(out_var)) {
+ stop("Model formula variables must be supplied as a character vector.")
+ }
+ if (!is.null(names(out_var)) && all(nzchar(names(out_var)))) {
+ out_var <- names(out_var)
+ }
if ("log2AspectRatio" %in% out_var) {
out_var[grep("log2AspectRatio", out_var)] <-
- "I(abs(log2AspectRatio) * as.numeric(dist_border<50))"
+ "I(abs(log2AspectRatio) * as.numeric(dist_border < 50))"
}
model_formula <- paste0("~(", paste(out_var, collapse = " + "),
")^2", sep = "")
+ if (verbose) {
+ message("Final formula used for Quality Score computation:")
+ message(model_formula)
+ }
+
return(model_formula)
}
+.qscoreSupportedPredictors <- function() {
+ c(
+ "log2SignalDensity",
+ "Area_um",
+ "log2AspectRatio",
+ "log2Ctrl_total_ratio"
+ )
+}
+
+.qscoreBorderTerm <- function() {
+ paste0(
+ "I(abs(log2AspectRatio)*",
+ "as.numeric(dist_border<50))"
+ )
+}
+
+.validateQScoreFormula <- function(modelFormula, dataNames, technology) {
+ model_formula <- tryCatch(
+ stats::as.formula(modelFormula),
+ error=function(e) {
+ stop("Unable to parse 'modelFormula': ", conditionMessage(e))
+ }
+ )
+ if (length(model_formula) != 2L) {
+ stop("'modelFormula' must be a one-sided formula.")
+ }
+
+ required_variables <- base::all.vars(model_formula)
+ supported <- .qscoreSupportedPredictors()
+ allowed_variables <- c(supported, "dist_border")
+ unsupported_variables <- setdiff(
+ required_variables,
+ allowed_variables
+ )
+ if (length(unsupported_variables) > 0L) {
+ stop(
+ "Unsupported Quality Score predictor(s): ",
+ paste(unsupported_variables, collapse=", "),
+ ". Formula structure may be customised, but predictors and ",
+ "transformations are restricted to: ",
+ paste(supported, collapse=", "),
+ "."
+ )
+ }
+
+ term_labels <- attr(stats::terms(model_formula), "term.labels")
+ if (length(term_labels) == 0L) {
+ stop("'modelFormula' must contain at least one supported predictor.")
+ }
+ normalise <- function(x) gsub("[[:space:]]+", "", x)
+ border_term <- .qscoreBorderTerm()
+ unsupported_terms <- character()
+ for (term_label in term_labels) {
+ components <- strsplit(
+ normalise(term_label),
+ ":",
+ fixed=TRUE
+ )[[1]]
+ valid_components <- components %in% supported |
+ components == border_term
+ if (!all(valid_components)) {
+ unsupported_terms <- c(unsupported_terms, term_label)
+ }
+ }
+ unsupported_terms <- unique(unsupported_terms)
+ if (length(unsupported_terms) > 0L) {
+ stop(
+ "Unsupported Quality Score formula term(s) or transformation(s): ",
+ paste(unsupported_terms, collapse=", "),
+ ". Formula structure and interactions may be customised, but ",
+ "predictors and transformations may not be extended. Supported ",
+ "predictors are: ",
+ paste(supported, collapse=", "),
+ "; the supported border-effect expression is ",
+ "'I(abs(log2AspectRatio) * as.numeric(dist_border < 50))'."
+ )
+ }
+
+ missing_variables <- setdiff(required_variables, dataNames)
+ if (length(missing_variables) > 0L) {
+ stop(
+ "Missing variables required by 'modelFormula': ",
+ paste(missing_variables, collapse=", "),
+ ". Run 'spatialPerCellQC()' before computing Quality Score."
+ )
+ }
+
+ is_cosmx <- technology %in%
+ c("Nanostring_CosMx", "Nanostring_CosMx_Protein")
+ uses_border_predictor <- any(
+ c("log2AspectRatio", "dist_border") %in% required_variables
+ )
+ if (uses_border_predictor && !is_cosmx) {
+ stop(
+ "'log2AspectRatio' and the border-effect expression are ",
+ "supported only for Nanostring CosMx datasets; detected ",
+ "technology: '", technology, "'."
+ )
+ }
+
+ formula_text <- if (is.character(modelFormula)) {
+ paste(modelFormula, collapse=" ")
+ } else {
+ paste(deparse(model_formula, width.cutoff=500L), collapse=" ")
+ }
+ training_metrics <- required_variables[
+ required_variables %in% supported
+ ]
+ list(
+ formula=model_formula,
+ text=formula_text,
+ variables=required_variables,
+ terms=term_labels,
+ training_metrics=training_metrics
+ )
+}
+
#' .computeXenMerTrainSet
#' @name dot-computeXenMerTrainSet
#' @rdname dot-computeXenMerTrainSet
@@ -1002,6 +1214,19 @@ computeOutliersQScore <- function(spe, metricList=c("log2SignalDensity","Area_um
"log2AspectRatio", "log2Ctrl_total_ratio")) {
stopifnot(is(spe, "SpatialExperiment"))
+ unsupported_metrics <- setdiff(
+ metricList,
+ .qscoreSupportedPredictors()
+ )
+ if (length(unsupported_metrics) > 0L) {
+ stop(
+ "Unsupported Quality Score metric(s): ",
+ paste(unsupported_metrics, collapse=", "),
+ ". Supported metrics are: ",
+ paste(.qscoreSupportedPredictors(), collapse=", "),
+ "."
+ )
+ }
cd <- colData(spe)
method <- .checkSkw(cd, metricList)
@@ -1118,8 +1343,9 @@ checkOutliers <- function(spe, verbose=FALSE) {
if (verbose) {
for (i in names(out_var)) {
message("Outliers found for ", i, ":")
- for (k in 1:length(names(table(cd[[out_var[i]]])))){
- message(paste0(names(table(cd[[out_var[i]]]))[k], ": ", table(cd[[out_var[i]]])[k]))
+ outlier_counts <- table(cd[[out_var[i]]])
+ for (k in seq_along(outlier_counts)) {
+ message(names(outlier_counts)[k], ": ", outlier_counts[k])
}
}
}
@@ -1190,8 +1416,8 @@ checkOutliers <- function(spe, verbose=FALSE) {
# remove zero-count cells once
zerocells <- spe$total==0
if (sum(zerocells) > 0) {
- warning(paste0(sum(zerocells),
- " cells with 0 counts were found. These cells will be removed."))
+ warning(sum(zerocells),
+ " cells with 0 counts were found. These cells will be removed.")
spe <- spe[, !zerocells]
}
if("log2CountArea" %in% names(colData(spe)))
@@ -1209,12 +1435,13 @@ checkOutliers <- function(spe, verbose=FALSE) {
return(list(df=df, out_var=out_var, tech=tech))
}
-#' .applyQScoreModel
-#' @name .applyQScoreModel
-#' @rdname dot-applyQScoreModel
+#' applyQScoreModel
+#' @name applyQScoreModel
+#' @rdname applyQScoreModel
#'
#' @description
-#' Internal: applies a previously trained Quality Score model to a new SpatialExperiment object.
+#' Applies a previously trained Quality Score model to a new SpatialExperiment
+#' object.
#' See details for important considerations when applying a model to a
#' different dataset.
#' @details
@@ -1255,16 +1482,20 @@ checkOutliers <- function(spe, verbose=FALSE) {
#'
#' ## Train the Quality Control (QC) score model on one dataset
#' spe_train <- computeQScore(spe_train)
-#' qc_model <- metadata(spe_train)$QScore_model
+#' qs_model <- metadata(spe_train)$QScore_model
#'
#' ## Apply the trained model to another dataset
#' spe_test <- applyQScoreModel(
#' spe=spe_test,
-#' qcModel=qc_model,
+#' qsModel=qs_model,
#' scoreName="QScore_transferred"
#' )
#'
#' summary(spe_test$QScore_transferred)
+applyQScoreModel <- function(spe, qsModel, scoreName="QScore") {
+ .applyQScoreModel(spe=spe, qsModel=qsModel, scoreName=scoreName)
+}
+
.applyQScoreModel <- function(spe, qsModel, scoreName="QScore") {
stopifnot(is(spe, "SpatialExperiment"))
@@ -1281,6 +1512,11 @@ checkOutliers <- function(spe, verbose=FALSE) {
)
df <- as.data.frame(colData(spe))
+ .validateQScoreFormula(
+ modelFormula=qsModel$model_formula,
+ dataNames=colnames(df),
+ technology=metadata(spe)$technology
+ )
ok <- .filterCompleteModelCases(
df=df,
@@ -1353,11 +1589,35 @@ checkOutliers <- function(spe, verbose=FALSE) {
return(modelMatrix)
}
+.getQScoreResponse <- function(trainDF) {
+ has_canonical <- "QScore_train" %in% colnames(trainDF)
+ has_legacy <- "qcscore_train" %in% colnames(trainDF)
+ if (!has_canonical && !has_legacy) {
+ stop(
+ "Training data must contain 'QScore_train' or ",
+ "'qcscore_train'."
+ )
+ }
+ if (has_canonical && has_legacy) {
+ canonical <- trainDF$QScore_train
+ legacy <- trainDF$qcscore_train
+ same <- (is.na(canonical) & is.na(legacy)) |
+ (!is.na(canonical) & !is.na(legacy) & canonical == legacy)
+ if (any(!same)) {
+ stop(
+ "'QScore_train' and 'qcscore_train' are both present ",
+ "but contain different values."
+ )
+ }
+ }
+ if (has_canonical) trainDF$QScore_train else trainDF$qcscore_train
+}
+
#' @importFrom stats complete.cases
.filterCompleteModelCases <- function(df, modelFormula, response=NULL,
context="cells") {
- vars <- all.vars(as.formula(modelFormula))
+ vars <- base::all.vars(stats::as.formula(modelFormula))
missing_vars <- setdiff(vars, colnames(df))
@@ -1386,3 +1646,143 @@ checkOutliers <- function(spe, verbose=FALSE) {
return(ok)
}
+
+## ---- Deprecated functions -----------------------------------------------
+
+#' Deprecated functions in SpaceTrooper
+#'
+#' @name SpaceTrooper-deprecated
+#' @description
+#' These functions are retained for compatibility with published SpaceTrooper
+#' APIs. They issue a deprecation warning and direct users to the canonical
+#' Quality Score functions.
+#'
+#' @details
+#' \itemize{
+#' \item `computeQCScore()` replaces canonical `QScore` and
+#' `metadata(spe)$QScore_model` with historical `QC_score` and
+#' `metadata(spe)$QCScore_model`.
+#' \item `computeQCScoreFlags()` consumes `QC_score` and creates
+#' `low_qcscore` and, when applicable, `low_threshold_qcscore`.
+#' \item `computeOutliersQCScore()` is the historical name for
+#' `computeOutliersQScore()`.
+#' \item `applyQCScoreModel()` retains the `qcModel` argument, the default
+#' `QC_score` output, and `metadata(spe)$QCScore_model_applied`.
+#' }
+#'
+#' Use `computeQScore()`, `computeQScoreFlags()`,
+#' `computeOutliersQScore()`, and `applyQScoreModel()` for canonical names.
+NULL
+
+#' computeQCScore (deprecated)
+#' @rdname SpaceTrooper-deprecated
+#'
+#' @param spe A `SpatialExperiment` object.
+#' @param bestLambda Passed to \code{\link{computeQScore}}.
+#' @param modelFormula Passed to \code{\link{computeQScore}}.
+#' @param verbose Passed to \code{\link{computeQScore}}.
+#' @return A `SpatialExperiment` object; see \code{\link{computeQScore}}.
+#' @export
+computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL,
+ verbose=FALSE) {
+ .Deprecated(
+ new="computeQScore",
+ package="SpaceTrooper",
+ msg=paste0(
+ "'computeQCScore' is deprecated.\n",
+ "Use 'computeQScore' instead.\n",
+ "See help('computeQScore') for details."
+ )
+ )
+ spe <- computeQScore(
+ spe,
+ bestLambda=bestLambda,
+ modelFormula=modelFormula,
+ verbose=verbose
+ )
+ spe$QC_score <- spe$QScore
+ spe$QScore <- NULL
+ metadata(spe)$QCScore_model <- metadata(spe)$QScore_model
+ metadata(spe)$QScore_model <- NULL
+ return(spe)
+}
+
+#' computeQCScoreFlags (deprecated)
+#' @rdname SpaceTrooper-deprecated
+#'
+#' @param spe A `SpatialExperiment` object.
+#' @param qsThreshold Passed to \code{\link{computeQScoreFlags}}.
+#' @param useQSQuantiles Passed to \code{\link{computeQScoreFlags}}.
+#' @return A `SpatialExperiment` object; see \code{\link{computeQScoreFlags}}.
+#' @export
+computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) {
+ .Deprecated(
+ new="computeQScoreFlags",
+ package="SpaceTrooper",
+ msg=paste0(
+ "'computeQCScoreFlags' is deprecated.\n",
+ "Use 'computeQScoreFlags' instead.\n",
+ "See help('computeQScoreFlags') for details."
+ )
+ )
+ stopifnot(is(spe, "SpatialExperiment"))
+ stopifnot("QC_score" %in% names(colData(spe)))
+
+ if (useQSQuantiles) {
+ spe$low_qcscore <- spe$QC_score <
+ quantile(spe$QC_score, probs=qsThreshold, na.rm=TRUE)
+ } else {
+ spe$low_qcscore <- spe$QC_score < qsThreshold
+ }
+
+ if ("threshold_flags" %in% names(colData(spe))) {
+ spe$low_threshold_qcscore <- spe$low_qcscore &
+ spe$threshold_flags
+ }
+ return(spe)
+}
+
+#' computeOutliersQCScore (deprecated)
+#' @rdname SpaceTrooper-deprecated
+#'
+#' @param spe A `SpatialExperiment` object.
+#' @param metricList Passed to \code{\link{computeOutliersQScore}}.
+#' @return A `SpatialExperiment` object; see
+#' \code{\link{computeOutliersQScore}}.
+#' @export
+computeOutliersQCScore <- function(spe, metricList=c(
+ "log2SignalDensity", "Area_um", "log2AspectRatio",
+ "log2Ctrl_total_ratio"
+)) {
+ .Deprecated(
+ new="computeOutliersQScore",
+ package="SpaceTrooper"
+ )
+ computeOutliersQScore(spe=spe, metricList=metricList)
+}
+
+#' applyQCScoreModel (deprecated)
+#' @rdname SpaceTrooper-deprecated
+#'
+#' @param spe A `SpatialExperiment` object with QC metrics already computed.
+#' @param qcModel A historical QC score model object, usually stored in
+#' `metadata(spe)$QCScore_model`.
+#' @param scoreName Name of the legacy output column in `colData`.
+#' @return A `SpatialExperiment` object with the applied score in `colData`
+#' and the model in `metadata(spe)$QCScore_model_applied`.
+#' @export
+applyQCScoreModel <- function(spe, qcModel, scoreName="QC_score") {
+ .Deprecated(
+ new="applyQScoreModel",
+ package="SpaceTrooper"
+ )
+ spe <- .applyQScoreModel(
+ spe=spe,
+ qsModel=qcModel,
+ scoreName=scoreName
+ )
+ metadata(spe)$QCScore_model_applied <-
+ metadata(spe)$QScore_model_applied
+ metadata(spe)$QScore_model_applied <- NULL
+ return(spe)
+}
diff --git a/R/spatialQCPlots.R b/R/spatialQCPlots.R
index c730dac..711fb8b 100644
--- a/R/spatialQCPlots.R
+++ b/R/spatialQCPlots.R
@@ -10,17 +10,23 @@
#' @param sampleId Character string identifying which sample to plot.
#' Default: `unique(spe$sample_id)`.
#' @param pointCol Color for the cell centroids. Default: `"firebrick"`.
-#' @param pointSize Numeric point size for the cell centroids. Default: `0.05`.
-#' @param pointAlpha Numeric transparency for the cell centroids. Default: `0.8`.
#' @param numbersCol Color for the FoV labels. Default: `"black"`.
-#' @param numberSize Numeric size for the FoV labels. Default: `1`.
-#' @param numbersAlpha Numeric transparency for FoV labels. Default: `0.8`.
+#' @param alphaNumbers Deprecated alias for `numbersAlpha`.
#' @param fovDim numeric with two named dimensions xdim, ydim. (Default is
#' metadata(spe)$fov_dim)
+#' @param size Deprecated alias for `pointSize`.
+#' @param alpha Deprecated alias for `pointAlpha`.
#' @param scaleBar A logical value indicating whether to add a scale bar to the
#' plot. (Default is `TRUE`)
#' @param micronConvFact Numeric conversion factor from pixels to microns.
#' Default is `0.12`.
+#' @param pointSize Numeric point size for the cell centroids. If supplied,
+#' takes precedence over `size`.
+#' @param pointAlpha Numeric transparency for the cell centroids. If supplied,
+#' takes precedence over `alpha`.
+#' @param numberSize Numeric size for the FoV labels. Default: `1`.
+#' @param numbersAlpha Numeric transparency for FoV labels. If supplied, takes
+#' precedence over `alphaNumbers`.
#'
#' @return A `ggplot` object showing cell centroids and FoV boundaries.
#'
@@ -62,16 +68,35 @@
#' g <- plotCellsFovs(spe)
#' print(g)
plotCellsFovs <- function(spe, sampleId=unique(spe$sample_id),
- pointCol="firebrick", pointSize=0.05,
- pointAlpha=0.8, numbersCol="black",
- numberSize= 1, numbersAlpha=0.8,
- fovDim=metadata(spe)$fov_dim,
- scaleBar=TRUE, micronConvFact = 0.12)
+ pointCol="firebrick", numbersCol="black",
+ alphaNumbers=0.8, fovDim=metadata(spe)$fov_dim,
+ size=0.05, alpha=0.8,
+ scaleBar=TRUE, micronConvFact=0.12,
+ pointSize=NULL, pointAlpha=NULL, numberSize=1,
+ numbersAlpha=NULL)
{
stopifnot(is(spe, "SpatialExperiment"))
stopifnot("fov" %in% names(colData(spe)))
stopifnot( all(names(fovDim) %in% c("xdim","ydim")) )
+ if (!is.null(pointSize) && !missing(size)) {
+ warning("Both 'size' and 'pointSize' were supplied; using 'pointSize'.")
+ }
+ if (!is.null(pointAlpha) && !missing(alpha)) {
+ warning(
+ "Both 'alpha' and 'pointAlpha' were supplied; using 'pointAlpha'."
+ )
+ }
+ if (!is.null(numbersAlpha) && !missing(alphaNumbers)) {
+ warning(
+ "Both 'alphaNumbers' and 'numbersAlpha' were supplied; ",
+ "using 'numbersAlpha'."
+ )
+ }
+ if (is.null(pointSize)) pointSize <- size
+ if (is.null(pointAlpha)) pointAlpha <- alpha
+ if (is.null(numbersAlpha)) numbersAlpha <- alphaNumbers
+
spd <- as.data.frame(spatialCoords(spe))
x_coord <- spatialCoordsNames(spe)[1]
y_coord <- spatialCoordsNames(spe)[2]
@@ -439,12 +464,10 @@ plotPolygons <- function(spe, colourBy="darkgrey", colourLog=FALSE,
#' plot. If `NULL`, no title is added. Default is `NULL`.
#' @param mapPointCol A character string specifying the color of the points
#' in the map. Default is `"darkmagenta"`.
-#' @param mapPointSize Numeric size for points in the map. Default: `0.5`.
-#' @param mapPointAlpha Numeric transparency for points in the map. Default: `0.8`.
-#' @param fovNumbersCol A character string specifying the color of the
-#' numbers on the FoV zoom-in. Default is `"black"`.
-#' @param fovNumberSize Numeric size for the FoV labels. Default: `1`.
-#' @param fovNumbersAlpha Numeric transparency for FoV labels. Default: `0.8`.
+#' @param mapNumbersCol Deprecated alias for `fovNumbersCol`.
+#' @param mapAlphaNumbers Deprecated alias for `fovNumbersAlpha`.
+#' @param csize Deprecated alias for `mapPointSize`.
+#' @param calpha Deprecated alias for `mapPointAlpha`.
#' @param scaleBars Logical or NULL. Default is `NULL`.
#' Master switch controlling the presence of scale bars in both panels.
#' If \code{TRUE}, scale bars are shown in both the map and polygon panels.
@@ -457,6 +480,15 @@ plotPolygons <- function(spe, colourBy="darkgrey", colourLog=FALSE,
#' These parameters are only used when \code{scaleBars} is \code{NULL};
#' otherwise they are overridden by \code{scaleBars}.
#' @param ... Additional arguments passed to `plotPolygons`.
+#' @param mapPointSize Numeric size for points in the map. If supplied, takes
+#' precedence over `csize`.
+#' @param mapPointAlpha Numeric transparency for points in the map. If supplied,
+#' takes precedence over `calpha`.
+#' @param fovNumbersCol Color for FoV labels. If supplied, takes precedence over
+#' `mapNumbersCol`.
+#' @param fovNumberSize Numeric size for the FoV labels. Default: `1`.
+#' @param fovNumbersAlpha Transparency for FoV labels. If supplied, takes
+#' precedence over `mapAlphaNumbers`.
#'
#' @return A combined plot showing a map of all FOVs with zoomed-in views of
#' the specified FOVs and their associated polygons.
@@ -475,17 +507,48 @@ plotPolygons <- function(spe, colourBy="darkgrey", colourLog=FALSE,
#' plotZoomFovsMap(spe, fovs=16, title="FOV 16")
plotZoomFovsMap <- function(spe, fovs=NULL, title=NULL,
mapPointCol="darkmagenta",
- mapPointSize=0.5, mapPointAlpha=0.8,
- fovNumbersCol="black", fovNumberSize=1,
- fovNumbersAlpha=0.8,
+ mapNumbersCol="black",
+ mapAlphaNumbers=0.8,
csize=0.05, calpha=0.8,
scaleBars=NULL,
scaleBarMap=TRUE,
scaleBarPol=TRUE,
- ...) {
+ ...,
+ mapPointSize=NULL, mapPointAlpha=NULL,
+ fovNumbersCol=NULL, fovNumberSize=1,
+ fovNumbersAlpha=NULL) {
stopifnot(is(spe, "SpatialExperiment"))
stopifnot("fov" %in% names(colData(spe)))
stopifnot(all(fovs %in% spe$fov))
+ if (!is.null(mapPointSize) && !missing(csize)) {
+ warning(
+ "Both 'csize' and 'mapPointSize' were supplied; ",
+ "using 'mapPointSize'."
+ )
+ }
+ if (!is.null(mapPointAlpha) && !missing(calpha)) {
+ warning(
+ "Both 'calpha' and 'mapPointAlpha' were supplied; ",
+ "using 'mapPointAlpha'."
+ )
+ }
+ if (!is.null(fovNumbersCol) && !missing(mapNumbersCol)) {
+ warning(
+ "Both 'mapNumbersCol' and 'fovNumbersCol' were supplied; ",
+ "using 'fovNumbersCol'."
+ )
+ }
+ if (!is.null(fovNumbersAlpha) && !missing(mapAlphaNumbers)) {
+ warning(
+ "Both 'mapAlphaNumbers' and 'fovNumbersAlpha' were supplied; ",
+ "using 'fovNumbersAlpha'."
+ )
+ }
+ if (is.null(mapPointSize)) mapPointSize <- csize
+ if (is.null(mapPointAlpha)) mapPointAlpha <- calpha
+ if (is.null(fovNumbersCol)) fovNumbersCol <- mapNumbersCol
+ if (is.null(fovNumbersAlpha)) fovNumbersAlpha <- mapAlphaNumbers
+
spefovs <- spe[, spe$fov %in% fovs]
if (!is.null(scaleBars)) {
stopifnot(is.logical(scaleBars), length(scaleBars) == 1L)
diff --git a/README.md b/README.md
index fd21d22..11d3d1f 100644
--- a/README.md
+++ b/README.md
@@ -122,13 +122,24 @@ spe <- computeQScoreFlags(spe, qsThreshold=0.5)
# 5. Visualization
## Visualize cells as dots in their centroid coordinates, colored by a column in `colData(spe)` (e.g., QS computed above).
-plotCentroids(spe, colourBy='QC_score')
+plotCentroids(spe, colourBy="QScore")
## Visualize cells using their polygon shapes, colored by a column in `colData(spe)` (e.g., QS computed above).
## To visualize polygons step 2 is mandatory.
## Polygons can be cumbersome to plot for large datasets (e.g., entire slides with more than 100,000 cells), hence centroids may be preferred.
-plotPolygons(spe, colourBy='QC_score')
+plotPolygons(spe, colourBy="QScore")
```
+
+Custom `modelFormula` values may select any subset of the supported predictors
+`log2SignalDensity`, `Area_um`, `log2AspectRatio`, and
+`log2Ctrl_total_ratio`. Additive formulas and selected interactions are
+preserved exactly; unsupported predictors and transformations are rejected.
+The border-effect expression is available only for CosMx data:
+`I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`.
+
+The former `computeQCScore()`, `computeQCScoreFlags()`,
+`computeOutliersQCScore()`, and `applyQCScoreModel()` APIs remain available as
+deprecated compatibility functions and retain their historical output names.
@@ -156,7 +167,8 @@ Refer to the first vignette if working with CosMx, Xenium or MERFISH spatial tra
Function-level help is available through standard R documentation, for example:
```r
-?computeQCScore
+?computeQScore
+?SpaceTrooper-deprecated
?plotPolygons
```
diff --git a/man/SpaceTrooper-deprecated.Rd b/man/SpaceTrooper-deprecated.Rd
new file mode 100644
index 0000000..f07cf57
--- /dev/null
+++ b/man/SpaceTrooper-deprecated.Rd
@@ -0,0 +1,74 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/QC.R
+\name{SpaceTrooper-deprecated}
+\alias{SpaceTrooper-deprecated}
+\alias{computeQCScore}
+\alias{computeQCScoreFlags}
+\alias{computeOutliersQCScore}
+\alias{applyQCScoreModel}
+\title{Deprecated functions in SpaceTrooper}
+\usage{
+computeQCScore(spe, bestLambda = NULL, modelFormula = NULL, verbose = FALSE)
+
+computeQCScoreFlags(spe, qsThreshold = 0.5, useQSQuantiles = FALSE)
+
+computeOutliersQCScore(
+ spe,
+ metricList = c("log2SignalDensity", "Area_um", "log2AspectRatio",
+ "log2Ctrl_total_ratio")
+)
+
+applyQCScoreModel(spe, qcModel, scoreName = "QC_score")
+}
+\arguments{
+\item{spe}{A `SpatialExperiment` object with QC metrics already computed.}
+
+\item{bestLambda}{Passed to \code{\link{computeQScore}}.}
+
+\item{modelFormula}{Passed to \code{\link{computeQScore}}.}
+
+\item{verbose}{Passed to \code{\link{computeQScore}}.}
+
+\item{qsThreshold}{Passed to \code{\link{computeQScoreFlags}}.}
+
+\item{useQSQuantiles}{Passed to \code{\link{computeQScoreFlags}}.}
+
+\item{metricList}{Passed to \code{\link{computeOutliersQScore}}.}
+
+\item{qcModel}{A historical QC score model object, usually stored in
+`metadata(spe)$QCScore_model`.}
+
+\item{scoreName}{Name of the legacy output column in `colData`.}
+}
+\value{
+A `SpatialExperiment` object; see \code{\link{computeQScore}}.
+
+A `SpatialExperiment` object; see \code{\link{computeQScoreFlags}}.
+
+A `SpatialExperiment` object; see
+ \code{\link{computeOutliersQScore}}.
+
+A `SpatialExperiment` object with the applied score in `colData`
+ and the model in `metadata(spe)$QCScore_model_applied`.
+}
+\description{
+These functions are retained for compatibility with published SpaceTrooper
+APIs. They issue a deprecation warning and direct users to the canonical
+Quality Score functions.
+}
+\details{
+\itemize{
+ \item `computeQCScore()` replaces canonical `QScore` and
+ `metadata(spe)$QScore_model` with historical `QC_score` and
+ `metadata(spe)$QCScore_model`.
+ \item `computeQCScoreFlags()` consumes `QC_score` and creates
+ `low_qcscore` and, when applicable, `low_threshold_qcscore`.
+ \item `computeOutliersQCScore()` is the historical name for
+ `computeOutliersQScore()`.
+ \item `applyQCScoreModel()` retains the `qcModel` argument, the default
+ `QC_score` output, and `metadata(spe)$QCScore_model_applied`.
+}
+
+Use `computeQScore()`, `computeQScoreFlags()`,
+`computeOutliersQScore()`, and `applyQScoreModel()` for canonical names.
+}
diff --git a/man/dot-applyQScoreModel.Rd b/man/applyQScoreModel.Rd
similarity index 87%
rename from man/dot-applyQScoreModel.Rd
rename to man/applyQScoreModel.Rd
index b126f31..f26f7b0 100644
--- a/man/dot-applyQScoreModel.Rd
+++ b/man/applyQScoreModel.Rd
@@ -1,10 +1,10 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/QC.R
-\name{.applyQScoreModel}
-\alias{.applyQScoreModel}
-\title{.applyQScoreModel}
+\name{applyQScoreModel}
+\alias{applyQScoreModel}
+\title{applyQScoreModel}
\usage{
-.applyQScoreModel(spe, qsModel, scoreName = "QScore")
+applyQScoreModel(spe, qsModel, scoreName = "QScore")
}
\arguments{
\item{spe}{A `SpatialExperiment` object with QC metrics already computed.}
@@ -18,7 +18,8 @@
A `SpatialExperiment` object with added Quality Score in `colData`.
}
\description{
-Internal: applies a previously trained Quality Score model to a new SpatialExperiment object.
+Applies a previously trained Quality Score model to a new SpatialExperiment
+object.
See details for important considerations when applying a model to a
different dataset.
}
@@ -52,12 +53,12 @@ spe_test <- spe[, -idx]
## Train the Quality Control (QC) score model on one dataset
spe_train <- computeQScore(spe_train)
-qc_model <- metadata(spe_train)$QScore_model
+qs_model <- metadata(spe_train)$QScore_model
## Apply the trained model to another dataset
spe_test <- applyQScoreModel(
spe=spe_test,
- qcModel=qc_model,
+ qsModel=qs_model,
scoreName="QScore_transferred"
)
diff --git a/man/computeLambda.Rd b/man/computeLambda.Rd
index e7becff..b79ee4c 100644
--- a/man/computeLambda.Rd
+++ b/man/computeLambda.Rd
@@ -4,17 +4,17 @@
\alias{computeLambda}
\title{computeLambda}
\usage{
-computeLambda(modelMatrix, trainDF)
+computeLambda(trainDF, modelFormula)
}
\arguments{
\item{trainDF}{`data.frame`
A data frame for QS model training that must include:
Predictor columns: All columns referenced in the formula returned by `getModelFormula()`.
- `qscore_train` A binary (0/1) response vector to be modeled.}
+ `QScore_train` or historical `qcscore_train`, containing the binary
+ (0/1) response to be modeled.}
-\item{modelFormula}{`character`
-A character string representing the model formula
- `~ log2SignalDensity + ...`, as returned by `getModelFormula()`.}
+\item{modelFormula}{A one-sided formula or character string describing the
+supported predictor terms.}
}
\value{
`numeric`
@@ -30,14 +30,14 @@ k-fold cross-validation to identify the optimal regularization parameter
\eqn{\lambda} for Quality Score (QS) model training.
}
\details{
-Internally, the function:
- constructs the design matrix via \code{model.matrix()},
- runs k-fold cross-validation of ridge logistic regression using `cv.glmnet` with `alpha = 0`,
- extracts and returns `ridge_cv$lambda.min`.
+Complete cases for the formula variables and response are retained. The
+function constructs a model matrix and passes it to an internal
+matrix-based implementation that runs ridge logistic-regression
+cross-validation using `cv.glmnet` with `alpha = 0`.
}
\examples{
example(computeTrainDF)
-modform <- getModelFormula(metadata(spe)$formula_variables)
+modform <- getModelFormula(names(metadata(spe)$formula_variables))
best_lambda <- computeLambda(df_train, modform)
print(best_lambda)
diff --git a/man/computeQScore.Rd b/man/computeQScore.Rd
index 730e87b..8b5ceac 100644
--- a/man/computeQScore.Rd
+++ b/man/computeQScore.Rd
@@ -11,10 +11,10 @@ computeQScore(spe, bestLambda = NULL, modelFormula = NULL, verbose = FALSE)
\item{bestLambda}{the best lambda typically computed using `computeLambda`.}
-\item{modelFormula}{a character string representing the formula to be used for
-training the model. If NULL, the formula is automatically generated
-based on the available metrics and their outliers in the dataset. See details
-for more information.}
+\item{modelFormula}{A one-sided formula or character string. If `NULL`, the
+default formula is generated from available supported metrics and their
+outliers. If supplied, its supported terms and interactions are fitted
+without rewriting. See Details.}
\item{verbose}{logical for having a verbose output. Default is FALSE.}
}
@@ -49,9 +49,9 @@ Note that the function is responsible for automatically including/excluding
metrics in the formula based on their availability in the `colData` of the
`SpatialExperiment` object.
-Inclusion of metrics in the formula depends also on the number of available
-outliers. If the number of outliers for each metric is less than 0.1\% out of the
-entire dataset, the metric will be excluded from the QS formula.
+For the default formula, inclusion of metrics also depends on the number of
+available outliers. If the number of outliers for a metric is less than
+0.1\% of the dataset, that metric is excluded from the default QS formula.
- Model fitting: ridge (L2) logistic regression is fitted (via `glmnet`) on
the balanced training set. The function uses `trainModel()` for fitting
@@ -63,19 +63,28 @@ be computed internally, just set a seed with `set.seed()` before running
`computeQScore`. Otherwise, a fixed value of lambda previously computed with
`computeLambda` preceeded by `computeTrainDF` and `getModelFormula` can be set.
-- Model formula details: the model formula is automatically generated as
-follows:
+- Model formula details: the most complete default formula is:
`~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2`.
- When user-provided, the formula must follow the same default syntax and
-removed (or added) terms should be written exactly as in the default formula,
-e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))` must have spaces
-around the `*` and `<` operators.
-In any case, metrics with insufficient outliers (less than 0.1\% of the dataset)
-will be excluded from the QS formula.
+A user-supplied formula is fitted as supplied: additive structures, selected
+interactions, `:` and `*` are preserved. The formula must include
+`log2SignalDensity` and may use any subset of the other supported predictors:
+`Area_um`,
+`log2AspectRatio`, and `log2Ctrl_total_ratio`. Arbitrary predictors and
+transformations are rejected. The supported CosMx-only border expression is
+`I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`.
+
+When a custom `modelFormula` is supplied, the supported base QC metrics in
+that formula are also used to compute outliers, define good and bad training
+cells, and construct `QScore_train`. Interactions do not create additional
+training metrics. In the supported border expression, `dist_border` is an
+auxiliary variable and `log2AspectRatio` is the corresponding training QC
+metric. Supported metrics omitted from the custom formula do not contribute
+to training-label construction. When `modelFormula=NULL`, the established
+default metric-selection and training-label workflow is unchanged.
The computed model output is stored in `metadata(spe)$QScore_model`.
For inspection and model coefficient transfer between datasets
-see also \code{\link{.applyQScoreModel}}. However, transferring models between
+see also \code{\link{applyQScoreModel}}. However, transferring models between
datasets is not recommended. As described in the paper, QS model training is
dataset-specific and does not generalize well across datasets. Moreover,
this step is computationally efficient. Please, refer to the paper for
diff --git a/man/computeTrainDF.Rd b/man/computeTrainDF.Rd
index 0a54223..3d0666b 100644
--- a/man/computeTrainDF.Rd
+++ b/man/computeTrainDF.Rd
@@ -27,7 +27,8 @@ If \code{TRUE}, prints the number of “bad” and “good” cells selected.}
}
\value{
A \code{data.frame} with one row per cell, including:
- \code{QScore_train} (0/1) indicating “bad” vs “good”,
+ \code{QScore_train} and its compatibility alias
+ \code{qcscore_train} (0/1) indicating “bad” vs “good”,
relevant \code{colData} columns used for modeling.
Deduplicates and down-samples “good” cells to match the number of “bad” cells.
}
diff --git a/man/getModelFormula.Rd b/man/getModelFormula.Rd
index 49e67fd..f312fae 100644
--- a/man/getModelFormula.Rd
+++ b/man/getModelFormula.Rd
@@ -4,21 +4,27 @@
\alias{getModelFormula}
\title{getModelFormula}
\usage{
-getModelFormula(metricList)
+getModelFormula(formulaVars, verbose = FALSE, ...)
}
\arguments{
-\item{formulaVars}{A named character vector mapping variable names
-(e.g. `"log2SignalDensity"`, `"Area_um"`, etc.) to their corresponding
-outlier label columns, typically from
-`metadata(spe)$formula_variables`.}
+\item{formulaVars}{A named character vector mapping metric names to outlier
+columns, as historically stored in
+`metadata(spe)$formula_variables`. An unnamed character vector of metric
+names is also accepted.}
+
+\item{verbose}{Logical. If `TRUE`, prints the generated formula.}
+
+\item{...}{May contain the named replacement `metricList`. If both
+`formulaVars` and `metricList` are supplied, `metricList` takes precedence
+with a warning.}
}
\value{
`character`
A one‐sided formula as a string (e.g. "~ log2SignalDensity + ...").
}
\description{
-Returns the right‐hand side of a model formula string based on formula
-variables found in the `metadata` of a `SpatialExperiment` object.
+Returns the right‐hand side of a model formula string based on a vector of
+metric names.
}
\examples{
example(checkOutliers)
diff --git a/man/plotCellsFovs.Rd b/man/plotCellsFovs.Rd
index d0377a7..e6972d5 100644
--- a/man/plotCellsFovs.Rd
+++ b/man/plotCellsFovs.Rd
@@ -8,14 +8,17 @@ plotCellsFovs(
spe,
sampleId = unique(spe$sample_id),
pointCol = "firebrick",
- pointSize = 0.05,
- pointAlpha = 0.8,
numbersCol = "black",
- numberSize = 1,
- numbersAlpha = 0.8,
+ alphaNumbers = 0.8,
fovDim = metadata(spe)$fov_dim,
+ size = 0.05,
+ alpha = 0.8,
scaleBar = TRUE,
- micronConvFact = 0.12
+ micronConvFact = 0.12,
+ pointSize = NULL,
+ pointAlpha = NULL,
+ numberSize = 1,
+ numbersAlpha = NULL
)
}
\arguments{
@@ -26,24 +29,33 @@ Default: `unique(spe$sample_id)`.}
\item{pointCol}{Color for the cell centroids. Default: `"firebrick"`.}
-\item{pointSize}{Numeric point size for the cell centroids. Default: `0.05`.}
-
-\item{pointAlpha}{Numeric transparency for the cell centroids. Default: `0.8`.}
-
\item{numbersCol}{Color for the FoV labels. Default: `"black"`.}
-\item{numberSize}{Numeric size for the FoV labels. Default: `1`.}
-
-\item{numbersAlpha}{Numeric transparency for FoV labels. Default: `0.8`.}
+\item{alphaNumbers}{Deprecated alias for `numbersAlpha`.}
\item{fovDim}{numeric with two named dimensions xdim, ydim. (Default is
metadata(spe)$fov_dim)}
+\item{size}{Deprecated alias for `pointSize`.}
+
+\item{alpha}{Deprecated alias for `pointAlpha`.}
+
\item{scaleBar}{A logical value indicating whether to add a scale bar to the
plot. (Default is `TRUE`)}
\item{micronConvFact}{Numeric conversion factor from pixels to microns.
Default is `0.12`.}
+
+\item{pointSize}{Numeric point size for the cell centroids. If supplied,
+takes precedence over `size`.}
+
+\item{pointAlpha}{Numeric transparency for the cell centroids. If supplied,
+takes precedence over `alpha`.}
+
+\item{numberSize}{Numeric size for the FoV labels. Default: `1`.}
+
+\item{numbersAlpha}{Numeric transparency for FoV labels. If supplied, takes
+precedence over `alphaNumbers`.}
}
\value{
A `ggplot` object showing cell centroids and FoV boundaries.
diff --git a/man/plotZoomFovsMap.Rd b/man/plotZoomFovsMap.Rd
index d19a97d..cc32161 100644
--- a/man/plotZoomFovsMap.Rd
+++ b/man/plotZoomFovsMap.Rd
@@ -9,17 +9,19 @@ plotZoomFovsMap(
fovs = NULL,
title = NULL,
mapPointCol = "darkmagenta",
- mapPointSize = 0.5,
- mapPointAlpha = 0.8,
- fovNumbersCol = "black",
- fovNumberSize = 1,
- fovNumbersAlpha = 0.8,
+ mapNumbersCol = "black",
+ mapAlphaNumbers = 0.8,
csize = 0.05,
calpha = 0.8,
scaleBars = NULL,
scaleBarMap = TRUE,
scaleBarPol = TRUE,
- ...
+ ...,
+ mapPointSize = NULL,
+ mapPointAlpha = NULL,
+ fovNumbersCol = NULL,
+ fovNumberSize = 1,
+ fovNumbersAlpha = NULL
)
}
\arguments{
@@ -35,16 +37,13 @@ plot. If `NULL`, no title is added. Default is `NULL`.}
\item{mapPointCol}{A character string specifying the color of the points
in the map. Default is `"darkmagenta"`.}
-\item{mapPointSize}{Numeric size for points in the map. Default: `0.5`.}
+\item{mapNumbersCol}{Deprecated alias for `fovNumbersCol`.}
-\item{mapPointAlpha}{Numeric transparency for points in the map. Default: `0.8`.}
+\item{mapAlphaNumbers}{Deprecated alias for `fovNumbersAlpha`.}
-\item{fovNumbersCol}{A character string specifying the color of the
-numbers on the FoV zoom-in. Default is `"black"`.}
+\item{csize}{Deprecated alias for `mapPointSize`.}
-\item{fovNumberSize}{Numeric size for the FoV labels. Default: `1`.}
-
-\item{fovNumbersAlpha}{Numeric transparency for FoV labels. Default: `0.8`.}
+\item{calpha}{Deprecated alias for `mapPointAlpha`.}
\item{scaleBars}{Logical or NULL. Default is `NULL`.
Master switch controlling the presence of scale bars in both panels.
@@ -60,6 +59,20 @@ These parameters are only used when \code{scaleBars} is \code{NULL};
otherwise they are overridden by \code{scaleBars}.}
\item{...}{Additional arguments passed to `plotPolygons`.}
+
+\item{mapPointSize}{Numeric size for points in the map. If supplied, takes
+precedence over `csize`.}
+
+\item{mapPointAlpha}{Numeric transparency for points in the map. If supplied,
+takes precedence over `calpha`.}
+
+\item{fovNumbersCol}{Color for FoV labels. If supplied, takes precedence over
+`mapNumbersCol`.}
+
+\item{fovNumberSize}{Numeric size for the FoV labels. Default: `1`.}
+
+\item{fovNumbersAlpha}{Transparency for FoV labels. If supplied, takes
+precedence over `mapAlphaNumbers`.}
}
\value{
A combined plot showing a map of all FOVs with zoomed-in views of
diff --git a/man/trainModel.Rd b/man/trainModel.Rd
index 769a52a..b5a913e 100644
--- a/man/trainModel.Rd
+++ b/man/trainModel.Rd
@@ -11,8 +11,8 @@ trainModel(modelMatrix, trainDF)
with `getModelFormula` and `model.matrix` functions.}
\item{trainDF}{`data.frame`
-A data frame containing at least the response column
-`qscore_train`, coded as 0/1.}
+A data frame containing `QScore_train` or historical `qcscore_train`,
+coded as 0/1. If both are present, they must agree.}
}
\value{
A \code{\link[glmnet]{glmnet}} model object fitted with
diff --git a/tests/testthat/test_QCScores.R b/tests/testthat/test_QCScores.R
index 9e59e9a..51880c8 100644
--- a/tests/testthat/test_QCScores.R
+++ b/tests/testthat/test_QCScores.R
@@ -1,19 +1,31 @@
library(testthat)
library(SpaceTrooper)
-# load the example SpatialExperiment
-spe0 <- example(readCosmxSPE)$value
+# Load the example SpatialExperiment directly so tests do not depend on an
+# installed help database.
+cosmx_path <- system.file(
+ "extdata",
+ "CosMx_DBKero_Tiny",
+ package="SpaceTrooper"
+)
+spe0 <- readCosmxSPE(cosmx_path, sampleName="DBKero_Tiny")
test_that("QC functions are exported", {
expect_true(exists("spatialPerCellQC", mode = "function"))
- expect_true(exists("computeQCScore", mode = "function"))
+ expect_true(exists("computeQScore", mode = "function"))
expect_true(exists("computeSpatialOutlier", mode = "function"))
- expect_true(exists("computeQCScoreFlags", mode = "function"))
+ expect_true(exists("computeQScoreFlags", mode = "function"))
expect_true(exists("computeThresholdFlags", mode = "function"))
+ ## deprecated wrappers should still be exported
+ expect_true(exists("computeQCScore", mode = "function"))
+ expect_true(exists("computeQCScoreFlags", mode = "function"))
+ expect_true(exists("computeOutliersQCScore", mode = "function"))
+ expect_true(exists("applyQScoreModel", mode = "function"))
+ expect_true(exists("applyQCScoreModel", mode = "function"))
})
-test_that("spatialPerCellQC adds per‐cell metrics to colData", {
+test_that("spatialPerCellQC adds per-cell metrics to colData", {
spe <- spatialPerCellQC(spe0, micronConvFact = 0.15)
expect_s4_class(spe, "SpatialExperiment")
cd <- colData(spe)
@@ -23,16 +35,46 @@ test_that("spatialPerCellQC adds per‐cell metrics to colData", {
expect_true(all(required %in% colnames(cd)))
})
-test_that("computeQCScore adds a flag_score between 0 and 1", {
+test_that("computeQScore adds a QScore between 0 and 1", {
spe <- spatialPerCellQC(spe0)
- spe2 <- computeQCScore(spe)
+ set.seed(42)
+ spe2 <- computeQScore(spe)
cd2 <- colData(spe2)
- expect_true("QC_score" %in% colnames(cd2))
- fs <- cd2$QC_score
+ expect_true("QScore" %in% colnames(cd2))
+ fs <- cd2$QScore
expect_true(is.numeric(fs))
- expect_true(all(fs >= 0 & fs <= 1))
+ expect_true(all(fs[!is.na(fs)] >= 0 & fs[!is.na(fs)] <= 1))
+})
+
+test_that("computeQCScore preserves legacy score and model names", {
+ spe <- spatialPerCellQC(spe0)
+ set.seed(42)
+ expect_warning(
+ spe2 <- computeQCScore(spe),
+ "deprecated"
+ )
+ expect_true("QC_score" %in% colnames(colData(spe2)))
+ expect_false("QScore" %in% colnames(colData(spe2)))
+ expect_true("QCScore_model" %in% names(metadata(spe2)))
+ expect_false("QScore_model" %in% names(metadata(spe2)))
})
+test_that("computeQCScoreFlags preserves legacy flag names", {
+ spe <- spe0
+ spe$QC_score <- seq(0, 1, length.out=ncol(spe))
+ spe$threshold_flags <- rep(c(TRUE, FALSE), length.out=ncol(spe))
+
+ expect_warning(
+ flagged <- computeQCScoreFlags(spe, qsThreshold=0.5),
+ "deprecated"
+ )
+
+ expect_true(all(c(
+ "low_qcscore",
+ "low_threshold_qcscore"
+ ) %in% colnames(colData(flagged))))
+ expect_false("low_QScore" %in% colnames(colData(flagged)))
+})
test_that("computeSpatialOutlier flags outliers for a chosen metric", {
spe <- spatialPerCellQC(spe0)
@@ -44,9 +86,10 @@ test_that("computeSpatialOutlier flags outliers for a chosen metric", {
})
-test_that("computeQCScoreFlags combines filters and returns filter_out", {
+test_that("computeQScoreFlags combines filters and returns filter_out", {
spe <- spatialPerCellQC(spe0)
- spe <- computeQCScore(spe)
+ set.seed(42)
+ spe <- computeQScore(spe)
ff <- computeThresholdFlags(spe,
totalThreshold = 10,
ctrlTotRatioThreshold = 0.2)
@@ -61,3 +104,120 @@ test_that("computeQCScoreFlags combines filters and returns filter_out", {
combined <- (ff$is_zero_counts & ff$is_ctrl_tot_outlier)
expect_identical(combined, cd_ff$threshold_flags)
})
+
+test_that("getModelFormula preserves historical calls", {
+ formula_vars <- c(
+ log2SignalDensity="log2SignalDensity_outlier_train",
+ Area_um="Area_um_outlier_sc"
+ )
+
+ expect_identical(
+ getModelFormula(formula_vars),
+ "~(log2SignalDensity + Area_um)^2"
+ )
+ expect_message(
+ getModelFormula(formulaVars=formula_vars, verbose=TRUE),
+ "Final formula"
+ )
+ expect_identical(
+ getModelFormula(metricList=names(formula_vars)),
+ "~(log2SignalDensity + Area_um)^2"
+ )
+ expect_warning(
+ getModelFormula(
+ formulaVars=formula_vars,
+ metricList="log2SignalDensity"
+ ),
+ "using 'metricList'"
+ )
+})
+
+test_that("plotCellsFovs old and canonical arguments have equal effects", {
+ old <- plotCellsFovs(
+ spe0,
+ size=2,
+ alpha=0.4,
+ alphaNumbers=0.3,
+ scaleBar=FALSE
+ )
+ canonical <- plotCellsFovs(
+ spe0,
+ pointSize=2,
+ pointAlpha=0.4,
+ numbersAlpha=0.3,
+ scaleBar=FALSE
+ )
+
+ expect_equal(old$layers[[1]]$aes_params, canonical$layers[[1]]$aes_params)
+ expect_equal(old$layers[[3]]$aes_params, canonical$layers[[3]]$aes_params)
+
+ historical_positional <- plotCellsFovs(
+ spe0,
+ unique(spe0$sample_id),
+ "firebrick",
+ "black",
+ 0.3,
+ metadata(spe0)$fov_dim,
+ 2,
+ 0.4,
+ FALSE,
+ 0.12
+ )
+ expect_equal(
+ historical_positional$layers[[1]]$aes_params$size,
+ 2
+ )
+ expect_warning(
+ plotCellsFovs(
+ spe0,
+ size=1,
+ pointSize=2,
+ scaleBar=FALSE
+ ),
+ "using 'pointSize'"
+ )
+})
+
+test_that("plotZoomFovsMap old and canonical arguments are equivalent", {
+ cosmx_polygons <- readCosmxSPE(
+ cosmx_path,
+ sampleName="DBKero_Tiny",
+ keepPolygons=TRUE
+ )
+ fov <- unique(cosmx_polygons$fov)[1]
+
+ old <- plotZoomFovsMap(
+ cosmx_polygons,
+ fovs=fov,
+ mapNumbersCol="blue",
+ mapAlphaNumbers=0.3,
+ csize=0.7,
+ calpha=0.4,
+ scaleBars=FALSE
+ )
+ canonical <- plotZoomFovsMap(
+ cosmx_polygons,
+ fovs=fov,
+ fovNumbersCol="blue",
+ fovNumbersAlpha=0.3,
+ mapPointSize=0.7,
+ mapPointAlpha=0.4,
+ scaleBars=FALSE
+ )
+
+ expect_true(isTRUE(all.equal(
+ old,
+ canonical,
+ check.environment=FALSE
+ )))
+ expect_warning(
+ plotZoomFovsMap(
+ cosmx_polygons,
+ fovs=fov,
+ csize=0.5,
+ mapPointSize=0.7,
+ scaleBars=FALSE
+ ),
+ "using 'mapPointSize'"
+ )
+})
diff --git a/tests/testthat/test_computeLambda.R b/tests/testthat/test_computeLambda.R
new file mode 100644
index 0000000..ce934b0
--- /dev/null
+++ b/tests/testthat/test_computeLambda.R
@@ -0,0 +1,119 @@
+library(testthat)
+library(SpaceTrooper)
+
+lambda_train_data <- function(response_name="QScore_train") {
+ set.seed(101)
+ train_df <- data.frame(
+ x1=rnorm(80),
+ x2=rnorm(80),
+ x3=rnorm(80)
+ )
+ response <- as.integer(
+ train_df$x1 - 0.5 * train_df$x2 + rnorm(80) > 0
+ )
+ train_df[[response_name]] <- response
+ train_df
+}
+
+test_that("computeLambda preserves historical positional and named calls", {
+ train_df <- lambda_train_data("qcscore_train")
+ model_formula <- "~ x1 + x2"
+
+ set.seed(11)
+ positional <- computeLambda(train_df, model_formula)
+ set.seed(11)
+ named <- computeLambda(
+ trainDF=train_df,
+ modelFormula=model_formula
+ )
+
+ expect_type(positional, "double")
+ expect_length(positional, 1L)
+ expect_true(is.finite(positional))
+ expect_identical(positional, named)
+})
+
+test_that("computeLambda accepts canonical and legacy responses", {
+ canonical <- lambda_train_data("QScore_train")
+ legacy <- canonical
+ names(legacy)[names(legacy) == "QScore_train"] <- "qcscore_train"
+
+ set.seed(22)
+ canonical_lambda <- computeLambda(canonical, ~ x1 + x2)
+ set.seed(22)
+ legacy_lambda <- computeLambda(legacy, ~ x1 + x2)
+ expect_identical(canonical_lambda, legacy_lambda)
+
+ both <- canonical
+ both$qcscore_train <- both$QScore_train
+ set.seed(22)
+ expect_identical(
+ computeLambda(both, ~ x1 + x2),
+ canonical_lambda
+ )
+
+ both$qcscore_train[1] <- 1L - both$QScore_train[1]
+ expect_error(
+ computeLambda(both, ~ x1 + x2),
+ "different values"
+ )
+})
+
+test_that("computeLambda filters complete cases", {
+ train_df <- lambda_train_data()
+ train_df$x2[c(2, 7)] <- NA_real_
+ complete_df <- train_df[complete.cases(
+ train_df[, c("x1", "x2", "QScore_train")]
+ ), ]
+
+ set.seed(33)
+ expect_warning(
+ with_missing <- computeLambda(train_df, ~ x1 + x2),
+ "2 training cells"
+ )
+ set.seed(33)
+ without_missing <- computeLambda(complete_df, ~ x1 + x2)
+
+ expect_identical(with_missing, without_missing)
+})
+
+test_that("computeLambda preserves additive and selected interaction formulas", {
+ train_df <- lambda_train_data()
+
+ set.seed(44)
+ additive <- computeLambda(train_df, ~ x1 + x2 + x3)
+ set.seed(44)
+ selected_interaction <- computeLambda(
+ train_df,
+ ~ x1 + x2 + x1:x3
+ )
+
+ expect_true(is.finite(additive))
+ expect_true(is.finite(selected_interaction))
+})
+
+test_that("public computeLambda agrees with matrix implementation and glmnet", {
+ train_df <- lambda_train_data()
+ model_formula <- ~ x1 + x2 + x1:x2
+ model_matrix <- stats::model.matrix(model_formula, data=train_df)
+ response <- train_df$QScore_train
+
+ set.seed(55)
+ public <- computeLambda(train_df, model_formula)
+ set.seed(55)
+ internal <- SpaceTrooper:::.computeLambda(
+ modelMatrix=model_matrix,
+ response=response
+ )
+ set.seed(55)
+ reference <- glmnet::cv.glmnet(
+ model_matrix,
+ response,
+ family="binomial",
+ alpha=0,
+ lambda=NULL
+ )$lambda.min
+
+ expect_identical(public, internal)
+ expect_identical(public, reference)
+})
diff --git a/tests/testthat/test_modelFormula.R b/tests/testthat/test_modelFormula.R
new file mode 100644
index 0000000..a07aa62
--- /dev/null
+++ b/tests/testthat/test_modelFormula.R
@@ -0,0 +1,326 @@
+library(testthat)
+library(SpaceTrooper)
+
+qscore_formula_spe <- function() {
+ path <- system.file(
+ "extdata",
+ "CosMx_DBKero_Tiny",
+ package="SpaceTrooper"
+ )
+ spatialPerCellQC(
+ readCosmxSPE(path, sampleName="DBKero_Tiny")
+ )
+}
+
+model_terms <- function(spe) {
+ attr(
+ stats::terms(stats::as.formula(
+ metadata(spe)$QScore_model$model_formula
+ )),
+ "term.labels"
+ )
+}
+
+captured_training_metrics <- function(spe, model_formula, seed=206) {
+ seen <- new.env(parent=emptyenv())
+ original_prep <- SpaceTrooper:::.prepQCContext
+ testthat::local_mocked_bindings(
+ .prepQCContext=function(spe, metricList, verbose=FALSE) {
+ seen$metricList <- metricList
+ original_prep(spe, metricList, verbose)
+ },
+ .package="SpaceTrooper"
+ )
+ set.seed(seed)
+ scored <- computeQScore(
+ spe,
+ bestLambda=0.01,
+ modelFormula=model_formula
+ )
+ list(metrics=seen$metricList, scored=scored)
+}
+
+test_that("supported Quality Score predictors are explicit and stable", {
+ expect_identical(
+ SpaceTrooper:::.qscoreSupportedPredictors(),
+ c(
+ "log2SignalDensity",
+ "Area_um",
+ "log2AspectRatio",
+ "log2Ctrl_total_ratio"
+ )
+ )
+})
+
+test_that("default formula keeps the established pairwise interactions", {
+ spe <- qscore_formula_spe()
+ set.seed(201)
+ scored <- computeQScore(spe)
+ terms <- model_terms(scored)
+
+ expect_true("log2SignalDensity" %in% terms)
+ expect_true("Area_um" %in% terms)
+ expect_true(any(grepl(":", terms, fixed=TRUE)))
+})
+
+test_that("additive custom formula remains additive", {
+ spe <- qscore_formula_spe()
+ supplied <- "~ log2SignalDensity + Area_um"
+ set.seed(202)
+ scored <- computeQScore(spe, modelFormula=supplied)
+
+ expect_identical(
+ metadata(scored)$QScore_model$model_formula,
+ supplied
+ )
+ expect_identical(
+ model_terms(scored),
+ c("log2SignalDensity", "Area_um")
+ )
+ expect_identical(
+ metadata(scored)$QScore_model$model_matrix_colnames,
+ c("(Intercept)", "log2SignalDensity", "Area_um")
+ )
+})
+
+test_that("custom formula preserves selected interactions and column order", {
+ spe <- qscore_formula_spe()
+ supplied <- paste(
+ "~ log2SignalDensity + Area_um +",
+ "log2SignalDensity:log2Ctrl_total_ratio"
+ )
+ set.seed(203)
+ scored <- computeQScore(spe, modelFormula=supplied)
+
+ expect_identical(
+ model_terms(scored),
+ c(
+ "log2SignalDensity",
+ "Area_um",
+ "log2SignalDensity:log2Ctrl_total_ratio"
+ )
+ )
+ expect_identical(
+ metadata(scored)$QScore_model$model_matrix_colnames,
+ c(
+ "(Intercept)",
+ "log2SignalDensity",
+ "Area_um",
+ "log2SignalDensity:log2Ctrl_total_ratio"
+ )
+ )
+})
+
+test_that("custom formula is not rebuilt by getModelFormula", {
+ spe <- qscore_formula_spe()
+ testthat::local_mocked_bindings(
+ getModelFormula=function(...) stop("getModelFormula was called"),
+ .package="SpaceTrooper"
+ )
+
+ set.seed(204)
+ expect_no_error(computeQScore(
+ spe,
+ modelFormula=~ log2SignalDensity + Area_um
+ ))
+})
+
+test_that("custom formulas define metrics passed to QScore preparation", {
+ spe <- qscore_formula_spe()
+
+ single <- captured_training_metrics(
+ spe,
+ ~ log2SignalDensity,
+ seed=206
+ )
+ expect_identical(single$metrics, "log2SignalDensity")
+
+ additive <- captured_training_metrics(
+ spe,
+ ~ log2SignalDensity + Area_um,
+ seed=207
+ )
+ expect_identical(
+ additive$metrics,
+ c("log2SignalDensity", "Area_um")
+ )
+
+ interaction <- captured_training_metrics(
+ spe,
+ ~ log2SignalDensity * Area_um,
+ seed=208
+ )
+ expect_identical(
+ interaction$metrics,
+ c("log2SignalDensity", "Area_um")
+ )
+ expect_identical(
+ model_terms(interaction$scored),
+ c(
+ "log2SignalDensity",
+ "Area_um",
+ "log2SignalDensity:Area_um"
+ )
+ )
+
+ selected <- captured_training_metrics(
+ spe,
+ ~ log2SignalDensity +
+ log2SignalDensity:log2Ctrl_total_ratio,
+ seed=209
+ )
+ expect_identical(
+ selected$metrics,
+ c("log2SignalDensity", "log2Ctrl_total_ratio")
+ )
+})
+
+test_that("default formula keeps the established training metrics", {
+ spe <- qscore_formula_spe()
+ default <- captured_training_metrics(spe, NULL, seed=210)
+
+ expect_identical(
+ default$metrics,
+ SpaceTrooper:::.qscoreSupportedPredictors()
+ )
+})
+
+test_that("supported CosMx border expression is accepted", {
+ spe <- qscore_formula_spe()
+ supplied <- paste0(
+ "~ log2SignalDensity + ",
+ "I(abs(log2AspectRatio) * as.numeric(dist_border < 50))"
+ )
+ set.seed(205)
+ scored <- computeQScore(spe, modelFormula=supplied)
+
+ expect_identical(
+ metadata(scored)$QScore_model$model_formula,
+ supplied
+ )
+ formula_info <- SpaceTrooper:::.validateQScoreFormula(
+ modelFormula=supplied,
+ dataNames=names(colData(spe)),
+ technology=metadata(spe)$technology
+ )
+ expect_identical(
+ formula_info$training_metrics,
+ c("log2SignalDensity", "log2AspectRatio")
+ )
+ expect_false("dist_border" %in% formula_info$training_metrics)
+})
+
+test_that("custom metric subsets change deterministic training labels", {
+ n_cells <- 100L
+ training_data <- data.frame(
+ cell_id=paste0("cell", seq_len(n_cells)),
+ log2SignalDensity=seq_len(n_cells),
+ Area_um=seq_len(n_cells),
+ log2SignalDensity_outlier_train=c(
+ "LOW",
+ rep("NO", n_cells - 1L)
+ ),
+ Area_um_outlier_sc=c(
+ "NO",
+ "HIGH",
+ rep("NO", n_cells - 2L)
+ )
+ )
+ outlier_columns <- c(
+ log2SignalDensity="log2SignalDensity_outlier_train",
+ Area_um="Area_um_outlier_sc"
+ )
+ formula_metrics <- function(model_formula) {
+ SpaceTrooper:::.validateQScoreFormula(
+ modelFormula=model_formula,
+ dataNames=names(training_data),
+ technology="10X_Xenium"
+ )$training_metrics
+ }
+
+ signal_metrics <- formula_metrics(~ log2SignalDensity)
+ subset_metrics <- formula_metrics(
+ ~ log2SignalDensity + Area_um
+ )
+
+ set.seed(211)
+ signal_train <- computeTrainDF(
+ training_data,
+ outlier_columns[signal_metrics],
+ tech="10X_Xenium"
+ )
+ set.seed(211)
+ subset_train <- computeTrainDF(
+ training_data,
+ outlier_columns[subset_metrics],
+ tech="10X_Xenium"
+ )
+
+ expect_identical(
+ sort(signal_train$cell_id[signal_train$QScore_train == 0]),
+ "cell1"
+ )
+ expect_identical(
+ sort(subset_train$cell_id[subset_train$QScore_train == 0]),
+ c("cell1", "cell2")
+ )
+})
+
+test_that("unsupported predictors and transformations are rejected", {
+ spe <- qscore_formula_spe()
+
+ expect_error(
+ computeQScore(
+ spe,
+ modelFormula=~ log2SignalDensity + customMetric
+ ),
+ "Unsupported Quality Score predictor.*customMetric"
+ )
+ expect_error(
+ computeQScore(
+ spe,
+ modelFormula=~ log2SignalDensity + log1p(Area_um)
+ ),
+ "Unsupported Quality Score formula term.*log1p"
+ )
+
+ spe$customMetric <- seq_len(ncol(spe))
+ expect_error(
+ computeOutliersQScore(
+ spe,
+ metricList=c("log2SignalDensity", "customMetric")
+ ),
+ "Unsupported Quality Score metric.*customMetric"
+ )
+})
+
+test_that("custom formulas retain the signal-density requirement", {
+ spe <- qscore_formula_spe()
+
+ expect_error(
+ computeQScore(spe, modelFormula=~ Area_um),
+ "'log2SignalDensity' is required.*training labels"
+ )
+})
+
+test_that("missing and technology-incompatible variables are reported", {
+ spe <- qscore_formula_spe()
+ spe$Area_um <- NULL
+ expect_error(
+ computeQScore(
+ spe,
+ modelFormula=~ log2SignalDensity + Area_um
+ ),
+ "Missing variables required.*Area_um"
+ )
+
+ xenium_like <- qscore_formula_spe()
+ metadata(xenium_like)$technology <- "10X_Xenium"
+ expect_error(
+ computeQScore(
+ xenium_like,
+ modelFormula=~ log2SignalDensity + log2AspectRatio
+ ),
+ "supported only for Nanostring CosMx"
+ )
+})
diff --git a/tests/testthat/test_modelTransfer.R b/tests/testthat/test_modelTransfer.R
new file mode 100644
index 0000000..446abc5
--- /dev/null
+++ b/tests/testthat/test_modelTransfer.R
@@ -0,0 +1,156 @@
+library(testthat)
+library(SpaceTrooper)
+
+qscore_transfer_data <- function() {
+ path <- system.file(
+ "extdata",
+ "CosMx_DBKero_Tiny",
+ package="SpaceTrooper"
+ )
+ spe <- spatialPerCellQC(
+ readCosmxSPE(path, sampleName="DBKero_Tiny")
+ )
+ split <- floor(ncol(spe) * 0.7)
+ list(
+ train=spe[, seq_len(split)],
+ query=spe[, seq.int(split + 1L, ncol(spe))]
+ )
+}
+
+test_that("applyQScoreModel transfers a supported custom model", {
+ data <- qscore_transfer_data()
+ supplied <- paste(
+ "~ log2SignalDensity + Area_um +",
+ "log2SignalDensity:log2Ctrl_total_ratio"
+ )
+ set.seed(301)
+ trained <- computeQScore(
+ data$train,
+ modelFormula=supplied
+ )
+ qscore_model <- metadata(trained)$QScore_model
+
+ applied <- applyQScoreModel(
+ data$query,
+ qsModel=qscore_model
+ )
+
+ expect_true("QScore" %in% colnames(colData(applied)))
+ expect_true(all(
+ applied$QScore[!is.na(applied$QScore)] >= 0 &
+ applied$QScore[!is.na(applied$QScore)] <= 1
+ ))
+ expect_identical(
+ metadata(applied)$QScore_model_applied$model_matrix_colnames,
+ qscore_model$model_matrix_colnames
+ )
+ expect_identical(
+ qscore_model$model_matrix_colnames,
+ c(
+ "(Intercept)",
+ "log2SignalDensity",
+ "Area_um",
+ "log2SignalDensity:log2Ctrl_total_ratio"
+ )
+ )
+})
+
+test_that("deprecated applyQCScoreModel preserves legacy output names", {
+ data <- qscore_transfer_data()
+ set.seed(302)
+ trained <- computeQScore(
+ data$train,
+ modelFormula=~ log2SignalDensity + Area_um
+ )
+ qscore_model <- metadata(trained)$QScore_model
+
+ expect_warning(
+ applied <- applyQCScoreModel(
+ data$query,
+ qcModel=qscore_model
+ ),
+ "deprecated"
+ )
+
+ expect_true("QC_score" %in% colnames(colData(applied)))
+ expect_false("QScore" %in% colnames(colData(applied)))
+ expect_true("QCScore_model_applied" %in% names(metadata(applied)))
+ expect_false("QScore_model_applied" %in% names(metadata(applied)))
+})
+
+test_that("model transfer rejects technology-incompatible border terms", {
+ data <- qscore_transfer_data()
+ border_formula <- paste0(
+ "~ log2SignalDensity + ",
+ "I(abs(log2AspectRatio) * as.numeric(dist_border < 50))"
+ )
+ set.seed(303)
+ trained <- computeQScore(
+ data$train,
+ modelFormula=border_formula
+ )
+
+ metadata(data$query)$technology <- "10X_Xenium"
+ expect_error(
+ applyQScoreModel(
+ data$query,
+ qsModel=metadata(trained)$QScore_model
+ ),
+ "supported only for Nanostring CosMx"
+ )
+})
+
+test_that("computeTrainDF and trainModel accept legacy training response", {
+ data <- qscore_transfer_data()
+ outliers <- computeOutliersQScore(data$train)
+ outliers <- checkOutliers(outliers)
+ set.seed(304)
+ train_df <- computeTrainDF(
+ colData(outliers),
+ metadata(outliers)$formula_variables,
+ metadata(outliers)$technology
+ )
+
+ expect_identical(
+ train_df$QScore_train,
+ train_df$qcscore_train
+ )
+
+ matrix <- stats::model.matrix(
+ ~ log2SignalDensity + Area_um,
+ data=train_df
+ )
+ legacy_only <- train_df
+ legacy_only$QScore_train <- NULL
+ expect_s3_class(
+ trainModel(matrix, legacy_only),
+ "glmnet"
+ )
+
+ conflicting <- train_df
+ conflicting$qcscore_train[1] <-
+ 1L - conflicting$QScore_train[1]
+ expect_error(
+ trainModel(matrix, conflicting),
+ "different values"
+ )
+})
+
+test_that("stored coefficient metadata contains the complete vector", {
+ data <- qscore_transfer_data()
+ set.seed(305)
+ trained <- computeQScore(
+ data$train,
+ modelFormula=~ log2SignalDensity + Area_um
+ )
+ model <- metadata(trained)$QScore_model
+
+ expect_equal(
+ nrow(model$coefficients),
+ length(model$model_matrix_colnames) + 1L
+ )
+ expect_equal(
+ nrow(model$coefficients_table),
+ nrow(model$coefficients)
+ )
+})
diff --git a/vignettes/SpaceTrooper_utilities.Rmd b/vignettes/SpaceTrooper_utilities.Rmd
index 8de8a96..b9e18b1 100644
--- a/vignettes/SpaceTrooper_utilities.Rmd
+++ b/vignettes/SpaceTrooper_utilities.Rmd
@@ -471,20 +471,17 @@ all pairwise interactions.
mandatory metric in the QS formula. If an insufficient number of outliers is
detected for this metric (fewer than 0.1% of the dataset after excluding
zero-count cells), QS computation cannot proceed using the remaining metrics.
-In such cases, the provided code will still add a`QScore` column to `colData`.
+In such cases, the provided code will still add a `QScore` column to `colData`.
This column is populated if the minimum requirement is met and contains `NA`
values otherwise.
```{r compute-QS-safe-run, message=TRUE}
-# safe run function
-safe_run <- function(expr) {
-tryCatch(
- list(result=expr, error=NULL),
+# safe run using tryCatch directly so errors are caught before evaluation
+out <- tryCatch(
+ list(result=computeQScore(spe), error=NULL),
error=function(e) list(result=NULL, error=e)
- )
-}
+)
-out <- safe_run(computeQScore(spe))
if (!is.null(out$error)) {
message("Failed: ", out$error$message)
colData(spe)$QScore <- NA
@@ -496,30 +493,37 @@ spe <- out$result
### Compute QS with custom formula
In `computeQScore`, the `modelFormula` parameter allows users to specify which
-metrics to include in the QS formula among the 4 currently supported, i.e.
-`log2SignalDensity`, `Area_um`, `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`, and `log2Ctrl_total_ratio`.
-Border effect still cannot be considered for Xenium and MERFISH/MERSCOPE datasets.
+metrics to include in the QS formula among the four currently supported:
+`log2SignalDensity`, `Area_um`, `log2AspectRatio`, and
+`log2Ctrl_total_ratio`. The supported border-effect expression is
+`I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`; aspect ratio and
+border effect are available only for CosMx datasets.
+
The default is `NULL` and the formula is automatically estimated based on the
available metrics, their outliers and the dataset technology. The most complete
formula is generated as follows:
`~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2`
- When user-provided, the formula must follow the same default syntax and the terms
- must be written exactly as they appear here (e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`
- must have spaces around the `*` and `<` operators). However, the terms do not need
- to follow any specific order. If the outliers for one or more metrics are not
- sufficient, they will not be considered for model training and the corresponding terms
- are automatically dropped from the formula.
- We provide interaction terms by default as they contribute to model flexibility.
- However, if users prefer to exclude interaction terms, they can do so by
- removing the `()^2` notation from the formula.
+When user-provided, the formula is fitted as supplied. It must include
+`log2SignalDensity`; users may select a subset of the other supported metrics,
+use an additive formula, or choose interactions with `:` and `*`. SpaceTrooper
+does not add or remove terms. Predictor and transformation
+validation is separate from formula structure: arbitrary metrics and
+transformations such as `log1p(Area_um)` are rejected. The supported base QC
+metrics present in the custom formula are also used to compute outliers, define
+good and bad training cells, and construct `QScore_train`; supported metrics
+omitted from the formula do not contribute to those labels. Interactions alter
+the fitted model but do not create additional training metrics. In the supported
+CosMx border expression, `log2AspectRatio` is the training QC metric and
+`dist_border` is only an auxiliary variable. With `modelFormula=NULL`, the
+established default metric-selection and training-label workflow is unchanged.
```{r compute-QS-different-formula, message=TRUE}
set.seed(1713)
-model_formula <- "~log2Ctrl_total_ratio + log2SignalDensity"
+model_formula <- "~ log2Ctrl_total_ratio + log2SignalDensity"
spe <- computeQScore(spe, modelFormula=model_formula, verbose=TRUE)
```
@@ -533,6 +537,27 @@ Compared to the plot displayed in [RNA](https://bioconductor.org/packages/devel/
QS values are higher for bigger cells and cells located on the FoV borders, as the
cell size and border effect contributions are not considered in the QS computation.
+### Apply a trained QS model
+
+The fitted model is stored in `metadata(spe)$QScore_model` and can be applied
+with `applyQScoreModel()` to an object containing the required QC metrics.
+Model transfer across datasets is generally discouraged because QS training is
+dataset-specific; when it is necessary, use datasets from the same experiment
+and technology, or train a supported formula that excludes incompatible terms.
+
+```{r apply-QS-model, message=TRUE}
+qs_model <- metadata(spe)$QScore_model
+spe <- applyQScoreModel(
+ spe,
+ qsModel=qs_model,
+ scoreName="QScore_transferred"
+)
+```
+
+The historical `applyQCScoreModel()` interface remains available with a
+deprecation warning and retains `qcModel`, `QC_score`, and
+`QCScore_model_applied` naming.
+
### Single metric flagging
The `SpaceTrooper` package allows the identification of numerically aberrant