From d62f72888614eb9d8c98569a088c5711b6a0315d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 24 Jul 2026 10:05:38 +0000 Subject: [PATCH 01/13] Initial plan From 7fd56ab75c0d0e46c38b745b43cc3a1039ac1675 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 24 Jul 2026 10:12:28 +0000 Subject: [PATCH 02/13] Fix code review issues: deprecated wrappers, formula handling, docs, tests, vignette --- NAMESPACE | 2 + R/QC.R | 110 ++++++++++++++++++++------ man/computeLambda.Rd | 16 ++-- man/computeQCScore-deprecated.Rd | 23 ++++++ man/computeQCScoreFlags-deprecated.Rd | 21 +++++ man/computeQScore.Rd | 6 +- man/dot-applyQScoreModel.Rd | 6 +- man/getModelFormula.Rd | 15 ++-- tests/testthat/test_QCScores.R | 35 +++++--- vignettes/SpaceTrooper_utilities.Rmd | 22 +++--- 10 files changed, 186 insertions(+), 70 deletions(-) create mode 100644 man/computeQCScore-deprecated.Rd create mode 100644 man/computeQCScoreFlags-deprecated.Rd diff --git a/NAMESPACE b/NAMESPACE index 5911763..ccdb3b0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,8 @@ export(computeLambda) export(computeMissingMetricsMerfish) export(computeMissingMetricsXenium) export(computeOutliersQScore) +export(computeQCScore) +export(computeQCScoreFlags) export(computeQScore) export(computeQScoreFlags) export(computeSpatialOutlier) diff --git a/R/QC.R b/R/QC.R index 7bd3e76..c54fed4 100644 --- a/R/QC.R +++ b/R/QC.R @@ -349,13 +349,13 @@ computeThresholdFlags <- function(spe, totalThreshold=0, #' k-fold cross-validation to identify the optimal regularization parameter #' \eqn{\lambda} for Quality Score (QS) model training. #' +#' @param modelMatrix `matrix` +#' The design matrix built from the training data, typically via +#' `model.matrix(as.formula(model_formula), data=trainDF)`. #' @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` A binary (0/1) response vector to be modeled. #' #' @return #' `numeric` @@ -364,14 +364,14 @@ computeThresholdFlags <- function(spe, totalThreshold=0, #' #' @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`. #' #' @examples #' example(computeTrainDF) -#' modform <- getModelFormula(metadata(spe)$formula_variables) -#' best_lambda <- computeLambda(df_train, modform) +#' modform <- getModelFormula(names(metadata(spe)$formula_variables)) +#' model_matrix <- model.matrix(as.formula(modform), data=df_train) +#' best_lambda <- computeLambda(model_matrix, df_train) #' print(best_lambda) #' #' @@ -432,9 +432,9 @@ computeLambda <- function(modelMatrix, trainDF) { #' 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 -#' 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. +#' removed (or added) terms should be written as in the default formula, +#' e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`. +#' Whitespace around operators is accepted. #' In any case, metrics with insufficient outliers (less than 0.1\% of the dataset) #' will be excluded from the QS formula. #' @@ -478,9 +478,10 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE 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) + ## Whitespace-tolerant detection of the border-effect interaction term + border_pat <- "I\\(abs\\(log2AspectRatio\\)\\s*\\*\\s*as\\.numeric\\(dist_border\\s*<\\s*50\\)\\)" + if (any(grepl(border_pat, metricList))) { + metricList <- gsub(border_pat, "log2AspectRatio", metricList) } } @@ -488,7 +489,9 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE ctx <- .prepQCContext(spe, metricList, 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)) + } if (verbose) { message("Using final model formula:") @@ -752,24 +755,23 @@ 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 metricList A character vector of metric names to include in the +#' formula (e.g. `"log2SignalDensity"`, `"Area_um"`, etc.), typically the +#' names of `metadata(spe)$formula_variables`. #' @return `character` #' A one‐sided formula as a string (e.g. "~ log2SignalDensity + ..."). #' @export #' @examples #' example(checkOutliers) -#' getModelFormula(metadata(spe)$formula_variables) +#' getModelFormula(names(metadata(spe)$formula_variables)) getModelFormula <- function(metricList) { out_var <- metricList 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 = "") @@ -1255,12 +1257,12 @@ 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_test <- .applyQScoreModel( #' spe=spe_test, -#' qcModel=qc_model, +#' qsModel=qs_model, #' scoreName="QScore_transferred" #' ) #' @@ -1386,3 +1388,61 @@ checkOutliers <- function(spe, verbose=FALSE) { return(ok) } + +## ---- Deprecated functions ----------------------------------------------- + +#' computeQCScore (deprecated) +#' @name computeQCScore +#' @rdname computeQCScore-deprecated +#' @description +#' **Deprecated.** Use \code{\link{computeQScore}} instead. +#' +#' \lifecycle{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." + ) + ) + computeQScore(spe, bestLambda=bestLambda, modelFormula=modelFormula, + verbose=verbose) +} + +#' computeQCScoreFlags (deprecated) +#' @name computeQCScoreFlags +#' @rdname computeQCScoreFlags-deprecated +#' @description +#' **Deprecated.** Use \code{\link{computeQScoreFlags}} instead. +#' +#' \lifecycle{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." + ) + ) + computeQScoreFlags(spe, qsThreshold=qsThreshold, + useQSQuantiles=useQSQuantiles) +} diff --git a/man/computeLambda.Rd b/man/computeLambda.Rd index e7becff..bf293a9 100644 --- a/man/computeLambda.Rd +++ b/man/computeLambda.Rd @@ -7,14 +7,14 @@ computeLambda(modelMatrix, trainDF) } \arguments{ +\item{modelMatrix}{`matrix` +The design matrix built from the training data, typically via +`model.matrix(as.formula(model_formula), data=trainDF)`.} + \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.} - -\item{modelFormula}{`character` -A character string representing the model formula - `~ log2SignalDensity + ...`, as returned by `getModelFormula()`.} + `QScore_train` A binary (0/1) response vector to be modeled.} } \value{ `numeric` @@ -31,14 +31,14 @@ k-fold cross-validation to identify the optimal regularization parameter } \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`. } \examples{ example(computeTrainDF) -modform <- getModelFormula(metadata(spe)$formula_variables) -best_lambda <- computeLambda(df_train, modform) +modform <- getModelFormula(names(metadata(spe)$formula_variables)) +model_matrix <- model.matrix(as.formula(modform), data=df_train) +best_lambda <- computeLambda(model_matrix, df_train) print(best_lambda) diff --git a/man/computeQCScore-deprecated.Rd b/man/computeQCScore-deprecated.Rd new file mode 100644 index 0000000..b844724 --- /dev/null +++ b/man/computeQCScore-deprecated.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/QC.R +\name{computeQCScore} +\alias{computeQCScore} +\title{computeQCScore (deprecated)} +\usage{ +computeQCScore(spe, bestLambda = NULL, modelFormula = NULL, verbose = FALSE) +} +\arguments{ +\item{spe}{A `SpatialExperiment` object.} + +\item{bestLambda}{Passed to \code{\link{computeQScore}}.} + +\item{modelFormula}{Passed to \code{\link{computeQScore}}.} + +\item{verbose}{Passed to \code{\link{computeQScore}}.} +} +\value{ +A `SpatialExperiment` object; see \code{\link{computeQScore}}. +} +\description{ +\strong{Deprecated.} Use \code{\link{computeQScore}} instead. +} diff --git a/man/computeQCScoreFlags-deprecated.Rd b/man/computeQCScoreFlags-deprecated.Rd new file mode 100644 index 0000000..42fcec2 --- /dev/null +++ b/man/computeQCScoreFlags-deprecated.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/QC.R +\name{computeQCScoreFlags} +\alias{computeQCScoreFlags} +\title{computeQCScoreFlags (deprecated)} +\usage{ +computeQCScoreFlags(spe, qsThreshold = 0.5, useQSQuantiles = FALSE) +} +\arguments{ +\item{spe}{A `SpatialExperiment` object.} + +\item{qsThreshold}{Passed to \code{\link{computeQScoreFlags}}.} + +\item{useQSQuantiles}{Passed to \code{\link{computeQScoreFlags}}.} +} +\value{ +A `SpatialExperiment` object; see \code{\link{computeQScoreFlags}}. +} +\description{ +\strong{Deprecated.} Use \code{\link{computeQScoreFlags}} instead. +} diff --git a/man/computeQScore.Rd b/man/computeQScore.Rd index 730e87b..28dc9d5 100644 --- a/man/computeQScore.Rd +++ b/man/computeQScore.Rd @@ -67,9 +67,9 @@ be computed internally, just set a seed with `set.seed()` before running 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 -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. +removed (or added) terms should be written as in the default formula, +e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`. +Whitespace around operators is accepted. In any case, metrics with insufficient outliers (less than 0.1\% of the dataset) will be excluded from the QS formula. diff --git a/man/dot-applyQScoreModel.Rd b/man/dot-applyQScoreModel.Rd index b126f31..4383d16 100644 --- a/man/dot-applyQScoreModel.Rd +++ b/man/dot-applyQScoreModel.Rd @@ -52,12 +52,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_test <- .applyQScoreModel( spe=spe_test, - qcModel=qc_model, + qsModel=qs_model, scoreName="QScore_transferred" ) diff --git a/man/getModelFormula.Rd b/man/getModelFormula.Rd index 49e67fd..3528f4c 100644 --- a/man/getModelFormula.Rd +++ b/man/getModelFormula.Rd @@ -7,20 +7,19 @@ getModelFormula(metricList) } \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{metricList}{A character vector of metric names to include in the +formula (e.g. `"log2SignalDensity"`, `"Area_um"`, etc.), typically the +names of `metadata(spe)$formula_variables`.} } \value{ `character` - A one‐sided formula as a string (e.g. "~ log2SignalDensity + ..."). + 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) -getModelFormula(metadata(spe)$formula_variables) +getModelFormula(names(metadata(spe)$formula_variables)) } diff --git a/tests/testthat/test_QCScores.R b/tests/testthat/test_QCScores.R index 9e59e9a..221b8fa 100644 --- a/tests/testthat/test_QCScores.R +++ b/tests/testthat/test_QCScores.R @@ -6,14 +6,17 @@ spe0 <- example(readCosmxSPE)$value 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")) }) -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,14 +26,25 @@ 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 (deprecated) still works and produces QScore", { + spe <- spatialPerCellQC(spe0) + set.seed(42) + expect_warning( + spe2 <- computeQCScore(spe), + "deprecated" + ) + expect_true("QScore" %in% colnames(colData(spe2))) }) @@ -44,9 +58,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) diff --git a/vignettes/SpaceTrooper_utilities.Rmd b/vignettes/SpaceTrooper_utilities.Rmd index 8de8a96..eac7cb2 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 @@ -506,11 +503,10 @@ 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. + should be written as they appear here (e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`). + 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. From b0b831a6f6058a2bbaef52a16a775d895b803add Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 24 Jul 2026 10:43:37 +0000 Subject: [PATCH 03/13] Fix R CMD check WARNING and NOTE: remove unused csize/calpha args, add stats::terms import --- NAMESPACE | 1 + R/QC.R | 2 +- R/spatialQCPlots.R | 1 - 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index ccdb3b0..0a0fa4d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -134,3 +134,4 @@ importFrom(stats,complete.cases) importFrom(stats,model.matrix) importFrom(stats,predict) importFrom(stats,quantile) +importFrom(stats,terms) diff --git a/R/QC.R b/R/QC.R index c54fed4..5508612 100644 --- a/R/QC.R +++ b/R/QC.R @@ -457,7 +457,7 @@ computeLambda <- function(modelMatrix, trainDF) { #' @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) diff --git a/R/spatialQCPlots.R b/R/spatialQCPlots.R index c730dac..71cdf75 100644 --- a/R/spatialQCPlots.R +++ b/R/spatialQCPlots.R @@ -478,7 +478,6 @@ plotZoomFovsMap <- function(spe, fovs=NULL, title=NULL, mapPointSize=0.5, mapPointAlpha=0.8, fovNumbersCol="black", fovNumberSize=1, fovNumbersAlpha=0.8, - csize=0.05, calpha=0.8, scaleBars=NULL, scaleBarMap=TRUE, scaleBarPol=TRUE, From 7abcb1083224f82f06d006a9f9dffbac3d0047e0 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 24 Jul 2026 11:05:17 +0000 Subject: [PATCH 04/13] fix stale plotZoomFovsMap Rd usage --- man/plotZoomFovsMap.Rd | 2 -- 1 file changed, 2 deletions(-) diff --git a/man/plotZoomFovsMap.Rd b/man/plotZoomFovsMap.Rd index d19a97d..453c7f2 100644 --- a/man/plotZoomFovsMap.Rd +++ b/man/plotZoomFovsMap.Rd @@ -14,8 +14,6 @@ plotZoomFovsMap( fovNumbersCol = "black", fovNumberSize = 1, fovNumbersAlpha = 0.8, - csize = 0.05, - calpha = 0.8, scaleBars = NULL, scaleBarMap = TRUE, scaleBarPol = TRUE, From b6a5134150bf5e4be824353577b7570fa4c3af74 Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 17:32:43 +0200 Subject: [PATCH 05/13] Restore deprecated QScore APIs and legacy outputs --- NAMESPACE | 4 +- R/QC.R | 100 ++++++++++++++++++++++++++++----- tests/testthat/test_QCScores.R | 36 ++++++++++-- 3 files changed, 122 insertions(+), 18 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 0a0fa4d..dc68eed 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) @@ -13,6 +14,7 @@ export(computeLambda) export(computeMissingMetricsMerfish) export(computeMissingMetricsXenium) export(computeOutliersQScore) +export(computeOutliersQCScore) export(computeQCScore) export(computeQCScoreFlags) export(computeQScore) diff --git a/R/QC.R b/R/QC.R index 5508612..212d511 100644 --- a/R/QC.R +++ b/R/QC.R @@ -1211,12 +1211,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 @@ -1260,13 +1261,17 @@ checkOutliers <- function(spe, verbose=FALSE) { #' qs_model <- metadata(spe_train)$QScore_model #' #' ## Apply the trained model to another dataset -#' spe_test <- .applyQScoreModel( +#' spe_test <- applyQScoreModel( #' spe=spe_test, #' 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")) @@ -1397,8 +1402,6 @@ checkOutliers <- function(spe, verbose=FALSE) { #' @description #' **Deprecated.** Use \code{\link{computeQScore}} instead. #' -#' \lifecycle{deprecated} -#' #' @param spe A `SpatialExperiment` object. #' @param bestLambda Passed to \code{\link{computeQScore}}. #' @param modelFormula Passed to \code{\link{computeQScore}}. @@ -1416,8 +1419,17 @@ computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, "See help('computeQScore') for details." ) ) - computeQScore(spe, bestLambda=bestLambda, modelFormula=modelFormula, - verbose=verbose) + 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) @@ -1426,8 +1438,6 @@ computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, #' @description #' **Deprecated.** Use \code{\link{computeQScoreFlags}} instead. #' -#' \lifecycle{deprecated} -#' #' @param spe A `SpatialExperiment` object. #' @param qsThreshold Passed to \code{\link{computeQScoreFlags}}. #' @param useQSQuantiles Passed to \code{\link{computeQScoreFlags}}. @@ -1443,6 +1453,70 @@ computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { "See help('computeQScoreFlags') for details." ) ) - computeQScoreFlags(spe, qsThreshold=qsThreshold, - useQSQuantiles=useQSQuantiles) + 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) +#' @name computeOutliersQCScore +#' @rdname SpaceTrooper-deprecated +#' @description +#' **Deprecated.** Use \code{\link{computeOutliersQScore}} instead. +#' +#' @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) +#' @name applyQCScoreModel +#' @rdname SpaceTrooper-deprecated +#' @description +#' **Deprecated.** Use \code{\link{applyQScoreModel}} instead. +#' +#' @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/tests/testthat/test_QCScores.R b/tests/testthat/test_QCScores.R index 221b8fa..66e9a4e 100644 --- a/tests/testthat/test_QCScores.R +++ b/tests/testthat/test_QCScores.R @@ -1,8 +1,14 @@ 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")) @@ -13,6 +19,9 @@ test_that("QC functions are exported", { ## 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")) }) @@ -37,16 +46,35 @@ test_that("computeQScore adds a QScore between 0 and 1", { expect_true(all(fs[!is.na(fs)] >= 0 & fs[!is.na(fs)] <= 1)) }) -test_that("computeQCScore (deprecated) still works and produces QScore", { +test_that("computeQCScore preserves legacy score and model names", { spe <- spatialPerCellQC(spe0) set.seed(42) expect_warning( spe2 <- computeQCScore(spe), "deprecated" ) - expect_true("QScore" %in% colnames(colData(spe2))) + 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) From 685e4f344db0679fae3e987846a9b8450a7986d9 Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 17:39:05 +0200 Subject: [PATCH 06/13] Restore published argument compatibility --- R/QC.R | 38 +++++++++-- R/spatialQCPlots.R | 102 ++++++++++++++++++++++------ tests/testthat/test_QCScores.R | 117 +++++++++++++++++++++++++++++++++ 3 files changed, 232 insertions(+), 25 deletions(-) diff --git a/R/QC.R b/R/QC.R index 212d511..7d41c1e 100644 --- a/R/QC.R +++ b/R/QC.R @@ -757,18 +757,39 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) { #' @description #' Returns the right‐hand side of a model formula string based on a vector of #' metric names. -#' @param metricList A character vector of metric names to include in the -#' formula (e.g. `"log2SignalDensity"`, `"Area_um"`, etc.), typically the -#' names of `metadata(spe)$formula_variables`. +#' @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 metricList Named replacement for `formulaVars`. If both arguments 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(names(metadata(spe)$formula_variables)) -getModelFormula <- function(metricList) +#' getModelFormula(metadata(spe)$formula_variables) +getModelFormula <- function(formulaVars, verbose=FALSE, metricList) { - out_var <- metricList + has_formula_vars <- !missing(formulaVars) + has_metric_list <- !missing(metricList) + 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) 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))" @@ -776,6 +797,11 @@ getModelFormula <- function(metricList) 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) } diff --git a/R/spatialQCPlots.R b/R/spatialQCPlots.R index 71cdf75..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,16 +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/tests/testthat/test_QCScores.R b/tests/testthat/test_QCScores.R index 66e9a4e..51880c8 100644 --- a/tests/testthat/test_QCScores.R +++ b/tests/testthat/test_QCScores.R @@ -104,3 +104,120 @@ test_that("computeQScoreFlags 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'" + ) +}) From 3c59e86c8bcc48db1a3f8fdecb13364e0e5a85bd Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 17:50:39 +0200 Subject: [PATCH 07/13] Restore public computeLambda contract --- R/QC.R | 80 +++++++++++++++---- tests/testthat/test_computeLambda.R | 119 ++++++++++++++++++++++++++++ 2 files changed, 184 insertions(+), 15 deletions(-) create mode 100644 tests/testthat/test_computeLambda.R diff --git a/R/QC.R b/R/QC.R index 7d41c1e..e993b09 100644 --- a/R/QC.R +++ b/R/QC.R @@ -349,13 +349,13 @@ computeThresholdFlags <- function(spe, totalThreshold=0, #' k-fold cross-validation to identify the optimal regularization parameter #' \eqn{\lambda} for Quality Score (QS) model training. #' -#' @param modelMatrix `matrix` -#' The design matrix built from the training data, typically via -#' `model.matrix(as.formula(model_formula), data=trainDF)`. #' @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. +#' `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: -#' 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(names(metadata(spe)$formula_variables)) -#' model_matrix <- model.matrix(as.formula(modform), data=df_train) -#' best_lambda <- computeLambda(model_matrix, df_train) +#' 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 @@ -513,7 +536,10 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE 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] @@ -1386,6 +1412,30 @@ applyQScoreModel <- function(spe, qsModel, scoreName="QScore") { 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") { 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) +}) From 56d7a25f90e5838e3ae40734154bc4f11c66b165 Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 17:53:57 +0200 Subject: [PATCH 08/13] Validate and preserve custom QScore formulas --- R/QC.R | 184 ++++++++++++++++++++++++--- tests/testthat/test_modelFormula.R | 196 +++++++++++++++++++++++++++++ 2 files changed, 360 insertions(+), 20 deletions(-) create mode 100644 tests/testthat/test_modelFormula.R diff --git a/R/QC.R b/R/QC.R index e993b09..e6491f4 100644 --- a/R/QC.R +++ b/R/QC.R @@ -493,27 +493,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") + supported_metrics <- .qscoreSupportedPredictors() + metric_list <- intersect( + supported_metrics, + names(colData(spe)) + ) + if (!"log2SignalDensity" %in% metric_list) { + stop( + "'log2SignalDensity' is required to construct Quality Score ", + "training labels." + ) + } + formula_info <- NULL if (!is.null(modelFormula)) { - model_formula <- modelFormula - metricList <- attr(terms(as.formula(modelFormula)), "term.labels") - metricList <- metricList[!grepl(":", metricList, fixed=TRUE)] - ## Whitespace-tolerant detection of the border-effect interaction term - border_pat <- "I\\(abs\\(log2AspectRatio\\)\\s*\\*\\s*as\\.numeric\\(dist_border\\s*<\\s*50\\)\\)" - if (any(grepl(border_pat, metricList))) { - metricList <- gsub(border_pat, "log2AspectRatio", metricList) - } + formula_info <- .validateQScoreFormula( + modelFormula=modelFormula, + dataNames=names(colData(spe)), + technology=metadata(spe)$technology + ) } - stopifnot("Not all required metrics in the colData.\nPlease run spatialPerCellQC first." = all(metricList %in% names(colData(spe)))) - ctx <- .prepQCContext(spe, metricList, verbose) + ## Training-label construction is intentionally based on all available, + ## supported metrics and remains independent of a user-selected fit formula. + ctx <- .prepQCContext(spe, metric_list, verbose) df <- ctx$df; out_var <- ctx$out_var; tech <- ctx$tech 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) { @@ -525,14 +536,17 @@ 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)) { @@ -542,7 +556,7 @@ computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE ) } - 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:") @@ -550,11 +564,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. @@ -831,6 +849,119 @@ getModelFormula <- function(formulaVars, verbose=FALSE, metricList) 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 <- stats::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=" ") + } + list( + formula=model_formula, + text=formula_text, + variables=required_variables, + terms=term_labels + ) +} + #' .computeXenMerTrainSet #' @name dot-computeXenMerTrainSet #' @rdname dot-computeXenMerTrainSet @@ -1056,6 +1187,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) @@ -1440,7 +1584,7 @@ applyQScoreModel <- function(spe, qsModel, scoreName="QScore") { .filterCompleteModelCases <- function(df, modelFormula, response=NULL, context="cells") { - vars <- all.vars(as.formula(modelFormula)) + vars <- stats::all.vars(stats::as.formula(modelFormula)) missing_vars <- setdiff(vars, colnames(df)) diff --git a/tests/testthat/test_modelFormula.R b/tests/testthat/test_modelFormula.R new file mode 100644 index 0000000..91305a5 --- /dev/null +++ b/tests/testthat/test_modelFormula.R @@ -0,0 +1,196 @@ +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" + ) +} + +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 fit formula does not redefine training-label metrics", { + spe <- qscore_formula_spe() + 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(206) + computeQScore( + spe, + modelFormula=~ log2SignalDensity + ) + + expect_identical( + seen$metricList, + 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 + ) +}) + +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("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" + ) +}) From 05bdf95ef5360f58e7a26d70f3acdc9b30821e82 Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 17:57:13 +0200 Subject: [PATCH 09/13] Complete QScore model transfer compatibility --- R/QC.R | 15 ++- tests/testthat/test_modelTransfer.R | 156 ++++++++++++++++++++++++++++ 2 files changed, 168 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/test_modelTransfer.R diff --git a/R/QC.R b/R/QC.R index e6491f4..65cf2e6 100644 --- a/R/QC.R +++ b/R/QC.R @@ -629,7 +629,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) } @@ -792,6 +793,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) } @@ -876,7 +880,7 @@ getModelFormula <- function(formulaVars, verbose=FALSE, metricList) stop("'modelFormula' must be a one-sided formula.") } - required_variables <- stats::all.vars(model_formula) + required_variables <- base::all.vars(model_formula) supported <- .qscoreSupportedPredictors() allowed_variables <- c(supported, "dist_border") unsupported_variables <- setdiff( @@ -1484,6 +1488,11 @@ applyQScoreModel <- function(spe, qsModel, scoreName="QScore") { ) df <- as.data.frame(colData(spe)) + .validateQScoreFormula( + modelFormula=qsModel$model_formula, + dataNames=colnames(df), + technology=metadata(spe)$technology + ) ok <- .filterCompleteModelCases( df=df, @@ -1584,7 +1593,7 @@ applyQScoreModel <- function(spe, qsModel, scoreName="QScore") { .filterCompleteModelCases <- function(df, modelFormula, response=NULL, context="cells") { - vars <- stats::all.vars(stats::as.formula(modelFormula)) + vars <- base::all.vars(stats::as.formula(modelFormula)) missing_vars <- setdiff(vars, colnames(df)) 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) + ) +}) From 173d11aefeb11d84b0b85b1225ada2428a94ef29 Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 18:04:45 +0200 Subject: [PATCH 10/13] Document compatible QScore APIs and behavior --- NAMESPACE | 2 +- NEWS.md | 33 +++++- R/QC.R | 100 +++++++++++------- README.md | 18 +++- man/SpaceTrooper-deprecated.Rd | 74 +++++++++++++ ...pplyQScoreModel.Rd => applyQScoreModel.Rd} | 13 +-- man/computeLambda.Rd | 22 ++-- man/computeQCScore-deprecated.Rd | 23 ---- man/computeQCScoreFlags-deprecated.Rd | 21 ---- man/computeQScore.Rd | 32 +++--- man/computeTrainDF.Rd | 3 +- man/getModelFormula.Rd | 21 ++-- man/plotCellsFovs.Rd | 36 ++++--- man/plotZoomFovsMap.Rd | 41 ++++--- man/trainModel.Rd | 4 +- vignettes/SpaceTrooper_utilities.Rmd | 48 ++++++--- 16 files changed, 322 insertions(+), 169 deletions(-) create mode 100644 man/SpaceTrooper-deprecated.Rd rename man/{dot-applyQScoreModel.Rd => applyQScoreModel.Rd} (89%) delete mode 100644 man/computeQCScore-deprecated.Rd delete mode 100644 man/computeQCScoreFlags-deprecated.Rd diff --git a/NAMESPACE b/NAMESPACE index dc68eed..a5e4a18 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,8 +13,8 @@ export(computeCenterFromPolygons) export(computeLambda) export(computeMissingMetricsMerfish) export(computeMissingMetricsXenium) -export(computeOutliersQScore) export(computeOutliersQCScore) +export(computeOutliersQScore) export(computeQCScore) export(computeQCScoreFlags) export(computeQScore) diff --git a/NEWS.md b/NEWS.md index f7aff73..9cd02e4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,35 @@ # 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 +* 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. +* 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 65cf2e6..a3aee26 100644 --- a/R/QC.R +++ b/R/QC.R @@ -437,9 +437,9 @@ computeLambda <- function(trainDF, modelFormula) { #' 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 @@ -451,19 +451,19 @@ computeLambda <- function(trainDF, modelFormula) { #' `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 as in the default formula, -#' e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`. -#' Whitespace around operators is accepted. -#' 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 may use any subset of +#' the four supported predictors: `log2SignalDensity`, `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))`. Formula selection +#' is independent of the supported metrics used to construct training labels. #' #' 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 @@ -472,10 +472,10 @@ computeLambda <- function(trainDF, modelFormula) { #' @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 @@ -610,8 +610,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. #' @@ -664,7 +664,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. #' @@ -810,18 +811,28 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) { #' `metadata(spe)$formula_variables`. An unnamed character vector of metric #' names is also accepted. #' @param verbose Logical. If `TRUE`, prints the generated formula. -#' @param metricList Named replacement for `formulaVars`. If both arguments are -#' supplied, `metricList` takes precedence with a warning. +#' @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(formulaVars, verbose=FALSE, metricList) +getModelFormula <- function(formulaVars, verbose=FALSE, ...) { has_formula_vars <- !missing(formulaVars) - has_metric_list <- !missing(metricList) + 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.") } @@ -831,7 +842,7 @@ getModelFormula <- function(formulaVars, verbose=FALSE, metricList) "using 'metricList'." ) } - out_var <- if (has_metric_list) metricList else formulaVars + 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.") } @@ -1625,11 +1636,33 @@ applyQScoreModel <- function(spe, qsModel, scoreName="QScore") { ## ---- Deprecated functions ----------------------------------------------- -#' computeQCScore (deprecated) -#' @name computeQCScore -#' @rdname computeQCScore-deprecated +#' Deprecated functions in SpaceTrooper +#' +#' @name SpaceTrooper-deprecated #' @description -#' **Deprecated.** Use \code{\link{computeQScore}} instead. +#' 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}}. @@ -1662,10 +1695,7 @@ computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, } #' computeQCScoreFlags (deprecated) -#' @name computeQCScoreFlags -#' @rdname computeQCScoreFlags-deprecated -#' @description -#' **Deprecated.** Use \code{\link{computeQScoreFlags}} instead. +#' @rdname SpaceTrooper-deprecated #' #' @param spe A `SpatialExperiment` object. #' @param qsThreshold Passed to \code{\link{computeQScoreFlags}}. @@ -1700,10 +1730,7 @@ computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { } #' computeOutliersQCScore (deprecated) -#' @name computeOutliersQCScore #' @rdname SpaceTrooper-deprecated -#' @description -#' **Deprecated.** Use \code{\link{computeOutliersQScore}} instead. #' #' @param spe A `SpatialExperiment` object. #' @param metricList Passed to \code{\link{computeOutliersQScore}}. @@ -1722,10 +1749,7 @@ computeOutliersQCScore <- function(spe, metricList=c( } #' applyQCScoreModel (deprecated) -#' @name applyQCScoreModel #' @rdname SpaceTrooper-deprecated -#' @description -#' **Deprecated.** Use \code{\link{applyQScoreModel}} instead. #' #' @param spe A `SpatialExperiment` object with QC metrics already computed. #' @param qcModel A historical QC score model object, usually stored in 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 89% rename from man/dot-applyQScoreModel.Rd rename to man/applyQScoreModel.Rd index 4383d16..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. } @@ -55,7 +56,7 @@ spe_train <- computeQScore(spe_train) qs_model <- metadata(spe_train)$QScore_model ## Apply the trained model to another dataset -spe_test <- .applyQScoreModel( +spe_test <- applyQScoreModel( spe=spe_test, qsModel=qs_model, scoreName="QScore_transferred" diff --git a/man/computeLambda.Rd b/man/computeLambda.Rd index bf293a9..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{modelMatrix}{`matrix` -The design matrix built from the training data, typically via -`model.matrix(as.formula(model_formula), data=trainDF)`.} - \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}{A one-sided formula or character string describing the +supported predictor terms.} } \value{ `numeric` @@ -30,15 +30,15 @@ k-fold cross-validation to identify the optimal regularization parameter \eqn{\lambda} for Quality Score (QS) model training. } \details{ -Internally, the function: - 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(names(metadata(spe)$formula_variables)) -model_matrix <- model.matrix(as.formula(modform), data=df_train) -best_lambda <- computeLambda(model_matrix, df_train) +best_lambda <- computeLambda(df_train, modform) print(best_lambda) diff --git a/man/computeQCScore-deprecated.Rd b/man/computeQCScore-deprecated.Rd deleted file mode 100644 index b844724..0000000 --- a/man/computeQCScore-deprecated.Rd +++ /dev/null @@ -1,23 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/QC.R -\name{computeQCScore} -\alias{computeQCScore} -\title{computeQCScore (deprecated)} -\usage{ -computeQCScore(spe, bestLambda = NULL, modelFormula = NULL, verbose = FALSE) -} -\arguments{ -\item{spe}{A `SpatialExperiment` object.} - -\item{bestLambda}{Passed to \code{\link{computeQScore}}.} - -\item{modelFormula}{Passed to \code{\link{computeQScore}}.} - -\item{verbose}{Passed to \code{\link{computeQScore}}.} -} -\value{ -A `SpatialExperiment` object; see \code{\link{computeQScore}}. -} -\description{ -\strong{Deprecated.} Use \code{\link{computeQScore}} instead. -} diff --git a/man/computeQCScoreFlags-deprecated.Rd b/man/computeQCScoreFlags-deprecated.Rd deleted file mode 100644 index 42fcec2..0000000 --- a/man/computeQCScoreFlags-deprecated.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/QC.R -\name{computeQCScoreFlags} -\alias{computeQCScoreFlags} -\title{computeQCScoreFlags (deprecated)} -\usage{ -computeQCScoreFlags(spe, qsThreshold = 0.5, useQSQuantiles = FALSE) -} -\arguments{ -\item{spe}{A `SpatialExperiment` object.} - -\item{qsThreshold}{Passed to \code{\link{computeQScoreFlags}}.} - -\item{useQSQuantiles}{Passed to \code{\link{computeQScoreFlags}}.} -} -\value{ -A `SpatialExperiment` object; see \code{\link{computeQScoreFlags}}. -} -\description{ -\strong{Deprecated.} Use \code{\link{computeQScoreFlags}} instead. -} diff --git a/man/computeQScore.Rd b/man/computeQScore.Rd index 28dc9d5..3cbdbae 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,19 @@ 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 as in the default formula, -e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`. -Whitespace around operators is accepted. -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 may use any subset of +the four supported predictors: `log2SignalDensity`, `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))`. Formula selection +is independent of the supported metrics used to construct training labels. 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 3528f4c..f312fae 100644 --- a/man/getModelFormula.Rd +++ b/man/getModelFormula.Rd @@ -4,22 +4,29 @@ \alias{getModelFormula} \title{getModelFormula} \usage{ -getModelFormula(metricList) +getModelFormula(formulaVars, verbose = FALSE, ...) } \arguments{ -\item{metricList}{A character vector of metric names to include in the -formula (e.g. `"log2SignalDensity"`, `"Area_um"`, etc.), typically the -names of `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 + ..."). + A one‐sided formula as a string (e.g. "~ log2SignalDensity + ..."). } \description{ -Returns the right-hand side of a model formula string based on a vector of +Returns the right‐hand side of a model formula string based on a vector of metric names. } \examples{ example(checkOutliers) -getModelFormula(names(metadata(spe)$formula_variables)) +getModelFormula(metadata(spe)$formula_variables) } 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 453c7f2..cc32161 100644 --- a/man/plotZoomFovsMap.Rd +++ b/man/plotZoomFovsMap.Rd @@ -9,15 +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{ @@ -33,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. @@ -58,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/vignettes/SpaceTrooper_utilities.Rmd b/vignettes/SpaceTrooper_utilities.Rmd index eac7cb2..dfde8d3 100644 --- a/vignettes/SpaceTrooper_utilities.Rmd +++ b/vignettes/SpaceTrooper_utilities.Rmd @@ -493,29 +493,32 @@ 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 - should be written as they appear here (e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))`). - 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. Users may select a +supported subset, 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 metrics used to +construct good and bad training labels remain governed by the established +SpaceTrooper rules and are not redefined by omitting a predictor from the fitted +formula. ```{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) ``` @@ -529,6 +532,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 From 9ffbaac6ec3728d9081a8001161d443939fb78f7 Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 18:13:35 +0200 Subject: [PATCH 11/13] Address QScore BiocCheck findings --- R/QC.R | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/R/QC.R b/R/QC.R index a3aee26..c60066b 100644 --- a/R/QC.R +++ b/R/QC.R @@ -1331,8 +1331,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]) } } } @@ -1403,8 +1404,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))) From 065d407d0770178672ea89cce56824567603d82c Mon Sep 17 00:00:00 2001 From: Dario Date: Wed, 29 Jul 2026 18:18:03 +0200 Subject: [PATCH 12/13] Bump development version to 1.1.9 --- DESCRIPTION | 4 ++-- NEWS.md | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) 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/NEWS.md b/NEWS.md index 9cd02e4..4a17db2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# Changes in version 1.1.8 +# Changes in version 1.1.9 * Added the canonical `computeQScore()`, `computeQScoreFlags()`, `computeOutliersQScore()`, and `applyQScoreModel()` APIs. From a61d732a45d3b11b2f15feaf1eb25f3f4f3c71b7 Mon Sep 17 00:00:00 2001 From: Dario Date: Fri, 31 Jul 2026 11:47:27 +0200 Subject: [PATCH 13/13] Align QScore training labels with custom formulas --- NEWS.md | 3 + R/QC.R | 52 +++++---- man/computeQScore.Rd | 17 ++- tests/testthat/test_modelFormula.R | 156 ++++++++++++++++++++++++--- vignettes/SpaceTrooper_utilities.Rmd | 19 ++-- 5 files changed, 203 insertions(+), 44 deletions(-) diff --git a/NEWS.md b/NEWS.md index 4a17db2..2e92b31 100644 --- a/NEWS.md +++ b/NEWS.md @@ -18,6 +18,9 @@ * 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 diff --git a/R/QC.R b/R/QC.R index c60066b..e63cde1 100644 --- a/R/QC.R +++ b/R/QC.R @@ -454,12 +454,21 @@ computeLambda <- function(trainDF, modelFormula) { #' - Model formula details: the most complete default formula is: #' `~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2`. #' A user-supplied formula is fitted as supplied: additive structures, selected -#' interactions, `:` and `*` are preserved. The formula may use any subset of -#' the four supported predictors: `log2SignalDensity`, `Area_um`, +#' 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))`. Formula selection -#' is independent of the supported metrics used to construct training labels. +#' `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 @@ -493,29 +502,28 @@ 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] } - supported_metrics <- .qscoreSupportedPredictors() - metric_list <- intersect( - supported_metrics, - names(colData(spe)) - ) - if (!"log2SignalDensity" %in% metric_list) { - stop( - "'log2SignalDensity' is required to construct Quality Score ", - "training labels." - ) - } - formula_info <- NULL - if (!is.null(modelFormula)) { + 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." + ) } - ## Training-label construction is intentionally based on all available, - ## supported metrics and remains independent of a user-selected fit formula. ctx <- .prepQCContext(spe, metric_list, verbose) df <- ctx$df; out_var <- ctx$out_var; tech <- ctx$tech @@ -969,11 +977,15 @@ getModelFormula <- function(formulaVars, verbose=FALSE, ...) } 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 + terms=term_labels, + training_metrics=training_metrics ) } diff --git a/man/computeQScore.Rd b/man/computeQScore.Rd index 3cbdbae..8b5ceac 100644 --- a/man/computeQScore.Rd +++ b/man/computeQScore.Rd @@ -66,12 +66,21 @@ be computed internally, just set a seed with `set.seed()` before running - Model formula details: the most complete default formula is: `~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2`. A user-supplied formula is fitted as supplied: additive structures, selected -interactions, `:` and `*` are preserved. The formula may use any subset of -the four supported predictors: `log2SignalDensity`, `Area_um`, +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))`. Formula selection -is independent of the supported metrics used to construct training labels. +`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 diff --git a/tests/testthat/test_modelFormula.R b/tests/testthat/test_modelFormula.R index 91305a5..a07aa62 100644 --- a/tests/testthat/test_modelFormula.R +++ b/tests/testthat/test_modelFormula.R @@ -21,6 +21,25 @@ model_terms <- function(spe) { ) } +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(), @@ -106,26 +125,62 @@ test_that("custom formula is not rebuilt by getModelFormula", { )) }) -test_that("custom fit formula does not redefine training-label metrics", { +test_that("custom formulas define metrics passed to QScore preparation", { spe <- qscore_formula_spe() - 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" + + single <- captured_training_metrics( + spe, + ~ log2SignalDensity, + seed=206 ) + expect_identical(single$metrics, "log2SignalDensity") - set.seed(206) - computeQScore( + additive <- captured_training_metrics( spe, - modelFormula=~ log2SignalDensity + ~ 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( - seen$metricList, + 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() ) }) @@ -143,6 +198,72 @@ test_that("supported CosMx border expression is accepted", { 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", { @@ -173,6 +294,15 @@ test_that("unsupported predictors and transformations are rejected", { ) }) +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 diff --git a/vignettes/SpaceTrooper_utilities.Rmd b/vignettes/SpaceTrooper_utilities.Rmd index dfde8d3..b9e18b1 100644 --- a/vignettes/SpaceTrooper_utilities.Rmd +++ b/vignettes/SpaceTrooper_utilities.Rmd @@ -505,14 +505,19 @@ 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 is fitted as supplied. Users may select a -supported subset, use an additive formula, or choose interactions with `:` and -`*`; SpaceTrooper does not add or remove terms. Predictor and transformation +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 metrics used to -construct good and bad training labels remain governed by the established -SpaceTrooper rules and are not redefined by omitting a predictor from the fitted -formula. +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}