--- title: "MSstats: Metabolomics workflow with MZMine" date: June 17th, 2026 --- ```{r style, echo = FALSE, results = 'asis'} BiocStyle::markdown() ``` ```{r global_options, include=FALSE} knitr::opts_chunk$set(fig.width=10, fig.height=7, warning=FALSE, message=FALSE) options(width=110) ``` ```{=html} ``` # __MSstats: Metabolomics workflow with MZMine__ Author: Swaraj Patil Date: June 17th, 2026 ## __Introduction__ `MSstats` supports differential analysis of metabolomics data acquired with LC-MS untargeted workflows. This vignette walks an end-to-end run: import MZMine feature quantifications and library annotations, layer in SIRIUS structure identifications, convert to the MSstats format, summarize features into compound-level abundance, and test for differences between conditions. Compound identification combines two evidence sources: * __MZMine compound names__ come from MS/MS spectral-library matching and correspond to MSI Level 2 putative identifications (Sumner et al., 2007). * __SIRIUS names__ come from in-silico structure prediction and correspond to MSI Level 3 identifications. The SIRIUS pass extends discovery coverage to features the MZMine spectral library does not cover. `MZMinetoMSstatsFormat` is re-exported from `MSstatsConvert`, so attaching `MSstats` alone is enough to run the full workflow. ## __1. Setup__ ```{r setup} library(MSstats) library(data.table) ``` ## __2. Load example data__ Example MZMine input, sample annotation, MZMine library annotations, and SIRIUS structure identifications ship with `MSstatsConvert` and are loaded via `system.file()`. ```{r load-data} input_path = system.file("tinytest/raw_data/MZMine/mzmine_input.csv", package = "MSstatsConvert") annotation_path = system.file("tinytest/raw_data/MZMine/annotation.csv", package = "MSstatsConvert") mzmine_ann_path = system.file("tinytest/raw_data/MZMine/mzmine_annotations.csv", package = "MSstatsConvert") sirius_path = system.file("tinytest/raw_data/MZMine/structure_identifications.tsv", package = "MSstatsConvert") mzmine_input = data.table::fread(input_path) annotation = data.table::fread(annotation_path) mzmine_annotations = data.table::fread(mzmine_ann_path) sirius_annotations = data.table::fread(sirius_path) head(mzmine_input, 5) head(annotation) head(mzmine_annotations) head(sirius_annotations) ``` The MZMine feature table is wide: one row per feature, columns `row ID`, `row m/z`, `row retention time`, and per-sample `" Peak area"` columns. The annotation table maps each MS run to its `Condition` and `BioReplicate`. `mzmine_annotations` is the spectral-library match table (`id`, `compound_name`, `score`, `adduct`); features with multiple library hits resolve to the highest-scoring compound. `sirius_annotations` is SIRIUS's `structure_identifications.tsv`; its `mappingFeatureId` joins to `row ID` in the MZMine input. ## __3. Convert with `MZMinetoMSstatsFormat`__ ```{r convert, message = FALSE} mzmine_msstats = MZMinetoMSstatsFormat( input = mzmine_input, annotation = annotation, mzmine_annotations = mzmine_annotations, sirius_annotations = sirius_annotations, use_log_file = FALSE ) head(mzmine_msstats) ``` `ProteinName` is assigned per feature in priority order: (1) the highest-scoring MZMine compound name when present, (2) the SIRIUS name when MZMine has no match, (3) an `m/z_RT` fallback identifier for features neither source identified. Every feature is retained -- discovery coverage is preserved at the cost of a wider multiple-testing burden in Section 5. Although the column is named `ProteinName` for compatibility with the rest of MSstats, here it holds the compound (analyte) name; a future release will expose it as `Analyte` for metabolomics data. ## __4. Summarize with `dataProcess`__ ```{r summarize, message = FALSE} summarized = dataProcess( mzmine_msstats, logTrans = 2, normalization = "equalizeMedians", featureSubset = "all", summaryMethod = "TMP", censoredInt = "NA", MBimpute = TRUE, use_log_file = FALSE ) head(summarized$FeatureLevelData) head(summarized$ProteinLevelData) ``` The settings above mirror a typical discovery proteomics workflow: log-2 transform, median-equalized normalization, all features used, and Tukey median polish summarization. Model-based imputation is enabled (`MBimpute = TRUE`), but no values are imputed in this small example. Caffeine is detected at two adducts (`[M+H]+` on feature 1, `[M+Na]+` on feature 6) and is summarized into a single compound-level abundance per run. ## __5. Test for differences with `groupComparison`__ We construct the contrast matrix by hand to show how a comparison is specified. The two-condition design admits a single Treatment-vs-Control contrast: ```{r contrast, message = FALSE} # A contrast matrix has one row per comparison and one column per condition. # Columns must match the condition levels (alphabetical here: Control, Treatment). # The -1 / +1 pair selects the two groups being compared. contrast_matrix = matrix(c(-1, 1), nrow = 1) colnames(contrast_matrix) = c("Control", "Treatment") rownames(contrast_matrix) = "Treatment vs Control" contrast_matrix comparison = groupComparison(contrast.matrix = contrast_matrix, data = summarized, use_log_file = FALSE) comparison$ComparisonResult ``` Each row of `ComparisonResult` is one compound (or `m/z_RT` fallback) tested against the contrast. Columns of interest: `log2FC`, `pvalue`, and `adj.pvalue`. The `issue` column flags compounds that could not be tested normally, for example one missing from an entire condition; in this small fixture the `issue` column is `NA` for every compound shown. ## __6. Visualization__ Profile plots show feature-level intensities alongside the compound-level summary. Caffeine is identified at two adducts in this dataset and is summarized into a single compound -- the profile plot makes that aggregation visible. ```{r profile, fig.width = 8, fig.height = 5} dataProcessPlots(summarized, type = "ProfilePlot", which.Protein = "Caffeine", address = FALSE) ``` For a study-wide view of fold-change versus significance, pass the `groupComparison` result to `groupComparisonPlots`. On a four-sample fixture the volcano is sparse. ```{r volcano, eval = FALSE} groupComparisonPlots(data = comparison$ComparisonResult, type = "VolcanoPlot", address = FALSE) ``` ## __References__ Sumner LW, Amberg A, Barrett D, et al. (2007). Proposed minimum reporting standards for chemical analysis: Chemical Analysis Working Group (CAWG) Metabolomics Standards Initiative (MSI). *Metabolomics* 3(3): 211-221. doi: [10.1007/s11306-007-0082-2](https://doi.org/10.1007/s11306-007-0082-2) ```{r session} sessionInfo() ```