The curated workbook

MRManalyzeR expects one .xlsx workbook holding the instrument export plus the two curation sheets that describe it:

Keeping the curation next to the data is what makes a run reproducible: the decisions about which compounds are reportable and which injections belong to the study are recorded in the workbook, not in a script.

Two practical rules govern whether a run will work at all:

• Names must match exactly. feature_metadata$Processing_name has to match the compound names in the acquisition method, and sample_metadata$Name has to match the injection names in the data sheet – including trailing spaces, which are the usual culprit.
• Exclusions are recorded, not deleted. An injection you do not want analysed gets Include = "NO" rather than being removed; a compound that should not be reported gets Report = "NO". Re-injections are the common case, and marking them keeps the workbook a faithful record of the run.

xlsx  <- system.file("extdata", "example_data.xlsx",
                     package = "MRManalyzeR")
fdata <- openxlsx::read.xlsx(xlsx, sheet = "feature_metadata")
meta  <- openxlsx::read.xlsx(xlsx, sheet = "sample_metadata")

c(injections = nrow(meta), compounds = nrow(fdata))

#> injections  compounds 
#>         79        164

The bundled workbook is a subset of the oxylipin LC-MS/MS data of Kolmert et al. (2018) (see References): 79 injections of mouse bronchoalveolar lavage over a calibration ladder, extraction blanks, pooled QC injections and 45 study samples.

DT::datatable(as.data.frame(table(Sample_type = meta$Sample_type)),
   rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

Because the panel is known in advance, feature_metadata can carry real biology alongside the reporting decision, for instance the enzymatic pathway a compound comes from, its precursor PUFA, its retention time. Those columns serve a functional role: they colour the PCA loadings, group the heatmap rows and define the metabolite ratios, and they travel with the data as variable_meta all the way to the output.

DT::datatable(fdata[fdata$Report == "YES",
         c("Compound", "Enzymatic_pathway", "PUFA", "RT", "Report")],
   rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

The TargetLynx export

TargetLynx does not export a matrix. Its “complete summary” (Edit > Copy > Complete Summary, pasted straight into the lcms_data sheet) is a stacked table – one block per compound, each block its own little table of injections, with the compound name in a header row between blocks. It is tedious to inspect visually and does not parse into R with existing tools, which is the first problem the package solves.

readTargetLynx() parses those blocks, pivots them into a wide sample x compound matrix of the requested datatype column, and (with snr) masks any peak whose signal-to-noise falls below the threshold:

m <- readTargetLynx(xlsx, datatype = "Area", snr = 3)
dim(m)

#> [1]  82 155

DT::datatable(round(m[, 1:6], 1), filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

Columns come back as the raw TargetLynx compound identifiers (Processing_name). feature_metadata maps those to the canonical Compound name used everywhere downstream, which is how an instrument method can be renamed without touching the analysis.

Three shapes have now appeared – 79 injections and 164 compounds in the metadata, a parsed matrix of 82 x 155 – and they describe different populations at different stages:

• feature_metadata lists the whole curated panel, 164 compounds, including everything acquired but not reported. Eighteen carry Report = "YES".
• sample_metadata lists 79 injections. The export holds 82: three were run but never entered in the metadata, and an injection the metadata does not describe cannot be analysed, so they are dropped on assembly.
• The parsed matrix is 82 x 155 because the export carries 155 compound blocks – internal standards and panel members outside the reported set included.

Assembly applies the reporting decision and the name matching, and nothing else:

de_real <- assembleDataset(m, fdata, meta)
dim(as.data.frame(de_real$data))

#> [1] 79 14

Fourteen compounds of the eighteen reported, and 79 injections rather than the 57 flagged Include = "YES": assembleDataset() keeps every injection the metadata describes, and Include is applied later, by processDataset() or by subsetDataset().

The four missing compounds are worth understanding, because they are a trap in the manual route only. assembleDataset() matches matrix columns against the Compound column, not Processing_name:

rep18 <- fdata[fdata$Report == "YES", c("Processing_name", "Compound")]
rep18[rep18$Processing_name != rep18$Compound, ]

#>           Processing_name             Compound
#> 13 12,13-DiHOME (erythro) erythro-12,13-DiHOME
#> 14   12,13-DiHOME (threo)   threo-12,13-DiHOME
#> 68             13-oxo-ODE              13-KODE
#> 72     LTB4 (12-epi-LTB4)                 LTB4

Exactly those four differ, and none of them appears in this export under either name, so here they are simply absent. But the mechanism is general: a matrix whose columns are raw Processing_name identifiers loses every compound whose two names differ, silently, because the column no longer matches a Compound value. The YAML route never hits this – processDataset() renames Processing_name to Compound before assembling. Calling assembleDataset() yourself means doing that rename yourself, or passing a matrix that already carries canonical names.

readSkyline() returns the same contract from a Skyline molecule x sample sheet, with signal_filter = "LOD" or "LOQ" in place of the S/N mask.

Output from any other software can be supplied as a plain sample-by-analyte matrix, which enters the workflow at assembleDataset() and is processed identically from there on. Readers for further vendor formats are planned for a later version; suitable test data is being sought, since a reader written without real exports to check it against is not one worth relying on.

The bundled dataset

The example above is a single chromatographic batch, injected from one reconstitution volume, with no per-sample amount to normalise to so normalisation, concentration adjustment and batch correction would all be incompatible. The package therefore bundles a second, entirely simulated workbook: a Skyline-layout export of 20 analytes across three enzymatic pathways, 54 injections split over two chromatographic batches, with a per-sample protein amount and a different reconstitution volume per batch.

The study design is balanced across the batches, so the batch effect and the treatment effect are not confounded.

sxlsx <- system.file("extdata", "example_synthetic.xlsx",
                     package = "MRManalyzeR")
sfeat <- openxlsx::read.xlsx(sxlsx, sheet = "feature_metadata")
smeta <- openxlsx::read.xlsx(sxlsx, sheet = "sample_metadata")

table(smeta$Chrom_Batch, smeta$Sample_type)

#>     
#>      Blank QC Sample
#>   B1     2  5     20
#>   B2     2  5     20

Assembling a DatasetExperiment

assembleDataset() aligns the matrix with its feature and sample metadata and returns a struct::DatasetExperiment. Every step from here on takes and returns that object, so the metadata travels with the data and no step has to be handed a matching sample table again.

sm <- readSkyline(sxlsx, "skyline_data", sfeat, signal_filter = "LOD")
de <- assembleDataset(sm, sfeat, smeta)
de

#> A "DatasetExperiment" object
#> ----------------------------
#> name:          
#> description:   
#> data:          54 rows x 20 columns
#> sample_meta:   54 rows x 16 columns
#> variable_meta: 20 rows x 9 columns

This works directly because the simulated panel uses one name for each compound – Processing_name and Compound are identical for all 20 – so the distinction drawn above never arises. On a real workbook, check that before handing a reader’s output straight to assembleDataset().

Blank filter

filterBlanks() finds the extraction blanks via sample_meta and sets any value below mean(blanks) * blank_filter to NA. A factor of 3 is a common choice – a peak has to be at least three times the blank level to be treated as real.

It needs blank values to work from, which makes it worth seeing on data that has them. Two compounds, three samples and two extraction blanks:

bl <- assembleDataset(
  data.frame(well_above = c(100, 120, 110, 10, 12),
             near_blank = c(30, 28, 26, 25, 24),
             row.names  = c("S1", "S2", "S3", "B1", "B2")),
  data.frame(Compound = c("well_above", "near_blank"), Report = "YES"),
  data.frame(Name        = c("S1", "S2", "S3", "B1", "B2"),
             Sample_type = c("Sample", "Sample", "Sample", "Blank", "Blank")))

before <- as.data.frame(bl$data)
bl     <- filterBlanks(bl, blank_filter = 3, blank_head = "Sample_type",
                       blank_name = "Blank")
cbind(before, after = as.data.frame(bl$data))

#>    well_above near_blank after.well_above after.near_blank
#> B1         10         25               NA               NA
#> B2         12         24               NA               NA
#> S1        100         30              100               NA
#> S2        120         28              120               NA
#> S3        110         26              110               NA

well_above has blanks at 10 and 12, so the threshold is 33 and all three samples clear it. near_blank has blanks at 25 and 24, so the threshold is 73.5 and nothing in the compound survives – which is the point of the step. Both are recorded:

as.data.frame(bl$variable_meta)[, c("blank_threshold", "n_masked_blank")]

#>            blank_threshold n_masked_blank
#> well_above            33.0              2
#> near_blank            73.5              5

n_masked_blank counts every cell masked, the blank injections themselves included: they sit below their own threshold by construction.

On the simulated workbook the step has nothing to do, and it is worth seeing why rather than running it and reading the number as an effect:

blank_rows <- rownames(de$data)[
  as.data.frame(de$sample_meta)$Sample_type == "Blank"]
c(missing_total = sum(is.na(de$data)),
  missing_in_blanks = sum(is.na(as.data.frame(de$data)[blank_rows, ])),
  blank_cells = length(blank_rows) * ncol(de$data))

#>     missing_total missing_in_blanks       blank_cells 
#>               154                80                80

Every blank cell is already NA: readSkyline(signal_filter = "LOD") masked the blanks along with everything else below the detection limit, which is what you would hope for in a clean extraction. filterBlanks() therefore has no blank values to average and leaves the matrix alone, and the missing values counted here are the LOD mask’s work, not the blank filter’s.

Detection filters

Detection filters remove poorly measured compounds and injections before imputation: filterFeatures() drops compounds, filterSamples() drops injections that missed most of the panel. Both are off by default.

What matters is which samples the detection fraction is computed over. Take a compound present throughout one treatment group and absent from the others:

demo <- assembleDataset(
  data.frame(everywhere    = c(10, 12, 11, 9, 10, 11),
             group_specific = c(NA, NA, NA, 20, 22, 21),
             sparse         = c(5, NA, NA, NA, NA, 6),
             row.names = paste0("S", 1:6)),
  # assembleDataset() keeps only rows flagged Report == "YES"
  data.frame(Compound = c("everywhere", "group_specific", "sparse"),
             Report   = "YES"),
  data.frame(Name        = paste0("S", 1:6),
             Sample_type = rep("Sample", 6),
             Treatment   = rep(c("PBS", "HDM"), each = 3)))

across <- filterFeatures(demo, min_frac = 0.6, method = "across")
within <- filterFeatures(demo, min_frac = 0.6, method = "within",
                         group_col = "Treatment")

list(across = colnames(across$data), within = colnames(within$data))

#> $across
#> [1] "everywhere"
#> 
#> $within
#> [1] "everywhere"     "group_specific"

"across" pools all six samples, so group_specific scores 3/6 = 0.5, falls short of 0.6 and is deleted because it is group-specific. Within-group filtering retains compounds consistently detected in at least one group: "within" scores each treatment separately, group_specific is 3/3 in HDM and survives, while sparse still goes.

The threshold matters as much as the method here. At min_frac = 0.5 the same compound scores exactly 0.5 under "across" and is kept, since the test is “at least”, so the two methods agree and nothing is demonstrated. A group-specific compound is only ever at risk once the pooled threshold is set above its group’s share of the samples.

A third mode, "QC", scores against the QC injections – whether a compound is detected consistently in the QCs, which is a prerequisite for trusting its precision rather than a measure of it. It falls back to "within" with a warning when a run has no QCs.

Small groups quantise the threshold: with n = 3 a fraction can only be 0, 1/3, 2/3 or 1, so every min_frac between 1/3 and 2/3 behaves identically.

QC and blank injections never count towards "within" or "across". They stay in the dataset – the blank filter and correctBatch() still need them – but a compound seen only in the QCs cannot keep itself alive. The fraction each feature was judged on is kept as detected_frac.

On the working dataset, where readSkyline() has already masked values below the LOD:

de <- filterFeatures(de, min_frac = 0.5, method = "within",
                     group_col = "Treatment")

filterSamples() then drops injections by missing-value fraction. Run it after filterFeatures(): judged against a feature list that still contains everything the instrument barely saw, every injection looks poor because the denominator is padded with noise. Only study samples are removed – a QC over the threshold is reported but kept, since dropping one can leave a batch with nothing for correctBatch() to anchor to.

de <- filterSamples(de, max_na = 0.8)
summary(as.data.frame(de$sample_meta)$na_frac)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.0000  0.0000  0.0500  0.1426  0.1000  1.0000

Both run once over the whole dataset rather than per batch: a detection fraction is a property of the study, and filtering per batch would leave each batch with a different set of columns. In a YAML config they are the filter_features: and filter_samples: blocks, applied after the blank filter and before imputation – imputation fills the very cells these measure.

Batch correction

Where an acquisition spans several days or columns, response can shift between batches. correctBatch() scales each batch so its per-feature median matches the grand median across the reference samples. It runs straight after the filters, on the instrument signal, and before normalisation, concentration adjustment and imputation. A missing value counts as zero: a compound missing from most of a batch’s QCs is described as such, thus impacting the median calculation, avoiding it being based on few injections where it was found. Where more than half of values for a compound are missing, the median is zero, and correction is not applied for this compound only. The failing batches are printed and kept in bc_failed_batches. To correct only part of the panel, pass feature_col – a feature_metadata column of YES / NO.

batch     <- stats::setNames(de$sample_meta$Chrom_Batch,
                             rownames(as.data.frame(de$data)))
de_uncorr <- de

de <- correctBatch(de, qc_label = "QC",
                    factor_name = "Sample_type",
                    batch_head  = "Chrom_Batch")
de_bc     <- de

The correction equalises the medians of the reference injections, so that is where to look at it. Taking the most completely measured analyte, and counting a missing value as zero exactly as the reference-median calculation does:

qc     <- as.data.frame(de_uncorr$sample_meta)$Sample_type == "QC"
j      <- which.min(colSums(is.na(as.data.frame(de_uncorr$data))))
zeroed <- function(v) median(replace(v, is.na(v), 0))

rbind(
  before = tapply(as.data.frame(de_uncorr$data)[qc, j], batch[qc], zeroed),
  after  = tapply(as.data.frame(de_bc$data)[qc, j],     batch[qc], zeroed))

#>             B1      B2
#> before 1662.70  837.20
#> after  1243.85 1243.85

The two QC medians land on the same value, which is the whole of what the step promises. The same table over all injections does not come out equal, and should not be read as a failure:

rbind(
  before = tapply(as.data.frame(de_uncorr$data)[[j]], batch, zeroed),
  after  = tapply(as.data.frame(de_bc$data)[[j]],     batch, zeroed))

#>             B1      B2
#> before 1670.40  837.20
#> after  1249.61 1243.85

Two limitations follow from that, and both matter when reading a corrected dataset:

• It removes a between-batch shift in the reference median, nothing else. Biological differences between the samples in each batch are left alone – which is the point, and also why a design confounded with batch cannot be rescued here.
• It does not correct drift within a batch. A signal falling steadily across a sequence still falls after correction; plotDrift() is where that shows up, and pmp’s spline-based QC-RSC is the tool for it.

Counting a missing value as zero applies only to the reference median – it is what stops a batch whose QCs mostly missed a compound from claiming a high reference. The missing measurements themselves stay missing, and are still missing after correction.

A compound whose reference median is zero in any batch is left as measured in all of them rather than corrected in some. In this dataset the two analytes simulated near their detection limit are below it in more than half of the batch-2 QCs, which puts that median at zero. They are listed here, and they still carry the batch difference in the PCA further down:

vm <- as.data.frame(de_bc$variable_meta)
vm[nzchar(vm$bc_failed_batches), c("Compound", "bc_failed_batches")]

#>              Compound bc_failed_batches
#> Analyte_07 Analyte_07                B2
#> Analyte_16 Analyte_16                B2

Use pooled QC injections as the reference (qc_label = "QC") whenever the design has enough of them in every batch; falling back to "Sample" assumes the batches are biologically comparable.

Normalise

normaliseMatrix() divides each sample by a known per-sample amount – protein, tissue weight, or the volume of biofluid extracted – so that concentrations are expressed per unit of starting material:

de <- normaliseMatrix(de, column = "protein_ug")
round(as.data.frame(de$data)[1:3, 1:3], 3)

#>         Analyte_01 Analyte_02 Analyte_03
#> SYN_001         NA         NA         NA
#> SYN_002      2.344     14.980     45.445
#> SYN_003      6.766      7.592     33.770

Concentration adjustment

Before the arithmetic, a word on units, because they change under your feet through this stage. The simulated workbook reports a simulated response, not a calibrated concentration – the generator plants signals, it does not model a calibration curve – so nothing below should be read as ng/mL. On a real panel the chain runs: the reported datatype (area, response or a concentration from the vendor’s curve), divided by the per-sample amount at normaliseMatrix() (so, per ug protein here), then rescaled by adjustConcentration(), which does not change the unit but does change what the number refers to – the amount in the original sample rather than in the vial. Whatever that chain ends at is what every axis label and every units: entry downstream should say.

A calibration curve returns the concentration in the vial the instrument sampled, and it does so based on the calibrant vial: the curve was built at the calibrant’s internal-standard concentration, and it assumes the sample carries the same. adjustConcentration() corrects that assumption (based on the same volume of calibration and study sample being injected into the instrument):

\[ \text{conc} \times \frac{\text{cal vol}}{\text{sample vol}} \times \frac{\text{sample IS}}{\text{cal IS}} \]

This uses the volumes of internal standard added as a proxy for the amounts of internal standard present, which holds only if the calibrant and the samples were spiked from the same stock at the same concentration. Where they were not, the ratio has to be built from amounts instead.

Enable it only when the correction has not already been applied upstream. Vendor software is often configured to do this at the point of quantification, and running it again here would apply the same factor twice – which is invisible in the output, since the result is a plausible number either way.

Ultimately leading to:

\[ \frac{\text{cal vol}}{\text{sample vol}} \times \frac{\text{sample IS}}{\text{cal IS}} \;=\; \frac{\text{sample IS} / \text{sample vol}}{\text{cal IS} / \text{cal vol}} \;=\; \frac{[\text{IS}] \text{ in the sample vial}} {[\text{IS}] \text{ in the calibrant vial}} \]

So the correction is the ratio of internal-standard concentrations. More vial volume per unit of IS means the IS is more dilute in your sample than in the calibrant; a smaller IS peak inflates the analyte/IS response ratio, the curve reads high, and this factor – below 1 – scales it back down.

Note what is not being corrected here. Reconstituting in more solvent dilutes the analyte and the internal standard equally, so the response ratio does not move and neither does the curve’s answer.

Setting starting_vol_col adds a second step, back to the concentration in the original, undried sample:

\[ \text{conc}_{\text{original}} = \text{conc}_{\text{vial}} \times \frac{\text{sample vol}}{\text{starting vol}} \]

de <- adjustConcentration(de,
                           sample_vol_col   = "sample_volume_uL",
                           cal_vol_col      = "cal_vol_uL",
                           sample_IS_col    = "sample_IS_vol_uL",
                           cal_IS_col       = "cal_IS_vol_uL",
                           starting_vol_col = "starting_vol_uL")
round(as.data.frame(de$data)[1:3, 1:3], 5)

#>         Analyte_01 Analyte_02 Analyte_03
#> SYN_001         NA         NA         NA
#> SYN_002    0.58607    3.74509   11.36127
#> SYN_003    1.69152    1.89799    8.44240

Impute missing values

imputeMissing() fills the remaining NAs in every injection that is not a blank – biological samples and QCs alike – with values below that compound’s lowest measurement, treating a non-detect as lying under the detection floor rather than at zero. Only the blanks are exempt, so frac_imputed below is a fraction of samples and QCs. The ceiling is scalar times the minimum – 0.2 here – and by default every cell gets its own random draw between 0 and that ceiling, from a normal centred mid-range. A single constant would give every non-detect the same value, which shrinks a compound’s variance and ties its ranks; method = "constant" does that if it is wanted. The draws use a fixed seed, so a rerun is identical. Blanks are left untouched.

The PCA pre-processing (mrmRunPCA(impute = )) offers the same random fill ("gaussian", its default) alongside "frac_min", "min", "half_min" and "none", and missing-value filters that drop features or samples above a missingness threshold instead of imputing them. For targeted panels a compound missing in most samples is usually better dropped than modelled.

Imputation comes last, after batch correction, normalisation and concentration adjustment. In a YAML-driven run it fills the exported matrix and the PCAs; the tests and plots use the measured values unless use_imputed: True is set. With replace_MVs off, the PCAs impute for themselves.

de <- imputeMissing(de, scalar = 0.2, blank_head = "Sample_type",
                     blank_name = "Blank")
sum(is.na(de$data))

#> [1] 80

Imputation is recorded in variable_meta: the ceiling (impute_max), how many cells it filled (n_imputed) and what share of the compound’s non-blank injections that is (frac_imputed).

Results for heavily imputed compounds are more sensitive to the imputation method. A high value does not by itself make a difference an artefact – a compound genuinely absent from one treatment carries one too – so read it next to the detection counts per group rather than on its own:

vm <- as.data.frame(de$variable_meta)
head(vm[order(-vm$frac_imputed),
        c("impute_max", "n_imputed", "frac_imputed")])

#>            impute_max n_imputed frac_imputed
#> Analyte_07 2.13364437        19         0.38
#> Analyte_16 0.06665625        19         0.38
#> Analyte_01 0.11721461         2         0.04
#> Analyte_02 0.15344484         2         0.04
#> Analyte_03 1.01602752         2         0.04
#> Analyte_04 0.01385953         2         0.04

The whole matrix in one call

processDataset() composes stages 1 and 2 – read, assemble, filter, batch-correct, normalise, adjust, then impute – driven by the same metadata sheets:

out <- processDataset(
  sfeat, smeta,
  xlsx_path        = sxlsx,
  data_source      = "skyline",
  data_tab_names   = "skyline_data",
  signal_filter    = "LOD",
  blank_filter     = 3,
  normalize        = "protein_ug",
  adjust_conc      = TRUE,
  starting_vol_col = "starting_vol_uL",
  replace_MVs      = 0.2,
  batch_correction = TRUE,
  bc_qc_label      = "QC",
  bc_header        = "Chrom_Batch",
  blank_head       = "Sample_type")

dim(out[[1]]$data)

#> [1] 54 20

out[[2]] holds the features dropped along the way with the reason recorded – the same table the data-quality report surfaces as “compounds excluded”, and the first place to look when a compound you expected is missing from the output.

The shape matches what the manual walkthrough produced, but the two are not the same calculation, and matching dimensions would not show that they were. The call above leaves filter_features and filter_samples at their defaults – off – where the walkthrough applied the detection filter, and processDataset() imputes per processing batch where imputeMissing() was called once over everything. Align the arguments if you want the two routes to agree cell for cell.

out$measured is the same dataset before imputation, gaps included – what the reports and statistics use in a YAML-driven run unless use_imputed: True is set. out[[1]] has its gaps filled everywhere except in the blanks, which are left as they were:

sum(is.na(out$measured$data))                         # gaps as measured

#> [1] 154

sum(as.data.frame(out[[1]]$variable_meta)$n_imputed)  # filled in out[[1]]

#> [1] 74

The rest of this vignette carries on from the run just made, using the measured copy – the one that still has its gaps:

syn <- out$measured
dim(as.data.frame(syn$data))

#> [1] 54 20

sum(is.na(as.data.frame(syn$data)))

#> [1] 154

That matters for everything that follows. CVs, group summaries and tests are meant to be read off measured values; computed on a filled matrix they would report invented numbers as observations, and a compound that was a non-detect in half of one group would show a full complement of samples. The PCA is the exception, and it handles itself: mrmRunPCA() fills internally for the decomposition, because a PCA cannot take a gap.

The same dataset also ships ready-made as system.file("extdata", "example_synthetic.RDS", package = "MRManalyzeR"), for examples that need a processed object without doing the processing. Note that the shipped copy is the filled one – out[[1]], not out$measured – so load it for convenience, not for a demonstration of measured values.

How well was each compound measured?

addCVMetrics() writes the per-compound quality metrics into variable_meta, matched to the feature they describe, so annotation and quality live in the same table. runMRManalyzeR() applies it to every run.

syn <- addCVMetrics(syn, qc_label = "QC")

cv <- as.data.frame(syn$variable_meta)[
  , c("Compound", "Enzymatic_pathway", "CV_QC", "CV_sample",
      "CV_sample_vs_QC")]

DT::datatable(cv[order(cv$CV_QC), ], rownames = FALSE, filter = "top",
              options = list(pageLength = 10, scrollX = TRUE))

• CV_QC is the coefficient of variation across the pooled QC injections. Every QC is the same material, so this is technical variability.

• CV_sample is the same statistic across the biological samples, so it carries the biological spread plus that same technical noise.

• CV_sample_vs_QC is CV_sample / CV_QC, a descriptive comparison of the two dispersions above. A large value says the samples vary much more than the QCs, which is what you want to see in a compound you intend to analyse; a value near 1 says the spread across the study is no wider than the repeatability of the method.

  Two cautions. It is a ratio of coefficients of variation, not of variances, and not a decomposition into biological and technical parts: a CV carries its own mean, so the comparison is only as meaningful as the similarity of the QC and sample means. And it is not the dispersion ratio of Broadhurst et al. (2018), which is defined on standard deviations and in the other direction, \(100 \times SD_{\mathrm{QC}} / SD_{\mathrm{samples}}\), with its commonly quoted thresholds attached to that definition. The two are not reciprocals of each other unless the means coincide, so do not carry a D-ratio threshold across to this column. Compute the published statistic explicitly if that is the one you need.

Analytical drift

plotDrift() puts one compound’s response against acquisition order, coloured by sample type so the pooled QCs can be tracked through the sequence. A trend, a step or a widening spread in the QCs is instrumental, not biological.

plotDrift(syn, "Analyte_02", value_label = "response / ug protein")

Normality

Q-Q plots help assess distributional assumptions within comparison groups. The parametric tests assume approximate normality within each group – or of the model residuals – not of the pooled measurements, and pooling two treatments that genuinely differ can bend the line while each group on its own is well behaved. So assess it per group:

qq_grp <- function(level)
  plotQQ(subsetDataset(syn, conditions = list(Treatment = level,
                                              Sample_type = "Sample")),
         "Analyte_02", transform = "log2") +
    ggplot2::labs(title = level)

cowplot::plot_grid(qq_grp("PBS"), qq_grp("HDM"), nrow = 1)

Targeted panels are frequently right-skewed, which is why the report renders a tab per transform rather than picking one. The displayed transformation does not change the statistical tests: runStats() works on the values it is handed, so analysing on a transformed scale means transforming the matrix.

PCA

mrmRunPCA() applies the missing-value filtering, imputation, transform, centring and scaling, and returns the fit together with a readable trace of what it did. mrmPlotPCA() draws the scores.

Throughout the package, transform = "log2" means log2(x + 1) – in the PCA and in both heatmaps. The pseudocount keeps zeros and near-zeros finite, but it is an absolute offset on a relative scale: on areas in the thousands it is negligible, while on values already near 1 it compresses the low end noticeably. Changing the units of the matrix therefore changes what this transform does to it, which is worth knowing before comparing a normalised run with an unnormalised one.

Blanks are dropped first. They are diagnostics rather than measurements, and leaving them in gives them a say in the components: a blank sits far from every real sample, so it can dominate the variance the PCA is trying to describe. mrmRunPCA() does not filter them for you, so the subset is explicit.

syn_bio <- subsetDataset(syn, conditions = list(Sample_type = c("Sample", "QC")))

pca <- mrmRunPCA(syn_bio, impute = "gaussian", impute_frac = 0.2,
               transform = "log2", feat_na_max = 0.6,
               sample_na_max = 0.6)
pca$steps

#> [1] "Feature missing filter (max 60% NA): dropped 0 feature(s)."                                                                                
#> [2] "Sample missing filter (max 60% NA): dropped 0 sample(s)."                                                                                  
#> [3] "Imputed NAs with random values between 0 and 0.2 x the per-feature minimum (normal, mean at half that ceiling, SD a sixth of it; seed 42)."
#> [4] "log2(x + 1) transform."                                                                                                                    
#> [5] "Mean-centred per feature."                                                                                                                 
#> [6] "Autoscaled (divide by per-feature SD)."

Colour scores by sample type, batch or injection order to investigate sources of variation. The report draws one plot per column listed under qc_pca.color_by.

mrmPlotPCA(pca, syn_bio, colour_by = "Sample_type")

Colouring by chromatographic batch is the sharpest of these: if the batches separate, batch-associated variation is driving the component and warrants investigation. It does not by itself make that variation analytical – a design where one batch holds more of one treatment, or where batches were run on material collected at different times, separates for reasons that are partly biological. The validation step warns when batch is confounded with a study factor for exactly this reason, and that warning is what tells the two readings apart.

This dataset had a batch shift planted in it, and the design is balanced, so here the separation really is the artefact – the same PCA before and after correctBatch():

keep <- list(Sample_type = c("Sample", "QC"))
before <- subsetDataset(de_uncorr, conditions = keep)
after  <- subsetDataset(de_bc,     conditions = keep)

pca_before <- mrmRunPCA(before, impute = "gaussian", impute_frac = 0.2,
                        transform = "log2")
pca_after  <- mrmRunPCA(after,  impute = "gaussian", impute_frac = 0.2,
                        transform = "log2")

cowplot::plot_grid(
  mrmPlotPCA(pca_before, before, colour_by = "Chrom_Batch") +
    ggplot2::labs(title = "Before batch correction"),
  mrmPlotPCA(pca_after, after, colour_by = "Chrom_Batch") +
    ggplot2::labs(title = "After batch correction"),
  nrow = 1)

Loadings

The scores say whether samples separate; plotLoadings() says which compounds drove it. Colouring by the enzymatic pathway is what turns a ranked list into something interpretable – the palette is shared with the heatmap annotation, so a class is the same colour in every figure.

plotLoadings(pca, syn_bio, colour_by = "Enzymatic_pathway")

What was reported, and what was dropped

DT::datatable(summariseFeatures(syn, out[[2]]), rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

Selecting the biological samples

subsetDataset() slices by any sample_meta condition, which is how the study samples are separated from the QCs and blanks before testing:

samples <- subsetDataset(syn,
                          conditions = list(Sample_type = "Sample"))

table(samples$sample_meta$Treatment, samples$sample_meta$Type)

#>      
#>       KO WT
#>   HDM 10 10
#>   PBS 10 10

Group summaries

summariseGroups() is the table the box and bar plots are drawn from, and it is written to the summary sheet of the stats workbook. Both SD and SE are returned: SD describes the spread of the animals, SE how well the mean is pinned down, and the choice belongs to the reader.

gs <- summariseGroups(samples, "Treatment")
DT::datatable(cbind(gs[, 1:3], round(gs[, 4:6], 3)), rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

Comparisons

runStats() runs the comparisons, correlations, linear models and metabolite ratios described by the stats_report: block of the YAML and returns tidy tables. The method is named explicitly per comparison: welch (the default), student or wilcoxon where a factor has two levels, anova or kruskal where it has three or more, the latter carrying Tukey and pairwise-Wilcoxon post-hoc tests.

Three things to fix before reading any of its output:

• Scale. The test runs on the matrix as handed to runStats() – here the measured, normalised values, untransformed. Nothing is logged on the way in.
• Direction. levels = c("PBS", "HDM") makes PBS group 1 and HDM group 2, and log2FC is log2(mean_HDM / mean_PBS): positive means higher in HDM.
• Families. FDR adjustment is applied within a family, and the family is named in each row’s adjust_family: one per comparison and method for the tests, one per correlation analysis, subset and method, one per linear-model term. A global test and its post-hoc results are therefore adjusted separately – and the post-hoc p-values arrive already adjusted within the feature, by Tukey for anova and by Holm for kruskal, so p_adj on those rows is a second, across-features correction of an already-corrected value.

Naming the test in the YAML, before the data is seen, keeps that choice visible and prevents a more convenient result being selected afterwards.

params <- list(
  sig_threshold = 0.05,
  comparisons = list(enabled = TRUE, entries = list(
    list(name    = "PBS_vs_HDM",
         method  = "welch",
         compare = list(factor = "Treatment",
                        levels = c("PBS", "HDM"))))),
  ion_ratios = list(enabled = TRUE, ratios = list(
    list(num = "Analyte_02", den = "Analyte_05",
         label = "A02/A05"))))

st <- runStats(samples, params)

Because the effects in this dataset were planted, the recovered fold changes can be checked against the truth:

tt <- st$stats[st$stats$method == "welch",
               c("feature", "n_g1", "n_g2", "log2FC", "p_value", "p_adj")]
truth <- as.data.frame(syn$variable_meta)[, c("Compound",
                                              "Simulated_log2FC")]
chk <- merge(tt, truth, by.x = "feature", by.y = "Compound")
chk <- chk[order(-abs(chk$Simulated_log2FC), chk$p_value), ]

chk$log2FC           <- round(chk$log2FC, 3)
chk$Simulated_log2FC <- round(chk$Simulated_log2FC, 3)
chk$p_value          <- signif(chk$p_value, 3)
chk$p_adj            <- signif(chk$p_adj, 3)

DT::datatable(chk, rownames = FALSE, filter = "top",
              options = list(pageLength = 8, scrollX = TRUE))

n_g1 and n_g2 are the non-missing counts the test actually used, which is why they are worth carrying: a fold change computed on three animals is not the same evidence as one computed on twenty.

Read the recovery directly. The five analytes with a non-zero Simulated_log2FC carry the largest recovered effects and the smallest p-values, in the right direction, and four of them clear FDR 0.05. Analyte_18 – the smallest planted effect, at -0.737 – does not, which is the expected outcome rather than a failure: an effect that size, in this number of samples, is not reliably recoverable at this threshold. The other fifteen scatter around zero.

Per-compound plots

plotBoxplot(samples, "Analyte_02", "Treatment",
             value_label = "response / ug protein")

type = "bar" draws the group means and error bars from summariseGroups(), so the figure and the summary sheet cannot drift apart:

plotBoxplot(samples, "Analyte_02", "Treatment", type = "bar",
             error_bar = "SE", value_label = "response / ug protein")

Volcano

Volcano plots show fold changes against statistical significance. The x axis is log2(HDM / PBS), so compounds to the right are higher in HDM.

The inference above used FDR-adjusted p-values, so the plot should too – use_adjusted = TRUE draws the same threshold the table was read against:

plotVolcano(st, "PBS_vs_HDM", use_adjusted = TRUE, sig_threshold = 0.05)

Whether a small fold change matters depends on the analyte and the question, not on the plot.

Heatmap

Where the comparisons test one compound at a time, the heatmap shows the whole panel at once – which is how you see that a treatment moved an entire pathway rather than a few compounds that happened to clear a threshold.

mrmPlotHeatmap(samples,
             color_samples_by  = "Treatment",
             group_features_by = "Enzymatic_pathway",
             transform         = "log2",
             title             = "All analytes")

plotGroupHeatmap(samples,
                   color_samples_by  = "Treatment",
                   group_features_by = "Enzymatic_pathway",
                   transform         = "log2",
                   title             = "All analytes, group means")

Feature correlations

Feature correlations identify compounds that vary together. Shared treatment responses, batch effects and normalisation can all produce them, so a block along the diagonal is a starting point for a hypothesis about shared pathway chemistry rather than evidence of it – correlating within a treatment group as well as across the pooled samples is the quickest way to see which it is.

cp <- list(correlations = list(enabled = TRUE, entries = list(
  list(name = "pairs", methods = "spearman",
       subsets = list(list())))))

plotCorrelation(runStats(samples, cp), samples,
                 group_by = "Enzymatic_pathway")

Metabolite ratios

Metabolite ratios describe relative abundances and can support hypotheses about pathway regulation (these are ratios between metabolites, not the qualifier-to-quantifier transition ratios). A product-to-precursor ratio cancels shared variation in how much substrate was present, which makes it a more stable quantity than either compound alone.

It is not a measure of enzyme activity or of flux. A concentration ratio at one timepoint is the net result of synthesis, degradation, transport and compartmentation; inferring a rate from it needs kinetic information this assay does not carry. runStats() returns the per-sample ratios with a group test attached, so they reach the ion_ratios sheet as well as the report.

ir <- st$ion_ratios[!is.na(st$ion_ratios$comparison), ]
DT::datatable(cbind(ir[, c("ratio", "comparison", "sample", "group")],
         round(ir[, c("value", "p_value")], 4)),
   rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))