--- title: "Introduction to quantMSImageR" author: "Matthew J. Smith" date: "`r Sys.Date()`" package: quantMSImageR output: BiocStyle::html_document: toc_float: true vignette: > %\VignetteIndexEntry{1. Introduction to quantMSImageR} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE) ``` # Introduction A targeted DESI-MRM imaging run produces one intensity per transition per pixel, across tens of thousands of pixels and several acquisitions. Turning that into a biological result runs into several challenges. A tissue section occupies part of the scanned area, and the surrounding slide still returns signal: chemical background, termed background pixels. Left in, those pixels drag down every summary statistic and make images look more uniform than they are. Which pixels are tissue has to be decided per acquisition, and the decision has to survive being re-applied to a different MRM panel of the same sample. Additionally, these pixels aid determining a pseudo noise level for filtering purposes. Targeted panels are limited in number, so often multiple scans of the same area need to be acquired and thus merged. Response depends on ionisation efficiency and local matrix suppression, so raw responses should not be interpreted as analyte amount, and their comparability between analytes, sections or instruments should not be assumed without appropriate normalisation and calibration. ## What this package does `quantMSImageR` addresses each of those: signal-to-noise (SNR) filtering against background pixels, tissue/background separation, per-feature ion images including `.txt` image outputs for overlaying with other imaging (H&E), quantile heatmaps -- and quantification against calibration standards for the third. ## Scope The package is currently validated for Waters DESI-MRM data. The software is extendable to targeted imaging data, which may be analysed after conversion to a compatible `MSImagingExperiment`. However, accessing raw data from other targeted imaging methods is not yet validated. Within that scope it is designed to handle: - background-referenced filtering against non-tissue pixels; - merging of targeted panels; - targeted ion-image visualisation and `.txt` export; - cross-sample summaries such as quantile heatmaps; - calibration-based amount estimation. It does not do the following, by design: | Not provided | Use instead | |---|---| | Untargeted peak detection | a general MSI processing package | | Feature annotation / molecular identification | dedicated annotation tools | | Vendor-neutral raw-data conversion | vendor or community converters | | Histological registration | [TissUUmaps](https://tissuumaps.github.io/); our own helper tools are still in development ([PMID 42438167](https://pubmed.ncbi.nlm.nih.gov/42438167/)) | | General-purpose image segmentation | image-analysis packages | | Complex spatial / colocalisation analysis | [corrMSI](https://github.com/MJS-708/corrMSI/tree/dev) | ## The workflow Functions are grouped into families, and each section below opens with the family it belongs to. The whole workflow is: ``` .raw acquisition ion library CSV | | +--------> readMRM() <----+ | MSImagingExperiment, fData joined to library v selectTissuePixels() --> tissue_pixels.csv (x, y, tissue/background) | labelROIs() --> roi_labels.csv (x, y, region label + id) | +-------------------+-------------------+ | ASSEMBLE alignFeatures -> combineMSIs / bindPanels | | stitchAcquisitions (pieces of one tissue) | | trimMSI | +-------------------+-------------------+ v quant_MSImagingExperiment +-------------------+-------------------+ | FILTER int2response (IS-normalise) | | zero2NA, removeBlankMzs | | int2SNR -> applySNR, back2NA | +-------------------+-------------------+ v +--------------+--------------+ v v VISUALISE QUANTIFY (needs a standards acquisition) imageR summariseCalLevels quantileHm v createCalCurve v int2conc -> "pg_pixel", "pg_mm2" v plotCalCoverage (does calibration cover the data?) WHOLE STUDY: runStudy() (YAML) . generateTxtImages() . runExample() ``` Per-feature images can also be exported as plain `.txt` intensity matrices, which load directly into ImageJ and similar tools. Useful for downstream multimodal analysis such as co-registration with histology or other imaging modalities. Quantification against calibration standards is covered separately in the *Quantification with calibration standards* vignette (`vignette("quantification", package = "quantMSImageR")`). The workflow implemented here is described and applied in [Smith *et al.*, *Analytical Chemistry* (2024)](https://pubs.acs.org/doi/10.1021/acs.analchem.4c02350). # Installation ```{r install, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("quantMSImageR") ``` # Quick start The fastest way to see the whole workflow is `runExample()`, which runs the imaging workflow on two synthetic groups (A, circular tissue, n = 3; B, square tissue, n = 3), renders the HTML report, and also demonstrates quantification on the bundled calibration standards: ```{r quickstart, eval = FALSE} library(quantMSImageR) res <- runExample() # both groups; imaging report + calibration demo res$calibrated # the calibrated object ``` The rest of this vignette reproduces the key steps by hand. # The data model Every acquisition is a `Cardinal::MSImagingExperiment`. `quantMSImageR` adds the `quant_MSImagingExperiment` subclass, which carries two extra slots for calibration and tissue metadata, so most functions coerce with `as()` first. The object is features x pixels: one row per MRM transition, one column per pixel. Calculated quantities -- `response`, `snr`, calibrated amount -- are stored as *additional* spectra layers. Masking functions are the exception: `applySNR()`, `back2NA()` and `zero2NA()` modify the selected `val_slot` in place, so keep a copy of the object when the unmodified values are needed afterwards. | Component | Contains | Accessor | |---|---|---| | spectra slots | `intensity` (raw), `response` (IS-normalised), `snr`, `pg_pixel` / `pg_mm2` (calibrated) | `spectra(x, "name")`, `spectraData(x)` | | feature metadata | `mz`, `name`, transition m/z, analyte type, ion-library columns | `fData(x)` | | pixel metadata | `x`, `y`, `run`, tissue/background label, stage position, region label | `pData(x)` | | calibration | fitted models, calibration levels, fit diagnostics | `calibrationModels(x)`, `calibrationLevels(x)`, `calibrationDiagnostics(x)` | | tissue summaries | pixel and ROI matrices | `tissueMatrix(x)`, `tissueData(x)` | Multiple acquisitions are held in one object, distinguished by the `run` column of `pData()` -- that is what `combineMSIs()` produces and what every per-sample summary groups on. Pixel metadata uses consistent vocabulary: `run` identifies the acquisition, `x`/`y` are grid coordinates, and the tissue mask lives in `sample_name` as `tissue_pixels` / `background_pixels`. The calibration path uses `sample_type` with `Cal` / `Tissue` / `Background`. ```{r data-model} library(quantMSImageR) p <- system.file("extdata", "example.raw", "section01.RDS", package = "quantMSImageR") obj <- as(readRDS(p), "quant_MSImagingExperiment") dim(obj) # features x pixels names(spectraData(obj)) # available layers ``` ```{r data-model-fdata} DT::datatable(as.data.frame(fData(obj)), rownames = FALSE, options = list(pageLength = 5, scrollX = TRUE)) ``` **Table 1.** Feature metadata: one row per MRM transition. ```{r data-model-pdata} DT::datatable(head(as.data.frame(pData(obj)), 100), rownames = FALSE, options = list(pageLength = 5, scrollX = TRUE)) ``` **Table 2.** Pixel metadata, first 100 pixels of one section. # Reading your own data The bundled example ships as `.RDS` so the vignette is reproducible, but real acquisitions enter through `readMRM()`, which reads a Waters DESI-MRM `.raw` folder and joins the ion library to the features: ```{r read-mrm, eval = FALSE} obj <- readMRM( name = "my_acquisition", # .raw folder name, without the suffix folder = "path/to/raw", # directory containing my_acquisition.raw lib_ion_path = "path/to/ion_library.csv" ) ``` The `.raw` folder must contain the exported per-transition text files (an `imaging/` subdirectory) or a previously cached `MSImagingExperiment.rds`. The ion library is a CSV with one row per transition; the columns used for matching are the precursor and product m/z, plus `transition_id` for the display name and `Type` to mark internal standards. Any additional columns -- pathway class, polarity, collision energy -- are carried through to `fData()` and become available for grouping, as shown in *Joining ion-library metadata* below. The result is a features x pixels object exactly like the one inspected above, ready for `selectTissuePixels()` and the rest of the workflow. # Loading the acquisitions *Families: `acquisition`, `combining acquisitions`.* | Function | Takes | Produces | |---|---|---| | `readMRM()` | `.raw` folder + ion-library CSV | `MSImagingExperiment`, `fData` joined to the library | | `trimMSI()` | object with a tissue mask | object with pure-background border rows/columns dropped | | `removeBlankMzs()` | object | object with features that have no data removed | | `alignFeatures()` | two objects | both reduced to common features, m/z-keyed | | `combineMSIs()` | two or more aligned objects | one object, each acquisition its own `run` | | `bindPanels()` | two objects on a shared `(x, y)` grid | features of both on the pixel-grid intersection | | `stitchAcquisitions()` | pieces of one tissue acquired separately | the pixels of every piece, placed where the stage recorded them, as one run | Note the difference between the last two. `combineMSIs()` stacks different sections side by side, adding pixels. `bindPanels()` merges two acquisitions of the same physical area -- typically a positive- and a negative-mode run of one section, or two MRM panels -- adding features to the same pixels: ```{r bind-panels, eval = FALSE} pos <- readMRM("acq_section1_pos", folder = "path/to/raw", lib_ion_path = "path/to/ion_library.csv") neg <- readMRM("acq_section1_neg", folder = "path/to/raw", lib_ion_path = "path/to/ion_library.csv") # Pixels are matched on (x, y), so the two runs need not have identical grids; # pixels present in only one are dropped. both <- bindPanels(pos, neg, label = "S1") ``` The bundled example is six sections: three of group A (circular tissue) and three of group B (square tissue). They ship as RDS, so `readMRM()` is not needed here. `combineMSIs()` binds them into a single object, keeping each section as its own `run`: ```{r load} library(quantMSImageR) library(ggplot2) sections <- sprintf("section%02d", 1:6) objs <- lapply(sections, function(s) readRDS(system.file("extdata", "example.raw", paste0(s, ".RDS"), package = "quantMSImageR"))) combined <- do.call(combineMSIs, objs) groups <- c("A", "A", "A", "B", "B", "B") # sections 01-03 = A, 04-06 = B combined ``` Each section contributed its pixels as a separate `run`, on the same feature axis: ```{r load-check} dim(combined) # features x (pixels from all six sections) table(pData(combined)$run) ``` Everything below operates on this one `combined` object. ## Pieces of one tissue A section too large for one acquisition is sometimes imaged in pieces -- a top and a bottom half, for instance. `stitchAcquisitions()` joins them. `readMRM()` keeps each pixel's stage position (`x_stage`, `y_stage`, in mm), so each piece is placed where it was acquired -- above, below or beside the others -- and the order they are passed does not matter. Pixels between the pieces are simply absent. Pieces that overlap stop with an error rather than averaging the shared pixels. The pieces must share a pixel size. Their grids rarely line up exactly, because each acquisition starts wherever the stage was set, so every pixel is moved to the nearest position on a common grid and the largest shift is reported. `bindPanels()` and `stitchAcquisitions()` do different jobs. Panels over the *same* area -- pos01 and pos02, or a negative- and a positive-mode run -- are merged pixel by pixel with `bindPanels()`. Pieces covering *different* areas are joined with `stitchAcquisitions()`. A sample can need both; the `pieces:` entry of the study YAML does them in the right order. # Selecting tissue pixels *Family: `acquisition`.* | Function | Takes | Produces | |---|---|---| | `selectTissuePixels()` | acquisition name, data path, ion library | `tissue_pixels.csv` (`x`, `y`, `tissue_pixels`, `background_pixels`) plus a preview of the saved mask | Before SNR filtering, each acquisition needs a tissue mask. `selectTissuePixels()` displays every feature's ion image so you can pick the one that best separates tissue from background, then opens an interactive ROI selector: ```{r select, eval = FALSE} selectTissuePixels( name = "my_acquisition", # .raw folder name (no suffix) data_path = "path/to/raw", lib_ion_path = "path/to/ion_library.csv" ) ``` It writes the resulting mask as `tissue_pixels.csv` inside the acquisition's `.raw` folder, where `int2SNR()` and `generateTxtImages()` read it. The mask is simply a per-pixel label -- tissue versus background -- which the bundled sections already carry in their `sample_name` column: ```{r mask-table} mask_tbl <- as.data.frame(pData(combined))[, c("run", "x", "y", "sample_name")] DT::datatable(head(mask_tbl, 500), rownames = FALSE, options = list(pageLength = 5, scrollX = TRUE)) ``` **Table 3.** Per-pixel tissue/background labelling used by `int2SNR()` (first 500 pixels). Plotting that column over the pixel coordinates shows the mask itself, which is the quickest way to confirm a selection before running the workflow -- the circular group A and square group B tissue should stand out cleanly from the surrounding background: ```{r mask-plot, fig.wide = TRUE, fig.cap = "Tissue masks for all six sections. Green pixels are labelled tissue_pixels and contribute to every downstream summary; grey pixels are background_pixels and are used to estimate the noise level in int2SNR()."} ggplot(mask_tbl, aes(x, -y, fill = sample_name)) + geom_tile() + facet_wrap(~ run, ncol = 3) + coord_equal() + scale_fill_manual(values = c(tissue_pixels = "#1b7837", background_pixels = "grey85")) + labs(x = NULL, y = NULL, fill = NULL) + theme_minimal() + theme(axis.text = element_blank(), panel.grid = element_blank(), strip.text = element_text(size = 7), legend.text = element_text(size = 7)) ``` # Labelling regions of interest *Family: `acquisition`.* | Function | Takes | Produces | |---|---|---| | `labelROIs()` | acquisition name, data path, ion library | `roi_labels.csv` (`x`, `y`, `roi_label`, `roi_id`) plus a map of the labelled regions | Once an acquisition has a tissue mask, regions within the tissue can be labelled: airways, parenchyma, a vessel. This is worth doing whenever a whole section is the wrong unit to compare: - **A transition marks a tissue type.** If one of the measured ions is characteristic of a structure -- an airway, a lesion, a particular cell layer -- the ion image itself tells you where that structure is, and the region can be drawn from it rather than from a stain on a neighbouring section. - **Dose defines the region.** In a toxicity study the toxicant is often one of the measured ions, so regions can be drawn where it is high and where it is low in the *same* tissue, and the endogenous species compared between them. That comparison is within one animal, so it is not carried by the differences between animals that a group comparison has to survive. `labelROIs()` shows the ion image of a selected transition, across tissue pixels only. Then draw a region, name it, and choose whether to add another: ```{r label-rois, eval = FALSE} labelROIs( name = "my_acquisition", data_path = "path/to/raw", lib_ion_path = "path/to/ion_library.csv" ) #> Draw a region on the image, then close the window. #> Label for this region (used in this study: airway, vessel): airway #> Saved as airway_03. #> Add another region? [y/n]: y ``` The same label can be used for several regions -- a section usually holds many airways. Each region is numbered as it is drawn (`airway_01`, `airway_02`, ...), so every region keeps its own identity while sharing its label with the others. Labels already used in the study are listed each time, so a region type carries the same name in every sample; the report compares regions by label. A pixel belongs to at most one region, and a region drawn over an earlier one takes those pixels. Tissue outside every region is `unassigned`; background has no region. The regions are saved as `roi_labels.csv` in the acquisition's `.raw` folder -- one file per acquisition, beside the tissue mask, which is left untouched. Running `labelROIs()` again adds to the existing regions; `overwrite = TRUE` starts again and keeps the old file as a backup. For a sample with several panels over the same area, label one of them and pass the others to `copy_to`; the regions are copied by pixel position. ## What the labels look like downstream Drawing regions needs a device and a person, so the bundled sections come with regions already: two called `a` near the middle of each tissue and two called `b` out by the rim, exactly as `labelROIs()` would have written them. Every labelled pixel carries the region's name and its numbered identity: ```{r roi-table} roi_tbl <- as.data.frame(pData(combined))[, c("run", "x", "y", "roi_label", "roi_id")] labelled <- roi_tbl[!is.na(roi_tbl$roi_id) & roi_tbl$roi_id != "unassigned", ] table(labelled$run, labelled$roi_id) ``` Two regions share the label `a` and two share `b`, which is the usual case: a section holds many airways, and each is its own region while the label is what the report compares. Tissue outside every region is `unassigned`, and background has no region at all (`NA`). ```{r roi-map, fig.wide = TRUE, fig.cap = "Regions of interest in all six sections. The two `a` regions sit near the middle of each tissue and the two `b` regions out by the rim; grey is tissue belonging to no region."} ggplot(subset(roi_tbl, !is.na(roi_label)), aes(x, -y, fill = roi_label)) + geom_tile() + facet_wrap(~ run, ncol = 3) + coord_equal() + scale_fill_manual(values = c(a = "#66C2A5", b = "#FC8D62", unassigned = "grey85")) + labs(x = NULL, y = NULL, fill = NULL) + theme_minimal() + theme(axis.text = element_blank(), panel.grid = element_blank(), strip.text = element_text(size = 7)) ``` The comparison the report draws is then one point per sample and region -- the median of that sample's pixels in regions of that label -- rather than one point per pixel: ```{r roi-summary, fig.cap = "One feature by region, one point per section. Each section contributes one value per region label, the median over its pixels there."} feature <- "12-HHTrE" # spectra() is a plain matrix with no dimnames, so a feature is taken by its # row number rather than by name. fi <- match(feature, fData(combined)$name) vals <- data.frame( run = as.character(pData(combined)$run), roi_label = as.character(pData(combined)$roi_label), intensity = as.numeric(spectra(combined, "intensity")[fi, ]) ) vals <- subset(vals, roi_label %in% c("a", "b")) per_sample <- aggregate(intensity ~ run + roi_label, data = vals, FUN = median) per_sample$group <- ifelse(per_sample$run %in% sections[1:3], "A", "B") ggplot(per_sample, aes(roi_label, intensity, fill = roi_label)) + geom_boxplot(outlier.shape = NA, alpha = 0.6, width = 0.6) + geom_point(shape = 21, size = 2) + facet_wrap(~ group) + scale_fill_manual(values = c(a = "#66C2A5", b = "#FC8D62"), guide = "none") + labs(x = NULL, y = paste(feature, "(median per section)")) + theme_classic() ``` Six points per panel, one per section, not the thousands of pixels behind them. The `a` regions are the brighter ones here because they sit where the signal is: nothing about the regions was made to differ, they are simply in different places. # Ion images *Family: `visualisation`.* | Function | Takes | Produces | |---|---|---| | `imageR()` | object, `feat_ind`, `val_slot`, `sample_lab`, `scale` | a `ggplot` ion image | | `quantileHm()` | object, `quant_val`, sample order/labels, `feature_split` | a `ComplexHeatmap::Heatmap` (rows = samples, columns = features) | | `contributionHm()` | the same arguments | the same `ComplexHeatmap::Heatmap`, with group mean as hue and per-sample contribution as opacity | | `quantPalettes()` | palette name, optional `n` | the hex colours used by both, for reuse in your own plots | Colours are consistent across the package and selectable per plot. `imageR()` takes `palette = "heatmap0"` (the default) or `"viridis"`; `quantileHm()` takes a diverging `palette` plus `group_palette` / `feature_palette` for the metadata bars, and `cell_border` for the hairline that keeps neighbouring cells of the same colour from merging into one block. In a YAML-driven study these are set once in the `colours:` block: ```yaml colours: ion_image: "heatmap0" # heatmap0 | viridis heatmap: "heatmap2" # heatmap2 | heatmap0 group: "hat" # hat | reading | heatmap0 feature: "reading" # reading | hat | heatmap0 region: "Set 2" # regions of interest; also any grDevices qualitative palette cell_border: "white" # any colour, or "none" ``` ```{r palettes} names(quantPalettes()) quantPalettes("heatmap0") ``` `region:` also accepts a `grDevices` qualitative palette, which is where its default comes from -- the region colours should not be mistaken for the group ones, and the package's own qualitative sets are already spoken for. To see what any of them look like before choosing: ```{r palette-preview, eval = FALSE} scales::show_col(quantPalettes("hat")) # a package palette grDevices::hcl.pals("qualitative") # the rest, by name scales::show_col(grDevices::hcl.colors(6, "Dark 3")) ``` The [colorspace palette reference](https://colorspace.r-forge.r-project.org/articles/hcl_palettes.html) shows every `grDevices` palette side by side, which is the quickest way to pick one. `imageR()` returns a `ggplot` ion image for a chosen feature. Because it is a plain `ggplot`, the layout can be adjusted -- here into a 3 x 2 grid: ```{r image, fig.wide = TRUE, fig.cap = "Ion image of the first feature across all six sections, before SNR filtering. Sections 01-03 are group A (circular tissue), 04-06 are group B (square)."} imageR(combined, feat_ind = 1, sample_lab = "run", scale = "suppress", palette = "heatmap0") + facet_wrap(~ sample, ncol = 3) + theme(strip.text = element_text(size = 7), legend.text = element_text(size = 7)) ``` ## Colour-scale treatment A single ion image can look very different depending on how the colour scale is built, and the right choice depends on what is being judged. `scale` offers three treatments: - `"suppress"` caps the scale at `percentile` (default 99), so a few hot pixels cannot flatten everything else. This is the default and the right starting point for comparing sections. - `"histogram"` equalises the distribution, which reveals structure in a feature whose signal spans a narrow range -- at the cost of a colour scale that no longer maps linearly to response. - `"sqrt"` compresses the top of the range more gently than histogram equalisation, keeping the scale monotonic in response. ```{r image-scales, fig.wide = TRUE, fig.height = 3, fig.cap = "The same feature and section under the three colour-scale treatments. Only `suppress` keeps the colour scale linear in response; the others trade that for visible structure."} .one <- combined[, pData(combined)$run == "section01"] .p <- lapply(c("suppress", "histogram", "sqrt"), function(sc) imageR(.one, feat_ind = 1, sample_lab = "run", scale = sc) + labs(title = sc) + theme(legend.position = "none", plot.title = element_text(size = 9, face = "bold"), strip.text = element_blank())) patchwork::wrap_plots(.p, nrow = 1) ``` `perc_scale = TRUE` additionally relabels the scale as a percentage of each image's maximum, which is useful when the absolute numbers are not comparable anyway. The HTML report produced by `runStudy()` shows a related set of three views per feature, built with Cardinal's own `image()` rather than `imageR()`: capped at the 95th quantile, histogram-normalised, and unscaled. The same per-feature images can be written to disk as `.txt` intensity matrices (`generateTxtImages(output_txt = TRUE)`, or `output_txt: true` in the YAML config). Those files open directly in ImageJ, which is the usual route into multimodal workflows. # Joining ion-library metadata `buildFeatureMeta()` joins ion-library annotation to features by rounded `(precursor, product)` m/z, so extra columns -- such as the pathway class `Met-1` -- become available for grouping. We do this before the heatmap because the heatmap uses `Met-1` to split its rows. The ion library is the CSV you supply as `lib_ion_path`, and it is worth looking at before anything else: it defines the panel. `transition_id`, `precursor_mz`, `product_mz`, `collision_eV`, `cone_V`, `Polarity` and `Type` are required; anything else is yours to add. Here `Met-1` carries the pathway class and `IS_norm` names the internal standard each analyte is normalised to (see [int2response()]). ```{r ion-lib} lib <- read.csv(system.file("extdata", "example_ion_library.csv", package = "quantMSImageR"), check.names = FALSE) DT::datatable(lib, rownames = FALSE, options = list(pageLength = 10, scrollX = TRUE)) ``` **Table 4.** The example ion library, as supplied on disk. ```{r meta} fm <- buildFeatureMeta(combined, lib) DT::datatable(fm[, c("name", "precursor_mz", "product_mz", "Type", "Met-1", "IS_norm")], rownames = FALSE, options = list(pageLength = 10, scrollX = TRUE)) ``` **Table 5.** The same metadata after `buildFeatureMeta()` has joined it to each measured feature by m/z. # Quantile heatmaps across samples `quantileHm()` summarises each feature to a per-sample quantile and returns a `ComplexHeatmap::Heatmap`. Rows are samples, grouped by study group; columns are features, split by the `Met-1` class joined above -- studies usually have more samples than features, so the long dimension runs vertically. Cells are drawn square (stretching to at most 1.5:1 when one dimension is much longer), so the heatmap grows with the data rather than distorting to fill the device. `cell_size` sets that size in millimetres and `fontsize` the label size, which is what to reach for when long transition names crowd the panel: enlarge the cells and shrink the text rather than letting the labels dictate the layout. The `quant_val` argument chooses which quantile summarises each feature -- a median (`0.5`) reflects the bulk of the tissue, while a high quantile (`0.95`) emphasises hot-spots: ```{r heatmap-median, fig.wide = TRUE, fig.cap = "Median (50th percentile) feature intensity per section. Rows are sections grouped by study group; columns are features split by Met-1 class."} feature_split <- fm[["Met-1"]] quantileHm(combined, quant_val = 0.5, heatmap_order = sections, heatmap_labs = groups, feature_split = feature_split, feature_split_name = "Met-1", cell_size = 10, fontsize = 7) ``` ```{r heatmap-95, fig.wide = TRUE, fig.cap = "The same summary at the 95th percentile, which emphasises hot-spots rather than the bulk of the tissue."} quantileHm(combined, quant_val = 0.95, heatmap_order = sections, heatmap_labs = groups, feature_split = feature_split, feature_split_name = "Met-1", cell_size = 10, fontsize = 7) ``` ## Seeing which sample carried a group effect `contributionHm()` draws the same matrix but encodes two things per cell: the **hue** is the sample's *group* mean for that feature, and the **opacity** is that *sample's* contribution to it. A uniformly solid block is a real group effect; one opaque tile among faded ones is a group mean resting on a single replicate, which a plain group-mean heatmap cannot show. It takes the same arguments as `quantileHm()` and returns the same kind of object, laid out identically -- same group and feature colour bars, same splits, same square cells. Only the cell fill differs, so the two can be read against each other directly. The HTML report draws one or the other, never both: `parameters: heatmap_style:` selects `per_sample`, `contribution`, or `auto`, which picks the contribution view once any group has three or more replicates -- the point at which a group mean can rest on one of them without saying so -- and the per-sample view otherwise. Two details matter for reading it. Values are z-scored per feature across *all* samples rather than within group, so a cell means "high or low for this feature". And the two channels get separate colour limits: a group mean of \(n\) replicates is roughly \(\sqrt{n}\) smaller than a single sample's z-score, so reusing one limit for both washes the panel out. ```{r heatmap-contribution, fig.wide = TRUE, fig.cap = "The same median summary as a group mean, with each sample's contribution to that mean drawn as opacity. A solid block is a group effect; a lone opaque tile is a mean resting on one section."} contributionHm(combined, quant_val = 0.5, heatmap_order = sections, heatmap_labs = groups, feature_split = feature_split, feature_split_name = "Met-1", cell_size = 10, fontsize = 7) ``` # Signal-to-noise filtering *Family: `filtering`.* | Function | Takes | Produces | |---|---|---| | `int2response()` | object + internal-standard name | intensities normalised to the IS (`response`) | | `int2SNR()` | object, `background`/`tissue` labels, `snr_thresh`, optional `snr_overrides` | adds an `snr` slot; below threshold and non-tissue become `NA` | | `applySNR()` | object carrying an `snr` slot | sub-threshold pixels set to `NA` in `val_slot` | | `back2NA()` | object + background label | background pixels set to `NA` | | `zero2NA()` | object, `val_slot` | zeros replaced with `NA` | `int2SNR()` computes a per-pixel signal-to-noise ratio for each feature against the background pixels, and `applySNR()` then sets every sub-threshold pixel to `NA` in the intensity slot. ```{r snr-apply} combined <- int2SNR(combined, val_slot = "intensity", snr_thresh = 3) combined <- applySNR(combined, val_slot = "intensity") ``` `int2SNR()` added an `snr` layer rather than replacing anything, and `applySNR()` then set the sub-threshold pixels of `intensity` to `NA`: ```{r snr-check} names(spectraData(combined)) # snr layer added mean(is.na(spectra(combined, "intensity"))) |> round(3) # fraction now masked ``` Compare the filtered image below with the unfiltered ion images in *Ion images* above: the low-signal background has been removed, leaving the tissue. ```{r snr-after, fig.wide = TRUE, fig.cap = "The same feature after applySNR(). Sub-threshold pixels are now NA, so the low-signal background has been removed and only tissue signal remains (compare to previous ion images)."} imageR(combined, feat_ind = 1, sample_lab = "run", val_slot = "intensity", scale = "suppress") + facet_wrap(~ sample, ncol = 3) + theme(strip.text = element_text(size = 7), legend.text = element_text(size = 7)) ``` # Running a full study from a YAML config *Family: `workflow`.* | Function | Takes | Produces | |---|---|---| | `runStudy()` | a YAML configuration file | runs the whole study: `.txt` images, optional calibration, one HTML report per SNR threshold | | `validateConfig()` | the same file | a report of every problem found in it, before anything is read or written | | `generateTxtImages()` | acquisition names, paths, ion library, SNR settings | list of combined/SNR-filtered objects; optional per-feature `.txt` images | | `runExample()` | nothing (all paths resolved internally) | the same list plus `calibrated` and `report`, having rendered the HTML report | Everything above runs inside `generateTxtImages()`, which `runStudy()` drives from a YAML file. In practice a whole study -- multiple acquisitions, SNR sweeps, ion ratios and optional calibration -- is configured in one place. A minimal config looks like: ```yaml study: "my_study" paths: data_path: "path/to/raw" out_path: "path/to/results" image_dir: "path/to/images" lib_ion_path: "path/to/ion_library.csv" samples: - pos: "acq_ctrl" run_id: "S1" group: "Ctrl" - pos: "acq_treated" run_id: "S2" group: "Treatment" parameters: snr_thresh: [3] output: render_report: true output_txt: true # write per-feature .txt images for ImageJ ``` The report's ion-image section opens with an animation of every feature -- one frame per metabolite, all sections at once, ordered by pathway, each frame labelled with a box in its pathway's colour and the intensity its colour scale is capped at. The same file is written beside the report as a `.gif`, which is what goes into a talk. Switch it off or narrow it with `output: ion_gif:`: ```yaml output: ion_gif: enabled: True features: "all" # or ["PGE2", "AA"] seconds: 0.4 ``` Every image in the report -- the animation, the dropdown's first rendering and the concentration images -- caps its colour scale at the same percentile, and each caption states the value it used: ```yaml output: ion_images: percentile: 95 ``` This is separate from `parameters: perc:`, which caps the `.txt` images written for other software: those are an input to further processing and want their dynamic range, while the report's images want a map that can be read. The shipped template documents every option (SNR overrides, ratios, calibration, ...) -- open it to use as a starting point: ```{r yaml-template, eval = FALSE} file.show(system.file("config_template.yaml", package = "quantMSImageR")) ``` Then point `runStudy()` at your filled-in config: ```{r run-study, eval = FALSE} runStudy("path/to/study_config.yaml") ``` The same runner is available from the command line, which is the usual way to run a large study unattended: ```sh Rscript -e 'quantMSImageR::runStudy("study_config.yaml")' ``` ## Samples made of several acquisitions A sample is often more than one `.raw` folder, and the entry says what the files are: ```yaml samples: # One acquisition. - neg: "M120-1" run_id: "Male-120" group: "Young" # Several panels over the SAME area -- neg + pos, or pos01 + pos02. # Merged pixel by pixel (bindPanels). - neg: "M316-PLs" pos: ["M316-pos01", "M316-pos02"] run_id: "Female-316" group: "Young" # One tissue acquired in PIECES, e.g. top and bottom. Panels are merged # within each piece, then the pieces are placed where the stage recorded # them and joined into one sample (stitchAcquisitions). - run_id: "Male-707" group: "Old" pieces: - neg: "M707-top" pos: ["M707-top_pos01", "M707-top_pos02"] - neg: "M707-btm" pos: ["M707_btm-pos01", "M707_btm-pos02"] ``` Two halves of one tissue are one sample, not two. Left as separate samples they count twice in every group summary, and because the halves agree with each other they look like a reproducible effect rather than one measurement. ## Comparing regions of interest When any acquisition has regions, every analysis section of the report is drawn at sample level first and then at region level, so the two are read the same way rather than the regions living in a section of their own: | Section | By sample | By region | |---|---|---| | 2. Heatmap | one row per sample | one row per sample and region, split by region, **each region z-scored against itself** | | 3. Comparison plots | box + violin over pixels | box plots with one point per sample or per region | | 4. Correlation | pixels of each sample | pixels of each region type | | 5. Regions of interest | -- | the region maps and how many pixels each holds | Scoring each region against itself is what makes the heatmap answer "which samples are high for this feature *in their airways*" rather than restating that airways and parenchyma differ. `heatmap: False` leaves that view out. The box plots are arranged by `compare`: - `"between_groups"` -- each region type on its own, groups side by side: airway, Ctrl v HDM; vessel, Ctrl v HDM; - `"between_rois"` -- each group on its own, region types side by side: Ctrl, airway v vessel; HDM, airway v vessel; - `"both"` -- both sets. An `ROI_summary` sheet in the tables workbook carries one row per point (see `unit`) and feature, with pixel quantiles and the pixel count. Correlation is computed at **one** level, never both, since it is the slowest part of the report: `output: colocalisation` takes `'roi'`, `'sample'`, `True` (the same as `'sample'`), `False`, or `auto`, which correlates within the regions when the study has them and is off otherwise. A region holds a fraction of a section's pixels, so the region level is the quicker one as well as the more specific question. What each point stands for is set by `unit`: - `"sample"` -- one point per sample and region type: the median over all of that sample's airways, for instance; - `"roi"` -- one point per region: `airway_01`, `airway_02`, ... each on its own. One or the other, not both: a comparison has a single unit, and offering two views of the same numbers invites reading whichever looks better. ```yaml roi: enabled: auto # auto (on when any acquisition has roi_labels.csv) | True | False include_unassigned: False # tissue outside every region as its own "unassigned" type | leave out compare: "between_groups" # between_groups | between_rois | both heatmap: True # the by-region heatmap view in section 2 | False unit: "sample" # rows/points in the region views: sample | roi output: colocalisation: auto # auto = within regions when there are any, else off | 'roi' | 'sample' | False ``` The points are samples or independent regions, not pixels. The pixels of one region are repeated measurements of the same tissue; a box plot over pixels would show thousands of points where there are only a handful of regions or animals. ## A second factor: blocking Some designs carry a factor that is not the comparison: the donor a paired set of treatments came from, the sex of an animal, the slide an acquisition sat on. `block:` names it per sample, beside the `group:` that is being compared. ```yaml samples: - neg: "acq_d1_ctrl" run_id: "D1_Ctrl" block: "D1" # the donor these two share group: "Ctrl" - neg: "acq_d1_treated" run_id: "D1_Treated" block: "D1" group: "Treated" - neg: "acq_d2_ctrl" run_id: "D2_Ctrl" block: "D2" group: "Ctrl" - neg: "acq_d2_treated" run_id: "D2_Treated" block: "D2" group: "Treated" parameters: block_name: "Donor" # what block: means, for legends and annotations block_centre: "center" # off | center | baseline ``` Declaring a block annotates and groups the report by it. Correcting for it is separate, and off by default, so adding `block:` on its own changes no number. ### What correcting does A block carries a level of its own: one donor's sections ionise better, one slide was run on a different day. That level sits on top of every comparison made inside the block. `block_centre` divides each feature by a reference computed within the block, which removes it: - `center` divides by the block's own geometric mean. It needs no privileged sample and pools every sample in the block, so it is the steadier choice. - `baseline` divides by the geometric mean of the block's `baseline_group` samples, so a value reads directly as a fold change against its own control. Clearer, and noisier: the reference's noise propagates into every other sample of that block. A block with no `baseline_group` sample falls back to its own mean and says so. The correction is multiplicative on intensity, which is additive in log space. A block-level gain is removed; differences in how far a treatment moves a feature are not. A donor who responds twice as strongly as another still reads twice as strongly afterwards, which is usually the point of the comparison. It reaches the heatmaps and the comparison plots, where the block's level gets in the way. The ion images are left uncorrected: an image should show the intensity that was measured, and its colour-scale cap should be a real number. A corrected report says so at the top and gains `blocked` in its filename, so it cannot overwrite an uncorrected run of the same study. `blockCentre()` applies the same correction outside the report. ### Ion images scaled within a block `scale_within` decides which image panels share a colour scale: ```yaml output: ion_images: percentile: 95 # where every image in the report caps its scale scale_within: "block" # all | block | group ``` `all` is one scale for the study, which is readable when the sections are of similar intensity and not when they are not: the weakest reads as empty and the strongest saturates. `block` gives each block its own, so each uses its full range; the grid is then one column per block, its groups as rows and replicates across, and every column states the value it is capped at. That number is what keeps the columns comparable once the colours no longer are. With a block the ion-image section draws four renderings per feature -- scaled within each block, then the three study-wide ones (capped at the percentile, histogram normalised, unscaled). Without one it draws the usual three. ### Pieces that overlap `stitchAcquisitions()` refuses two pieces that claim the same grid position, because averaging them or dropping one silently would invent a pixel. Passes that deliberately share a strip are a real case, so the behaviour is settable: ```yaml parameters: stitch_overlap: "first" # error (default) | first | last ``` `first` keeps the earlier piece's pixels where they collide, `last` the later one's, and either reports how many were dropped from each piece. ## The configuration as a record of the analysis Driving a study from a file rather than from a session means a complete human-readable configuration file to document how the data were processed. Every choice that affects the results (which acquisitions belong to which group, the SNR thresholds and any per-analyte overrides, whether internal-standard normalisation ran and against which standard, the calibration design and its weighting, how much extrapolation was tolerated) is stored for future reference. The file is small enough to deposit alongside the data, cite in a methods section, or track in version control in line with FAIR principles (Wilkinson *et al.*, 2016). To ensure accuracy of this documentation `validateConfig()` runs every check rather than stopping at the first failure, so one call reports every problem in the file: ```{r validate, eval = FALSE} validateConfig("path/to/study_config.yaml") ``` ``` quantMSImageR validation: 2 error(s), 1 warning(s), 26 check(s) Errors: - parameters$average_method = 'geometric'; use 'mean' or 'median'. - parameters$is_name = 'PGE2-d4' is not a value of the ion library's 'Type' column. Values present: Analyte, IS. Note this is a type label, not a transition name. ``` `runStudy()` calls it first and stops if there are errors, so a configuration that cannot be trusted never produces an output. The checks are deliberately cross-file for every enabled functionality (the ion library, the calibration metadata and the config have to agree with each other, not merely be individually well-formed). # 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). Wilkinson MD, *et al.* The FAIR Guiding Principles for scientific data management and stewardship. *Scientific Data* 3, 160018 (2016). # Session information ```{r session} sessionInfo() ```