Metabodeconplus — An R Package for Automated Deconvolution and Alignment of Metabolomics Data

1. Introduction

1.1. Background

Nuclear magnetic resonance (NMR) spectroscopy is a widely used analytical technique for detecting and quantifying metabolites in biological specimens such as urine, plasma, and tissue extracts [1]. By exploiting the magnetic properties of atomic nuclei, NMR provides structural and quantitative information about the molecular composition of a sample. Alongside mass spectrometry, this makes it a key tool in metabolomics, where researchers seek to characterize the metabolic profile of biological systems.

1.2. Challenges in 1D NMR Spectra Analysis

The goal of a metabolomics experiment is to identify and quantify the metabolites present in a sample. However, this information cannot be measured directly. Instead, certain chemical properties are recorded, allowing inferences about the contained metabolites. In NMR spectroscopy, the measured properties are the resonance frequencies of atomic nuclei in a magnetic field, which are influenced by their chemical environment, and the intensities at which those resonance frequencies appear. These resonance frequencies, often expressed as chemical shifts, together with their intensities act as a “fingerprint” that can be used to infer the presence and concentrations of specific metabolites [2]. However, despite its advantages, NMR spectroscopy faces challenges that complicate the direct extraction of metabolite concentrations from raw spectra. Two major issues that hinder accurate quantification are peak overlap and peak shifts.

1.2.1. Peak Overlap

Ideally, the resonance frequencies in an NMR spectrum would appear as sharp peaks, allowing straightforward identification and quantification. However, in reality, these peaks possess distinct linewidths and, therefore, may overlap with one another, making it difficult to distinguish individual signals. Several factors contribute to peak broadening, including T₂ relaxation effects [3], magnetic field inhomogeneities [4], chemical exchange processes [5], sample viscosity, molecular interactions, which both contribute to relaxation [6], and the measurement process itself.

The process of reconstructing the individual metabolite contributions from a complex spectrum of overlapping peaks is known as deconvolution. This is one of the key challenges that metabodeconplus seeks to address.

1.2.2. Peak Shifts

Another major challenge in NMR-based metabolomics is that the resonance frequencies of the same metabolite do not always appear at exactly the same positions across samples. Several factors can cause these peak shifts, including variations in temperature [2], differences in sample pH [5], and variations in ionic strength and solvent composition [6].

These inconsistencies complicate direct comparisons between samples and hinder both identification and quantification of metabolites. To address this issue, signal alignment techniques are used to shift peaks across spectra to concordant positions, allowing for more accurate quantitative comparisons.

1.3. Existing Solutions

Available open-source solutions include the R packages BATMAN [7], rDolphin [8], Speaq 2.0 [9], ASICS [10], SigMa [11], and MetaboDecon1D [12] (the predecessor of this package), the Python program Decon1d [13], the MATLAB applications icoshift [14] and COW [15], and the web applications BAYESIL [16], DEEP Picker1D [17], and NMRProcFlow [18]. Proprietary software suites available for NMR spectral analysis include NMR Workbook Suite by ACD/Labs, AMIX by Bruker, Chenomx NMR Suite by Chenomx Inc., and MNova NMR by Mestralab Research. Furthermore, Bruker TopSpin 4.2 includes the program mldecon [19].

A summary of the main features of each tool is provided in Table 1; detailed descriptions of individual solutions are given in Appendix A.

1.4. The Role of Metabodeconplus

The aim of metabodeconplus is to deconvolute 1D NMR spectra of complex biofluids such as plasma, urine, or CSF that contain many overlapping peaks, in order to determine accurate integrals of the underlying peaks regardless of whether the corresponding metabolites are known. As a result, unidentified signals can also be used for classification tasks. We previously developed an R package for the deconvolution and integration of 1D NMR data (MetaboDecon1D [12]). This enabled the deconvolution of overlapping NMR signals. For each spectrum, the output is a list of deconvoluted signals characterized by their spectral positions together with accurate integral values. However, to make full use of these data, additional signal alignment is required for subsequent statistical analyses. metabodeconplus is not a new version of MetaboDecon1D but a new and substantially broader package: whereas MetaboDecon1D performed deconvolution of a single spectrum, metabodeconplus provides a complete, automated pipeline that takes a set of raw 1D NMR spectra all the way to an aligned, classification-ready table of signal integrals. First, the CluPA algorithm [20] of the speaq 2.0 package [9] was included. To this end, CluPA was adapted to the output format of MetaboDecon1D. The CluPA algorithm was chosen because it offers reliable alignment between target and reference spectra and has been implemented in R, which facilitated its addition to the metabodeconplus package. In addition, the underlying deconvolution algorithm was optimised for speed and accuracy.

The contribution of metabodeconplus is thus a single, automated end-to-end workflow for 1D NMR fingerprinting: from raw spectra it produces a two-dimensional table of aligned signal integrals that feeds directly into built-in sample classification. It achieves this by integrating three previously separate steps — deconvolution, peak alignment via the adapted CluPA algorithm, and classification — into one package, and by automating the parameter search these steps require. A substantially faster deconvolution backend makes the complete workflow practical on datasets of several hundred spectra. For a user, this means the entire path from raw spectra to a classification-ready feature matrix is handled by one package, replacing the ad-hoc combination of separate deconvolution and alignment tools that this analysis previously required.

2. Materials and Methods

2.1. Package Availability

The R package is available on GitHub at https://github.com/spang-lab/metabodeconplus and will be available on CRAN via https://cran.r-project.org/package=metabodeconplus.

2.2. Study Cohorts

Four datasets were used throughout this study. The Sim2 and Urine datasets are distributed with the package, whereas the Blood and AKI datasets can be obtained via download_example_datasets() or downloaded manually from GitHub.

The Sim2 dataset consists of 100 simulated spectra, each containing 2048 data points covering a chemical shift range of 3.59 to 3.28 ppm (datapoint spacing $0.00015$ ppm). The spectra are evenly split into two groups, A and B (50 each). Each spectrum contains the same 25 underlying Lorentzian peaks: a fixed set of reference peak positions, peak intensities, and half-widths is reused across all spectra, with per-spectrum variability introduced by (i) a scalar global ppm shift drawn from $\mathcal{N}(0,0.{00120}^2)$ ppm (sd $=8$ datapoints), (ii) a per-peak ppm jitter drawn from $\mathcal{N}(0,0.{00060}^2)$ ppm (sd $=4$ datapoints), (iii) per-peak amplitude scaling drawn from $\mathcal{U}[0.4,1.6]$, (iv) per-peak half-width scaling drawn from $\mathcal{U}[0.9,1.1]$, and (v) additive Gaussian noise. In group A, the amplitudes of 3 of the 25 peaks (the discriminating peaks) are additionally scaled by fixed factors of $1.25$, $1.25$, and $1/1.25$ ($\approx 0.80$); group B is left unmodified. The first spectrum is generated without any positional shift and serves as a clean unshifted reference. The shared reference parameters were derived from an initial deconvolution of the Blood dataset with manually optimized settings, so that the simulated spectra closely resemble real blood plasma spectra while providing exact ground truth for peak parameters, group labels, and noise levels.

The Blood dataset consists of 16 human blood plasma spectra measured with a 1D ¹H Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence.

The Urine dataset consists of 2 human urine spectra measured with a 1D ¹H NOESY pulse sequence. Specimens from both the Blood and Urine data sets were obtained from specimens of the German Chronic Kidney Disease (GCKD) study collected at the baseline time point[21,22]. The GCKD study was carried out in accordance with the Declaration of Helsinki, registered in the German Register of Clinical Trials (DRKS 00003971), and approved by the ethics committees of the participating institutions.

The AKI dataset consists of 106 urinary 1D ¹H NMR spectra acquired with a 1D ¹H NOESY pulse sequence in the context of a study on acute kidney injury (AKI) following cardiac surgery with cardiopulmonary bypass (CPB) use [23]. Urinary specimens were collected 24 h after surgery. These spectra had been previously collected with written informed consent of the patients upon ethical approval from the University Clinic Erlangen, Erlangen, Germany. The study was approved by an institutional review board at the Faculty of Medicine at the Friedrich-Alexander-University (FAU) Erlangen-Nuremberg under the identifier #4010 [24].

All experimental measurements for these datasets were performed on a 600 MHz Bruker Avance III spectrometer (Bruker BioSpin GmbH, Rheinstetten, Germany) employing a helium cooled cryogenic probe and a cooled sample changer. The measurements were performed using optimized pulse sequences for the suppression of water signals.

An overview of all datasets is given in Table 2.

2.3. The Deconvolution Method

The deconvolution algorithm used in metabodeconplus was first described in Koh et al. (2009) and implemented as part of MetaboDecon1D v0.2.2 [12]. The implementation in metabodeconplus was completely restructured into a modular framework with four distinct stages: (i) smoothing, (ii) peak detection, (iii) peak filtering, and (iv) Lorentzian function fitting. A concise overview of the major functional differences between MetaboDecon1D v0.2.2 and metabodeconplus is provided in Table A1. The overall algorithm is visualized for a small, simulated spectrum in Figure 1.

2.3.1. Smoothing

Smoothing is required to suppress spurious detections caused by noise. Following Koh et al. (2009), we used a moving average filter with user-defined window size and iteration count before peak detection.

2.3.2. Peak Detection

Peak detection operates on the smoothed intensity vector using the curvature-based method introduced by Koh et al. (2009) [25]. Peaks are identified based on local minima in the second derivative of the signal. Each peak is represented by a left, center, and right point ${\omega }_i$, with $i=1,2,3$, respectively. Each peak is assigned the score introduced by Koh et al. (2009) [25]: $\mathrm{score}=min\left({\sum }_{i={\omega }_1}^{{\omega }_2}|S^{\prime \prime }\left(w_i\right)|,{\sum }_{i={\omega }_2}^{{\omega }_3}\left|S^{\prime \prime }\left(w_i\right)\right|\right)$, i.e. the smaller of the two cumulative absolute second-derivative sums over the peak’s left and right halves. Sharper, more prominent peaks therefore receive higher scores.

2.3.3. Peak Filtering

Peaks in the signal-free region are used to estimate the noise distribution. Only peaks with a score exceeding $\mu +\delta ·\sigma$ are retained, where $\mu$ and $\sigma$ are the mean and standard deviation of the noise scores and $\delta$ is a user-tunable multiplier (default 6.4) that controls peak-detection sensitivity. Peaks inside user-specified ignore regions, such as the water-artifact region, are also removed.

2.3.4. Lorentzian Function Fitting

Each detected peak provides a triplet of observed points $({\omega }_i,y_i)$, $i=1,2,3$, which uniquely determines the center ${\omega }_0$, half-width $\lambda$, and amplitude A of a Lorentzian function passing exactly through those three points. These are the initial parameter estimates. When the resulting Lorentzian curves are superimposed, the sum tends to exceed the observed spectrum due to spill-over from neighboring peaks. Following Koh et al. [25], each peak’s height is, therefore, iteratively reduced by a rule of proportion until the superposition matches the raw spectrum. Compared to the original MetaboDecon1D implementation, metabodeconplus uses a more efficient strategy for computing the Lorentzian parameters during both initial estimation and iterative refinement; mathematical details are provided in Appendix C.

2.3.5. Rust Implementation

In addition to the native R implementation, metabodeconplus provides an optional Rust backend that implements the same deconvolution algorithm with further refinements for numerical stability and computational efficiency. The backend is accessible through the standard R interface via the use_rust parameter of deconvolute() and yields results that are numerically very close, though not identical, to the R implementation. The substantially reduced runtimes allow for a systematic grid search of optimal deconvolution parameters (Section 2.5.1) even for large datasets consisting of several hundred spectra.

2.3.6. Scoring of Deconvolution Quality

For scoring of deconvolution quality we devised a new measure called PRARPX, which scores both the number of correctly identified peaks and the approximation of the measured spectra. It is defined as $\mathrm{PRARPX}=\mathrm{PRX}·(1-\mathrm{AR})$, where the Extended Peak Ratio $\mathrm{PRX}=c/(t+i)$ compares correctly identified peaks c to true peaks t plus incorrectly identified peaks i, and the Area Ratio $\mathrm{AR}=R/S$ is the total absolute residual R over the total absolute observed signal S. Both factors lie in $[0,1]$ and reach 1 only when the correct peaks have been recovered and the reconstruction matches the observed spectrum, so PRARPX approaches 1 for an ideal deconvolution. Because PRARPX requires ground-truth peak parameters, it is only applicable to simulated data or controlled measurements of known compounds; Appendix D motivates it in detail and contrasts it with using PR or AR alone.

2.4. The Alignment Method

Proper signal alignment is crucial to correct for chemical-shift variations across samples and to ensure that subsequent statistical analyses compare the same signals from the same metabolites. Following deconvolution, the obtained peak lists are passed through two consecutive steps: a coarse CluPA (cluster-based peak alignment) shift, followed by a refinement step snap_to_ref(). Both operate on the discrete peak grid: CluPA shifts whole spectrum segments by the same integer-datapoint offset, so all peaks within one segment move together — one peak may end up well-aligned while another in the same segment is still slightly off. The snap_to_ref() step removes this residual mismatch by snapping each remaining peak individually onto the nearest peak column of the reference spectrum. The two steps are exposed individually as clupa() and snap_to_ref() and are chained into a single call by the wrapper align(x, maxShift, maxCombine). The combined procedure is visualized for a set of simulated spectra in Figure 2.

2.4.1. CluPA

The first stage is based on the hierarchical cluster-based peak alignment (CluPA) approach from the speaq package [9,20]. CluPA aligns a target spectrum to a reference spectrum in two phases. First, it shifts the spectrum globally. Then, it iteratively splits the spectrum into smaller local segments that are individually aligned. The segment boundaries are defined by hierarchical clustering applied to the combined peak list of reference and target spectrum: peaks are grouped into a tree based on their distances, and at each iteration all segments are bisected according to the dendrogram. For each segment, the optimal integer datapoint shift is determined by fast Fourier transform cross-correlation, bounded by the user-supplied tolerance maxShift. Iterations stop when a segment contains only peaks from the same spectrum or fewer than three peaks. CluPA changes peak centers but does not change the number of peaks per spectrum.

2.4.2. snap_to_ref

CluPA shifts whole spectrum segments, not individual peaks, so peaks within the same segment all move equally in terms of direction and distance. Following CluPA, the alignment of spectra is substantially improved, but individual peaks are still not guaranteed to exhibit the exact same chemical shift across spectra, and a refinement step is needed to enforce a common peak grid. The snap_to_ref() step does this by collapsing (snapping) each spectrum’s peak list onto the peak grid of the reference spectrum chosen by CluPA. For each target peak, the nearest reference peak within maxCombine columns is its target; peaks farther than maxCombine from every reference peak are dropped. This may lead to the loss of rare peaks present in only few of the analyzed spectra and not in the reference spectrum. The interpretation of maxCombine is, therefore, the residual positional uncertainty that remains after CluPA. Figure 2 visualises the cumulative effect of the two steps: the top row overlays the six reconstructed spectra after each step, the middle row renders them as intensity heatmaps, and the bottom row reduces each spectrum to its sparse peak-position marks.

Both clupa() and snap_to_ref() align/snap each spectrum towards a fixed reference spectrum. This makes application to new samples trivial: each new sample is aligned and snapped to the same reference.

2.5. Parameter Optimization

2.5.1. Unsupervised Parameter Optimization

metabodeconplus’s deconvolution exposes four tunable parameters: nfit, smit, smws and delta. nfit sets the number of iterative Lorentzian-refinement passes, smit the number of moving-average smoothing passes, smws the smoothing window size (in datapoints), and delta the noise-threshold multiplier from Section 2.3.3. Without class labels, these can be optimized for each spectrum individually in an automated fashion by minimizing the area ratio (AR) between observed intensities $y_i$ and reconstructed intensities ${\widehat{y}}_i$:

$$\mathrm{AR}={\displaystyle \frac{{\sum }_i|y_i-{\widehat{y}}_i|}{{\sum }_i\left|y_i\right|}}.$$

(1)

For evenly sampled spectra this is the mean absolute residual divided by the mean absolute signal, and smaller values indicate a better fit.

Minimizing the AR alone is not always optimal as the AR keeps decreasing as more Lorentzians are fitted, so the optimum may drift towards over-fit reconstructions with many spurious peaks. We, therefore, added as a user-defined parameter the maximum number of expected peaks npmax. The automated parameter optimization of the deconvolution step then sweeps a 60-cell parameter grid $\mathtt{smit}\in \{1,2,3\}$, $\mathtt{smws}\in \{3,5,7,9\}$, $\mathtt{delta}\in \{1.6,3.2,4.8,6.4,8.0\}$, $\mathtt{nfit}=10$ and picks the smallest-AR cell whose peak count does not exceed npmax. As a consequence the user has to set only npmax manually. A suitable npmax can be set from domain knowledge (e.g., expected maximum peak numbers in blood or urine spectra) or chosen by maximizing downstream classification performance in case that class labels are available (Section 2.5.2).

The 60-cell grid above was chosen by deconvoluting representative blood and urine spectra over a broader candidate grid, visually inspecting the resulting Lorentzian reconstructions, and trimming the grid to the values that consistently produced plausible fits. The delta step of $1.6$ was picked so that the previous hand-tuned default $\mathtt{delta}=6.4$ falls on a grid point. Users who wish to use a finer control can pass an explicit deg argument.

2.5.2. Supervised Parameter Optimization

When class labels are available, npmax, maxShift and maxCombine can be optimized jointly by maximizing classification performance on the resulting feature matrix. Because the same metabolites may not be present in all spectra, zeros in the feature matrix encode “missing peak” rather than “zero area”. Linear models confound these zeros with small observed values and become unstable; we therefore use a non-linear random forest classifier from the ranger package [26], which can partition on presence as well as magnitude. A further advantage of random forests is the out-of-bag (OOB) error [27], a generalization-error estimate obtained during training itself, which we use as our score without an extra held-out split or inner cross-validation loop.

The two alignment parameters — maxShift (CluPA tolerance) and maxCombine (snap-to-ref window width) — are governed by different sources of positional error and their empirical optima do not necessarily coincide, so tuning them jointly is in general not redundant. All three parameters can therefore in principle be swept as a three-dimensional grid via the supervised search in fit_mdm().

In practice, we do not recommend tuning all three at once: the runtime grows multiplicatively with the size of each axis, and the deconvolution and CluPA stages dominate that cost. We instead recommend the following two-step procedure:

1.

Pick npmax from domain knowledge (e.g., urine spectra typically contain more metabolites than blood) or by visual inspection of a representative deconvolution, as in Section 2.5.1.

2.

Rely on CluPA’s internal optimum by passing the sentinel maxShift = -1. This sentinel triggers an adaptive sweep that doubles maxShift through $\{1,2,4,8,\dots \}$, runs CluPA at each step, computes the average pairwise Pearson correlation of the aligned Lorentzian superpositions, and stops one step before that correlation first decreases.

The supervised search is then reduced to a one-dimensional sweep over maxCombine, which is by far the cheapest of the three stages to rerun.

fit_mdm() relies on caching and multiprocessing to keep its runtime reasonable. Three optimizations are combined:

1.

Each npmax value has an associated set of per-spectrum deconvolution parameters $(\mathtt{nfit},\mathtt{smit},\mathtt{smws},\mathtt{delta})$ that give the lowest reconstruction error. Finding them requires an internal grid search per spectrum, which is run once at function entry and attached to each spectrum. Whenever npmax changes the optimal parameters can be looked up instead of repeatedly recomputed.

2.

The grid is traversed in (npmax, maxShift, maxCombine) order, so npmax varies slowest. If npmax is unchanged between two rows, the deconvolution of the previous row is reused instead of recomputed; if npmax and maxShift are both unchanged, the CluPA alignment is reused as well.

3.

Deconvolution, alignment and fitting each use several workers. We parallelize over spectra rather than over grid rows to keep memory low: parallelizing over grid rows would force every worker to hold a copy of all spectra, whereas parallelizing over spectra means each worker only holds the spectra it is currently processing.

3. Results

3.1. Deconvolution Quality on Sim2: Metabodeconplus vs. MetaboDecon1D and Grid Search

The evaluations that follow quantify the payoff of the integrated metabodeconplus workflow for the user: automatically tuned deconvolution and alignment, fast enough to scale to hundreds of spectra, at no cost in downstream classification performance. We scored deconvolution quality with PRARPX as defined in the materials and methods section. To this end, we devised a simulated dataset Sim2 (see Table 2 for details). Note that this simulated dataset also contained noise to approximate a realistic experimental data set. Figure A2 shows an example spectrum of this data set. Each of the 100 Sim2 spectra was deconvoluted with MetaboDecon1D (default parameters), metabodeconplus (default parameters), and metabodeconplus with unsupervised grid-search parameter selection at $\mathtt{npmax}\in \{10,20,30,40,50\}$. The signal-free region was auto-derived per spectrum so that all methods saw equivalent inputs. Note that PRARPX is computed from ground-truth peak parameters, whereas the general grid search parameter optimization has to operate on the area-ratio alone. Therefore, a PRARPX of 1 is unreachable in practice.

Table 3 summarises PRARPX (mean, sd, min, max) across all spectra for each configuration, and Figure 3 plots the per-spectrum trace. The default metabodeconplus configuration already improved on MetaboDecon1D, and the grid-search cutoff at moderate npmax closed most of the remaining gap to the metabodeconplus optimal bound.

Even with the metabodeconplus optimal parameters an ideal value of $\mathrm{PRARPX}=1$ is not reached for two reasons. First, the superposition of the estimated noiseless Lorentzian signals cannot fully reproduce the noisy observation. Figure A2A illustrates this. Second, small peaks completely shadowed by larger neighbours produce no detectable curvature minimum and are therefore not recovered, so $\mathrm{PRX}<1$ for ordinary spectra (Figure A2B).

3.2. Alignment Quality on Sim2: CluPA and snap_to_ref

To quantify alignment quality directly we exploited the fact that every Sim2 spectrum contains the same 25 simulated peaks named A–Y at known positions, with the first spectrum simulated without any per-spectrum jitter or global shift so that it serves as the unshifted reference. For each extracted peak we recorded to which of these 25 simulated peaks it corresponded to: a peak whose fitted center lay within $\pm 2$ data points of a simulated position inherited that peak’s name; a peak farther away from any simulated position was labeled Z (a spurious detection). Labels are assigned to the raw deconvolution output, before any alignment, so each peak carries the identity of the simulated source it came from through the rest of the pipeline.

After CluPA each peak is potentially shifted but still positioned on the continuous chemical-shift axis. The snap_to_ref() step then collapses these onto the discrete peak grid of the unshifted reference spectrum: each peak is moved to the nearest reference column within a snap window of maxCombine datapoints, and peaks outside of this window are dropped. With labels in hand we asked, for each spectrum and each true-peak, whether the reference column belonging to a certain peak ended up populated by the right peak. The fraction of correct pairs over the full $100\times 25$ grid is the snap purity reported in Figure 4.

Results show that the snap purity is highest where CluPA has already brought peaks within a small maxCombine window of their reference columns. On Sim2 the maxCombine axis shows no valley because the 25 simulated peaks are spaced far enough apart (tens of data points) that even the widest snap windows rarely sweep in a wrong neighbour; on a more crowded spectrum, widening maxCombine past the optimum would start trading correct snaps for wrong-label intrusions and a clear valley would appear (data not shown).

3.3. Supervised Parameter Optimization on Sim2

In the following, it is demonstrated that the supervised search described in Section 2.5.2 where class labels are known is feasible. To this end, we ran a 5000-tree ranger probability forest over a small grid. The grid spans all three parameters ($\mathtt{npmax}\in \{10,20,30,40,50\}$, $\mathtt{maxShift}\in \{2,4,8,16,32\}$, $\mathtt{maxCombine}\in \{1,5,10\}$). This Sim2 sweep is the only instance where all three parameters were tuned at once: the small spectrum count and coarse grid keep the runtime under a few minutes. This allowed to analyze whether the supervised optimum of parameters obtained here is consistent with the unsupervised one. The grid search selected $\mathtt{npmax}=30$, $\mathtt{maxShift}=4$, $\mathtt{maxCombine}=5$ as the best combination (OOB accuracy $90\%$, OOB AUC $94.6\%$). Applied to the held-out 50 test spectra, this configuration achieved a test accuracy of $86\%$ and AUC of $81.9\%$. Of the top-10 features ranked by ranger permutation importance, 3 fall within 3 datapoints of a discriminative peak.

The supervised optimum $\mathtt{npmax}=30$ coincides with the unsupervised optimum of Section 3.1. All three discriminating peak positions are recovered in most spectra. Note, the supervised search optimises classification accuracy rather than reconstruction fidelity as in the unsupervised search, so the two criteria need not select the same $\mathtt{npmax}$ in general; in contexts where peak recovery rather than classification accuracy is the criterion of interest, $\mathtt{npmax}$ is more naturally set to a prior upper bound on the expected peak count than tuned by the supervised search.

The out-of-bag (OOB)-AUC and OOB-accuracy heatmaps obtained for the different parameter combinations are shown in Figure 5 a and b, respectively. A heatmap of all features selected in the training data set by the ranger classification (Figure 5 c). Data are shown for the optimal set of parameters obtained in the supervised search. Each column is labeled with its corresponding chemical shift position in ppm. Columns to the left of the center are features whose mean is higher in Group A (positive two-sample t-score) and vice versa for columns to the right. Within each half , columns are sorted by ranger feature importance [26,27] so that the most informative features are located at the left or right borders of the plot. A feature is counted as matching a discriminative peak when its center lies within $\pm 3$ data points of one of the 3 discriminatory peaks. Figure 5 d shows a superposition of 5 group-A and 5 group-B training spectra after alignment and snapping, zoomed on the range of ranger features (plus $0.05$ ppm of padding), with one vertical line per feature shown in (c). Green lines mark features within $\pm 3$ datapoints of a discriminative peak.

Data show that good classification performance is obtained with automated parameter optimization. As can be seen from Figure 5 a and b classification performance is relatively robust with respect to parameter settings with a broad maximum in the parameter matrix. Figure 5 c shows that the true discriminatory features are amongst the most informative features picked by the classification algorithm. That also other features are considered by the classification algorithm is due to the fact that the three discriminatory features of group A were up and down scaled only very moderately. Furthermore, a per-peak amplitude scaling and the addition of Gaussian noise were performed on the Sim2 data set as described in the Study cohorts section. As can be seen from Figure 5 d the three discriminatory peaks are of only moderate intensity and show partial overlap with other signals. Therefore, the Sim2 data set provides a challenging realistic example.

3.4. End-to-End Prediction Performance on the AKI Dataset

We evaluated end-to-end predictive performance on the AKI dataset (34 AKI versus 72 control urinary spectra, creatinine-normalised) by comparing two models. The first was an equidistant-binning baseline that mirrors the fixed-bin feature representation originally used by Zacharias et al. [23], who introduced the AKI dataset: a 700-bin feature matrix consisting of 300 bins covering $6.5$–$9.5$ ppm and 400 bins covering $0.5$–$4.5$ ppm (a fixed bin width of $0.01$ ppm). Whereas Zacharias et al. [23] classified these features with a support-vector machine (SVM), we replaced the SVM with a ranger-based random forest [26,27,28] run in probability mode (one class probability per sample), matching the classifier used by the metabodeconplus pipeline and ensuring a fair comparison. The second was the full metabodeconplus pipeline: deconvolute → CluPA → snap_to_ref → ranger, with $\mathtt{npmax}=1000$ fixed.

Following the recommendation of Section 2.5.2, we fix $\mathtt{npmax}=1000$ on the AKI dataset rather than tuning it. The choice was made by visual inspection of the deconvolution output on the two Urine spectra introduced in Table 2: at $\mathtt{npmax}=1000$, the fitted Lorentzian superposition reconstructs every visually-resolvable peak of those reference urine spectra while keeping the residual free of obvious unfit signal, so the same setting was expected to be a sensible default for the AKI urine cohort, which was acquired with the same NMR pulse sequence. We further rely on CluPA’s internal optimum for $\mathtt{maxShift}$ by passing the sentinel $\mathtt{maxShift}=-1$ (see Section 2.5.2), so that only $\mathtt{maxCombine}$ is swept on the training fold of each cross-validation split via the inner grid-search in fit_mdm().

Predictive performance was estimated by stratified 10-fold cross-validation repeated over three random seeds (so 30 folds in total per model). Each fold’s training portion was independently re-fit before scoring on the held-out fold, and standard errors are taken across folds. The binning baseline reached $75.4\pm 2.33$ % accuracy and $\mathrm{AUC}=0.810\pm 0.027$; metabodeconplus reached $74.3\pm 2.05$ % and $\mathrm{AUC}=0.818\pm 0.023$, thus matching the binning baseline with no loss in accuracy or AUC, while operating on a higher-resolution feature matrix whose columns correspond to individually localised peaks rather than fixed-width bins.

To compare what the two models had learnt, we refitted both on the full AKI dataset, extracted ranger permutation importances, and overlaid the top 20 features of each model on the alignment reference spectrum across the $0.5$–$4.5$ ppm aliphatic region (Figure 6).

On its higher-resolution, per-peak feature matrix, metabodeconplus performed on par with the binning model. One main advantage of using a deconvolution approach is that each selected feature is a single deconvoluted peak that contributes to only one metabolite. This allows for an unambiguous metabolite assignment of selected features. This is of special importance in case the measurement of predictive signatures should be transferred to other methods from, for example, clinical chemistry. In addition, this allows for a biological interpretation of obtained signatures. In this context it should be noted that predictive signatures only provide correlations between metabolites and groups of specimens; they do not provide causal relationships. Whereas, as can be seen from Figure 6, for the binning model a bin may contain contributions from multiple peaks and metabolites, which substantially hinders assignment of predictive signatures to corresponding metabolites.

3.5. Runtime Performance and Parallel Scaling

We compared the per-spectrum, single-core runtime of the original MetaboDecon1D implementation (Figure 7 A), the R implementation of the new metabodeconplus package (Figure 7 B), and the Rust backend of metabodeconplus as activated from within the R environment (Figure 7 C). To this end, spectra with different numbers of data points ($2^{11}$ to $2^{17}$) and peak counts ($2^4$ to $2^{12}$) were simulated. Typical 1D ¹H spectra of biofluids such as urine and plasma consist of $2^{16}$ to $2^{17}$ data points with approximately 1000 peaks. For the original MetaboDecon1D implementation, this resulted in an average runtime of around 40 s at $2^{17}$ data points (Figure 7 A). In the R implementation of the new metabodeconplus package, the runtime is substantially reduced to well below one second per spectrum in this regime (Figure 7B). The Rust backend (Figure 7 C) provides a further speed-up that becomes most pronounced for densely-peaked spectra: at $2^{17}$ data points and 4096 peaks the deconvolution takes around $3.7$ s in pure R but only around $2.4$ s with Rust, while at the more typical $\approx 1000$ peaks the two backends are essentially indistinguishable on a single spectrum. Panels A–C of Figure 7 are single-core single-spectrum measurements; parallelisation across spectra is shown separately in Figure 7 D showing the deconvolution times of all 106 AKI spectra with respect to the degree of parallelization ($\mathtt{nworkers}\in \{1,\dots ,10\}$) for both the R and Rust backend. Runtime scales close to the ideal $1/k$ reference throughout the measured range, with the Rust backend completing the full 106-spectrum batch in under 6 s at $\mathtt{nworkers}=10$ and the R implementation in under 8 s. These short runtimes make it tractable to determine the optimal deconvolution parameters for each dataset individually by systematic grid search.

4. Discussion

With metabodeconplus we provide a single end-to-end workflow for the analysis of 1D NMR fingerprinting data of complex biofluids, requiring no prior knowledge of the individual molecules present. The workflow integrates three steps that previously required separate tools: deconvolution of overlapping signals into individual Lorentzian lines, alignment of the deconvoluted peaks, and sample classification, tied together by an automated parameter search. Results showed that, with the new implementation of the deconvolution algorithm, a substantial improvement in runtime could be obtained. With this the deconvolution of large datasets containing several hundreds of spectra becomes feasible. This is especially true when the Rust implementation of the deconvolution part is used. This substantial improvement in runtime permits a systematic search of optimal deconvolution and alignment parameters. To analyze the performance of metabodeconplus in terms of peak detection and deconvolution, a carefully designed challenging simulated data set was used, where the ground truth is known. We could show that metabodeconplus facilitates a reliable peak detection and deconvolution of overlapping signals in complex data. The residual gap to $\mathrm{PRARPX}=1$ on the simulated data is dominated by the additive-noise floor (raising AR above zero) and by small peaks shadowed by larger neighbours (suppressing PRX), as detailed in Section 3.1 and illustrated in Figure A2. Deconvolution quality depends on appropriate parameter selection. We devised a grid search procedure to search parameter space by optimizing the area ratio. However, especially for complex experimental spectra it is not guaranteed that this will always lead to an optimal selection. One critical parameter is the expected maximal number of peaks npmax that ensures that not a large number of spurious signals is detected. For experimental data such as urine and plasma spectra we recommend setting the the maximum number of expected peaks based on prior knowledge e.g. to 1000 for human urine. To make full use of the deconvoluted data in subsequent statistical analysis proper peak alignment is required to correct for small variations in peak positions across spectra. We employ the well established CluPA algorithm that was extended by the snap_to_ref function to ensure that corresponding peaks are put to the exact same position. As we could show, this in general provides a good alignment. However, alignment quality depends on the two parameters maxShift and maxCombine. The parameters npmax, maxShift and maxCombine may be manually set by the user after careful spectra inspection. We also devised an automated supervised grid search via fit_mdm() that can in principle tune all three jointly when two well-separable classes with known class labels are available. In practice we recommend, as mentioned above, fixing npmax from domain knowledge or visual inspection, leaving maxShift to the adaptive CluPA-internal sweep selected by maxShift = -1, and tuning only maxCombine on the labelled data — this avoids the multiplicative runtime of the full three-dimensional grid without sacrificing classification performance. The higher spectral resolution and interpretability of the peak-based feature matrix were achieved without loss of predictive performance: classification matched the equidistant-binning baseline, while each discriminatory feature could be assigned to an individual metabolite — which a coarse bin, mixing several peaks, does not permit. In summary, metabodeconplus makes 1D NMR fingerprinting a single automated workflow — from deconvolution through alignment to classification — with automated parameter optimization and a fast, parallelised backend that make this end-to-end pipeline practical at the scale of modern metabolomics cohorts.

Limitations Direct file-format support is restricted to Bruker and JCAMP-DX; other vendor formats are not parsed natively. Users can however hand spectra in as a chemical-shift vector and a corresponding signal-intensity vector, so any format that can be loaded into the R session is usable. Native readers for further vendor formats remain future work. No routines are provided to group deconvoluted peaks into multiplet structures, and no automatic peak-to-metabolite assignment (e.g. matching against reference spectra) is included. Finally, PRARPX has a hard upper bound below 1 on noisy spectra: additive measurement noise keeps $\mathrm{AR}>0$ and peaks completely overlapped by larger neighbours produce no detectable curvature minimum and are therefore not recovered (Figure A2).