Scope, and where this sits among related packages

MRManalyzeR begins with reviewed vendor result tables. It does not perform chromatographic peak detection, peak integration or calibration-curve fitting from raw mass-spectrometry data; those judgements stay in TargetLynx or Skyline, where the chromatograms are.

That is the line which separates it from the Bioconductor packages nearest to it, and from the existing targeted-MRM tools:

The design assumes a validated targeted panel: roughly 30 MRM transitions upward, known compounds carrying class or pathway annotation, authentic standards, and an experiment designed to test a hypothesis rather than to discover features. Feature metadata is treated as biology, not bookkeeping – compound class and enzymatic pathway group and colour the heatmap annotations, PCA loadings and correlation blocks throughout.

Assumptions and limitations

Worth knowing before you rely on the output.

• Quantification is internal-standard based. MRManalyzeR deliberately provides conservative normalisation and batch-adjustment methods suited to validated targeted assays, where an isotope-labelled standard has already absorbed most of the run-order and matrix variation. More extensive signal-drift modelling – QC-RSC spline correction and similar – remains available through complementary packages such as pmp, and is the right choice when the assay does not carry standards of that kind.
• Batch correction equalises batch medians. It assumes the reference samples (QC, or biological samples if the design is balanced across batches) are comparable between batches. correct_batch() refuses to run without enough reference injections per batch, and warns when the batch variable is confounded with a study factor – but it cannot detect a design where that confounding is total.
• Concentration adjustment rescales vendor-reported values. It applies the internal-standard concentration ratio and the sample-volume correction; it does not fit or refit a calibration curve.
• Imputation is a fixed fraction of each compound’s observed minimum. That is a convention, not an estimate. Values below the detection limit and values missing for other reasons are not currently distinguished.
• Statistical methods are chosen by you, per comparison. Nothing is selected automatically from the data, and no test is run in parallel with another so that the more convenient answer can be picked afterwards.
• Multiple-testing correction is applied within a family, never across the whole workbook – see the adjust_family column of the stats output.

How the package is organised

The workflow falls into four stages. Everything from stage 1 onwards is carried in a single struct::DatasetExperiment – the assay matrix, the sample metadata and the feature metadata kept aligned in one object – so every function takes that object and returns it, and nothing has to be handed a matching metadata table a second time.

   project_config.yml            workbook.xlsx
           |                           |
           |                  read_targetlynx() / read_skyline()   [1]
           |                           |
           |                  wide sample x compound matrix
           |                           |
           |                    assemble_dataset()
           |                           |
           |                           v
           |                   DatasetExperiment
           |                           |
           |          filter_blanks -> normalise_matrix ->         [2]
           +--------> adjust_concentration -> impute_missing
                      then correct_batch
                                       |
                                       v
                          processed DatasetExperiment
                            |          |          |
             [3] QC check   |          |          |   [4] Stats
             add_cv_metrics |          |          |   summarise_groups
             run_pca        |          |          |   run_stats
             summarise_     |          |          |
               features     |          |          |
                            |          |          |
              tables + plots|          |          |tables + plots
                            v          v          v
              data_quality       .xlsx  .RDS      stats_report
                _report.html                      .html
                                       |
                                  _stats.xlsx

Every number in a figure comes from the same table that is written to the spreadsheet – run_stats(), summarise_groups() and add_cv_metrics() compute once, the plot_*() functions consume. A figure and the sheet beside it cannot disagree.

Each plot is an exported function rather than code buried in a report template, and that is deliberate. Most users will never call one: the YAML route renders both reports and the figures appear. But when a figure needs to go into a paper with different labels, or one compound needs looking at on its own, or a reviewer asks for the same plot on a subset, the choice should not be between editing an R Markdown template and starting from the raw matrix. plot_volcano(st, "PBS_vs_HDM") is one line, takes the objects the workflow already produced, and returns an ordinary ggplot you can add layers to, restyle or save. The reports call exactly the same functions, so what you get interactively is what the report shows.

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(as.data.frame(table(Sample_type = meta$Sample_type)),
   rownames = FALSE)

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

dt(fdata[fdata$Report == "YES",
         c("Compound", "Enzymatic_pathway", "PUFA", "RT", "Report")],
   rownames = FALSE)

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 awkward to read by eye and worse to read programmatically, which is the first problem the package solves.

read_targetlynx() 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 <- read_targetlynx(xlsx, datatype = "Area", snr = 3)
dim(m)

#> [1]  82 155

dt(round(m[, 1:6], 1))

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.

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

The simulated dataset used from here on

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 no-ops on it. 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.

Nothing in it is measured data. It is generated from a fixed seed by data-raw/make_example_synthetic.R (which ships with the package source), and what it contains was deliberately planted so each step can be seen recovering something known: a response shift between the two batches, a treatment effect on 5 of the 20 analytes – recorded in feature_metadata$Simulated_log2FC – and QC injections drawn with about 8 percent noise against 35 percent between animals. 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

assemble_dataset() 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 <- read_skyline(sxlsx, "skyline_data", sfeat, signal_filter = "LOD")
de <- assemble_dataset(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

Blank filter

filter_blanks() 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.

de <- filter_blanks(de, blank_filter = 3, blank_head = "Sample_type",
                    blank_name = "Blank")
sum(is.na(de$data))

#> [1] 154

Normalise

normalise_matrix() 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 <- normalise_matrix(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      3.134     22.430     68.330
#> SYN_003      9.044     11.367     50.775

Concentration adjustment

A calibration curve returns the concentration in the vial the instrument sampled, and it does so on the calibrant’s terms: the curve was built at the calibrant’s internal-standard concentration, and it silently assumes your sample carries the same one. adjust_concentration() corrects that assumption:

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

Those two ratios are one thing, not two. Multiply them out and they collapse to a single quantity:

\[ \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. The reconstitution volume matters only through the IS concentration, which is exactly what the expression above captures. (This reasoning depends on quantification being IS-ratio based, as it is for a targeted assay with stable-isotope standards. Pointing adjust_conc at a raw Area datatype with no IS correction would break the assumption.)

Setting starting_vol_col adds a second step, back to the concentration in the original, undried sample – the plasma or lavage volume actually taken into extraction – so results are reported per unit of biofluid rather than per vial. Following one analyte amount A through both steps:

  curve gives     conc_reported = (A / IS_sample) x [IS] in calibrant
  step 1 gives    conc_vial     = A / sample_vol
  step 2 gives    conc_original = A / starting_vol

sample_vol cancels between the two steps, and it should: the final number depends only on how much analyte was there and the volume of biofluid it came from, never on how much solvent you happened to reconstitute in.

de <- adjust_concentration(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.78343    5.60746   17.08246
#> SYN_003    2.26112    2.84183   12.69373

Impute missing values

impute_missing() replaces the remaining NAs in the biological samples with a fraction of that compound’s minimum observed value, treating a non-detect as lying just below the lowest measurement rather than as a zero or a hole. scalar = 0.5 is the familiar half-minimum rule; a smaller fraction (0.2) is more conservative if the panel has compounds sitting near the detection limit. Blanks are left untouched – they are diagnostics, not measurements.

Other strategies are available where they matter most, in the PCA pre-processing (run_pca(impute = )), which offers "none", "min", "half_min" and "frac_min" alongside missing-value filters that drop features or samples above a missingness threshold instead of imputing them. For targeted panels the choice is rarely critical: a compound missing in most samples is usually better dropped than modelled.

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

#> [1] 80

Batch correction

Where an acquisition spans several days or columns, response can shift between batches. correct_batch() scales each batch so its per-feature median matches the grand median across the reference samples:

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

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

Taking the first analyte as an example, the two batch medians are pulled together:

med <- function(x) tapply(as.data.frame(x$data)[[1]], batch, median)

rbind(before = med(de_uncorr),
      after  = med(de))

#>        B1 B2
#> before NA NA
#> after  NA NA

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.

The whole matrix in one call

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

out <- process_dataset(
  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.5,
  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 processed dataset also ships ready-made, which is what the rest of this vignette uses:

syn <- load_dataset(system.file("extdata", "example_synthetic.RDS",
                                package = "MRManalyzeR"))
dim(syn$data)

#> [1] 54 20

How well was each compound measured?

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

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

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

dt(cv[order(cv$CV_QC), ], rownames = FALSE, page = 10)

• CV_QC is the coefficient of variation across the pooled QC injections. Every QC is the same material, so this is technical variability – the measurement noise of that compound in that run. It is the primary quantification filter: a common rule of thumb keeps compounds under 20-30 percent CV_QC and reports the rest as detected but not reliably quantified. It is NA when a run contains no QC injections.
• 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 ratio comfortably above 1 means the biological signal is larger than the noise it is measured against – the compound has room to show a real difference. A ratio near 1 means a compound is varying no more between animals than between replicate injections of identical material: whatever passed the CV_QC filter, there is little biology left to find in it. Sorting the panel by this ratio is a quick way to see which analytes are worth interpreting before running any test.

Analytical drift

plot_drift() 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.

plot_drift(syn, "Analyte_02", value_label = "ng/mL")

Normality

plot_qq() asks whether a compound’s distribution supports the parametric tests. Points along the line mean the t-test is on safe ground; a curve means the Wilcoxon result reported alongside it is the one to trust. Targeted panels are frequently right-skewed, which is why the report renders a tab per transform rather than picking one.

plot_qq(syn, "Analyte_02", transform = "log2", levels = "Sample")

PCA

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

pca <- run_pca(syn, impute = "frac_min", 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 per-feature 0.2 x minimum."               
#> [4] "log2(x + 1) transform."                                    
#> [5] "Mean-centred per feature."                                 
#> [6] "Autoscaled (divide by per-feature SD)."

The colouring is the whole point: the same scores answer a different question depending on what you colour them by. The report draws one plot per column listed under qc_pca.color_by.

plot_pca(pca, syn, colour_by = "Sample_type")

Colouring by chromatographic batch is the sharpest of these – if the batches separate, the variance driving component 1 is analytical rather than biological. This dataset had a batch shift planted in it, so here is the same PCA before and after correct_batch():

pca_by_batch <- function(x, title){
  p <- run_pca(x, impute = "frac_min", impute_frac = 0.2,
               transform = "log2")
  plot_pca(p, x, colour_by = "Chrom_Batch") +
    ggplot2::labs(title = title)
}

cowplot::plot_grid(
  pca_by_batch(de_uncorr, "Before batch correction"),
  pca_by_batch(de,        "After batch correction"),
  nrow = 1)

Loadings

The scores say whether samples separate; plot_loadings() 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.

plot_loadings(pca, syn, colour_by = "Enzymatic_pathway")

What was reported, and what was dropped

dt(summarise_features(syn, out[[2]]), rownames = FALSE)

Selecting the biological samples

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

samples <- subset_dataset(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

summarise_groups() 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 <- summarise_groups(samples, "Treatment")
dt(cbind(gs[, 1:3], round(gs[, 4:6], 3)), rownames = FALSE)

Comparisons

run_stats() runs the comparisons, correlations, linear models and ion ratios described by the stats_report: block of the YAML and returns tidy tables. The method is chosen from the design rather than by the user: two levels get a t-test and a Wilcoxon test, three or more get ANOVA and Kruskal-Wallis with Tukey and pairwise-Wilcoxon post-hoc tests, and p-values are FDR-adjusted.

Reporting a parametric and a non-parametric test side by side is deliberate: targeted panels are frequently run on small groups where normality cannot be taken on trust, and agreement between the two is itself a piece of evidence.

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 <- run_stats(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", "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), ]

dt(cbind(chk[, 1], round(chk[, -1], 4)), rownames = FALSE, page = 8)

The five analytes with a non-zero Simulated_log2FC should carry the largest recovered effects and the smallest p-values; the other fifteen should scatter around zero.

Per-compound plots

plot_boxplot(samples, "Analyte_02", "Treatment",
             value_label = "ng/mL")

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

plot_boxplot(samples, "Analyte_02", "Treatment", type = "bar",
             error_bar = "SE", value_label = "ng/mL")

Volcano

Effect size against significance. A large fold change on three noisy animals is not a finding, and a tiny but reproducible shift is rarely interesting either; the compounds worth attention sit in the upper corners.

plot_volcano(st, "PBS_vs_HDM", sig_threshold = 0.05)

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.

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

Feature correlations

Which compounds move together, which in a targeted panel usually means they share a branch of the pathway – so blocks along the diagonal are the pathway structure showing itself.

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

plot_correlation(run_stats(samples, cp), samples,
                 group_by = "Enzymatic_pathway")

Ion ratios

An ion ratio is an enzyme-activity surrogate: the ratio of a product to its precursor tracks flux through a pathway more directly than either compound alone, because it cancels the shared variation in how much substrate was there. run_stats() 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(cbind(ir[, c("ratio", "comparison", "sample", "group")],
         round(ir[, c("value", "p_value")], 4)),
   rownames = FALSE)

Merging several acquisition panels

Where a study is split across methods – a polar-acid and a polar-basic panel, say – each is processed on its own with run_MRManalyzeR() first, and the resulting .RDS files are merged:

file.copy(
  system.file("extdata", "example_combine_config.yml",
              package = "MRManalyzeR"),
  "combine_config.yml")

res <- run_MRManalyzeR_combine("combine_config.yml")

Two things decide whether a merge succeeds: sample_id_col must name an identifier that is stable across the panels (not the acquisition Name, which usually carries a method-specific prefix), and each dataset’s tag prefixes its feature names so that a compound measured on two methods does not collide.