Skip to contents

Exporting to standard and community formats

Source: vignettes/articles/tutorial-08-export-formats.Rmd
tutorial-08-export-formats.Rmd

Tutorial Prerequisites: Full workflow

Sharing results or running downstream statistics usually means handing the data to another tool. MRMhub offers two export routes: mzTab-M, the HUPO-PSI community standard that repositories expect, and a Bioconductor SummarizedExperiment for downstream analysis with limma (differential abundance) or lipidr (lipid-specific analysis).

Both examples start from a processed object. Here we build one from the bundled lipidomics_dataset; in practice this would be your own.

library(mrmhub)

mexp <- lipidomics_dataset |>
  normalize_by_istd() |>
  quantify_by_istd()

1. mzTab-M

mzTab-M is the HUPO-PSI community standard for reporting metabolomics and lipidomics quantification. It is a plain, tab-delimited text format that opens in Excel yet is fully machine-readable, and is what repositories such as MetaboLights expect on submission.

save_dataset_mztab() writes a processed experiment to an mzTab-M 2.0.0-M file:

out_dir <- tempdir()

save_dataset_mztab(mexp, file.path(out_dir, "experiment.mzTab"))

By default the final concentrations (feature_conc) are written as the per-sample abundances, with the concentration unit declared in the header. If the experiment has not been quantified, the exporter falls back to the raw feature_intensity and declares an “Arbitrary quantification unit”. Choose a different abundance variable with variable:

save_dataset_mztab(
  mexp, file.path(out_dir, "raw_areas.mzTab"),
  variable = "area")

The full dataset is exported, every analysis (including QC, blank and calibration samples) and every feature:

mzTab-M section mrmhub source
MTD metadata title, units, one ms_run/assay per analysis, one study_variable per qc_type
SMF (feature) one row per feature_id (quantifiers, qualifiers and ISTDs); abundance_assay[n] = chosen variable; ISTDs flagged via opt_global_is_internal_standard
SML (summary) one row per analyte, grouping its features; the quantifier drives the summary abundance and per-group mean / %CV
SME (evidence) a minimal identification stub per feature

Enrich the metadata header with optional arguments:

save_dataset_mztab(
  mexp, file.path(out_dir, "experiment.mzTab"),
  instrument = "Agilent 6495C QqQ",
  contact = "Jane Doe",
  publication = "doi:10.1234/example")

mzTab-M is a quantification report, not a full processing record: internal-standard relationships, QC and calibration metrics, and drift/batch state are not part of the model and are not reproduced on round-trip. Keep the MRMhubExperiment (or the Excel report from save_report_xlsx()) as the definitive record.

Validating the file

The output targets mzTab-M 2.0.0-M. To confirm conformance, upload the file to the HUPO-PSI / LIFS web validator at https://apps.lifs-tools.org/mztabvalidator/, or parse it back with the reference R package rmzTabM:

m <- rmzTabM::readMzTab(file.path(out_dir, "experiment.mzTab"))
rmzTabM::extractSmallMoleculeFeatures(m)

MRMhub has no runtime dependency on rmzTabM; the writer is self-contained.

Importing mzTab-M

import_data_mztab() ingests mzTab-M produced by other tools (for example Lipid Data Analyzer, MS-DIAL or MZmine) into an MRMhubExperiment:

mexp_in <- MRMhubExperiment(title = "Imported lipidomics")
mexp_in <- import_data_mztab(mexp_in, "LDA_export.mzTab")

Each Small Molecule Feature (SMF) becomes a feature and each assay an analysis. The per-assay abundances are imported as feature_intensity, and feature identities (name, formula, neutral mass, m/z, retention time) come from the SMF/SML sections. Where one analyte is reported as several features (e.g. different adducts), the adduct is appended to keep feature_id unique (Cer d18:1/16:0 | [M-H]-).

Import is partial by nature: mzTab-M carries a single abundance per feature, so internal-standard relationships, QC-type assignments, and calibration metadata are absent and must be supplied with add_metadata(). study_variable groups are imported best-effort as batch_id.

2. SummarizedExperiment

SummarizedExperiment is the Bioconductor container many statistical and lipidomics tools consume. save_dataset_summarizedexperiment() produces one, with features as rows and analyses as columns:

library(SummarizedExperiment)

se <- save_dataset_summarizedexperiment(mexp)
se
#> class: SummarizedExperiment
#> dim: 29 499
#> metadata(11): title analysis_type ... var_batch_corrected mrmhub_version
#> assays(9): rt area ... pmol_total conc
#> rownames(29): CE 18:1 CE 18:1 d7 (ISTD) ... TG 48:2 [-18:1] TG 48:2 [SIM]
#> rowData names(18): feature_id feature_class ... remarks feature_label
#> colnames(499): Longit_BLANK-01 (Eluent A) ... Longit_BLANK-07 (Eluent A)
#> colData names(13): analysis_order analysis_id ... annot_order_num remarks

Install the package once with BiocManager::install("SummarizedExperiment"). Pass a path to also write the object to disk as an .rds; without one it is only returned:

save_dataset_summarizedexperiment(mexp, "experiment.rds")

The mapping from experiment slots is:

Component mrmhub source
assays() one matrix per feature variable, named without the feature_ prefix
rowData() annot_features: one row per feature_id
colData() annot_analyses: one row per analysis_id
metadata() title, analysis type, processing status, is_* flags, concentration unit, MRMhub version

The parallel feature variables become parallel assays, so raw, normalized and quantified values live side-by-side in one object:

assayNames(se)
#> [1] "rt"     "area"   "height" "fwhm"   "width"  "intensity"
#> [7] "norm_intensity" "pmol_total" "conc"

assay(se, "conc")[1:3, 1:2]

Export a subset of variables with variable, or QC-filtered data by filtering first and passing filter_data = TRUE:

se_conc <- save_dataset_summarizedexperiment(mexp, variable = "conc")

mexp_filt <- mexp |>
  calc_qc_metrics() |>
  filter_features_qc(
    include_qualifier = FALSE,
    include_istd = FALSE,
    max.cv.conc.bqc = 25)

se_filt <- save_dataset_summarizedexperiment(mexp_filt, filter_data = TRUE)

Subsetting to study samples

Everything is exported and flagged rather than dropped, since downstream tools need it (lipidr, for example, requires the internal-standard annotation). But nothing downstream reads qc_type, so subset to study samples before a PCA, normalization, or differential test; otherwise blanks and QCs are folded in with the samples:

se_spl <- se[!rowData(se)$is_istd, se$qc_type == "SPL"]
dim(se_spl)
#> [1]  20 374

QC metrics are deliberately not written to rowData(): nothing downstream reads them, and QC filtering belongs in MRMhub where it is tested. Use filter_data = TRUE as above, or save_feature_qc_metrics() for the metrics.

Differential abundance with limma

limma works directly on continuous data such as concentrations. voom() is not used: it models count data and is invalid here.

The bundled lipidomics_dataset carries no phenotype: every study sample is plasma from the same longitudinal series. The grouping below is simulated to demonstrate the interface only; the p-values are meaningless by construction. With your own data the group would come from your sample metadata and already be in colData().

library(limma)

set.seed(1)
se_spl$group <- factor(
  sample(c("ctrl", "trt"), ncol(se_spl), replace = TRUE))

# concentrations are right-skewed; limma assumes normality
y <- log2(assay(se_spl, "conc"))

design <- model.matrix(~ 0 + group, data = colData(se_spl))
colnames(design) <- levels(se_spl$group)

fit <- lmFit(y, design)
fit <- contrasts.fit(fit, makeContrasts(trt - ctrl, levels = design))
fit <- eBayes(fit)

topTable(fit, number = 5)

Because colData() carries the full analysis annotation, covariates are available without extra joins, for example blocking on batch with ~ 0 + group + batch_id.

Lipid-specific analysis with lipidr

lipidr works on a LipidomicsExperiment, a subclass of SummarizedExperiment. Produce one directly:

le <- save_dataset_summarizedexperiment(
  mexp,
  variable = "intensity",
  as = "LipidomicsExperiment")

MRMhub fills in what lipidr requires: Molecule from feature_id, Class from feature_class, istd from is_istd, plus the summarized / logged / normalized flags lipidr reads but does not validate; an object missing them constructs cleanly and then misbehaves.

Use a peak-area scale variable with lipidr, not concentrations. lipidr log-transforms with log = TRUE by default and clamps values below 1 to 1 first, which assumes Skyline peak-area magnitudes. Concentrations in µmol/L are mostly below 1, so they are silently flattened to log2(1) = 0; eBayes() then reports zero residual variances and every fold change collapses. MRMhub warns when you export a mostly-sub-1 assay this way.

The trigger is absolute scale, not concentration as such: the clamp is a fixed threshold at 1. If you need concentrations, pass log = FALSE to lipidr, or log-transform them yourself beforehand.

library(lipidr)

le_spl <- le[!rowData(le)$istd, le$qc_type == "SPL"]
set.seed(1)
le_spl$group <- factor(
  sample(c("ctrl", "trt"), ncol(le_spl), replace = TRUE))

le_norm <- normalize_pqn(le_spl, measure = "intensity", log = TRUE)

de_results <- de_analysis(
  le_norm,
  trt - ctrl,
  measure = "intensity",
  group_col = "group")
plot_results_volcano(de_results, show.labels = FALSE)

Next steps