 Research article
 Open Access
 Published:
Linear MALDIToF simultaneous spectrum deconvolution and baseline removal
BMC Bioinformatics volume 19, Article number: 123 (2018)
Abstract
Background
Thanks to a reasonable cost and simple sample preparation procedure, linear MALDIToF spectrometry is a growing technology for clinical microbiology. With appropriate spectrum databases, this technology can be used for early identification of pathogens in body fluids. However, due to the low resolution of linear MALDIToF instruments, robust and accurate peak picking remains a challenging task. In this context we propose a new peak extraction algorithm from raw spectrum. With this method the spectrum baseline and spectrum peaks are processed jointly. The approach relies on an additive model constituted by a smooth baseline part plus a sparse peak list convolved with a known peak shape. The model is then fitted under a Gaussian noise model. The proposed method is well suited to process low resolution spectra with important baseline and unresolved peaks.
Results
We developed a new peak deconvolution procedure. The paper describes the method derivation and discusses some of its interpretations. The algorithm is then described in a pseudocode form where the required optimization procedure is detailed. For synthetic data the method is compared to a more conventional approach. The new method reduces artifacts caused by the usual twosteps procedure, baseline removal then peak extraction. Finally some results on real linear MALDIToF spectra are provided.
Conclusions
We introduced a new method for peak picking, where peak deconvolution and baseline computation are performed jointly. On simulated data we showed that this global approach performs better than a classical one where baseline and peaks are processed sequentially. A dedicated experiment has been conducted on real spectra. In this study a collection of spectra of spiked proteins were acquired and then analyzed. Better performances of the proposed method, in term of accuracy and reproductibility, have been observed and validated by an extended statistical analysis.
Background
Linear matrixassisted laser desorption/ionization timeofflight mass spectrometry (MALDIToF MS) has now revolutionized identification of bacteria, yeasts and molds in clinical microbiology [1]. The technology is simple, accurate, fast, and for large laboratories less expensive than conventional methods. Despite a lower resolution than other analyzers used in modern proteomics, linear ToF are preferred in microbiology because of a better sensitivity in the 220 kDa mass range, where proteins contain phylogenic information. Moreover, the lower cost of linear instruments favored a wider adoption by health institutions. In essence, identifications are performed in minutes by simply acquiring an experimental spectrum of the whole microorganism cells and comparing the resulting peak list with a database [2, 3]. In this context we propose a new method for peak extraction especially adapted to linear MALDIToF spectra. The usual approach for MALDI mass spectra processing generally consists of chaining several procedures. Most of the times we have a smoothing step, a baseline correction step and only then the final the peak extraction [4]. The main idea of our new method is to jointly perform these steps with the aim of reducing the potential unrecoverable artifacts introduced by a sequential processing. In the next sections we briefly present the three main steps of the usual approaches. We then describe our method and its detailed derivation.
Smoothing
A popular [5] smoothing technique in the spectrometry community is the use of SavitzkyGolay 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 SavitzkyGolay for instance) of the spectrum and to extract local minima [4]. The use of second order derivative minima instead of the zerocrossing 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 sparsityinducing norms, like the l_{1} norm) [20–23]. Further generalizations can be obtained in case of blinddeconvolution [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.
Contributions
Computing the baseline correction and finding peaks are two strongly linked problems, it is thus natural to perform these two operations jointly. In this work we propose such an approach.
In the first part of the paper we describe a direct model with an additive noise where the spectrum is modelized by a smooth baseline plus a sparse peak list convolved by a given peak shape function. We describe how we chose our priors to enforce baseline smoothness and sparsity of the peak list. We then assume a Gaussian distribution for the noise. This allows us to use Euclidean distance to quantify the error between our model and the measured spectrum.
Next we show how the unknown baseline can be eliminated from the model. This manipulation leads to a modified problem very close to the classical deconvolution one. We underline this similarity and rigorously describe the two limiting cases, zero or infinite penalization for the baseline smoothness. As a by product we can interpret that our new deconvolution method is equivalent in some way to deconvolve a regularized second order derivative of the initial spectrum.
We then carefully examine the behavior of the computed baseline at the spectrum boundaries. We observed that when the smoothness penalty is too strong, the computed baseline can become overly flat. To avoid this effect a correction allowing to define baseline values at boundaries is proposed. With this modification the behavior of the baseline at the boundaries is no more affected by strong baseline smoothness penalty.
An effective optimization method to compute the solution of the deconvolution problem is exposed. This optimization algorithm is used twice in our twopasses deconvolution procedure. In the first pass a sparsity prior is used and a first optimization problem is solved to find peak centers. In the second pass the sparsity prior is replaced by the previously found peak positions. This second optimization problem is solved to compute peak height values.
Finally the new method is compared to one instance of the smoothing/baseline correction/deconvolution classical approach. The advantage of the joint baseline computation and peak deconvolution is demonstrated on synthetic data. An example on “real” data is shown and a reference to a more detailed comparison between our method and classical ones is given.
Method
Problem definition
The proposed developments are founded on a natural model for the observed spectrum y as a spiky signal x_{ p } convolved with a peak shape p superimposed onto a smooth baseline x_{ b }
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.
The problem at stake is to recover both the signal of interest x_{ p } and the baseline x_{ b }. It is a difficult task at least for three reasons.

It is underdetermined since the number of unknowns is twice the number of data.

The convolution reduces the resolution due to peak enlargement and possible overlap.

Measurement noise and possible model inadequacy induce additional uncertainties.
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.
x _{ b } smoothness
A simple way to account for smoothness of x_{ b } is to penalize its fluctuations through
where D is a finite differences matrix of size (n−1)×n (given in Appendix “Smoothness and convolution matrix”) and μ>0.
x _{ p } sparsity and positivity
In order to favor sparsity for x_{ p } (spiky property) an elasticnet penalty is introduced
The degree of sparsity is controlled by λ_{1}. The coefficient λ_{2} is generally set to zero or to a very small value. The reason why we have introduced this extra regularization is that a small positive value can sometimes improve convergence speed of the algorithm. In practice this only happens for spectra of several thousand of m/z channels and always has a limited impact on the obtained solution.
In addition, the solution also requires positivity of x_{ p }, i.e. positivity for each component x_{ p }[ k]. The l_{1} norm then simplifies
where \(\mathbb {1}_{n}\) is the n dimensional column vector with each entry to 1. This substitution turns the convex, but nonsmooth, penalty \(\mathcal {P}_{p}(\mathbf {x}_{p})\) into a smooth quadratic one. This idea has already been used in [23, 28].
Datafidelity term
A common and natural way to quantify data  model discrepancy founded on Eq. 1 relies on a squared Euclidean norm:
where the n×n band matrix L represents the convolution with the peak shape p.
Complete objective
Using Eqs. 2, 3 and 4, we get a first expression of the complete objective:
that is to say:
The minimization of this objective function
gives the desired solution \(\left (\widehat {\mathbf {x}}_{p},\widehat {\mathbf {x}}_{b}\right)\).
Elimination of x _{ b }
To find the solution of problem Eq. 7 we begin by solving it for the x_{ b } vector. This is an unconstrained quadratic problem and an analytical solution can be found. We shall first write down the gradient of \(\mathcal {J}\) with respect to x_{ b }:
Solving \(\nabla _{\mathbf {x}_{b}} \mathcal {J} =0\) yields the minimizer:
where \(\mathbf {B}_{\mu }=\mathbf {I}_{b}+\mu \mathbf {D}^{t}\mathbf {D}\). This operation is always possible since B_{ μ } is invertible (sum of the identity matrix and a semipositive matrix).
It is then possible to substitute x_{ b } for \(\widehat {\mathbf {x}}_{b}\) in objective Eq. 5 to obtain a reduced objective. After some algebra (proof is given in Appendix “Derivation of reduced objective”), we have
where C is a constant, and \(\widetilde {\mathcal {J}}(\mathbf {x}_{p})\) is defined by:
The matrix A_{ μ } is defined by:
and its interpretation is discussed in detail in “Analysis of the A_{ μ } matrix” section).
This quadratic form \(\widetilde {\mathcal {J}}(\mathbf {x}_{p})\) is our objective function and the solution \(\widehat {\mathbf {x}}_{p}\) is its minimizer subject to positivity:
Once this constrained quadratic problem is solved, we can retrieve \(\widehat {\mathbf {x}}_{b}\) using Eq. 9:
then the initial problem Eq. 7 is solved.
The gradient of \(\widetilde {\mathcal {J}}(\mathbf {x}_{p})\) is easily deduced form Eq. 10 and reads:
Interpretation of the reduced criterion
We can make a pause here to interpret Eq. 12. If there was no baseline a natural way to deconvolve the spectrum would be to solve:
If we use the positivity constraint and expand this objective function we get the following equivalent problem:
Now if we look Eq. 10 and set λ_{2}=0 we get:
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 μ.
For small μ, the A_{ μ } operator acts like a regularized second order derivation. For a sufficiently small μ (such that the spectral radius satisfies ρ(μD^{t}D)<1) a Taylor expansion gives:
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.
The regularization strength increases as the baseline parameter μ increases. At the limit \(\mu \rightarrow \infty \) (proof given in Appendix “Expression of \( \underset{\mu \to \infty }{\lim }{\mathbf{A}}_{\mu } \)”) we get:
When applied to any vector y we get
Thus the action of A_{ μ } centers the signal y by subtracting the constant \(\bar {\mathbf {y}}\mathbb {1}_{n}\) vector, where the scalar \(\bar {\mathbf {y}}\) is the mean value of the y vector’s components. The same centering holds for the peak shape function p. It follows that in this limiting case the solved problem is the usual deconvolution procedure applied on the centered spectrum. The computed baseline is simply a constant equal to the spectrum mean.
Debiasing
The l_{1} 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_{1} 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.
First stage: peak support selection
Given μ, λ_{1} and λ_{2}, we solve the minimization problem Eq. 12. The obtained solution \(\widehat {\mathbf {x}}_{p}\) is a sparse vector containing peak intensities. The intensities are biased if the λ hyperparameters are not null. However we can use \(\widehat {\mathbf {x}}_{p}\) to define the peak support Ω. In the present work the peak support is simply defined by keeping the local maxima of \(\widehat {\mathbf {x}}_{p}\).
The peak intensities are going to be corrected in the second stage, but their final positions are defined by the condition Eq. 16. More complex procedures can be used to find the peak support. One such example is the procedure presented in [22], Postprocessing and thresholding. These more refined methods can be introduced in a straightforward way in our approach by using them instead of the basic condition Eq. 16.
Second stage: peak intensity correction
In the second step, as we have peak support Ω, we do not need the l_{1} regularization anymore. Hence we simply ignore it and solve the simplified optimization problem:
Despite its more complex appearance, this problem is no more complicated than Eq. 12. Solving this problem will correct peak intensities by removing the bias induced by the previously used l_{1} penalty.
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
After boundary correction, the modified Eq. 9 is
and the modified Eq. 10 is
The solution \(\widehat {\mathbf {x}}_{p}\) is the unique minimizer of \(\widetilde {\mathcal {J}}_{2}\) subject to positivity constraint:
The gradient is obtained by a straightforward computation:
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.
Effective minimization
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 BarzilaiBorwein method described in [30]. The BarzilaiBorwein method [31] dramatically improves the classical steepest descent with Cauchy step.
For a convex quadratic form
the steepest descent direction is defined by
with
The associated Cauchy step is defined by
A straightforward calculation leads to:
The BarzilaiBorwein method uses the same steepest descent direction but replaces the Cauchy step by one (or an alternating sequence) of the two socalled BarzilaiBorwein steps. These step lengths are defined by [30]:
Despite its simple update formula the BarzilaiBorwein method works surprisingly well [32, 33]. This method has no line search and exhibits a nonmonotonic convergence. It has been shown to be globally convergent for the strictly convex quadratic case [34]. A direct extension to the constrained case is obtained by replacing the iterate x^{(k+1)} by its projection on the feasible domain:
For bound constraints x∈ [ l,u] the projection operator P is simply defined (in a component wise fashion) by
Without any line search [30] gives a counterexample where this method is not convergent. However in practice the method is generally successful, especially for badly conditioned problems where it can outperform the constrained conjugate gradients algorithm [35].
For its simplicity and good behavior in practice we have chosen to use this method. Our implementation is given in pseudocode in Algorithm 2.
Variants of this method with nonmonotone 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
Algorithm 2 stopping criterion needs to be defined. We use a rigorous one based on KarushKuhnTucker (KKT) conditions [38]. For a smooth convex optimization problem
necessary and sufficient conditions are:

1
Stationarity \(\nabla _{\mathbf {x}_{p}} \mathcal {F}+\mathbf {\lambda }^{\mathbf {u}}\mathbf {\lambda }^{\mathbf {l}}=\mathbf {0}\)

2
Primal feasibility \(\mathbf {x}_{p} \in \left [\mathbf {l},\mathbf {u}\right ]\)

3
Dual feasibility

Lower bounds λ^{l}≥0

Upper bounds λ^{u}≥0


4
Complementary slackness

Lower bounds \(\forall i,\ \mathbf {\lambda }^{\mathbf {l}}\left [i\right ](\mathbf {l}\left [i\right ]\mathbf {x}_{p}\left [i\right ])=0\)

Upper bounds \(\forall i,\ \mathbf {\lambda }^{\mathbf {u}}\left [i\right ](\mathbf {x}_{p}\left [i\right ]\mathbf {u}\left [i\right ])=0\)

Our stopping criterion reflects KKT conditions violation. Enumerating all possible cases this criterion is computed as follow
with s[ i] defined by:
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 nonmonotonic convergence behavior of the BarzilaiBorwein 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.
Results and discussion
Synthetic data
The presented algorithm has been tested on synthetic data. The synthetic spectra consist of y vectors of n=500 components. Each y is generated by summing contributions from a synthetic baseline x_{ b }, a synthetic peak list x_{ p } convolved by a Gaussian peak shape p and a Gaussian noise ε:
The detailed description of these three contributions follows.

Synthetic baseline
The analytic expression of the baseline is
$$\mathbf{x}_{b}\left[i\right] = C(s)+s\exp{\left(3\frac{i}{n}\right)}2\frac{i}{n} $$where s can take one of the following values s=−1,0,or 1.
For s=0 the baseline is a straight line, for s=−1 the baseline is concave, for s=+1 it is convex.
The constant C(s) insure a positive spectrum. Its value is C(s)=5 for s=1 and C(s)=2 otherwise.

Synthetic peak list
The analytic expression of the peak list contribution is:
$$\left(\mathbf{x}_{p}\ast p\right)[\!i]=\sum\limits_{k=1}^{n_{p}}\alpha_{k} p\left(\mu_{k},\sigma_{p},i\right) $$where we have taken n_{ p }=10 and
$$p\left(\mu_{k},\sigma_{p},i\right)= \exp{\left(\frac{1}{2}\left(\frac{i\mu_{k}}{\sigma_{p}}\right)^{2}\right)} $$These n_{ p } Gaussian peaks are defined by a constant shape factor σ_{ p }=10. The individual heights α_{ k } and centers μ_{ k } are tabulated in Table 1.

Simulated zero mean Gaussian noise
The noise ε follows a Normal law of zero mean and σ_{ noise } standard deviation:
$$\mathbf{\epsilon}[\!i] \sim \mathcal{N}\left(0,\sigma_{noise}\right) $$
Algorithm implementation
We provide a reference implementation^{Footnote 1} 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 computationpeak deconvolution”, “Comparison with the usual sequential approach” and “Real data” sections. Typical runtimes 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 computationpeak deconvolution
As explained in “Final model formulation” section, page 8, the joint baseline computation and peak deconvolution is performed in two steps:

1
a first resolution of Eq. 20 is performed with a strong λ_{1} penalty. The role of this step, is to obtain peak centers from the noisy spectrum. This set of peak centers is denoted by Ω. This step is illustrated in Figs. 4 and 5.

2
a second resolution, is then performed with a null λ_{1} penalty. This penalty is replaced by the restricted support Ω for peak centers we found in the first step. The role of this step is to correct peak heights and baseline values which were biased by the presence of the strong λ_{1} penalty of the first step. Illustrations are given in Figs. 6 and 7.
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 λ_{1} 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 λ_{1} 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. xaxis 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
In this second part we provide a comparison between our new method and a common procedure used to extract peaks. The objective is to validate the approach and to see if there are some advantages to use the joint baseline computation and peak deconvolution over more traditional approaches. Our method essentially depends on 3 hyperparameters λ_{1},λ_{2} and μ. The parameter λ_{2} is of minor importance and is constantly set to 0.1. The two other hyperparameters are:

1
the λ_{1} parameter enforcing the sparsity of the solution x_{ p },

2
the μ parameter enforcing the smoothness of the baseline x_{ b }.
In order to perform a fair comparison with another approach we tried to stick to an approach that also only uses two hyperparameters. The usual approaches in spectra processing generally chain at least two procedures [4]. In peculiar, the baseline subtraction procedure is followed by a peak picking procedure.

1
Baseline subtraction: the Statisticssensitive Nonlinear Iterative Peakclipping (SNIP) algorithm [40–42] is an efficient algorithm to compute spectrum baseline. It generally gives good results and is easy to implement. It uses only one parameter, the window width m_{SNIP}, but requires a smoothed spectrum. To smooth the spectrum we use a SavitzkyGolay filter [6, 7]. The idea is to locally fit the spectrum by a polynomial. This leastsquares fitted polynomial of degree d is computed using m_{SG} points of a running window. For each window position, the spectrum value at the window center is then replaced by the polynomial value. In order to have only one parameter, a fixed window width m_{SG}=39 is used for the smoothing step.

2
Peak picking: once the baseline has been subtracted from the original spectrum, there are a large panel of methods to extract peaks [4]. We can mention simple thresholding [42], second derivative computation using SavitzkyGolay filters, extraction from wavelet coefficients or deconvolution. For our comparison we are more interested by deconvolutionlike methods [21, 22] which are closer to our approach. By consequence we have decided to extract peaks thanks to the usual sparse deconvolution problem:
$$ \arg\min\limits_{\mathbf{x}_{p}\ge0} \frac{1}{2}\big\\mathbf{y}  \mathbf{L}\mathbf{x}_{p}\big\^{2}_{2} + \lambda_{1} \big\\mathbf{x}_{p}\_{1} + \frac{\lambda_{2}}{2} \big\\mathbf{x}_{p}\big\_{2}^{2} $$(26)As for our method, the coefficient λ_{2} is not critical and we set it to a constant value λ_{2}=0.1. The remaining λ_{1} is our second free parameter. To solve Eq. 26 we use the same BarzilaiBorwein solver, detailed in Algorithm 2 and the same twosteps procedure which consists in peak selection (high λ_{1}) and peak heights debiasing (λ_{1}=0).
To compare the two approaches we perform a systematic grid search for the two hyperparameters. The grid for our new method is given Table 2, whereas Table 3 is used by the usual method, that sequentially removes the baseline and then performs the peak deconvolution.
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 matrixvector 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:
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.
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 noisefree 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 SavitzkyGolay 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 MALDIToF 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}=1, the second run used λ_{1}=0.5.
The developed algorithm can also be used with high resolution MALDIToF 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 λ_{1}=0.5, the second run used λ_{1}=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 λ_{1} 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 BHIPRO 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.
Conclusions
We have introduced a new method for peak deconvolution. This new approach jointly performs baseline computation and peak deconvolution. The baseline equation can be solved in a closed form and is substituted into the deconvolution equation. The new deconvolution equation contains a linear operator A_{ μ } acting as a smoothed second order derivator. The problem of boundaries is exposed and tackled rigorously leading to a modified equation. The problem is efficiently solved by the projected BarzilaiBorwein algorithm. A comparison with a traditional approach relying on a sequential baseline removal and peak picking is detailed. The benefits of our new approach are put in evidence, including a better baseline approximation avoiding ascent below the large peaks, and a better reconstruction of the deconvolved peaks. Finally the method has been tested on real data.
Perspectives
There are two main directions we want to explore. The first direction is to go one step further in a joint processing approach by allowing automatic adjustment of the peak shape function. Ideally, with such an approach, the baseline, the peak shape and the deconvolved peaks would be computed jointly. The second direction would be to devise a procedure for automatic tuning of the two main hyperparameters. The resulting algorithm would be of great value to process a large number of spectra in batch mode.
Appendix
Smoothness and convolution matrix
To enforce baseline x_{ b } smoothness the first order finite differences matrix D of dimension (n−1)×n
has been introduced. For ease of reading we also provide the D^{t}D matrix expression. Its dimension is n×n and its components are:
The matrices D^{t}D and \(\mathbf {B}_{\mu }=\mathbf {I}_{b}+\mu \mathbf {D}^{t}\mathbf {D}\) are tridiagonal matrices. Hence the associated matrixvector products can be computed efficiently. The Thomas algorithm [43], which is a Gaussian elimination optimized for tridiagonal matrices, also provides an efficient way to compute the matrixvector products associated to the \(\mathbf {A}_{\mu }= \mathbf {I}_{b}\mathbf {B}_{\mu }^{1}\) operator.
In the presentation the convolution product between the peak list and peak shape function, x_{ p }∗p, is often represented by the Lx_{ p } matrixvector 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.
Let φ be a function which takes as input a vector z of size n and produces as output the scalar quantity
where q_{0} is a scalar, q is a vector of size n and Q is a strictly positivedefinite symmetric matrix of size n×n. The minimizer of φ is \(\hat {\mathbf {z}} = \mathbf {Q}^{1}\mathbf {q}\) and the minimum is
If we introduce the gradient ∇φ(z)=Qz+q and the Hessian ∇^{2}φ(z)=Q we notice that the previous relation takes the following form:
If we consider Eq. 6 as a function of x_{ b }
where
and
then application of Eq. 27 directly leads to Eq. 10:
with a constant term \(\mathbf {C}=\frac {1}{2}\mathbf {y}^{t} \mathbf {A}_{\mu } \mathbf {y}\).
Expression of \(\underset {\mu \rightarrow \infty }{\lim } \mathbf {A}_{\mu }\)
The matrix D^{t}D is real and symmetric, hence it admits an eigen decomposition
where Δ is a diagonal matrix of eigenvalues ν_{ k } and V is an orthonormal matrix. Each column k of V is an eigenvector v_{ k }.
Let’s write down two basic statements about eigenpair \(\left (\nu,\mathbf {v}\right)\) of any square matrix M:
 the first one, where I denotes the identity matrix, is:
 the second one, that requires M to be invertible (hence ν≠0), is:
Now observe that like \(\mathbf {I}_{b}+\mu \mathbf {D}^{t}\mathbf {D}\) is a strictly diagonally dominant matrix it is invertible [44]. Using the two previously enumerated relations it is easy to show that this matrix admits the following eigen decomposition:
By construction (first order finite differences) the product of the matrix D with a constant vector is the null vector. It follows that we have \(\mathbf {D}^{t}\mathbf {D}\,\mathbf {v}_{1}=\mathbf {0}\) where \(\mathbf {v}_{1}=\frac {1}{\sqrt {n}}\mathbb {1}_{n}\). Thus v_{1} is the normalized eigen vector associated to the ν_{1}=0 eigenvalue. Moreover, as D^{t}D is an unreduced^{Footnote 2} real symmetric tridiagonal matrix, all its eigenvalues are simples ([45], Lemma 7.7.1). By consequence we can write:
with ν_{ k }≠0 for all k≥2. It is then trivial to compute the limit:
which is a square matrix of dimensions n×n with all entries equal to 1.
From the definition of A_{ μ } in Eq. 11 we get the expected result:
with δ_{i,j} equals to 1 if i=j and 0 otherwise.
Boundary correction
Our boundary correction approach exploits the opportunity to convert an equality constrained problem
into a unconstrained modified problem
Let begin by writing the Lagrangian [38] associated to Eq. 28.
The KKT conditions are then easily obtained by differentiating the Lagrangian and by imposing constraint satisfaction:
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:
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:
By identification we can write down the announced result:
which is an unconstrained system for which the minimizer satisfies \(\mathbf {x}\left [k\right ]=\bar {y}_{k}\). For presentation clarity we only have considered one constraint, but this scheme can be used sequentially for any set of constraints
Now we can see how to use this technique to our deconvolution problem. We must go back to Eq. 8:
By identification the role of the matrix Q is played by the B_{ μ } matrix. Furthermore all the modifications done on the vector q can be transposed to modifications associated to the −y vector.
As an illustrative example, we consider the most usual case that consists of imposing the left \(\bar {y}_{1}\) and right \(\bar {y}_{n}\) baseline values at boundaries. Direct application of Eq. 30 provides the result. The modified B_{ μ } matrix is:
The modified y vector is:
With these modifications, the modified Eq. 9 that allows to retrieve the baseline x_{ b } from the deconvolved peaks x_{ p } is now:
The modified A_{ μ } operator Eq. 11 is computed as before
The final quadratic form Eq. 13 must be replaced by:
Notes
 1.
 2.
No offdiagonal element is null
Abbreviations
 BB:

BarzilaiBorwein
 BHIPRO:

Bayesian hierarchical inversion for mass spectrometry. Application to discovery and validation of new PROtein biomarkers
 CWT:

Continuous wavelet transform
 DWT:

Discrete wavelet transform
 KKT:

KarushKuhnTucker
 LC:

Liquid chromatography
 MALDIToF:

Matrixassisted laser desorption/ionizationtime of flight
 MS:

Mass spectrometry
 SG:

Savitzkygolay
 SNIP:

Statisticssensitive nonlinear iterative peakclipping
 UDWT:

Undecimated wavelet transform
References
 1
van Belkum A, Welker M, Pincus D, Charrier JP, Girard V. Matrixassisted laser desorption ionization timeofflight mass spectrometry in clinical microbiology: What are the current issues?Ann Lab Med. 2017; 37(6):475.
 2
Ilina EN, Borovskaya AD, Malakhova MM, Vereshchagin VA, Kubanova AA, Kruglov AN, Svistunova TS, Gazarian AO, Maier T, Kostrzewa M, Govorun VM. Direct bacterial profiling by matrixassisted laser desorptionionization timeofflight mass spectrometry for identification of pathogenic neisseria. J Mol Diagn. 2009; 11(1):75–86.
 3
Singhal N, Kumar M, Kanaujia PK, Virdi JS. Malditof mass spectrometry: an emerging technology for microbial identification and diagnosis. Front Microbiol. 2015; 6:791. https://doi.org/10.3389/fmicb.2015.00791.
 4
Yang C, He Z, Yu W. Comparison of public peak detection algorithms for maldi mass spectrometry data analysis. BMC Bioinformatics. 2009; 10(1):4.
 5
Larive CK, Sweedler JV. Celebrating the 75th anniversary of the acs division of analytical chemistry: A special collection of the most highly cited analytical chemistry papers published between 1938 and 2012. Anal Chem. 2013; 85(9):4201–2.
 6
Savitzky A, Golay MJ. Smoothing and differentiation of data by simplified least squares procedures. Anal Chem. 1964; 36(8):1627–39.
 7
Gorry PA. General leastsquares smoothing and differentiation by the convolution (savitzkygolay) method. Anal Chem. 1990; 62(6):570–3.
 8
Browne M, Mayer N, Cutmore TR. A multiscale polynomial filter for adaptive smoothing. Digit Signal Process. 2007; 17(1):69–75.
 9
Barak P. Smoothing and differentiation by an adaptivedegree polynomial filter. Anal Chem. 1995; 67(17):2758–62.
 10
Shensa MJ. The discrete wavelet transform: wedding the a trous and mallat algorithms. IEEE Trans Signal Process. 1992; 40(10):2464–82.
 11
Nason GP, Silverman BW. The stationary wavelet transform and some statistical applications In: Antoniadis A, Oppenheim G, editors. Wavelets and Statistics. New York: Springer: 1995. p. 281–99. https://doi.org/10.1007/9781461225447_17.
 12
Coombes KR, Tsavachidis S, Morris JS, Baggerly KA, Kuerer HM. Improved peak detection and quantification of mass spectrometry data acquired from surfaceenhanced laser desorption and ionization by denoising spectra with the undecimated discrete wavelet transform. Proteomics. 2005; 5:4107–17.
 13
Antoniadis A, Bigot J, LambertLacroix S. Peaks detection and alignment for mass spectrometry data. J Société Française Stat. 2010; 151(1):17–37.
 14
PerezPueyo R, Soneira MJ, RuizMoreno S. Morphologybased automated baseline removal for raman spectra of artistic pigments. Appl Spectrosc. 2010; 64(6):595–600.
 15
Morháč M, Matoušek V. Peak clipping algorithms for background estimation in spectroscopic data. Appl Spectrosc. 2008; 62(1):91–106.
 16
Mazet V, Carteret C, Brie D, Idier J, Humbert B. Background removal from spectra by designing and minimising a nonquadratic cost function. Chemometr Intell Lab Syst. 2005; 76(2):121–33.
 17
Ruckstuhl AF, Jacobson MP, Field RW, Dodd JA. Baseline subtraction using robust local regression estimation. J Quant Spectrosc Radiative Transf. 2001; 68(2):179–93.
 18
Li Z, Zhan DJ, Wang JJ, Huang J, Xu QS, Zhang ZM, Zheng YB, Liang YZ, Wang H. Morphological weighted penalized least squares for background correction. Analyst. 2013; 138(16):4483–92.
 19
Dubrovkin J. Evaluation of the peak location uncertainty in secondorder derivative spectra. Case study: symmetrical lines. J Emerg Technol Comput Appl Sci. 2014; 3:9.
 20
MohammadDjafari A, Giovannelli JF, Demoment G, Idier J. Regularization, maximum entropy and probabilistic methods in mass spectrometry data processing problems. Int J Mass Spectrom. 2002; 215(1):175–93.
 21
Renard BY, Kirchner M, Steen H, Steen JA, Hamprecht FA. Nitpick: peak identification for mass spectrometry data. BMC Bioinformatics. 2008; 9(1):355.
 22
Slawski M, Hussong R, Tholey A, Jakoby T, Gregorius B, Hildebrandt A, Hein M. Isotope pattern deconvolution for peptide mass spectrometry by nonnegative least squares/least absolute deviation template matching. BMC Bioinformatics. 2012; 13(1):291.
 23
Giovannelli JF, Coulais A. Positive deconvolution for superimposed extended source and point sources. Astron Astrophys. 2005; 439:401–12.
 24
Hirsch M, Schölkopf B, Habeck M. A blind deconvolution approach for improving the resolution of cryoem density maps. J Comput Biol. 2011; 18(3):335–46.
 25
Lange E, Gröpl C, Reinert K, Kohlbacher O, Hildebrandt A. HighAccuracy Peak Picking Of Proteomics Data Using Wavelet Techniques. In: Biocomputing 2006. World Scientific: 2012. p. 243–54. https://www.worldscientific.com/doi/abs/10.1142/9789812701626_0023.
 26
Du P, Kibbe WA, Lin SM. Improved peak detection in mass spectrum by incorporating continuous wavelet transformbased pattern matching. Bioinformatics. 2006; 22(17):2059–65.
 27
Mallat S, Zhong S. Characterization of signals from multiscale edges. IEEE Trans Pattern Anal Mach Intell. 1992; 7:710–32.
 28
Figueiredo MA, Nowak RD, Wright SJ. Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems. IEEE J Sel Top Sign Process. 2007; 1(4):586–97.
 29
Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B Methodol. 1996; 58(1):267–88. http://www.jstor.org/stable/2346178.
 30
Dai YH, Fletcher R. Projected barzilaiborwein methods for largescale boxconstrained quadratic programming. Numer Math. 2005; 100(1):21–47.
 31
Barzilai J, Borwein JM. Twopoint step size gradient methods. IMA J Numer Anal. 1988; 8(1):141–8.
 32
Fletcher R. On the barzilaiborwein method In: Qi L, Teo K, Yang X, editors. Optimization and Control with Applications. Boston: Springer US: 2005. p. 235–56.
 33
Raydan M. The barzilai and borwein gradient method for the large scale unconstrained minimization problem. SIAM J Optim. 1997; 7(1):26–33.
 34
Raydan M. On the barzilai and borwein choice of steplength for the gradient method. IMA J Numer Anal. 1993; 13(3):321–6.
 35
Moré JJ, Toraldo G. On the solution of large quadratic programming problems with bound constraints. SIAM J Optim. 1991; 1(1):93–113.
 36
Kim D, Sra S, Dhillon IS. A nonmonotonic method for largescale nonnegative least squares. Optim Methods Softw (OMS). 2012.
 37
Birgin EG, Martínez JM, Raydan M. Nonmonotone spectral projected gradient methods on convex sets. SIAM J Optim. 2000; 10(4):1196–211.
 38
Boyd S, Vandenberghe L. Convex optimization.Cambridge university press; 2009.
 39
Mercier C, Klich A, Truntzer C, Picaud V, Giovannelli JF, Ducoroy P, Grangeat P, MaucortBoulch D, Roy P. Variance component analysis to assess protein quantification in biomarker discovery. application to MALDITOF mass spectrometry. Biom J. 2017; 60(2):262–74. https://doi.org/10.1002/bimj.201600198.
 40
Ryan C, Clayton E, Griffin W, Sie S, Cousens D. Snip, a statisticssensitive background treatment for the quantitative analysis of pixe spectra in geoscience applications. Nucl Inst Methods Phys Res Sec B: Beam Interactions with Mater Atoms. 1988; 34(3):396–402.
 41
Morháč M, Kliman J, Matoušek V, Veselskỳ M, Turzo I. Background elimination methods for multidimensional coincidence γray spectra. Nucl Inst Methods Phys Res Sec A: Accelerators, Spectrometers, Detectors and Assoc Equip. 1997; 401(1):113–32.
 42
Gibb S, Strimmer K. Maldiquant: a versatile R package for the analysis of mass spectrometry data. Bioinformatics. 2012; 28(17):2270–1.
 43
Conte SD, Boor CWD. Elementary numerical analysis: an algorithmic approach: McGrawHill High Educ.1980.
 44
Varga RS. Geršgorin and his circles: Springer Berlin Heidelberg; 2004. https://doi.org/10.1007/9783642177989.
 45
Parlett BN. The symmetric eigenvalue problem. SIAM. 1980; 7:134.
Acknowledgment
We thank Bruno Lacroix (bioMérieux) and Pierre Mahe (bioMérieux) for their contribution during the ANR Grant submission and execution. We thank Laurent Gerfault (CEA) for his advises for the development of the BHIPRO MALDI model, data processing and data analysis. We thank Amna Klich, Delphine MaucortBoulch, Pascal Roy (Biostatistique, Hospices Civils de Lyon) and Patrick Ducoroy (CLIPP) for their advises and comments during the development and tests of the method.
Funding
The BHIPRO project has been partially funded by the Agence Nationale de la Recherche under grant ANR 2010 BLAN 0313 and by Commisariat à l’Energie Atomique et aux Energies Alternatives (CEA) for CEA authors contribution (VP, PG).
Availability of data and materials
The datasets used to generate the figures and a basic implementation of the proposed algorithm are avalaible under the https://github.com/vincentpicaud/Joint_Baseline_PeakDeconvGitHub repository.
Author information
Affiliations
Contributions
VP initial idea of the method, theoretical developments and C++ implementation. The contributions of JFG and AG regard (a) the design of the criterion and its optimization and (b) the writing of the paper. CT read and approved the final manuscript and proposed valuable comments about preprocessing of spectra in the context of mass spectrometry. CM participated in the problem definition and evaluated the methods on real spectra. JPC provided input in the MALDIToF field and the result interpretation. He acquired and provided real MALDIToF spectra to test the algorithm. PG was the BHIPRO project manager. He has coordinated the conception and design of the processing algorithms, in particular the design of the acquisition chain MALDI and MRM models and inversion algorithms, and the interpretation of the data. He has revised the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable.
Competing interests
JPC is employed by bioMérieux. The other authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Picaud, V., Giovannelli, J., Truntzer, C. et al. Linear MALDIToF simultaneous spectrum deconvolution and baseline removal. BMC Bioinformatics 19, 123 (2018). https://doi.org/10.1186/s1285901821163
Received:
Accepted:
Published:
Keywords
 Mass spectrometry
 Peak picking
 Deconvolution
 Baseline