Linear MALDI-ToF simultaneous spectrum deconvolution and baseline removal

Smoothing

A popular [5] smoothing technique in the spectrometry community is the use of Savitzky-Golay linear filters [6, 7]. These moving average filters perform a least squares fit of a small set of consecutive data points to a polynomial. The value of this fitted polynomial at the window central point is the filter output. One can also compute a smoothed derivative by using the derivative of the fitted polynomial to compute the central point value. This smoothed derivative can also be used by peak picking algorithms [4]. This method is versatile and efficient. However, its main drawbacks are that we have to manually choose the polynomial degree and the window length. Some studies to automatically choose the former [8] or later [9] exist but to the best of our knowledge they are not used as often as the original approach.

A more recent approach to smooth spectra is the use the wavelet transform. The Undecimated Wavelet Transform (UDWT) [10, 11] is generally preferred to the Discrete Wavelet Transform (DWT) as it produces less artifacts after coefficient thresholding. The UDWT is equivalent to an averaged DWT computed for all integer shifts of the signal and is thus a redundant and shift invariant transform. In applications it has been reported to yield better qualitative denoising [12].

Baseline correction

Baseline correction is a difficult problem that potentially also introduces artifacts [13]. There are at least two kinds of approaches for baseline correction. One category of methods is close to mathematical morphology. In these methods a lower envelope of the spectrum [14, 15] is computed. Methods of this category generally need a smoothed signal (see “Smoothing” section). The other category contains methods using an asymmetric loss function to fit spectrum baseline without being biased by peaks [16, 17]. Finally some other methods mix the two previous approaches [18].

Peak picking

After baseline removal the next step is generally a peak picking procedure. Several approaches are possible. Perhaps the most intuitive approach is to compute a regularized second order derivative (using Savitzky-Golay for instance) of the spectrum and to extract local minima [4]. The use of second order derivative minima instead of the zero-crossing of the first order derivative allows, to some extend, to detect overlapping peaks [19].

A second kind of approach, especially useful in case of overlapping peaks, is peak deconvolution. Overlapping of complex peak patterns can be deconvolved if one uses specially tuned point spread function and judicious regularizations (positivity constraint and sparsity-inducing norms, like the l₁ norm) [20–23]. Further generalizations can be obtained in case of blind-deconvolution [24]. However these kinds of approaches are much more computationally intensive and are not widely used in mass spectrometry.

Finally we can mention the Continuous Wavelet Transform (CWT) [25, 26] which can be efficiently computed using the Fast Fourier Transform. The idea is to follow wavelet modulus maximum. Theses ridges characterize the regularity of the signal [27] and can be used to detect peaks.

Problem definition

where e includes measurement and model errors. It is a common linear model with additive uncertainties. These quantities are represented by vectors of size n, where n is the number of m/z channels of the original spectrum y.

As a consequence, information must be accounted for regarding the expected signals x _p and x _b . In the following developments, x _b is expected to be smooth while x _p is expected to be spiky and positive. This knowledge will be included in the next sections.

Interpretation of the reduced criterion

The two equations are very similar apart from the fact that the operator A _μ is applied to the reconstructed peaks Lx _p and to the raw spectrum y. It could be interpreted as the precision (inverse of the covariance) matrix in a correlated noise framework. Instead of this classical approach and to better understand its role in our deconvolution context we study the evolution of A _μ in the two limiting cases \(\mu \rightarrow 0\) and \(\mu \rightarrow \infty \).

Analysis of the A _μ matrix

Figure 1 shows the column \(j=\lfloor {\frac {n}{2}}\rfloor \) of A _μ for two different values of μ.

Hence compared to Eq. 14, one can interpret Eq. 15 as an usual deconvolution applied on the second order derivative A _μ y of the initial spectrum y. The used point spread function is also the second order derivative A _μ L of the initial peak shape p introduced in Eq. 1. The overall effect of the A _μ operator is to cancel the slow varying component of the signal. In another terms, the baseline is removed thanks to a derivation.

Debiasing

The l₁ penalty acts as a soft threshold to select peaks ([29], Section 10). This leads to a bias in peak intensity estimation. These intensities are artificially reduced when the l₁ penalty increases. We use the ideas introduced in [28] to get corrected peak intensities. This yields a resolution procedure involving two stages. The first stage selects the peaks, the second one corrects their intensities.

Boundary conditions

There is a possible improvement of the method concerning boundary conditions. As Eq. 2 suggests, a strong μ penalty forces the baseline to be constant. In Fig. 2 we see that this phenomenon is especially present at the domain boundaries.

We solved this problem by imposing baseline values at boundaries. This corrected solution is also shown in Fig. 2 and we can see that the corrected solution does not suffer from boundary effect anymore. Appendix “Boundary correction”, page 11, provides all the details on how to modify Eqs. 9 and 13 to introduce some constraints on the baseline values x _b . These modified equations will constitute our final model formulation.

Final model formulation

As before, \(\widehat {\mathbf {x}}_{b}\) is computed from \(\widehat {\mathbf {x}}_{p}\) using Eq. 18. The explicit forms of \(\widetilde {\mathbf {B}}_{\mu }\), \(\widetilde {\mathbf {y}}\) and \(\widetilde {\mathbf {A}}_{\mu }\) are given in Appendix “Boundary correction”, respectively by Eqs. 31, 32 and 33.

Algorithm summary

For ease of reading Algorithm 1 recapitulates the main steps of the proposed method in its final formulation. The optimization procedure used to solve efficiently the two minimization problems will be described in detail in the next section.

Quadratic programming with bound constraints

To solve the optimization problems Eq. 12 or Eq. 20 and their associated debiasing step Eq. 17 we use the projected Barzilai-Borwein method described in [30]. The Barzilai-Borwein method [31] dramatically improves the classical steepest descent with Cauchy step.

For its simplicity and good behavior in practice we have chosen to use this method. Our implementation is given in pseudo-code in Algorithm 2.

Variants of this method with non-monotone line search [36, 37] have been tested. Theoretically this allows to prove global convergence, but in practice we have observed a performance degradation compared to the simpler method presented in [30].

Stopping criterion

Illustration of the method

We give in Fig. 3 typical convergence behavior for the proposed method applied to problem Eq. 12 or Eq. 17.

We observe the non-monotonic convergence behavior of the Barzilai-Borwein method. As reported by [30], we observe that it is more effective to alternate between BB1 Eq. 22 and BB2 Eq. 23 steps than using only one type of step update.

Synthetic data

Algorithm implementation

We provide a reference implementation¹ that can be used to reproduce the results of the following sections. This implementation is coded in C++ and runs under Linux. The main page of the project details all the steps to reproduce the results of the “Joint baseline computation-peak deconvolution”, “Comparison with the usual sequential approach” and “Real data” sections. Typical run-times are one second for the synthetic data and three seconds for the examples using real spectra data. To keep this implementation as simple as possible a basic projected gradient descent is used instead of the more effective Algorithm 2. Also note that this implementation requires CSV input files. A more versatile version of our algorithm, using Algorithm 2, is used in [39]. However, this implementation is not yet publicly available.

Joint baseline computation-peak deconvolution

The baseline approximation obtained at the end of the first step is given in Fig. 4. This figure shows a common flaw of most of the baseline removal methods which is an ascent of the baseline under the peaks.

The deconvolved peaks obtained after the first step are given in Fig. 5. We notice the negative impact of a strong λ₁ penalty which leads to an underestimation of the peak heights. The figure also shows that the deconvolved solution is less spiky in regions of strong peak overlaps (right part of the spectrum).

The baseline computed after the second step is shown in Fig. 6. During this second stage, the λ₁ penalty is removed and replaced by a restricted support Ω computed using Eq. 16. We see that the ascent of the baseline under the peaks has been corrected (c.g. x-axis ranges from 150 to 250 and from 350 to 400).

The new peak heights are shown in Fig. 7. The main role of this second stage is to debias peak heights. Compared to Fig. 5 we can see that this objective is quite well fulfilled.

Comparison with the usual sequential approach

For each processed spectrum we perform a systematic parameter grid search and report the best result obtained by the two methods. To find the best solution we compare the deconvolved peaks with the ground truth of our synthetic data. We use the following procedure to compare these two peak lists. As defined in our model, a peak list is stored in a sparse vector x _p . To each nonzero component \(\mathbf {x}_{p}\left [i\right ]\) corresponds a peak at position i and height \(x_{p}\left [i\right ]\). Like computed peak centers and ground truth peak centers might not be perfectly aligned we reconstruct a continuous signal by convolving the two peak lists by the known peak shape p. This blurring is required later to compute the Euclidean distance between these two peak lists. This reconstruction is performed by the following matrix-vector products: \(L\mathbf {x}_{p}^{\star }\) for the ground truth peak list and \(L\hat {\mathbf {x}}_{p}\) for the peaks extracted by the algorithm. Using these vectors we can now compute a normalized Euclidean distance defined by:

$$\mathcal{E}_{e} =\frac{\big\|\mathbf{L}\left(\widehat{\mathbf{x}}_{p}-\mathbf{x}_{p}^{*}\right)\big\|_{2}}{\big\|\mathbf{L}\mathbf{x}_{p}^{*}\big\|_{2}} $$

For each noise level and each algorithm the computations are performed 10 times, each time with a different noise realization. From these 10 replica we report the obtained mean and standard deviation of the \(\mathcal {E}_{e}\) merit factor. The results are given in Table 4 and plotted in Fig. 8.

σ noise	Sequential processing	Proposed method
	\(\mathcal {E}_{e}\times 10^{-3}\)	\(\mathcal {E}_{e}\times 10^{-3}\)
0.	15.	27.6
0.2	47.8±8.76	45.2 ±9.12
0.4	88.6±15.6	81.6 ±11.9
0.6	131.±21.4	118. ±18.
0.8	176.±26.	154. ±20.1
1.	217.±33.7	184. ±21.3
1.2	262.±45.7	215. ±25.5
1.4	295.±48.2	245. ±29.4
1.6	330.±52.5	275. ±32.7
1.8	360.±59.8	300. ±38.2
2.	390.±70.4	322. ±44.3

\(\mathcal {E}_{e}\) merit factor, lower is better. Apart from the noise free signal, the joint approach outperforms the sequential one for all noise levels

The italic means the best result (when comparing the two methods)

On this example the joint evaluation of the baseline and of the deconvolved peaks outperforms the sequential approach in nearly all configurations except for the noise-free case. In this condition, as we are using synthetic data, a perfect reconstruction of the “ground truth” spectrum is possible. By consequence, the reconstruction error is very small. To explain this result, we think that our iterative solver had stopped its iterative resolution too early.

To get more intuition on what is happening we provide Fig. 9 illustrating the baseline computation using the SNIP algorithm.

As explained the SNIP algorithm needs to work on a previously smoothed spectrum. This smoothing is performed using a Savitzky-Golay filter. This sequential approach, smoothing then SNIP baseline, can introduce an important bias on the computed baseline. This error is then transferred to the peak picking algorithm. The result can be a poor peak extraction.

Baselines obtained by our method and by the classical SNIP algorithm are shown in Fig. 10.

We see that the simultaneous baseline and deconvolution approach allows a nearly perfect reconstruction of the baseline. Unlike our approach, the SNIP algorithm suffers from initial filtering of the spectrum and presents baseline ascent below peaks.

Real data

We present some results obtained by our algorithm on real MALDI-ToF mass spectra. The Fig. 11 shows the peaks obtained when the low resolution spectrum is deconvolved with a tight peak shape. We have used a Gaussian peak shape with σ_peak=0.4 m/z. The goal is to evaluate the behavior of the algorithm when we try to deconvolve isotopic motifs. The notable result here is to see that the returned peak centers are approximately spaced by 1 m/z which is the expected value. This result is encouraging as this information was not provided to the algorithm to perform its deconvolution.

However isotopic motif deconvolution without using any extra information (like an expected 1 m/z spacing between peaks) can lack robustness. That is the reason why it is certainly safer to use a wider peak shape to model the unresolved isotopic motif as a whole. This is illustrated by Fig. 12. We have run the algorithm twice. The first run used a sparsity promoting regularization of λ₁=1, the second run used λ₁=0.5.

The developed algorithm can also be used with high resolution MALDI-ToF mass spectrum (reflectron mode). Figure 13 gives such an example. A Gaussian peak shape with σ_peak=0.15 m/z is used for the deconvolution. As before, we have used two values for the sparsity promoting regularization. The first run used λ₁=0.5, the second run used λ₁=0.2. We clearly see the impact on the results, in the first case only the main peaks are extracted, in the second one very small peaks are also extracted. For the moment this λ₁ parameter has to be manually tuned by the user.

It is generally quite difficult to evaluate peak picking methods on real spectra because we do not have the ground truth at our disposal. However, in our case we have expensively used the proposed approach in the BHI-PRO project. The study [39] compares our method against a more usual one on mass spectra of spiked proteins. This study quantifies the algorithmic part of the variance of the measured protein abundances and shows a clear gain in favor of our algorithm.

Smoothness and convolution matrix

In the presentation the convolution product between the peak list and peak shape function, x _p ∗p, is often represented by the Lx _p matrix-vector product. However, in practical computations the L matrix is never explicitly formed and the convolution product is computed directly using a specialized subroutine.

Derivation of reduced objective

This appendix focuses on the reduced criterion Eq. 10.

Expression of \(\underset {\mu \rightarrow \infty }{\lim } \mathbf {A}_{\mu }\)

Let’s write down two basic statements about eigenpair \(\left (\nu,\mathbf {v}\right)\) of any square matrix M:

Boundary correction

Let begin by writing the Lagrangian [38] associated to Eq. 28.

$$\mathcal{L}(\mathbf{x},\lambda_{k})=\frac{1}{2} \mathbf{x}^{t} \widetilde{\mathbf{Q}} \mathbf{x} + \mathbf{x}^{t} \widetilde{\mathbf{q}}+\lambda_{k}\mathbf{e}_{k}^{t}\left(\mathbf{x}\left[k\right]-\bar{y}_{k}\right) $$

The KKT conditions are then easily obtained by differentiating the Lagrangian and by imposing constraint satisfaction:

$$\left(\nabla_{\mathbf{x}} \mathcal{L}\left(\mathbf{x},\lambda_{k}\right)=\mathbf{Q}\mathbf{x}+\mathbf{q}+\lambda_{k}\mathbf{e}_{k}=0\right)\wedge\left(\mathbf{x}\left[k\right]=\bar{y}_{k}\right) $$

where λ _k is the Lagrangian multiplier associated to the constraint \(\mathbf {x}[k]=\bar {y}_{k}\) and e _k the vector of component δ _k (i). We can expand the \(\nabla _{\mathbf {x}} \mathcal {L}(\mathbf {x},\lambda _{k})=0\) equation to get:

It is then obvious that if we take:

$$\lambda_{k}=-\left(\mathbf{v}^{t}\mathbf{x}^{-} + \bar{y}_{k}\mathbf{Q}(k,k) + \mathbf{w}^{t}\mathbf{x}^{+} + \mathbf{q}[\!k]\right) $$

then the constraint \(\mathbf {x}[k]=\bar {y}_{k}\) is fulfilled. Substituting this into the \(\nabla _{\mathbf {x}} \mathcal {L}(\mathbf {x},\lambda _{k})=0\) equation we get:

This equation is usable but its drawback is that we have lost the symmetry of the modified Q matrix. This can be fixed by another modification of the vector term q. Exploiting the fact that by construction \(\mathbf {x}[\!k]=\bar {y}_{k}\) we have the following equivalent system:

(30)