## ----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)

## ----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"))
}

## ----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))

## ----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))

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("quantMSImageR")

## ----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")

## ----pix-table----------------------------------------------------------------
DT::datatable(pix, rownames = FALSE,
              options = list(dom = "t", scrollX = TRUE))

## ----cal-meta-----------------------------------------------------------------
DT::datatable(cal_metadata, rownames = FALSE,
              options = list(pageLength = 10, scrollX = TRUE))

## ----fit----------------------------------------------------------------------
cal <- summariseCalLevels(cal, cal_metadata, val_slot = "intensity",
                            cal_label = "Cal", id = "identifier")
cal <- createCalCurve(cal, cal_type = "cal")

## ----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)]

## ----eqn-table----------------------------------------------------------------
DT::datatable(eqn, rownames = FALSE,
              options = list(dom = "t", scrollX = TRUE)) |>
  DT::formatSignif(c("slope", "intercept", "r2"), digits = 4)

## ----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")

## ----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

## ----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")

## ----oor----------------------------------------------------------------------
calibrationDiagnostics(study)

## ----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")

## ----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")

## ----run-study, eval = FALSE--------------------------------------------------
# runStudy("path/to/study_config.yaml")

## ----session------------------------------------------------------------------
sessionInfo()