Dataset 4 — Steroid Assay
Dataset 4 is a fully quantitative steroid panel used to compare MRMhub with a vendor pipeline (Agilent MassHunter). Fifteen steroids were measured in human serum with the IBL International Steroid Panel LC-MS kit (cat. 30191875) (IBL International GmbH 2022) according to the instructions for use, with modifications. The kit acquires its analytes over two chromatographic panels; only Panel 1 (15 analytes) is used here, which is why every sample identifier carries the _P1 suffix (e.g. Cal A_P1, QC_High1_P1).
Unlike Datasets 1 and 3 (relative lipidomics), this is an externally calibrated absolute assay. The analytical design comprises an external calibration curve (Cal 0–Cal F), low and high quality-control samples (QC_Low1, QC_High1), and external quality-assessment (EQA) samples: SKML ring-trial materials with reference-method-assigned target concentrations (Jansen et al. 2014).
The workflow compares MRMhub and MassHunter at both pipeline stages:
The Dataset 4 data are small enough to be included in the GitHub repository under data/dataset-4/, so this notebook renders from a plain git clone with no further downloads.
The same files, together with those of the other datasets, are also deposited in the Zenodo record 10.5281/zenodo.15370293 (MRMhub-workflows) as mrmhub-workflows.zip. MRMhub-Dataset4.zip in that record holds the raw mzML files together with the MRMhub-INTEGRATOR application and all of its input files, so the peak integration compared below can be re-run as published or with new parameters. It is not required here: this notebook starts from INTEGRATOR’s long.csv output.
Raw .d files were converted to mzML and integrated with MRMhub-INTEGRATOR, as described for Datasets 1 and 3. The batch was integrated in two subsets, exported as separate long-format tables and combined downstream: the system-suitability injections (SST, ten replicate injections) and the assay samples (ASSAY: blanks, calibrators, QCs, EQA samples). INTEGRATOR reports peak areas, retention times and peak metadata in long.csv, and writes per-transition integration plots (the shaded, integrated region) to the by_transition folder.
Figure 1 shows six of these plots: three clean, well-resolved QC peaks, and beneath them three cases where the analyte sits next to a co-eluting neighbour, so that the integration boundaries are less obvious.
The same raw files were independently re-integrated in Agilent MassHunter Quantitative Analysis, providing the reference peak areas and concentrations used throughout this document.
QC_Low1): 11-deoxycortisol, Aldosterone, Dihydrotestosterone. Bottom row — three demanding cases, each with a co-eluting neighbour: 21-deoxycortisol in the high QC (QC_High1, large neighbour to the right), Aldosterone in EQA sample SKML2026_1Bb (neighbour to the right, peak at 0.3 × Cal A) and 11-deoxycorticosterone in SKML2026_1A (flanked on both sides; the largest MRMhub–MassHunter area difference in the batch).
The example chromatograms are cropped from the INTEGRATOR by_transition PDFs by scripts/dataset4-chromatograms/dataset4-chromatograms.R (run once; requires the magick and pdftools R packages) and embedded here as static images from images/dataset4/. The notebook render itself needs neither.
Unlike Datasets 1 and 3, this is a short pipeline: a 15-analyte, single-batch panel needs no drift/batch correction or feature filtering. The QUANT module (library(mrmhub)) imports the INTEGRATOR results and metadata, normalizes to the internal standards, fits external calibration curves, and quantifies.
This notebook was rendered with mrmhub 0.9.9 (QUANT module of MRMhub); install with pak::pak("SLINGhub/MRMhub"). The coloured console output additionally uses the fansi package (pak::pak("fansi")).
library(mrmhub)
library(mirai) # parallel processing (picked up automatically by mrmhub)
# Colour mrmhub's cli console feedback in the rendered HTML (requires the fansi package)
mrmhub_enable_cli_color()
# Set number of cores for parallel processing (all-but-one, so the system stays responsive)
n_cores <- {
n <- parallel::detectCores()
if (is.na(n)) n <- 4L
max(1L, n - 1L)
}
if (mirai::status()$daemons == 0) mirai::daemons(n_cores)The two INTEGRATOR subsets share columns and hold disjoint analyses, so they are combined by row-binding. The combined table is written to output/ because it is re-imported later (Figure 12) with MassHunter areas substituted in.
sst_long <- read_csv(
file.path(d4, "Dataset4_MRMhub-INTEGRATOR_SST.csv"),
show_col_types = FALSE
)
assay_long <- read_csv(
file.path(d4, "Dataset4_MRMhub-INTEGRATOR_ASSAY.csv"),
show_col_types = FALSE
)
long_combined <- bind_rows(sst_long, assay_long)
long_combined_path <- file.path("output", "Dataset4_long_combined.csv")
write_csv(long_combined, long_combined_path)
# Analysis ids per subset (import strips the ".mzML" extension).
sst_analysis_ids <- unique(str_remove(
sst_long$raw_data_filename,
"[.]mzML$"
))
assay_analysis_ids <- unique(str_remove(
assay_long$raw_data_filename,
"[.]mzML$"
))
meta_path <- file.path(d4, "Dataset4_Metadata.xlsx")
mexp <- MRMhubExperiment(title = "Dataset 4 — Steroid Assay")
mexp <- import_data_mrmhub(
mexp,
path = long_combined_path,
import_metadata = TRUE
)✔ Imported 30 analyses with 64 features.ℹ feature_area selected as default feature intensity. Modify with `set_intensity_var()`.✔ Analysis metadata associated with 30 analyses.✔ Feature metadata associated with 64 features.mexp <- import_metadata_msorganiser(
mexp,
path = meta_path,
excl_unmatched_analyses = TRUE, ignore_warnings = TRUE
)
#> Found no errors, 2 warnings, and 2 notes in the metadata.
#> ----------------------------------------------------------------
#> Type Table Column Issue Count
#> 1 W* Analyses analysis_id Analyses not in analysis data 2
#> 2 W* Features feature_id Feature(s) without metadata 34
#> 3 N Analyses sample_id Not defined for all analyses 16
#> 4 N Features analyte_id Not defined for all features 15
#>
#> ----------------------------------------------------------------
#> E = Error, W = Warning, W* = Suppressed Warning, N = Note
#> ----------------------------------------------------------------✔ Analysis metadata associated with 30 analyses.✔ Feature metadata associated with 30 features.✔ Internal Standard metadata associated with 15 ISTDs.✔ QC concentration metadata associated with 13 samples and 15 analytesEach feature is normalized to its stable-isotope-labelled internal standard, then a quadratic, 1/x-weighted external calibration curve is fitted per analyte from the Cal 0–Cal F levels.
mexp <- normalize_by_istd(mexp)✔ 15 features normalized with 15 ISTDs in 30 analyses.mexp <- calc_calibration_results(
mexp,
fit_overwrite = TRUE,
fit_model = "quadratic",
fit_weighting = "1/x"
)✔ Calibration curve fits calculated for all 15 quantifier features. Average r²: 0.9978.The fitted curves and their metrics are inspected before quantification.
plot_calibrationcurves(
mexp,
fit_overwrite = FALSE,
fit_model = "quadratic",
fit_weighting = "1/x",
include_istd = FALSE,
include_qualifier = FALSE,
show_progress = FALSE
)✔ Calibration curve fits calculated for all 15 quantifier features. Average r²: 0.9978.✔ Done
get_calibration_metrics(mexp, include_qualifier = FALSE, summary_table = TRUE)Concentrations are computed from the normalized intensities via the calibration curves. Qualifier transitions are excluded; analytes whose curve failed to fit are skipped rather than aborting the run.
mexp <- quantify_by_calibration(
mexp,
fit_overwrite = FALSE,
include_qualifier = FALSE,
ignore_failed_calibration = TRUE,
fit_model = "quadratic",
fit_weighting = "1/x"
)✔ Calibration curve fits calculated for all 15 quantifier features. Average r²: 0.9978.ℹ 77 concentration values fall outside the calibrated range (retained, flagged in feature_conc_out_of_range).✔ Concentrations calculated for 15 features in 30 analyses.✔ Concentrations are given in nmol/L.
# Export + reload the concentration table (wide: analysis_id, qc_type, analytes).
conc_path <- file.path("output", "Dataset4_quant_results.csv")
save_dataset_csv(
mexp,
path = conc_path,
variable = "conc",
add_qctype = TRUE
)✔ Concentration values for 30 analyses and 15 features have been exported to 'output/Dataset4_quant_results.csv'.d_conc_raw <- read_csv(conc_path, show_col_types = FALSE)The MassHunter Quantitative Analysis export is imported as a second experiment and reshaped to one tidy table — the single MassHunter source for all comparisons below. The trailing (ISTD) suffix on labelled internal standards is stripped so feature identifiers match MRMhub.
mh_reint <- import_data_masshunter(
MRMhubExperiment(title = "MassHunter re-integration"),
path = file.path(d4, "Dataset4_IntegrationResults_MassHunter.csv"),
import_metadata = TRUE, silent = TRUE
)ℹ feature_area selected as default feature intensity. Modify with `set_intensity_var()`.✔ Analysis metadata associated with 30 analyses.✔ Feature metadata associated with 30 features.mh_data <- get_analyticaldata(mh_reint, annotated = FALSE) |>
transmute(
analysis_id,
feature_id = str_remove(feature_id, "[ ]*\\(ISTD\\)$"),
rt_mh = feature_rt, area_mh = feature_area,
conc_mh = feature_conc_final, mi_mh = feature_manual_integration
)Study samples (SPL = the SKML EQA samples) whose quantifier peak area falls below the corresponding Cal A area are below the lowest non-zero calibrator and treated as not reliably quantifiable. The rule is applied per analyte and independently per method (MRMhub SPL area vs MRMhub Cal A area, MassHunter SPL area vs MassHunter Cal A area), yielding the flags excl_mrm / excl_mh that downstream comparisons apply to their SPL rows. Internal standards are never filtered.
CAL_A_ID <- "Cal A_P1"
# MRMhub quantifier areas (quantifiers only, excluding the ISTDs).
mrm_q_area <- get_analyticaldata(mexp, annotated = TRUE) |>
filter(is_quantifier, !is_istd) |>
transmute(
analysis_id,
analyte_id = feature_id, qc_type, area_mrm = feature_area
)
cala_mrm <- mrm_q_area |>
filter(analysis_id == CAL_A_ID) |>
select(analyte_id, cala_area_mrm = area_mrm)
cala_mh <- mh_data |>
filter(analysis_id == CAL_A_ID) |>
transmute(analyte_id = feature_id, cala_area_mh = area_mh)
# One row per SPL (analysis, analyte); flag below-Cal-A per method (NA-safe).
spl_area_excl <- mrm_q_area |>
filter(qc_type == "SPL") |>
left_join(
select(mh_data, analysis_id, analyte_id = feature_id, area_mh),
by = c("analysis_id", "analyte_id")
) |>
left_join(cala_mrm, by = "analyte_id") |>
left_join(cala_mh, by = "analyte_id") |>
mutate(
excl_mrm = !is.na(area_mrm) & !is.na(cala_area_mrm) &
area_mrm < cala_area_mrm,
excl_mh = !is.na(area_mh) & !is.na(cala_area_mh) &
area_mh < cala_area_mh
) |>
select(
analysis_id, analyte_id,
area_mrm, area_mh, cala_area_mrm, cala_area_mh,
excl_mrm, excl_mh
)Quantified concentrations of the low and high QC samples, one row per analyte.
d_conc_raw |>
filter(qc_type %in% c("LQC", "HQC")) |>
pivot_longer(
!c(analysis_id, qc_type),
names_to = "analyte", values_to = "conc"
) |>
mutate(qc = factor(qc_type, levels = c("LQC", "HQC"))) |>
select(analyte, qc, conc) |>
pivot_wider(names_from = qc, values_from = conc) |>
mutate(across(where(is.numeric), ~ signif(.x, 3)))Figure 7 compares MRMhub-INTEGRATOR peak areas with the MassHunter re-integration of the same raw files. MRMhub quantifier and ISTD peaks are joined to MassHunter on analysis_id + feature_id, dropping near-noise peaks (area < 100), empty CAL rows, and below-Cal-A SPL points.
# Manuscript house theme + palettes (shared by both figures).
theme_set(
theme_bw(base_size = 8, base_family = "sans") +
theme(
panel.grid.minor = element_blank(),
legend.position = "bottom",
plot.title = element_text(size = 8, face = "plain"),
legend.key.size = unit(3.2, "mm"), legend.margin = margin(1, 1, 1, 1)
)
)
col_method <- c(MRMhub = "#2C7FB8", MassHunter = "#D95F0E")
smp_cols <- c(
CAL = "#0072B2",
HQC = "#D55E00", LQC = "#E69F00", SPL = "grey30", SST = "#009E73",
SBLK = "#CC79A7", PBLK = "#999999"
)
# Manuscript composites are saved as vector PDF + 300-dpi PNG at 180 mm width
# with `save_plot()`, which selects the cairo and ragg devices (both render the
# Unicode in axis labels such as "µmol/L").mrm_peaks_q <- get_analyticaldata(mexp, annotated = TRUE) |>
filter(is_quantifier | is_istd) |>
transmute(
analysis_id,
feature_id, qc_type, sample_id,
feature_type = factor(
if_else(is_istd, "ISTD", "Quantifier"),
levels = c("Quantifier", "ISTD")
),
rt_mrm = feature_rt, area_mrm = feature_area
)
# Peak-level MRMhub↔MassHunter comparison table, built up by a filter cascade so
# only well-integrated, comparable peaks survive.
d_peakcmp_reint <- mrm_peaks_q |>
# Keep only peaks both integrators reported for the same analysis + feature.
inner_join(mh_data, by = c("analysis_id", "feature_id")) |>
filter(
!is.na(rt_mrm),
!is.na(rt_mh), !is.na(area_mrm), !is.na(area_mh)
) |>
# Drop near-noise peaks (either method below area 100).
filter(area_mrm >= 100, area_mh >= 100) |>
# Drop empty (unassigned) calibrator rows; ISTDs are always kept.
filter(!(qc_type == "CAL" & is.na(sample_id) & feature_type != "ISTD")) |>
left_join(
select(spl_area_excl, analysis_id, analyte_id, excl_mrm, excl_mh),
by = c("analysis_id", "feature_id" = "analyte_id")
) |>
# Drop below-Cal-A (below-LOQ) study samples flagged by either method.
filter(
!(qc_type == "SPL" &
(coalesce(excl_mrm, FALSE) | coalesce(excl_mh, FALSE)))
) |>
select(-excl_mrm, -excl_mh) |>
mutate(
area_diff_pct = (area_mh - area_mrm) / area_mrm * 100,
rt_diff = rt_mh - rt_mrm,
# Tag each analysis as an SST (system-suitability) or ASSAY run.
subset = case_when(
analysis_id %in% sst_analysis_ids ~ "SST",
analysis_id %in% assay_analysis_ids ~ "ASSAY", TRUE ~ NA_character_
)
)Log–log peak areas against the 1:1 line (Pearson r on log10 values). CAL/QC/SST samples span four orders of magnitude.
smp_lv <- c("CAL", "HQC", "LQC", "SST")
smp_cols4 <- c(
CAL = "#0072B2",
HQC = "#D55E00", LQC = "#E69F00", SST = "#009E73"
)
ft_shapes <- c(Quantifier = 16, ISTD = 17)
scale_smp <- scale_colour_manual(
values = smp_cols4,
limits = smp_lv, drop = FALSE
)
scale_ft <- scale_shape_manual(values = ft_shapes)
d_cmp <- d_peakcmp_reint |>
filter(qc_type %in% smp_lv) |>
mutate(qc_type = factor(qc_type, levels = smp_lv))
r_log <- cor(log10(d_cmp$area_mrm), log10(d_cmp$area_mh))
ax_lims <- range(c(d_cmp$area_mrm, d_cmp$area_mh))
p_corr <- d_cmp |>
ggplot(aes(area_mrm, area_mh)) +
geom_abline(
slope = 1,
intercept = 0, linetype = "dashed", colour = "grey40"
) +
geom_point(
aes(colour = qc_type, shape = feature_type),
size = 1.1, alpha = 0.7
) +
# Ring the peaks MassHunter integrated manually (MI flag), keeping the
# sample-type colour / feature shape encoding underneath.
geom_point(
data = ~ filter(.x, mi_mh),
aes(alpha = "Manual integration (MassHunter)"), shape = 1,
colour = "grey15", size = 2.3, stroke = 0.5
) +
annotate(
"text",
x = ax_lims[1], y = ax_lims[2], hjust = 0, vjust = 1, size = 2.4,
label = sprintf("Pearson r = %.4f", r_log)
) +
scale_x_log10(limits = ax_lims) +
scale_y_log10(limits = ax_lims) +
scale_smp +
scale_ft +
scale_alpha_manual(
name = NULL,
values = c("Manual integration (MassHunter)" = 1)
) +
guides(
colour = guide_legend(order = 1, nrow = 2),
shape = guide_legend(order = 2, nrow = 2),
alpha = guide_legend(
order = 3,
override.aes = list(
shape = 1,
colour = "grey15", size = 2.3, stroke = 0.5
)
)
) +
labs(
x = "MRMhub peak area",
y = "MassHunter peak area", colour = "Sample type", shape = "Feature",
title = "Peak-area agreement"
)
p_corr
Per-feature %CV over the ten replicate SST injections, MRMhub vs MassHunter. The two methods show comparable precision, with points scattering around the 1:1 line.
sst_cv <- d_peakcmp_reint |>
filter(subset == "SST") |>
group_by(feature_type, feature_id) |>
summarise(
cv_mrmhub = sd(area_mrm) / mean(area_mrm) * 100,
cv_masshunter = sd(area_mh) / mean(area_mh) * 100, .groups = "drop"
)
cv_lims <- range(c(sst_cv$cv_mrmhub, sst_cv$cv_masshunter)) + c(-0.3, 0.3)
p_cv <- sst_cv |>
ggplot(aes(cv_mrmhub, cv_masshunter, shape = feature_type)) +
geom_abline(
slope = 1,
intercept = 0, linetype = "dashed", colour = "grey60"
) +
geom_point(size = 1.6, alpha = 0.85, colour = "#009E73") +
coord_cartesian(xlim = cv_lims, ylim = cv_lims) +
scale_ft +
guides(shape = "none") +
labs(
x = "MRMhub %CV (10 SST)",
y = "MassHunter %CV", title = "Precision (%CV) equivalence"
)
p_cv
Percent difference (MassHunter − MRMhub) against the mean area, excluding SST and near-noise peaks (mean area < 10³). The three highlighted analytes are among those with the largest real-sample differences (see Table 3); their chromatograms are shown below.
highlight_analytes <- c(
"21-deoxycortisol",
"11-deoxycortisol", "Aldosterone"
)
hl_cols <- setNames(c("#D81B60", "#1E88E5", "#8E44AD"), highlight_analytes)
d_ba <- d_peakcmp_reint |>
filter(qc_type != "SST") |>
mutate(
mean_area = (area_mrm + area_mh) / 2,
pct_diff = (area_mh - area_mrm) / mean_area * 100,
hl = factor(
if_else(
feature_id %in% highlight_analytes,
feature_id, NA_character_
),
levels = highlight_analytes
)
) |>
filter(mean_area >= 1e3)
ba_bias <- mean(d_ba$pct_diff)
ba_loa <- ba_bias + c(-1.96, 1.96) * sd(d_ba$pct_diff)
ba_xr <- max(d_ba$mean_area)
p_ba <- d_ba |>
ggplot(aes(mean_area, pct_diff)) +
geom_hline(yintercept = 0, colour = "grey70") +
geom_hline(yintercept = ba_bias, colour = "grey20") +
geom_hline(yintercept = ba_loa, linetype = "dashed", colour = "grey50") +
geom_point(
data = ~ filter(.x, is.na(hl)),
colour = "grey80", size = 0.9, alpha = 0.6
) +
geom_point(
data = ~ filter(.x, !is.na(hl)),
aes(colour = hl), size = 1.7, alpha = 0.9
) +
# Ring the peaks MassHunter integrated manually (MI flag).
geom_point(
data = ~ filter(.x, mi_mh),
aes(alpha = "Manual integration (MassHunter)"), shape = 1,
colour = "grey15", size = 2.0, stroke = 0.5
) +
annotate(
"text",
x = ba_xr, y = ba_bias, hjust = 1, vjust = -0.6, size = 2.1,
colour = "grey20", label = sprintf("mean bias %.1f%%", ba_bias)
) +
annotate(
"text",
x = ba_xr, y = ba_loa[2], hjust = 1, vjust = 1.4, size = 2,
colour = "grey45", label = "+95% LoA"
) +
annotate(
"text",
x = ba_xr, y = ba_loa[1], hjust = 1, vjust = -0.7, size = 2,
colour = "grey45", label = "-95% LoA"
) +
scale_colour_manual(
values = hl_cols,
name = "Highlighted analyte", na.translate = FALSE
) +
scale_alpha_manual(
name = NULL,
values = c("Manual integration (MassHunter)" = 1)
) +
scale_x_log10(limits = c(1e3, NA)) +
coord_cartesian(ylim = c(-20, 20)) +
guides(
alpha = guide_legend(
order = 2,
override.aes = list(
shape = 1,
colour = "grey15", size = 2.0, stroke = 0.5
)
)
) +
labs(
x = "Mean peak area of the two methods (>= 10^3)",
y = "Percent difference (MH - MRMhub) %",
title = sprintf(
"Bland-Altman (excl. SST, area >= 10^3): mean bias %.1f%%, 95%% LoA [%.0f, %.0f]%%",
ba_bias, ba_loa[1], ba_loa[2]
)
)
p_ba
The largest real-sample area differences occur where a peak sits on an elevated baseline or beside a co-eluting neighbour, where the two integrators place the integration boundaries differently. The candidate table ranks these cases; the panel crops the three chosen examples (one per highlighted analyte) from the MRMhub by_transition chromatograms.
chrom_candidates <- d_peakcmp_reint |>
filter(
feature_type == "Quantifier",
qc_type %in% c("CAL", "HQC", "LQC", "SPL", "SST"), area_mrm >= 20000
) |>
arrange(desc(abs(area_diff_pct))) |>
transmute(
feature_id,
analysis_id, qc_type,
area_mrm = round(area_mrm),
area_mh = round(area_mh), area_diff_pct = round(area_diff_pct, 1)
) |>
head(10)
chrom_candidates# Panel-D chromatograms are the pre-generated crops (see the INTEGRATOR section);
# subtitles carry the data-derived area difference for each example.
chrom_examples <- tibble::tribble(
~name,
~feature_id, ~analysis_id, "fig1d_21-deoxycortisol", "21-deoxycortisol",
"SKML2025_6B_P1", "fig1d_11-deoxycortisol", "11-deoxycortisol",
"SKML2026_2B_P1", "fig1d_Aldosterone", "Aldosterone", "SKML2026_2A_P1"
)
chrom_meta <- chrom_examples |>
left_join(
distinct(d_peakcmp_reint, feature_id, analysis_id, area_diff_pct),
by = c("feature_id", "analysis_id")
) |>
mutate(
subtitle = sprintf(
"%s (dArea %+.1f%%)",
str_remove(analysis_id, "_P1$"), area_diff_pct
)
)
chrom_plots <- Map(
chrom_gg,
chrom_meta$name, chrom_meta$feature_id, chrom_meta$subtitle
)
p_chrom <- wrap_elements(full = wrap_plots(chrom_plots, nrow = 1))
p_chrom
The four panels above are composed into the publication figure (Figure 1) and saved as a vector PDF and a 300-dpi PNG (180 mm double-column width).
row_ab <- (p_corr + p_cv) +
plot_layout(guides = "collect") &
theme(legend.position = "bottom")
fig1 <- row_ab /
p_ba /
p_chrom +
plot_layout(heights = c(1, 0.72, 0.58)) +
plot_annotation(tag_levels = "A")
save_plot(
fig1, "output/fig1_integration_correspondence", format = c("pdf", "png"),
width = 180, height = 230)✔ Saved plot (180 x 230 mm) to output/fig1_integration_correspondence.pdf, output/fig1_integration_correspondence.png
output/fig1_integration_correspondence.{pdf,png} (180 mm).
Figure 12 compares the two pipelines at the quantification stage. To separate the post-processing (calibration + quantification) from the peak integration, MassHunter’s own peak areas are also run through the MRMhub pipeline, with identical areas and an identical calibration model, so that any concentration difference arises from the post-processing alone.
# Rebuild the INTEGRATOR long table with MassHunter's peak areas swapped in for
# MRMhub's (matched on analysis + feature). Running the identical pipeline on this
# table — same areas, same calibration model — makes any concentration difference
# vs MassHunter attributable to post-processing alone, not to peak integration.
mh_area <- mh_data |>
transmute(analysis_id, feature_name = feature_id, area_mh)
long_mh <- read_csv(long_combined_path, show_col_types = FALSE) |>
mutate(.aid = str_remove(raw_data_filename, "[.]mzML$")) |>
inner_join(mh_area, by = c(".aid" = "analysis_id", "feature_name")) |>
mutate(area = area_mh) |>
select(-.aid, -area_mh)
# The pipeline reads from disk, so write the area-substituted table to a temp file.
tmp_long <- tempfile(fileext = ".csv")
write_csv(long_mh, tmp_long)
mexp_mh <- MRMhubExperiment(title = "MassHunter areas via MRMhub pipeline")
mexp_mh <- import_data_mrmhub(
mexp_mh,
path = tmp_long,
import_metadata = TRUE
)✔ Imported 30 analyses with 30 features.ℹ feature_area selected as default feature intensity. Modify with `set_intensity_var()`.✔ Analysis metadata associated with 30 analyses.✔ Feature metadata associated with 30 features.mexp_mh <- import_metadata_msorganiser(
mexp_mh,
path = meta_path,
excl_unmatched_analyses = TRUE, ignore_warnings = TRUE
)
#> Found no errors, 1 warning, and 2 notes in the metadata.
#> ----------------------------------------------------------------
#> Type Table Column Issue Count
#> 1 W* Analyses analysis_id Analyses not in analysis data 2
#> 2 N Analyses sample_id Not defined for all analyses 16
#> 3 N Features analyte_id Not defined for all features 15
#>
#> ----------------------------------------------------------------
#> E = Error, W = Warning, W* = Suppressed Warning, N = Note
#> ----------------------------------------------------------------✔ Analysis metadata associated with 30 analyses.
✔ Feature metadata associated with 30 features.✔ Internal Standard metadata associated with 15 ISTDs.✔ QC concentration metadata associated with 13 samples and 15 analytesmexp_mh <- normalize_by_istd(mexp_mh)
#> Warning: [1m[22m[33m
The same comparison as a table, with three fits per analyte (quadratic, y = a·x² + b·x + c, 1/x weighting): MassHunter’s own fit, the MRMhub end-to-end fit (MRMhub-INTEGRATOR areas), and the MRMhub pipeline run on MassHunter’s own peak areas. The gap between the two MRMhub rows is the peak-integration effect; the distance between MRMhub (on MH areas) and MassHunter isolates the post-processing (normalization, weighting, regression) alone.
library(gt)
# Measured LQC / HQC concentrations (nmol/L) per engine, one column each.
qc_wide <- function(df) {
df |>
filter(qc_type %in% c("LQC", "HQC")) |>
pivot_longer(
!c(analysis_id, qc_type),
names_to = "analyte_id", values_to = "conc"
) |>
select(analyte_id, qc_type, conc) |>
pivot_wider(names_from = qc_type, values_from = conc)
}
tmp_mhconc <- tempfile(fileext = ".csv")
save_dataset_csv(
mexp_mh,
path = tmp_mhconc,
variable = "conc",
add_qctype = TRUE
)✔ Concentration values for 30 analyses and 15 features have been exported to '/var/folders/3r/ywcsb3896zj4_0xlb_70yvlh0000gn/T//Rtmpe7n07o/file702e6e482fd4.csv'.conc_e2e <- qc_wide(d_conc_raw)
conc_onmh <- qc_wide(read_csv(tmp_mhconc, show_col_types = FALSE))
conc_mhun <- mh_data |>
filter(analysis_id %in% c("QC_Low1_P1", "QC_High1_P1")) |>
transmute(
analyte_id = feature_id,
qc_type = if_else(analysis_id == "QC_Low1_P1", "LQC", "HQC"),
conc = conc_mh
) |>
pivot_wider(names_from = qc_type, values_from = conc)
# Assemble one row per analyte for each of the three fits. mrmhub reports
# coefficients intercept-first (coef_a = x^0 … coef_c = x^2), so the two MRMhub
# blocks remap to MassHunter's descending order (a = x^2, b = x, c = intercept);
# MassHunter's export is already in that order.
mh_own <- read_csv(
file.path(d4, "Dataset4_MassHunter-CalibrationFits.csv"),
show_col_types = FALSE
) |>
transmute(
analyte_id,
method = "MassHunter", a = coef_a, b = coef_b, c = coef_c, r2 = r2
) |>
left_join(conc_mhun, by = "analyte_id")
# MRMhub pipeline fit on MassHunter's peak areas (isolates post-processing).
mrm_mh <- get_calibration_metrics(mexp_mh, include_qualifier = FALSE) |>
transmute(
analyte_id = feature_id,
method = "MRMhub (on MH areas)", a = coef_c, b = coef_b, c = coef_a,
r2 = r2
) |>
left_join(conc_onmh, by = "analyte_id")
# MRMhub end-to-end fit on MRMhub-INTEGRATOR areas.
mrm_e2e <- get_calibration_metrics(mexp, include_qualifier = FALSE) |>
transmute(
analyte_id = feature_id,
method = "MRMhub (end-to-end)", a = coef_c, b = coef_b, c = coef_a,
r2 = r2
) |>
left_join(conc_e2e, by = "analyte_id")
cal_lv <- c("MassHunter", "MRMhub (on MH areas)", "MRMhub (end-to-end)")
fmt_conc <- function(x) {
ifelse(
is.na(x),
"—",
formatC(
signif(x, 4),
format = "fg", big.mark = ",", drop0trailing = TRUE
)
)
}
cal_coef_tbl <- bind_rows(mh_own, mrm_mh, mrm_e2e) |>
mutate(method = factor(method, levels = cal_lv)) |>
arrange(analyte_id, method) |>
# Stripe by analyte group (3 rows each), not by row.
mutate(.grp = match(analyte_id, unique(analyte_id))) |>
transmute(
.grp,
Analyte = if_else(as.character(method) == cal_lv[1], analyte_id, ""),
Method = as.character(method), `a (x²)` = sprintf("%.4g", a),
`b (x)` = sprintf("%.4g", b), `c` = sprintf("%.4g", c),
`R²` = sprintf("%.5f", r2), `LQC Concentration` = fmt_conc(LQC),
`HQC Concentration` = fmt_conc(HQC)
)
gt(cal_coef_tbl) |>
cols_hide(".grp") |>
cols_align(
"right",
columns = c(
"a (x²)",
"b (x)", "c", "R²", "LQC Concentration", "HQC Concentration"
)
) |>
# Zebra shading applied per analyte, so each analyte's three method rows
# share one stripe and alternate against the next analyte. Both parities get
# an explicit inline fill (odd = white) so Quarto's per-row table striping
# can't leak through the unfilled groups.
tab_style(
cell_fill(color = "#FFFFFF"),
cells_body(rows = .grp %% 2 == 1)
) |>
tab_style(
cell_fill(color = "#F2F2F2"),
cells_body(rows = .grp %% 2 == 0)
) |>
# MRMhub end-to-end fit distinguished by blue text (see caption).
tab_style(
cell_text(color = "#2C7FB8"),
cells_body(rows = Method == "MRMhub (end-to-end)")
) |>
tab_options(table.font.size = px(12), data_row.padding = px(2))| Analyte | Method | a (x²) | b (x) | c | R² | LQC Concentration | HQC Concentration |
|---|---|---|---|---|---|---|---|
| 11-deoxycorticosterone | MassHunter | 0.01256 | 1.471 | 0.04589 | 0.99815 | 0.6529 | 10.32 |
| MRMhub (on MH areas) | 0.01256 | 1.471 | 0.0459 | 0.99815 | 0.6529 | 10.32 | |
| MRMhub (end-to-end) | 0.01313 | 1.427 | 0.04811 | 0.99791 | 0.6723 | 10.19 | |
| 11-deoxycortisol | MassHunter | 0.005944 | 0.9796 | 0.06228 | 0.99929 | 1.458 | 22.95 |
| MRMhub (on MH areas) | 0.005944 | 0.9796 | 0.06228 | 0.99929 | 1.458 | 22.95 | |
| MRMhub (end-to-end) | 0.005598 | 1.022 | 0.0644 | 0.99970 | 1.451 | 23.48 | |
| 17-hydroxyprogesterone | MassHunter | 0.001839 | 0.9685 | 0.03238 | 0.99971 | 1.488 | 26.33 |
| MRMhub (on MH areas) | 0.001839 | 0.9685 | 0.03238 | 0.99971 | 1.488 | 26.33 | |
| MRMhub (end-to-end) | 0.001363 | 0.9671 | 0.03846 | 0.99982 | 1.513 | 27.05 | |
| 21-deoxycortisol | MassHunter | 0.006042 | 0.7758 | 0.008574 | 0.99935 | 1.698 | 22.46 |
| MRMhub (on MH areas) | 0.006042 | 0.7758 | 0.008574 | 0.99935 | 1.697 | 22.46 | |
| MRMhub (end-to-end) | 0.005988 | 0.7606 | 0.002312 | 0.99938 | 1.684 | 22.37 | |
| Aldosterone | MassHunter | 0.01507 | 0.5417 | 0.1124 | 0.99656 | 0.9199 | 9.798 |
| MRMhub (on MH areas) | 0.01507 | 0.5417 | 0.1124 | 0.99656 | 0.9199 | 9.798 | |
| MRMhub (end-to-end) | 0.01406 | 0.518 | 0.08996 | 0.99623 | 0.9878 | 9.864 | |
| Androstenedione | MassHunter | 0.005103 | 0.6132 | 0.04322 | 0.99955 | 1.81 | 28.02 |
| MRMhub (on MH areas) | 0.005103 | 0.6132 | 0.04322 | 0.99955 | 1.81 | 28.02 | |
| MRMhub (end-to-end) | 0.004824 | 0.6072 | 0.0427 | 0.99945 | 1.842 | 27.68 | |
| Corticosterone | MassHunter | -6.449e-05 | 0.4008 | -0.01636 | 0.99695 | 3.758 | 81.6 |
| MRMhub (on MH areas) | -6.45e-05 | 0.4008 | -0.01636 | 0.99695 | 3.757 | 81.6 | |
| MRMhub (end-to-end) | -6.664e-05 | 0.4109 | -0.01616 | 0.99722 | 3.601 | 82.82 | |
| Cortisol | MassHunter | 1.227e-05 | 0.1237 | -0.06246 | 0.99942 | 24.37 | 504.8 |
| MRMhub (on MH areas) | 1.227e-05 | 0.1237 | -0.06246 | 0.99942 | 24.37 | 504.8 | |
| MRMhub (end-to-end) | 1.526e-05 | 0.1212 | -0.03548 | 0.99964 | 24.89 | 504.6 | |
| Cortisone | MassHunter | -0.0003965 | 0.4315 | -0.1158 | 0.99967 | 6.238 | 115.2 |
| MRMhub (on MH areas) | -0.0003965 | 0.4315 | -0.1158 | 0.99967 | 6.238 | 115.2 | |
| MRMhub (end-to-end) | -0.0003793 | 0.428 | -0.09925 | 0.99970 | 6.269 | 115.6 | |
| DHEA | MassHunter | -3.25e-05 | 0.03091 | 0.1772 | 0.98415 | 8.355 | 129.2 |
| MRMhub (on MH areas) | -3.25e-05 | 0.03091 | 0.1772 | 0.98415 | 8.355 | 129.2 | |
| MRMhub (end-to-end) | -4.076e-05 | 0.03256 | 0.1832 | 0.98458 | 8.937 | 131.9 | |
| DHEAS | MassHunter | -7.313e-09 | 0.003939 | -0.0933 | 0.99903 | 524.2 | 9,931 |
| MRMhub (on MH areas) | -7.313e-09 | 0.003939 | -0.0933 | 0.99903 | 524.2 | 9,931 | |
| MRMhub (end-to-end) | -1.016e-08 | 0.003949 | -0.09189 | 0.99918 | 521.2 | 10,090 | |
| Dexamethasone | MassHunter | -1.546e-05 | 0.1036 | -0.01331 | 0.99797 | 6.046 | 110.4 |
| MRMhub (on MH areas) | -1.546e-05 | 0.1036 | -0.01331 | 0.99797 | 6.046 | 110.4 | |
| MRMhub (end-to-end) | -1.714e-05 | 0.1036 | -0.007976 | 0.99826 | 5.983 | 111.7 | |
| Dihydrotestosterone | MassHunter | -0.03613 | 0.5642 | -0.1789 | 0.97656 | 1.224 | 6.929 |
| MRMhub (on MH areas) | -0.03613 | 0.5642 | -0.1789 | 0.97656 | 1.224 | 6.929 | |
| MRMhub (end-to-end) | -0.01988 | 0.4645 | -0.07173 | 0.99770 | 1.212 | 5.861 | |
| Progesterone | MassHunter | -0.002326 | 0.9786 | 0.03361 | 0.99961 | 1.349 | 31.82 |
| MRMhub (on MH areas) | -0.002326 | 0.9786 | 0.03361 | 0.99961 | 1.349 | 31.82 | |
| MRMhub (end-to-end) | -0.001442 | 0.9495 | 0.05202 | 0.99967 | 1.322 | 31.06 | |
| Testosterone | MassHunter | 0.003025 | 1.22 | 0.07571 | 0.99918 | 0.87 | 31.79 |
| MRMhub (on MH areas) | 0.003025 | 1.22 | 0.07571 | 0.99918 | 0.87 | 31.79 | |
| MRMhub (end-to-end) | 0.002567 | 1.216 | 0.0762 | 0.99922 | 0.8558 | 32.85 |
MassHunter areas run through the MRMhub pipeline give the same concentrations as MassHunter end-to-end (Pearson r on log10 values below), so the post-processing steps of the two pipelines give the same result on identical input.
tmp_conc <- tempfile(fileext = ".csv")
save_dataset_csv(
mexp_mh,
path = tmp_conc,
variable = "conc",
add_qctype = TRUE
)✔ Concentration values for 30 analyses and 15 features have been exported to '/var/folders/3r/ywcsb3896zj4_0xlb_70yvlh0000gn/T//Rtmpe7n07o/file702e7f335eb3.csv'.d_qonly <- read_csv(tmp_conc, show_col_types = FALSE) |>
pivot_longer(
!c(analysis_id, qc_type),
names_to = "analyte_id", values_to = "conc_mrmpipe"
) |>
left_join(
mh_data |> select(analysis_id, analyte_id = feature_id, conc_mh),
by = c("analysis_id", "analyte_id")
) |>
filter(!is.na(conc_mrmpipe), !is.na(conc_mh)) |>
filter(!qc_type %in% c("SBLK", "PBLK")) |>
mutate(qc_type = factor(qc_type, levels = names(smp_cols)))
posq <- d_qonly |> filter(conc_mrmpipe > 0, conc_mh > 0)
q_r <- cor(log10(posq$conc_mrmpipe), log10(posq$conc_mh))
axq <- range(c(posq$conc_mrmpipe, posq$conc_mh))
p_qonly <- posq |>
ggplot(aes(conc_mh, conc_mrmpipe, colour = qc_type)) +
geom_abline(
slope = 1,
intercept = 0, linetype = "dashed", colour = "grey40"
) +
geom_point(size = 1.1, alpha = 0.6) +
annotate(
"text",
x = axq[1], y = axq[2], hjust = 0, vjust = 1, size = 2.4,
label = sprintf("Pearson r = %.5f", q_r)
) +
scale_x_log10(limits = axq) +
scale_y_log10(limits = axq) +
scale_colour_manual(
values = smp_cols,
drop = TRUE, name = "Sample type"
) +
coord_equal() +
labs(
x = "MassHunter end-to-end (nmol/L)",
y = "MRMhub-QUANT on MH areas (nmol/L)",
title = "Concentration equivalence (same areas)"
) +
theme(
legend.position = "inside",
legend.position.inside = c(0.98, 0.02),
legend.justification.inside = c(1, 0),
legend.background = element_rect(
fill = alpha("white", 0.7),
colour = NA
),
legend.key.size = unit(2.8, "mm"),
legend.title = element_text(size = 6.5),
legend.text = element_text(size = 6)
)
p_qonly
Independent end-to-end concentrations (each method’s own integration and calibration) for the LQC, HQC and SPL samples.
d_full <- d_conc_raw |>
filter(qc_type %in% c("LQC", "HQC", "SPL")) |>
pivot_longer(
!c("analysis_id", "qc_type"),
names_to = "analyte_id", values_to = "conc"
) |>
left_join(
mh_data |> select(analysis_id, analyte_id = feature_id, conc_mh),
by = c("analysis_id", "analyte_id")
) |>
left_join(
select(spl_area_excl, analysis_id, analyte_id, excl_mrm, excl_mh),
by = c("analysis_id", "analyte_id")
) |>
mutate(
conc = if_else(coalesce(excl_mrm, FALSE), NA_real_, conc),
conc_mh = if_else(coalesce(excl_mh, FALSE), NA_real_, conc_mh)
) |>
select(-excl_mrm, -excl_mh) |>
filter(!is.na(conc), !is.na(conc_mh), conc > 0, conc_mh > 0) |>
mutate(qc_type = factor(qc_type, levels = names(smp_cols)))
f_r <- cor(log10(d_full$conc), log10(d_full$conc_mh))
axf <- range(c(d_full$conc, d_full$conc_mh))
p_full <- d_full |>
ggplot(aes(conc_mh, conc, colour = qc_type)) +
geom_abline(
slope = 1,
intercept = 0, linetype = "dashed", colour = "grey40"
) +
geom_point(size = 1.3, alpha = 0.7) +
annotate(
"text",
x = axf[1], y = axf[2], hjust = 0, vjust = 1, size = 2.4,
label = sprintf("Pearson r = %.4f", f_r)
) +
scale_x_log10(limits = axf) +
scale_y_log10(limits = axf) +
scale_colour_manual(
values = smp_cols,
drop = TRUE, name = "Sample type"
) +
coord_equal() +
labs(
x = "MassHunter end-to-end (nmol/L)",
y = "MRMhub end-to-end (nmol/L)",
title = "Concentration agreement (full pipeline)"
) +
theme(legend.position = "none")
p_full
Per-analyte percent bias against the reference concentrations (the SKML EQA reference values and the nominal LQC/HQC target concentrations), one point per sample per method (dashed lines: ±15 %; LQC and HQC are grouped to the left of the EQA samples). EQA points whose quantifier area falls below the lowest calibrator (Cal A, i.e. below the limit of quantification) are drawn as open symbols rather than excluded. The reference concentrations are read from the metadata workbook’s QC Concentrations sheet (imported above with the rest of the metadata), keyed by sample and analyte. DHEA and DHEAS are excluded (large biases that would compress the other panels).
# Reference (target) concentrations come from the metadata's `QC Concentrations`
# sheet (imported above), keyed by `sample_id` + `analyte_id`, and surface via the
# `annot_qcconcentrations` slot (all in nmol/L). Both the SKML EQA reference values
# and the nominal LQC/HQC targets are kept, so accuracy is assessed for all three.
acc_target <- mexp@annot_qcconcentrations |>
filter(str_starts(sample_id, "SKML") | sample_id %in% c("LQC", "HQC")) |>
transmute(sample_id, analyte_id, target = concentration)
# analysis_id -> sample_id map (groups the `…b` replicates onto one EQA sample).
sid_map <- get_analyticaldata(mexp, annotated = TRUE) |>
distinct(analysis_id, sample_id)
d_acc <- d_conc_raw |>
filter(qc_type %in% c("LQC", "HQC", "SPL")) |>
pivot_longer(
!c("analysis_id", "qc_type"),
names_to = "analyte_id", values_to = "conc_mrm"
) |>
left_join(
mh_data |> select(analysis_id, analyte_id = feature_id, conc_mh),
by = c("analysis_id", "analyte_id")
) |>
left_join(
select(spl_area_excl, analysis_id, analyte_id, excl_mrm, excl_mh),
by = c("analysis_id", "analyte_id")
) |>
left_join(sid_map, by = "analysis_id") |>
inner_join(acc_target, by = c("sample_id", "analyte_id")) |>
# Compact axis label: SKML EQA -> 2026.1A; LQC/HQC keep their name.
mutate(
sample_lbl = str_replace(
str_remove(sample_id, "^SKML"),
"^([0-9]{4})_", "\\1."
)
) |>
# Below-Cal-A (below-LOQ) values are kept rather than masked; a per-method flag
# (`blo_*`) lets the panel draw them as open symbols. LQC/HQC carry no exclusion
# flag (never below-Cal-A here), so they count as quantifiable.
transmute(
analyte_id,
sample_lbl, analysis_id, target, conc_mrm, conc_mh,
bias_mrmhub = (conc_mrm - target) / target * 100,
bias_masshunter = (conc_mh - target) / target * 100,
blo_mrmhub = coalesce(excl_mrm, FALSE),
blo_masshunter = coalesce(excl_mh, FALSE)
) |>
filter(!(is.na(conc_mrm) & is.na(conc_mh)))
# x-axis order: LQC, HQC first, then the EQA samples (sorted by label).
eqa_lbls <- sort(setdiff(unique(d_acc$sample_lbl), c("LQC", "HQC")))
x_levels <- c("LQC", "HQC", eqa_lbls)
p_eqa <- d_acc |>
filter(!analyte_id %in% c("DHEA", "DHEAS")) |>
# Reshape bias and the below-LOQ flag together (one row per method).
pivot_longer(
c(bias_mrmhub, bias_masshunter, blo_mrmhub, blo_masshunter),
names_to = c(".value", "method"),
names_pattern = "(bias|blo)_(mrmhub|masshunter)"
) |>
filter(!is.na(bias)) |>
mutate(
method = recode(method, mrmhub = "MRMhub", masshunter = "MassHunter"),
method = factor(method, levels = c("MRMhub", "MassHunter")),
sample_lbl = factor(sample_lbl, levels = x_levels),
loq = factor(
if_else(blo, "Below LOQ (< Cal A)", "Quantifiable (>= Cal A)"),
levels = c("Quantifiable (>= Cal A)", "Below LOQ (< Cal A)")
)
) |>
ggplot(aes(
sample_lbl, bias,
colour = method, shape = loq, group = method
)) +
geom_hline(yintercept = 0, colour = "grey40") +
geom_hline(
yintercept = c(-15, 15),
linetype = "dashed", colour = "grey75"
) +
# Divider between the LQC/HQC group (left) and the EQA samples (right).
geom_vline(xintercept = 2.5, colour = "grey85", linewidth = 0.3) +
geom_point(
size = 1.5,
alpha = 0.85, position = position_dodge(width = 0.4)
) +
scale_colour_manual(values = col_method, name = "Method") +
# Open symbols flag EQA values below the lowest calibrator (below LOQ).
scale_shape_manual(
values = c("Quantifiable (>= Cal A)" = 16, "Below LOQ (< Cal A)" = 1),
name = NULL
) +
facet_wrap(vars(analyte_id), scales = "free_y", ncol = 3) +
guides(
colour = guide_legend(override.aes = list(shape = 16)),
shape = guide_legend(override.aes = list(colour = "grey30"))
) +
labs(
x = NULL,
y = "Bias vs reference (%)",
title = "Accuracy - bias per analyte (QC & EQA; excl. DHEA, DHEAS)"
) +
theme(
aspect.ratio = 1,
panel.grid.minor = element_blank(),
strip.text = element_text(size = 6),
axis.text.x = element_text(angle = 45, hjust = 1, size = 5),
axis.text.y = element_text(size = 5.2)
)
p_eqa
The four panels above are composed into the publication figure (Figure 2).
# Middle row (B, C) uses coord_equal(); give it enough height that the two square
# panels grow to fill the full width, so their L/R borders align with A and D.
# The bias panel is a square, 3-column grid as a standalone figure; inside the
# full-width composite it fills better as a wide 4-column grid without the forced
# square aspect, so override those two settings for Figure 2 only.
p_eqa_fig2 <- p_eqa +
facet_wrap(vars(analyte_id), scales = "free_y", ncol = 4) +
theme(aspect.ratio = NULL)
fig2 <- p_r2 /
(p_qonly + p_full) /
p_eqa_fig2 +
plot_layout(heights = c(0.68, 1.3, 1.2)) +
plot_annotation(tag_levels = "A")
save_plot(
fig2, "output/fig2_quantification_comparison", format = c("pdf", "png"),
width = 180, height = 236)✔ Saved plot (180 x 236 mm) to output/fig2_quantification_comparison.pdf, output/fig2_quantification_comparison.png
output/fig2_quantification_comparison.{pdf,png} (180 mm).
At the integration stage, MRMhub-INTEGRATOR and MassHunter peak areas follow the 1:1 line over four orders of magnitude, with a mean bias below 1 % but 95 % limits of agreement of roughly ±19 %, and comparable precision over the ten replicate SST injections (median %CV 6.3 % for both) (Figure 7). The larger differences are confined to real samples where the peak sits on an elevated baseline or beside a co-eluting neighbour, and the two integrators set the boundaries differently; a number of those peaks were integrated manually in MassHunter.
At the quantification stage, MassHunter areas run through the MRMhub pipeline give the same concentrations as MassHunter end-to-end, so the two post-processing chains agree on identical input. Run end-to-end, each with its own integration and calibration, the two pipelines show comparable bias against the nominal LQC/HQC targets and the SKML EQA reference values (Figure 12). The remaining differences between them therefore enter at the peak-integration step rather than in the QUANT post-processing.