Label-free quantitative phosphoproteomics with novel pairwise abundance normalization reveals synergistic RAS and CIP2A signaling

Abstract

Hyperactivated RAS drives progression of many human malignancies. However, oncogenic activity of RAS is dependent on simultaneous inactivation of protein phosphatase 2A (PP2A) activity. Although PP2A is known to regulate some of the RAS effector pathways, it has not been systematically assessed how these proteins functionally interact. Here we have analyzed phosphoproteomes regulated by either RAS or PP2A, by phosphopeptide enrichment followed by mass-spectrometry-based label-free quantification. To allow data normalization in situations where depletion of RAS or PP2A inhibitor CIP2A causes a large uni-directional change in the phosphopeptide abundance, we developed a novel normalization strategy, named pairwise normalization. This normalization is based on adjusting phosphopeptide abundances measured before and after the enrichment. The superior performance of the pairwise normalization was verified by various independent methods. Additionally, we demonstrate how the selected normalization method influences the downstream analyses and interpretation of pathway activities. Consequently, bioinformatics analysis of RAS and CIP2A regulated phosphoproteomes revealed a significant overlap in their functional pathways. This is most likely biologically meaningful as we observed a synergistic survival effect between CIP2A and RAS expression as well as KRAS activating mutations in TCGA pan-cancer data set, and synergistic relationship between CIP2A and KRAS depletion in colony growth assays.

Figure 1. A schematic effect of a normalization bias caused by manipulation of RAS and PP2A phosphoproteomes.

(a) Protein phosphatase 2A (PP2A) participates in the regulation of a large part of phosphoproteome, including major serine/threonine kinases AKT and ERK that are also key downstream effectors of the RAS oncoproteins. RNAi mediated depletion of RAS, PP2A activation by depletion of CIP2A protein, and PP2A inhibition by OA were used as model perturbations, to study the influence of global phosphorylation changes on the performance of different normalization methods in label-free quantitative phosphoproteomics. (b) Centering normalization is often used in quantitative proteomics and phosphoproteomics data (upper panel). However, a global phosphorylation change shifts the distribution of the phosphorylation ratios (middle panel). In such cases, centering leads to normalization bias, which introduces false positive phosphorylations in the opposite direction from the global change and also false negatives in the direction of the global change (lower panel).

Figure 2. Pairwise normalization developed for label-free quantitative phosphoproteomics.

(a) HeLa cells with different treatments were subjected to cell lysis, spiking α-casein standard, and tryptic digestion. Peptides with and without TiO₂ phosphopeptide enrichment were analyzed by LC-MS/MS. Peptides were identified by Mascot database search, followed by phosphorylation site validation by phosphoRS. Phosphopeptide identification was supplemented by SimSpectraST spectral library search. Following label-free quantification, peptide abundance was normalized with different methods, including the pairwise normalization for TiO₂ data developed in this study. (b) The principle of the pairwise normalization method. Fifty-two phosphopeptides were quantified in both the non-enriched digests and TiO₂-enriched samples (i.e. 52 digest-TiO₂ pairs). Abundance profiles of two hypothetical phosphopeptides are illustrated as examples. An abundance ratio was calculated by pairwise comparison (digest/TiO₂) for each phosphopeptide. Eleven pairs were excluded as outliers (see the criteria in Supplementary Fig. 3). The median of normalized abundance ratios was then calculated for the remaining 41 pairs and used as a pairwise normalization factor for the TiO₂ data. The TiO₂ data were pre-normalized with the global centering method, whereas the digest data were normalized with the global centering or quantile centering method (i.e. global pairwise and quantile pairwise, respectively).

Figure 3. Fold change distributions of phosphorylations after different normalizations.

(a) Fold changes for each phosphopeptide ion feature was calculated for the CIP2A, RAS, or OA samples compared to the control 1 samples (log-transformed). The abundance of the features was normalized with global centering and quantile centering methods. Median and mean levels are marked with a solid and dashed line on the box plots, respectively, and whiskers represent 1.5 × interquartile range. (b) Ratio of up- and down-regulated phosphosites (differentially regulated phosphosites compared to the control 1 samples; t-test, p < 0.01) is shown for both normalization methods and non-normalized data. Abundances of the features with identical protein phosphorylations were summed up for calculating phosphosite abundance. The centering normalizations resulted in similar ratios of up- and downregulated phophosites in contrast to the expected phosphoproteome changes (i.e. increase in protein phosphorylation after OA treatment and dephosphorylation after CIP2A or RAS depletion, refer to Fig. 1a). (c) Fold changes of phosphopeptide features and (d) ratio of up- and down-regulated phosphosites (t-test, p < 0.01) after pairwise normalizations. Global pairwise normalization of the data resulted in the best agreement with the expected global phosphoproteome changes (see Fig. 1a).

Figure 4. Hierarchical clustering of the samples after different normalizations.

(a) The log-transformed, normalized phosphosite data was clustered using a variety of distance metrics and clustering strategies. Euclidean distance-based Ward’s minimum variance clustering for the global-pairwise-normalized data is shown here as an example. CIP2A and RAS formed a tight cluster that was clearly separated from OA, and also distinguished from the control sample cluster. (b) Various cuts on the clustering distance height were applied (horizontal lines 1, 2 or 3 in panel a) to produce subclusters of different sizes. Here, clustering solutions with 2, 3 or 6 clusters are shown. (c) The sample clusters at various height cuts were compared to the original sample groups using the adjusted Rand index computed for each of the 5 normalization methods, and AUC was used to compare between the methods. The AUC values for different clustering parameter combinations are shown in Supplementary Table 7. (d) PCA plots for the quantile-centering-normalized and quantile-pairwise-normalized data. Variance among the OA samples led to sub-optimal grouping in the quantile centering normalization (left panel).

Figure 5. Western blot validation of phosphorylations.

(a) Western blotting was performed on the cell lysates used for LC-MS/MS analysis. Representative western blots for each antibody are shown. (See Supplementary Fig. 5 for different exposure times). (b) Quantitative results of the phosphorylation regulations obtained by western blotting were compared with LC-MS/MS results with different normalizations. Fold-changes (average of triplicates) compared to the control 1 samples are shown. The directions of phosphosite regulations (i.e. up or down) in the CIP2A, RAS, and OA samples (individual replicates) were also compared to the average of control 1 samples. The agreement with western blot was compared between different normalizations using Fisher’s exact test. (c) Average correlation coefficients for phosphosites were calculated between the western blotting and LC-MS/MS results on log-transformed data. As the OA samples significantly skewed the data dominating the Pearson’s correlation coefficients, they were excluded from the calculations. Global pairwise normalization led to the highest correlation. (d) Spearman’s ρ and Kendall’s τ rank correlation coefficients were also calculated for phosphosites in all samples (i.e. the OA samples included). WB: western blotting, GP: global pairwise, QC: quantile centering, QP: quantile pairwise, GC: global centering, NN: non-normalized, and Ca: casein.

Figure 6. Pathway analysis for protein and phosphorylation regulations.

(a) The protein and phosphosite fold changes (compared to control 1) were calculated from global-centering-normalized non-enriched data and global-pairwise-normalized TiO₂ data, respectively. In Ingenuity Pathway Analysis, core analysis was performed for differentially regulated proteins and phosphosites (t-test, p < 0.05), followed by comparison analysis between the CIP2A, RAS, and OA core analyses. The top hits from the category “Diseases and Bio functions” are shown. (b) The phosphosite data was filtered for those regulated by both CIP2A and RAS depletions (t-test, p < 0.05), and the core analysis was performed. Upstream regulator analysis restricted to kinases is shown. (c) AKT and ERK target sites were predicted by NetworKIN and GPS tools or retrieved from the PhosphositePlus database (see Supplementary Fig. 8). The average fold changes for AKT and ERK target sites are presented for the global-centering (left) and global-pairwise (right) normalized data. The expected regulations of AKT and ERK mediated phosphorylations were clearly observed by global pairwise normalization. The error bars represent standard error of the mean (SEM). The asterisks represent level of statistical significance for up-/down-regulations (one sample t-test, *p < 0.05, **p < 0.01, ***p < 0.001).

Figure 7. CIP2A and KRAS show synergy in TCGA data and in colony formation assay.

(a) The Kaplan-Meier plot of TCGA pan-cancer survival profiles. The patients were split into groups by pan-cancer normalized gene-expression data for KRAS, NRAS and CIP2A (downregulated < 0, upregulated ≥ 0), and by KRAS mutational status. Combination of low CIP2A expression with non-mutated or low expression of RAS was associated with best survival. Log-rank test was used for comparing the survival distributions between the group with the best survival and the other groups. (b) Colony formation assay following CIP2A and KRAS depletions as well as CIP2A + KRAS co-depletion. Assay was performed 3 times using different siRNAs for CIP2A in HeLa, CW-2, HCA7, and NCI H747 cell lines. The average result is shown. See also Supplementary Fig. 10 for the triple-RAS depletion. (c) Averages of 3 colony formation experiments for each cell line. The error bars represent SEM. The asterisks represent level of statistical significance (t-test, *p = 0.0162, **p = 0.0012).