--- title: "Quantification with calibration standards" author: "Matthew J. Smith" date: "`r Sys.Date()`" package: quantMSImageR output: BiocStyle::html_document: toc_float: true vignette: > %\VignetteIndexEntry{2. Quantification with calibration standards} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE) # The schematics in the next section are drawn before the worked example # attaches its packages, so ggplot2 is needed this early. library(ggplot2) ``` # Introduction Mass spectrometry imaging measures analyte response at each pixel rather than analyte amount directly. Using calibration standards, `quantMSImageR` can convert this response into an estimated amount per pixel (pg/pixel) and, when pixel dimensions are known, an estimated real amount (pg/mm²). These estimates are conditional on the calibration design, internal standard choice and strategy, matrix matching, acquisition conditions, and validity of the fitted model. ## On-slide versus on-tissue calibration ```{r cal-schematic-helper, echo = FALSE} # Two small schematics per calibration design: where the standards sit, and how # amount is recovered from the fit. Illustrative numbers throughout -- these are # diagrams, not measurements. .slide_plot <- function(on_tissue) { tis <- data.frame(xmin = 0.4, xmax = 3.6, ymin = 0.6, ymax = 3.4) sp <- if (on_tissue) expand.grid(x = c(0.9, 1.7, 2.5, 3.3), y = c(1.2, 2.0, 2.8)) else expand.grid(x = c(4.3, 4.9, 5.5), y = c(1.2, 2.0, 2.8)) ggplot() + geom_rect(data = tis, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), fill = "#7db954", colour = "grey40", alpha = 0.55) + geom_point(data = sp, aes(x, y), shape = 21, size = 4, fill = "#e8351e", colour = "grey20", stroke = 0.4) + annotate("text", x = 2.0, y = 0.25, label = "tissue section", size = 3, colour = "grey25") + annotate("text", x = if (on_tissue) 2.1 else 4.9, y = if (on_tissue) 3.75 else 0.25, label = "standards", size = 3, colour = "#cd023d") + coord_fixed(xlim = c(0, 6), ylim = c(0, 4)) + theme_void(base_size = 11) + theme(panel.border = element_rect(colour = "grey70", fill = NA)) } .fit_plot <- function(on_tissue) { add <- c(0, 10, 20, 30, 40) b0 <- if (on_tissue) 45 else 5 # response at zero added amount pts <- data.frame(amount = add, response = b0 + 2 * add) g <- ggplot() + geom_hline(yintercept = 0, colour = "grey30", linewidth = 0.4) + geom_vline(xintercept = 0, colour = "grey30", linewidth = 0.4) + geom_line(data = data.frame(x = c(0, 45), y = b0 + 2 * c(0, 45)), aes(x, y), colour = "#0f8096", linewidth = 0.7) + geom_point(data = pts, aes(amount, response), shape = 21, size = 2.6, fill = "#e8351e", colour = "grey20") if (on_tissue) { # The extrapolation back to the x-intercept is what standard addition adds. g <- g + geom_line(data = data.frame(x = c(-25, 0), y = b0 + 2 * c(-25, 0)), aes(x, y), colour = "#0f8096", linewidth = 0.7, linetype = "dashed") + geom_point(aes(x = -22.5, y = 0), shape = 4, size = 3.4, stroke = 1.1, colour = "#cd023d") + annotate("text", x = -22.5, y = 16, label = "endogenous amount", hjust = 0.5, vjust = 0, size = 2.9, lineheight = 0.9, colour = "#cd023d") } g + scale_x_continuous(limits = if (on_tissue) c(-30, 48) else c(-5, 48)) + labs(x = "amount deposited (pg / pixel)", y = "response") + theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) } ``` On-slide (external) calibration deposits standards on the slide adjacent to the tissue. This captures instrumental response during the same acquisition but does not reproduce tissue-specific matrix effects. Use `cal_type = "cal"`. ```{r cal-onslide, echo = FALSE, fig.wide = TRUE, fig.height = 2.4, fig.cap = "On-slide calibration. Standards are deposited on the slide beside the section, so they share the acquisition but not the tissue matrix. Amount is read straight off the fitted line, whose intercept is instrument background. Schematic."} patchwork::wrap_plots(.slide_plot(FALSE), .fit_plot(FALSE), widths = c(1, 1.1)) ``` On-tissue calibration deposits standards onto tissue, or onto a representative tissue homogenate. Matrix matching is better, but endogenous analyte signal and spatial heterogeneity still need careful treatment. Use `cal_type = "std_addition"`. Standard addition: `createCalCurve()` fits response against deposited amount, takes the x-intercept of that line as the endogenous amount already present in the tissue, shifts every deposited amount by it so the axis expresses *total* amount (endogenous + added), drops the background level and refits. The estimate therefore accounts for analyte the tissue already contained. ```{r cal-ontissue, echo = FALSE, fig.wide = TRUE, fig.height = 2.4, fig.cap = "On-tissue calibration by standard addition. Standards are deposited onto the section, so they experience the same matrix as the analyte. The response at zero added amount already contains endogenous analyte, so the line is extrapolated back (dashed) to the x-intercept, and that amount is added to every deposited amount. Schematic."} patchwork::wrap_plots(.slide_plot(TRUE), .fit_plot(TRUE), widths = c(1, 1.1)) ``` `cal_type` has no default: the two modes do materially different things, so it must be chosen explicitly. ## Scope and assumptions Amount estimates are only meaningful where the following hold. Each is the user's responsibility: - the standard is the same chemical species as the endogenous analyte; - transitions and acquisition settings match between standards and sample; - calibration levels expand the response range of the signal actually measured in tissue (checked by `plotCalCoverage()`); - the fitted model is adequate, with acceptable residual behaviour; - the calibration substrate is a reasonable matrix match for the tissue; - any internal standard is analogous to the endogenous analyte and is consistently measured throughout the acquisition; - pixel dimensions are known, for conversion to pg/mm². # Installation ```{r install, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("quantMSImageR") ``` # Calibration workflow overview When a calibration standards acquisition is available, `quantMSImageR` fits a response-versus-amount model per analyte and converts tissue-pixel responses to `pg_pixel` and `pg_mm2`. The workflow is three steps: 1. `summariseCalLevels()` -- average each calibration level per standard spot, 2. `createCalCurve()` -- fit the response-vs-amount model per analyte, 3. `int2conc()` -- apply those models to convert tissue-pixel responses to estimated amounts, then `plotCalCoverage()` to confirm the standards actually bracket the signal measured in tissue. *Family: `calibration`.* | Function | Takes | Produces | |---|---|---| | `summariseCalLevels()` | object with `Cal` pixels + `identifier`; `cal_metadata` (`identifier`, `analyte`, `amount_pg`, `level`) | fills `cal_response_data` and the per-level pixel counts | | `createCalCurve()` | that object; `cal_type = "cal"` or `"std_addition"` | fills `cal_list` (one linear model per analyte) and `r2_df` | | `int2conc()` | object carrying `cal_list`; `pixel_header` / `pixels`; `experimentData$pixelSize` for the mm² layer | adds `pg_pixel` and `pg_mm2`; features without a curve are dropped | | `plotCalCoverage()` | a calibrated object (after `int2conc()`) | a `patchwork` figure: pixel histogram above the calibration curve | This vignette demonstrates the workflow on a small synthetic calibration dataset bundled with the package. It applies the resulting models to a study section, and then shows how the same steps are driven from a study YAML configuration. `runExample(calibrate = TRUE)` runs the chain end-to-end. The core imaging workflow that produces the object these operate on -- SNR filtering, ion images, heatmaps -- is covered in the *Introduction to quantMSImageR* vignette. The approach is described and applied in Smith *et al.* (2024); see *References*. # The bundled calibration dataset The synthetic calibration object (generated by `inst/scripts/generate_cal_data.R`) is a `quant_MSImagingExperiment` whose pixels are labelled `Cal` (standard spots), `Tissue` or `Background` in `sample_type`, with a unique `identifier` per standard spot. Its analytes are named to match features in the imaging example, so the fitted curves can be applied to a study section later. ```{r load} library(quantMSImageR) library(ggplot2) cal_dir <- system.file("extdata", "cal_example.raw", package = "quantMSImageR") cal <- as(readRDS(file.path(cal_dir, "cal_MSI.RDS")), "quant_MSImagingExperiment") cal_metadata <- read.csv(file.path(cal_dir, "calibration_metadata.csv")) pix <- as.data.frame(table(pData(cal)$sample_type)) names(pix) <- c("sample_type", "n_pixels") ``` ```{r pix-table} DT::datatable(pix, rownames = FALSE, options = list(dom = "t", scrollX = TRUE)) ``` **Table 1.** How the calibration acquisition's pixels are labelled. `Cal` are the standard spots used to build the curves; `Tissue` and `Background` are the surrounding sample and background. The companion metadata table is the link between a physical standard spot and what was deposited on it: each `identifier` maps to an `analyte`, the amount spotted (`amount_pg`) and which dilution `level` it belongs to. ```{r cal-meta} DT::datatable(cal_metadata, rownames = FALSE, options = list(pageLength = 10, scrollX = TRUE)) ``` **Table 2.** Calibration metadata: amount of each analyte deposited at every standard spot. # Fitting and diagnosing calibration models `summariseCalLevels()` averages the signal at each calibration spot, and `createCalCurve()` fits a linear response-vs-amount model per analyte (here `cal_type = "cal"` for on-slide standards; use `"std_addition"` for standard addition on tissue). ```{r fit} cal <- summariseCalLevels(cal, cal_metadata, val_slot = "intensity", cal_label = "Cal", id = "identifier") cal <- createCalCurve(cal, cal_type = "cal") ``` Each fit is a straight line, `response = slope x (pg/pixel) + intercept`, which `int2conc()` later inverts to go from response back to amount. The slope is the analyte's sensitivity. ```{r eqn} cal_list <- calibrationModels(cal) r2_df <- calibrationDiagnostics(cal) eqn <- data.frame( analyte = names(cal_list), slope = vapply(cal_list, function(m) unname(coef(m)[2]), numeric(1)), intercept = vapply(cal_list, function(m) unname(coef(m)[1]), numeric(1)), row.names = NULL ) eqn$r2 <- r2_df$r2[match(eqn$analyte, r2_df$feature)] ``` ```{r eqn-table} DT::datatable(eqn, rownames = FALSE, options = list(dom = "t", scrollX = TRUE)) |> DT::formatSignif(c("slope", "intercept", "r2"), digits = 4) ``` **Table 3.** Fitted calibration model per analyte: the slope, intercept and R-squared of the line that `int2conc()` inverts. Whilst a useful metric R² reports only how much of the variance in the standards the straight line accounts for. Treat a low R² as a definite warning and a high R² as necessary but not sufficient. Before relying on the estimates, examine at minimum: - the calibration range, and whether the tissue signal falls inside it (`plotCalCoverage()`, below); - residual patterns against fitted amount; - precision across calibration replicates; - the lower and upper quantifiable limits implied by the design; - whether the fit needs weighting. `createCalCurve()` currently applies inverse-amount (`1/amount`) weighting. This reduces the influence of the highest calibration levels, but its suitability should be assessed from the residual structure of your analytical method. Plotting those fits gives the calibration curves. The line must come from the stored model: ```{r curve-plot, fig.wide = TRUE, fig.cap = "Calibration curves. Each point is one standard spot's summarised response; the line is predicted from the stored, weighted model that int2conc() inverts."} cal_data <- calibrationLevels(cal) models <- calibrationModels(cal) # Predict from each stored model across the amounts it was fitted to. cal_lines <- lapply(names(models), function(analyte) { m <- models[[analyte]] xx <- seq(min(m$model$pg_perpixel), max(m$model$pg_perpixel), length.out = 100) data.frame(analyte = analyte, pg_perpixel = xx, response_perpixel = predict(m, newdata = data.frame(pg_perpixel = xx))) }) cal_lines <- do.call(rbind, cal_lines) ggplot(cal_data, aes(pg_perpixel, response_perpixel)) + geom_point() + geom_line(data = cal_lines, colour = "steelblue", linewidth = 0.8) + facet_wrap(~ analyte, scales = "free") + labs(x = "amount (pg / pixel)", y = "response per pixel") ``` # Applying the models to a study section The models are fitted on the standards, then carried across to the study data. Here we load one section of the imaging example (the circular tissue of group A) and attach the fitted models before converting its tissue pixels: ```{r apply} study <- as(readRDS(system.file("extdata", "example.raw", "section01.RDS", package = "quantMSImageR")), "quant_MSImagingExperiment") # carry the fitted curves over, then convert this section's tissue pixels calibrationData(study) <- calibrationData(cal) study <- int2conc(study, val_slot = "intensity", pixel_header = "sample_name", pixels = "tissue_pixels") analyte <- fData(study)$name[1] analyte ``` ## Checking calibration coverage Converting a response into an amount estimate is only trustworthy where the standards constrain the fit. `plotCalCoverage()` enables this to be evaluated. The lower panel is the calibration itself -- the standard spots (red) and the fitted model (dashed) that `int2conc()` inverts. The upper panel is the distribution of the section's own pixels on exactly the same x-axis, with the calibration levels marked in red. When the histogram sits inside the span of the red lines, every converted pixel is interpolated between real standards. ```{r dyn-range, fig.wide = TRUE, fig.height = 5, fig.cap = "Calibration coverage of the measured data. Pixels falling inside the marked range are interpolated between real standards rather than extrapolated."} plotCalCoverage(study, features = analyte, val_slot = "intensity") ``` Called without `features` it draws one facet per calibrated analyte, which is the usual way to screen a whole panel at once for standards that do not bracket the data. ## Handling out-of-range pixels Either the calibration needs re-acquiring such that tissue pixels are interpolated (recheck with `plotCalCoverage()` as above). Or a small % of extrapolation of extreme pixels is acceptable for your purpose and is set by `max_out_of_range` (`1` accepts any). The proportion is recorded per feature, so it can be carried into a report alongside R-squared: ```{r oor} calibrationDiagnostics(study) ``` Here the synthetic standards span 1.25 to 100 pg/pixel and every tissue pixel falls inside that, so `out_of_range` is zero throughout. The same section can now be imaged either way. First in raw response, where the colour scale is arbitrary instrument signal: ```{r img-intensity, fig.small = TRUE, fig.cap = "Study section in raw intensity. The colour scale is arbitrary instrument response and is not comparable to any other acquisition."} imageR(study, feat_ind = 1, sample_lab = "run", val_slot = "intensity", scale = "suppress") ``` Then the same data in calibrated units, where the colour scale is an estimated amount per unit area: ```{r img-conc, fig.small = TRUE, fig.cap = "The same section after int2conc(). The colour scale is now an estimated amount per unit area rather than instrument response, which is interpretable on its own terms and comparable across acquisitions where standards, matrix and acquisition conditions are matched."} imageR(study, feat_ind = 1, sample_lab = "run", val_slot = "pg_mm2", value = "pg/mm2", scale = "suppress") ``` # Driving calibration from a study YAML In a full study, calibration is enabled through the `calibration:` block of the configuration file (`system.file("config_template.yaml", package = "quantMSImageR")`), positioned just before the `ratios:` block: ```yaml calibration: enabled: true cal_acquisition: "cal_run_slide1" # standards .raw acquisition cal_roi_csv: "cal_ROIs.csv" # pixel -> sample_type = "Cal" + identifier cal_metadata: "calibration_metadata.csv" cal_type: "std_addition" # std_addition | cal val_slot: "intensity" background_level: "background" quantify_pixels: ["tissue_pixels"] ``` With `enabled: true`, the study runner builds the models from the standards acquisition and writes a `_calibrated.RDS` alongside the reports. The rendered report gains a *Quantification* section carrying the fits, the coverage plot and calibrated ion images: ```{r run-study, eval = FALSE} runStudy("path/to/study_config.yaml") ``` # Interpretation and limitations The values produced here are calibrated amount estimates, not volumetric concentrations and not traceable reference measurements. Reporting them as `pg/pixel` or `pg/mm²` is a statement about the calibration, not a claim about the absolute biological content of a voxel. External or on-slide calibration in particular does not automatically account for: - tissue-specific ion suppression; - regional variation in desorption and extraction efficiency; - differences in section thickness; - analyte recovery; - source geometry and day-to-day instrumental changes; - mismatch between the calibration substrate and the biological tissue. On-tissue calibration improves matrix matching but does not remove these concerns, and adds the problem of separating added from endogenous analyte. In practice this means amount estimates are most defensible when used to compare like with like -- the same analyte, prepared and acquired the same way, with standards that bracket the observed signal -- and progressively weaker as any of those conditions relax. Where a comparison matters, report the calibration diagnostics (range, coverage, R², out-of-range proportion) alongside the values. # References Smith MJ, Nie M, Adner M, Säfholm J, Wheelock CE. Development of a Desorption Electrospray Ionization–Multiple-Reaction-Monitoring Mass Spectrometry (DESI-MRM) Workflow for Spatially Mapping Oxylipins in Pulmonary Tissue. *Analytical Chemistry* (2024). # Session information ```{r session} sessionInfo() ```