Abstract
Background
Highthroughput proteomics techniques, such as mass spectrometry (MS)based approaches, produce very highdimensional datasets. In a clinical setting one is often interested in how mass spectra differ between patients of different classes, for example spectra from healthy patients vs. spectra from patients having a particular disease. Machine learning algorithms are needed to (a) identify these discriminating features and (b) classify unknown spectra based on this feature set. Since the acquired data is usually noisy, the algorithms should be robust against noise and outliers, while the identified feature set should be as small as possible.
Results
We present a new algorithm, Sparse Proteomics Analysis (SPA), based on the theory of compressed sensing that allows us to identify a minimal discriminating set of features from mass spectrometry datasets. We show (1) how our method performs on artificial and realworld datasets, (2) that its performance is competitive with standard (and widely used) algorithms for analyzing proteomics data, and (3) that it is robust against random and systematic noise. We further demonstrate the applicability of our algorithm to two previously published clinical datasets.
Availability of data and methods
The method sourcecode can be downloaded from our homepage: http://software.medicalbioinformatics.de. The used data will be made available through BioMed Central.
List of abbreviations
AIC  Akaike Information Criterion; BIC  Bayesian Information Criterion; CS  Compressed Sensing; MALDITOF: MatrixAssisted Laser Desorption Ionization TimeOfFlight; ML  Machine Learning; MS  Mass Spectrometry; SPA  Sparse Proteomics Analysis; SVM  Support Vector Machine; TP  True Positive; TN  True Negative; FP  False Positive; FN  False Negative
1 Introduction
During the last decade, highthroughput assays systems^{1}^{1}1Assays, e.g. immunoassays, are used in molecular diagnostics to detect concentrations of specific molecules even in low concentrations from a biological sample, such as blood[1]. for measuring a variety of different biological sources have become standard in modern laboratories. This allows for the quick and cheap creation of very large datasets which characterize for example the status of a cell by its billions of constituents, e.g. nucleotides, RNAs, contained proteins, or metabolites. Ideally, analyzing these massive datasets leads to a better understanding of the underlying biological processes. Especially in the context of characterizing—and ultimately understanding—diseases, a first step is often to find significant differences in the data between samples from healthy and diseased individuals. There are many successful examples where this approach based on omics data (e.g., genomics, proteomics, or metabolomics) led to the identification of biological markers, enabling a new type of molecular diagnostics. We call a collection of biological markers that represents the differences on the data level a disease fingerprint.
Many diseaserelevant mechanisms are controlled by proteins (e.g. hormones), which can be detected in biological samples (blood, urine, etc.) using mass spectrometry (MS). This technique allows (potentially) for monitoring the entire set of proteins—the socalled proteome—in a given sample. Due to its wide availability in hospitals, MSbased proteomics can bring the next wave of progress in diagnostics, since even subtle changes in the proteome can be detected and linked to disease onset and progression [2, 3, 4, 5].
Disease Fingerprints: The main idea of the identification of disease fingerprints using MSbased proteomics is sketched in Fig. 1:
(a) A mass spectrum is generated reflecting the constitution of a given (blood)sample with respect to contained molecules. (b) Based on mass spectra from two sample groups (representing a healthy control group and a group having a particular disease) differences are detected. This set of differences precisely corresponds to a disease fingerprint, since it represents a trace caused by a particular disease in the proteome. Several studies have shown that this approach works well in practice and found differences do indeed reflect correlations between changes in the mass spectrum, the proteome, and phenotypic changes ([6, 7, 8, 9, 10]). Panels of proteomic markers (fingerprints) have been shown to be more sensitive and specific than conventionally biomarker approaches [3], for example when diagnosing cancer [11, 12, 13]. However, a single proteomics dataset can contain tens of millions of signals which is many orders of magnitudes larger than the number of available observations in a typical study.
Our ultimate goal is therefore to build a library of proteomics disease fingerprints which are extracted from highthroughput MS experiments. These would enable to diagnose diseases based on their proteomic fingerprints—just by analyzing an individual’s proteome. Ideally, these fingerprints are of lowcomplexity allowing easy interpretation by experts, e.g. medical doctors, and the implementation of medical assays for routine diagnostics, e.g. in an hospital environment. Clearly, the less components an assay is composed of, the easier it is to implement and interpret.
Thus, a fingerprint should only consist of a minimal collection of proteins specific for a particular disease and should be robust against noisy measurements. On the other hand, the acquired data from the highthroughput experiments is very highdimensional and contains large amounts of random and systematic noise which makes an automatic analysis of mass spectra a very challenging task. Hence, the discovery of biomarkers is still a widely open research topic and there are several analytic problems that hinder reproduction of results (see [14] for example).
Despite these challenges there is indeed hope that these disease specific, lowcomplexity fingerprints exist: It has been shown for several cancer types that a small numbers of genes and proteins can be identified that serve as biomarkers (e.g. for lung cancer [15], breast cancer [16] or pancreas cancer [17]). This means that only a few signals in a mass spectrum can be used to derive a sparse classifier.
MS1 Data: In this work we consider mass spectrometry (MS) data acquired from a standard MALDITOF instrument because it is easy to obtain using comparatively cheap MSinstruments which are widely available, e.g. in hospitals. Opposed to other approaches such as tandem mass spectrometry (MS/MS), we directly work on the raw data acquired in profile mode and do not aim for identification. Thus, each mass spectrum (sample) always has the same number of dimensions (number of entries). ^{2}^{2}2The datasets used in this paper contain dimensions in each MS1 spectrum but our approach is not limited by that. Recall, that the entries in a mass spectrum are a weightordered list of ioncounts of the respective ionmasses. (See also Fig. 1.)
One of the reasons for this is that standard approaches for MS data analysis usually convert the MS data to peak lists as a first step and work on the converted data. However, signals can be missed by this conversion step due to noise or missing values in the raw data which hinders peak detection. Opposed to this, our approach does not rely on any peak identification but works on the raw data. This allows for a more robust analysis in presence of noise which is a typical challenge in MS data analysis.
1.1 Problem Definition
In this article, we will focus on the following problem setting:
We assume that we are given data of mass spectra derived from biological samples (e.g. from blood of individual patients) in form of pairs . Here, represents the mass spectrum of the th sample (e.g. the th patient) and its respective class, e.g., healthy or diseased. Thus, each (representing an individual mass spectrum) contains entries. The goal is to identify a (small) set of features, i.e. indices in the mass spectrum, separating these two classes. Thus, a feature represents a specific position (or mass) in a mass spectrum in which the two groups (e.g. healthy vs. diseased) differ. This corresponds to the well known problem of feature selection^{3}^{3}3In feature selection, one is interested in identifying relevant dimensions of the data (features) which can be used to distinguish between two (or more) classes within a dataset. and leads to a potential disease fingerprint for the given data.
Mathematically, this can be formulated as the identification of a feature vector such that^{4}^{4}4Here, denotes the sign function, i.e., if and if .
(1) 
with a linear decision function .
From a geometric perspective, this means that the hyperplane with normal vector
appropriately separates the datapoints of the respective classes. This means that can be used as a linear classifier where each entry of corresponds to a specific position in a spectrum and the nonzero entries (which we call features) indicate their significance. Our goal is therefore to learn a sparse for which Eqn. 1 holds. As a particular consequence, a classifier based on such will yield good prediction accuracy.In most realistic scenarios for feature selection, unfortunately, the number of features is much larger than available samples () and the data suffers from noisy measurements. For these reasons, the number of feasible classifiers can become extremely large, so that the problem of overfitting can occur. In order to allow interpretability and generalization of the classifier, it is in fact inevitable to restrict the solution space for . In this paper, we focus on very sparse^{5}^{5}5We call a vector sparse if the number of nonzero entries is small. vectors satisfying (1), which precisely reflects our wish for a minimal disease fingerprint.
At this point, it should be emphasized that (1) does not need to hold for all samples but rather for most of them. Allowing for such a small “mismatch” in the model, we incorporate the crucial fact that a simple binary output model, such as (1), might describe the disease label only with high accuracy but not necessarily exactly. In turn, this asks for a certain robustness of the used method against wrong predictions with regard to (1).
We will approach this challenge by formulating the feature selection problem as a constrained (or regularized) optimization problem:
(2) 
where is a loss (error) function, is a regularization (cost) function that encourages a particular structure of (e.g., sparsity), and the parameter controls the degree of model complexity. Given any potential feature vector and the (true) output label
, the loss function
measures the discrepancy between the actual and the desired prediction.As already pointed out, we are particularly interested in a method that produces optimal and robust solutions in the following situation:

The input data are noisy,

the number of data dimensions is large (typically: ),

the number of samples is relatively small (typically: ), and

the set of highlyrelevant features is small (i.e., a minimal disease fingerprint indeed exists), which corresponds to a small number of nonzero elements in (typically: ).
On the contrary, we are not mainly interested in the methods’ overall classification performance. Measures of classification performance such as accuracy are indicators whether a learned classifier accurately separates the data into classes. In our case, we assume that the data can be characterized well by a sparse classifier whose nonzero entries are those used for classification and are therefore of medical relevance. That means, if is sparse and leads to good classification accuracy then only a few entries contribute and medical interpretation becomes feasible. However, if there does not exists a sparse such that Eqn. 1 holds, there is strong evidence that no sparse (simple) characterization is possible. This indicates that the underlying biological mechanisms are too complex to be captured by a sparse (simple) model. If this is the case, every sparsityencouraging method will fail, meaning that a sparse classifier will always give poor classification. As a consequence, an important assumption of this work is that a sparse (groundtruth) exists.
As we will see later it is often possible to find nonsparse classifiers which achieve better classification accuracy. This might be favorable in some situations in which the main focus is indeed on overall classification accuracy. However, in these situations overfitting becomes an issue and the identification of interpretable, highlydiscriminative features might be extremely difficult. In the context of MSdata analysis such a classifier would be especially hard to interpret because of the very high dimensionality of the data.
1.2 State of the Art in Sparse Feature Selection
There are numerous approaches for feature selection which mainly fall into three categories:

Filters
: Using some score or correlation function (e.g., based on Fisher’s, ttest, information theoretic criteria) evaluating the importance of each feature in a
univariate way and taking the toprated features. 
Wrappers
: Using machinelearning algorithms to evaluate and choose features using some search strategy (e.g. simulated annealing or genetic algorithms).

Embedded methods
: Selecting variables by directly optimizing an objective function (usually in a multivariate way) with respect to: goodnessoffit and (optionally) number of features. This could be achieved with algorithms like leastsquare regression, support vector machines (SVM), or decision trees.
In this paper, we will mainly focus on embedded methods. Regarding this category, the literature contains several wellknown options for choosing combinations of loss and regularization functions (cf. (2)), some of which are exemplarily listed in Table 1.
Name  Loss function ()  Regularizer () 

AIC/BIC  
Lasso  
Elastic Net  +  
Regularized Least Absolute  
Deviations Regression  
Classic SVM  ^{*}  
SVM  ^{*}  
Logistic Regression  
^{*}This is the so called Hinge loss. 
Prominent options for choosing loss function and regularizer in feature extraction algorithms. The
 and norm of a vector are defined by and , respectively. The “norm” , simply counts the number of nonzero entries of .Different combinations can influence the results dramatically: Fig. 2 demonstrates the effect of sparsity by comparing a  and regularized version.
In this example, a proteomics dataset was created that contains three discriminant features between the two subgroups. It can be easily seen how the results differ: While the based result is optimized for selecting only a few features, the variant selects much more features which in turn results in a better fit of the observation model. In this paper, we are interested in developing a method that selects as few features as possible while achieving the best possible fit under this constraint. This is in contrast to methods that aim at only achieving the best possible fit. A lowcomplexity model is of particular interest in biological applications because each selected feature is usually analyzed in subsequent experiments, which creates additional costs.
Various approaches can be used to assess the outcome of a feature selection method, when appropriate training and test data are available. We will use the following three measures of quality: (i) correctness of the selected features, (ii) size of the selected feature set, (iii) performance of classifying an unknown test set (specificity, sensitivity, accuracy). Obviously, (i) can only be used if the correct features are known, which is the case in our benchmark datasets (for more details see Section 4).
1.3 Contribution
As already mentioned above, the major challenge of sparse feature extraction is to robustly identify a small set of variables (nonzero components of ) that can be used to accurately classify unknown proteomics data (e.g. healthy or diseased) according to (1). This paper introduces Sparse Proteomics Analysis (SPA), a novel framework for feature selection and classification. The key step of our method is based on 1bit compressed sensing (cf. Section 2) and solves the following optimization problem:^{6}^{6}6Here, again denotes the Euclidean scalar product.
(3) 
where the regularization is now defined by two inequality constraints on the feature vector .^{7}^{7}7For the sake of convenience, we formulate our algorithm as in (3), but with some slight modifications, it could be equivalently stated in the form of (2). The above approach is motivated by the general theory of compressed sensing, which was originally introduced by Donoho as well as by Candès, Romberg, and Tao (cf. [18, 19, 20]) and provides a modern framework for efficiently acquiring and processing highdimensional (nearly) sparse signals (for more details see Section 2).
We shall verify the competitiveness of our method by applying it to several synthetic and realworld datasets and comparing the results to those of other widelyused algorithms in this field. Although the core of the algorithm (3) is surprisingly simple, we will observe that SPA (including pre and postprocessing steps) finds optimal feature vectors which are extremely sparse, allow for highly accurate classification, and are robust against noise. In particular, for “verysparse” situations, it even turns out that SPA outperforms the standard methods listed in Table 1.
Note that computational solutions to (2) or (3) are usually based on solving a convex program by standard optimization techniques, such as interior point methods. However, these methods sometimes scale poorly with increasing number of samples and data dimension , as it is typically the case for omics data analysis. Several strategies have been proposed in the literature to speed up the calculations, e.g., using stochastic decent ([21, 22, 23, 24, 25]). In this article, we shall not focus on such computational issues but rather on providing a novel way of formalizing and solving the feature selection problem, namely in the context of compressed sensing.
Apart from the specific approach of (3), it is a general concern of this work to promote the benefit of sparse embedded methods. In contrast to classical (univariate) approaches, such as statistical tests, the process of variable selection takes place in an automatic fashion here. In this way, a costly preprocessing (e.g., peak detection) as well as subsequent feature assessments can be avoided as much as possible. Especially in a situation where only a very few samples are available, those additional steps may cause further instability and their success strongly relies on the specific data structure. In fact, it was already succinctly emphasized by Vapnik in [26, p. 12] that
“If you possess a restricted amount of information for solving some problem, try to solve the problem directly and never solve the more general problem as an intermediate step. It is possible that the available information is sufficient for a direct solution but is insufficient for solving a more general intermediate problem.”
This fundamental principle is precisely reflected by our viewpoint, which only makes a few (generic) assumptions on the underlying data model. Finally, we would like to mention that recently, rigorous theoretical guarantees for sparse feature selection from MS data were shown in [27]. Using the novel idea of optimal problem representations, the mathematical framework of [27] even goes beyond the binary output scheme of (1) and allows for a unified treatment of general observation and data models.
1.4 Outline of this Paper
We start by shortly reviewing the background of compressed sensing in Section 2, and then describe our novel feature selection approach SPA in detail (Section 3). Finally, we present several benchmark results in Sections 4 and 5 for simulated and real datasets and compare them to current stateoftheart algorithms.
2 Background: Compressed Sensing
2.1 Compressed SensingBased Data Analysis
In its most simple form, compressed sensing (CS) studies the recovery of an unknown vector from linear measurements . Here, is an matrix and the entries of contain the measurements. The major challenge is now to design the measurement process in such a way that the number of measurements is as small as possible and, at the same time, is still (uniquely) recoverable from . Thus, we are asking for the maximal compressibility of by linear measurements.
Obviously, when , we require some additional information to obtain a unique solution of . The prior information on which is studied in compressed sensing is the assumption of sparsity, i.e., most coefficients of are assumed to be zero, or at least very small. One naive approach to incorporate this additional property is to search for the sparsest solution of :^{8}^{8}8Here, simply counts the number of nonzero elements of .
(4) 
Unfortunately, this problem is nonconvex and cannot be efficiently solved in general. Therefore, one usually replaces (4) by its convex relaxation, which is also known as the basis pursuit ([28]):
(5) 
One of the first key results in compressed sensing states that, if is chosen randomly, e.g., with independent and identically distributed Gaussian entries, and
, then (with “high probability”) every
sparse vector (i.e., ) can be uniquely recovered from (5). The most surprising fact is that the number of required measurements only logarithmically depends on the (possibly large) dimension of the ambient space. Hence, random measurement processes indeed allow for a very strong compression of sparse vectors (see also [18, 19, 20] for more details).In order to consider more complicated situations, the stability and robustness of the basis pursuit algorithm was extensively studied. Various theoretical results and numerical experiments show that this algorithmic approach can also be applied for the stable recovery of vectors which are only nearly sparse, as well as to noisy measurements of the form . To obtain a robust version of (5), one may replace its equality constraint by for some appropriate noise level . Not very surprisingly, this approach is also closely related to the Lasso introduced by Tibshirani in [29] (see also (2) and Table 1).
2.2 1Bit Compressed Sensing
In many practical scenarios, especially when working with computers, there is no way to represent real numbers exactly. Thus, it is reasonable to assume that the measurement vector is acquired in a quantized (and therefore nonlinear) fashion. The most extreme form directly leads to 1bit measurements, i.e., only the signs of are known:
(6) 
where are the rows of the measurement matrix . As in classical compressed sensing, we are asking for an appropriate recovery of from (6) using as few measurements as possible. This challenge was originally considered in [30] as 1bit compressed sensing, and has been extensively studied in [31, 32].
A surprisingly simple convex recovery approach was proposed by Plan and Vershynin in [32]:
(7) 
where denotes the sparsitycontrolling parameter. To get some intuition, we first note that we have if and only if holds. Hence, maximizing the sum in (7) will ensure the consistency of many measurements , according to (6). However, the total consistency is not enforced so that (7) indeed allows for noisy inputs that do not satisfy (6). On the other hand, the constraint of (7) promotes sparsity of the final outcome. To see this, we may consider the set and observe that (cf. [32, Sec. III])^{9}^{9}9Here, denotes the convex hull of the set .
This means that (7) optimizes over a convex relaxation of the set which contains all sparse vectors in the unit ball. For more details, see also [32]. The main statement of [32] proves that the robust 1bit compressed sensing algorithm (7) indeed allows for an appropriate recovery of sparse vectors, using only measurements. Moreover, it is surprisingly robust against several types of noise, including (random) bitflips of the labels.
2.3 Why Compressed Sensing?
At a first sight, the main challenges of compressed sensing and machine learning (ML) seem to be very different. In compressed sensing, we intend to design a measurement process in order to compress a vector , whereas in machine learning, the training data is already contained in the rows of and we are rather willing to explain the observations by some appropriate vector . However, in both areas we are asking for a (sparse) recovery from a certain type of measurement. Indeed, a linear regression in ML exactly corresponds to classical CS model (see Subsection 2.1), and a classification problem is actually equivalent to 1bit CS (see Subsection 2.2).
Therefore, it is not very surprising that the applied algorithms for compressed sensing and machine learning resemble each other, and that theoretical results in both fields rely on the same mathematical foundations (concentration of measure, convex geometry, etc.). Unfortunately, both communities only rarely interacted with each other. In this paper, we would like to emphasize the viewpoint of compressed sensing, in particular, because it is still not very common for the classification tasks that we deal with.
With the recent progress in compressed sensing and related areas as lowrank matrix recovery or quantized CS, also new algorithms like nuclear norm minimization or 1bit CS have been proposed. Although these methods are typically motivated by theoretical studies, they perform also very well for realworld data. In general, we believe that these alternative perspectives allow for deeper theoretical insights, finally leading to the improvement of the classical (based) tools from machine learning.
For an extensive introduction to compressed sensing, we refer to [34, 35]
. As we already mentioned above, comparing this text to literature from statistical learning theory (see
[36] for example), the reader will quickly notice many interesting connections between both fields.3 Sparse Proteomics Analysis (SPA)
In this section, we present the details of our novel framework which is based on the ideas of 1bit compressed sensing presented in the previous section. The first part provides a mathematical formulation of the feature selection problem as well as a brief overview of the steps that are performed in SPA. The rest of this section is then devoted to a detailed description and discussion of the single steps.
3.1 Setting and Overview
As already mentioned in the introduction, we assume that our learning process is supervised, i.e., we know which spectrum belongs to the class of healthy () and diseased () samples in advance. If the data vectors , are mass spectra, the indices of correspond to the values^{10}^{10}10 is the unit for the masstocharge ratio. and its entries represent the intensities. The nonzero entries of the feature vector describe the location of the disease fingerprints and its respective values the significance of these features.
In the setting of classical learning theory, we are asking for a hyperplane which correctly separates most of the data points labeled by . More precisely, this means^{11}^{11}11Compared to Section 2, we are now using the standard notations from learning theory. In particular, the measurement vectors are denoted by (instead of ) and the feature vector is (instead of ).
(8) 
Equivalently, we can view (8) as a problem from 1bit compressed sensing (cf. Section 2.3), i.e., we have acquired noisy 1bit measurements and are now looking for a sparse recovery.
In the development of SPA, we have primarily focused on the latter perspective, and therefore, the 1bit recovery program (7) forms the key step of our algorithm:
Algorithm 1 (SPA at a Glance).
Input: Raw data samples Output: Sparse feature vector Preprocessing: Normalize data to make the spectra comparable. Perform smoothing by a convolution with Gaussian density. Standardize data. Sparse Feature Selection: Perform 1bit CS optimization (7) to find feature vector . Postprocessing: Detect the connected components of to obtain a sparsified version . (Optional) Reduce dimension by projecting data onto the feature space.3.2 Algorithmic Details
In the following, we are going to specify and discuss the single steps of Algorithm 1.
Step 1: Normalization of the Data
This step heavily depends on the underlying acquisition method of the data. Every spectrum is normalized by a certain scaling factor , i.e., for . The individual scalars should be chosen such that the resulting data vectors are “comparable.”
For example, when we assume that the data are acquired by MALDITOFMS as described in Fig. 1, it seems to be quite natural to normalize them by the total ion count. Mathematically, this means that we would divide every spectrum by its norm, i.e., we choose .
Step 2: Smoothing by Gaussian Density
We already pointed out that one major challenge is the strong noise within the raw data. Therefore, it is crucial to perform some noise reduction before trying to extract features. For this purpose, we suggest a simple smoothing strategy by a Gaussian density:
Let denote the (centered) Gaussian density function
with fixed standard deviation
, that is,The smoothed spectra are then obtained by a discrete convolution
(9) 
Using the fast Fourier transform (FFT), this computation can be performed quickly with
operations. In a very simplified scenario, a spectrum can be written as the sum of Gaussianshaped peaks plus some baseline noise in each mass channel. Since the convolution of two Gaussian densities is again Gaussian, the original (local) structure of the spectra is essentially preserved in , whereas the noise of is significantly reduced. Note that the deviation serves as a tuning parameter of the algorithm. A good choice of clearly depends on the nature of the data; usually it depends on the noise level as well as on the (average) width of the peaks.Finally, we would like to emphasize another interesting interpretation of the above smoothing approach: The convolution in (9) can be written as a scalar product of with the shifted Gaussian density (note that is symmetric), that is, . Thus, the entries of are actually the analysis coefficients of the Gaussian dictionary . The perspective of analyzing data by a dictionary offers several opportunities for generalization. For instance, one could also consider (redundant) dictionaries with more than one standard deviation or more sophisticated functions than .
Step 3: Standardizing the Data
The 1bit optimization of (7) does not incorporate a bias term. Hence, it is necessary to center the data first. For this, we compute the mean spectrum^{12}^{12}12Actually, we use the smoothed data vectors from Step 2 as input for this computation. But in order to keep the notation simple, we still write . This convention holds also for all forthcoming steps.
i.e., contains the average of the th entry of all spectra. The spectra are further scaled by dividing the nonconstant features by their standard deviation
The standardized spectra are then obtained by
In this way, all feature variables are centered and have an empirical standard deviation equal to , so that they get equally weighted in the selection process.
Step 4: Sparse Feature Selection
We are now ready to perform the actual feature extraction step, using the 1bit recovery method presented in Subsection 2.2:
Algorithm 2 (1Bit Compressed Sensing).
Input: Samples , sparsity parameter , threshold Output:Estimated feature vector Compute:The second part (in (11)) is a simple hard thresholding that tries to eliminate computational inaccuracies by setting almost zero entries of to ( is usually very small, e.g., ).
The actual feature selection takes place in (10). Recalling the observation model from (8), we conclude that the th sample is correctly classified by a vector if and only if . Hence, the objective functional of (10) will be particularly large if sufficiently many samples are correctly classified by . However, a consistent prediction of all measurements (i.e., for all ) is not strictly enforced, and therefore, our strategy enjoys a certain robustness against (random) perturbation of the model (8). This could occur in practice, for example, when a training sample was wrongly classified from the very beginning. On the other hand, the constraint of (10) guarantees that the maximizer will be “effectively” sparse (depending on the choice of the sparsity parameter ). This intuition indicates that the estimator will be indeed a sparse vector allowing for an appropriate separation of the two classes.
Step 5: Detecting the Connected Components
One advantage of Algorithm 2 is that it does not make any assumptions on the structure of the data vectors . Hence, it might be even suited for much more general types of data. However, its “universality” comes with the drawback that the characteristic peak structure of MS data is not captured at all. In fact, a spectrum does not consist of sharp spikes but rather widespread Gaussian shaped peaks. Hence, if the algorithm finds a significant feature position, say at the maximum of some peak, it usually tends to select also those features which are close to this position. Such a behavior is not very surprising, because nearby features are highly correlated to the maximum of the peak, and therefore, they may allow for a good separation as well.
Empirical results have shown that this process of selection “evolves” in a continuous fashion when changing the sparsity level . As a consequence, the support of a feature vector from Algorithm 2 typically consists of a few connected “intervals” (consecutive sequences of indices) which are centered around the selected peaks (see also Fig. 3). The actual sparsity of should be therefore measured by means of its connected intervals and not by simply counting its nonzero entries.
For this reason, we may easily improve the sparsity of by reducing every interval to its most significant entry:^{13}^{13}13Here denotes the support of , i.e., the set of indices corresponding to its nonzero entries.
Algorithm 3 (Sparsification of ).
Input: (Sparse) feature vector Output: Sparsified version Compute: Find the connected components of . For every do the following: Set all entries of in to , except from . The resulting vector is .Step 6: Dimension Reduction
This final (optional) step does not involve any further computations but shows how to proceed with our result . As mentioned before, the main purpose of SPA is not just to classify (unknown) samples, but rather to reduce the data to its significant entries (dimensions). Indeed, we may use for a dimension reduction: Let be some (possibly unknown) data vector. Then, we can project onto the selected feature positions of . More precisely, all entries that do not belong to are set to :
(12) 
The resulting data vector is now trivially embedded into a lowdimensional space of dimension .^{14}^{14}14In practice, one would simply reject all indices that are not contained in . But it still contains the most important information which has been found by the above algorithm. Note that we have not made any use of the actual values of but merely of its support.
By this projection, we may reduce the danger of overfitting. In particular, by working in a lowdimensional space, a large tool set from machine learning is now available for classification and clustering. But how to explicitly proceed with the data heavily depends on the specific application and is therefore not part of SPA.
4 Experimental Results: Feature Selection from Simulated DataSets
In this section, we assess our framework of SPA with regard to a typical situation in massspectrometry analysis: We would like to extract discriminating features from MS data with respect to two groups (e.g., healthy and diseased patients). A major difficulty is usually that only a small number of measurements (observations) is available. Building on this, we ask for the following: Given a simulated dataset for which the position and number of discriminating peaks are known (this will be called below), how many samples are needed to identify these features with high accuracy?
We shall compare our results to the widely used stateoftheart algorithms LIBLINEAR (regularized SVM) and the standard MATLAB implementation of Lasso.
4.1 Creating a Simulated DataSet
We assume that our sample set follows a certain joint random distribution , where each sample is independently drawn. In order to make the problem tractable, let us make two model assumptions on and . First, the mass spectra are generated as follows:
where determines the (random) amplitude of the th peak, specifies its position and shape, and represents the lowamplitude baseline noise. We shall assume that the amplitudes and the noise are Gaussian, that is, with positive definite and with . Note that the generated data might have negative components. This does not mimic the structure of realworld mass spectra which is always nonnegative. However, since centering is part of our preprocessing anyway (cf. Step 3 in Subsection 3.2), the assumption of meanzero amplitudes is quite natural. The (disease) labels are then simply modeled as 1bit observations (see also (8))
(13) 
where is the sparse groundtruth feature vector, which we intend to estimate. In the following, each nonzero entry of is located at the center of a specific peak (see Fig. 4(d)–(f)), so that actually determines all biologically relevant peaks (molecular structures). Since is invertible (i.e., the features are linearly independent), this collection of peaks is an optimal fingerprint in the sense that removing or adding any feature variable would decrease the prediction accuracy (with respect to the “perfect” model of (13)).
In our experiments, we create datasets , each one consisting of equidistant peaks (atoms ) shaped like Gaussian density function of width . The vector is chosen to have five nonzero components, which means that only five prechosen peaks were used to generate the labels . Hereafter, we will refer to these as condition positive peaks. Fig. 4 shows three different data instances magnifying only the first seven peaks, generated in the described way. In order to verify our method, we will use two types of datasets DS1 and DS2 which only differ in their correlation matrix . For DS1,
is chosen to be the identity matrix. This implies that the heights of all of the 200 peaks are standard Gaussian random variables. For DS2, we have chosen three pairs of negative peaks to be positively correlated and in addition, one condition positive peak was chosen to be positively correlated with one of the negative peaks. Thus, there are a few entries of value
off the main diagonal in . To test the algorithm’s performance increasing amount of Gaussian noise with was added to DS1 and DS2. These corresponds to signaltonoise (SNR) ratio of 10, 3.33 repectively ^{15}^{15}15Signaltonoise ratio was calculated as .. The values of SNR are chosen to represent the behaviour of the algorithm up to the levels of noise that are normally found in MS data.4.2 Setup and Evaluation Criteria
Let us recall the essential question of our experiments: Can we recover the support of , and if so, how many samples do we need for that? For this purpose, we shall successively increase the number of available samples in the (training) dataset and examine whether SPA (or Lasso, or SVM) succeeds in recovering . Since each of the considered algorithms involves a variable parameter, we have decided to perform an adaptive tuning for each problem instance. In fact, the sparsity parameter was chosen such that the resulting classifier matches the sparsity level of . But of course, this does not automatically imply that the supports of and completely coincide.^{16}^{16}16Due to the redundancy of the peakassociated feature variables (cf. Step 5 in Subsection 3.2), an estimated feature vector is considered to be equal to the groundtruth vector with some tolerance, which particularly depends on the width of the peaks. For each problem instance, the smallest sparsity parameter which resulted in a classifier with five nonzero entries was chosen in the following way: The initial value of the sparsity parameter for SPA and SVM (Lasso) was set to the value which corresponds to the classifier with less than (more than) five nonzero values^{17}^{17}17This difference arises from the implementation of Lasso.. For SPA and SVM (Lasso), the sparsity parameter was increased for a preset step size until the outcome had five or more (five or fewer) nonzero entries. If the previous step provided a sparse classifier with strictly more than (strictly less than) five nonzero entries, the bisection method was used on the interval between the two last sparsity parameter values. The bisection method was used until the optimal sparsity parameter was found or the difference between the two consecutive parameters became smaller than a preset tolerance.
We will use a measure based on sensitivity. Sensitivity, defined as^{18}^{18}18TP  true positives, i.e. correctly identified peaks
FP  false positives, i.e. incorrectly identified peaks
TN  true negatives, i.e. correctly rejected peaks
FN  false negatives, i.e. incorrectly rejected peaks
is an appropriate measure for our objectives because it represents an algorithm’s ability to detect the relevant features. Note that ideally, the number of condition positives () is equal to predicted condition positives (). In such a situation, the precision, given by is equal to the sensitivity. However, in the presence of noise it is possible that the final selection encompasses several features which are associated with a single peak. This could lead to a precision value equal to if all of the selected values are declared as true positives, though some other true features remain undetected. Since for us, it is equally important to penalize both false positives and false negatives, we have chosen the sensitivity to be the main point of reference. A measure of similar importance is the specificity, which is defined by
Finally, due to the possibly imbalanced number of relevant features, we shall also take into account the socalled balanced accuracy
4.3 Results
Datasets of sample sizes between and were generated as described above and each of the methods was performed for standardized input data. Note that the hard thresholding step described in (11) was also applied to the classifiers obtained from Lasso or SVM. Otherwise, any computational inaccuracy would completely destroy the sparsity structure of the results.
For the sake of statistical stability, each experiment was repeated times. The averaged results are presented in the Fig. 5. We can see that SPA () performs better than the SVM or Lasso with regard to the capability of recognizing the true positive features (sensitivity in Fig. 5). In our setting, if one method fails to select 5 condition positive peaks because one of them was selected twice, and the other method selects exactly the same 4 peaks and one false positive in addition, the specificity penalizes only the latter one. But effectively, both cases are suboptimal, since only the 5 positive peaks together can predict the class correctly. This effect is reflected by a smaller value of specificity of the 1bit approach comparing to the specificity of other two methods for datasets with less than 300 spectra (column 2 in Fig. 5). However, this also implies that SPA performs sligtly worse in rejecting true negatives than the other two approaches. The average results for balanced accuracy are visualized in the third column of Fig. 5. We observe that SPA outperforms the other two methods and even achieves 100% accuracy with relatively few observations. With further decreasing SNR the performance of the three algorithms becomes more similar. Fig. 6 shows the numerical outcomes for the dataset DS2. The nontrivial correlation structure of DS2 eventually leads to a slight drop of sensitivity and accuracy for SPA (compared to DS1), whereas the performance of the other two methods essentially remains unaffected. As before with further decreasing SNR the performance of the three algorithms becomes more similar in terms of sensitivity and balanced accuracy.
5 Experimental Results: Analyzing RealWorld MALDITOF MS Data
In this section, we present results of SPA, Lasso, and SVM for analyzing realworld massspectrometry data and compare them to the MALDIquant proteomics analysis workflow [37]. All data was acquired in our earlier studies [38, 11]. It was approved by the local ethics committees and fulfils the requirements of the Helsinki declaration. All subjects gave informed consent to participate in the study. We will demonstrate the performance of our method on two datasets:

Spiked Data: The spiked dataset is a labelled groundtruth dataset containing control (e.g. healthy) and case (e.g. diseased) mass spectra where the true labels are known. It is created from human blood samples^{19}^{19}19Blood serum of apparently healthy individuals from a clinical study ([38]) was used. which were either unchanged (control group) or in which a proteinmix has been mixed (spiked) into (case group). In order to simulate different strength of an effect caused e.g. by a disease, we further subdivided the case group into five subgroups where the amount of spikedin proteins is increasing. The five volumes in the case subgroups were spiked with the following concentrations of the protein mix^{20}^{20}20Protein calibration standard mix Part No.: 206355 & 206196) from Bruker Daltronics (Leipzig, Germany): 0.075pMol/L, 3.03pMol/L, 0.30nMol/L, 0.76nMol/L and 121.21nMol/L. This mix contains the hormones Angiotensin, ACTH, clip 1839, Substance P and the cell protein Ubiquitin. The peptide mix was added before sample pretreatment and 64 spectra were measured due to 4fold spotting (technical replicates). Mass spectra were acquired using the protocol described in the supplementary material (S1). Each volume corresponds to a dataset. What differentiates the datasets are the amplitudes of the 6 spikes resulting from the added substances. The signaltonoise ratio of the spikedin peaks is shown in the Fig. 7^{21}^{21}21The power of noise for each of the 5 analyzed datasets is estimated as an average of intensity of noise of the observations using median absolute deviation..

Pancreas Cancer Data (P. CA): A total of 120 patients with pancreatic cancer and controls were recruited for this study [11]
. For the discovery study sera were obtained from two different clinical centres (University Hospital Leipzig (UHL, set L) and Heidelberg (UHH, set H)) as described in the supplementary material (S1). Note that each acquired spectrum has been assigned a classlabel, i.e., healthy or diseased. So, the health status of the training samples is known in advance (supervised learning).
Baseline removal was performed on the raw MS data using TopHat filtering ([39]). In particular, no additional calibration or noise reduction steps have been applied. More information on the data and sample preparation can be found in the supplementary material (S1).
The missingdata problem
When dealing with data coming from measurements of, say, a Mass Spectrometer instrument, the so called missingdata problem usually occurs. This means that the instrument failed to give measurements for some of the measured masses, usually due to the stochastic nature of the process happening inside the device. Due to the smoothing step in our algorithm and the arguments of e.g. Rubin et al. ([40]) this problem can be mainly ignored in our case for identifying the relevant features. However, this does not necessarily hold for the classification step, i.e. applying the identified sparse classifier to an unknown dataset. In this scenario, where data is missing in an unknown sample, there are basically two options: (1) applying a method for inferring the missing data or (2) stopping the classification and return an error message to the user. In this work we decided to follow the latter approach, since inferring missing data is not in the scope of this paper^{22}^{22}22The interested reader might find a good starting point about this topic in these two reviews [41, 42] but is an unarguable crucial point in any data analysis pipeline and should depend on the actual usecase.
Accuracy vs. Number of Features
We performed these experiments with respect to the same evaluation categories as in the case of simulated data. Note that the normalization and standardization as described in Section 3.2 were applied as preprocessing steps in each of the methods. Similarly, a hard thresholding as described in (11) was applied to all classifiers estimated by the examined algorithms.
For the each of the algorithms, we are testing the performance of the obtained classifiers learned on the pure dataset which corresponds to the condition negative class and one spiked dataset at a time corresponding to the condition positive class.
The results of the classifier with nonzeros on the spiked dataset are shown in Table 2. The main question in these experiments is how successful each of the algorithms is in detecting the peaks that were initially spiked (see the dataset description above). We can see that the values of sensitivity for SPA are at least as high as those of the other methods, which implies that the approach of bit CS is very competitive in this situation and mostly achieves the best detection rate. However, the relatively poor performance of all the algorithms on the spiked dataset can be explained by the nature of the data. Since the peptide mix was added to the blood samples before acquiring the mass spectra, the spiked peaks are not always present in all the resulting mass spectra in the positions where we expect to find them. There exist datasets for which all the mass spectra failed to exhibit certain spiked peaks at their expected locations as can be seen in the Fig. 7. Thus, we cannot expect any of the algorithms to find these missing peaks. Nonetheless, there is still a chance to build a reliable fingerprint out of the remaining spikes while there is no chance to detect the missing spikes because the dataset is not rich enough to represent it. On the other hand, this spiked dataset combines the advantages of both simulated and clinical data, since the positions of the desired biomarkers are known in advance while their representative behavior in the spectra is quite realistic.
In contrast to that in the case of pancreas cancer datasets, we do not know the truepositive feature positions. Consequently, we can only rely on the classification performance of the obtained sparse classifiers by each of the algorithms. To evaluate the reliability of our results, for each of the methods, we have employed the crossvalidation scheme as described in the Algorithm 4 with the number of folds set to 5. In order to ensure statistical stability, each experiment was repeated 10times. Fig.8 shows the average results over 10 repetitions.
Algorithm 4 (CrossValidation of Classification Performance).
Input: Raw data Number of CVfolds ; Output: Classification accuracy ; Average sparsity (number of selected features) Split the sample set randomly into disjoint folds of (almost) equal size. For each fold perform the following steps 24: Compute the feature vector employing the desired method on the samples of Dimension reduction as described in Sec. 3.2 (step 6). Project all spectra onto and put Classification of : Use the projected samples of to predict the labels of the spectra by an ordinary SVM^{23}^{23}23Here, the standard MATLAB implementation of SVM was used.. Denote the prediction accuracy by . Compute the average accuracy and average sparsityIn order to ensure statistical stability, each experiment was repeated times. Fig. 8 shows the results. Note that our results show that accurate predictions are already possible with a very few features, so that the assumption of small disease fingerprint seems to hold for this dataset. Furthermore, it can be seen that SPA is especially well suited for situations where a sparse classifier (containing only very few features) is preferred. This is appealing because fewer features enable an easier interpretation of the actual components of a potential disease fingerprint. Moreover, followup experiments that often involve an individual treatment of each component (e.g., potential biomarkers) would become much less costly. Note that in the nonsparse region with more than features selected, it is not meaningful to relate the achieved accuracy to the quality of the learned feature vector due to the small sample size. The considered algorithms assume the underlying fingerprint to be sparse. This assumption usually does not fully hold in practice. Therefore, we cannot expect that a learned feature vector achieves perfect classification. The classification accuracy should be therefore considered as an indicator of how well our model assumption of the sparse fingerprint fits to the unknown groundtruth. If we let the algorithms operate out of the region for which they have been designed for, we may achieve indeed a higher accuracy, but this is probably a consequence of overfitting. And even more importantly, the learned feature vector (model) is not reliable anymore.
Best Classifier
Apart from that, we are interested in the performance of the best sparse classifier (i.e. small number of features) found by each of the algorithms (SPA, Lasso, SVM). For all learned classifiers with to nonzero components, Table 3 presents those with the best classification accuracy. Furthermore, we also considered a typical analysis pipeline (MALDIQuant) to see how the “purelydatabased” approaches (SPA, Lasso, SVM) compare to a modelbased approach^{24}^{24}24By “modelbased” we mean that specific model assumptions on the data are made and exploited, such as noisestructure for denoising or Gaussianshaped structures for peak detection.. In Table 3, it can be seen that SPA provides the sparsest solutions while achieving competitive results with respect to sensitivity and specificity at the same time. Lasso and SVM select almost the same features and therefore perform similarly. On the other hand, MALDIQuant selects the features based on a prior modelbased peak detection followed by a feature selection based on shrinkage diagonal discriminant analysis ([43]). But however, it still performs worst on the UHL dataset.
Medical Interpretation of Results
Pancreatic cancer is not only a common and increasingly frequent [44], but also still a fatal disease, with a survival rate of 35% five years after diagnosis [45]. The conventional tumor marker, Carbohydrate Antigen 199 (CA199), as a blood group antigen not present in a significant proportion of the patients [46], shows insufficient diagnostic sensitivity and specificity (AUC 0.71), even in combination with the secondline tumor marker Carcinoembryonic Antigen (CEA, combined AUC 0.75) [47]. The need for better markers for screening and differential diagnosis is evident, as panceratic carcinoma would be principally curable if detected and identified very early in the course of the disease. Along with the emerging “omics”technologies great hope was rised to find tumorspecific peptides or metabolic alterations to increase sensitivity and specificity of early and differential diagnostics, and several combinatory marker models could be identified by proteomics [11] and metabolomics [48]. Pancreatic carcinoma is a complex disease  it affects the metabolism as a whole (e.g. the socalled Warburg effect) [49], but also alters proteolytic activity [50]. Therefore, it might be naïve to expect a single marker capable to indicate presence, progression and exact type of the malignancy at once [51] 9– it might even be overly reductionistic to attribute these capabilities to a single model, even if it consists of several entities measured by different “ omics” technologies [46]. As Raftery states “basing inferences on a single “best” model as if the single selected model were true ignores model uncertainty, which can result in underestimating uncertainty about quantities of interest” [52], and the larger the “omics” datasets grow, the larger is the ‘probability, that there is not one “single best” predictive marker model, but instead several with comparable selectivity [51]. And it is very reasonable to assume that, even on the same dataset, different algorithms might favor different models consisting of different feature sets and bring forth completely different results, when only the best differentiating models are regarded. For an indepth comparison of the validity of the results of different algorithms, the underlying peak features should also be taken into account, and similarities in the selected features corroborate the algorithms superimposed on them. In the case of our study, we have the great advantage, that the same dataset was evaluated in three different studies: the principal one by Fiedler et al. [11], a subsequent BinDAalgorithmbased manuscript by Gibb and Strimmer published recently [53], and the present one. Fiedler et al. [11] identified one discriminating peptide, Platelet Factor 4 (m/z 3884, identified in italics, double hits in bold) within four discriminating peaks (m/z 3194, 3884, 4055, and 5959). The 30 most differential peaks in Gibb et al. [53] were m/z 4495, 8868, 8989, 1855, 4468, 8937, 2023, 1866, 5864, 5946, 1780, 2093, 5906, 5960, 8131, 1207, 4236, 2953, 9181, 1021, 1466, 4092, 4251, 5005, 8184, 1897, 3264, 2756, 6051, and 1264, with m/z 8937 identified as pancreatic progenitor cell differentiation and proliferation factorlike protein. m/z 3884 could not be identified as discriminating marker (while it might play a role in pancreatic carcinoma nonetheless [54]), whereas m/z 1466 can be attributed to a fragment of fibrinopeptide A (DSGEGDFLAEGGGVR), as previously described in tumor samples [55]. In the present study, the peaks m/z 1464, 1546, 1944, 5904, 1619, 4209, and 2662 could be identified as discriminating features. The slight mass shift of about 2 Da for m/z 1464 / 1466 and 5904 / 5906 is probably arising from different peak preprocessing procedures, peaks are wide enough to tolerate this deviation. Further investigations and the application of further algorithms on the same dataset are highly likely to yield a similar, partially overlapping set of features, each with a comparable discriminating power (Fiedler et al. [56] AUC 1.0; Gibb et al. [53] in a 5feature model: accuracy of 0.96, sensitivity of 0.96, specificity of 0.97, positive predictive value of 0.97 and negative predictive value of 0.95; the present study accuracy 0.96, sensitivity 0.97, specificity 0.95 and accuracy 0.98, sensitivity 0.99, specificity 0.97. This also corresponds to a recently published comparable study investigating a glycoprotein marker panel (AUC 0.95) [57]. Biomarkers for clinical diagnostics comprise a wide field of applications (e.g. populationwide screening, early diagnostics, characterization, treatment guidance, efficacy and toxicity monitoring, prognosis, susceptibility estimation and many more) [46], each with special requirements for sensitivity and specificity, that are only partially condensed in the AUC as an overall selectivity measure [51]. Especially for screening purposes, sensitivity is extremely important [48], and clinically applied tests e.g. for newborn screening frequently surpass the 0.99 hallmark [56]. Compared with the conventional, “notforscreening” marker CA199, the SPAbased model shows considerable improvement, however there is still a big gap to screening suitability, which in the next years might be bridged by improved sensitivity of new instrumentation, refined algorithms (as the SPA), and combination with other “markers” from the “big data” field, enabling a more holistic view – not only of the disease, but also of the affected patient[46].
6 Conclusion
In this work, we introduced Sparse Proteomics Analysis (SPA), a new framework for the analysis of proteomics data, generated by mass spectrometry measurements. The framework solves the problem of selecting a minimum set of features from highdimensional data in a situation where only relatively few measurements are available. Our approach particularly allows for a biomedical interpretation and enables a classification of unknown samples. This is done by formulating and solving a regularized optimization problem, using ideas from bit compressed sensing combined with several generic pre and postprocessing steps. We have shown by several numerical experiments that SPA can indeed compete with standard (and widely used) algorithms as well as with a specifically adapted method for the analysis of proteomics data (MALDIQuant). In the “verysparse” situation, it has even turned out that SPA outperforms the other approaches with respect to prediction accuracy.
Competing Interests
The authors declare that they have no competing interests.
Author’s Contributions
Conceived and designed the experiments: TC, MG, JV, GK, and CS. Performed the experiments: NC, TC, and NW. All authors contributed to writing the paper. All authors read and approved the final manuscript.
Acknowledgments
JV was supported by the ERC CZ grant LL1203 of the Czech Ministry of Education. TC, MG, NC, JV, GK and CS were supported by the Einstein Center for Mathematics Berlin (ECMath), project grant CH2, and by the DFG Research Center Matheon Mathematics for key technologies, Berlin. TC and CS are supported by the German Ministry of Research and Education (BMBF) project Grant 3FO18501 (Forschungscampus MODAL). GK acknowledges support by the Einstein Foundation Berlin, by the Deutsche Forschungsgemeinschaft (DFG), and by the DFG Collaborative Research Center TRR 109 Discretization in Geometry and Dynamics.
The authors are thankful to Irena Bojarovska for fruitful discussions and help conducting the experiments.
Supporting Information
S1  Mass Spectrometry Data Generation
Chemicals, Standards and Consumables
Gradient grade acetonitrile, ethanol, and HPLCwater were obtained from J.T. Baker (Phillipsburg, NJ, USA); p.a. trifluoroacetic acid (TFA) and acetone were purchased from SigmaAldrich (Taufkirchen, Germany). The peptide and protein MALDITOF calibration standards I and cyano4hydroxycinnamic acid (HCCA) were purchased from Bruker Daltonics (Bremen, Germany). Automated magnetic bead preparations were performed using 96 well plates, TubePlates from Biozym (Hessisch Oldendorf, Germany), polypropylene tubes (low profile) from Abgene (Surrey, UK), and modular reservoir quarter modules from Beckman (Fullerton, USA). For sample storage 450 L CryoTubesTM were purchased from Sarstedt (Nmbrecht, Germany). Multifly needle sets and polypropylene serum monovettes with clotting activators were also obtained from Sarstedt.
Peptidome Separation
All serum samples of the discovery set were processed at one time and analyzed simultaneously to avoid proceduredependent errors. The external validation set was prepared, processed and analyzed separately. Peptidome separation of the samples was performed using the ClinPro Tools profiling purification kits from Bruker Daltonics. Magnetic particles with defined surface functionalities (magnetic beadimmobilized metal ion affinity chromatography (MBIMAC Cu), magnetic beadhydrophobic interaction (MBHIC C8) and weak cation exchange (MBWCX)) were processed by the ClinPro Tools liquid handling robot according to the manufacturers protocol (Bruker Daltonics). Serum specimens were thawed on ice for 30 min and immediately processed according to our standardized protocol for serum peptidomics [58].
Mass Spectrometry
A linear MALDITOF mass spectrometer (Autoflex I, Bruker Daltonics) was used for the peptidome profiling. Daily mass calibration was performed using the standard calibration mixture of peptides and proteins in a mass range of 110 kDa. Mass spectra were recorded and processed using AutoXecute tool of the flexControl acquisition software (Ver. 2.0; Bruker Daltonics). The settings were applied as follows: Ion source 1: 20 kV; ion source 2, 18.50 kV; lens, 9.00 kV; pulsed ion extraction, 120 ns; nitrogenpressure, 2500 mbar. Ionization was achieved by a nitrogen laser (=337 nm) operating at 50 Hz. For matrix suppression a high gating factor with signal suppression up to 500 Da was used. Mass spectra were detected in linear positive mode.
Baseline Removal
The baseline is an exponential like offset dependent on the value (masstocharge; xvalue). A baseline correction is performed to remove this rather lowfrequency noise from the spectrum. We use a morphological TopHat filter to eliminate certain spatial structures within the signal, in our case the baseline. Note that this technique does not produce negative intensity values.
References
 1. Rissin D, Kan C, Campbell T, Howes S, Fournier D, Song L, et al. Singlemolecule enzymelinked immunosorbent assay detects serum proteins at subfemtomolar concentrations. Nat Biotechnol. 2010;28.
 2. Aebersold R, Mann M. Mass spectrometrybased proteomics. Nature. 2003;422(6928):198–207.
 3. Petricoin EF, Belluco C, Araujo RP, Liotta LA. The blood peptidome: a higher dimension of information content for cancer biomarker discovery. Nat Rev Cancer. 2006 Dec;6(12):961–967. Available from: http://dx.doi.org/10.1038/nrc2011.
 4. Rai AJ, Chan DW. Cancer proteomics: Serum diagnostics for tumor marker discovery. Ann N Y Acad Sci. 2004 Jun;1022:286–294. Available from: http://dx.doi.org/10.1196/annals.1318.044.
 5. Coombes KR, Morris JS, Hu J, Edmonson SR, Baggerly KA. Serum proteomics profiling–a young technology begins to mature. Nat Biotechnol. 2005 Mar;23(3):291–292. Available from: http://dx.doi.org/10.1038/nbt0305291.
 6. Liotta LA. Clinical proteomics: written in blood. Nature. 2003;425.
 7. Phizicky E, Bastiaens PIH, Zhu H, Snyder M, , Fields S. Protein analysis on a proteomic scale. Nature. 2003;422.
 8. Issaq HJ, Xiao Z, Veenstra TD. Serum and plasma proteomics. Chem Rev. 2007 Aug;107(8):3601–3620. Available from: http://dx.doi.org/10.1021/cr068287r.
 9. Stühler K, Meyer HE. MALDI: more than peptide mass fingerprints. Curr Opin Mol Ther. 2004 Jun;6(3):239–248.
 10. Sitek B, WalderaLupa DM, Poschmann G, Meyer HE, Stühler K. Application of labelfree proteomics for differential analysis of lung carcinoma cell line A549. Methods Mol Biol. 2012;893:241–248. Available from: http://dx.doi.org/10.1007/9781617798856_16.
 11. Fiedler GM, Leichtle A, Kase J, Baumann S, Ceglarek U, Felix K, et al. Serum Peptidome Profiling Revealed Platelet Factor 4 as a Potential Discriminating Peptide Associated With Pancreatic Cancer. Clin Cancer Res. 2009 June;15(11):3812–3819. Available from: http://publications.mi.fuberlin.de/155/.
 12. Strenziok R, Hinz S, Wolf C, Conrad TOF, Krause H, Miller K, et al. Surfaceenhanced laser desorption/ionization timeofflight mass spectrometry: serum protein profiling in seminoma patients. World J of Urology. 2009 June;28(2):193–197. Available from: http://publications.mi.fuberlin.de/156/.
 13. Leichtle A, Nuoffer JM, Ceglarek U, Kase J, Conrad TOF, Witzigmann H, et al. Serum amino acid profiles and their alterations in colorectal cancer. Metabolomics. 2011 September;Available from: http://publications.mi.fuberlin.de/1081/.
 14. Diao L, Clarke CH, Coombes KR, Hamilton SR, Roth J, Mao L, et al. Reproducibility of SELDI Spectra Across Time and Laboratories. Cancer Inform. 2011;10:45–64. Available from: http://dx.doi.org/10.4137/CIN.S6438.
 15. Marrugal A, Ojeda L, PazAres L, MolinaPinelo S, Ferrer I. ProteomicBased Approaches for the Study of Cytokines in Lung Cancer. Disease Markers. 2016;.
 16. Tang S, Zhou F, Sun Y, Wei L, Zhu S, Yang R, et al. CEA in breast ductal secretions as a promising biomarker for the diagnosis of breast cancer: a systematic review and metaanalysis. Breast Cancer. 2016;p. 1–7. Available from: http://dx.doi.org/10.1007/s1228201606809.
 17. Le N, Sund M, Vinci A, Beyer G, Javed MA, Krug S, et al. Prognostic and predictive markers in pancreatic adenocarcinoma. Digestive and Liver Disease. 2016;48(3):223 – 230. Available from: http://www.sciencedirect.com/science/article/pii/S1590865815006908.
 18. Donoho DL. Compressed sensing. IEEE Trans Inform Theory. 2006;52:1289–1306.

19.
Candés EJ, Tao T.
Decoding by linear programming.
IEEE Trans Inform Theory. 2005;51:4203–4215.  20. Candés EJ, Romberg J, Tao T. Stable signal recovery from incomplete and inaccurate measurements. Comm Pure Appl Math. 2006;59:1207–1223.
 21. Genkin A, Lewis D, Madigan D. Largescale Bayesian logistic regression for text categorization. Technometrics. 2007;49:291–304.
 22. Friedman J, Hastie T, Tibshirani R. Regularized paths for generalized linear models via coordinate descent. Department of Statistics, Stanford University; 2008.
 23. Efron B, Hastie T, Johnstone I, Tibshirani R. Least angle regression. Ann Statist. 2004;32:407–499.
 24. Koh K, Kim S, Boyd S. An InteriorPoint Method for LargeScale l1Regularized Least Squares. Selected Topics in Signal Processing. 2007;1(4):606–617. Available from: http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=4407767.
 25. Wu TT, Lange K. Coordinate descent algorithms for lasso penalized regression. Ann Appl Stat. 2008;2:224–244.
 26. Vapnik VN. Statistical Learning Theory. John Wiley & Sons; 1998.
 27. Genzel M, Kutyniok G. Towards a Mathematical Theory of Feature Selection from RealWorld Data with NonLinear Observations. preprint. 2016;.
 28. Chen SS, Donoho DL, Saunders MA. Atomic decomposition by basis pursuit. SIAM J Sci Comput. 1998;20:33–61.
 29. Tibshirani R. Regression shrinkage and selection via the Lasso. J Royal Statist Soc B. 1996;58:267–288.
 30. Boufounos PT, Baraniuk RG. 1Bit compressive sensing. In: Proceedings of the 42nd Annual Conference on Information Sciences and Systems; 2008. p. 16–21.
 31. Plan Y, Vershynin R. Onebit compressed sensing by linear programming. Comm Pure Appl Math. 2013;66:1275–1297.
 32. Plan Y, Vershynin R. Robust 1bit compressed sensing and sparse logistic regression: a convex programming approach. IEEE Trans Inf Theory. 2013;59(1):482–494.
 33. Zou H, Hastie T. Regularization and variable selection via the elastic net. J Royal Statist Soc B. 2005;67(2):301–320.
 34. Davenport MA, Duarte MF, Eldar YC, Kutyniok G. Introduction to compressed sensing. Cambridge Univ. Press; 2012.
 35. Foucart S, Rauhut H. A mathematical introduction to compressive sensing. Birkhäuser; 2013.
 36. Bühlmann P, Van De Geer S. Statistics for highdimensional data: methods, theory and applications. Springer; 2011.
 37. Gibb S, Strimmer K. MALDIquant: a versatile R package for the analysis of mass spectrometry data. Bioinformatics. 2012;28(17):2270–2271.
 38. Kratzsch J, Fiedler GM, Leichtle A, Brügel M, Buchbinder S, Otto L, et al. New reference intervals for thyrotropin and thyroid hormones based on National Academy of Clinical Biochemistry criteria and regular ultrasonography of the thyroid. Clin Chem. 2005 Aug;51(8):1480–1486. Available from: http://dx.doi.org/10.1373/clinchem.2004.047399.
 39. Sauve AC, Speed TP. Normalization, baseline correction and alignment of highthroughput mass spectrometry data. In: Proceedings of the Data Procedings Gensips; 2004. .
 40. Rubin DB. Inference and missing data. Biometrika. 1976;63(3):581–592.
 41. Pigott TD. A Review of Methods for Missing Data. Educational Research and Evaluation. 2001;7(4):353–383.

42.
Schafer JL, Olsen MK.
Multiple Imputation for Multivariate MissingData Problems: A Data Analyst’s Perspective.
Multivariate Behavioral Research. 1998;33(4):545–571. PMID: 26753828.  43. Ahdesmäki A, Strimmer K. Feature selection in omics prediction problems using cat scores and false nondiscovery rate control. Ann Appl Stat. 2010;4.
 44. Yeo TP, Lowenfels AB. Demographics and epidemiology of pancreatic cancer. Cancer J. 2012;18(6):477–484. Available from: http://dx.doi.org/10.1097/PPO.0b013e3182756803.
 45. Michl P, Pauls S, Gress TM. Evidencebased diagnosis and staging of pancreatic cancer. Best Pract Res Clin Gastroenterol. 2006 Apr;20(2):227–251. Available from: http://dx.doi.org/10.1016/j.bpg.2005.10.005.
 46. Leichtle A. Biomarker – vom Sein und Wesen. J Lab Med. 2015;39.
 47. Reitz D, Gerger A, Seidel J, Kornprat P, Samonigg H, Stotz M, et al. Combination of tumour markers CEA and CA199 improves the prognostic prediction in patients with pancreatic cancer. J Clin Pathol. 2015 Mar;Available from: http://dx.doi.org/10.1136/jclinpath2014202451.
 48. Leichtle A, Ceglarek U, Weinert P, Nakas CT, Nuoffer JM, Kase J, et al. Pancreatic carcinoma, pancreatitis, and healthy controls  metabolite models in a threeclass diagnostic dilemma. Metabolomics. 2013 May;9(3):677–687. Available from: http://publications.mi.fuberlin.de/1165/.
 49. Zhou W, Capello M, Fredolini C, Racanicchi L, Piemonti L, Liotta LA, et al. Proteomic analysis reveals Warburg effect and anomalous metabolism of glutamine in pancreatic cancer cells. J Proteome Res. 2012 Feb;11(2):554–563. Available from: http://dx.doi.org/10.1021/pr2009274.
 50. Brand RE, Nolen BM, Zeh HJ, Allen PJ, Eloubeidi MA, Goldberg M, et al. Serum biomarker panels for the detection of pancreatic cancer. Clin Cancer Res. 2011 Feb;17(4):805–816. Available from: http://dx.doi.org/10.1158/10780432.CCR100248.
 51. Leichtle AB, Dufour JF, Fiedler GM. Potentials and pitfalls of clinical peptidomics and metabolomics. Swiss Med Wkly. 2013;143:w13801. Available from: http://dx.doi.org/10.4414/smw.2013.13801.
 52. Raftery AE, Madigan D, Hoeting JA. Bayesian model averaging for linear regression models. JASA. 1997;92(437):179–191.
 53. Gibb S, Strimmer K. Differential Protein Expression and Peak Selection in Mass Spectrometry Data by Binary Discriminant Analysis. Bioinformatics. 2015;31(19).
 54. Poruk KE, Firpo MA, Huerter LM, Scaife CL, Emerson LL, Boucher KM, et al. Serum platelet factor 4 is an independent predictor of survival and venous thromboembolism in patients with pancreatic adenocarcinoma. Cancer Epidemiol Biomarkers Prev. 2010 Oct;19(10):2605–2610. Available from: http://dx.doi.org/10.1158/10559965.EPI100178.
 55. Villanueva J, Shaffer DR, Philip J, Chaparro CA, ErdjumentBromage H, Olshen AB, et al. Differential exoprotease activities confer tumorspecific serum peptidome patterns. J Clin Invest. 2006 Jan;116(1):271–284.
 56. Ceglarek U, Leichtle A, Brügel M, Kortz L, Brauer R, Bresler K, et al. Challenges and developments in tandem mass spectrometry based clinical metabolomics. Mol Cell Endocrinol. 2009 Mar;301(12):266–271. Available from: http://dx.doi.org/10.1016/j.mce.2008.10.013.
 57. Nie S, Lo A, Wu J, Zhu J, Tan Z, Simeone DM, et al. Glycoprotein biomarker panel for pancreatic cancer discovered by quantitative proteomics analysis. J Proteome Res. 2014 Apr;13(4):1873–1884. Available from: http://dx.doi.org/10.1021/pr400967x.
 58. Baumann S, Ceglarek U, Fiedler GM, Lembcke J, Leichtle A, Thiery J. Standardized approach to proteome profiling of human serum based on magnetic bead separation and matrixassisted laser desorption/ionization timeofflight mass spectrometry. Clin Chem. 2005 Jun;51(6):973–80.