Analytical model of peptide mass cluster centres with applications

Background The elemental composition of peptides results in formation of distinct, equidistantly spaced clusters across the mass range. The property of peptide mass clustering is used to calibrate peptide mass lists, to identify and remove non-peptide peaks and for data reduction. Results We developed an analytical model of the peptide mass cluster centres. Inputs to the model included, the amino acid frequencies in the sequence database, the average length of the proteins in the database, the cleavage specificity of the proteolytic enzyme used and the cleavage probability. We examined the accuracy of our model by comparing it with the model based on an in silico sequence database digest. To identify the crucial parameters we analysed how the cluster centre location depends on the inputs. The distance to the nearest cluster was used to calibrate mass spectrometric peptide peak-lists and to identify non-peptide peaks. Conclusion The model introduced here enables us to predict the location of the peptide mass cluster centres. It explains how the location of the cluster centres depends on the input parameters. Fast and efficient calibration and filtering of non-peptide peaks is achieved by a distance measure suggested by Wool and Smilansky.


Background
The mass spectrometric (MS) technique is widely used to identify proteins in biological samples [1][2][3][4]. The proteins are cleaved into peptides by a residue specific protease, e.g. trypsin. The resulting cleavage products can then be analysed by Peptide Mass Fingerprinting (PMF) [5] or subjected to MS/MS fragment ion analysis [6,7], which both rely on the comparison of peptide or peptide fragment ion spectra with spectra simulated from protein sequence databases [8].
The sensitivity and specificity of the peptide identification can be increased by various post-processing methods, for example calibration [9][10][11][12] and identification of non-peptide peaks [10,13,14]. The fact that peptide masses are not uniformly distributed across the mass range but form equidistantly spaced clusters [15] is employed by some of these methods. In dependence on the atomic composition of the peptide, the monoisotopic mass would emerge below (e.g. cystein rich peptides) or above (e.g. lysine rich peptides) the cluster centres. The deviation from the cluster centre is a result of the mass defect, which is the differ-ence between the nominal mass and the monoisotopic mass (Table 1). The mass defect is a result of atom fusion [16,17].

Calibration
Mass spectrometric peptide peak-lists of peptide mass finger print experiments [18] can be calibrated by comparing the location of measured peptide masses with the location of the peptide mass cluster centres. Gras et al. [19] suggested the use of maximum likelihood methods in order to determine the calibration coefficients a and b. They defined the likelihood function by: where m i is the i-th mass in the peak-list, and Δm is a search window. P(m, Δm) is the probability to find a mass in [m, m + Δm] given the theoretical distribution of peptide masses. The parameters a, b for arg max ∑ i P(am i + b, Δm) can then be used to calibrate the peak-lists. The authors, however, do not provide information on whether P(m, Δm) was determined from the exact distribution of the peptide masses or if a model approximating the distribution was used. They also do not mention which algorithm was used to maximise the likelihood. They reported that a mass measurement accuracy of 0.2Da and better was obtained after calibration.
Wool and Smilansky [10] have used Discrete Fourier Transformation (DFT) to determine the frequency λ and phase ϕ of a peak-list or mass spectrum. By comparing the experimental λ and ϕ with the theoretical λ = 1.000495 and ϕ = 0, they determined the slope and intercept of the calibration function. The authors reported a 40 -60% reduction of the mass measurement error. Furthermore, they presented a scoring scheme for sequence database searches. This scoring scheme approximates the probability P(m, Δm) to observe a peptide peak of mass m with given measurement error Δm.

Matrix noise filtration
The most widely used MALDI matrices for the analysis of peptides are 3,5-Dimethoxy-4-hydroxycinnamic acid (synapic acid), alpha-Cyano-4-hydroxycinnamic acid (alpha cyano) [20] and 2,5-dihydroxybenzoic acid (DHB) [21]. Unfortunately, clusters of matrix molecules can be ionised and cause peaks in the same mass range where peptide peaks are measured. Matrix aggregate formation can be minimised but not eliminated by adding ammonium acetate [21].
Some of the database search scoring schemes incorporate the number of signals (peaks) not assigned to a protein when computing the identification scores [22]. Therefore, the presence of matrix signals in MS spectra decreases the sensitivity of the MS spectra interpretation. Hence, the removal of peaks strongly deviating from the cluster centres is applied [21,23]. The measure of deviation from cluster centres introduced here provides a simple tool to filter non-peptide peaks.

Data reduction
A further application which employs the property of peptide mass clustering is the binning of the mass measurement range. By applying this technique the amount of data is reduced, thus increasing the speed with which the pairwise comparison of spectra can be made [24,25].
All these applications require us to know the exact location of or the distance between the peptide mass cluster centres. The distance between the cluster centres, which we will henceforth call wavelength λ, is commonly computed by first generating an in silico digest of the database. Afterwards, the linear dependence between the decimal point and the integer part is determined by regression analysis, for a relatively small mass range of 500 to 1000Da [23]. Various authors report different values of the distance between clusters: Wool and Smilansky reported 1.000495 [10], Gay et al. 1.000455 [15], while Tabb et al. used a wavelength of 1.00057 [24].
In this work we present an analytical model allowing us to predict the mass of the peptide cluster centres. The parameters of the model include: the frequencies of the amino acids in the sequence database [26], the average protein length of the proteins in the database, the cleavage sites of the proteolytic enzyme and the cleavage probability. Based on this model we introduced a measure of deviation of peptide masses from the nearest cluster centre, which is a refinement of a measure proposed by Wool and Smilansky [10]. Using this distance measure, we developed a calibration procedure which employs least squares linear regression in order to determine the affine model of the mass measurement error and subsequently to calibrate the spectra. Using this method we reached higher calibration accuracy as reported by Wool and Smilansky [10], and Gras et al [19]. We used the same distance measure to identify and remove non-peptide peaks prior to database searches performed by the Mascot search engine [22].

Results and discussion
A simple way to predict the peptide mass cluster centres of a protein database Figure 1 shows the mass defect, the difference of the monoisotopic (m (M) ) and nominal (m (N) )masses of peptides of a sequence specific in silico protein sequence data-base digest [27], as a function of m (N) . The peptides were produced with the restriction that no missed cleavages were allowed. A strong linear dependence of the mass defect on m (N) can be observed.
The peptide mass rule  The first model of this dependence which we examined was m (M) -m (N) = c 1 ·m (N) . We fixed the intercept at 0, because a hypothetical peptide with a nominal mass of 0 must have a monoisotopic mass equal to 0. The slope coefficient c 1 , determined by linear regression (cf. Methods) equalled 4.98·10 -4 ( Figure 1, Panel A -red dashed line), which is a value similar to the values 4.95·10 -4 reported by Wool and Smilansky [10].
We were interested in determining the dependence between monoisotopic and nominal mass analytically.
For example, the monoisotopic mass (m (M) ) of hypothetical peptides built only of one amino acid i can be pre-

The influence of the digestion enzyme on the wavelength of peptide mass clusters
In case of a complete sequence specific cleavage of proteins, the number of generated peptides is C P + 1 peptides, given that C P is the number of cleavage sites per protein.
The peptides generated from the terminus of the protein (further called terminal) will not bear a cleavage site residue R C at their end. All the other peptides, which we call internal, will have such a residue at their end. The fraction of the internal peptides f c,n is given by where n is the number of missed cleavages per protein. We approximate C P , for a sequence database, by: where are the relative frequencies of the cleavage sites and |P| is the average protein length in the database. The fraction of the terminal peptides in case of n missed cleavages is given by 1 -f c,n . The fraction of cleavage site residues R C in a internal peptide of mass m pep , with n missed cleavage sites is denoted f m,n and approximated by: where is the average mass of an amino acid residue. A more accurate model of f m,n is provided in the Appendix. In the case of terminal peptides the fraction of cleavage site residues R C equals f m,n -1 . The fraction of all the other amino acid residues R\R C equals 1 -f m,n or 1 -f m,n -1 respectively. Table 2 summarises these results.
In the case of internal peptides, the average contribution of the amino acid residues to the peptide mass is the weighted sum: where is the average mass of non cleavage residues, and: is the average mass of the cleavage site residues R C . Finally, the wavelength of internal peptides is presented as: The wavelength of terminal peptides was determined by: Cleavage probability p c In practice, the cleavage probability will depend on various factors, for example on the incubation time and the efficiency of the protease used. The probability to generate a peptide with n ∈ 0...∞ missed cleavage sites, given the cleavage probability p c can be modelled using the geometric distribution: holds. Hence, given the cleavage probability is p c and cleavage residues R C , we express the peptide mass by: where S n = (f c,n f m,n + f m,(n-1) -f c,n f m,(n-1) ). (20) Therefore, the wavelength λ of peptides if the cleavage probability is p c is given by: The monoisotopic mass as a function of the nominal mass can be expressed by: This equation represents our final model of the peptide mass cluster centres. To illustrate the accuracy of the prediction we computed the residuals Δ between the monoisotopic masses of the in silico database digest and the cluster centres predicted by Equation 24. Figure 3 shows the relative residuals Δ ppm (m) = Δ(m)/m·10 6 , in parts per million. The grey line shows the moving average of the residuals Δ ppm (m) computed for a window of 15Da.
R cleavage -frequencies of cleavage site residues; R non-cleavage -frequencies of non-cleavage site residues; f m,n -see Equation 6; f c,n -see Equation 4.
These coefficients are in good agreement with the slope and intercept determined by linear regression for the in silico sequence database digest ( Figure 1).
Furthermore, we observed that the intercept c 0 will be pos-

Error of the model
Combinatorial restrictions may cause significant differences between the linear prediction of the model (Equation 24) introduced and the actual location of the cluster centre. To asses this error we first computed the location of the cluster centres (average of all monoisotopic masses in cluster) of the in silico database digest, and afterwards determined the difference to the cluster centre location predicted by model of Equation 24. This difference (cluster) is shown in Figure 5.
For a moving window of 100Da we computed the maximum and minimum (orange), third and first quartile (red), median (blue) and mean(gree) of (cluster). The combinatorial restriction decreases with increasing mass and for peptide masses greater than 1000Da it is negligible. However, (cluster) increases again for masses greater than 2500Da because peptide masses may deviate more strongly from the cluster centres and furthermore much fewer long peptides are generated.

The type of distribution around the cluster centres
In order to remove non-peptide peaks prior to database search, filtering thresholds have to be chosen. In  To determine the type of distribution we use qqplots [28] shown in Figure 6. We compared the distribution of the residues Δ ppm (m), observed for four different mass windows (m ∈ (500 -530), m ∈ (1000 -1110), m ∈ (2000 -2200) and m ∈ (3400 -3700)) with the normal distribution and t-distributions with various degrees of freedom. The t-distribution with degrees of freedom μ ∈ (15, 25) is a good approximation of the empirical distribution of Δ ppm for masses > 2000,.

Sensitivity analysis
The input parameters to the model of the peptide mass cluster centres included: • f i -frequencies of the amino acids.
• cleavage specificity of the protease R C • |P| -Protein length • p c -cleavage probability To examine how the output of the model is influenced by these factors we varied the protein length |P| in steps of 100 from 300 to 800 amino acids per protein. We determined the amino acid frequencies f i for 9 sequence databases (cf. Methods) and used them as inputs to the model. Furthermore, six cleavage specificities (shown in Table 3 where c 0 , c 1 can be obtained using the Equations 27 and 26, m N is the nominal mass (an integer).
Therefore, for a given m M , c 0 and c 1 we can determine the deviation Δ from the closest cluster centre of smaller mass by using the modulo operator as suggested by Wool and Smilansky [10]: However, in order to determine the distance to the closest cluster centre we considered two cases: The units of Δ λ (m i , 0) are in [m/z]. The magenta dot dashed curves in Figure 3 indicate the maximum detectable distance from cluster centres in ppm (±0.5Da/ m·10 6 [ppm]). Deviations from the cluster centres outside the range enclosed by these two curves are assigned to the wrong cluster. In case of theoretical peptide masses and experimental masses calibrated to high precision, such distances are observed only for masses greater than 2500Da. Fortunately, the majority of tryptic peptide masses detected in a mass spectrometric peptide fingerprint experiment are below this mass.

Linear regression on peptide mass rule LR/PR
The limitations of calibration methods based on the property of peptide mass clustering are a mass accuracy of only 0.2Da, its sensitivity to non-peptide peaks in the spectra, and that it completely fails if the number of peptide peaks in the peak list is small [10,14,19]. Hence, in practice, the method is used to confirm the results of internal calibration only [14,29]. However, the advantage of the calibration methods based on the property of peptide mass clustering, over other calibration methods [12], is that no internal or external calibrants are required in order to calibrate the peptide mass lists.
We propose here a novel method for the calibration of PMF data, based on robust linear regression and the distance measure introduced in the Equation 30. To determine the slope of the mass measurement error we computed the deviation from the peptide mass rule for every pair of peak masses (m i , m j ) within a peak-list, employing the following equation:    To determine the intercept coefficient of the mass measurement error we subsequently computed Δ λ (m corrected , 0) (using Equation 30), for all peak-list masses. Strip-charts of the data-set for a mass range of 2210 -2212Da and 842 -843Da, including the tryptic autolysis peaks 842.508Da and 2211.100Da. Gray triangles -raw data; blue "+" -Wool Smilansky algorithm (cf. Appendix); red "o" -LR/RP algorithm for tryptic peaks .
brate mass spectrometric peak-lists to an accuracy of 0.1Da. This measurement accuracy surpasses the other published calibration methods [10,19] at least two-fold.
Filtering of non-peptide peaks using the peptide mass rule Non-peptide peaks can be recognised according to their deviation from the cluster centres. The amino acids that Schema of non-peptide mass filtering where λ DB = 1.000511 (Equation 2). We used the relative deviation of Δ ppm from the cluster centre in parts per million instead of using absolute values. Figure 3 shows that only very short peptides approach the lower bound of -413ppm. This is due to the low frequency of Cysteine (C). The high frequencies of K, L, I (whose λ ≈ 1.00074) mean that the theoretical upper bound of 241ppm can indeed be reached by some peptides with a mass of ≈ l000Da. Peptides of higher mass never approach the upper and lower theoretical bound due to the rapidly decreasing probability to consist of K, L or I, or of C only. The lines for the standard deviation of S N (orange lines) and of the 1% and 99% quantile (green lines) in Figure 3 indicate that it is an exceedingly rare event to encounter a peptide mass for which (m, 0) will deviate more Δ λ ppm Δ λ ppm Scatter plot : abscissae -peptide mass m i , ordinate -m i modλ with λ = 1.000495

Figure 10
Scatter plot : abscissae -peptide mass m i , ordinate -m i modλ with λ = 1.000495. In red are highlighted peaks removed from the dataset because of their high frequencies. In green, peaks removed due to the strong deviation from the peptide mass cluster centres.
than 200ppm from the peptide cluster centre predicted by our model. Therefore, we use 200ppm as a filtering threshold. An essential requirement, to apply this filtering method successfully is that peak-list must be calibrated to high precision [12].  [12] are shown in red.
We studied how the non-peptide peak filtering influences the Probability Based Mascot Score (PBMS) [22]. In theory, for example one cystein rich peptide strongly deviating from the peptide mass rule and with a unique mass in the database digest, if properly assigned is sufficient to identify the protein unambiguously [10]. In case of PBMS, which requires multiple matches to peptide masses, a single match of a unique peptide mass, even if properly assigned, will not give a score indicating reliable identification of the protein. Furthermore, this scoring scheme takes into account the number of non-matching peaks. If many unassigned peaks are observed, the score is decreased and the assignment is interpreted as insignificant. Therefore, the removal of non-peptide peaks should increase the identification sensitivity. Table 4 demonstrates that an increase of 2.5% in the number of identified samples can be obtained by removing all peaks with a distance (m, 0) > 200ppm from the peptide peaklists. Row 8 of Table 4 shows that non-peptide peak filter-ing increases the PBMS score in 30 -55% of cases.
Removal of peptide peaks due to filtering caused a decrease of the PBMS score in less than 1% of samples.
We concluded that non-peptide peak filtering increases the sensitivity of protein identification if using the PBMS scoring schema. However, to which extend these results can be reproduced is dependent on the database search algorithm used.

Conclusion
We introduced here a simple model to predict the cluster centres of peptide masses. The input parameters of the model can be easily determined for the sequence databases. We studied how these parameters influence the location of cluster centres, concluding that the cleavage specificity of the enzyme used for peptide digestion and the cleavage probability are the main factors. The change of the cluster centre location due to changes in average protein length or due to variability of amino acid frequencies among the databases is relatively small. However, our analysis also illustrates that, due to combinatorial constraints, the location of the cluster centres for masses smaller than l000Da can differ from the average location. Based on the model of the peptide mass cluster centres we derived a measure to determine the deviation of an experimental peptide mass from the nearest cluster centre. We used this distance measure to calibrate the peptide peaklists and to recognise non-peptide peaks. The calibration method, linear regression on peptide rule, is a robust and accurate method to calibrate single peak lists without resorting to internal calibrants. With this method higher calibration precision was obtained in comparison to other calibration methods, which also employ the property of peptide mass clustering. Columns: Arabidopsis t., Rhodopirelulla b., Mus musculus -peptide mass fingerprint datasets (cf. Methods). Row 1 -number of samples with a significant PBMS score prior to filtering of non-peptide peak masses. Row 2 -number of samples with a significant PBMS score for peak-lists with non-peptide removed. Row 3 -relative change of the identification rate (Row 2 -Row 1)/Row1 100. Row 4 -Total number of samples which produced a PBMS score. Row 5 -number of samples for which an increase of the PBMS score due to non peptide peak filtering was observed. Row 6 -number of samples for which no change of the PBMS score due to non-peptide peak filtering was observed. Row 7 -number of samples for which a decrease of the PBMS score due to non-peptide peak filtering was observed. Row 8-9 -relative increase and decrease of the PBMS score, respectively.
The same distance measure was used to recognise nonpeptide peaks and to remove them from the peak-lists. Due to their removal, an increase of the identification rate of up to 2.5% for the PBMS scoring schema was observed.

Data sets
In this study, we used three data sets generated in different proteome analyses: All PMF MS spectra derive from tryptic protein digests of individually excised protein spots. For this purpose, the whole tissue/cell protein extracts of the aforementioned organisms were separated by two-dimensional (2D) gel electrophoresis [33] and visualised with MS compatible Coomassie brilliant blue G250 [32]. The MALDI-TOF MS analysis was performed using a delayed ion extraction and by employing the MALDI AnchorChip ™targets (Bruker Daltonics, Bremen, Germany). Positively charged ions in the m/z range of 700 -4, 500m/z were recorded. Subsequently, the SNAP algorithm of the XTOF spectrum analysis software (Bruker Daltonics, Bremen, Germany) detected the monoisotopic masses of the measured peptides. The sum of the detected monoisotopic masses constitutes the raw peak-list.

Calibration
In order to perform filtering of non-peptide peaks the dataset must be calibrated to high mass measurement accuracy. To align the dataset we used a calibration sequence [12] consisting of several calibration procedures.
First calibration using external calibration samples was performed in order to remove higher order terms of the mass measurement error [11]. Next, the affine mass measurement error of all samples on the sample support was determined by linear regression on the peptide mass rule introduced here. Subsequently, the thin plate splines were used to model the mass measurement error in dependence of the sample support positions to calibrate the spectra. Finally, the spectra were aligned using a modified spanning tree algorithm [12].

Mascot database search
Processed peak-lists were then used for the protein database searches with the Mascot search software (Version 1.8.1) [22], employing a mass accuracy of ± 0.1Da. Methionine oxidation was set as a variable and carbamidomethylation of cysteine residues as fixed modification. We allowed only one missed proteolytic cleavage site in the analysis.

Sequence databases
We determined the amino acid frequencies of the nine protein sequence databases listed in Table 5. Seven of these databases are organism specific subsets of the NCBI non-redundant protein database [34].

In silico protein digestion
The theoretical digestion of the protein databases was done with ProtDigest [35], a command line program taking a protein sequence database file in fasta format and cleavage specificities as input. Other optional input parameters included fixed as well as variable modifications and number of missed cleavages. The output file contains all theoretically resulting peptides with their corresponding masses.

Regression analysis
The complete tryptic insilico digest of the SwissProt [27] database generated more than 7 million peptides. In order to compute the slope coefficient we were sampling 500 times 10000 monoisotopic and corresponding nominal masses. For each sample we fitted the affine linear model with and without fixed intercept using linear regression. The slope and intercept coefficients in Figure 1 are the medians of these 500 samples.