diff --git a/DESCRIPTION b/DESCRIPTION index 75e58a1..85f3527 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.7 +Version: 1.1.8 Authors@R: c(person("Dario", "Righelli", email="dario.righelli@gmail.com", diff --git a/NAMESPACE b/NAMESPACE index cc92c05..5911763 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,10 +1,10 @@ # Generated by roxygen2: do not edit by hand +export(.applyQScoreModel) export(.getActiveGeometryName) export(.renameGeometry) export(.setActiveGeometry) export(addPolygonsToSPE) -export(applyQCScoreModel) export(checkOutliers) export(computeAreaFromPolygons) export(computeAspectRatioFromPolygons) @@ -12,9 +12,9 @@ export(computeCenterFromPolygons) export(computeLambda) export(computeMissingMetricsMerfish) export(computeMissingMetricsXenium) -export(computeOutliersQCScore) -export(computeQCScore) -export(computeQCScoreFlags) +export(computeOutliersQScore) +export(computeQScore) +export(computeQScoreFlags) export(computeSpatialOutlier) export(computeThresholdFlags) export(computeTrainDF) diff --git a/NEWS.md b/NEWS.md index 0887798..f7aff73 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,7 @@ # Changes in version 1.1.8 -* fixing naming of of spacetrooper utilities vignette -* adding functions for QC model transfer across datasets +* implementing Quality Score computation with custom formula +* fixing naming of spacetrooper utilities vignette * adding citation file with biorxiv paper # Changes in version 1.1.7 diff --git a/R/QC.R b/R/QC.R index a265b95..7bd3e76 100644 --- a/R/QC.R +++ b/R/QC.R @@ -4,10 +4,11 @@ #' @description #' Computes quality‐control metrics for each cell and adds them to `colData`. #' -#' @param spe A `SpatialExperiment` object containing spatial data. +#' @param spe A `SpatialExperiment` object containing spatial omics data. #' @param micronConvFact Numeric factor to convert pixels to microns. Default -#' `0.12`. -#' @param rmZeros logical for removing zero counts cells (default is TRUE). +#' is 0.12. +#' @param rmZeros logical for removing zero total counts (protein intensities) +#' cells (default is TRUE). #' @param negProbList Character vector of patterns to identify negative probes. #' Defaults include: #' Nanostring CosMx: `"NegPrb"`, `"Negative"`, `"SystemControl"` @@ -26,47 +27,57 @@ #' #' Key behaviours and expectations: #' - Feature detection: negative‑probe patterns supplied in `negProbList` are -#' used to build `subsets_*` groups passed to `scater::addPerCellQC()` -#' (via `use_altexps` when requested). `addPerCellQC()` must be able to find -#' matching feature names (`rownames` of spe) and will create `subsets_*_sum` -#' and `subsets_*_detected` columns used below. -#' - Required columns: the function expects `sum`, `detected` and `total` -#' (from `addPerCellQC()` and the SPE assays) to be present; these are used to -#' compute `control_sum`, `control_detected`, `target_sum` and -#' `target_detected.` +#' passed to `scater::addPerCellQC()` (via `use_altexps` when requested). +#' `addPerCellQC()` must be able to find matching control feature names +#' (`rownames` of spe) and will create `subsets_*_sum` and `subsets_*_detected` +#' columns used below for each of control feature types. +#' +#' The function appends `sum`, `detected` and `total` to the SPE colData by +#' calling internally `addPerCellQC()`; these are used to compute `control_sum`, +#' `control_detected`, `target_sum` and `target_detected.` +#' #' - Control metrics: `control_sum` and `control_detected` are computed by #' summing matching `subsets_*` columns; `target_*` metrics are computed as -#' the complement vs sum / detected. -#' - Ratios and logs: `ctrl_total_ratio` (`control` / 9) and its stabilized -#' log2 transform `log2Ctrl_total_ratio` are added. +#' sum(detected) - control_sum(_detected). +#' +#' - Ratios and logs: `ctrl_total_ratio` (`control_sum`/`sum`) and its stabilized +#' log2-transformed `log2Ctrl_total_ratio` are added. +#' #' - Coordinate and area handling: #' - For CosMx technologies (Nanostring_CosMx and Nanostring_CosMx_Protein) #' spatial coordinates are converted from pixels to microns using #' `micronConvFact` and appended to `colData` (column names have px -> um). +#' #' - For CosMx, `Area_um` is derived from an existing Area column scaled by #' micronConvFact^2 and `.computeBorderDistanceCosMx()` is invoked to #' compute `dist_border`. +#' #' - For Nanostring_CosMx_Protein, a legacy `Area.um2` column (if present) #' is renamed to `Area_um` to standardize naming. +#' #' - For Xenium (10X_Xenium), if `Area_um` is missing the function will #' attempt to use `cell_area` as a fallback and issue a warning. +#' #' - Aspect ratio: if `AspectRatio` exists it is logged (`log2AspectRatio`); if #' missing a warning is emitted. +#' #' - Signal density: `SignalDensity` is computed as `sum` / `Area_um` for most -#' technologies; for Nanostring_CosMx_Protein it is set to total. A -#' log transform `log2SignalDensity` is also added. -#' - Zero‑count removal: when `rmZeros = TRUE` cells with `sum == 0` are -#' removed from the returned SpatialExperiment (message printed). +#' technologies; for Nanostring_CosMx_Protein it coincides with `total`. A +#' log2-transformed `log2SignalDensity` is also added. +#' +#' - Zero‑count (protein intensity) removal: when `rmZeros = TRUE` cells with +#' `sum == 0` are removed from the returned SpatialExperiment (message printed). +#' #' - Side effects: the function modifies `colData(spe)` (adds multiple new -#' columns), may add spatial coordinates into `colData` if missing, and may +#' columns). It may add spatial coordinates into `colData` if missing, and may #' subset the SPE to remove zero‑count cells. It issues warnings when expected #' inputs (e.g. area, aspect ratio, polygon‑derived fields) are missing or #' when fallbacks are used. #' #' Use this information to ensure the input SpatialExperiment contains the #' necessary assays and fields (feature names, `sum`/`detected`/`total`, -#' `Area`/`cell_area` when available) so metrics are computed and assigned -#' correctly. +#' `Area`/`cell_area`, log2AspectRatio, dist_border when available) so +#' metrics are computed and assigned correctly. #' #' @importFrom SummarizedExperiment colData #' @importFrom scater addPerCellQC @@ -109,7 +120,7 @@ spatialPerCellQC <- function(spe, micronConvFact=0.12, rmZeros=TRUE, eps <- 0.0001 spe$log2Ctrl_total_ratio <- log2(spe$ctrl_total_ratio+eps) if(metadata(spe)$technology == "Nanostring_CosMx_Protein") { - # Only for proteins will be included in QCScore + # Only for proteins will be included in QScore idx <- which(names(colData(spe)) == "Area.um2") if(length(idx)!=0) { names(colData(spe))[idx] <- "Area_um" } } @@ -155,14 +166,14 @@ spatialPerCellQC <- function(spe, micronConvFact=0.12, rmZeros=TRUE, #' @name .computeBorderDistanceCosMx #' @rdname dot-computeBorderDistanceCosMx #' @description -#' Calculates the minimum distance of each cell to the field‐of‐view border -#' and adds it to `colData`. +#' Calculates the minimum distance of each cell to the Field‐Of‐View (FOV) +#' border and adds it to `colData`. #' #' @param spe A `SpatialExperiment` object with CosMx data. -#' @param xwindim Width of FOV in x (default from `metadata(spe)$fov_dim`). -#' @param ywindim Height of FOV in y (default from `metadata(spe)$fov_dim`). +#' @param xwindim Width of FOV along x axis (default from `metadata(spe)$fov_dim`). +#' @param ywindim Height of FOV along y axis (default from `metadata(spe)$fov_dim`). #' -#' @return A `SpatialExperiment` object with `dist_border` columns in +#' @return A `SpatialExperiment` object with `dist_border` column in #' `colData`. #' #' @importFrom dplyr left_join @@ -191,37 +202,33 @@ spatialPerCellQC <- function(spe, micronConvFact=0.12, rmZeros=TRUE, #' @name computeSpatialOutlier #' @rdname computeSpatialOutlier #' @description -#' Computes outliers based on the Area (in micron) of the experiment. -#' It gives the possibility to choose between the medcouple (mc method argument) -#' and the MADs (scuttle method argument). +#' Computes outliers for quantitative cell metadata based on a statistical test: +#' medcouple (`mc` method argument) or the Median Absolute Deviation (MAD, `scuttle` +#' method argument), when the metric distribution is respectively asymmetric or symmetric. #' #' @details #' The medcouple method is a measure for the skeweness of univariate #' distribution as described in Hubert M. et al. (2008). #' In particular, the computed medcouple value must be in a range between -0.6 -#' and 0.6 to computed adjusted boxplots and perform the outlier detection. -#' For median absolute deviations (MADs) method we just wrap the isOutlier -#' function in the scuttle package. Please see McCarthy DJ et al (2017) -#' for further details. -#' -#' @param spe a SpatialExperiment object with target_counts, area in micron -#' and log2 of the aspect ratio in the `colData`. -#' @param computeBy character indicating a `colData` column name on which -#' compute the outlier. +#' and 0.6 to compute adjusted boxplots and perform the outlier detection. +#' For MAD method we just wrap the isOutlier function in the scuttle package. +#' Please see McCarthy DJ et al (2017) for further details. +#' +#' @param spe a SpatialExperiment object with quantitative cell metrics in the `colData`. +#' @param computeBy character indicating a `colData` column name to compute the outliers on. #' @param method one of `mc`, `scuttle`, `both`. -#' Use `mc` for medcouple, `scuttle` for median absolute deviations as computed -#' in `scuttle`, `both` for computing both of them. -#' @param mcDoScale logical indicating if the values to compute the medcouple -#' for the outlier detection should be scaled (default is FALSE, as suggested -#' by the original Medcouple authors.). See \link[robustbase]{mc} for further -#' readings. +#' Use `mc` for medcouple, `scuttle` for MAD as computed in `scuttle`, `both` +#' for computing both of them. +#' @param mcDoScale logical indicating if the metric values should be scaled +#' to compute medcouple (default is FALSE, as suggested by the original medcouple +#' authors). See \link[robustbase]{mc} for further readings. #' @param scuttleType One of `"both"`, `"lower"`, `"higher"` for scuttle method. #' #' @return a SpatialExperiment object with additional column(s) (named as -#' the column name indicated in `column_by` followed by the outlier_sc/mc -#' nomenclature) with the outlier detection as `outlier.filter` logical class -#' object. This allows to store the thresholds as attributes of the column. -#' use attr(,"thresholds") to retrieve them. +#' the column name indicated in `computeBy` followed by the outlier_sc/mc). +#' Resulting columns are `outlier.filter` logical class objects, allowing to store +#' the thresholds as attributes of the columns. +#' Use attr(,"thresholds") to retrieve them. #' #' @export #' @importFrom robustbase mc adjbox @@ -288,24 +295,25 @@ computeSpatialOutlier <- function(spe, computeBy=NULL, #' @name computeThresholdFlags #' @rdname computeThresholdFlags #' @description -#' Compute Flagged cells using fixed thresholds for SpatialExperiment. +#' Computes flagged cells using fixed thresholds for SpatialExperiment. #' -#' This function calculates flagged cells only for total counts and control on -#' total probe counts ratio using fixed thresholds for a `SpatialExperiment` +#' This function calculates flagged cells only for total counts (protein intensities) +#' and control-to-total probe counts ratio using fixed thresholds for a `SpatialExperiment` #' object. #' -#' @param spe A `SpatialExperiment` object with spatial transcriptomics data. +#' @param spe A `SpatialExperiment` object with spatial omics data. #' @param totalThreshold A numeric value for the threshold of total counts to #' identify cells with low counts. Default is `0`. #' @param ctrlTotRatioThreshold A numeric value for the threshold of #' control-to-total ratio to flag cells over a certain threshold. Default is #' `0.1`. #' -#' @return The `SpatialExperiment` object with added filter flags in `colData`. +#' @return The `SpatialExperiment` object with added flagged cells as TRUE/FALSE +#' in `colData`. #' #' @details The function flags cells basing on zero counts and control-to-total #' ratio to identify junk cells. -#' It also combines these flags into a single filter flag. +#' It also combines these two flags into a single flag metric named `threshold_flags`. #' #' @importFrom SummarizedExperiment colData #' @export @@ -334,15 +342,15 @@ computeThresholdFlags <- function(spe, totalThreshold=0, #' computeLambda #' @description -#' Compute Optimal Ridge Regularization Parameter \eqn{\lambda} via -#' Cross-Validation +#' Computes optimal ridge regularization parameter \eqn{\lambda} via +#' k-fold cross-validation. #' #' \code{computeLambda} performs ridge (L2) logistic regression with -#' cross-validation to identify the optimal regularization parameter -#' \eqn{\lambda} for a binary response. +#' k-fold cross-validation to identify the optimal regularization parameter +#' \eqn{\lambda} for Quality Score (QS) model training. #' #' @param trainDF `data.frame` -#' A data frame for training that must include: +#' 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` @@ -356,9 +364,9 @@ computeThresholdFlags <- function(spe, totalThreshold=0, #' #' @details #' Internally, the function: -#' Constructs the design matrix via \code{model.matrix()}, -#' Runs ridge logistic regression cross-validation using `cv.glmnet` with `alpha = 0`, -#' Extracts and returns `ridge_cv$lambda.min`. +#' 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) @@ -368,92 +376,84 @@ computeThresholdFlags <- function(spe, totalThreshold=0, #' #' #' @export -computeLambda <- function(trainDF, modelFormula) { - # model_formula <- getModelFormula(technology) - train_ok <- .filterCompleteModelCases( - df=trainDF, - modelFormula=modelFormula, - response=trainDF$qcscore_train, - context="training cells for lambda selection" - ) - - trainDF <- trainDF[train_ok, , drop=FALSE] - model_matrix <- model.matrix(as.formula(modelFormula), data=trainDF) - model_matrix <- .dropModelIntercept(model_matrix) - ridge_cv <- cv.glmnet(model_matrix, trainDF$qcscore_train, +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) } -#' computeQCScore -#' @name computeQCScore -#' @rdname computeQCScore +#' computeQScore +#' @name computeQScore +#' @rdname computeQScore #' @description -#' Compute QC score and automatically define weights for QC score -#' through glm training. This function computes QC score with a formula -#' defined on the metrics as "log2SignalDensity", "Area_um", +#' Computes Quality Score (QS) and automatically define weights of the QS model +#' terms through glm training. This function computes QS with a formula +#' defined on the metrics "log2SignalDensity", "Area_um", #' "log2AspectRatio", "log2Ctrl_total_ratio", as computed from #' the `spatialPerCellQC` function. -#' It automatically computes the number of available outliers for each available -#' metric, as they are needed for the glm training. +#' It automatically computes the number of outliers for each available +#' metric, which are needed for the glm training. #' See Details for further information. #' -#' @details For all the techologies, the QC Score formula depends on the follow -#' metrics: +#' @details For all the technologies, the QS formula can contain +#' only the following terms at most: #' -#' QC score ~ count density - aspect ratio - control-total ratio - size +#' Quality Score ~ signal density - size - border effect - background signal #' -#' Where count density is the total counts-to-size ratio, aspect ratio -#' represents the ratio between the width and the height of the cell -#' (computed from the provided polygons if not already present in the experiment -#' metadata) and control-total ratio is the aspecific signal; -#' size is the area for CosMx and Xenium, while it is the volume for Merfish. +#' Where signal density is the total counts/protein intensities-to-size ratio, +#' border effect is represented by the aspect ratio (i.e. width/height of the cell) +#' and is considered only for cells with a distance from the FOV border less than +#' 50 px and only for CosMx datasets. Background signal is represented by the ratio +#' of control probe counts to total counts (protein intensities) and cell size is +#' the area for CosMx and Xenium, while it is the volume for Merfish. #' For each couple of variables interaction terms are computed. #' -#' Additionally, for CosMx datasets, the distance from the border of the FOV is -#' also included in the formula as a metric to take into account. -#' For Xenium and Merscope datasets, QC score cannot depend on FoV border effect, -#' as no FOV border effect was captured through this metric. -#' #' 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 < 0.1% out of the -#' entire dataset, the metric will be excluded from the QC score formula. +#' 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. #' #' - Model fitting: ridge (L2) logistic regression is fitted (via `glmnet`) on #' the balanced training set. The function uses `trainModel()` for fitting -#' and `computeLambda()` (cross‑validation) to select lambda unless -#' `bestLambda` is supplied. +#' and `computeLambda()` to select lambda unless `bestLambda` is supplied. #' #' - Lambda details: because of the randomness in the training set selection, -#' results may vary so that it is possible to set a fixed lambda value -#' previously computed with `computeLambda` preceeded by `computeTrainDF` and -#' `getModelFormula`. -#' This is useful for reproducibility across different runs. -#' Otherwise, an easier way is to let be lambda computed internally, just set -#' a seed with `set.seed()` before running `computeQCScore`. -#' -#' The computed model output is stored in `metadata(spe)$QCScore_model` -#' for inspection and reuse (see also \code{\link{applyQCScoreModel}}). +#' results may vary. The easiest way to ensure reproducibility is to let lambda +#' 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: +#' `~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2`. +#' When user-provided, the formula must follow the same default syntax and +#' removed (or added) terms should be written exactly as in the default formula, +#' e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))` must have spaces +#' around the `*` and `<` operators. +#' In any case, metrics with insufficient outliers (less than 0.1\% of the dataset) +#' will be excluded from the QS formula. +#' +#' 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 +#' 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 +#' details on scalability and coefficient transfer. #' -#' @param spe A `SpatialExperiment` object with spatial transcriptomics data. +#' @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 model formula to be -#' used for training the model. If NULL, the formula is automatically generated -#' based on the available metrics and outliers in the dataset. -#' See Details for more information. -#' Note that the automatically generated formula will include interaction -#' terms between the metrics, and will exclude metrics with insufficient -#' outliers (< 0.1% of the dataset). If a custom `modelFormula` is provided, -#' it will be used as is without modification or checks for outlier counts. -#' -#' @return The `SpatialExperiment` object with added QC score in `colData`. +#' @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. +#' @return The `SpatialExperiment` object with added Quality Score in `colData`. #' @export #' @importFrom dplyr case_when filter mutate distinct pull #' @importFrom glmnet glmnet cv.glmnet @@ -461,9 +461,9 @@ computeLambda <- function(trainDF, modelFormula) { #' @examples #' example(spatialPerCellQC) #' set.seed(1998) -#' spe <- computeQCScore(spe) -#' summary(spe$QC_score) -computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE) { +#' spe <- computeQScore(spe) +#' summary(spe$QScore) +computeQScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALSE) { stopifnot(is(spe, "SpatialExperiment")) if (dim(spe[,spe$total == 0])[2] != 0) { warning(paste0(dim(spe[,spe$total == 0])[2], @@ -471,40 +471,49 @@ computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALS spe <- spe[,spe$total > 0] } metricList <- c("log2SignalDensity", "Area_um", - "log2AspectRatio", "log2Ctrl_total_ratio") + "log2AspectRatio", + "log2Ctrl_total_ratio") + + if (!is.null(modelFormula)) { + model_formula <- modelFormula + metricList <- attr(terms(as.formula(modelFormula)), "term.labels") + metricList <- metricList[!grepl(":", metricList, fixed=TRUE)] + if("I(abs(log2AspectRatio) * as.numeric(dist_border < 50))" %in% metricList) { + metricList <- gsub("I\\(abs\\((log2AspectRatio)\\) \\* as\\.numeric\\((dist_border) < 50\\)\\)", + "log2AspectRatio", metricList) + } + } stopifnot("Not all required metrics in the colData.\nPlease run spatialPerCellQC first." = all(metricList %in% names(colData(spe)))) ctx <- .prepQCContext(spe, metricList, verbose) df <- ctx$df; out_var <- ctx$out_var; tech <- ctx$tech - train_df <- computeTrainDF(df, out_var, tech, verbose) + model_formula <- getModelFormula(names(out_var)) - if (is.null(modelFormula)) { - model_formula <- getModelFormula(out_var, verbose) - } else { - model_formula <- modelFormula - if (verbose) { - message("Using user-supplied model formula:") - message(model_formula) - } + if (verbose) { + message("Using final model formula:") + message(model_formula) } + + train_df <- computeTrainDF(df, out_var, tech, verbose) + train_ok <- .filterCompleteModelCases( df=train_df, modelFormula=model_formula, - response=train_df$qcscore_train, + 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 <- .dropModelIntercept(model_matrix) model <- trainModel(model_matrix, train_df) + if(is.null(bestLambda)) { - bestLambda <- computeLambda(train_df, model_formula) + bestLambda <- computeLambda(model_matrix, train_df) } - coefs <- coef(model, s=bestLambda) + coefs <- coef(model, s=bestLambda)[coef(model, s=bestLambda)[,1]!= 0,,drop=FALSE] if (verbose) { message("Model coefficients for every term used in the formula:") @@ -518,20 +527,17 @@ computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALS full_matrix <- model.matrix(as.formula(model_formula), data=df[full_ok, , drop=FALSE]) - full_matrix <- .dropModelIntercept(full_matrix) - full_matrix <- full_matrix[, colnames(model_matrix), 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. - ## Predict only complete cases; incomplete cells keep QC_score = NA. - qc_score <- rep(NA_real_, nrow(df)) - qc_score[full_ok] <- as.vector( + ## Predict only complete cases; incomplete cells keep QScore = NA. + qscore <- rep(NA_real_, nrow(df)) + qscore[full_ok] <- as.vector( predict(model, s=bestLambda, newx=full_matrix, type="response") ) - spe$QC_score <- qc_score + spe$QScore <- qscore - metadata(spe)$QCScore_model <- list( + metadata(spe)$QScore_model <- list( model=model, bestLambda=bestLambda, model_formula=model_formula, @@ -551,7 +557,7 @@ computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALS #' @name trainModel #' @rdname trainModel #' @description -#' Fit a Ridge Logistic Regression Model +#' Fits a ridge logistic regression model. #' #' \code{trainModel} fits an L2-regularized (ridge) logistic regression #' using \pkg{glmnet}, given a design matrix and a training data frame. @@ -576,7 +582,7 @@ computeQCScore <- function(spe, bestLambda=NULL, modelFormula=NULL, verbose=FALS #' coef(fit, s = 0.01) trainModel <- function(modelMatrix, trainDF) { - model <- glmnet(x=modelMatrix, y=trainDF$qcscore_train, + model <- glmnet(x=modelMatrix, y=trainDF$QScore_train, family="binomial", lambda=NULL, alpha=0) return(model) } @@ -586,22 +592,22 @@ trainModel <- function(modelMatrix, trainDF) #' @rdname computeTrainDF #' #' @description -#' Build a Balanced Training Data Frame from a SpatialExperiment +#' Builds a balanced training data frame from a SpatialExperiment. #' #' \code{computeTrainDF} takes a \code{SpatialExperiment} object #' and assembles a balanced training set of “good” vs “bad” cells for #' subsequent model fitting. #' @param colData A per-cell metadata table. Typically #' `as.data.frame(colData(spe))`. Must include at least: -#' `cell_id`, raw metric columns named in `formulaVars` (e.g. +#' `cell_id`, raw metric columns named as in `formulaVars` (e.g. #' `log2SignalDensity`, `Area_um`, `log2Ctrl_total_ratio`, -#' optionally `log2AspectRatio`), and the corresponding outlier-label +#' optionally `log2AspectRatio` and `dist_border`), and the corresponding outlier-label #' columns referenced by `formulaVars`. #' @param formulaVars A named character vector mapping variable name to #' its outlier label column name, e.g. #' `c(log2SignalDensity="log2SignalDensity_outlier_train", ...)`. #' @param tech Character string with the acquisition technology. Used to -#' enable CosMx-specific handling for `log2AspectRatio`. Expected +#' enable CosMx-specific inclusion of border effect. Expected #' values include `"Nanostring_CosMx"` or `"Nanostring_CosMx_Protein"`. #' @param verbose Logical; print progress messages. #' @@ -610,7 +616,7 @@ trainModel <- function(modelMatrix, trainDF) #' #' @return #' A \code{data.frame} with one row per cell, including: -#' \code{qcscore_train} (0/1) indicating “bad” vs “good”, +#' \code{QScore_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. #' @@ -624,11 +630,11 @@ trainModel <- function(modelMatrix, trainDF) #' #' @examples #' example(spatialPerCellQC) -#' spe <- computeOutliersQCScore(spe) +#' spe <- computeOutliersQScore(spe) #' spe <- checkOutliers(spe) #' df_train <- computeTrainDF(colData(spe), metadata(spe)$formula_variables, #' metadata(spe)$technology) -#' table(df_train$qcscore_train) +#' table(df_train$QScore_train) #' #' @export computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) { @@ -637,7 +643,7 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) { train_bad_var <- character() train_good_var <- character() - stopifnot("log2SignalDensity is not included in the QC score formula.\nQC score cannot be computed"="log2SignalDensity" %in% names(out_var)) + stopifnot("log2SignalDensity is not included in the Quality Score formula.\nQuality Score cannot be computed"="log2SignalDensity" %in% names(out_var)) cfg <- list( log2SignalDensity=list(bad="LOW", good=c(0.90,0.99)), @@ -702,12 +708,12 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) { train_bad <- df |> dplyr::filter(cell_id %in% train_bad_var) |> - dplyr::mutate(qcscore_train=0) + dplyr::mutate(QScore_train=0) train_good <- df |> dplyr::filter(cell_id %in% train_good_var) |> dplyr::mutate( - qcscore_train=1, + QScore_train=1, is_a_bad_boy=cell_id %in% train_bad$cell_id ) @@ -742,8 +748,6 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) { return(train_df) } - - #' getModelFormula #' @name getModelFormula #' @rdname getModelFormula @@ -754,28 +758,22 @@ computeTrainDF <- function(colData, formulaVars, tech, verbose=FALSE) { #' (e.g. `"log2SignalDensity"`, `"Area_um"`, etc.) to their corresponding #' outlier label columns, typically from #' `metadata(spe)$formula_variables`. -#' @param verbose Logical. If `TRUE`, prints the final formula used for QC score #' @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) +getModelFormula <- function(metricList) { - out_var <- formulaVars - if ("log2AspectRatio" %in% names(out_var)) { - names(out_var)[grep("log2AspectRatio_outlier", out_var)] <- + out_var <- metricList + if ("log2AspectRatio" %in% out_var) { + out_var[grep("log2AspectRatio", out_var)] <- "I(abs(log2AspectRatio) * as.numeric(dist_border<50))" } - model_formula <- paste0("~(", paste(names(out_var), collapse = " + "), + model_formula <- paste0("~(", paste(out_var, collapse = " + "), ")^2", sep = "") - if (verbose) { - message("Final formula used for QC score computation:") - message(model_formula) - } - return(model_formula) } @@ -783,7 +781,7 @@ getModelFormula <- function(formulaVars, verbose=FALSE) #' @name dot-computeXenMerTrainSet #' @rdname dot-computeXenMerTrainSet #' @description -#' Internal: Build Training Set for Xenium & MERFISH +#' Internal: builds Training Set for Xenium & MERFISH. #' Splits a SpatialExperiment into “bad” vs “good” cells based on #' pre-computed outlier labels on log2SignalDensity. #' @param spe \code{SpatialExperiment} @@ -806,9 +804,9 @@ getModelFormula <- function(formulaVars, verbose=FALSE) #' @name dot-computeCosmxTrainSet #' @rdname dot-computeCosmxTrainSet #' @description -#' Internal: Build Training Set for CosMx +#' Internal: builds Training Set for CosMx. #' Splits a SpatialExperiment into “bad” vs “good” cells based on -#' outliers in aspect ratio near tissue border or low count area. +#' outliers in aspect ratio near FOV border or low signal density. #' @param spe \code{SpatialExperiment} #' @return #' A list with elements \code{bad} and \code{good}, each a data.frame @@ -837,9 +835,9 @@ getModelFormula <- function(formulaVars, verbose=FALSE) #' @name dot-computeCosmxProteinTrainSet #' @rdname dot-computeCosmxProteinTrainSet #' @description -#' Internal: Build Training Set for CosMx-Protein +#' Internal: builds Training Set for CosMx-Protein. #' Splits a SpatialExperiment into “bad” vs “good” cells based on -#' outliers in aspect ratio near tissue border or low count area. +#' outliers in aspect ratio near FOV border or low signal density. #' @param spe \code{SpatialExperiment} #' @return #' A list with elements \code{bad} and \code{good}, each a data.frame @@ -867,47 +865,46 @@ getModelFormula <- function(formulaVars, verbose=FALSE) } -#' computeQCScoreFlags -#' @name computeQCScoreFlags -#' @rdname computeQCScoreFlags +#' computeQScoreFlags +#' @name computeQScoreFlags +#' @rdname computeQScoreFlags #' @description -#' Compute flagged cells based on a manually chosen threshold on quality score +#' Computes flagged cells based on a manually chosen threshold on Quality Score (QS). #' -#' This function Compute flagged cells based on a manually chosen threshold on -#' quality score stored in `SpatialExperiment` object. +#' This function computes flagged cells based on a manually chosen threshold on +#' QS stored in `SpatialExperiment` object. #' -#' @param spe A `SpatialExperiment` object with spatial transcriptomics data. -#' @param qsThreshold Numeric threshold or quantile for quality score. Default -#' `0.5`. -#' @param useQSQuantiles Logical; if `TRUE`, treat `qsThreshold` as a -#' percentile. +#' @param spe A `SpatialExperiment` object with spatial omics data. +#' @param qsThreshold Numeric threshold or quantile for QS. Default is `0.5`. +#' @param useQSQuantiles Logical; if `TRUE`, treat `qsThreshold` as a percentile. #' -#' @return The `SpatialExperiment` object with added filter flags in `colData`. +#' @return The `SpatialExperiment` object with added low-quality flagged cells as +#' TRUE/FALSE in `colData`. #' #' @importFrom SummarizedExperiment colData #' @export #' @examples -#' example(computeQCScore) -#' spe <- computeQCScoreFlags(spe) -#' table(spe$low_qcscore) +#' example(computeQScore) +#' spe <- computeQScoreFlags(spe) +#' table(spe$low_QScore) #' # if fixed filters are defined we have an additional column #' spe <- computeThresholdFlags(spe) -#' spe <- computeQCScoreFlags(spe) -#' table(spe$low_threshold_qcscore) -computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { +#' spe <- computeQScoreFlags(spe) +#' table(spe$low_threshold_QScore) +computeQScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { stopifnot(is(spe, "SpatialExperiment")) - stopifnot("QC_score" %in% names(colData(spe))) + stopifnot("QScore" %in% names(colData(spe))) if(useQSQuantiles) { - spe$low_qcscore <- ifelse( - spe$QC_score < quantile(spe$QC_score, probs=qsThreshold), + spe$low_QScore <- ifelse( + spe$QScore < quantile(spe$QScore, probs=qsThreshold), TRUE, FALSE) } else { - spe$low_qcscore <- spe$QC_score < qsThreshold + spe$low_QScore <- spe$QScore < qsThreshold } if("threshold_flags" %in% names(colData(spe))) { - spe$low_threshold_qcscore <- (spe$low_qcscore & + spe$low_threshold_QScore <- (spe$low_QScore & spe$threshold_flags) } return(spe) @@ -917,11 +914,11 @@ computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { #' @name dot-checkSkw #' @rdname dot-checkSkw #' @description -#' Check skewness of metrics to choose outlier detection method. +#' Internal: checks skewness of metrics to choose outlier detection method. #' #' @param cd colData of `SpatialExperiment` object. #' @param metricList A character vector specifying the metrics to include in -#' the QC score formula. Defaults are "log2SignalDensity", "Area_um", +#' the Quality Score formula. Defaults are "log2SignalDensity", "Area_um", #' "log2AspectRatio", "log2Ctrl_total_ratio". #' #' @return @@ -954,11 +951,11 @@ computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { return(method) } -#' computeOutliersQCScore -#' @name computeOutliersQCScore -#' @rdname computeOutliersQCScore +#' computeOutliersQScore +#' @name computeOutliersQScore +#' @rdname computeOutliersQScore #' @description -#' Compute outlier cells for each metric that can be used in QC score formula +#' Computes outlier cells for each metric that can be used in Quality Score formula #' for SpatialExperiment. #' #' This function calculates outlier cells for each variable specified @@ -966,16 +963,18 @@ computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { #' the `colData` of the `SpatialExperiment` object as a minimum requirement. #' The user can choose which metrics to include among the following: Area_um, #' log2Ctrl_total_ratio, log2AspectRatio. For Xenium and Merfish datasets, -#' log2AspectRatio is automatically removed from the formula. -#' +#' even if outliers for log2AspectRatio are detected, neither the metric will be +#' included in the Quality Score formula nor nor the outliers will be used for +#' training. +#' #' @param spe A `SpatialExperiment` object with spatial omics data. #' @param metricList A character vector specifying the metrics to include in -#' the QC score formula. Default is `c("log2SignalDensity", "Area_um", +#' the Quality Score formula. Default is `c("log2SignalDensity", "Area_um", #' "log2AspectRatio", "log2Ctrl_total_ratio")`. #' #' @return The `SpatialExperiment` object with added outlier variables in -#' `colData` and the temporary QCScore metric variables that in the -#' `metadata`. +#' `colData` and the temporary metrics considered for Quality Score computation +#' stored in the object `metadata`. #' #' @details The function computes outliers for each specified metric after #' automatically choosing the appropriate method according to the skewness of @@ -983,10 +982,10 @@ computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { #' Internally the function: #' \enumerate{ #' \item Calls \code{.checkSkw()} to choose the proper outlier detection -#' method according to the variable skewness, +#' method according to the variable skewness; #' \item Calls \code{computeSpatialOutlier()} on each included metric to get -#' fences, -#' \item Labels cells as “LOW”/“HIGH” outliers or “NO” +#' fences; +#' \item Labels cells as “LOW”/“HIGH” outliers or “NO”. #' } #' #' @importFrom SummarizedExperiment colData @@ -997,9 +996,9 @@ computeQCScoreFlags <- function(spe, qsThreshold=0.5, useQSQuantiles=FALSE) { #' @examples #' example(readCosmxSPE) #' spe <- spatialPerCellQC(spe) -#' spe <- computeOutliersQCScore(spe) +#' spe <- computeOutliersQScore(spe) #' table(spe$log2SignalDensity_outlier_train) -computeOutliersQCScore <- function(spe, metricList=c("log2SignalDensity","Area_um", +computeOutliersQScore <- function(spe, metricList=c("log2SignalDensity","Area_um", "log2AspectRatio", "log2Ctrl_total_ratio")) { stopifnot(is(spe, "SpatialExperiment")) @@ -1024,7 +1023,7 @@ computeOutliersQCScore <- function(spe, metricList=c("log2SignalDensity","Area_u idx <- which(is.na(method)) for (i in idx) { message("Metric ", names(method)[i], " produced NA skewness values.", - "It will be removed from the QC score formula.") + "It will be removed from the Quality Score formula.") method <- method[-i] } } @@ -1083,25 +1082,25 @@ computeOutliersQCScore <- function(spe, metricList=c("log2SignalDensity","Area_u #' @rdname checkOutliers #' @description #' Checks if computed outliers meet the minimum numerical requirement, being -#' at least 0.1% of total cells for each metric to be used in QC score formula. +#' at least 0.1\% of total cells for each metric to be used in Quality Score formula. #' If the requirement is not met, the variable is removed from the formula. #' #' @param spe A `SpatialExperiment` object with spatial omics data. #' @param verbose Logical. If `TRUE`, prints how many outliers were found for #' each metric. #' -#' @return The `SpatialExperiment` object with added QCScore metric variables -#' in the `metadata`. +#' @return The `SpatialExperiment` object with added metric variables that +#' will be used for Quality Score computation stored in the object` metadata`. #' #' @details The function checks if computed outliers for each metric meet -#' the minimum number to get the metric included in the QC score formula. +#' the minimum number to get the metric included in the Quality Score formula. #' If verbose is TRUE, it also prints how many outliers were found for each #' metric. #' #' @importFrom SummarizedExperiment colData #' @export #' @examples -#' example(computeOutliersQCScore) +#' example(computeOutliersQScore) #' spe <- checkOutliers(spe, verbose = TRUE) #' metadata(spe)$formula_variables checkOutliers <- function(spe, verbose=FALSE) { @@ -1110,7 +1109,7 @@ checkOutliers <- function(spe, verbose=FALSE) { m1 <- paste0("Not enough outlier cells for ", var, ".\n") m2 <- switch(warnstop, s="In this version of SpaceTrooper, presence of outliers - for log2SignalDensity is required to compute QC score.", + for log2SignalDensity is required to compute Quality Score.", w="This variable will not be used in the final formula") return(paste0(m1, m2)) } @@ -1119,12 +1118,14 @@ checkOutliers <- function(spe, verbose=FALSE) { if (verbose) { for (i in names(out_var)) { message("Outliers found for ", i, ":") - message(paste(table(cd[[out_var[i]]]), collapse = " ")) + 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])) + } } } stopifnot( - "log2SignalDensity is not included in the QC score formula.\n - QC score cannot be computed"= + "log2SignalDensity is not included in the Quality Score formula.\n + Quality Score cannot be computed"= "log2SignalDensity" %in% names(out_var) ) cfg <- list( @@ -1155,8 +1156,8 @@ checkOutliers <- function(spe, verbose=FALSE) { pat <- paste0(var, "_outlier") is_cosmx <- metadata(spe)$technology %in% c("Nanostring_CosMx", "Nanostring_CosMx_Protein") + idx <- grep(pat, out_var) if (is_cosmx && (var %in% names(out_var))) { - idx <- grep(pat, out_var) if (length(idx)) { col <- out_var[idx] labs <- factor(cd[[col]], levels=c("LOW","HIGH","NO")) @@ -1168,11 +1169,13 @@ checkOutliers <- function(spe, verbose=FALSE) { if (low < nmin && high < nmin) { warning(warnstopmsg(var, "w")) - out_var <- out_var[-grep(pat, out_var)] + out_var <- out_var[-idx] } } } else { - out_var <- out_var[-grep(pat, out_var)] + if (length(idx)) { + out_var <- out_var[-idx] + } } metadata(spe)$formula_variables <- out_var return(spe) @@ -1196,7 +1199,7 @@ checkOutliers <- function(spe, verbose=FALSE) { stop("log2CountArea found in colData.\n", "Please updated the object with the latest version of spatialPerCellQC.") } - spe1 <- computeOutliersQCScore(spe, metricList) + spe1 <- computeOutliersQScore(spe, metricList) spe1 <- checkOutliers(spe1, verbose) out_var <- metadata(spe1)$formula_variables @@ -1206,29 +1209,35 @@ checkOutliers <- function(spe, verbose=FALSE) { return(list(df=df, out_var=out_var, tech=tech)) } -#' applyQCScoreModel -#' +#' .applyQScoreModel +#' @name .applyQScoreModel +#' @rdname dot-applyQScoreModel +#' #' @description -#' Apply a previously trained QC score model to a new SpatialExperiment object. +#' Internal: 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 -#' The authors do not reccommend applying a QC score model trained on one +#' The authors do not recommend applying a Quality Score (QS) model trained on one #' dataset to a different dataset, but this could be useful in some cases, -#' for example to transfer a model trained on a dataset core and the applied to -#' other cores from the same experiment. -#' The QC score model should have been trained on a dataset with the same -#' set of QC metrics and similar data distribution to have better -#' predictions. The function will check for the presence of required model -#' variables in the new dataset and will handle missing or extra variables -#' accordingly. +#' for example to transfer a model trained on a dataset sourced from a tissue +#' slide sample and then applied to other slides from the same experiment. +#' When transferring the model from a reference dataset of the same technology +#' of the query dataset, there is no need to set a common formula. However, +#' when transferring from a different technology, reference QS computation +#' should be performed providing a formula that excludes untransferable or unwanted +#' variables (e.g. computing QS on a reference CosMx dataset omitting the border +#' effect term to transfer the QS model to a Xenium query dataset. This is done by +#' providing a custom formula in computeQScore, see link computeQScore). The function will check +#' for the presence of required model variables in the new dataset and will handle +#' missing or extra variables accordingly. #' #' @param spe A `SpatialExperiment` object with QC metrics already computed. -#' @param qcModel A QC score model object, usually stored in -#' `metadata(spe)$QCScore_model`. +#' @param qsModel A Quality Score model object, usually stored in +#' `metadata(spe)$QScore_model`. #' @param scoreName Name of the output column in `colData`. #' -#' @return A `SpatialExperiment` object with added QC score in `colData`. +#' @return A `SpatialExperiment` object with added Quality Score in `colData`. #' #' @export #' @examples @@ -1245,18 +1254,18 @@ checkOutliers <- function(spe, verbose=FALSE) { #' spe_test <- spe[, -idx] #' #' ## Train the Quality Control (QC) score model on one dataset -#' spe_train <- computeQCScore(spe_train) -#' qc_model <- metadata(spe_train)$QCScore_model +#' spe_train <- computeQScore(spe_train) +#' qc_model <- metadata(spe_train)$QScore_model #' #' ## Apply the trained model to another dataset -#' spe_test <- applyQCScoreModel( +#' spe_test <- applyQScoreModel( #' spe=spe_test, #' qcModel=qc_model, -#' scoreName="QC_score_transferred" +#' scoreName="QScore_transferred" #' ) #' -#' summary(spe_test$QC_score_transferred) -applyQCScoreModel <- function(spe, qcModel, scoreName="QC_score") { +#' summary(spe_test$QScore_transferred) +.applyQScoreModel <- function(spe, qsModel, scoreName="QScore") { stopifnot(is(spe, "SpatialExperiment")) required <- c( @@ -1267,27 +1276,25 @@ applyQCScoreModel <- function(spe, qcModel, scoreName="QC_score") { ) stopifnot( - "qcModel is missing required fields" = - all(required %in% names(qcModel)) + "qsModel is missing required fields" = + all(required %in% names(qsModel)) ) df <- as.data.frame(colData(spe)) ok <- .filterCompleteModelCases( df=df, - modelFormula=qcModel$model_formula, + modelFormula=qsModel$model_formula, context="cells" ) new_matrix <- model.matrix( - as.formula(qcModel$model_formula), + as.formula(qsModel$model_formula), data=df[ok, , drop=FALSE] ) - new_matrix <- .dropModelIntercept(new_matrix) - missing_cols <- setdiff( - qcModel$model_matrix_colnames, + qsModel$model_matrix_colnames, colnames(new_matrix) ) @@ -1304,7 +1311,7 @@ applyQCScoreModel <- function(spe, qcModel, scoreName="QC_score") { extra_cols <- setdiff( colnames(new_matrix), - qcModel$model_matrix_colnames + qsModel$model_matrix_colnames ) if (length(extra_cols) > 0L) { @@ -1315,14 +1322,14 @@ applyQCScoreModel <- function(spe, qcModel, scoreName="QC_score") { ) } - new_matrix <- new_matrix[, qcModel$model_matrix_colnames, drop=FALSE] + new_matrix <- new_matrix[, qsModel$model_matrix_colnames, drop=FALSE] score <- rep(NA_real_, nrow(df)) score[ok] <- as.vector( predict( - qcModel$model, - s=qcModel$bestLambda, + qsModel$model, + s=qsModel$bestLambda, newx=new_matrix, type="response" ) @@ -1332,7 +1339,7 @@ applyQCScoreModel <- function(spe, qcModel, scoreName="QC_score") { cd[[scoreName]] <- score colData(spe) <- cd - metadata(spe)$QCScore_model_applied <- qcModel + metadata(spe)$QScore_model_applied <- qsModel return(spe) } diff --git a/R/readMerfish.R b/R/readMerfish.R index d7d0dd8..b629640 100644 --- a/R/readMerfish.R +++ b/R/readMerfish.R @@ -191,7 +191,7 @@ computeMissingMetricsMerfish <- function(polFile, coldata, warning("Volume column not found in colData.\nComputing area from polygons instead.") area <- computeAreaFromPolygons(polygons) } else { - warning("Volume is used to compute QC score for MERFISH technology. + warning("Volume is used to compute Quality Score for MERFISH technology. For simplicity, it is renamed as Area_um.") area <- cd$volume } diff --git a/R/spatialQCPlots.R b/R/spatialQCPlots.R index ad6a231..c730dac 100644 --- a/R/spatialQCPlots.R +++ b/R/spatialQCPlots.R @@ -10,12 +10,13 @@ #' @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 alphaNumbers Numeric transparency for FoV labels. Default: `0.8`. +#' @param numberSize Numeric size for the FoV labels. Default: `1`. +#' @param numbersAlpha Numeric transparency for FoV labels. Default: `0.8`. #' @param fovDim numeric with two named dimensions xdim, ydim. (Default is #' metadata(spe)$fov_dim) -#' @param size Numeric point size for the cell centroids. Default: `0.05`. -#' @param alpha Numeric transparency for the cell centroids. Default: `0.8`. #' @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. @@ -29,7 +30,7 @@ #' - `metadata(spe)$fov_positions`: a matrix or data.frame #' (or list with named elements) containing at minimum `x_global_px`, `y_global_px`, #' and `fov`. Values `x_global_px`/`y_global_px` are in pixels and represent -#' the origin (top-left) of each FoV. +#' the origin (bottom-left) of each FoV. #' - `metadata(spe)$fov_dim` (or the `fovDim` argument): a named numeric with #' `xdim` and `ydim` giving FoV width/height in pixels. #' @@ -61,9 +62,10 @@ #' g <- plotCellsFovs(spe) #' print(g) plotCellsFovs <- function(spe, sampleId=unique(spe$sample_id), - pointCol="firebrick", numbersCol="black", - alphaNumbers=0.8, fovDim=metadata(spe)$fov_dim, - size=0.05, alpha=0.8, + 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) { stopifnot(is(spe, "SpatialExperiment")) @@ -78,7 +80,7 @@ plotCellsFovs <- function(spe, sampleId=unique(spe$sample_id), y=.data[[y_coord]]), colour=pointCol, fill=pointCol, - size=size, alpha=alpha) + + size=pointSize, alpha=pointAlpha) + annotate("rect", xmin=metadata(spe)$fov_positions["x_global_px"][ , , drop=TRUE], xmax=metadata(spe)$fov_positions["x_global_px"][ , , drop=TRUE] + @@ -92,7 +94,8 @@ plotCellsFovs <- function(spe, sampleId=unique(spe$sample_id), y=metadata(spe)$fov_positions["y_global_px"][,,drop=TRUE]+ fovDim[["ydim"]]/2, label=metadata(spe)$fov_positions["fov"][,,drop=TRUE]), - color=numbersCol, fontface="bold", alpha=alphaNumbers) + + color=numbersCol, size = numberSize, alpha=numbersAlpha, + fontface="bold") + ggtitle(sampleId) + .fov_image_theme(backColor="white", backBorder="white", titleCol="black") + ggplot2::coord_fixed() @@ -436,14 +439,12 @@ 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 mapNumbersCol A character string specifying the color of the -#' numbers on the map. Default is `"black"`. -#' @param mapAlphaNumbers A numeric value specifying the transparency of the -#' numbers on the map. Default is `0.8`. -#' @param csize A numeric value specifying the size of the points in the map. -#' Default is `0.05`. -#' @param calpha A numeric value specifying the transparency of the points in -#' the map. Default is `0.8`. +#' @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 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. @@ -474,8 +475,9 @@ plotPolygons <- function(spe, colourBy="darkgrey", colourLog=FALSE, #' plotZoomFovsMap(spe, fovs=16, title="FOV 16") plotZoomFovsMap <- function(spe, fovs=NULL, title=NULL, mapPointCol="darkmagenta", - mapNumbersCol="black", - mapAlphaNumbers=0.8, + mapPointSize=0.5, mapPointAlpha=0.8, + fovNumbersCol="black", fovNumberSize=1, + fovNumbersAlpha=0.8, csize=0.05, calpha=0.8, scaleBars=NULL, scaleBarMap=TRUE, @@ -491,8 +493,9 @@ plotZoomFovsMap <- function(spe, fovs=NULL, title=NULL, scaleBarPol <- scaleBars } map <- plotCellsFovs(spefovs, pointCol=mapPointCol, - numbersCol=mapNumbersCol, alphaNumbers=mapAlphaNumbers, - sampleId=NULL, size=csize, alpha=calpha, scaleBar=scaleBarMap) + pointSize=mapPointSize, pointAlpha=mapPointAlpha, + numbersCol=fovNumbersCol, numberSize=fovNumberSize, + numbersAlpha=fovNumbersAlpha, sampleId=NULL, scaleBar=scaleBarMap) g2 <- plotPolygons(spefovs, sampleId=NULL, scaleBar=scaleBarPol, ...) final_plot <- ggpubr::ggarrange(map, g2, ncol=2) if (!is.null(title)) { diff --git a/README.md b/README.md index 1627303..fd21d22 100644 --- a/README.md +++ b/README.md @@ -117,8 +117,8 @@ spe <- readAndAddPolygonsToSPE(spe, boundariesType="csv") spe <- spatialPerCellQC(spe) # 4. Compute QS and, optionally, flag cells with a score higher than 'qsThreshold'. -spe <- computeQCScore(spe) -spe <- computeQCScoreFlags(spe, qsThreshold=0.5) +spe <- computeQScore(spe) +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). diff --git a/man/applyQCScoreModel.Rd b/man/applyQCScoreModel.Rd deleted file mode 100644 index 3e67195..0000000 --- a/man/applyQCScoreModel.Rd +++ /dev/null @@ -1,61 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/QC.R -\name{applyQCScoreModel} -\alias{applyQCScoreModel} -\title{applyQCScoreModel} -\usage{ -applyQCScoreModel(spe, qcModel, scoreName = "QC_score") -} -\arguments{ -\item{spe}{A `SpatialExperiment` object with QC metrics already computed.} - -\item{qcModel}{A QC score model object, usually stored in -`metadata(spe)$QCScore_model`.} - -\item{scoreName}{Name of the output column in `colData`.} -} -\value{ -A `SpatialExperiment` object with added QC score in `colData`. -} -\description{ -Apply a previously trained QC score model to a new SpatialExperiment object. -See details for important considerations when applying a model to a -different dataset. -} -\details{ -The authors do not reccommend applying a QC score model trained on one -dataset to a different dataset, but this could be useful in some cases, -for example to transfer a model trained on a dataset core and the applied to -other cores from the same experiment. -The QC score model should have been trained on a dataset with the same -set of QC metrics and similar data distribution to have better -predictions. The function will check for the presence of required model -variables in the new dataset and will handle missing or extra variables -accordingly. -} -\examples{ - -example(readCosmxSPE) - -## Compute per-cell quality control metrics -spe <- spatialPerCellQC(spe) - -## Split the object only for example purposes -set.seed(1998) -idx <- sample(seq_len(ncol(spe)), floor(ncol(spe) / 2)) -spe_train <- spe[, idx] -spe_test <- spe[, -idx] - -## Train the Quality Control (QC) score model on one dataset -spe_train <- computeQCScore(spe_train) -qc_model <- metadata(spe_train)$QCScore_model - -## Apply the trained model to another dataset -spe_test <- applyQCScoreModel( - spe=spe_test, - qcModel=qc_model, - scoreName="QC_score_transferred" -) - -summary(spe_test$QC_score_transferred) -} diff --git a/man/checkOutliers.Rd b/man/checkOutliers.Rd index d29bffa..0cb22d7 100644 --- a/man/checkOutliers.Rd +++ b/man/checkOutliers.Rd @@ -13,22 +13,22 @@ checkOutliers(spe, verbose = FALSE) each metric.} } \value{ -The `SpatialExperiment` object with added QCScore metric variables - in the `metadata`. +The `SpatialExperiment` object with added metric variables that +will be used for Quality Score computation stored in the object` metadata`. } \description{ Checks if computed outliers meet the minimum numerical requirement, being -at least 0.1% of total cells for each metric to be used in QC score formula. +at least 0.1\% of total cells for each metric to be used in Quality Score formula. If the requirement is not met, the variable is removed from the formula. } \details{ The function checks if computed outliers for each metric meet -the minimum number to get the metric included in the QC score formula. +the minimum number to get the metric included in the Quality Score formula. If verbose is TRUE, it also prints how many outliers were found for each metric. } \examples{ -example(computeOutliersQCScore) +example(computeOutliersQScore) spe <- checkOutliers(spe, verbose = TRUE) metadata(spe)$formula_variables } diff --git a/man/computeLambda.Rd b/man/computeLambda.Rd index 6b4060e..e7becff 100644 --- a/man/computeLambda.Rd +++ b/man/computeLambda.Rd @@ -4,11 +4,11 @@ \alias{computeLambda} \title{computeLambda} \usage{ -computeLambda(trainDF, modelFormula) +computeLambda(modelMatrix, trainDF) } \arguments{ \item{trainDF}{`data.frame` -A data frame for training that must include: +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.} @@ -22,18 +22,18 @@ A character string representing the model formula `cv.glmnet` that minimizes the cross-validation error. } \description{ -Compute Optimal Ridge Regularization Parameter \eqn{\lambda} via -Cross-Validation +Computes optimal ridge regularization parameter \eqn{\lambda} via +k-fold cross-validation. \code{computeLambda} performs ridge (L2) logistic regression with -cross-validation to identify the optimal regularization parameter -\eqn{\lambda} for a binary response. +k-fold cross-validation to identify the optimal regularization parameter +\eqn{\lambda} for Quality Score (QS) model training. } \details{ Internally, the function: - Constructs the design matrix via \code{model.matrix()}, - Runs ridge logistic regression cross-validation using `cv.glmnet` with `alpha = 0`, - Extracts and returns `ridge_cv$lambda.min`. + 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) diff --git a/man/computeOutliersQCScore.Rd b/man/computeOutliersQScore.Rd similarity index 68% rename from man/computeOutliersQCScore.Rd rename to man/computeOutliersQScore.Rd index 21cfd69..844bf39 100644 --- a/man/computeOutliersQCScore.Rd +++ b/man/computeOutliersQScore.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/QC.R -\name{computeOutliersQCScore} -\alias{computeOutliersQCScore} -\title{computeOutliersQCScore} +\name{computeOutliersQScore} +\alias{computeOutliersQScore} +\title{computeOutliersQScore} \usage{ -computeOutliersQCScore( +computeOutliersQScore( spe, metricList = c("log2SignalDensity", "Area_um", "log2AspectRatio", "log2Ctrl_total_ratio") @@ -14,16 +14,16 @@ computeOutliersQCScore( \item{spe}{A `SpatialExperiment` object with spatial omics data.} \item{metricList}{A character vector specifying the metrics to include in -the QC score formula. Default is `c("log2SignalDensity", "Area_um", +the Quality Score formula. Default is `c("log2SignalDensity", "Area_um", "log2AspectRatio", "log2Ctrl_total_ratio")`.} } \value{ The `SpatialExperiment` object with added outlier variables in -`colData` and the temporary QCScore metric variables that in the -`metadata`. +`colData` and the temporary metrics considered for Quality Score computation +stored in the object `metadata`. } \description{ -Compute outlier cells for each metric that can be used in QC score formula +Computes outlier cells for each metric that can be used in Quality Score formula for SpatialExperiment. This function calculates outlier cells for each variable specified @@ -31,7 +31,9 @@ in `metricList` for a `SpatialExperiment`. log2SignalDensity must be present in the `colData` of the `SpatialExperiment` object as a minimum requirement. The user can choose which metrics to include among the following: Area_um, log2Ctrl_total_ratio, log2AspectRatio. For Xenium and Merfish datasets, -log2AspectRatio is automatically removed from the formula. +even if outliers for log2AspectRatio are detected, neither the metric will be +included in the Quality Score formula nor nor the outliers will be used for +training. } \details{ The function computes outliers for each specified metric after @@ -40,15 +42,15 @@ the distribution. Internally the function: \enumerate{ \item Calls \code{.checkSkw()} to choose the proper outlier detection - method according to the variable skewness, + method according to the variable skewness; \item Calls \code{computeSpatialOutlier()} on each included metric to get - fences, - \item Labels cells as “LOW”/“HIGH” outliers or “NO” + fences; + \item Labels cells as “LOW”/“HIGH” outliers or “NO”. } } \examples{ example(readCosmxSPE) spe <- spatialPerCellQC(spe) -spe <- computeOutliersQCScore(spe) +spe <- computeOutliersQScore(spe) table(spe$log2SignalDensity_outlier_train) } diff --git a/man/computeQCScore.Rd b/man/computeQCScore.Rd deleted file mode 100644 index edc0af6..0000000 --- a/man/computeQCScore.Rd +++ /dev/null @@ -1,85 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/QC.R -\name{computeQCScore} -\alias{computeQCScore} -\title{computeQCScore} -\usage{ -computeQCScore(spe, bestLambda = NULL, modelFormula = NULL, verbose = FALSE) -} -\arguments{ -\item{spe}{A `SpatialExperiment` object with spatial transcriptomics data.} - -\item{bestLambda}{the best lambda typically computed using `computeLambda`.} - -\item{modelFormula}{a character string representing the model formula to be -used for training the model. If NULL, the formula is automatically generated -based on the available metrics and outliers in the dataset. -See Details for more information. -Note that the automatically generated formula will include interaction -terms between the metrics, and will exclude metrics with insufficient -outliers (< 0.1% of the dataset). If a custom `modelFormula` is provided, -it will be used as is without modification or checks for outlier counts.} - -\item{verbose}{logical for having a verbose output. Default is FALSE.} -} -\value{ -The `SpatialExperiment` object with added QC score in `colData`. -} -\description{ -Compute QC score and automatically define weights for QC score -through glm training. This function computes QC score with a formula -defined on the metrics as "log2SignalDensity", "Area_um", -"log2AspectRatio", "log2Ctrl_total_ratio", as computed from -the `spatialPerCellQC` function. -It automatically computes the number of available outliers for each available -metric, as they are needed for the glm training. -See Details for further information. -} -\details{ -For all the techologies, the QC Score formula depends on the follow -metrics: - -QC score ~ count density - aspect ratio - control-total ratio - size - -Where count density is the total counts-to-size ratio, aspect ratio -represents the ratio between the width and the height of the cell -(computed from the provided polygons if not already present in the experiment -metadata) and control-total ratio is the aspecific signal; -size is the area for CosMx and Xenium, while it is the volume for Merfish. -For each couple of variables interaction terms are computed. - -Additionally, for CosMx datasets, the distance from the border of the FOV is -also included in the formula as a metric to take into account. -For Xenium and Merscope datasets, QC score cannot depend on FoV border effect, -as no FOV border effect was captured through this metric. - -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 < 0.1% out of the -entire dataset, the metric will be excluded from the QC score formula. - -- Model fitting: ridge (L2) logistic regression is fitted (via `glmnet`) on -the balanced training set. The function uses `trainModel()` for fitting -and `computeLambda()` (cross‑validation) to select lambda unless -`bestLambda` is supplied. - -- Lambda details: because of the randomness in the training set selection, -results may vary so that it is possible to set a fixed lambda value -previously computed with `computeLambda` preceeded by `computeTrainDF` and -`getModelFormula`. -This is useful for reproducibility across different runs. -Otherwise, an easier way is to let be lambda computed internally, just set -a seed with `set.seed()` before running `computeQCScore`. - -The computed model output is stored in `metadata(spe)$QCScore_model` -for inspection and reuse (see also \code{\link{applyQCScoreModel}}). -} -\examples{ -example(spatialPerCellQC) -set.seed(1998) -spe <- computeQCScore(spe) -summary(spe$QC_score) -} diff --git a/man/computeQCScoreFlags.Rd b/man/computeQCScoreFlags.Rd deleted file mode 100644 index 41100a4..0000000 --- a/man/computeQCScoreFlags.Rd +++ /dev/null @@ -1,35 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/QC.R -\name{computeQCScoreFlags} -\alias{computeQCScoreFlags} -\title{computeQCScoreFlags} -\usage{ -computeQCScoreFlags(spe, qsThreshold = 0.5, useQSQuantiles = FALSE) -} -\arguments{ -\item{spe}{A `SpatialExperiment` object with spatial transcriptomics data.} - -\item{qsThreshold}{Numeric threshold or quantile for quality score. Default -`0.5`.} - -\item{useQSQuantiles}{Logical; if `TRUE`, treat `qsThreshold` as a -percentile.} -} -\value{ -The `SpatialExperiment` object with added filter flags in `colData`. -} -\description{ -Compute flagged cells based on a manually chosen threshold on quality score - -This function Compute flagged cells based on a manually chosen threshold on -quality score stored in `SpatialExperiment` object. -} -\examples{ -example(computeQCScore) -spe <- computeQCScoreFlags(spe) -table(spe$low_qcscore) -# if fixed filters are defined we have an additional column -spe <- computeThresholdFlags(spe) -spe <- computeQCScoreFlags(spe) -table(spe$low_threshold_qcscore) -} diff --git a/man/computeQScore.Rd b/man/computeQScore.Rd new file mode 100644 index 0000000..730e87b --- /dev/null +++ b/man/computeQScore.Rd @@ -0,0 +1,89 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/QC.R +\name{computeQScore} +\alias{computeQScore} +\title{computeQScore} +\usage{ +computeQScore(spe, bestLambda = NULL, modelFormula = NULL, verbose = FALSE) +} +\arguments{ +\item{spe}{A `SpatialExperiment` object with spatial omics data.} + +\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{verbose}{logical for having a verbose output. Default is FALSE.} +} +\value{ +The `SpatialExperiment` object with added Quality Score in `colData`. +} +\description{ +Computes Quality Score (QS) and automatically define weights of the QS model +terms through glm training. This function computes QS with a formula +defined on the metrics "log2SignalDensity", "Area_um", +"log2AspectRatio", "log2Ctrl_total_ratio", as computed from +the `spatialPerCellQC` function. +It automatically computes the number of outliers for each available +metric, which are needed for the glm training. +See Details for further information. +} +\details{ +For all the technologies, the QS formula can contain +only the following terms at most: + +Quality Score ~ signal density - size - border effect - background signal + +Where signal density is the total counts/protein intensities-to-size ratio, +border effect is represented by the aspect ratio (i.e. width/height of the cell) +and is considered only for cells with a distance from the FOV border less than +50 px and only for CosMx datasets. Background signal is represented by the ratio +of control probe counts to total counts (protein intensities) and cell size is +the area for CosMx and Xenium, while it is the volume for Merfish. +For each couple of variables interaction terms are computed. + +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. + +- Model fitting: ridge (L2) logistic regression is fitted (via `glmnet`) on +the balanced training set. The function uses `trainModel()` for fitting +and `computeLambda()` to select lambda unless `bestLambda` is supplied. + +- Lambda details: because of the randomness in the training set selection, +results may vary. The easiest way to ensure reproducibility is to let lambda +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: +`~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2`. + When user-provided, the formula must follow the same default syntax and +removed (or added) terms should be written exactly as in the default formula, +e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))` must have spaces +around the `*` and `<` operators. +In any case, metrics with insufficient outliers (less than 0.1\% of the dataset) +will be excluded from the QS formula. + +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 +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 +details on scalability and coefficient transfer. +} +\examples{ +example(spatialPerCellQC) +set.seed(1998) +spe <- computeQScore(spe) +summary(spe$QScore) +} diff --git a/man/computeQScoreFlags.Rd b/man/computeQScoreFlags.Rd new file mode 100644 index 0000000..78054a9 --- /dev/null +++ b/man/computeQScoreFlags.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/QC.R +\name{computeQScoreFlags} +\alias{computeQScoreFlags} +\title{computeQScoreFlags} +\usage{ +computeQScoreFlags(spe, qsThreshold = 0.5, useQSQuantiles = FALSE) +} +\arguments{ +\item{spe}{A `SpatialExperiment` object with spatial omics data.} + +\item{qsThreshold}{Numeric threshold or quantile for QS. Default is `0.5`.} + +\item{useQSQuantiles}{Logical; if `TRUE`, treat `qsThreshold` as a percentile.} +} +\value{ +The `SpatialExperiment` object with added low-quality flagged cells as +TRUE/FALSE in `colData`. +} +\description{ +Computes flagged cells based on a manually chosen threshold on Quality Score (QS). + +This function computes flagged cells based on a manually chosen threshold on +QS stored in `SpatialExperiment` object. +} +\examples{ +example(computeQScore) +spe <- computeQScoreFlags(spe) +table(spe$low_QScore) +# if fixed filters are defined we have an additional column +spe <- computeThresholdFlags(spe) +spe <- computeQScoreFlags(spe) +table(spe$low_threshold_QScore) +} diff --git a/man/computeSpatialOutlier.Rd b/man/computeSpatialOutlier.Rd index 929208d..f4ae52e 100644 --- a/man/computeSpatialOutlier.Rd +++ b/man/computeSpatialOutlier.Rd @@ -13,43 +13,39 @@ computeSpatialOutlier( ) } \arguments{ -\item{spe}{a SpatialExperiment object with target_counts, area in micron -and log2 of the aspect ratio in the `colData`.} +\item{spe}{a SpatialExperiment object with quantitative cell metrics in the `colData`.} -\item{computeBy}{character indicating a `colData` column name on which -compute the outlier.} +\item{computeBy}{character indicating a `colData` column name to compute the outliers on.} \item{method}{one of `mc`, `scuttle`, `both`. -Use `mc` for medcouple, `scuttle` for median absolute deviations as computed -in `scuttle`, `both` for computing both of them.} +Use `mc` for medcouple, `scuttle` for MAD as computed in `scuttle`, `both` +for computing both of them.} -\item{mcDoScale}{logical indicating if the values to compute the medcouple -for the outlier detection should be scaled (default is FALSE, as suggested -by the original Medcouple authors.). See \link[robustbase]{mc} for further -readings.} +\item{mcDoScale}{logical indicating if the metric values should be scaled +to compute medcouple (default is FALSE, as suggested by the original medcouple +authors). See \link[robustbase]{mc} for further readings.} \item{scuttleType}{One of `"both"`, `"lower"`, `"higher"` for scuttle method.} } \value{ a SpatialExperiment object with additional column(s) (named as -the column name indicated in `column_by` followed by the outlier_sc/mc -nomenclature) with the outlier detection as `outlier.filter` logical class -object. This allows to store the thresholds as attributes of the column. -use attr(,"thresholds") to retrieve them. +the column name indicated in `computeBy` followed by the outlier_sc/mc). +Resulting columns are `outlier.filter` logical class objects, allowing to store +the thresholds as attributes of the columns. +Use attr(,"thresholds") to retrieve them. } \description{ -Computes outliers based on the Area (in micron) of the experiment. -It gives the possibility to choose between the medcouple (mc method argument) -and the MADs (scuttle method argument). +Computes outliers for quantitative cell metadata based on a statistical test: +medcouple (`mc` method argument) or the Median Absolute Deviation (MAD, `scuttle` +method argument), when the metric distribution is respectively asymmetric or symmetric. } \details{ The medcouple method is a measure for the skeweness of univariate distribution as described in Hubert M. et al. (2008). In particular, the computed medcouple value must be in a range between -0.6 -and 0.6 to computed adjusted boxplots and perform the outlier detection. -For median absolute deviations (MADs) method we just wrap the isOutlier -function in the scuttle package. Please see McCarthy DJ et al (2017) -for further details. +and 0.6 to compute adjusted boxplots and perform the outlier detection. +For MAD method we just wrap the isOutlier function in the scuttle package. +Please see McCarthy DJ et al (2017) for further details. } \examples{ example(spatialPerCellQC) diff --git a/man/computeThresholdFlags.Rd b/man/computeThresholdFlags.Rd index fc27cdd..150e512 100644 --- a/man/computeThresholdFlags.Rd +++ b/man/computeThresholdFlags.Rd @@ -7,7 +7,7 @@ computeThresholdFlags(spe, totalThreshold = 0, ctrlTotRatioThreshold = 0.1) } \arguments{ -\item{spe}{A `SpatialExperiment` object with spatial transcriptomics data.} +\item{spe}{A `SpatialExperiment` object with spatial omics data.} \item{totalThreshold}{A numeric value for the threshold of total counts to identify cells with low counts. Default is `0`.} @@ -17,19 +17,20 @@ control-to-total ratio to flag cells over a certain threshold. Default is `0.1`.} } \value{ -The `SpatialExperiment` object with added filter flags in `colData`. +The `SpatialExperiment` object with added flagged cells as TRUE/FALSE + in `colData`. } \description{ -Compute Flagged cells using fixed thresholds for SpatialExperiment. +Computes flagged cells using fixed thresholds for SpatialExperiment. -This function calculates flagged cells only for total counts and control on -total probe counts ratio using fixed thresholds for a `SpatialExperiment` +This function calculates flagged cells only for total counts (protein intensities) +and control-to-total probe counts ratio using fixed thresholds for a `SpatialExperiment` object. } \details{ The function flags cells basing on zero counts and control-to-total ratio to identify junk cells. -It also combines these flags into a single filter flag. +It also combines these two flags into a single flag metric named `threshold_flags`. } \examples{ example(readCosmxSPE) diff --git a/man/computeTrainDF.Rd b/man/computeTrainDF.Rd index 858a664..0a54223 100644 --- a/man/computeTrainDF.Rd +++ b/man/computeTrainDF.Rd @@ -9,9 +9,9 @@ computeTrainDF(colData, formulaVars, tech, verbose = FALSE) \arguments{ \item{colData}{A per-cell metadata table. Typically `as.data.frame(colData(spe))`. Must include at least: -`cell_id`, raw metric columns named in `formulaVars` (e.g. +`cell_id`, raw metric columns named as in `formulaVars` (e.g. `log2SignalDensity`, `Area_um`, `log2Ctrl_total_ratio`, -optionally `log2AspectRatio`), and the corresponding outlier-label +optionally `log2AspectRatio` and `dist_border`), and the corresponding outlier-label columns referenced by `formulaVars`.} \item{formulaVars}{A named character vector mapping variable name to @@ -19,7 +19,7 @@ its outlier label column name, e.g. `c(log2SignalDensity="log2SignalDensity_outlier_train", ...)`.} \item{tech}{Character string with the acquisition technology. Used to -enable CosMx-specific handling for `log2AspectRatio`. Expected +enable CosMx-specific inclusion of border effect. Expected values include `"Nanostring_CosMx"` or `"Nanostring_CosMx_Protein"`.} \item{verbose}{`logical(1)` (default \code{FALSE}) @@ -27,12 +27,12 @@ If \code{TRUE}, prints the number of “bad” and “good” cells selected.} } \value{ A \code{data.frame} with one row per cell, including: - \code{qcscore_train} (0/1) indicating “bad” vs “good”, + \code{QScore_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. } \description{ -Build a Balanced Training Data Frame from a SpatialExperiment +Builds a balanced training data frame from a SpatialExperiment. \code{computeTrainDF} takes a \code{SpatialExperiment} object and assembles a balanced training set of “good” vs “bad” cells for @@ -44,10 +44,10 @@ in the `metadata` of the `SpatialExperiment` object. } \examples{ example(spatialPerCellQC) -spe <- computeOutliersQCScore(spe) +spe <- computeOutliersQScore(spe) spe <- checkOutliers(spe) df_train <- computeTrainDF(colData(spe), metadata(spe)$formula_variables, metadata(spe)$technology) -table(df_train$qcscore_train) +table(df_train$QScore_train) } diff --git a/man/dot-applyQScoreModel.Rd b/man/dot-applyQScoreModel.Rd new file mode 100644 index 0000000..b126f31 --- /dev/null +++ b/man/dot-applyQScoreModel.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/QC.R +\name{.applyQScoreModel} +\alias{.applyQScoreModel} +\title{.applyQScoreModel} +\usage{ +.applyQScoreModel(spe, qsModel, scoreName = "QScore") +} +\arguments{ +\item{spe}{A `SpatialExperiment` object with QC metrics already computed.} + +\item{qsModel}{A Quality Score model object, usually stored in +`metadata(spe)$QScore_model`.} + +\item{scoreName}{Name of the output column in `colData`.} +} +\value{ +A `SpatialExperiment` object with added Quality Score in `colData`. +} +\description{ +Internal: 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{ +The authors do not recommend applying a Quality Score (QS) model trained on one +dataset to a different dataset, but this could be useful in some cases, +for example to transfer a model trained on a dataset sourced from a tissue +slide sample and then applied to other slides from the same experiment. +When transferring the model from a reference dataset of the same technology +of the query dataset, there is no need to set a common formula. However, +when transferring from a different technology, reference QS computation +should be performed providing a formula that excludes untransferable or unwanted +variables (e.g. computing QS on a reference CosMx dataset omitting the border +effect term to transfer the QS model to a Xenium query dataset. This is done by +providing a custom formula in computeQScore, see link computeQScore). The function will check +for the presence of required model variables in the new dataset and will handle +missing or extra variables accordingly. +} +\examples{ + +example(readCosmxSPE) + +## Compute per-cell quality control metrics +spe <- spatialPerCellQC(spe) + +## Split the object only for example purposes +set.seed(1998) +idx <- sample(seq_len(ncol(spe)), floor(ncol(spe) / 2)) +spe_train <- spe[, idx] +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 + +## Apply the trained model to another dataset +spe_test <- applyQScoreModel( + spe=spe_test, + qcModel=qc_model, + scoreName="QScore_transferred" +) + +summary(spe_test$QScore_transferred) +} diff --git a/man/dot-checkSkw.Rd b/man/dot-checkSkw.Rd index a5bc0fb..865cdcd 100644 --- a/man/dot-checkSkw.Rd +++ b/man/dot-checkSkw.Rd @@ -15,7 +15,7 @@ \item{cd}{colData of `SpatialExperiment` object.} \item{metricList}{A character vector specifying the metrics to include in -the QC score formula. Defaults are "log2SignalDensity", "Area_um", +the Quality Score formula. Defaults are "log2SignalDensity", "Area_um", "log2AspectRatio", "log2Ctrl_total_ratio".} } \value{ @@ -23,6 +23,6 @@ A vector containing the list of chosen outlier detection method for each metric. } \description{ -Check skewness of metrics to choose outlier detection method. +Internal: checks skewness of metrics to choose outlier detection method. } \keyword{internal} diff --git a/man/dot-computeBorderDistanceCosMx.Rd b/man/dot-computeBorderDistanceCosMx.Rd index a5f27d8..a2f872a 100644 --- a/man/dot-computeBorderDistanceCosMx.Rd +++ b/man/dot-computeBorderDistanceCosMx.Rd @@ -13,16 +13,16 @@ \arguments{ \item{spe}{A `SpatialExperiment` object with CosMx data.} -\item{xwindim}{Width of FOV in x (default from `metadata(spe)$fov_dim`).} +\item{xwindim}{Width of FOV along x axis (default from `metadata(spe)$fov_dim`).} -\item{ywindim}{Height of FOV in y (default from `metadata(spe)$fov_dim`).} +\item{ywindim}{Height of FOV along y axis (default from `metadata(spe)$fov_dim`).} } \value{ -A `SpatialExperiment` object with `dist_border` columns in +A `SpatialExperiment` object with `dist_border` column in `colData`. } \description{ -Calculates the minimum distance of each cell to the field‐of‐view border -and adds it to `colData`. +Calculates the minimum distance of each cell to the Field‐Of‐View (FOV) +border and adds it to `colData`. } \keyword{internal} diff --git a/man/dot-computeCosmxProteinTrainSet.Rd b/man/dot-computeCosmxProteinTrainSet.Rd index 013259b..446c9ec 100644 --- a/man/dot-computeCosmxProteinTrainSet.Rd +++ b/man/dot-computeCosmxProteinTrainSet.Rd @@ -15,8 +15,8 @@ A list with elements \code{bad} and \code{good}, each a data.frame with \code{qscore_train} and (for “good”) an \code{is_a_bad_boy} flag. } \description{ -Internal: Build Training Set for CosMx-Protein +Internal: builds Training Set for CosMx-Protein. Splits a SpatialExperiment into “bad” vs “good” cells based on -outliers in aspect ratio near tissue border or low count area. +outliers in aspect ratio near FOV border or low signal density. } \keyword{internal} diff --git a/man/dot-computeCosmxTrainSet.Rd b/man/dot-computeCosmxTrainSet.Rd index e51f559..b5aae7e 100644 --- a/man/dot-computeCosmxTrainSet.Rd +++ b/man/dot-computeCosmxTrainSet.Rd @@ -15,8 +15,8 @@ A list with elements \code{bad} and \code{good}, each a data.frame with \code{qscore_train} and (for “good”) an \code{is_a_bad_boy} flag. } \description{ -Internal: Build Training Set for CosMx +Internal: builds Training Set for CosMx. Splits a SpatialExperiment into “bad” vs “good” cells based on -outliers in aspect ratio near tissue border or low count area. +outliers in aspect ratio near FOV border or low signal density. } \keyword{internal} diff --git a/man/dot-computeXenMerTrainSet.Rd b/man/dot-computeXenMerTrainSet.Rd index cfbed78..c524a74 100644 --- a/man/dot-computeXenMerTrainSet.Rd +++ b/man/dot-computeXenMerTrainSet.Rd @@ -15,7 +15,7 @@ A list with elements \code{bad} and \code{good}, each a data.frame with \code{qscore_train} and (for “good”) an \code{is_a_bad_boy} flag. } \description{ -Internal: Build Training Set for Xenium & MERFISH +Internal: builds Training Set for Xenium & MERFISH. Splits a SpatialExperiment into “bad” vs “good” cells based on pre-computed outlier labels on log2SignalDensity. } diff --git a/man/getModelFormula.Rd b/man/getModelFormula.Rd index e4379d7..49e67fd 100644 --- a/man/getModelFormula.Rd +++ b/man/getModelFormula.Rd @@ -4,15 +4,13 @@ \alias{getModelFormula} \title{getModelFormula} \usage{ -getModelFormula(formulaVars, verbose = FALSE) +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{verbose}{Logical. If `TRUE`, prints the final formula used for QC score} } \value{ `character` diff --git a/man/plotCellsFovs.Rd b/man/plotCellsFovs.Rd index b320190..d0377a7 100644 --- a/man/plotCellsFovs.Rd +++ b/man/plotCellsFovs.Rd @@ -8,11 +8,12 @@ plotCellsFovs( spe, sampleId = unique(spe$sample_id), pointCol = "firebrick", + pointSize = 0.05, + pointAlpha = 0.8, numbersCol = "black", - alphaNumbers = 0.8, + numberSize = 1, + numbersAlpha = 0.8, fovDim = metadata(spe)$fov_dim, - size = 0.05, - alpha = 0.8, scaleBar = TRUE, micronConvFact = 0.12 ) @@ -25,17 +26,19 @@ 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{alphaNumbers}{Numeric transparency for FoV labels. Default: `0.8`.} +\item{numberSize}{Numeric size for the FoV labels. Default: `1`.} + +\item{numbersAlpha}{Numeric transparency for FoV labels. Default: `0.8`.} \item{fovDim}{numeric with two named dimensions xdim, ydim. (Default is metadata(spe)$fov_dim)} -\item{size}{Numeric point size for the cell centroids. Default: `0.05`.} - -\item{alpha}{Numeric transparency for the cell centroids. Default: `0.8`.} - \item{scaleBar}{A logical value indicating whether to add a scale bar to the plot. (Default is `TRUE`)} @@ -55,7 +58,7 @@ metadata in metadata(spe): - `metadata(spe)$fov_positions`: a matrix or data.frame (or list with named elements) containing at minimum `x_global_px`, `y_global_px`, and `fov`. Values `x_global_px`/`y_global_px` are in pixels and represent -the origin (top-left) of each FoV. +the origin (bottom-left) of each FoV. - `metadata(spe)$fov_dim` (or the `fovDim` argument): a named numeric with `xdim` and `ydim` giving FoV width/height in pixels. diff --git a/man/plotZoomFovsMap.Rd b/man/plotZoomFovsMap.Rd index 7903d96..d19a97d 100644 --- a/man/plotZoomFovsMap.Rd +++ b/man/plotZoomFovsMap.Rd @@ -9,8 +9,11 @@ plotZoomFovsMap( fovs = NULL, title = NULL, mapPointCol = "darkmagenta", - mapNumbersCol = "black", - mapAlphaNumbers = 0.8, + mapPointSize = 0.5, + mapPointAlpha = 0.8, + fovNumbersCol = "black", + fovNumberSize = 1, + fovNumbersAlpha = 0.8, csize = 0.05, calpha = 0.8, scaleBars = NULL, @@ -32,17 +35,16 @@ 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{mapNumbersCol}{A character string specifying the color of the -numbers on the map. Default is `"black"`.} +\item{mapPointSize}{Numeric size for points in the map. Default: `0.5`.} -\item{mapAlphaNumbers}{A numeric value specifying the transparency of the -numbers on the map. Default is `0.8`.} +\item{mapPointAlpha}{Numeric transparency for points in the map. Default: `0.8`.} -\item{csize}{A numeric value specifying the size of the points in the map. -Default is `0.05`.} +\item{fovNumbersCol}{A character string specifying the color of the +numbers on the FoV zoom-in. Default is `"black"`.} -\item{calpha}{A numeric value specifying the transparency of the points in -the map. Default is `0.8`.} +\item{fovNumberSize}{Numeric size for the FoV labels. Default: `1`.} + +\item{fovNumbersAlpha}{Numeric transparency for FoV labels. Default: `0.8`.} \item{scaleBars}{Logical or NULL. Default is `NULL`. Master switch controlling the presence of scale bars in both panels. diff --git a/man/spatialPerCellQC.Rd b/man/spatialPerCellQC.Rd index 6c34445..82c343a 100644 --- a/man/spatialPerCellQC.Rd +++ b/man/spatialPerCellQC.Rd @@ -14,12 +14,13 @@ spatialPerCellQC( ) } \arguments{ -\item{spe}{A `SpatialExperiment` object containing spatial data.} +\item{spe}{A `SpatialExperiment` object containing spatial omics data.} \item{micronConvFact}{Numeric factor to convert pixels to microns. Default -`0.12`.} +is 0.12.} -\item{rmZeros}{logical for removing zero counts cells (default is TRUE).} +\item{rmZeros}{logical for removing zero total counts (protein intensities) +cells (default is TRUE).} \item{negProbList}{Character vector of patterns to identify negative probes. Defaults include: @@ -44,47 +45,57 @@ may subset the returned SpatialExperiment. Key behaviours and expectations: - Feature detection: negative‑probe patterns supplied in `negProbList` are -used to build `subsets_*` groups passed to `scater::addPerCellQC()` -(via `use_altexps` when requested). `addPerCellQC()` must be able to find -matching feature names (`rownames` of spe) and will create `subsets_*_sum` -and `subsets_*_detected` columns used below. -- Required columns: the function expects `sum`, `detected` and `total` -(from `addPerCellQC()` and the SPE assays) to be present; these are used to -compute `control_sum`, `control_detected`, `target_sum` and -`target_detected.` +passed to `scater::addPerCellQC()` (via `use_altexps` when requested). +`addPerCellQC()` must be able to find matching control feature names +(`rownames` of spe) and will create `subsets_*_sum` and `subsets_*_detected` +columns used below for each of control feature types. + +The function appends `sum`, `detected` and `total` to the SPE colData by +calling internally `addPerCellQC()`; these are used to compute `control_sum`, +`control_detected`, `target_sum` and `target_detected.` + - Control metrics: `control_sum` and `control_detected` are computed by summing matching `subsets_*` columns; `target_*` metrics are computed as -the complement vs sum / detected. -- Ratios and logs: `ctrl_total_ratio` (`control` / 9) and its stabilized -log2 transform `log2Ctrl_total_ratio` are added. +sum(detected) - control_sum(_detected). + +- Ratios and logs: `ctrl_total_ratio` (`control_sum`/`sum`) and its stabilized +log2-transformed `log2Ctrl_total_ratio` are added. + - Coordinate and area handling: - For CosMx technologies (Nanostring_CosMx and Nanostring_CosMx_Protein) spatial coordinates are converted from pixels to microns using `micronConvFact` and appended to `colData` (column names have px -> um). + - For CosMx, `Area_um` is derived from an existing Area column scaled by micronConvFact^2 and `.computeBorderDistanceCosMx()` is invoked to compute `dist_border`. + - For Nanostring_CosMx_Protein, a legacy `Area.um2` column (if present) is renamed to `Area_um` to standardize naming. + - For Xenium (10X_Xenium), if `Area_um` is missing the function will attempt to use `cell_area` as a fallback and issue a warning. + - Aspect ratio: if `AspectRatio` exists it is logged (`log2AspectRatio`); if missing a warning is emitted. + - Signal density: `SignalDensity` is computed as `sum` / `Area_um` for most -technologies; for Nanostring_CosMx_Protein it is set to total. A -log transform `log2SignalDensity` is also added. -- Zero‑count removal: when `rmZeros = TRUE` cells with `sum == 0` are -removed from the returned SpatialExperiment (message printed). +technologies; for Nanostring_CosMx_Protein it coincides with `total`. A +log2-transformed `log2SignalDensity` is also added. + +- Zero‑count (protein intensity) removal: when `rmZeros = TRUE` cells with +`sum == 0` are removed from the returned SpatialExperiment (message printed). + - Side effects: the function modifies `colData(spe)` (adds multiple new -columns), may add spatial coordinates into `colData` if missing, and may +columns). It may add spatial coordinates into `colData` if missing, and may subset the SPE to remove zero‑count cells. It issues warnings when expected inputs (e.g. area, aspect ratio, polygon‑derived fields) are missing or when fallbacks are used. Use this information to ensure the input SpatialExperiment contains the necessary assays and fields (feature names, `sum`/`detected`/`total`, -`Area`/`cell_area` when available) so metrics are computed and assigned -correctly. +`Area`/`cell_area`, log2AspectRatio, dist_border when available) so +metrics are computed and assigned correctly. } \examples{ example(readCosmxSPE) diff --git a/man/trainModel.Rd b/man/trainModel.Rd index 5b8d3a5..769a52a 100644 --- a/man/trainModel.Rd +++ b/man/trainModel.Rd @@ -20,7 +20,7 @@ A \code{\link[glmnet]{glmnet}} model object fitted with \eqn{\lambda} values. } \description{ -Fit a Ridge Logistic Regression Model +Fits a ridge logistic regression model. \code{trainModel} fits an L2-regularized (ridge) logistic regression using \pkg{glmnet}, given a design matrix and a training data frame. diff --git a/vignettes/Protein_tutorial_vignette.Rmd b/vignettes/Protein_tutorial_vignette.Rmd index e192fca..fa99299 100644 --- a/vignettes/Protein_tutorial_vignette.Rmd +++ b/vignettes/Protein_tutorial_vignette.Rmd @@ -18,7 +18,7 @@ editor_options: ```{r setup, include=FALSE} # Set chunk options: suppress echo, messages, and warnings in code output -knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE) +knitr::opts_chunk$set(echo=TRUE, message=FALSE, warning=FALSE) ``` ## Introduction @@ -49,7 +49,7 @@ library(SpaceTrooper) library(ggplot2) protfolder <- system.file( "extdata", "S01_prot", package="SpaceTrooper") -spe <- readCosmxProteinSPE(protfolder, sampleName = "CosMx_Protein_Tonsil") +spe <- readCosmxProteinSPE(protfolder, sampleName="CosMx_Protein_Tonsil") spe ``` @@ -66,7 +66,7 @@ any spatial shift, as such misalignment directly affects the QS computation. ```{r plot-fovs} # to check misalignment -plotCellsFovs(spe, size = 3, alpha = 0.7) +plotCellsFovs(spe, pointSize=3, pointAlpha=0.7, numberSize=6) ``` If any misalignment is observed, it can be corrected by adjusting the FoV @@ -84,7 +84,7 @@ Rerun the FoV plot to check whether the shift was corrected. ```{r plot-fovs-2} # check shift correction -plotCellsFovs(spe, size = 3, alpha = 0.7) +plotCellsFovs(spe, pointSize=3, pointAlpha=0.7, numberSize=6) ``` Cell centroids, shown in dark red, now are all contained by FoV boundaries, @@ -105,7 +105,7 @@ only for visualization purposes. They are stored as an `sf` object within This step is not mandatory for CosMx, because the pipeline can be executed even without them. -```{r load-poly, message = TRUE} +```{r load-poly, message=TRUE} # polygon loading spe <- readAndAddPolygonsToSPE(spe, boundariesType="csv") ``` @@ -132,7 +132,7 @@ custom vector, as shown below. By default, the function automatically removes 0 count cells, but this can be handled with the `rmZeros` parameter. -```{r cosmx-analysis-qc, message = TRUE} +```{r cosmx-analysis-qc, message=TRUE} spe <- spatialPerCellQC(spe, rmZeros=TRUE, negProbList=c("Ms IgG1", "Rb IgG")) @@ -158,7 +158,7 @@ the [SpaceTrooper utilities](https://bioconductor.org/packages/devel/bioc/vignet ### Compute Quality Score -`computeQCScore` function calculates QS per each cell. QS combines several +`computeQScore` function calculates QS per each cell. QS combines several metrics into a formula: **log2SignalDensity**, corresponding to signal density, **Area_um**, i.e. cell size, **log2Ctrl_total_ratio**, namely background signal and **log2AspectRatio** combined with **dist_border**, which jointly correspond @@ -173,30 +173,30 @@ formula contributes to cell quality. examples (fewer than 0.1% of the total number of cells), it is excluded from the formula and therefore not used in the QS computation. -The QS (stored as `QC_score` in `coldata`) ranges from 0 to 1, with 0 meaning +The QS (stored as `QScore` in `coldata`) ranges from 0 to 1, with 0 meaning low-quality and 1 high-quality. We are setting a seed to ensure reproducibility in this tutorial, because there are stochastic processes underlying QS computation. -```{r cosmx-analysis-score, message = TRUE} +```{r cosmx-analysis-score, message=TRUE} set.seed(1713) -spe <- computeQCScore(spe) +spe <- computeQScore(spe) -format(summary(spe$QC_score), scientific=FALSE, digits = 4) +format(summary(spe$QScore), scientific=FALSE, digits=4) ``` In this case, all the terms were used as no warning appeared. Subsequently, it is possible also to assess which cells have a QS lower than a certain threshold (default is 0.5) with the following function. It creates a -new column called `low_qcscore` inside of `colData`. +new column called `low_QScore` inside of `colData`. -```{r cosmx-analysis-score2, message = TRUE} -spe <- computeQCScoreFlags(spe, qsThreshold=0.5) +```{r cosmx-analysis-score2, message=TRUE} +spe <- computeQScoreFlags(spe, qsThreshold=0.5) -table(spe$low_qcscore) +table(spe$low_QScore) ``` Using this threshold, 274 cells are flagged as low-quality. We do not suggest a @@ -212,7 +212,7 @@ in handy. ```{r plot-hist} # view quantitative metric distribution -plotMetricHist(spe, metric = "QC_score") +plotMetricHist(spe, metric="QScore") ``` The QS distribution exhibits a left tail starting around 0.75. @@ -221,17 +221,17 @@ Cell visualization can be obtained by using either centroids (recommended when the dataset has a large number of cells) or polygons. `plotCentroids` plots cell centroids that can be colored by a certain metric -contained in the `colData` slot, by using the `colour_by` parameter. +contained in the `colData` slot, by using the `colourBy` parameter. Additionally, if you have a palette column in `colData`, containing colors for each cell, it can be given to `palette` parameter, so that it automatically -matches the column passed in `colour_by`. As an example, we are using the +matches the column passed in `colourBy`. As an example, we are using the cell types obtained as described in the [paper](https://www.biorxiv.org/content/10.64898/2025.12.24.696336v1) with their own color palette. ```{r plot-centroids-labels} labf <- system.file(file.path("extdata", "S01_prot", "labels_tiny.tsv"), package="SpaceTrooper") -labs <- read.table(file=labf, sep="\t", header=TRUE, comment.char = "") +labs <- read.table(file=labf, sep="\t", header=TRUE, comment.char="") spe$labels <- as.factor(labs[match(spe$cell_id, labs$cell_id),]$label) spe$labels_colors <- as.factor(labs[match(spe$cell_id, labs$cell_id),]$lab_color) @@ -265,11 +265,11 @@ Since all plotting functions are based on `ggplot2`, you can easily customize the graphical outputs by adding standard ggplot2 components. ```{r plot-polygons-fov-2} -plotPolygons(spe, colourBy="QC_score") + - scale_fill_viridis_c(option = "plasma") +plotPolygons(spe, colourBy="QScore") + + scale_fill_viridis_c(option="plasma") ``` -We can see that the QS is able to detect both the aspects highlighted by +We can see that the QS is able to detect all the aspects highlighted by `log2SignalDensity`,`Area_um`, `log2Ctrl_total_ratio` or `log2AspectRatio`. Cells that showed either lower signal density, bigger size, higher background signal or border effect also display @@ -279,24 +279,24 @@ It's up to the user to choose an appropriate threshold to flag cells according to the observed QS distribution. ```{r plot-polygons-fov-3} -plotPolygons(spe, colourBy="low_qcscore") + - scale_fill_manual(values=c("TRUE"="red", "FALSE" = "#c0c8cf")) +plotPolygons(spe, colourBy="low_QScore") + + scale_fill_manual(values=c("TRUE"="red", "FALSE"="#c0c8cf")) ``` -You can rerun`computeQCScoreFlags` to check how many cells would be flagged +You can rerun`computeQScoreFlags` to check how many cells would be flagged using another threshold. -```{r cosmx-analysis-score3, message = TRUE} -spe <- computeQCScoreFlags(spe, qsThreshold=0.75) +```{r cosmx-analysis-score3, message=TRUE} +spe <- computeQScoreFlags(spe, qsThreshold=0.75) -table(spe$low_qcscore) +table(spe$low_QScore) ``` Using this threshold, 527 cells are flagged as low-quality. ```{r plot-polygons-fov-4} -plotPolygons(spe, colourBy="low_qcscore") + - scale_fill_manual(values=c("TRUE"="red", "FALSE" = "#c0c8cf")) +plotPolygons(spe, colourBy="low_QScore") + + scale_fill_manual(values=c("TRUE"="red", "FALSE"="#c0c8cf")) ``` The threshold is more stringent and more cells are flagged compared to the diff --git a/vignettes/RNA_tutorial_vignette.rmd b/vignettes/RNA_tutorial_vignette.rmd index a195646..913e314 100644 --- a/vignettes/RNA_tutorial_vignette.rmd +++ b/vignettes/RNA_tutorial_vignette.rmd @@ -18,7 +18,7 @@ editor_options: ```{r setup, include=FALSE} # Set chunk options: suppress echo, messages, and warnings in code output -knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE) +knitr::opts_chunk$set(echo=TRUE, message=FALSE, warning=FALSE) ``` ## Introduction @@ -83,7 +83,7 @@ any spatial shift, as such misalignment directly affects the QS computation. ```{r plot-fovs} # to check misalignment -plotCellsFovs(spe, size = 3, alpha = 0.7) +plotCellsFovs(spe, pointSize=3, pointAlpha=0.7, numberSize=6) ``` The dataset is a subset with just a single FoV, whose number is displayed at the @@ -112,7 +112,7 @@ only for visualization purposes. They are stored as an `sf` object within This step is not mandatory for CosMx, because the pipeline can be executed even without them. -```{r load-poly, message = TRUE} +```{r load-poly, message=TRUE} # polygon loading spe <- readAndAddPolygonsToSPE(spe, boundariesType="csv") ``` @@ -139,7 +139,7 @@ custom vector, as shown below. By default, the function automatically removes 0 count cells, but this can be handled with the `rmZeros` parameter. -```{r cosmx-analysis-qc, message = TRUE} +```{r cosmx-analysis-qc, message=TRUE} spe <- spatialPerCellQC(spe, rmZeros=TRUE, negProbList=c("NegPrb", "Negative", "SystemControl")) @@ -166,7 +166,7 @@ the [SpaceTrooper utilities](https://bioconductor.org/packages/devel/bioc/vignet ### Compute Quality Score -`computeQCScore` function calculates QS per each cell. QS combines several +`computeQScore` function calculates QS per each cell. QS combines several metrics into a formula: **log2SignalDensity**, corresponding to signal density, **Area_um** or **volume**, i.e. cell size, **log2Ctrl_total_ratio**, namely background signal and **log2AspectRatio** combined with **dist_border**, which @@ -181,30 +181,30 @@ formula contributes to cell quality. examples (fewer than 0.1% of the total number of cells), it is excluded from the formula and therefore not used in the QS computation. -The QS (stored as `QC_score` in `coldata`) ranges from 0 to 1, with 0 meaning +The QS (stored as `QScore` in `coldata`) ranges from 0 to 1, with 0 meaning low-quality and 1 high-quality. We are setting a seed to ensure reproducibility in this tutorial, because there are stochastic processes underlying QS computation. -```{r cosmx-analysis-score, message = TRUE} +```{r cosmx-analysis-score, message=TRUE} set.seed(1713) -spe <- computeQCScore(spe) +spe <- computeQScore(spe) -summary(spe$QC_score) +summary(spe$QScore) ``` In this case, all the terms were used as no warning appeared. Subsequently, it is possible also to assess which cells have a QS lower than a certain threshold (default is 0.5) with the following function. It creates a -new column called `low_qcscore` inside of `colData`. +new column called `low_QScore` inside of `colData`. -```{r cosmx-analysis-score2, message = TRUE} -spe <- computeQCScoreFlags(spe, qsThreshold=0.5) +```{r cosmx-analysis-score2, message=TRUE} +spe <- computeQScoreFlags(spe, qsThreshold=0.5) -table(spe$low_qcscore) +table(spe$low_QScore) ``` Using this threshold, 97 cells are flagged as low-quality. We do not suggest a @@ -220,7 +220,7 @@ in handy. ```{r plot-hist} # view quantitative metric distribution -plotMetricHist(spe, metric = "QC_score") +plotMetricHist(spe, metric="QScore") ``` The QS distribution exhibits a left tail starting around 0.75. @@ -229,11 +229,11 @@ Cell visualization can be obtained by using either centroids (recommended when the dataset has a large number of cells) or polygons. `plotCentroids` plots cell centroids that can be colored by a certain metric -contained in the `colData` slot, by using the `colour_by` parameter. +contained in the `colData` slot, by using the `colourBy` parameter. Additionally, if you have a palette column in `colData`, containing colors for each cell, it can be given to `palette` parameter, so that it automatically -matches the column passed in `colour_by`. As an example, we are using the +matches the column passed in `colourBy`. As an example, we are using the cell types obtained as described in the [paper](https://www.biorxiv.org/content/10.64898/2025.12.24.696336v1) with their own color palette. ```{r plot-centroids-labels} @@ -273,11 +273,11 @@ Since all plotting functions are based on `ggplot2`, you can easily customize the graphical outputs by adding standard ggplot2 components. ```{r plot-polygons-fov-2} -plotPolygons(spe, colourBy="QC_score") + - scale_fill_viridis_c(option = "plasma") +plotPolygons(spe, colourBy="QScore") + + scale_fill_viridis_c(option="plasma") ``` -We can see that the QS is able to detect both the aspects highlighted by +We can see that the QS is able to detect all the aspects highlighted by `log2SignalDensity`,`Area_um`, `log2Ctrl_total_ratio` or `log2AspectRatio`. Cells that showed either lower signal density, bigger size, higher background signal or border effect also display @@ -287,24 +287,24 @@ It's up to the user to choose an appropriate threshold to flag cells according to the observed QS distribution. ```{r plot-polygons-fov-3} -plotPolygons(spe, colourBy="low_qcscore") + - scale_fill_manual(values=c("TRUE"="red", "FALSE" = "#c0c8cf")) +plotPolygons(spe, colourBy="low_QScore") + + scale_fill_manual(values=c("TRUE"="red", "FALSE"="#c0c8cf")) ``` -You can rerun`computeQCScoreFlags` to check how many cells would be flagged +You can rerun`computeQScoreFlags` to check how many cells would be flagged using another threshold. -```{r cosmx-analysis-score3, message = TRUE} -spe <- computeQCScoreFlags(spe, qsThreshold=0.75) +```{r cosmx-analysis-score3, message=TRUE} +spe <- computeQScoreFlags(spe, qsThreshold=0.75) -table(spe$low_qcscore) +table(spe$low_QScore) ``` Using this threshold, 193 cells are flagged as low-quality. ```{r plot-polygons-fov-4} -plotPolygons(spe, colourBy="low_qcscore") + - scale_fill_manual(values=c("TRUE"="red", "FALSE" = "#c0c8cf")) +plotPolygons(spe, colourBy="low_QScore") + + scale_fill_manual(values=c("TRUE"="red", "FALSE"="#c0c8cf")) ``` The threshold is more stringent and more cells are flagged compared to the diff --git a/vignettes/SpaceTrooper_utilities.Rmd b/vignettes/SpaceTrooper_utilities.Rmd index d44e414..8de8a96 100644 --- a/vignettes/SpaceTrooper_utilities.Rmd +++ b/vignettes/SpaceTrooper_utilities.Rmd @@ -18,7 +18,7 @@ editor_options: ```{r setup, include=FALSE} # Set chunk options: suppress echo, messages, and warnings in code output -knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE) +knitr::opts_chunk$set(echo=TRUE, message=FALSE, warning=FALSE) ``` ## Introduction @@ -69,7 +69,7 @@ spe <- readCosmxSPE(dirName=cosmxFolder, metadataFPattern="metadata_file.csv", polygonsFPattern="polygons.csv", fovPosFPattern="fov_positions_file.csv", - fovdims=c(xdim = 4256, ydim = 4256) + fovdims=c(xdim=4256, ydim=4256) ) ``` @@ -126,7 +126,7 @@ prot <- readCosmxProteinSPE(dirName=protfolder, metadataFPattern="metadata_file.csv", polygonsFPattern="polygons.csv", fovPosFPattern="fov_positions_file.csv", - fovdims=c(xdim = 4256, ydim = 4256) + fovdims=c(xdim=4256, ydim=4256) ) # code line for shift correction @@ -141,7 +141,7 @@ corresponding file name patterns. **XENIUM** -```{r load-xenium, message = TRUE, warning = TRUE} +```{r load-xenium, message=TRUE, warning=TRUE} xeniumFolder <- system.file( "extdata", "Xenium_small", package="SpaceTrooper") xen <- readXeniumSPE(dirName=xeniumFolder, @@ -196,7 +196,7 @@ file is available, you should provide a vector containing only that format. `computeMissingMetrics` is set to TRUE by default and enables the computation of `Area_um` and `AspectRatio` cell metrics from boundary polygons. **IMPORTANT**: it is strongly recommended to keep this parameter set to `TRUE` -as both metrics are required to run `computeQCScore()`, even if they are not used +as both metrics are required to run `computeQScore()`, even if they are not used in the final formula. `keepPolygons` controls whether polygons are retained in the `SpatialExperiment` @@ -225,9 +225,9 @@ information is extracted. This option is not mandatory and is therefore set to `SpatialExperimentIO::readXeniumSXE`, and multiple probe identifier sets could be available, Gene Symbols are selected and used as row names. -**MERFISH** +**MERFISH/MERSCOPE** -```{r load-merfish, message = TRUE, warning= TRUE} +```{r load-merfish, message=TRUE, warning= TRUE} merfishFolder <- system.file("extdata", "Merfish_Tiny", package="SpaceTrooper") mer <- readMerfishSPE(dirName=merfishFolder, @@ -243,7 +243,7 @@ mer <- readMerfishSPE(dirName=merfishFolder, ``` As with the other data-reading functions, `readMerfishSPE` requires the path to -the directory containing the MERFISH output files. The minimum set of inputs +the directory containing the MERFISH output files. The minimum set of inputs required to run the `SpaceTrooper` pipeline includes: - expression matrix @@ -274,7 +274,7 @@ That said, cell volume appears to be a more appropriate measure of cell size for QS computation with MERFISH data. In the current version of `SpaceTrooper`, when `computeMissingMetrics=TRUE`, the -`volume` is renamed as `Area_um` for compatibility with the `computeQCScore` +`volume` is renamed as `Area_um` for compatibility with the `computeQScore` function, whereas cell area is not computed. ### Load polygons @@ -292,7 +292,7 @@ associated with them; - stores the resulting `sf` object in a `colData` column whose name can be specified via the `polygonsCol` argument; -```{r load-poly, message = TRUE} +```{r load-poly, message=TRUE} # polygon loading spe <- readAndAddPolygonsToSPE(spe, polygonsCol="polygons", @@ -312,17 +312,17 @@ By default, the `keepMultiPol` parameter is set to `TRUE`. Under this setting, cells associated with multi-polygons are retained in the dataset; however, their `AspectRatio` is `NA`. -```{r load-poly-xen, message = TRUE} +```{r load-poly-xen, message=TRUE} # xenium polygon loading xen <- readAndAddPolygonsToSPE(xen, boundariesType=c("parquet")) ``` -```{r load-poly-mer, message = TRUE} +```{r load-poly-mer, message=TRUE} # merfish polygon loading mer <- readAndAddPolygonsToSPE(mer, boundariesType=c("parquet")) ``` -```{r load-poly-prot, message = TRUE} +```{r load-poly-prot, message=TRUE} # protein polygon loading prot <- readAndAddPolygonsToSPE(prot, boundariesType=c("csv")) ``` @@ -417,17 +417,17 @@ involves underlying stochastic processes. **IMPORTANT**: when working with recent MERSCOPE datasets, or with any dataset that uses numeric cell identifiers, it is strongly recommended to convert cell -IDs to character strings before calling the `computeQCScore` function. This +IDs to character strings before calling the `computeQScore` function. This prevents errors during the selection of training examples when estimating the coefficients of the QS formula. -```{r compute-QS, message = TRUE} +```{r compute-QS, message=TRUE} # MERSCOPE cell_id conversion # spe$cell_id <- as.character(spe$cell_id) set.seed(1713) -spe <- computeQCScore(spe, verbose=TRUE) +spe <- computeQScore(spe, verbose=TRUE) ``` When `verbose=TRUE`, additional information regarding the underlying steps is @@ -471,40 +471,80 @@ 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`QC_score` 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} +```{r compute-QS-safe-run, message=TRUE} # safe run function safe_run <- function(expr) { tryCatch( - list(result = expr, error = NULL), - error = function(e) list(result = NULL, error = e) + list(result=expr, error=NULL), + error=function(e) list(result=NULL, error=e) ) } -out <- safe_run(computeQCScore(spe)) +out <- safe_run(computeQScore(spe)) if (!is.null(out$error)) { message("Failed: ", out$error$message) -colData(spe)$QC_score <- NA +colData(spe)$QScore <- NA } else { 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. +The default is `NULL` and the formula is automatically estimated based on the +available metrics, their outliers and the dataset technology. The most complete +formula is generated as follows: + + `~(log2SignalDensity + Area_um + I(abs(log2AspectRatio) * as.numeric(dist_border < 50)) + log2Ctrl_total_ratio)^2` + + When user-provided, the formula must follow the same default syntax and the terms + must be written exactly as they appear here (e.g. `I(abs(log2AspectRatio) * as.numeric(dist_border < 50))` + must have spaces around the `*` and `<` operators). However, the terms do not need + to follow any specific order. If the outliers for one or more metrics are not + sufficient, they will not be considered for model training and the corresponding terms + are automatically dropped from the formula. + We provide interaction terms by default as they contribute to model flexibility. + However, if users prefer to exclude interaction terms, they can do so by + removing the `()^2` notation from the formula. + +```{r compute-QS-different-formula, message=TRUE} + +set.seed(1713) + +model_formula <- "~log2Ctrl_total_ratio + log2SignalDensity" + +spe <- computeQScore(spe, modelFormula=model_formula, verbose=TRUE) +``` + +```{r plot-polygons-dbkero} +plotPolygons(spe, colourBy="QScore") + + scale_fill_viridis_c(option="plasma") +``` + +Compared to the plot displayed in [RNA](https://bioconductor.org/packages/devel/bioc/vignettes/SpaceTrooper/inst/doc/RNA_tutorial_vignette.html) vignette, +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. + ### Single metric flagging The `SpaceTrooper` package allows the identification of numerically aberrant cells through outlier detection for each quantitative metric, using the same -procedures applied internally by the `computeQCScore` function. This enables +procedures applied internally by the `computeQScore` function. This enables metric-specific cell flagging in a manner analogous to standard scRNA-seq quality control workflows. -```{r compute-outliers-1, message = TRUE, warning = TRUE} +```{r compute-outliers-1, message=TRUE, warning=TRUE} # detecting outliers for cell area in Xenium dataset with both Medcouple # and MAD -xen <- computeSpatialOutlier(xen, computeBy="Area_um", method = "both") +xen <- computeSpatialOutlier(xen, computeBy="Area_um", method="both") ``` Pay attention to warnings: setting the method to `both`, applies both the @@ -514,8 +554,8 @@ and `_outlier_sc` to the variable name. If a warning indicates that the variable is symmetric, only the `_outlier_sc` column should be used for downstream analyses. -```{r compute-outliers-3, message = TRUE} -spe <- computeSpatialOutlier(spe, computeBy="Mean.CD68", method = "both") +```{r compute-outliers-3, message=TRUE} +spe <- computeSpatialOutlier(spe, computeBy="Mean.CD68", method="both") table(spe$Mean.CD68_outlier_mc) ``` @@ -534,7 +574,7 @@ corresponding to the upper and lower thresholds, used to define `HIGH` and `LOW` outliers, respectively. ```{r plot-hist-1} -spe <- computeSpatialOutlier(spe, computeBy="log2SignalDensity", method = "both") +spe <- computeSpatialOutlier(spe, computeBy="log2SignalDensity", method="both") plotMetricHist(spe, metric="log2SignalDensity", useFences="log2SignalDensity_outlier_mc") @@ -549,7 +589,8 @@ by `plotCellsFovs` function. ```{r plot-zoom-fov-1} # plotting FoVs map and zoom in of selected FoV colored by QS -plotZoomFovsMap(spe, fovs=16, colourBy="QC_score", csize=3, scaleBars=TRUE) +plotZoomFovsMap(spe, fovs=16, colourBy="QScore", mapPointSize=1, + fovNumberSize=7, scaleBars=TRUE) ``` For this plot, `scaleBars` parameter controls whether scale bars are displayed @@ -557,17 +598,17 @@ in both plots simultaneously. To show the scale bar in only one of the two plots, set either `scaleBarMap=FALSE` or `scaleBarPol=FALSE`. ```{r plot-zoom-fov-2} -plotZoomFovsMap(prot, fovs = 344, colourBy="Area_um", csize=3, scaleBarMap=TRUE, - scaleBarPol=FALSE) +plotZoomFovsMap(prot, fovs=344, colourBy="Area_um", mapPointSize=1, + fovNumberSize=7, scaleBarMap=TRUE, scaleBarPol=FALSE) ``` For the other plots, the `scaleBar` parameter can be set to `FALSE` to hide the scale bar. ```{r plot-centroids} -plotCentroids(spe, colourBy="QC_score", size=3, scaleBar=FALSE) + - scale_fill_viridis_c(option = "plasma") + - scale_color_viridis_c(option = "plasma") +plotCentroids(spe, colourBy="QScore", size=3, scaleBar=FALSE) + + scale_fill_viridis_c(option="plasma") + + scale_color_viridis_c(option="plasma") ``` ```{r plot-polygons}