Comparative Assessment of Quantification Methods for Tumor Tissue Phosphoproteomics.

With increasing sensitivity and accuracy in mass spectrometry, the tumor phosphoproteome is getting into reach. However, the selection of quantitation techniques best-suited to the biomedical question and diagnostic requirements remains a trial and error decision as no study has directly compared their performance for tumor tissue phosphoproteomics. We compared label-free quantification (LFQ), spike-in-SILAC (stable isotope labeling by amino acids in cell culture), and tandem mass tag (TMT) isobaric tandem mass tags technology for quantitative phosphosite profiling in tumor tissue. Compared to the classic SILAC method, spike-in-SILAC is not limited to cell culture analysis, making it suitable for quantitative analysis of tumor tissue samples. TMT offered the lowest accuracy and the highest precision and robustness toward different phosphosite abundances and matrices. Spike-in-SILAC offered the best compromise between these features but suffered from a low phosphosite coverage. LFQ offered the lowest precision but the highest number of identifications. Both spike-in-SILAC and LFQ presented susceptibility to matrix effects. Match between run (MBR)-based analysis enhanced the phosphosite coverage across technical replicates in LFQ and spike-in-SILAC but further reduced the precision and robustness of quantification. The choice of quantitative methodology is critical for both study design such as sample size in sample groups and quantified phosphosites and comparison of published cancer phosphoproteomes. Using ovarian cancer tissue as an example, our study builds a resource for the design and analysis of quantitative phosphoproteomic studies in cancer research and diagnostics.

Introduction

Phosphorylation events mediated by oncogenic kinases are widely recognized as major drivers of tumorigenesis. Despite more than 50 kinase inhibitors approved for cancer treatment by the Food and Drug Administration (FDA) and European Medicines Evaluation Agency (EMEA),[1] immediate and acquired drug resistance remain a key medical challenge. Toward personalized treatments, genome and transcriptome analysis based detection of somatic mutations and mRNA changes has led to breakthroughs such as the identification of HER2-positivity that guides therapies with HER-targeting antibodies[2,3] or BCR-ABL fusion indicating a response to imatinib and its analogues.[4] Although these therapies significantly improve the survival rates, e.g., in breast cancer and leukemia, they suffer from significant initial and acquired resistance.[5−7] For most cancers and kinase-directed drug therapies reliable predictive markers of initial drug response and acquired resistance are largely missing. Pathway analysis in cancer cell lines highlights the potential to predict drug sensitivity[8] and calls for translation to patient tumor tissues. Measuring the phosphorylation events and networks directly targeted by kinase inhibitors in the patient’s tumor may offer direct access to tumor- and patient-specific alterations governing drug resistance. Comprehensive quantitative coverage of the cancer- and patient-specific tumor phosphoproteome may identify patient subgroups with common response or escape mechanisms and may open new avenues to precision oncology by revealing a specific kinase signature in individual patients. Antibody-based techniques such as immunohistochemistry (IHC) and reverse-phase protein arrays (RPPA) are used in diagnostics and research, but they are limited by the availability of reliable antibodies and the low quantitative accuracy and robustness of these assays.[9,10] Mass spectrometry (MS)-based phosphoproteomic methods that combine the latest generation of instruments[11−14] with advanced computational tools[15−18] allow us to detect and quantify thousands of phosphosites. Label-free quantification (LFQ), stable isotope labeling by amino acids in cell culture (SILAC),[19] and chemical labeling e.g. with isobaric TMT[20] or isobaric tags for relative and absolute quantitation (iTRAQ)[21] are widely used for quantitative (phospho)proteomics and have been characterized in detail regarding their performance for studies in in vitro systems such as yeast or mammalian cell cultures.[22−27] However, we do not yet know whether the results translate to clinical phosphosproteomics in tumors. Tumor tissue differs from cell cultures in many aspects that are relevant for sample preparation and MS analysis. For instance, tumor tissue is typically snap-frozen and must be powderized prior to lysis, whereas cultured cells can be directly taken up in lysis buffer. Tumor tissue also contains different cell types as well as extracellular matrix and is thus more heterogeneous than cultured cells.[28] In the last 5 years, a series of studies have used LFQ, spike-in-SILAC, and TMT-based quantification separately for tumor proteomics[29−34] or phosphoproteomics[35−37] highlighting the need for a comparative assessment of quantitative phosphoproteomic methodologies in tumor tissue. The classical SILAC approach requires full metabolic labeling of the entire proteome and has been mainly limited to cell culture analysis with a few exceptions.[38] Mann and colleagues introduced spike-in-SILAC,[39] whereby a SILAC-labeled reference cell line is used to generate thousands of isotopically labeled peptides as internal standards for tissue proteome quantification. Ovarian cancer is well-accessible to clinical proteomics and presents a high unmet medical need due to limited treatment options and short survival.[40] Several studies have compared the performance of different quantitative proteome methods. Hogrebe et al.[26] and Stepath et al.[41] systematically compared LFQ, SILAC, and TMT workflows in human cell lines but not tissue samples, while Itzhak et al.[42] applied LFQ, SILAC, and TMT to the proteome analysis in both cell lysate and mouse tumor tissue lysate; however, they did not analyze phosphoproteomes. Using ovarian cancer as an example, we compared the performance of LFQ-, spike-in-SILAC-, and TMT-based phosphoproteomics regarding accuracy, precision, and robustness toward variation of the matrix and the phosphosite abundance. We use the term robustness in this article as the relative independence of phosphosite quantification from the number of replicates, phosphorylation levels, and matrix effects. The SKOV3 cell line, derived from ascites of a human ovarian adenocarcinoma[43] was used to compare our results from tumor tissue to a matched cell line.

Experimental Section

SKOV3 Cell Culture and SILAC Labeling

Experimental details of the cell culture conditions can be found in the Supplemental Experimental Section.

SKOV3 Cell Line and Ovarian Cancer Tissue Lysis

The detailed procedure of sample lysis can be found in the Supplemental Experimental Section.

Dephosphorylation and Sample Pooling

For dephosphorylation of the samples, it was required to exchange the protein lysis buffer with a phosphatase reaction buffer (10 mM Tris-HCl, 5 mM MgCl2, 100 mM KCl, 0.02% Triton X-100, pH 8.0), as alkaline phosphatases are inactive in the urea-containing protein lysis buffer. For this purpose, lysates of light labeled cells and tissues were loaded on a centrifugal filter unit (AMICON ULTRA-15 15 ML-10 kDa, catalog no. UFC901024) and centrifuged at 3 000g, RT for 1 h. Subsequently, the centrifugal filter unit was washed five times with 5 mL of phosphatase reaction buffer. To remove the phosphatase groups of the proteins,[44,45] alkaline phosphatase (protein, TSAP, 100:1, w/w, FastAP Thermosensitive Alkaline Phosphatase, catalog no. EF0651) was added to the lysates. Samples were incubated at 37 °C for 1 h. After inactivation of the reaction at 74 °C for 15 min, the buffer was changed back to the protein lysis buffer, as described before using the centrifugal filter unit. Subsequently, the protein concentrations were determined using the BCA assay kit. The intensities of the phosphopeptides before and after phosphatase treatment were measured using LC-MS/MS. The dephosphorylation efficiency was calculated as the ratio of the total intensity of the phosphopeptides before and after phosphatase treatment. To achieve three sample groups with different phosphorylation levels, the samples with and without phosphatase treatment were pooled according to the scheme in Figure . The original nontreated lysates were mixed with the dephosphorylated samples lysate at ratios of 1:5 (1×), 1:1 (3×), and 1:0 (6×). Next, each sample was divided into six aliquots containing 1 mg for LFQ, 500 μg for spike-in-SILAC, and 300 μg for TMT, each. A total of 500 μg of the heavy SILAC labeled cell lysate was mixed 1:1 with 500 μg of the 1×, 3×, or 6× tissue samples or the 1×, 3×, or 6× of the light SILAC-labeled cell line samples, respectively (Figure , panel C).

Figure 1

Study design: Sample preparation and MS measurement. (Panel A) Proteins extracted from the tumor tissue and cell line were divided into two aliquots. One aliquot was dephosphorylated by alkaline phosphatase. To generate the spike-in-SILAC standard, the cells were cultured in heavy labeled SILAC media for 5 passages. (Panel B) To obtain samples with known phosphosite quantities for evaluating the quantitative performance, the original and dephosphorylated aliquots were mixed in a ratio of 1:5, 1:1, and 1:0, resulting in three sample groups with different phosphosite quantities (1×, 3×, and 6×). The three sample groups allow the comparison of three different fold changes: 2, 3, and 6 (FC2, FC3, FC6). (Panel C) The samples were analyzed by LFQ, spike-in-SILAC, and TMT 10 plex, with 6 technical replicates per sample. The TMT 10 plex included one reference sample for normalizing the batch variation. The reference consisted of a mixture of sample 6×, 3×, and 1× in a ratio of 1:1:1. The scheme of sample labeling is shown in the table. (Panel D) The acquired raw data were processed with the MaxQuant suite, and the statistical analysis was done in R.

Protein Reduction, Alkylation, and Digestion

Protein lysates were subjected to protein reduction, alkylation, and then trypsin digestion. The detailed experimental procedure can be found in the Supplemental Experimental Section.

TMT Labeling

For TMT labeling, the dried peptides were resolubilized in 849 μL of TEAB/ACN buffer (80% H2O, 10% 1 M TEAB, 10% ACN). The TMT (TMT, Thermo Scientific, TMT10plex Isobaric Label Reagent Set, cat. no. 90406) reagents were solubilized in 100% ACN to a final concentration of 100 mM. A volume of 150 μL of the TMT stock reagents were added to the peptides according to Figure and incubated for 1 h at RT. To prevent side reactions, hydroxylamine (Thermo Scientific, 50% hydroxylamine for TMT experiments, catalog no. 90115) was added to a final concentration of 0.25% [w/v], followed by incubation for 15 min at RT. Next, samples were pooled according to the scheme in Figure and incubated for an additional 15 min. To reduce the concentration of ACN below 5%, the samples were diluted 1:1 with 2% TFA, before dilution with H2O. Samples were desalted using SepPak tC18 cartridges (Sep-Pak tC18 1 cc Vac Cartridge, 200 mg Sorbent per Cartridge, 37–55 μm, cat. no. WAT054925). The desalting steps were performed as described before increasing the used buffer volumes to 2 mL. A total of 50 μg of the eluent was used for the LC-MS/MS analysis and 2950 μg for the following IMAC enrichment. All the samples were dried using a SpeedVac.

Phosphopeptides Enrichment

The phosphopeptides were enriched using the immobilized metal affinity chromatography (IMAC) method (High-Select Fe-NTA Phosphopeptide Enrichment Kit, catalog no. A32992) according to the manufacturer protocol. In brief, the lyophilized peptide samples were suspended in 200 μL of IMAC binding/wash buffer. After removing the bottom closure of the IMAC spin columns and unscrewing the screw caps, the columns were centrifuged at 1000g for 30 s to remove the storage buffer. For equilibration, the columns were washed twice with IMAC binding/wash buffer. A volume of 200 μL of the suspended peptide samples were loaded on the equilibrated spin columns, and the columns were closed with the screw caps. The resin was mixed with the sample gently until the resin was in suspension. The suspension was incubated for 30 min, with gentle mixing every 10 min. Columns were washed three times with binding/wash buffer and one time with HPLC grade H2O. The phosphopeptides were eluted by adding 100 μL of elution buffer to the column two times and centrifuging at 1000g for 30 s. Samples were dried in a SpeedVac.

LC-MS/MS Measurement and Raw Data Processing

For LC-MS/MS analysis, samples were injected on an ultrahigh performance nano liquid chromatography system (Dionex UltiMate 3000 RSLCnano, Thermo Scientific, Bremen, Germany) coupled to an Orbitrap mass spectrometer (Fusion, Thermo Fisher Scientific) with a nano electrospray source. All LC-MS/MS data were processed with MaxQuant[15] version 1.6.5 and Proteome Discoverer suits (Sequest, version HT). Detailed information on the LC-MS/MS separation and MS parameters can be found in the Supplemental Experimental Section. All the raw and MaxQuant processed data is available at ProteomeXchange under a PXD030450 identifier. The identified phosphosite and phosphopeptides with quantitative values are available in the Supporting Information (Table S1 and Table S2).

Data Preprocessing and Statistical Analysis

The preprocessed LC-MS/MS data were further processed and statistically analyzed using R 3.6.5.[46] Median scale normalization was performed for LFQ and spike-in-SILAC. The intensities of phosphosites and peptides quantified in the enriched and nonenriched samples were log2-transformed. For each quantified phosphosite/peptide, the median intensity of the corresponding phosphosite/peptide data set was subtracted, and the median intensity of whole corresponding phosphosite/peptide data set was added. The phosphosites intensity of TMT samples was normalized to the reference sample to remove the variance from batches,[47] before further normalization on the total protein intensity. For the correlation matrix heatmap, the correlation matrix of normalized intensity between each sample were first calculated and then plotted by ggplot2.[48] Linear regression model were applied to classify the sample groups and the performance of regression models are visualized by the ROC curve in ggplot2.[48] To visualize kinase substrate sites that were enriched among the quantified phosphosites, enrichment analysis of predictive kinases was performed. The kinase-substrate relationship was extracted from the PhosphoSitePlus[49] database (downloaded on 2020-09-16). For the sequence motif analysis of the quantified phosphorylation sites, the ggseqlogo[50] package was applied to generate the sequence logos. The probability of the amino acid around the identified phosphosites is shown in the sequence motif plot, with annotating the residue physicochemical properties by color. The probability was calculated from the frequency of the amino acids surrounding the phosphosites in the phosphopeptide primary sequence.[50] The barplot, distribution plot, and violin plots were plotted using the ggplot2[48] package. The VennDiagram[51] package was used to plot Venn diagrams. Since the fold-change (ground truth) is known, we are able to determine the accuracy and precision. Accuracy describes how closely the median fold-change meets the ground truth fold-change. Precision describes the variability of the fold-change and is expressed as standard deviation or variance. R script used for data analysis and visualization is available on GitHub at https://github.com/functional-proteo-metabolomics/quantitative_phosphoproteomics.

Results

Sample Preparation and MS Measurements

The snap-frozen ovarian tumor tissue was powderized and taken up in ice-cold lysis buffer. SKOV3 cells were cultured in full medium. After a wash step, ice-cold lysis buffer was added directly to the tissue culture plates, and the cells were scraped off and transferred to tubes. Samples were sonicated and centrifuged, and the protein concentrations were adjusted to 1 mg/mL. All steps were performed on ice. To obtain samples with known ratios of phosphopeptide quantities (Figure A, panel A), we dephosphorylated tissue and cell lysates containing 20 mg of protein with alkaline phosphatase. The dephosphorylation efficiency was assessed by LC-MS/MS. After alkaline phosphatase treatment, the overall intensity of identified phosphopeptides was reduced to 1.6% for SKOV3 cells and 3% for ovarian tumor tissue (Figure S1A). We mixed the nontreated lysate with the dephosphorylated lysate at ratios of 1:5 (termed in the following as sample 1×), 1:1 (3×), and 1:0 (6×) (Figure B). Thus, the sample with the lowest phosphosite quantity was 1×. The phosphosite quantity in the 3× sample was thrice as high, and it was 6 times as high in the 6× sample. This resulted in a known fold change of 2 (FC2) for the 6× versus the 3× sample. The 1× versus 3× samples exhibited a known FC of 3 (FC3), and the 1× versus 6× samples had a known FC of 6 (FC6). The known fold changes served as the ground truth based on which we assessed in the following the performance of spike-in-SILAC, TMT-, and LFQ-based phosphosite quantification. In order to assess the technical variability, we performed the experiments in six technical replicates, which were prepared separately from the lysate samples (Figure ). For this purpose, each sample was divided into six aliquots of 1 mg for LFQ, 500 μg for SILAC, and 300 μg for TMT.

To generate the spike-in-SILAC standard, the SKOV3 cells were cultured for five passages in medium containing heavy lysine and arginine.[52] For the cell lysate to be compared with tissue sample, SKOV3 cells were cultured in light SILAC media. The labeling efficiency was >98% and the arginine-to-proline conversion was 1% or below (Figure S1B,C). 500 μg of the heavy SILAC-labeled cell lysate was mixed 1:1 with 500 μg of the 1×, 3×, or 6× tissue lysates or of the 1×, 3×, or 6× light SILAC-labeled cell lysates, respectively (Figure B,C), resulting in a total amount of 1 mg for each sample. All samples were reduced, alkylated, digested with trypsin for 16 h, and desalted using reversed phase solid phase extraction. TMT 10-plex was used to label six technical replicates of the 6×, 3×, and 1× tissue and cell line samples, which we assigned to two batches (Figure C). Each batch contained in addition a reference sample that was composed of 100 μg of the 6×, the 3×, and the 1× sample. For TMT analysis, two TMT batches were prepared, each consisting of 10 TMT-channels composed of 9 samples and 1 reference sample Figure C (bottom). For each channel, 300 μg of TMT-derivatized sample was used. This resulted in a total sample amount of 3 mg per TMT batch. Thus, we stayed below the maximum binding capacity of the Fe-NTA (nitrilotriacetic acid)-based IMAC columns (5 mg). Before IMAC enrichment, 50 μg of each sample was set apart for total protein concentration analysis. The rest of the samples were subjected to phosphopeptide enrichment by IMAC. All samples were dried and analyzed by LC-MS/MS using an Orbitrap Fusion instrument. The gradient length was kept the same for all samples (2 h). The LFQ and SILAC samples were analyzed by data dependent acquisition (DDA)[53] tandem MS (MS2). The TMT samples were analyzed by synchronous precursor selection MS3 (SPS-MS3),[54,55] to minimize reporter ion cross-contamination from peptides coisolated for fragmentation in the same isolation window. The data were analyzed using the MaxQuant[15] suite and SequestHT. MaxQuant is one of the most broadly used search engines in the community, which is why we focused on MaxQuant for our analysis. Shotgun proteomics by DDA is a widely used technology.[53] Although DDA achieves fast and accurate large-scale proteome analysis, missing values of 50% or higher are inherent characteristics of the stochastic precursor selection.[56,57] Data preprocessing with match between runs (MBR) is a common approach to reduce the number of missing values in DDA and perform identification transfer by matching nonidentified single-stage MS (MS1) features to identified ones.[27,58] We compared data processing without and with the MBR[59,60] option for LFQ and spike-in-SILAC, including the requantify (req) option that is typically applied in SILAC experiments to extract the signal of nondetected peaks, which are complementary to successfully identified light or heavy SILAC singleton peaks.[15]

Phosphosite and Phosphopeptide Identifications

We compared the total measurement time, the number of recorded MS2 spectra and MS3 spectra, the number of identified MS2 spectra, and the number of identified phosphosites among LFQ, spike-in-SILAC, and TMT (Table S3). In TMT, the MS3 spectra were used for quantification and the MS2 spectra were used for identification. Therefore, we compared the number of acquired and identified MS2 spectra to assess the identification performance. The threshold for the phosphosite localization probability was recommended to 0.75.[59] We calculated the number of phosphosites identified with different probability thresholds of phosphosite localization (Table S4). As expected, the number of phosphosites identified increased with lower phosphosite localization probability thresholds, and the trends were similar for the different methods. However, the use of low thresholds has an impact on the error rate, and it must be critically weighed which error rate can be tolerated. We used the PTM-score for phoshosite localization implemented in MaxQuant, with the widely used default parameter of 0.75 phosphosite localization probability threshold.[61] Note that other phosphosite localization tools (Ascore, phosphoRS, pSite, or tools) based on MS/MS spectra prediction with mass shift introduction at phosphorylated amino acids may provide different phosphosite localization results.[62]

For tumor tissue, spike-in-SILAC provided the highest number of detected (nspectra = 613 676) and identified (nid = 240 688) MS2 spectra, followed by LFQ (nspectra = 471 216, nid = 171 752) and TMT (nspectra= 63 421 and nid = 8 122). However, the number of phosphopeptide identifications was higher with LFQ (8 932) than with spike-in-SILAC (6 923). The discrepancy between the number of MS2 spectra and the number of identified phosphopeptides may be explained by the spike-in-SILAC pair-induced increased complexity at the MS1 level. The isolation and fragmentation of the light and heavy phosphopeptide lead to a redundant identification, which reduces the number of identified phosphopeptides compared to LFQ.[63] TMT resulted in the lowest number of identified phosphopeptides (n = 4 281), which can be assigned to the longer duty cycle with SPS-MS3, resulting in a lower number of fragmented and identified MS1 peaks. However, comparison of the number of identified phosphosites per minute shows that TMT identifies over 10-fold more phosphosites per minute LC-MS/MS analysis time than LFQ and spike-in-SILAC. Reducing the sample complexity in TMT by prefractionation, e.g., by high-pH reversed-phase chromatography, followed by SPS-MS3 likely increases the number of identified phosphosites. Similar observations were made for cell lysate (Table S3).

In line with the number of identified phosphopeptides, LFQ provided the highest number of identified phosphosites (n = 5 812) followed by spike-in-SILAC (n = 4 282) and TMT (n = 3 578) for tumor tissue. Similar results were obtained for cell lysate (Table S3). In total, 1 640 (20%) and 1 205 (18.6%) of the phosphosites were identified by all three methods in the tumor and cell line samples, respectively (Figure A,B). Spike-in-SILAC and LFQ showed an overlap of approximately 40% (40.2% for tumor tissue, 43.2% for cell line). Lower overlap was observed between TMT and LFQ or spike-in-SILAC for both tissue (TMT-LFQ, 25.9%; TMT-spike-in-SILAC, 21%) and cell line (TMT-LFQ, 20.7%; TMT-spike-in-SILAC, 21.7%). This may be explained by TMT-derivatization due to which the resulting phosphopeptides differ in their chemical composition and physicochemical properties[20] from phosphopeptides analyzed by LFQ and spike-in-SILAC that are chemically equivalent and differ only in their isotopic composition. We compared the frequency of b- and y-fragment ions of all identified phosphosites (Figure C,D). LFQ and spike-in-SILAC exhibited similar fragment ion frequencies in both tissue and cell line samples, with a higher frequency of y-ions (up to 75%) than b-ions (up to 25%), as expected for nonderivatized tryptic peptides.[64] Similar relative intensity distributions of b- and y-ions were obtained in SPS-MS3 for TMT-derivatized phosphopeptides and as well as for the nonphosphorylated peptides (Figure S2A,B). The increased proportion of b-ions is likely related to the CID fragmentation in the ion trap.[65] Also the basic 2,6-dimethylpiperidine residue of the TMT-tag increases N-terminal fragmentation and b-ion formation according to the mobile proton model.[66] We assessed this idea by reanalyzing label-free, TMT-MS2, and TMT-SPS-MS3 phosphopeptide data from the study of Hogrebe et al.[26]Figure S3 depicts the y and b ion relative intensity distribution. We report that both the fragmentation method and the TMT label influence the y and b ion relative intensity distributions. The TMT label has a larger effect than the fragmentation types (HCD TMT-MS2 versus CID/HCD TMT-SPS-MS3). The higher frequency of the b-ions may be explained by the introduction of a tertiary amine at the N-terminus as a result of TMT labeling. The tertiary amine of the TMT-tag has a higher gas phase basicity than the primary amine of the N-terminus of a tryptic peptide. As a result, the proton at the piperidine ring of the TMT-tag linked to the peptide N-terminus is less mobile, and b- and y-fragment ions are generated at similar frequencies.[66] The TMT-label impacts both the ionization and the fragmentation of the derivatized phosphopeptide,[67] which may explain the low overlap of identifications in TMT with LFQ and spike-in-SILAC. Our results illustrate that a derivatization strategy can shift the identifications to different subpopulations of phosphosites and peptides, which needs to be considered when comparing data sets acquired without and with chemical labeling. To further characterize the different identified phosphosites, we assessed the number and percentage of identified phospho-serine, phospho-threonine, and phospho-tyrosine sites (Table S5). In tumor lysate, a higher percentage of phospho-threonine and phospho-tyrosine were identified in the TMT samples. In cell lysate, the proportions of phospho-serine, phospho-threonine, and phospho-tyrosine sites were similar for the three methods. In addition to MaxQuant, we analyzed the LC-MS/MS data with Proteome Discoverer (version 2.4) using SequestHT and Percolator to evaluate the number of phosphopeptides and phosphosites identifications. The number of phosphopeptides and phosphosites identified was similar to the results obtained with MaxQuant (Figure S4 and Table S6), suggesting that our findings can be replicated with other search algorithms. We used the data generated with MaxQuant for further comparison of the performance of different quantitative phospho-proteomic methods.

Figure 2

Comparison of phosphosite identifications by LFQ, spike-in-SILAC, and TMT in tumor tissue lysate or cell lysate samples. (A,B) Venn diagram showing the number of identified phosphosites for the three different quantification methods (LFQ, spike-in-SILAC, and TMT) in the tumor tissue lysates (A) and cell lysates (B). (C,D) Density plots showing for each quantification method the distribution of b- and y-ions of the identified phosphopeptides in tumor tissue lysates (C) and cell lysates (D). (E,F) Bar plots showing the number of phosphosites identified in each sample group (6×, 3×, and 1×) for the different quantification methods used in tumor tissue lysates (E) and cell lysates (F). The color intensity indicates the number of replicates in which the phosphosites were identified. Only phosphosites with a localization probability of at least 0.75 were considered for the analysis.

We next compared the reproducibility of phosphosite identifications across the six technical replicates in tumor tissue (Figure E). The LFQ and spike-in-SILAC data sets were processed with or without activation of the MBR feature in MaxQuant. The reproducibility of the phosphosite identification is reflected by the distribution of the numbers of identified phosphosites across the replicates. MBR can be considered as an identification transfer step and therefore influences the reproducibly of phosphosite identification. MBR has recently been implemented for TMT.[68] However, MBR cannot match MS1 peaks that have not been submitted to fragmentation but only fragmented but not identified MS1 peaks. Such peptides generally yield low-quality MS2 spectra with low-abundance fragments. The MS2 fragments are difficult to fragment further by SPS-MS3, resulting in low intensity reporter ions and poor quantification. Since the current version of MaxQuant does not provide sample specific reporter ion information for nonfragmented high-abundance MS1 signals and thus no sample specific quantification values, we performed the TMT data analyses without MBR and opted for SPS (synchronous precursor selection)-MS3 technology to enhance the quantification accuracy.[54,55] As expected, MBR increased the number of reproducibly identified phosphosites for both LFQ and spike-in-SILAC. The number of phosphosites identified across six replicates decreased with decreasing phosphosite quantity (i.e., 6× vs 3× vs 1×) for LFQ and spike-in-SILAC with and without MBR. This indicates that with LFQ and spike-in-SILAC, reproducible identification increases with phosphosite quantity and the resulting signal intensity. This was in contrast to TMT for which the number of phosphosites reproducibly identified across six replicates was constant over the entire range of phosphorylation quantity (n = 1593 in 6×, n = 1591 in 3×, and n = 1582 in 1×). This likely comes from the multiplexing of the 1×, 3×, and 6× samples due to which low intensity signals in the 1× sample can be identified based on the higher average intensity from the pooled 1×, 3×, and 6× samples.[20] The higher average phosphosite quantity and signal intensity in the pooled TMT sample also explains the better reproducibility of identifications for samples with low phosphosite quantity (i.e., 1× and 3×) as compared to LFQ and spike-in-SILAC. For the high phosphosite quantity (6×) TMT yielded also more reproducibly identified phosphosites than spike-in-SILAC, but the MBR option increased this number beyond that of TMT. LFQ without and with MBR yielded higher numbers of reproducible identifications than TMT for high phosphosite quantities (6×). The reproducible identification of phosphosites in six or three replicates in TMT (with other replicate numbers largely missing) arises from the fact that the TMT samples were measured in two batches with three replicates per batch (Figure C). Therefore, a given phosphosite was identified in either of the batches (i.e., in 3 replicates), or in both batches (i.e., in 6 replicates). Comparable results were obtained for the reproducibility of the identification of phosphopeptides, with higher overall numbers as phosphosites can be assigned to several peptides (Figure S5A).

The overall numbers of reproducibly identified phosphosites and phosphopeptides from cell lysate were on average lower than from tissue lysate (Figure F and Figure S4B). Also here, TMT yielded a similar number of reproducible identifications across the different phosphorylation quantities. TMT yielded a higher number of reproducible identifications than spike-in-SILAC without and with MBR and LFQ across all phosphosite quantities. Only LFQ with MBR in a sample with a high phosphosite quantity (6×) performed better than TMT in terms of reproducibly identified phosphorylations. In both matrices, tumor tissue and cultured cells, TMT exhibited better robustness than LFQ and spike-in-SILAC regarding reproducibility of identifications across the whole range of phosphosite quantities. In tumor tissue with high phosphorylation quantity, spike-in-SILAC and LFQ outperformed TMT regarding the number of reproducible identifications, but this advantage was lost for lower phosphosite quantities.

Reproducibility, Accuracy, and Precision of Phosphosite Quantification

We assessed the reproducibility of phosphosite quantification by calculating the coefficient of variation (CV) distribution (Figures S6A,B and S7A–D). We observed no major difference in CV distribution in data sets analyzed with or without MBR both for LFQ (Figure S7A,B) or spike-in-SILAC (Figure S7C,D). TMT provided the lowest CVs, followed by spike-in-SILAC+MBRreq and LFQ+MBR (Figure S6A,B). For TMT, the CV distribution was similar across phosphosite quantities (1×, 3×, 6×). Thus, TMT exhibits the highest reproducibility for quantification of low-abundance phosphosites. TMT multiplexing reduces the technical variability as multiple samples are pooled and processed together during sample preparation (e.g., phosphopeptide enrichment and desalting by solid-phase extraction) and LC-MS/MS measurements. A reference channel in each TMT-batch further reduced the measurement variability between the two batches. In spike-in-SILAC, the spiked-in heavy-labeled phosphopeptides serve a similar purpose as the TMT reference channel, namely, to reduce variability introduced during sample preparation and LC-MS/MS analysis. This likely explains the lower CV values with spike-in-SILAC as compared to LFQ, for which all experimental and measurement steps were conducted independently and no reference sample was spiked-in. TMT shows lower CVs compared to spike-in-SILAC, although mixing of the multiplexed samples in spike-in-SILAC occurs earlier than in TMT. The smaller CVs of TMT are most likely due to the peak quantification in smoothed MS3-orbitrap mass spectra, whereas spike-in-SILAC used MS1-based 2D peak quantification in an LC-MS map. In these, the chromatographic peaks have a higher noise level than the signals in smoothed MS3-orbitrap mass spectra.[69] This may enable a more precise quantification in TMT data as compared to spike-in-SILAC data.

The CV distribution reflects the variability based on the mean of six replicates for each phosphosite. To assess the agreement between data sets in a pairwise manner, we analyzed the linear correlation between two data sets by calculating the pairwise correlation of phosphosite quantities for phosphopeptides identified and quantified in both samples (Figures S6C–H and S7E–H). Similar to the CV plots, no major difference was observed in the correlation heatmaps with and without MBR/MBRreq for LFQ and spike-in-SILAC (Figure S7E–H). For tissue and cell lysates, both TMT (Figure S6G,H) and spike-in-SILAC (Figure S6E,F) showed Pearson correlation coefficients above 0.9 within a sample group for all phosphosite quantities (1×, 3×, 6×). In TMT, the correlation was always above 0.9 and did not change between different phosphosite quantities, although a batch effect was visible that broadened the correlation coefficient distribution. For spike-in-SILAC, lower coefficients were observed for low phosphosite quantities (1×). This was even more pronounced for LFQ, resulting in a Pearson correlation coefficient below 0.9 for tumor lysate with low-abundance phosphosites (1×, Figure S6C). For all comparisons between sample groups (i.e., 6× vs 3×, 3× vs 1×, 6× vs 1×), TMT yielded Pearson correlation coefficients above 0.9 (Figure S6G,H), reflecting the high reproducibility of TMT across different phosphosite quantities. For comparisons between sample groups, spike-in-SILAC and LFQ showed considerably lower Pearson correlation coefficients than TMT. The correlation coefficients were particularly low when samples with low phosphosite quantity (1×) were part of the comparison, reflecting the drop in reproducibility and higher error in the quantification of low-abundance peaks for spike-in-SILAC and LFQ as the phosphosite quantity decreased. In summary, TMT exhibited the lowest variability and highest correlation for tissue and cell lysates, especially when the phosphosites were at low abundance.

We evaluated the precision and accuracy of the determination of the fold-changes between samples with known relative phosphosite quantities (1×, 3×, 6×) (Figure ). We assessed the distribution of detected fold changes relative to the theoretical fold changes (FC2, FC3, and FC6) by violin plots (Figure A,B) and by mean squared errors (MSE)[41] of the sum of positive deviation or variance (Figure C,D), indicative of the quantification error in accuracy and precision, respectively. LFQ and spike-in-SILAC determined all fold-changes with high accuracy and yielded more accurate quantifications than TMT in tumor and tissue lysates (Figure A–D). The use of MBR decreased the accuracy (Figure C,D). MBR also reduced the precision of all determined fold changes, and this effect was more pronounced in spike-in-SILAC than in LFQ. Matched features with low abundance or false feature matching may be the reasons for lower precision and accuracy in MBR.[69] This shows that MBR trades off a smaller number of missing values for accuracy and precision. TMT showed the highest precision and the lowest accuracy, consistently underestimating the expected fold change (Figure A–D). Our results are consistent with the study from Hogrebe et al., showing that TMT performs best in terms of precision while showing lower accuracy.[26]

Figure 3

Evaluation of precision and accuracy of phosphosite quantification. (A,C) tumor tissue lysates. (B,D) cell lysates. (A,B) Violin plots showing log2-transformed fold changes to evaluate quantification precision and accuracy errors of the methods. The boxes show the first, second (median), and third quartile. Whiskers show the minimum/maximum value within the 1.5 interquartile range. Expected log2-transformed fold changes were highlighted by colored lines. (C,D) Bar plots showing mean squared errors of accuracy and precision. Mean squared errors were calculated as the sum of the square of positive deviation and variance for each method and all replicates.

We also investigated the performance in detecting significant differences in the abundance of phosphosites between the different fold changes. Phosphosites had to be quantified in at least three out of six replicates. Statistical analysis was performed using a two-sided t test and the Benjamini–Hochberg procedure to correct for multiple testing. Across all fold-changes, TMT determined the highest number of phosphosites as being significantly different for both tissue (Figure A) and cell (Figure B) lysates. The next-best results were yielded by LFQ+MBR and spike-in-SILAC+MBR. LFQ and spike-in-SILAC performed worst. As multiplexing in TMT reduces the issue of missing values, almost all identified phosphosites (>98% for all sample groups) could be used for differential analysis. As stated earlier, also quantification based on smoothed MS3-orbitrap mass spectra yields more precise quantification and lower CVs, likely enhancing the discrimination of small differences. Spike-in-SILAC and LFQ resulted in a considerable proportion of phosphosites identified in less than three replicates, which prevented them from being used for differential analysis. The use of MBR in LFQ and spike-in-SILAC increased the number of significantly different phosphosites, in particular for the higher fold changes (FC3, FC6), which could be detected despite the lower precision (Figure C,D).

Figure 4

Performance in detecting quantitative phosphosite differences. (A,B) Barplots showing the number of identified and quantified phosphosites in tumor tissue lysates (A) and cell lysates (B). Identified only: phosphosites identified in less than three replicates. Quantified nonsignificant: phosphosites identified and quantified in at least three out of six replicates, with FDR >0.05. Quantified significant: phosphosites identified and quantified in at least three out of six replicates, FDR ≤0.05. (C,D) Receiver operating characteristic (ROC) curves for evaluating the ability to diagnose the different phosphorylation quantities in tumor tissue lysates (C) and cell lysates (D). Zoomed-in ROC curves are presented in the lower right corner of each graph.

We calculated true-positive-rates (TPR) and false-positive-rates (FPR) of the differentially quantified phosphosites for each quantification method and visualized them in a receiver operating characteristic (ROC) curve[70] (Figure C,D). We also calculated the area under the ROC (AUROC), which provide classification performance of TPR and FPR for whole threshold ranges (Table S7). At an FPR threshold of 0.05, the TPR was above 0.95 for all quantification approaches and the largest fold-change (FC6), indicating that all methods were able to correctly determine a fold-change of six (FC6) between the phosphorylation quantities (Figure C, Table S7). These results were also reflected by the AUROC values showing value higher than 0.98 for all methods for both tissue and cell lysate (Table S7). In general, TMT showed the highest TPRs and AUROCs for all fold-changes (Table S7). For the tissue lysate and the smallest fold-change (FC2), TMT (TPR, 0.99) and LFQ (TPR, 0.99) showed a TPR higher than 0.95. In contrast, TMT (TPR, 0.97) and spike-in-SILAC (TPR, 0.95) showed the highest TPRs for cell lysate, whereas LFQ yielded TPRs below 0.95 for FC2 (TPR, 0.59) and FC3 (TPR, 0.92). AUROC values showed similar trends to TPR at the FPR level of 0.05 (Table S7). MBR or MBRreq generally resulted in lower TPRs and AUROCs, further highlighting that MBR/MBRreq traded the numbers of identified phosphosites (Figure E,F) and of significantly different phosphosites (Figure A,B), for accuracy, precision (Figure C,D), and robustness of the quantification (Figure C,D).

Taken together, TMT showed the lowest accuracy but the highest precision and the highest number of significantly differential phosphosites for all fold-changes in tumor and cell lysate. TMT also provided the highest TPR regardless of sample matrix. In contrast, the matrix affected the TPRs for LFQ and spike-in-SILAC, with spike-in-SILAC exhibiting a TPR below 0.95 for tissue and LFQ yielding a TPR below 0.95 for the cell lysate.

Coverage of Kinase Targets by the Quantified Phosphosites

As TMT differs from LFQ and spike-in-SILAC regarding the identified phosphosite profiles (Figure A,B), we asked whether different quantification methods introduce a bias in the coverage of quantified kinase substrate sites. For this purpose, an enrichment analysis[71] for kinase substrates assigned to major oncogenic signaling pathways was conducted using the kinase-substrates database from PhosphoSitePlus[49] (Figure S8A,B). TMT showed a higher coverage of most pathways, as compared to LFQ and spike-in-SILAC. This was mitigated by MBR and MBRreq, likely due to the higher number of quantified phosphosites through identification transfer from samples with higher phosphorylation quantities (Figure A,B). However, only 10% of the quantified phosphosites were annotated in the kinase substrate database, limiting the power of this analysis and raising the possibility that further biases may have been missed. We therefore analyzed the coverage of kinase target motifs among the significantly quantified phosphosite. Whereas LFQ and spike-in-SILAC without or with MBR/MBRreq covered a similar array of kinase target motifs, their coverage differed for TMT (Figure S9A,B). We conclude that the profiles of both identified and quantified phosphosites differ between LFQ, spike-in-SILAC, and TMT.

Discussion

We report the first comprehensive comparative analysis of the quantitative performance of LFQ, spike-in-SILAC, and TMT for the analysis of the tumor tissue phosphoproteome. In summary, LFQ yielded the highest number of phosphosite identifications. MBR and MBRreq increased the number of identified phosphosites for LFQ and spike-in-SILAC (Figure ). The lower number of phosphosites identified in TMT is likely explained by the SPS-MS3 approach where MS2 was used for phosphopeptide identification and phosphosite localization, while MS3 for quantification of the reporter ions requires a longer duty cycle and ion trap resonance CID. In ion-trap resonance CID used for TMT analysis, the depth of the potential prevents ion ejection and precursors can only be excited to a few electron volts, requiring long activation times to build up sufficiently high internal energy for fragmentation.[72] Gas-phase rearrangement reactions can occur prior to dissociation and have been reported for phosphate moieties of phosphoserine- and phosphothreonine-containing peptides.[73] Furthermore, only the precursor is excited in ion-resonance CID.[74] For phosphopeptides containing phosphoserines and phosphothreonines, this can result in abundant nonsequence informative fragment ions corresponding to the neutral loss of the phosphate moiety from the precursor.[75,76] In contrast, in beam-type HCD, which was used for fragmentation in LFQ and spike-in-SILAC, all ions are activated, and fragments including the phosphate neutral loss can undergo several consecutive fragmentation events.[77] Therefore, beam-type HCD spectra contain more sequence informative fragment ions than resonance CID spectra, and the shorter activation times reduced potential gas-phase rearrangements. As beam-type HCD spectra are better suited for accurate phosphosite identification than resonance CID,[78] TMT yields in total less identifications than LFQ or spike-in-SILAC.

Nevertheless, TMT provided the highest number of reproducible identifications across all fold-changes and phosphosite quantities, whereas for LFQ and spike-in-SILAC, the number of reproducibly identified phosphosites decreased with decreasing phosphosite quantities, which is mainly related to DDA precursor selection stochasticity. This suggests that in these approaches, reproducible identification is highly dependent on phosphosite quantities. The isobaric and multiplexing nature of TMT tags increases the signal intensity of phosphorylated peptides at the MS1 level, reducing the number of missing values. Although the number of reproducibly identified phosphosites could be increased with MBR and MBRreq for LFQ and spike-in-SILAC, TMT outperformed the other approaches in this respect for the analysis of both tissue and cell lysates. The TMT method requires a smaller initial sample amount as samples are pooled after trypsin digestion and TMT-labeling. This represents a considerable advantage for the study of clinical samples whose availability is often limiting. Regarding phosphosite quantification, LFQ without MBR provided the highest accuracy and the lowest precision, whereas TMT showed the lowest accuracy and the highest precision, irrespective of the matrix (Figure ). MBR and MBRreq sacrifice precision and accuracy for a higher number of quantified phosphosites, with spike-in-SILAC suffering more severely from this issue than LFQ. In both tissue and cell lysate, TMT exhibited the highest TPRs for the differentially quantified phosphosites (Figure , Table S7). In contrast, we found matrix effects for both LFQ and spike-in-SILAC. LFQ yielded high TPRs for the phosphoproteome in tissue but not cell lysate, whereas spike-in-SILAC showed high TPRs in cell lysate but not for tumor tissue. MBR decreased the TPRs and could not mitigate these matrix-effects. However, MBR increased the number of significantly different quantified phosphosites in LFQ and spike-in-SILAC, albeit not to the same amount as reached by TMT. This further highlights that MBR trades the number of phosphosites identified and quantified for quantification accuracy, precision, and robustness.

We analyzed the tumor phosphoproteome side-by-side with a matched cancer cell line to ensure comparability of our data with previous studies on cell lines and assess the robustness of the quantitative methods toward matrices as different as tumor and cell lysate. Our analysis of the cell lysate reflects overall the findings of Hogrebe et al.,[26] although there are some differences in study design that need to be considered. Hogrebe et al.[26] used the protease Lys-C, whereas we used trypsin as the most widely used protease for bottom-up proteomics.[79,80] TMT labeling after titanium dioxide (TiO2) enrichment, as done by Hogrebe et al.,[26] may reduce the advantage of TMT as multiplexing occurs late in the protocol and does not compensate technical variability introduced in the enrichment step. Therefore, we performed TMT labeling and multiplexing directly after trypsin digestion and used IMAC instead of TiO2 due to its higher selectivity, identification numbers, and quantitative reproducibility.[81]

The quantification method significantly determines the study result in terms of identified, quantified, and differentially assigned phosphosites as well as quantification accuracy and last but not least the tumor related kinase phosphosite target profiles covered. For study design, the choice of the quantification method depends on the study aims. Although TMT outperformed spike-in-SILAC regarding quantification precision and was least amenable to matrix effects, it exhibited the lowest accuracy. If accuracy is a key requirement, TMT is therefore not the method of choice. Spike-in-SILAC provided the best compromise between accuracy, precision, and robustness toward differences in phosphosite quantities. However, spike-in-SILAC suffered from a low phosphosite coverage and did not reliably detect low fold changes between phosphosite quantities in tissue. Its performance should therefore be assessed for the matrix to be analyzed (i.e., the tumor tissue) in order to determine the suitability of spike-in-SILAC regarding the TPR of low fold changes. LFQ was susceptible to matrix effects but exhibited a higher TPR in tumor tissue than in cell lysate. However, LFQ offered low precision, which was not counterbalanced by MBR. Rather, MBR degraded the precision and robustness of quantification. Recently, Bekker-Jensen et al. showed that data-independent acquisition (DIA) LFQ may overcome some of the limitations in DDA LFQ, which we observed in our study.[82] Also, practical considerations such as instrument time may influence the choice of method. TMT multiplexing enabled the concomitant analysis of 9 samples in one measurement, reducing the analysis time by a factor of 9 compared to LFQ and spike-in-SILAC. This represents an inherent advantage of TMT regarding time and costs. In the meantime, the TMTpro reagent set became available, allowing multiplexed analysis of up to 18 samples,[83] further enhancing the efficient use of instrument time. To compensate for the longer duty cycle of the SPS-MS3 approach, several studies[84,85] found that real-time search improved the number of identified and quantified peptides and proteins and resulted in higher reproducibility and accuracy of quantification and a larger dynamic quantification range of bottom-up proteomics analysis. Similarly, FAIMS and hrMS2/SPS-MS3 methods[86,87] have been reported to reduce sample complexity and coisolation of interfering small peptides in MS1 DDA precursor selection, which improved the number of identified peptides and proteins and the accuracy and precision of their quantification. These methods are therefore attractive to further improve the performance of SPS-MS3 TMT, for instance for applications in large-scale clinical studies.

In conclusion, different quantification methods offer the highest accuracy, precision, and phosphosite coverage in tumor tissue proteomics and thus the choice of quantification method is critical. The different behavior of TMT, LFQ, and spike-in-SILAC as well as the influence of MBR should also be taken into account when comparing published tumor tissue phosphoproteomes. We advocate careful annotation of tissue phosphoproteomes to enable meaningful comparison.