---
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 = -21, y = 0, label = "endogenous amount",
               hjust = 0, vjust = 2, size = 3, 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**: `create_cal_curve()` 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 bracket the signal actually measured in tissue
  (checked by `plot_cal_coverage()`);
- the fitted model is adequate, with acceptable residual behaviour;
- the calibration substrate is a reasonable matrix match for the tissue;
- any internal standard behaves consistently across the acquisition;
- pixel dimensions are known, for conversion to pg/mm².

Keeping the analytical steps distinct from the interpretive one:

1. raw intensity;
2. internal-standard-normalised response (`int2response()`);
3. background-referenced response / SNR (`int2snr()`);
4. calibrated amount estimate (`int2conc()`);
5. biological interpretation.

Steps 1--4 are analytical transformations performed by the package. Step 5
remains conditional on the assumptions above.

# 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. `summarise_cal_levels()` -- average each calibration level per standard spot,
2. `create_cal_curve()` -- fit the response-vs-amount model per analyte,
3. `int2conc()` -- apply those models to convert tissue-pixel responses to
   estimated amounts,

then `plot_cal_coverage()` to confirm the standards actually bracket the signal
measured in tissue.

*Family: `calibration`.*

| Function | Takes | Produces |
|---|---|---|
| `summarise_cal_levels()` | object with `Cal` pixels + `identifier`; `cal_metadata` (`identifier`, `analyte`, `amount_pg`, `level`) | fills `cal_response_data` and the per-level pixel counts |
| `create_cal_curve()` | 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 |
| `plot_cal_coverage()` | 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. `run_example(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

`summarise_cal_levels()` averages the signal at each calibration spot, and
`create_cal_curve()` 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 <- summarise_cal_levels(cal, cal_metadata, val_slot = "intensity",
                            cal_label = "Cal", id = "identifier")
cal <- create_cal_curve(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.

## R² is not enough

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
  (`plot_cal_coverage()`, 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. `create_cal_curve()` 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. `plot_cal_coverage()` 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."}
plot_cal_coverage(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 `plot_cal_coverage()` 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 `<study>_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}
run_study("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).
<https://doi.org/10.1021/acs.analchem.4c02350>

# Session information

```{r session}
sessionInfo()
```