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
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().
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)
## row ID row m/z row retention time sampleA.mzML Peak area sampleB.mzML Peak area sampleC.mzML Peak area
## <int> <num> <num> <int> <int> <int>
## 1: 1 123.056 1.23 1000 1100 1200
## 2: 2 245.129 3.45 5000 4800 5200
## 3: 3 367.201 5.67 800 0 750
## 4: 4 489.334 7.89 2000 2100 1900
## 5: 5 555.447 9.10 100 0 0
## sampleD.mzML Peak area
## <int>
## 1: 1300
## 2: 4900
## 3: 820
## 4: 2050
## 5: 0
head(annotation)
## Run Condition BioReplicate
## <char> <char> <int>
## 1: sampleA.mzML Control 1
## 2: sampleB.mzML Control 2
## 3: sampleC.mzML Treatment 3
## 4: sampleD.mzML Treatment 4
head(mzmine_annotations)
## id compound_name score adduct
## <int> <char> <num> <char>
## 1: 1 Caffeine 0.95 [M+H]+
## 2: 2 GlucoseLow 0.72 [M+H]+
## 3: 2 GlucoseHigh 0.91 [M-H]-
## 4: 3 Lactate 0.88 [M+H]+
## 5: 6 Caffeine 0.80 [M+Na]+
head(sirius_annotations)
## mappingFeatureId name ConfidenceScoreExact ConfidenceScoreApproximate SiriusScore
## <int> <char> <num> <num> <num>
## 1: 1 DuplicateFromSirius 0.30 0.40 5.5
## 2: 4 Caffeic acid 0.85 0.88 22.1
## 3: 5 0.10 0.12 1.0
## 4: 99 Ghost 0.50 0.55 8.0
The MZMine feature table is wide: one row per feature, columns row ID,
row m/z, row retention time, and per-sample "<run> 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
mzmine_msstats = MZMinetoMSstatsFormat(
input = mzmine_input,
annotation = annotation,
mzmine_annotations = mzmine_annotations,
sirius_annotations = sirius_annotations,
use_log_file = FALSE
)
## INFO [2026-08-05 19:08:29] ** Raw data from MZMine imported successfully.
## INFO [2026-08-05 19:08:29] ** MZMine ProteinName assignment: MZMine compound: 4 feature(s); SIRIUS name: 1 feature(s); m/z-RT fallback: 1 feature(s).
## INFO [2026-08-05 19:08:29] ** Raw data from MZMine cleaned successfully.
## INFO [2026-08-05 19:08:29] ** Using provided annotation.
## INFO [2026-08-05 19:08:29] ** Run labels were standardized to remove symbols such as '.' or '%'.
## INFO [2026-08-05 19:08:29] ** The following options are used:
## - Features will be defined by the columns: PeptideSequence, PrecursorCharge, FragmentIon, ProductCharge
## - Shared peptides will not be removed.
## - Proteins with single feature will not be removed.
## - Features with less than 3 measurements across runs will be kept.
## INFO [2026-08-05 19:08:29] ** Features with all missing measurements across runs are removed.
## INFO [2026-08-05 19:08:29] ** Multiple measurements in a feature and a run are summarized by summaryforMultipleRows: max
## INFO [2026-08-05 19:08:29] ** Features with all missing measurements across runs are removed.
## INFO [2026-08-05 19:08:29] ** Run annotation merged with quantification data.
## INFO [2026-08-05 19:08:29] ** Updated quantification data to make balanced design. Missing values are marked by NA
## INFO [2026-08-05 19:08:29] ** Finished preprocessing. The dataset is ready to be processed by the dataProcess function.
head(mzmine_msstats)
## ProteinName PeptideSequence PrecursorCharge FragmentIon ProductCharge IsotopeLabelType Condition
## 1 Caffeine 1 NA <NA> NA Light Control
## 2 Caffeine 1 NA <NA> NA Light Control
## 3 Caffeine 1 NA <NA> NA Light Treatment
## 4 Caffeine 1 NA <NA> NA Light Treatment
## 5 GlucoseHigh 2 NA <NA> NA Light Control
## 6 GlucoseHigh 2 NA <NA> NA Light Control
## BioReplicate Run Fraction Intensity
## 1 1 sampleAmzML 1 1000
## 2 2 sampleBmzML 1 1100
## 3 3 sampleCmzML 1 1200
## 4 4 sampleDmzML 1 1300
## 5 1 sampleAmzML 1 5000
## 6 2 sampleBmzML 1 4800
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
summarized = dataProcess(
mzmine_msstats,
logTrans = 2,
normalization = "equalizeMedians",
featureSubset = "all",
summaryMethod = "TMP",
censoredInt = "NA",
MBimpute = TRUE,
use_log_file = FALSE
)
## INFO [2026-08-05 19:08:29] ** Log2 intensities under cutoff = 7.93 were considered as censored missing values.
## INFO [2026-08-05 19:08:29] ** Log2 intensities = NA were considered as censored missing values.
## INFO [2026-08-05 19:08:29] ** Use all features that the dataset originally has.
## INFO [2026-08-05 19:08:29]
## # proteins: 5
## # peptides per protein: 1-2
## # features per peptide: 1-1
## INFO [2026-08-05 19:08:29] Some proteins have only one feature:
## 555.447_9.1,
## Caffeic acid,
## GlucoseHigh,
## Lactate ...
## INFO [2026-08-05 19:08:29]
## Control Treatment
## # runs 2 2
## # bioreplicates 2 2
## # tech. replicates 1 1
## INFO [2026-08-05 19:08:29] Some features are completely missing in at least one condition:
## 5_NA_NA_NA,
## NA ...
## INFO [2026-08-05 19:08:29] == Start the summarization per subplot...
##
|
| | 0%
|
|==================== | 20%
|
|======================================== | 40%
|
|============================================================ | 60%
|
|================================================================================ | 80%
|
|====================================================================================================| 100%
## INFO [2026-08-05 19:08:29] == Summarization is done.
head(summarized$FeatureLevelData)
## PROTEIN PEPTIDE TRANSITION FEATURE LABEL RUN GROUP SUBJECT FRACTION INTENSITY ABUNDANCE
## 1 555.447_9.1 5_NA NA_NA 5_NA_NA_NA L 1 Control 1 1 100 7.125593
## 2 Caffeic acid 4_NA NA_NA 4_NA_NA_NA L 1 Control 1 1 2000 11.447521
## 3 Caffeine 1_NA NA_NA 1_NA_NA_NA L 1 Control 1 1 1000 10.447521
## 4 Caffeine 6_NA NA_NA 6_NA_NA_NA L 1 Control 1 1 600 9.710556
## 5 GlucoseHigh 2_NA NA_NA 2_NA_NA_NA L 1 Control 1 1 5000 12.769449
## 6 Lactate 3_NA NA_NA 3_NA_NA_NA L 1 Control 1 1 800 10.125593
## originalRUN censored newABUNDANCE predicted
## 1 sampleAmzML TRUE NA NA
## 2 sampleAmzML FALSE 11.447521 NA
## 3 sampleAmzML FALSE 10.447521 NA
## 4 sampleAmzML FALSE 9.710556 NA
## 5 sampleAmzML FALSE 12.769449 NA
## 6 sampleAmzML FALSE 10.125593 NA
head(summarized$ProteinLevelData)
## RUN Protein LABEL LogIntensities originalRUN GROUP SUBJECT TotalGroupMeasurements
## 1 1 Caffeic acid L 11.447521 sampleAmzML Control 1 2
## 2 2 Caffeic acid L 10.753000 sampleBmzML Control 2 2
## 3 3 Caffeic acid L 10.949522 sampleCmzML Treatment 3 2
## 4 4 Caffeic acid L 10.943670 sampleDmzML Treatment 4 2
## 5 1 Caffeine L 10.079039 sampleAmzML Control 1 4
## 6 2 Caffeine L 9.440618 sampleBmzML Control 2 4
## NumMeasuredFeature MissingPercentage more50missing NumImputedFeature
## 1 1 0 FALSE 0
## 2 1 0 FALSE 0
## 3 1 0 FALSE 0
## 4 1 0 FALSE 0
## 5 2 0 FALSE 0
## 6 2 0 FALSE 0
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:
# 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
## Control Treatment
## Treatment vs Control -1 1
comparison = groupComparison(contrast.matrix = contrast_matrix,
data = summarized, use_log_file = FALSE)
## INFO [2026-08-05 19:08:29] == Start to test and get inference in whole plot ...
##
|
| | 0%
|
|========================= | 25%
|
|================================================== | 50%
|
|=========================================================================== | 75%
|
|====================================================================================================| 100%
## INFO [2026-08-05 19:08:30] == Comparisons for all proteins are done.
comparison$ComparisonResult
## Protein Label log2FC SE Tvalue DF pvalue adj.pvalue issue
## 1 Caffeic acid Treatment vs Control -0.15366483 0.34727290 -0.4424901 2 0.70138788 0.9068237 NA
## 2 Caffeine Treatment vs Control 0.09860066 0.32162317 0.3065720 2 0.78814165 0.9068237 NA
## 3 GlucoseHigh Treatment vs Control -0.05611639 0.42400990 -0.1323469 2 0.90682374 0.9068237 NA
## 4 Lactate Treatment vs Control -0.51047981 0.01148012 -44.4664311 1 0.01431445 0.0572578 NA
## MissingPercentage ImputationPercentage
## 1 0.00 0
## 2 0.00 0
## 3 0.00 0
## 4 0.25 0
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.
dataProcessPlots(summarized,
type = "ProfilePlot",
which.Protein = "Caffeine",
address = FALSE)
##
|
| | 0%
|
|====================================================================================================| 100%
##
##
|
| | 0%
##
|
|====================================================================================================| 100%
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.
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
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_GB
## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] data.table_1.18.4 MSstats_4.21.1 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.60 ggplot2_4.0.3 htmlwidgets_1.6.4
## [5] caTools_1.18.4 ggrepel_0.9.8 lattice_0.22-9 vctrs_0.7.3
## [9] tools_4.6.1 Rdpack_2.6.6 bitops_1.1-0 generics_0.1.4
## [13] parallel_4.6.1 tibble_3.3.1 pkgconfig_2.0.3 Matrix_1.7-6
## [17] KernSmooth_2.23-26 checkmate_2.3.4 RColorBrewer_1.1-3 S7_0.2.2
## [21] lifecycle_1.0.5 compiler_4.6.1 farver_2.1.2 gplots_3.3.0
## [25] statmod_1.5.2 litedown_0.10 htmltools_0.5.9 yaml_2.3.12
## [29] preprocessCore_1.75.0 marray_1.91.0 plotly_4.12.1 pillar_1.11.1
## [33] nloptr_2.2.1 tidyr_1.3.2 MASS_7.3-66 limma_3.69.2
## [37] reformulas_0.4.4 boot_1.3-32 nlme_3.1-170 commonmark_2.0.0
## [41] gtools_3.9.5 tidyselect_1.2.1 digest_0.6.39 stringi_1.8.9
## [45] dplyr_1.2.1 purrr_1.2.2 labeling_0.4.3 splines_4.6.1
## [49] fastmap_1.2.0 grid_4.6.1 cli_3.6.6 magrittr_2.0.5
## [53] dichromat_2.0-1 survival_3.8-9 withr_3.0.3 backports_1.5.1
## [57] scales_1.4.0 rmarkdown_2.31 httr_1.4.8 otel_0.2.0
## [61] lme4_2.0-6 evaluate_1.0.5 knitr_1.51 rbibutils_2.4.1
## [65] MSstatsConvert_1.23.3 viridisLite_0.4.3 markdown_2.0 rlang_1.3.0
## [69] Rcpp_1.1.2 glue_1.8.1 BiocManager_1.30.27 minqa_1.2.8
## [73] jsonlite_2.0.0 R6_2.6.1