--- title: "Processing targeted LC-MS/MS lipidomics and metabolomics with MRManalyzeR" author: - name: Matthew J. Smith affiliation: Karolinska Institutet, Stockholm, Sweden email: mattyjsmith123@gmail.com date: "`r Sys.Date()`" package: MRManalyzeR output: BiocStyle::html_document: toc_float: true number_sections: true vignette: > %\VignetteIndexEntry{Processing targeted LC-MS/MS lipidomics and metabolomics with MRManalyzeR} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup-knitr, include = FALSE} knitr::opts_chunk$set(comment = "#>", fig.width = 6.5, fig.height = 4, dpi = 100, message = FALSE, warning = FALSE) # Tables are shown the way the HTML reports show them: searchable, # sortable and scrollable rather than truncated with head(). dt <- function(x, rownames = TRUE, page = 5){ DT::datatable(x, rownames = rownames, filter = "top", options = list(pageLength = page, scrollX = TRUE, dom = "ftip")) } ``` # Introduction Targeted LC-MS/MS is analytically less complex than untargeted profiling. The compound panel is defined in advance, every analyte has an authentic standard, and the questions that remain are about **data quality** and the **accurate reporting of known metabolites** rather than about feature discovery. Normalisation to stable-isotope internal standards is commonplace, and reported values are true concentrations only once they have been adjusted for the internal-standard concentration each vial actually carried and for the volume of sample the extract came from. MRManalyzeR is built around exactly that: assessing data quality on the way to accurate quantification, and giving a basic interpretation of the analytes that come out. The output is deliberately visual -- quality control and statistics are presented as plots with simple statistical summaries beside them, not as a modelling framework. Everything is driven from a curated `.xlsx` workbook and a single YAML file, so a run is reproducible and the workflow stays accessible to analysts who do not write R. The package reads Waters **TargetLynx** exports (the primary format), **Skyline** exports, and plain sample-by-analyte matrices from any other software, and covers signal filtering, blank filtering, normalisation, concentration adjustment, missing-value imputation and batch correction, followed by quality control, statistics and two self-contained HTML reports. ## 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: | Task | Tool | |---|---| | Raw untargeted LC-MS feature detection and alignment | `xcms` | | General metabolomics peak-matrix preprocessing | `pmp` | | Metabolite annotation-table workflows | `MetMashR` | | Chromatogram processing, peak integration, automated identification | MRMAnalyzer, MRMkit, MRMQuant | | Post-acquisition targeted-assay processing, QC and reporting | **MRManalyzeR** | 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. | Stage | Question it answers | Functions | |---|---|---| | 1. Data parse | What did the instrument report? | `read_targetlynx()`, `read_skyline()`, `assemble_dataset()`, `load_config()`, `load_dataset()`, `combine_datasets()` | | 2. Peak-matrix processing | What is the true concentration? | `filter_blanks()`, `normalise_matrix()`, `adjust_concentration()`, `impute_missing()`, `correct_batch()`, `subset_dataset()`, `process_dataset()` | | 3. QC check | Can these numbers be trusted? | `add_cv_metrics()`, `summarise_features()`, `run_pca()`, `plot_drift()`, `plot_qq()`, `plot_pca()`, `plot_loadings()` | | 4. Stats | What do these numbers mean? | `run_stats()`, `summarise_groups()`, `plot_boxplot()`, `plot_volcano()`, `plot_heatmap()`, `plot_correlation()` | ``` 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. # Installation ```{r install, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("MRManalyzeR") ``` The development version: ```{r install-dev, eval = FALSE} BiocManager::install("MJS-708/MRManalyzeR") ``` ```{r load} library(MRManalyzeR) ``` # Quick start `run_example()` runs all four stages on the bundled dataset and writes its outputs to a temporary directory. `render = FALSE` builds the peak matrix and the xlsx/RDS outputs but skips the (slower) HTML report rendering; drop it to get the reports as well. Running this successfully is the recommended check that an installation works before pointing the package at real data. ```{r run-example} res <- run_example(render = FALSE, open = FALSE) list.files(res$output_directory) ``` Everything after this section unpacks what that single call did. # Stage 1: data parse ## The curated workbook MRManalyzeR expects one `.xlsx` workbook holding the instrument export plus the two curation sheets that describe it: | Sheet | One row per | Purpose | |---|---|---| | `lcms_data` (or `skyline_data`) | peak | the instrument export itself | | `feature_metadata` | compound | which compounds to report, and their annotation | | `sample_metadata` | injection | which injections to include, and their design | 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. ```{r dataset} 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)) ``` The bundled workbook is a subset of the oxylipin LC-MS/MS data of Kolmert *et al.* (2018) (see [References](#references)): 79 injections of mouse bronchoalveolar lavage over a calibration ladder, extraction blanks, pooled QC injections and 45 study samples. ```{r dataset-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. ```{r dataset-features} 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: ```{r read} m <- read_targetlynx(xlsx, datatype = "Area", snr = 3) dim(m) ``` ```{r read-table} 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. ```{r syn-data} 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) ``` ## 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. ```{r step-assemble} sm <- read_skyline(sxlsx, "skyline_data", sfeat, signal_filter = "LOD") de <- assemble_dataset(sm, sfeat, smeta) de ``` # Stage 2: peak-matrix processing ## 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. ```{r step-blanks} de <- filter_blanks(de, blank_filter = 3, blank_head = "Sample_type", blank_name = "Blank") sum(is.na(de$data)) ``` ## 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: ```{r step-normalise} de <- normalise_matrix(de, column = "protein_ug") round(as.data.frame(de$data)[1:3, 1:3], 3) ``` ## 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. ```{r step-adjust} 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) ``` ## Impute missing values `impute_missing()` replaces the remaining `NA`s 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. ```{r step-impute} de <- impute_missing(de, scalar = 0.5, blank_head = "Sample_type", blank_name = "Blank") sum(is.na(de$data)) ``` ## 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: ```{r step-batch} 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: ```{r step-batch-medians} med <- function(x) tapply(as.data.frame(x$data)[[1]], batch, median) rbind(before = med(de_uncorr), after = med(de)) ``` 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: ```{r step-process} 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) ``` `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: ```{r syn-load} syn <- load_dataset(system.file("extdata", "example_synthetic.RDS", package = "MRManalyzeR")) dim(syn$data) ``` # Stage 3: QC check ## 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. ```{r qc-cv} 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. ```{r qc-drift} 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. ```{r qc-qq, fig.width = 5, fig.height = 4} 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. ```{r qc-pca} 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 ``` 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`. ```{r qc-pca-type, fig.height = 4.2} 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()`: ```{r qc-pca-batch, fig.height = 3.4} 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. ```{r qc-loadings, fig.width = 7, fig.height = 5} plot_loadings(pca, syn, colour_by = "Enzymatic_pathway") ``` ## What was reported, and what was dropped ```{r qc-features} dt(summarise_features(syn, out[[2]]), rownames = FALSE) ``` # Stage 4: statistics ## 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: ```{r stats-subset} samples <- subset_dataset(syn, conditions = list(Sample_type = "Sample")) table(samples$sample_meta$Treatment, samples$sample_meta$Type) ``` ## 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. ```{r stats-summary} 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. ```{r stats-run} 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: ```{r stats-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 ```{r stats-box, fig.width = 5, fig.height = 4} 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: ```{r stats-bar, fig.width = 5, fig.height = 4} 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. ```{r stats-volcano, fig.width = 5.5, fig.height = 5} 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. ```{r stats-heatmap, fig.width = 7.5, fig.height = 6} 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. ```{r stats-corr, fig.width = 6.5, fig.height = 6} 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. ```{r stats-ratios} 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) ``` # Driving it from a YAML config In routine use none of the above is typed out. One YAML file describes the data source, every processing step and both reports, and one call runs it: ```{r run-yaml, eval = FALSE} res <- run_MRManalyzeR("path/to/your_config.yml") ``` The bundled reference configuration is a commented template to copy and edit: ```{r config-path} system.file("extdata", "example_config.yml", package = "MRManalyzeR") ``` Its `PeakMatrixProcessing:` block is a one-to-one map of stage 2 -- each key switches a step on or names the metadata column it should use: ```yaml PeakMatrixProcessing: execute: True # False = skip, reuse existing RDS name_col: Name # sample id (matches the matrix rows) include_col: Include # sample include-filter column include_value: "YES" tl_data: # TargetLynx source (or skyline_data:) enabled: True data_tab_names: [lcms_data] tl_headers: ["ID", "Name", "Area", "ng/mL", "Response", "S/N"] snr: 3 # S/N threshold | False blank_filter: 3 # n x blank mean | False blank_head: 'Sample_type' blank_name: 'Blank' normalize: protein_ug # False | sample_metadata column adjust_conc: enabled: True starting_vol_col: starting_vol_uL 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 replace_MVs: 0.5 # False | scalar x per-feature min batch_correction: True bc_qc_label: "QC" bc_factor_name: "Sample_type" bc_header: "Chrom_Batch" ``` `data_quality_report:` and `stats_report:` configure stages 3 and 4 in the same declarative way -- which column identifies QCs, which factors to compare, which plots to draw. See the bundled file for the fully commented set. Outputs are written next to your data: the processed matrix as `.xlsx` and `.RDS`, a `_stats.xlsx` carrying `stats`, `correlations`, `linear_models`, `ion_ratios` and `summary`, the two self-contained HTML reports, and a text file recording the exact parameters used. The stats workbook opens on a `key` sheet defining every column heading it contains -- what `adjust_family` covers, how `log2FC` is signed, which effect size `effect_type` names. The workbook usually outlives the session that produced it, so it carries its own glossary rather than relying on this vignette. ## 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: ```{r run-combine, eval = FALSE} 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. # References {#references} Kolmert J. *et al.* (2018). Prominent release of lipoxygenase generated mediators in a murine house dust mite-induced asthma model. *Prostaglandins & Other Lipid Mediators*. doi:[10.1016/j.prostaglandins.2018.05.005](https://doi.org/10.1016/j.prostaglandins.2018.05.005) -- source of the bundled example data. Lloyd G.R., Jankevics A., Weber R.J.M. (2020). struct: an R/Bioconductor-based framework for standardized metabolomics data analysis. *Bioinformatics*. doi:[10.1093/bioinformatics/btaa227](https://doi.org/10.1093/bioinformatics/btaa227) -- the `DatasetExperiment` container used throughout. # Session info ```{r session-info} sessionInfo() ```