MSstats: Metabolomics workflow with MZMine
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.
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## 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## 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]+## 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.0The 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 22:01:55] ** Raw data from MZMine imported successfully.
## INFO [2026-08-05 22:01:55] ** MZMine ProteinName assignment: MZMine compound: 4 feature(s); SIRIUS name: 1 feature(s); m/z-RT fallback: 1 feature(s).
## INFO [2026-08-05 22:01:55] ** Raw data from MZMine cleaned successfully.
## INFO [2026-08-05 22:01:55] ** Using provided annotation.
## INFO [2026-08-05 22:01:55] ** Run labels were standardized to remove symbols such as '.' or '%'.
## INFO [2026-08-05 22:01:55] ** 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 22:01:55] ** Features with all missing measurements across runs are removed.
## INFO [2026-08-05 22:01:55] ** Multiple measurements in a feature and a run are summarized by summaryforMultipleRows: max
## INFO [2026-08-05 22:01:55] ** Features with all missing measurements across runs are removed.
## INFO [2026-08-05 22:01:55] ** Run annotation merged with quantification data.
## INFO [2026-08-05 22:01:55] ** Updated quantification data to make balanced design. Missing values are marked by NA
## INFO [2026-08-05 22:01:55] ** Finished preprocessing. The dataset is ready to be processed by the dataProcess function.## 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 4800ProteinName 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 22:01:55] ** Log2 intensities under cutoff = 7.93 were considered as censored missing values.
## INFO [2026-08-05 22:01:55] ** Log2 intensities = NA were considered as censored missing values.
## INFO [2026-08-05 22:01:55] ** Use all features that the dataset originally has.
## INFO [2026-08-05 22:01:55]
## # proteins: 5
## # peptides per protein: 1-2
## # features per peptide: 1-1
## INFO [2026-08-05 22:01:55] Some proteins have only one feature:
## 555.447_9.1,
## Caffeic acid,
## GlucoseHigh,
## Lactate ...
## INFO [2026-08-05 22:01:55]
## Control Treatment
## # runs 2 2
## # bioreplicates 2 2
## # tech. replicates 1 1
## INFO [2026-08-05 22:01:55] Some features are completely missing in at least one condition:
## 5_NA_NA_NA,
## NA ...
## INFO [2026-08-05 22:01:55] == Start the summarization per subplot...
## | | | 0% | |==================== | 20% | |======================================== | 40% | |============================================================ | 60% | |================================================================================ | 80% | |====================================================================================================| 100%
## INFO [2026-08-05 22:01:55] == Summarization is done.## 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## 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 0The 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 1comparison = groupComparison(contrast.matrix = contrast_matrix,
data = summarized, use_log_file = FALSE)## INFO [2026-08-05 22:01:56] == Start to test and get inference in whole plot ...
## | | | 0% | |========================= | 25% | |================================================== | 50% | |=========================================================================== | 75% | |====================================================================================================| 100%
## INFO [2026-08-05 22:01:56] == Comparisons for all proteins are done.## 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 0Each 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.
## | | | 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.
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
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 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: Etc/UTC
## 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 bslib_0.12.0 ggplot2_4.0.3
## [5] htmlwidgets_1.6.4 caTools_1.18.4 ggrepel_0.9.8 lattice_0.22-9
## [9] vctrs_0.7.3 tools_4.6.1 Rdpack_2.6.6 bitops_1.1-0
## [13] generics_0.1.4 parallel_4.6.1 tibble_3.3.1 pkgconfig_2.0.3
## [17] Matrix_1.7-6 KernSmooth_2.23-26 checkmate_2.3.4 RColorBrewer_1.1-3
## [21] S7_0.2.2 lifecycle_1.0.5 compiler_4.6.1 farver_2.1.2
## [25] gplots_3.3.0 statmod_1.5.2 htmltools_0.5.9 sys_3.4.3
## [29] buildtools_1.0.0 sass_0.4.10 yaml_2.3.12 preprocessCore_1.75.0
## [33] marray_1.91.0 plotly_4.12.1 pillar_1.11.1 nloptr_2.2.1
## [37] jquerylib_0.1.4 tidyr_1.3.2 MASS_7.3-66 cachem_1.1.0
## [41] limma_3.69.2 reformulas_0.4.4 boot_1.3-32 nlme_3.1-170
## [45] gtools_3.9.5 tidyselect_1.2.1 digest_0.6.39 stringi_1.8.9
## [49] dplyr_1.2.1 purrr_1.2.2 labeling_0.4.3 maketools_1.3.2
## [53] splines_4.6.1 fastmap_1.2.0 grid_4.6.1 cli_3.6.6
## [57] magrittr_2.0.5 survival_3.8-9 withr_3.0.3 backports_1.5.1
## [61] scales_1.4.0 rmarkdown_2.31 httr_1.4.8 otel_0.2.0
## [65] lme4_2.0-6 evaluate_1.0.5 knitr_1.51 rbibutils_2.4.1
## [69] MSstatsConvert_1.23.3 viridisLite_0.4.3 rlang_1.3.0 Rcpp_1.1.2
## [73] glue_1.8.1 BiocManager_1.30.27 minqa_1.2.8 jsonlite_2.0.0
## [77] R6_2.6.1