if (!require("remotes")) install.packages("remotes")
#remotes::install_github("SLINGhub/MRMhub", subdir = "quant")
#install.packages("here") # Used for parallel processingMRMhub Data Processing Workflow for Dataset 3
Supplementary Code 2
Overview Dataset
This example workflow uses a targeted lipidomics dataset measuring 829 features across 4,591 samples (Chen et al. (Chen et al. 2025)) and demonstrates the complete MRMhub data processing pipeline, including peak integration, quantitation, quality control, and reporting. This document covers many elements of the MRMhub workflow for demonstration purposes, some of which may not be applicable or necessary when processing other datasets.
Important
All datasets used and generated within this workflow have been deposited in the Zenodo repository https://zenodo.org/records/15370293. Within this repository, the raw mass spectrometry files and the MRMhub-INTEGRATOR input files required for peak integration are available under the record “Dataset 3”. The data and metadata files used and produced by the code in this Quarto notebook can be found under the”mrmhub-workflows” record. Please note that due to size limitations, the corresponding GitHub repository does not include any of these datasets.
Raw Data Processing: Peak Picking and Integration
This section describes key aspects of raw data processing workflow used for this dataset. For more a detailed description on the workflow, see Workflow for Dataset 1 and the INTEGRATOR Manual.
Conversion of vendor raw files to mzML
The original raw data files (Agilent .d) were converted to mzML using msconvert from the ProteoWizard software (https://proteowizard.sourceforge.io) with the following settings: output format = mzML; binary encoding precision = 32 bit; write index = true.
Preparation of the INTEGRATOR input files
The final input files used for this dataset are available from the Zenodo repository (see above). The following subsections summarize the key settings used for this example.
Global Settings File (param.txt)
The parameter mz_tol, defining the maximum tolerance in m/z values for identifying transitions, was set to 0.06. The RT_tol parameter, defining the window around the expected RT (see Feature/Transition table below) in which features are searched, was set to 0.04 min. This relatively narrow RT window was chosen to reduce selection of incorrect peaks in complex chromatograms, as observed in this dataset. The global peak_width parameter was set to [0.15, 0.08, 0.08, 0.18], which matched the majority of the peaks. For specific features with broader or convoluted peaks, feature-specific peak widths were applied, overwriting the global settings (see Feature Table below). The maximum allowed RT shift correction was set to (-0.1, 0.4), and the maximum allowed sample-to-sample RT shift was 0.05 min. The full set of settings in param.txt is shown in the figure below.
Sample Table
In the sample table (run_order_20251009.csv) the samples to be processed were listed by their .mzML file names. The first two BQC (batch quality control) samples were selected as retention time reference anchors. The sample type for each sample was defined in column C (sample_type) using the MRMhub QC type nomenclature. This is optional and INTEGRATOR distinguishes only between blanks (identified by the presence of the text BLK), and non-blanks.
Feature/Transition (Feature) Table
The feature/transition table (transition_list_20251010.csv) defined the transitions to be extracted from the .mzML files by precursor and product m/z values. Each feature (peak or peak group) was specified by its expected retention time (RT, in minutes). For some transitions, multiple features were defined. The corresponding internal standard feature for each entry was included in the table, however, these internal standards were not used by INTEGRATOR but were exported in the integration results.
Feature-specific ‘peak_width’ settings were applied for multiple features in this dataset when peaks were broader, exhibited increased tailing, or consisted of convoluted peaks that had to be co-integrated (see Figure 4 a-d). Fixed borders were assigned for features for which automatic integration failed or produced inconsistent results due to being poorly separated from adjacent peaks or having noisy, and consequently were applied when necessary (see Figure 4 e-h). The final version of the feature/transition table (transition_list_20251010.csv) is the result from several rounds of parameter optimization with INTEGRATOR.
Running MRMhub-INTEGRATOR and Review of Results
The most recent release OF INTEGRATOR was obtained from https://github.com/SLINGhub/MRMhub/releases. After preparing all input files, the mrmhub application is started. Steps 1 to 4 then run consecutively without further input to perform peak detection, peak picking, integration, and export of data and PDF results. The integration results are reported in the generated file ‘long.csv’, a long-format table that includes peak areas, actual retention time, peak width, and other metadata. This data file will be used for subsequent data postprocessing (see the following section). A wide-format table with peak areas is available as ‘quant_raw.csv’. PDFs of the integrated transitions are available in the by_* folders. See the INTEGRATOR Manual for more details.
The peak integration results were inspected for all features across samples using the generated PDFs (see above). In this dataset, RT shifts across the entire run of 4,591 samples were less than 0.04 min, except for ~ 20 samples at the start of a batch towards the end of this analysis, which had increased RT shifts.
Data Postprocessing and QC
Setting up a Quarto Notebook and Installation of {mrmhub}
A Quarto notebook was created in RStudio the {mrmhub} R package was installed by running following code in the R console, which also installs two other R packages used in this workflow.
The Quarto notebook code generated in this example workflow is available from https://github.com/SLINGhub/MRMhub-workflows (file Dataset3.qmd). In the subsequent section, the data postprocessing is encoded and documented step-by-step into the notebook as individual code chunks using {mrmhub} functions. Certain code chunks include additional R code to showcase specific examples for illustrative purposes. In real-world applications, this supplementary may not be necessary.
Postprocessing workflow
Load mrmhub and other required R packages
To improve performance multi-threading was used for some of the calculations and plottings, for which the R package {mirai} needs to be installed and loaded.
library(mrmhub)
library(mirai) # Used for parallel processingImport the MRMhub-INTEGRATOR results
The results from the peak integration performed with the INTEGRATOR workflow, as described in the previous section, were imported into a MRMhubExperiment data object. This object represents the central data container used in this postprocessing workflow. Nogte that the original result file from INTEGRATOR (‘long.csv’) has been renamed to ‘Dataset3_MRMhub-INTEGRATOR_20251010.csv’
data_path <- "./data/dataset-3/Dataset3_MRMhub-INTEGRATOR_20251010.csv"
mexp <- MRMhubExperiment()
mexp <- import_data_mrmhub(mexp, data_path, import_metadata = TRUE)
#> ✔ Imported 4591 analyses with 828 features
#> ℹ `feature_area` selected as default feature intensity. Modify with `set_intensity_var()`.
#> ✔ Analysis metadata associated with 4591 analyses.
#> ✔ Feature metadata associated with 828 features.Import the analysis metadata
Detailed analysis metadata, describing the analytes, samples, features, internal standards (ISTDs), and response curves required for subsequent steps, were imported from the Excel workbook ‘Dataset3_Metadata_20251010.xlsx’. This workbook is available within the previously mentioned Zenodo archive (record Dataset 3).
file_path <- "./data/dataset-3/Dataset3_Metadata_20251010.xlsx"
mexp <- import_metadata_msorganiser(mexp, file_path, ignore_warnings = TRUE)
#> ! Metadata has following warnings and notifications:
#> --------------------------------------------------------------------------------------------
#> Type Table Column Issue Count
#> 1 W* Features feature_id Feature(s) not in analysis data 31
#> 2 W* Features feature_id Feature(s) without metadata 1
#> --------------------------------------------------------------------------------------------
#> E = Error, W = Warning, W* = Supressed Warning, N = Note
#> --------------------------------------------------------------------------------------------
#> ✔ Analysis metadata associated with 4591 analyses.
#> ✔ Feature metadata associated with 827 features.
#> ✔ Internal Standard metadata associated with 26 ISTDs.
#> ✔ Response curve metadata associated with 36 annotated analyses.Analytical design and timeline
An overview of the analysis structure, detailing the types and sequence of QC samples analyzed is provided by the plot below. It also shows information on the start and end dates, total duration, and median run time for each sample. Setting show_timestamp = TRUE in the code below will plot the absolute dates as the x-axis, showing three short interruptions, and the last 4 batches were measured approximately one year later.
plot_runsequence(
mexp,
show_batches = TRUE,
qc_types = c("SPL", "BQC", "TQC", "PBLK", "UBLK", "RQC", "SBLK", "LTR", "NIST"),
batch_zebra_stripe = TRUE, base_font_size = 6,
batch_fill_color = "#fffbdb",
segment_linewidth = 0.25,
show_timestamp = FALSE) +
theme(plot.title = element_text(size = 5))
Overview Chromatographic Separation
The following plot shows the retention time distribution of all detected lipid species. The ceramides (Cer) species show little chromatographic separation in contrast to other lipid classes, with the majority eluting between 10.0 to 10.5 min. The triglycerides (TG) also elute within a narrower retention time range at around 11 min.
plot_abundanceprofile(
data = mexp,
log_scale = FALSE,
variable = "rt",
density_strip = TRUE,
qc_types = "SPL",
analysis_range = c(1,4000),
show_sum = FALSE,
#x_lim = c(6.5, 7.5),
x_label = NA,
feature_map = "lipidomics")
Overall trends and check for possible outliers
To examine overall technical trends and issues affecting most analytes (features), the Relative Log Abundance (RLA) plot is useful (De Livera et al., 2015 (Livera et al. 2015)). In an RLA plot, each feature is normalized to the across-sample or within-batch median and the resulting values are shown as boxplots for each sample. This visualization helps identify technical problems such as pipetting errors, sample loss, or changes in injection volume or instrument sensitivity. First, an overview of all analyses (runs) is plotted.
First, an overview of all analyses (samples) or samples is plotted, whereby on the signals from the spiked-in internal standards signals (ISTDs) are shown. The red horizontal lines indicate 1.5X IQR outlier limits.
fig2a <- plot_rla_boxplot(
data = mexp,
rla_type_batch = c("within"),
variable = "intensity",
qc_types = c("BQC", "TQC", "SPL", "NIST", "LTR", "SBLK"),
#plot_range = c(1, 4900),
rla_limit_to_range = FALSE,
filter_data = FALSE,
min_feature_intensity = 1000,
include_feature_filter = "ISTD",
#y_lim = c(-4,4),
show_timestamp = FALSE,
outlier_method = "fold",
outlier_k = c(-0.4, 0.3),
outlier_exclude = FALSE,
x_gridlines = FALSE,
show_plot = TRUE,
batch_zebra_stripe = TRUE,
linewidth = 0.1
)
#> ℹ Found 148 outliers in the 4545 shown analyses
Several outliers among study and QC samples were detected using this RLA analysis (visible also in the plot above)
# Get outlier in ISTD total signal excluding 3 batches with overall lower ISTD
istd_outlier <- fig2a$outliers |>
filter(!(batch_id %in% c("P-01", "P-02", "P-43") & val_res_median > -3))
print(istd_outlier)
#> # A tibble: 79 × 5
#> analysis_id analysis_order batch_id val_res_median qc_type
#> <chr> <int> <chr> <dbl> <chr>
#> 1 LV1017832619_MS-3 113 P-08 -11.8 SPL
#> 2 LV1017835020_MS-3 804 P-14 -11.8 SPL
#> 3 BQC-18_0_MS-4 1157 P-18 -2.69 BQC
#> 4 LV1017831110_MS-4 1185 P-18 -11.7 SPL
#> 5 LV1017832095_MS-4 1213 P-18 -12.1 SPL
#> 6 LV1017831999_MS-4 1328 P-19 -12.0 SPL
#> 7 LV1017834759_MS-4 1990 P-25 -11.9 SPL
#> 8 LV1017834663_MS-5 2218 P-27 -11.6 SPL
#> 9 TQC_217_MS-5 2320 P-28 -12.0 TQC
#> 10 TQC_218_MS-5 2334 P-28 -12.9 TQC
#> # ℹ 69 more rowsTo confirm these potential outliers, the peak areas of all ISTDs in different QC types is plotted against the run order (RunScatter plot). The outliers seen\ in these plots (with values close to zero) correspond to the outliers seen in the RLA plot above.
plot_runscatter(mexp, variable = "intensity",
include_qualifier = FALSE,
qc_types = c("SPL", "BQC", "TQC", "PBLK", "RQC", "SBLK"),
include_feature_filter = "IS",
#y_min = 0.00, y_max = 0.15,
#plot_range = c(0, 910),
point_size = .2,
point_border_width = 0.1,
point_transparency = .7,
base_font_size = 4,
cols_page = 3,
rows_page = 11,
show_progress = FALSE,
cap_outliers = TRUE)
#> Generating plots (1 page)