Abstract
Quantitative strategies relying on stable isotope labeling and isotope dilution mass spectrometry have proven to be a very robust alternative to the well established gelbased techniques for the study of the dynamic proteome. Postdigestion ^{18}O labeling is becoming very popular mainly due to the simplicity of the enzymecatalyzed exchange reaction, the peptide handling and storage procedures, and the flexibility and versatility introduced by decoupling protein digestion from peptide labeling. Despite recent progresses, peptide quantification by postdigestion ^{18}O labeling still involves several computational problems. In this work we analyzed the behavior of large collections of peptides when they were subjected to postdigestion labeling and concluded that this process can be explained by a universal kinetic model. On the basis of this observation, we developed an advanced quantification algorithm for this kind of labeling. Our method fits the entire isotopic envelope to parameters related with the kinetic exchange model, allowing at the same time an accurate calculation of the relative proportion of peptides in the original samples and of the specific labeling efficiency of each one of the peptides. We demonstrated that the new method eliminates artifacts produced by incomplete oxygen exchange in subsets of peptides that have a relatively low labeling efficiency and that may be considered indicative of false protein ratio deviations. Finally using a rigorous statistical analysis based on the calculation of error rates associated with false expression changes, we showed the validity of the method in the practice by detecting significant expression changes, produced by the activation of a model preparation of T cells, with only 5 μg of protein in three proteins among a pool of more than 100. By allowing a full control over potential artifacts, our method may improve automation of the procedures for relative protein quantification using this labeling strategy.
Beyond the control of synthesis and degradation of mRNA, living organisms possess a wide variety of regulatory mechanisms not affecting mRNA levels, such as protein synthesis, degradation, posttranslational modification or control of sub cellular location. Largescale quantitative measurements at the protein level, also called quantitative proteomics, are thus complementary to the well established mRNAbased gene expression profiling techniques and will probably become one of the cornerstones of systems biology in the near future.
In recent years, alternative strategies to the classical twodimensional gel electrophoresis approaches based on the chromatographic separation of complex mixtures of proteasegenerated peptides, also referred to as “shotgun proteomics,” have been developed. Peptides eluting from an RP^{1}HPLC column are sprayed into a mass spectrometer that detects, isolates, and fragments specific peptide ions to obtain structural information, which is correlated against protein databases to identify the peptide sequence. Recently these approaches have incorporated stable isotopic labeling techniques to produce relative quantitative information besides protein identification (1). When a peptide sample labeled with a heavy isotope is equimolarly mixed with an unlabeled sample, peptide ions appear as doublets showing the same intensity separated by a masstocharge offset determined by the mass of the label and the charge state of the peptide. Peptides differentially represented in each sample show departures from the expected 1:1 ratio, allowing the quantification of changes in the corresponding protein levels. Stable isotopic labeling can be achieved chemically by modifying reactive peptide groups (for example, cysteine side chains) with sitespecific reactants (2, 3), metabolically by growing cells in culture media enriched in heavy isotopes such as ^{15}N or ^{13}Clabeled amino acids (4–6), or enzymatically by adding [^{18}O]water molecules into peptide C termini during or after proteolytic cleavage (7–10). Among these techniques, ^{18}O labeling is emerging as a powerful labeling strategy for quantitative proteomics applications. Unlike other strategies that are restricted to peptides containing specific amino acids, any peptide generated by proteolytic digestion can be labeled with H_{2}^{18}O (9–12).
Peptide labeling with ^{18}O tags is performed by digesting the proteins with a proper endoprotease, originally trypsin, in the presence of [^{18}O]water; this produces incorporation of two ^{18}O atoms at the Cterminal end of peptides. It has been shown that proteolytic ^{18}O labeling can be decoupled from protein digestion so that the endoprotease can be used in a separate step to label peptides after proteolysis has taken place. This procedure has the advantage that proteins can be kept in solution prior to digestion using adequate chaotropic buffers or surfactants, that labeling can be performed with a limited volume of H_{2}^{18}O, and that digestion and labeling conditions can be optimized separately (10). This postlabeling strategy is becoming very popular mainly due to the simplicity of the exchange reaction and the availability of fast and efficient peptide desalting procedures and has been adopted by most researchers using this technique to make relative quantification of complex protein mixtures in solution.
The ^{18}O labeling process introduces either 2 or 4Da mass tags depending on the number of Cterminal oxygen atoms exchanged. Because the efficiency of the exchange reaction is not always complete with all peptides, a rather complex isotopic envelope pattern is usually obtained due to the overlap of the natural envelopes of unlabeled, singly labeled, and doubly labeled peptides. This has limited the use of threedimensional ion trap mass spectrometers, which produce low resolution survey scans, for peptide quantification by this labeling method. However, ion traps can perform high resolution mass spectra (or “ZoomScans”) of selected ions over a limited m/z range, and we have demonstrated previously that by using these scanning modes in a linear ion trap and taking advantage of its high scanning speed it is possible to make accurate quantitative measurements without compromising the ability of this machine to perform high throughput peptide identification (13).
Described software applications for quantification from ^{18}O labeling data either use simple schemes assuming complete exchange to the doubly labeled species so that the ^{16}O/^{18}O ratio is directly computed as the intensity ratio of the monoisotopic peak and the peak located at +4 Da (14, 15) or implement quantification algorithms that basically follow the formulation of Yao et al. (8) and Zang et al. (16). These algorithms compute the ^{16}O/^{18}O ratio either by sequentially subtracting the contributions of the peak heights of the three overlapping isotopic clusters (17, 18) or by directly measuring the heights of the first, third, and fifth isotopic peaks and then applying an analytical formula that corrects the intensities for the cluster overlap (17, 19). Here we present an advanced quantification algorithm specifically designed to deal with peptides labeled by postdigestion ^{18}O exchange. Our method fits the entire isotopic envelope to parameters related with a kinetic exchange model, allowing at the same time an accurate calculation of the relative proportion of peptides in the original samples and of the specific labeling efficiency of each one of the peptides. By improving the quantification procedure and allowing a full control over potential artifacts, we think that our method is a significant step toward the complete automation of the method for relative protein quantification using this popular labeling strategy.
EXPERIMENTAL PROCEDURES
Preparation of a Model Mixture of Proteins—
A mixture of standard proteins purchased from SigmaAldrich was prepared by mixing the following proteins: bovine serum albumin (4 pmol/μl), chicken egg lysozyme (59 pmol/μl), horse cytochrome c (20.2 pmol/μl), horse heart apomyoglobin (5.9 pmol/μl), and bovine αlactalbumin (0.5 pmol/μl). 5 μl of this mixture were used for enzymatic digestion.
Preparation of Endothelial Cell Extracts—
Cells were grown at 37 °C at 5% CO_{2}. The EA.hy926 cell line (kindly provided by Dr. Antonio MartinezRuiz) was cultured in Dulbecco's modified Eagle's medium with hypoxanthine, aminopterin, and thymidine supplement, 20% fetal bovine serum, 100 units/ml penicillin, 100 μg/ml streptomycin, and 5 μg/ml gentamicin. For preparing cell protein extracts, confluent cells were scraped and resuspended in nondenaturing lysis solution (50 mm TrisHCl, pH 7.4, 300 mm NaCl, 5 mm EDTA, 0.1 mm neocuproine, and 1% Triton X100 plus protease inhibitor mixture), incubated in ice for 15 min, and centrifuged at 10,000 × g for 15 min at 4 °C. Supernatants were collected, and the protein content was quantified using the Bradford reagent (BioRad).
Preparation of T Cell Extracts—
These extracts were kindly provided by Drs. Montse Carrascal and Joaquín Abián and were prepared as follows. Mononuclear cells were purified from PBSdiluted blood using the Ficoll gradient centrifugation method as described by Boyum (20). To this end, FicollPaque (Amersham Biosciences) was used according to the instructions of the manufacturer. The purity of the mononuclear fraction was then verified by checking under the microscope the fraction of mononuclear cells among all cells observed; it was around the expected value of 95%. Activation was carried according to the following procedure: either 5 ml of PBS (negative control) or 5 ml of 10 μg/μl PBSdiluted αCD3 antibody (BD Biosciences) were added to T75 vials to coat the surface with antibody. After a 5h incubation, flasks were washed twice with PBS, and 55 × 10^{6} cells, diluted in 10% FCSsupplemented RPMI 1640 medium, were added to each vial. After 15 min of incubation, cells were washed twice with PBS and centrifuged at 500 × g for 10 min. Cell pellets were stored at −80 °C until use. To extract proteins, lysis buffer (6 m urea, 50 mm Tris, 10 mm DTT) was then added. After alkylation with iodoacetamide (55 mm), protein extracts were acetoneprecipitated prior to digestion.
Protein Digestion and Postdigestion Labeling—
Total cell extracts from endothelial cells were subjected to cold acetone precipitation, lyophilized to dryness, and suspended in 8 m urea, 25 mm ammonium bicarbonate. 200 μg of starting protein material were used for cellular extracts, and 5 μl were used for the model mixture of proteins. Proteins were reduced for 1 h in 10 mm DTT, 25 mm ammonium bicarbonate, pH 8.25, at 37 °C and then alkylated for 45 min in 50 mm iodoacetamide. Digestion was carried out overnight at 37 °C with trypsin (1:50 protease/protein). In the case of the T cell extract, it was digested upon reconstitution in 25 mm ammonium bicarbonate, 2 m urea without further treatment. After acidification with 1% TFA, the digest was lyophilized. The dried pellet was then resuspended in 50 μl of 0.1% TFA and desalted in a Vydac C_{8} RP column (2.1 × 25 mm), using a Beckman Gold HPLC system, by onestep elution with 80% acetonitrile containing 0.1% TFA. The clean peptide pool was then dried down and dissolved in 40 μl of 20% acetonitrile (v/v); 40 μl of 100 mm ammonium acetate, pH 6.75, was added together with 40 μl of 50 mm calcium chloride and 10 μg of sequence grade trypsin; and the resulting mixture was dried in a vacuum centrifuge. The sample was then resuspended in 40 μl of either H_{2}^{16}O or H_{2}^{18}O (95%, SigmaAldrich) containing 20% acetonitrile and incubated at 37 °C for 48 h. Labeling was stopped with 1% formic acid.
Mass Spectrometry—
Protein digests were analyzed by LCESIlinear ion trapMS/MS using a Surveyor LC system coupled to a linear ion trap mass spectrometer model LTQ (ThermoFinnigan, San Jose, CA). The separation column was a 0.18 × 150mm Biobasic RP column (ThermoHypersilKeystone) operating at 1.5 μl/min. Peptides were eluted using 300min gradients from 5 to 40% solvent B (solvent A: 0.1% formic acid; solvent B: 0.1% formic acid, 80% acetonitrile). The linear ion trap was operated in datadependent ZoomScan and MS/MS switching mode using the three most intense precursors detected in a survey scan from 400 to 1600 amu (three microscans). ZoomScan mass windows were set to 12 Da to allow monitoring of the entire ^{16}O/^{18}O isotopic envelope of doubly and triply charged peptides irrespective of which isotopic peak was chosen as the precursor. Singly charged ions were excluded for MS/MS analysis. ZoomScan settings (maximum injection time, 50 ms; zoom target parameter, 3000 ions; and the number of microscans, 10) were exactly as optimized in a previous work (13). Normalized collision energy was set to 35%, and dynamic exclusion was applied during 3min periods to avoid fragmenting each ion more than twice.
Database Search—
Protein identification in human International Protein Index (IPI) version 3.12 database, downloaded from the European Bioinformatics Institute website, was carried out as described using SEQUEST. Database search parameters included as variable modifications Met oxidation and Lys and Arg + 4Da labeling and a fixed carbamidomethylation of Cys. Calculation of error rates of peptide identification was carried out as we have proposed earlier (21). Only ZoomScan spectra corresponding to peptide matches with a false discovery rate lower than 5% were used for quantification.
Kinetic Model for Trypsincatalyzed ^{18}O Labeling—
Let's consider a general case where peptides are subjected to a certain labeling reaction where two labeled species are produced using a reagent whose concentration is in large excess. Let's also assume that the labeled species may return to the unlabeled state and that the labeling reaction takes place according to the following kinetic mechanism, where B_{0} stands for the concentration of nonlabeled peptide; B_{1} and B_{2} stand for the concentration of mono and dilabeled peptide, respectively; k is the kinetic constant; and α and p are constants. A kinetic analysis of this equation (see the supplemental information) reveals that the concentrations of the three species along the labeling reaction are mathematically related in a manner that does not depend on the time or the kinetic constant.
In the case of trypsincatalyzed H_{2}^{18}O incorporation into peptides, the parameter p is equal to the fraction of water molecules containing ^{18}O, i.e. to the purity of H_{2}^{18}O. Besides oxygen exchange takes place in one of the two atoms at the C termini; hence only half the fraction of monolabeled species will give rise to a dilabeled species or return to the nonlabeled state. Therefore, α = 0.5, and the relation between the three species also becomes independent on the constant p (see the supplemental information). This relation may be expressed as a function of the total peptide concentration (B = B_{0} + B_{1} + B_{2})) and the labeling efficiency f.
Note that it would also be possible to derive Equations 1–3 by assuming that the exchange of the two oxygen atoms is independent of each other. In other words, when α = 0.5, the labeled oxygen atoms become homogeneously distributed among the mono and dilabeled peptide species even though a true equilibrium between these species has not actually been reached.
ZoomScan Spectra Preprocessing and Initial Quantification Step—
ZoomScan spectra are base linecorrected by subtracting from all points the median intensity of all local minima. A smoothened spectrum is then obtained by applying a moving average window, which is weighted using a tricubic function to better preserve peak heights. An initial estimate of relative isotope amounts is carried out on the basis of the formulation described by Zang et al. (16). The amounts of nonlabeled peptide (A) and of mono (B_{1}) and dilabeled peptide (B_{2}), expressed in the same units as the ion intensities determined by the MS detector, are determined from where I_{0}, I_{2}, and I_{4} are the raw intensities of the first, third, and fifth peak maxima of the isotopic cluster and M_{0}, M_{2}, and M_{4} are the relative proportions of the first (monoisotopic), third, and fifth peaks in the natural isotopic envelope. They are calculated from the elemental composition using a combinatorial subroutine determined from the peptide sequence. Isotopic abundance values are those reported by Chapman (22). The monoisotopic mass is also determined from the peptide sequence. Peptide sequence and charge state are taken from SEQUEST outputs once they are statistically analyzed to determine true peptide assignations.
Refinement of Quantification Using an Envelope Shape Model and Calculation of Standard Errors—
The initial estimates of the relative content of the three isotopic species mentioned above (A, B_{1}, and B_{2}) are used as starting parameters for a procedure where the isotopic envelope profile of the ZoomScan spectrum is fitted to a theoretical curve. Fitting is performed by a nonlinear, NewtonGauss unweighted least squares iterative method (23). The fitted function is the area covered by the first eight peaks of the signal doublet, which is modeled as a sum of three peptide isotopic envelopes, each one composed of a set of peaks evenly spaced 1/z units. To model peptide peaks, we observed that Gaussian distributions did not explain adequately the behavior of ions upon resonance ejection, which showed a variable degree of leptokurtosis. For this reason, peaks were modeled using mixed Gaussian/double exponential distributions where I(x) is the intensity at m/z x for a peak whose area is the unity, μ is the location parameter (the mean in the Gaussian distribution), σ is the scale parameter (the standard deviation in the Gaussian distribution), and β evaluates the proportion of double exponential component within the mixed distribution; the two last parameters are individually fitted for each spectrum. Five variable parameters are used to fit the function to the isotope envelope: A, B, f, σ, and β. Standard errors for these parameters were calculated from the matrix of partial derivatives as described elsewhere (23).
A Modified Algorithm Including a Correction for Incomplete Labeling Efficiency—
The corrected method was identical to the original one except that instead of fitting the parameters A, B_{1}, and B_{2} the procedure fits the parameters A, B, and f. ZoomScan spectra were modeled, as in the original protocol, as a sum of three isotopic envelopes; the first one was computed as the sum of A and B_{0} species (see Fig. 1), and B_{0}, B_{1}, and B_{2} were calculated from B and f using Equations 1–3. The modified algorithm used f = 1 as starting value for the curve fitting. This procedure allows a direct determination of corrected A and B proportions and of the labeling efficiency of each one of the peptides.
The theoretical relationship between the ratios determined by the standard (R_{S}) and corrected (R_{C}) algorithms and the labeling efficiency may be calculated considering the following. and Applying Equations 1–3 gives Equation 10. Using this equation it is possible to estimate the effect of labeling efficiency on the ratio determined using the standard method.
Statistical Analysis of Differential Expression Events—
Expression ratios were statistically analyzed by a method similar to that described by Li et al. (24). Briefly a ratio histogram in a log_{2} scale was constructed, and the distribution was fitted by least squares to a Gaussian curve, determining the mean and the S.D. A standard twotailed z test was then used to determine probability associated to peptides showing an expression ratio significantly higher or lower than the mean, which is usually, but not always, centered on zero. By introducing a sampledependent normalization, this method corrects for any systematic errors introduced during sample handling (24). To calculate the final p values, the standard error of the fitted parameters (i.e. the experimental error performed in the calculation of the ratio) was also taken into account by using error propagation theory (25, 26); this compensated for the fact that a noticeable departure from the center of distribution may not be statistically significant if the error made during the determination of the ratio (the fitting procedure) is sufficiently high.
The final statistical significance of potential differential expression events was assessed by determining the false discovery rates (FDRs) (27, 28), defined as the proportion of peptides expected to pass the p threshold, calculated from the fitted Gaussian distribution and the total number of determinations, among the observed number of peptides actually passing this threshold.
RESULTS
A Standard Method for Quantification from ZoomScan Spectra Obtained by Linear Ion Trap Mass Spectrometry—
In a previous report we described a stable isotopic labeling quantitative proteomics method based on ^{16}O/^{18}O labeling and linear ion trap mass spectrometry that takes advantage of the high scanning speed of this machine. Peptides identified in survey scans are subjected to a high resolution scanning mode (ZoomScan), which is used to quantify the isotopic partners, and then subjected to MS/MS analysis. This allowed both identification and relative quantification of peptides (13). Quantification was performed by decomposing the peak heights of the isotopic envelope into the relative contributions of unlabeled (A), singly labeled (B_{1}), and doubly labeled (B_{2}) peptide species using a formula described earlier (16). In this method all nonlabeled species are assumed to come from the nonlabeled sample (A), and the proportion of peptide coming from the labeled one is determined by adding up the proportion of singly labeled and doubly labeled species (i.e. B = B_{1} + B_{2}). The ratio is then computed as A/(B_{1} + B_{2}).
In this work we developed an improved algorithm that determines the same parameters by curve fitting of the isotopic envelope obtained in ZoomScan spectra. This method is based on the same assumption, uses the results obtained by our previous method as initial estimates, and makes a curve fitting to the peak envelope by using a peak shape model as explained under “Experimental Procedures” and exemplified in Fig. 1, A–C. Among other advantages, this algorithm, to which we will refer to as the standard method, allowed a more accurate quantification of proportions and also calculation of errors of estimates from curvefitting residuals.
Quantification Bias Associated to Incomplete ^{18}O Exchange—
To test the performance of the standard method in the practice, a negative control was prepared by making a comparative expression analysis of two identical aliquots of a peptide digest from a crude proteome extract from endothelial cells. The two samples were labeled with H_{2}^{16}O and H_{2}^{18}O water, respectively, and mixed in equal proportion, and their relative proportions were analyzed. The distribution of ratio values was analyzed using a base 2 logarithmic scale, which is a common practice for gene expression data and is expected to produce a symmetric distribution tightly centered on zero. We observed that the distribution of log_{2} ratios was not correctly centered on the expected value but was significantly biased toward the nonlabeled sample. Besides it was significantly asymmetric and showed an extended right tail (Fig. 2A).
After a close inspection of isotopic envelopes of peptides in the labeled sample, we noticed that this distortion was due to the presence of a small proportion of peptides having a low but significant proportion of nonlabeled species. To analyze a possible oxygen backexchange, labeled BSA tryptic digests were left at room temperature for several days in the presence of nonlabeled water containing 1% formic acid; in these experiments we did not observe any significant increase in the proportion of nonlabeled species. However, the proportion of nonlabeled species in the labeled proteome samples mentioned above was found to be higher than 50% for some peptides even when a fast analysis was performed just after the labeling reaction. We concluded that the observed isotopic proportions are not due to oxygen backexchange but to an incomplete incorporation of [^{18}O]water in a subset of peptides during the labeling reaction.
This was an unexpected observation given the reported good labeling behavior of tryptic peptide digests of BSA and other widely used purified proteins (10, 13, 29). In further experiments with different proteomes, we obtained similar results, corroborating the fact that the performance of the labeling process cannot be extrapolated from simple protein mixtures to very complex samples, such as those derived from the analysis of whole proteomes. Using higher enzyme/substrate ratios or higher concentrations of organic solvents or diluting the peptides in these solvents before adding water and enzyme failed to avoid the presence of residual amounts of peptides labeled at a low extent. These results suggest that this resistance toward Cterminal oxygen exchange is more strongly dependent on an unknown set of peptidespecific structural constraints, which are difficult to be controlled by the user, than on the particular labeling conditions. This is particularly relevant in the case of complex samples where, among a large number of structural patterns, there is a certain probability of finding a subset displaying very slow exchange kinetics. Incomplete labeling of these particular peptides may produce deviations in the expected mean ratio, suggesting true differential expression events.
A Kinetic Model for ^{18}O Exchange—
In an attempt to overcome this problem, we analyzed the kinetic behavior of the labeling process. As described under “Experimental Procedures,” a mathematical analysis of the kinetic process revealed that the relative proportions of nonlabeled and mono and dilabeled species from any peptide and at any time along the labeling reaction are expected to relate to each other according to Equations 1–3 even when the purity of [^{18}O]water is low. According to this equation, it would be possible to predict the relative proportion of the three isotopic species for any peptide if the labeling efficiency f is known. Conversely it would be possible to estimate the proportion of nonlabeled species (B_{0}) from that of the other two species (B_{1} and B_{2}) and hence to correct for the lack of a complete incorporation of at least one ^{18}O atom.
To check the validity of the kinetic model and the predictions of Equations 1–3, a digested proteome extract from a preparation of T cells was dried down and subjected to trypsincatalyzed ^{18}O labeling under incomplete labeling conditions. This was done by performing the labeling incubation for 24 h only, a time that we have observed previously not to be sufficient to achieve a complete exchange for many peptides (13). The labeled peptides were analyzed by linear ion trap mass spectrometry, and the ZoomScan spectra were used for quantification using the standard method. Because in this experiment no nonlabeled sample was mixed (i.e. A = 0), this analysis allowed us to determine the relative amounts of the three isotopic species (B_{0}, B_{1}, and B_{2}) and from these the labeling efficiency f for each one of the peptides. In Fig. 3, the observed proportions of the different species (circles) were compared with the theoretical predictions given by Equations 1–3 (solid lines). As shown, the majority of the peptides had a low proportion of nonlabeled species, around or less than 10% (open circles); the mean labeling efficiency was around 0.7. As expected from our previous observations, a noticeable number of peptides showed a significant proportion of nonlabeled species, which in some cases reached 70% (Fig. 3). The relative proportion of the three species was in perfect agreement with the predictions of Equations 1–3 in all the peptides without any exception. The same results have been obtained in our laboratory after the analysis of other samples (not shown), strongly suggesting that the kinetic model is universally valid for this kind of labeling process, and hence predictions given by Equations 1–3 may be potentially used to correct for an incomplete labeling efficiency.
A Corrected Method to Calculate Labeling Efficiency and Compensate for Incomplete ^{18}O Labeling—
As indicated earlier, all established quantification methods assume that the proportion of peptides in the nonlabeled sample (A) corresponds to those of nonlabeled species, whereas those in the labeled sample (B) come from the mono (B_{1}) and dilabeled peptide species (B_{2}) (Fig. 1, A–C); this implicitly assumes that the proportion of nonlabeled B species (B_{0}) is negligible. However, from a kinetics viewpoint, peptide labeling proceeds dynamically so that the proportion of B_{0} decreases along the oxygen exchange process but actually never takes a completely zero value. It is then difficult to ascertain at what point of the kinetics can B_{0} be neglected in the quantification. Obviously it would be desirable to take always into account the proportion of nonlabeled B species in the calculation of peptide ratios. With this idea in mind, we considered that the proportions of the three B species follow Equations 1–3 and hence may be calculated from a single parameter f and devised a modified algorithm where instead of fitting the ZoomScan isotope envelopes to a function constructed with the parameters A, B_{1}, and B_{2} as in the standard method (Fig. 1, A–C) they are fitted to a function containing the parameters A, B, and f. The proportion of the three B components is calculated using Equations 1–3, and the intensity of the first isotopic cluster is fitted to the sum of A and B_{0} (Fig. 1, D and E). This method, referred to here as the corrected method, makes it unnecessary to introduce additional fitting parameters, thus maintaining the degrees of freedom and hence the goodness of fit of the previous method. And despite its simplicity, it allows correcting for the presence of nonlabeled B species and permits a direct and simultaneous calculation of the proportions of A and B and also of the specific efficiency with which each one of the peptides is labeled.
Performance of the Corrected Method for Equimolar Protein Mixtures—
The performance of the corrected method was tested in the practice by applying it to the analysis of simple, equimolar protein mixtures subjected to incomplete labeling conditions. A tryptic digest of a mixture containing five different proteins was divided into two identical aliquots, and they were subjected to trypsincatalyzed enzymatic labeling in either H_{2}^{16}O or H_{2}^{18}O for 24 h. The two samples were then mixed and analyzed by linear ion trap mass spectrometry. The resulting spectra were then used for relative quantification using the standard and corrected methods. In Fig. 4A, the ratios determined from peptides originating from the same protein were averaged and compared. As shown, by using the standard method only two proteins could be accurately quantified, and the ratios of the other three were clearly overestimated; besides these ratios, when considered at the peptide level, showed a large dispersion. In clear contrast, the corrected method took into account the bias produced by low labeling efficiencies, yielding ratios that in all the cases were very close to the unity; this was reflected in the experimental error of the average ratios, which were much lower than those obtained with the standard method.
To analyze to what extent the corrected method was able to compensate for low labeling efficiencies in the practice, a tryptic digest from a complex protein mixture was again divided into two identical aliquots and subjected to labeling in the same conditions of incomplete ^{18}O incorporation. The relative quantification results obtained by using the standard (empty circles) and the corrected methods (filled circles) were then compared by plotting the labeling efficiencies of each one of the quantified peptides versus the calculated ratios in a base 2 log scale (Fig. 4B). As shown, the ratios calculated using the standard method were gradually shifted toward increasing values as the labeling efficiency of the peptides decreased; this shift was in good agreement with that theoretically expected from the kinetic model according to Equation 10 (Fig. 4B, black line), suggesting that the labeling efficiency was accurately estimated by using the corrected method. Consistently the log_{2} ratios determined using the corrected method were centered on zero and showed no bias even at labeling efficiencies as low as 0.2. These results suggested that the corrected method was able to compensate for the labeling efficiency effect in a wide range of values when the ratios were around the unity.
The corrected algorithm was then applied to the analysis of results obtained with the endothelial cell peptide extract, which was labeled using the optimized conditions and was analyzed previously using the standard method (Fig. 2A). As shown in Fig. 2B, the corrected algorithm eliminated the bias produced by incomplete incorporation, generating a symmetric distribution of ratios centered on the 1:1 value, which was satisfactorily fitted by a Gaussian envelope. Statistical analysis of these results failed to provide any peptide ratio showing a significant deviation from the average, thus confirming the internal consistency of the method. The distribution of labeling efficiencies calculated by the corrected method together with a plot of efficiencies versus log2 ratios is shown in Fig. 2C.
Performance of the Corrected Method When Peptide Ratios Take Extreme Values—
In a further set of experiments we tested the performance of the method in the extreme situations of peptide ratios. For this purpose we used the same dataset used to test the kinetic model (Fig. 3), containing labeled peptides only (i.e. A = 0 and B = 1), and analyzed the results obtained by subjecting this sample to quantification using the corrected method. This served to determine to what extent the modified algorithm was able to correctly assign the signal coming from the first isotopic cluster to the B_{0} species. As shown in Fig. 5A, where the estimated labeling efficiencies (i.e. the efficiencies computed by the method) are plotted as a function of the real ones (i.e. the efficiencies calculated in Fig. 3), the parameter f was calculated with good accuracy but only when it was higher than 0.4; below this value, it was overestimated so that part of the intensity of the first isotopic cluster is erroneously assigned to the A component. Consistently when the estimated labeling efficiencies were plotted as a function of the ratios calculated using the standard method, the shift toward higher values was more marked than that predicted according to Equation 10 when labeling efficiencies were lower than 0.4 (Fig. 5B, compare empty circles with the curve). Similarly the corrected method was only able to make an efficient correction of the effect of labeling efficiency when it was higher that 0.4 (Fig. 5B, filled circles). These results indicated that in these conditions a minimum labeling efficiency was actually needed in the practice to make an efficient correction; when a corrected peptide ratio is low, it should only be trusted provided that the labeling efficiency estimated by the algorithm is at least higher than 0.5.
The algorithm was then applied to the analysis of a nonlabeled peptide digest (i.e. A = 1 and B = 0). In these conditions, errors in estimation of labeling efficiency did not have any appreciable effect on the calculation of B using the corrected method because the estimated value of B took in all cases a very low value close to zero. These results are illustrated in the histogram of Fig. 5C; as shown, 96% of the peptides were assigned an A value equal to or higher than 0.9, and the lowest value of A determined from this experiment was 0.8. Similar results were also obtained when analyzing peptides at a ratio of 76:1 (not shown). These results showed that the corrected method was also able to produce accurate values for peptides displaying a high ratio.
Although the labeling efficiency has no significant effect on the calculation of corrected ratios when A ≫ B, it should be noted that in the extreme case of a completely nonlabeled peptide it is mathematically impossible to distinguish between a very high ratio (A ≈ 1) and a very low labeling efficiency (f ≈ 0). In this situation, these two solutions produce the same curve fitting to the isotope envelope. The reason that the corrected algorithm produces consistent and reproducible results with A ≈ 1 is that the starting parameters for the curve fitting are estimated from a model that implicitly assumes that there are no nonlabeled B species (B_{0} ≈ 0), i.e. that labeling efficiency is high; hence the corrected method converges to the first solution. Although this may in theory confound the corrected algorithm, in the practice it only happens when the labeling efficiency is very close to zero, and this is a situation that we have never observed experimentally even when using incomplete labeling conditions. Besides this situation is easily avoided in the practice by checking that labeling efficiencies do not take values that are too low. In a situation like that of Fig. 4B, for instance, it would be enough to keep the efficiencies above 0.2, whereas in conditions with more extreme ratios, like those of Fig. 5, B and C, labeling efficiencies must be cautiously kept above 0.5. We conclude that our analysis gives solid evidence that the corrected method makes very efficient corrections for labeling efficiency in any range of ratios provided that a reasonable level of labeling is experimentally achieved.
Application of the Corrected Method to the Analysis of Differential Expression Events Induced by Stimulation of T Cells—
After the encouraging results obtained with controlled experiments where the peptide ratios were a priori known, we applied this method to analyze changes in the protein pattern of small amounts of a crude proteome from a preparation of T cells when they were stimulated with antiCD3 antibody. For this end, samples from control and stimulated preparations were adjusted for the total protein content (5 μg each), trypsindigested, desalted, and subjected to enzymatic labeling for 48 h. The stimulated sample was labeled with H_{2}^{18}O, whereas the control was incubated with nonlabeled water. The mixtures were then mixed and analyzed by linear ion trap mass spectrometry in only one HPLC run using a long gradient. MS/MS spectra were then subjected to a database search using SEQUEST, and the FDR of peptide identification was calculated using a previously published empirical method (21). 849 MS/MS spectra were encountered having a FDR of 5%; they allowed the identification of 570 unique peptides belonging to 235 unique proteins. Among the ZoomScan spectra corresponding to these peptides, 495 were of sufficient quality to allow quantification; these spectra corresponded to 317 unique peptides belonging to 119 different proteins.
When quantification data obtained using the standard method was analyzed, despite the long labeling times used, the distribution of log_{2} ratio values was again found to be significantly biased toward the control sample and showed an extended right tail (Fig. 6A). Analyzing the data using the corrected method, it was possible to estimate the labeling efficiency f of each one of the quantified peptides. A semilog plot of f versus the ratios obtained by the standard method (Fig. 6C) revealed that the vast majority of peptides were labeled with an efficiency around 0.75; this is a typical result using the labeling conditions optimized in a previous study (13), and at this efficiency the proportion of nonlabeled sample (B_{0}) only accounts for 6% of the total amount of B. However, Fig. 6C also revealed the presence of a subpopulation of peptides labeled with a lower efficiency whose ratios are shifted toward higher values. These peptides were responsible for the bias of the distribution toward higher ratios and the presence of a right tail. In contrast, when the corrected method was used, the bias in the distribution was efficiently corrected and appeared centered at the 1:1 ratio, indicating that the two samples were in fact correctly adjusted for their protein content, and the log_{2} ratio histogram could be fitted to a symmetric Gaussian distribution (Fig. 6, B and D) with no evidence of a right tail, indicating that it was produced by peptides with a low labeling efficiency and not by true differential expression events.
A representative example of how the correction for labeling efficiency is made in poorly labeled peptides by using the modified method is shown in Fig. 7, A–D. This figure shows the ZoomScan spectra of four different peptides corresponding to the same protein. The ratios obtained using the standard method were 3.19, 1.54, 1.24, and 1.16, respectively. As it can be seen, the first peptide, which shows an unusually low incorporation of ^{18}O with a labeling efficiency of 0.37 (Fig. 7A), yields a shifted ratio measurement that would have indicated a significant expression change. However, the ratio calculated using the corrected method is in much better accordance with that of the other ones (Fig. 7, B–D) despite that its isotopic envelope might resemble that of an unlabeled peptide. The corrected ratios were 1.48, 1.29, 1.07, and 1.05, respectively, and none of them reflects a significant change in expression. Another example is found in Fig. 1 where a triply charged peptide with a rather low labeling efficiency was quantified using the standard and the corrected methods.
Statistical significance analysis of expression changes was then applied taking into account both the Gaussian ratio distribution obtained using the modified method and the errors associated to ratio calculation by curve fitting by using error propagation theory. We used a criterion based on the FDR defined in this context as the proportion of peptides that by chance alone are expected to display a similar increase or decrease in their ratios among the population of peptides with a given p value. We believe that this criterion, similar to that used for peptide identification (21), is much more adequate than conventional ones based on plain p values because it takes effectively into account the total number of quantified peptide pairs and estimates the number of expected peptide pairs that by chance alone deviate from the distribution. As shown in Table I, accepting a FDR not higher than 5%, only six spectra corresponding to tryptic peptides were detected as displaying significant expression changes; they belonged to three proteins, actin, histone H4, and thrombospondin1. Two further peptides with an FDR below 12% and belonging to the same set of proteins were also considered significant. These proteins were identified by three additional peptides; although they did not pass the rigorous statistical criterion, all of them still displayed expression changes that were in good agreement with that of the other peptides belonging to the same protein (Table I). One of the histone H4 peptides was quantified in its oxidized form. During the analysis of ZoomScan spectra, we also observed one that suggested a very clear expression change but whose MS/MS spectra did not give a significant match in databases. This spectrum was subjected to de novo sequencing using the program DeNovoX 1.0, and the resulting sequence tags were subjected to a string database search yielding an actinlike sequence; manual inspection of the MS/MS spectra and SEQUEST database search using nonenzyme constraints finally confirmed that this spectrum belonged to a nontryptic peptide from actin (Table I). Representative ZoomScan spectra of some of these peptides are shown in Fig. 7, E–H; the spectra of the remaining peptides can be found in the supplemental information, which also contains the complete list of peptides quantified in this experiment.
Just for illustrative purposes, Table I also shows the results that would have been obtained using a common but less rigorous criterion; using as a threshold a p value lower than 0.01 (i.e. 99% confidence), eight additional peptides would have been detected as reflecting a significant change. Interestingly all of them belonged to different proteins, and four of these belonged to proteins that were also quantified by other peptides not showing significant expression changes. This observation supports the validity of using the FDR criterion to detect significant expression changes.
The expression changes detected by this method appear to reflect true differential expression events because the three proteins identified have been implicated previously in T cell activation and other functions. Thus, thrombospondin1, an adhesive secreted glycoprotein that mediates celltocell and celltomatrix interactions, has been shown to play a role in T cell activation by antiCD3 (30); interaction of this protein with CD47 has also been shown to mediate expansion of inflammatory T cells (31). Histone H4 gene activation has been related to cell cycle activation (32), and it has also been implied in chromatin remodeling in response to T cell activation (33). Similarly actin cytoskeletal rearrangement is known to occur upon T cell stimulation (34). In this regard, the nontryptic actin peptide is probably the result of a proteolytic event. Besides actin presents a rather large number of isoforms in databases, and the four identified peptides are not homogeneously distributed among them, making it impossible to assign unequivocally the observed changes at the protein level; in addition, each one of the peptides is present in at least one isoform that does not contain the other ones. Cytoskeletal rearrangements probably affect the several actin types in a different manner, thus explaining the differences in expression ratios observed at the peptide level that were not observed in any other protein in this experiment (see supplemental information). In conclusion, the method was able to detect specific expression changes in three proteins, which are consistent with the process of T cell activation, among a pool of more than 100 proteins that did not exhibit significant variation in their cellular content.
DISCUSSION
In the original methods for ^{18}O labeling, it was assumed that trypsincatalyzed exchange reactions proceed completely to the doubly labeled (+4 Da) species (14, 15). Later it was acknowledged that although complete exchange is not always achieved, complete incorporation of at least one ^{18}O atom was considered feasible in the practice. This allowed using quantification algorithms where the isotopic envelope of the peptide pairs was assumed to be composed by one isotopic cluster deriving from the nonlabeled species, or sample A, and another two clusters deriving from the mono (B_{1}) and dilabeled (B_{2}) species, which together compose the sample B (B = B_{1} + B_{2}) (16–19). To achieve a complete incorporation of at least one oxygen, several labeling protocols have been described. A number of groups have reported efficient stable incorporation through a variety of ways, each one pointing out a set of recommendations to ensure the efficiency of the labeling process. The variability in the protocols spans the postdigestion buffer pH, the initial solubilization with strong organic solvents (acetonitrile or DMSO), the addition of calcium chloride, or the inactivation of trypsin after labeling to avoid backexchange among others (10, 29, 35). To date, none of these methods is still considered a standard. Postdigestion labeling strategies relying on other proteases (9, 11, 12) add further complexity to the situation. Besides and due to the small difference in the masses of the two labeled and the nonlabeled species, to obtain accurate results, a minimum mass resolution is needed, limiting the use of popular mass spectrometers like ion traps. These technical problems have limited the universal application of this quantification method.
In a previous work we reported an optimized labeling protocol and demonstrated that this method may be used for peptide quantification by using linear ion trap mass spectrometry and performing high resolution scans in narrow mass ranges (ZoomScan) (13). In the present work, we tried to further automate this method by developing a peakfitting algorithm that takes into account all the information contained in the whole isotopic envelope and not just the heights of the first, third, and fifth peaks (Fig. 1, left panels). This algorithm was more robust than previous ones, and by assessing the goodness of fit, it allowed the evaluation of the accuracy of fit and hence the introduction of statistically relevant information about the estimation of ratios. However, when this improved algorithm was applied in the practice to the analysis of a complex mixture of peptides derived from real proteomes, we noticed that a small but significant proportion of the peptides had no complete incorporation of at least one ^{18}O atom even when optimized protocols were used, producing a significant bias in the ratio distribution and false positive changes in expression ratios. Apparently when labeling very complex peptide pools, the probability of finding specific peptides showing a particularly low rate of ^{16}O/^{18}O exchange increased with the number of different peptides present in the sample so that some of them may be spuriously found not to be completely labeled with at least one ^{18}O atom.
To deal with this situation, we considered that the labeling process should actually be treated as a kinetic reaction where each one of the peptides had its own specific rate and reached a certain labeling efficiency after the exchange reaction step. Consistently the peptide population is expected to display a certain distribution of labeling efficiencies so that the proportion of completely nonlabeled species should never be neglected; it may be more or less close to zero depending on the labeling efficiency of each peptide. When peptide labeling behavior was analyzed in the practice, we found that labeling of all peptides analyzed followed precisely a kinetic exchange model where the fraction of nonlabeled and mono and dilabeled species could be accurately predicted as a function of the labeling efficiency. This relation, given by Equations 1–3, was independent of kinetic exchange rates and hence of peptide sequence or catalytic rate of the enzyme. The existence of a fixed relation between the three B species (B_{0}, B_{1}, and B_{2}) allowed the development of a modified algorithm that by direct curve fitting simultaneously determined the correct proportion of the A and B species and also of the labeling efficiency f that better explained the isotope distribution (Fig. 1). This procedure maintains the degrees of freedom and overcomes the limitation of current methods; it is in theory of potential application in any point of the exchange reaction process even when low purity [^{18}O]water is used.
Our results also indicate that the modified method is accurate in the practice even when peptide ratios take extreme values provided that a minimum labeling efficiency, as calculated by the algorithm itself, is reached. Although peptides labeled with lower efficiencies could be accurately quantified, a conservative criterion is to use a minimum f value of about 0.4 for a trustable quantification. This means that the modified algorithm is able to efficiently compensate for a proportion of nonlabeled peptide species as high as 36%. Expression ratios for peptides with lower f values should be considered with caution or even rejected. Direct calculation of labeling efficiency is also extremely useful in the practice to determine the goodness of the experiment and the general confidence to interpret quantification of peptide pairs. For this purpose the efficiencies estimated by the algorithm itself may be used to ascertain the degree of labeling achieved in the experiment; the efficiency versus log_{2} ratios plots are particularly informative to inspect the validity of the corrected ratios. Some examples of these plots, obtained under different conditions of labeling, are presented in this work where it is shown how experimental bias related to labeling efficiency is effectively removed, and there is a complete absence of false positives when performing the quantification of complex peptide samples against themselves.
The validity of the modified method was tested in the practice by determining significant expression changes produced by the activation of T cells with only 5 μg of protein extracts. More than 300 peptide pairs, corresponding to 100 proteins, could be quantified, yielding a symmetric Gaussian distribution of ratios closely centered around 1:1. Although the width of the ratio distribution suggested that statistically significant changes in expression ratios even lower than 2fold could be accurately detected at the 95% confidence level, we used a more rigorous statistical criterion. First, we took advantage of the curvefitting algorithm, which allowed estimation of the error associated to calculation of the ratio in each determination; this error was taken into account in the calculation of the p value with respect to the null hypothesis. Second, we considered all the expression data obtained in the experiment as a whole and calculated the FDR, or proportion of peptides that were expected to display just by chance a similar or higher change in expression ratio, at the same p threshold. The importance of taking into account all these factors becomes evident by the following considerations. Using the conventional statistical significance level of 95%, i.e. a p threshold of 0.05, 43 peptide pairs would have been detected as having a significant deviation in the ratio, corresponding to the 19 unique peptides shown in Table I. However, the 5% of the total population of quantified peptide pairs used to construct the Gaussian distribution amounts to 25 peptides; this means that 58% of the observed changes are expected to occur by chance alone (Table I). In addition, if the propagation of fit errors and the labeling efficiency correction are also omitted, the number of spectra in the set of differential expression candidates rises to 87 among which 66.3% are expected to occur by chance alone. Obviously and despite that it is widely used in this kind of quantitative experiments, using a prefixed p value is not an acceptable criterion by itself and leads to a false impression of high sensitivity at the expense of error rates that may become very large. In this particular experiment, we found it necessary to lower the p threshold down to 0.003, i.e. to use a 99.7% confidence level, until we achieved a satisfactory maximum error rate of 5%. The validity of this procedure is validated by the fact that all the expression changes detected at the peptide level are coherent among peptides belonging to the same proteins and are assigned to proteins that are known to be associated with T cell activation.
In conclusion, our results show that by introducing a correction for labeling efficiency in the quantification algorithm, ^{18}O labeling is a good quantitative method for shotgun proteomics analysis of highly complex peptide mixtures. The method described here is not only more robust and accurate that others described previously, but it also provides a means to determine the labeling performance obtained in the experiment. We think that this is an important advantage in the practice because stable isotope dilution techniques for quantitative proteomics are still not widely implemented, and skilled users with experience in controlling all the experimental factors related to labeling incorporation and quantification software are required. We also think that our approach makes the ^{18}O labeling method particularly attractive due to its simplicity, its almost universal applicability, the reasonable tradeoff between performance and costs, and the fact that this may be performed in the practice without high resolution machines using linear ion traps.
Acknowledgments
We are grateful to M. Carrascal, J. Abián, and A. MartínezRuiz for kindly providing the cell protein extracts and to F. MartínMaroto and X. Du for helpful comments on differential equations.
Footnotes

Published, MCP Papers in Press, February 23, 2007, DOI 10.1074/mcp.T600029MCP200

^{1} The abbreviations used are: RP, reverse phase; FDR, false discovery rate.

* This work was supported by Comisión Interministerial de Ciencia y Tecnología Grants BIO200301805 and BIO200610085, Comunidad Autonoma de Madrid Grants GR/SAL/0141/2004 and S2006/BIO0194, by the Spanish Cardiovascular Research Network (RECAVA), and by an institutional grant by Fundación Ramón Areces to Centro de Biología Molecular Severo Ochoa. The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

S The online version of this article (available at http://www.mcponline.org) contains supplemental material.

‡ Present address: Centro Nacional de Biotecnología, Universidad Autónoma de Madrid, 28049 Cantoblanco, Madrid, Spain.

§ Present address: Biological Sciences Division and Environmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, Richland, WA 99352.
 Received June 8, 2006.
 Revision received February 16, 2007.
 © 2007 The American Society for Biochemistry and Molecular Biology