diff --git a/NAMESPACE b/NAMESPACE index 5911763..0a0fa4d 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) @@ -132,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 7bd3e76..5508612 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. #' @@ -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) @@ -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/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, 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/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, 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.