Tissue-Specific Signaling Networks Rewired by Major Somatic Mutations in Human Cancer Revealed by Proteome-Wide Discovery

More than 46 million somatic mutations have been accumulated through several large-scale cancer genome sequencing projects such as The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium projects (1, 2). Currently, how to assess the impact of these somatic mutations in the process of tumorigenesis and disease progression is a main challenge. Considering the existing observations that many somatic mutations promote tumorigenesis by rewiring protein signaling networks (3), one potential strategy is to incorporate somatic mutations with the protein structure information and investigate them at functional sites (e.g., phosphorylation sites; refs. 4–8). Phosphorylation-dependent signaling network is fundamental in cellular physiology and its dysfunction plays a critical role in tumorigenesis. Recently, the rapid advancements of technologies (e.g., mass spectrometry) have made a wealth of kinase–substrate interactions with the specific phosphorylation sites available (9). This provides us with a unique opportunity to interrogate the impact of somatic mutations on kinase–substrate phosphorylation sites to determine their pathophysiologic roles in cancer and prioritize potentially actionable mutations (10).

Previous studies have suggested that many infrequent but specific somatic mutations may affect cancer hallmark processes and lead to tumorigenesis (11). However, a systematic investigation of the somatic mutations in phosphorylation networks and pathways governing cell signaling has not been done yet. Meanwhile, the human kinome has become one of the most important classes of drug targets in the pharmaceutical industry (12, 13). So far, more than 20 drugs targeting one or more kinases have been approved for clinical use in a variety of cancers, including lung, breast, melanoma, colorectal, pancreatic, and prostate cancers (14). Moreover, as of 2012, more than 500 kinase inhibitors have been investigated as therapeutic agents, approximately a third of which are undergoing clinical trials (15, 16). Frustratingly, while effective, patients treated with these kinase inhibitors eventually develop resistance, and the prolong survivals are typically only a few months (17, 18). One primary reason for resistance is that kinases are extensively involved in complex biological mechanisms through adaptive crosstalk or feedback within cellular networks (19). This suggests that it is imperative not only to investigate somatic mutations on individual kinase–substrate protein relationships but also to acquire in-depth knowledge of the interconnectivity of these networks and the biochemical specificity of phospho-catalytic events (20, 21). So far, a systematic investigation of signaling networks disturbed by somatic mutations at the kinome-wide has not been done yet. There is an urgent need to develop novel computational approaches for kinome-wide identification of their rewiring signaling networks altered by somatic mutations at functional sites (e.g., phosphorylation sites) so that the findings, especially the top candidates, can expedite cancer pharmacogenomics studies in the emerging era of precision oncology (10, 22).

In this study, we incorporated the somatic missense mutations into the protein phosphorylation sites to decipher the functional role of somatic mutations. We found a significant higher mutation rate at phosphorylation sites in comparison with the whole protein sequences in all the 16 cancer types that we examined. Our analysis revealed that mutations at protein phosphorylation sites were more likely deleterious than those at nonphosphorylation sites. We examined over 740,000 missense mutations in nearly 5,000 tumor-normal matching samples across 16 cancer types from TCGA. We found 47 statistically significantly mutated proteins (SMP) that have enriched missense mutations at their phosphorylation sites in this pan-cancer analysis. We further characterized the distribution of phosphorylation site mutations in kinase-substrate network and developed a new algorithm to detect significantly tissue-specific kinase-substrate interaction modules altered by somatic mutations. Functional annotation analysis of these modules revealed the tissue-specific phosphorylation-related mechanisms of somatic mutations in cancer. We found that the significantly mutated phosphorylation sites or proteins and their perturbed network modules were highly associated with patient survival and anticancer drug responses. Finally, we reported a kinome-wide landscape of pharmacogenomic interactions through integrating somatic mutation-rewired signaling networks in approximately 5,000 primary tumors and approximately 1,000 cancer cell lines using a statistical framework.

The detailed descriptions for Experimental Procedures are provided in the Supplementary Data.

We downloaded the tumor-normal pairwise somatic mutation data for patients from three sources: (i) Elledge Lab website at Harvard University (http://elledgelab.med.harvard.edu/?page_id=689; accessed in March 2014); (ii) Sanger website: (ftp://ftp.sanger.ac.uk/pub/cancer/AlexandrovEtAl, December 16, 2013); and (iii) Catalogue of Somatic Mutations in Cancer (v69). To reduce redundancy and ensure the quality of somatic mutation data, in this study, we only focused on the somatic mutations in TCGA tumor-normal matched samples from aforementioned three datasets (see Supplementary Materials and Methods).

We collected over 170,000 carefully curated nonredundant phosphorylation sites covering 18,610 proteins from the PhosphoSitePlus (23) and dbPTM3 databases (24).

We compiled high-quality kinase–substrate interactions (KSI) from two sources: (i) literature-derived low-throughput experiments, and (ii) high-throughput experiments (e.g., protein microarrays and mass spectrometry-based approaches; ref. 25). As a result, we compiled 7,346 unique KSI pairs connecting 379 kinases and 1,961 nonkinase substrate proteins (see Supplementary Methods and Materials).

We downloaded putative somatic mutations for 1,001 cancer cell lines and natural log IC₅₀ and the area under the dose–response curve values for all screened cell line/drug combinations from the Genomics of Drug Sensitivity in Cancer (http://www.cancerrxgene.org/).

All flanking regions with 7-residue window of phosphorylation sites as described previously (5) were tested for the recurrence of mutations. We assumed that the observed number of mutations for a given region followed a binomial distribution, binomial (⁠|T$|⁠,|{\rm{\ }}{p_{{g_i}}}$|⁠), where |T$| was the total number of mutations observed in one gene and |{p_{{g_i}}}$| was the estimated mutation rate for the region of interest in gene |{g_i}$| under the null hypothesis that the region was not recurrently mutated. Using |length\ ( {{g_i}} )$| to represent the length of protein product of gene |{g_i}$|⁠, for each phosphorylation site, we computed the P value, the probability of observing more than |k$| mutations around this phosphorylation site out of |T$| total mutations observed in this gene, using the equation below:

in which |{p_{{g_i}}} = {{\rm{\ }}\frac{{{\rm{window}}\_{\rm{size}}}}{{{\rm{length}}( {{g_i}} )}}$|⁠. Finally, we set the minimal P value across all the phosphorylation sites in a specific protein as the representative P value of its coding gene |{g_i}$|⁠, denoted |P( {{g_i}} )$|⁠.

We aim to integrate information from both kinase–substrate interaction network (KSIN) and coexpression networks to identify cancer-specific modules that are enriched for phosphorylation site mutations and highly coexpressed in the according tissue. First, for each cancer type, we used the expression data of matched normal tissues to calculate the coexpression value for each kinase–substrate interaction and used cut-off P < 0.05 to construct coexpressed KSIN. For each gene, we assigned a score |S( {{g_i}} ) = - {\rm{log}}( {P( {{g_i}} )} )$| based on its respective P value. Then, given one module |M = \{ {{g_1},{g_2}, \cdots, {g_n}} \}$|⁠, we could define the objective scoring function: |{\rm{Obj}}( M ) = {\frac{{\mathop \sum \nolimits_{g \epsiv M} ( {S( g ) - \mu } )}}{{\sqrt n }}$|⁠, in which μ is the mean value of the scores for all the genes considered.

We proposed four following steps to detect densely mutated kinase-substrate modules in coexpressed KSIN:

• Start from one of the substrates with significantly mutated phosphorylation sites and its upstream kinase and downstream substrate as the seed.

• Select one gene from the neighbors of the current module, and use it to expand the current module. Continue this expansion until the objective score stops growing.

• Repeat step (ii) 10,000 times and 10,000 modules are generated accordingly.

• Construct a consensus module for each cancer type by merging the top 5% modules with the highest objective scores.

For this ensemble of the genes in the consensus module, we calculated how many times each gene appears in different suboptimal solutions and finally assigned a “confidence score” for each gene.

For each drug, we constructed a drug-response vector consisting of IC₅₀ values from treatment of cell lines. Then, drug-response vector was modeled as a linear combination of the tissue of origin of the cell lines, screening medium, growth properties, and the status of a genomic feature. In this study, considering the data sparsity, we only performed pan-cancer analysis. A genomic feature–drug pair was tested only if the final drug-response vector contained at least three positive cell lines and at least three negative cell lines. Effect size was quantified through the Cohen's |d$| using the difference between two means divided by a pooled SD for the data. The resulting P values were corrected by Benjamini–Hochberg method (26).

All the statistical analyses were performed by the R package (v3.2.3, http://www.r-project.org/).

In this study, we developed a novel computational oncoproteomics approach, namely kinome-wide network module for cancer pharmacogenomics (KNMPx), to identify novel actionable mutations that have rewired signaling networks and to further characterize tumorigenesis and anticancer drug responses. As shown in Fig. 1, KNMPx includes three components: (i) identify SMPs that harbor the enriched missense mutations at their phosphorylation sites (Fig. 1A); (ii) detect significantly tissue-specific rewired signaling network modules through integrating somatic mutation–altered phosphorylation sites into tissue-specific KSIN using a novel network-based statistical approach (Fig. 1B); (ii) build a kinome-wide landscape of pharmacogenomic interactions by integrating somatic mutation-rewired signaling networks in 5,000 primary tumors and 1,001 cancer cell lines using a statistical framework (Fig. 1C). In the current KNMPx approach, we investigated a total of 746,631 somatic mutations in 4,997 tumor samples across 16 major cancer types/subtypes from TCGA. In this study, we only used missense mutations—the single nucleotide base substitutions leading to amino acid changes. Among the 16 cancer types, acute myeloid leukemia (LAML) had the lowest mean mutation rate (8 missense mutations per sample), whereas uterine corpus endometrial carcinoma (UCEC) had the highest mutation rate (497 per sample). Overall, the average mutation rate across the 16 cancer types is 160.

We mapped the 746,631 missense mutations into 15,401 human proteins harboring phosphorylation sites. Then, we obtained 8,057 unique missense mutations directly at the phosphorylation sites across 4,331 unique human proteins. After considering the somatic mutations at the seven immediate flanking residues of each phosphorylation site similarly applied in a previous study, we obtained 97,277 phosphorylation site mutations affecting 12,552 proteins (including 5,547 protein having recurrently mutated phosphorylation sites—at least 3 mutations affecting the same site).

We first examined the 4,331 proteins bearing mutations directly at the phosphorylation sites. We compared mutation frequency in the whole protein sequences with that at the phosphorylation sites across 16 cancer types. We found a significantly higher mutation rate at the direct phosphorylation sites than that of the whole protein sequence in nine cancer types among the total 16 cancer types (Fig. 2A). For the instance of breast invasive carcinoma, the average mutation rate is 13.8 per million amino acids at phosphorylation sites. This compared with 5.5 mutations per million amino acids at the whole protein level. The mutation rate at the phosphorylation sites is nearly three times of that at the whole protein level (P = 1.0 × 10⁻⁶¹, Wilcoxon rank-sum test). On the sample level, in pan-cancer analysis, more than 50% of samples harbored somatic mutation at phosphorylation site in at least one protein (Supplementary Table S1). Especially in lung squamous cell carcinoma (LUSC) and skin cutaneous melanoma (SKCM), the percentage of samples having phosphorylation site mutations is up to 80%. These observations indicated that phosphorylation sites were significantly altered by somatic missense mutations.

We tested the hypothesis that somatic missense mutations at phosphorylation sites have more deleterious impact than missense mutations at other sites. We measured the functional impact scores of mutations using two popular and complementary tools: SIFT (27) and PolyPhen-2 (28). We examined the cumulative distribution of SIFT and PolyPhen-2 scores for phosphorylation site mutations (direct position and seven flanking residues) and nonbinding site mutations (amino acid positions excluding the binding sites and seven flanking residues). The distribution is shown in Fig. 2B and C. We defined the mutations with SIFT scores <0.95 or PolyPhen-2 score >0.909 as deleterious mutations based on previous studies (4). We found that phosphorylation site mutations were more likely to be deleterious than nonphosphorylation site mutations when they were evaluated using both SIFT (P = 1.47 × 10⁻⁶, Fisher exact test; Fig. 2B) and PolyPhen-2 (P = 0.0038; Fig. 2C) scores.

By analyzing 746,631 missense mutations from 4,997 tumors across 16 cancer types from TCGA using KNMPx, we found that more than 50% tumor samples harbored phosphorylation-associated single nucleotide variants (SNV), and further identified 47 proteins that harbored significantly mutated phosphorylation sites (false discovery rate q < 0.3) including several well-known cancer proteins: BRAF: p.T599 (q = 2.2 × 10⁻³⁰⁸), TP53: p.S269 (q = 2.2 × 10⁻⁴⁵), and EGFR: p.T290 (q = 3.2 × 10⁻⁷; Fig. 3A). We next evaluated the predictions using two benchmark cancer gene sets. Specifically, we examined the enrichment of these 47 SMPs using two well-curated cancer gene datasets: (i) 547 experimentally validated cancer genes from Cancer Gene Census (accessed on November 5, 2014, denoted as CGC genes; ref. 29), and (ii) a curated list of 693 significantly mutated genes (SMG) in cancer from TCGA publications and several other publications (details were provided in Supplementary Table S2 and described in our previous study; ref. 7). We found that the genes harboring significantly mutated phosphorylation sites were significantly enriched in CGC (P = 1.1 × 10⁻¹⁹, observed overlap = 29, expected = 3, hypergeometric test) and SMGs (P = 1.3 × 10⁻²⁵, observed overlap = 37, expected = 3), respectively.

Next, we prioritized potential SMPs for each individual cancer type (Fig. 3B; Supplementary Fig. S1; Supplementary Table S3). In total, we identified 74 SMPs with q < 0.3 by integrating results in 16 cancer types, including 29 CGC proteins (P = 4.7 × 10⁻²⁶, hypergeometric test) and 37 published SMG products (P = 1.5 × 10⁻³⁷). Using the benchmark classification from Vogelstein (30), these 74 SMPs included nine tumor suppressor proteins (TP53, APC, PIK3R1, ATM, SMAD4, SOX9, VHL, CREBBP, and ARID1A) and five oncogene proteins (BRAF, CTNNB1, EGFR, KIT, and ERBB2). In the case of BRCA, we found 2 SMPs: ERBB2: Y772 (q = 0.00043), TP53: S269 (q = 0.0071). Both are well-known breast cancer–related proteins. For ERBB2, which encodes HER2, four mutations (V777L, D769Y, D769H, and I767M) may alter its phosphorylation site of p.Y772. Three of them (V777L, D769Y, and D769H) were functionally characterized as activating mutations in HER2 gene amplification-negative breast cancer (31). Interestingly, R273H on TP53 altering its phosphorylation site p.S269 was recently reported as an oncogenic mutation in the mouse model (32). In the analysis of UCEC, we found nine SMPs with q < 0.3 including CTNNB1: p.S37 (q = 1.7 × 10⁻⁹³), CCND1: p.T286 (q = 0.00049), PIK3R1: p.T576 (q = 0.0024), and TP53: p.S269 (q = 0.0047). CTNNB1 (Catenin beta-1) is the key downstream component of the canonical Wnt signaling pathway, which plays crucial roles in the regulation of diverse cell behaviors, including cell fate, proliferation, survival, differentiation, migration, and polarity (33). In BLCA, there were 4 SMPs: AHNAK: p.S1943 (q = 0.0012), TP53: p.T155 (q = 0.012), ILF3: p.S482 (q = 0.044), and ARID1A: p.S607 (q = 0.1). AHNAK is originally identified as a nuclear phosphoprotein in human neuroblastomas and skin epithelial cells. A recent study reported that AHNAK could mediate negative regulation of cell growth and act as novel tumor suppressor through potentiation of TGFβ signaling (34). ILF3 promotes breast tumor progression by regulating sustained urokinase-type plasminogen activator expression (35). Altogether, we demonstrated that KNMPx is a useful approach for identifying cancer-associated proteins harboring enriched somatic mutations at their phosphorylation sites.

We next investigated tissue-specific phosphorylation site mutations in kinase–substrate network. For each cancer type, the Pearson correlation coefficient (PCC) of each kinase–substrate pair in KSIN was calculated based on the normalized gene expression data from RNASeqV2 of the matched normal tissue in TCGA. Here, we defined a putative somatic mutation-altered kinase–substrate interaction if the corresponding substrate harbored any phosphorylation site mutations. Accordingly, we found that in nine cancer types the kinase–substrate interactions affected by phosphorylation site mutations were significantly shifted toward higher correlation coefficients when compared with the distributions expected for the randomly selected kinase–substrate interactions that are not affected by phosphorylation site mutations (P < 0.05 in 9 of 12 cancer types that had matched normal tissues, Kolmogorov–Smirnov test). The most significant patterns were observed in LUSC (P = 5.9 × 10⁻⁶), UCEC (P = 0.0026), and lung adenocarcinoma (LUAD; P = 0.0028), which is consistent with the reports that phosphorylation site mutations have been frequently observed in these three cancer types. Altogether, the clustering of phosphorylation site mutations in highly correlated kinase-substrate interactions that we observed might suggest tissue-specific patterns of phosphorylation site mutations in KSIN. Standing on this feature, we next integrated tissue-specific coexpression information to search for the highly mutated modules in kinase–substrate network.

We hypothesized that phosphorylation site mutations rewired kinase–substrate interaction networks and further associated with tumorigenesis or anticancer drug responses (36). As shown in KNMPx framework (Fig. 1B), we searched the consensus mutant modules via 4 steps as described in Materials and Methods. We identified one consensus module for each cancer type, resulting in a total of 16 consensus modules (Supplementary Table S4). We then combined 16 consensus mutant modules, which included 330 genes for follow up pan-cancer analysis. We found that these 330 genes were significantly enriched in CGC (P = 8.4 × 10⁻²⁴, hypergeometric test) and the previously reported SMGs (P = 3.1 × 10⁻³¹, hypergeometric test). In addition, genes in the consensus mutant modules for all the 16 individual cancer types were also significantly enriched in CGC (P values ranging from 2.8 × 10⁻¹¹ to 0.017, hypergeometric test) and previously reported SMGs (P values ranging from 1.3 × 10⁻²¹ to 0.0035, hypergeometric test). These enrichment analyses revealed the functional relevance of consensus mutant modules identified by our KNMPx approach. We further examined the correlation of expression of genes in the consensus modules with the patient survival information using LUAD and ovarian serous cystadenocarcinoma (OV) as two examples. We found that the expression of genes in consensus modules was highly correlated with the patient survival (Supplementary Fig. S2). For example, high expression of several genes in consensus modules was significantly associated with poor survival, such as CD44 (P = 0.01) and LCK (P = 0.044) in OV, and BAD (P = 0.012) and RGPD4 (P = 0.048) in LUAD.

In UCEC, the final consensus module (mutated in 90/248 samples) consisted of 11 genes including two well-known cancer genes (CTNNB1 and TP53) but also the genes with mutations that were too rare to be significant by single-gene tests (Fig. 4A). Most of these genes have literature evidence supporting their roles in cancer. For instance, we found Bruton tyrosine kinase (BTK), a key component of B-cell receptor (BCR) signaling and functions as an important regulator of cell proliferation and cell survival in various B-cell malignancies, has three somatic mutations around its phosphorylation site p.S21. A previous study has reported that Btk phosphorylation at p.S21 creates a binding site for the prolyl isomerase Pin1, which modulates Btk activity in a cell-cycle–dependent manner (37). We observed significant mutual exclusivity across the 11 genes in the module, which was largely due to the previously reported mutual exclusivity between TP53 and CTNNB1 mutations in UCEC (Fig. 4B; ref. 38). For example, 66 missense mutations on CTNNB1 in UCEC around site p.S37 (Fig. 4C and D) might alter β-catenin signaling in colon cancer (39). A previous study has suggested that the mutations within a pathway tended to exhibit mutual exclusivity (40). Enrichment analysis also showed that seven of them (GSK3B, FYN, RASSF1, TP53, BTK ABL1, and MAP2) belong to developmental process annotated in Gene Ontology (P = 7.21 × 10⁻³, http://geneontology.org/). These observations suggested that the phosphorylation site mutations in these 11 genes might define a subtype of UCEC.

Although it is well known that phosphorylation is frequently involved in the binding of drug inhibitors to target proteins, there has not been a systematic survey of the mutations at the phosphorylation sites against drug sensitivity. Here, we analyzed putative somatic mutations across 1,001 human cancer cell lines and how their phosphorylation site mutations could correlate with drug sensitivity or resistance to the 265 drugs examined.

First, we examined the phosphorylation site mutation profiles from cancer cell lines and for each cancer cell line, we evaluated the extent of recapitulating the mutational profile observed in the corresponding primary tumor genomes. To this end, we applied the same method in analyzing primary tumor data to map the somatic mutations in cell lines to protein structures. Because kidney chromophobe (KICH) did not have the corresponding cell line data, we analyzed 15 cancer types only. Among the 277,073 putative somatic mutations revealed in these cancer cell lines, 29,821 were at the phosphorylation sites (locating at the direct phosphorylation sites or their 7-residue flanking regions) and half of them (15,626/29,821) also occurred in primary tumor. When comparing the 47 SMPs that we identified in primary tumor pan-cancer analysis, the correlation between the frequency of phosphorylation site mutations in cell lines and primary tumors was greater than 0.5 (from 0.53 in LAML to 0.97 in skin cutaneous melanoma [SKCM]) for all the 13 cancer types that had at least 10 samples with both variant and drug response data in primary tumors and the corresponding cell lines with the similar origins (Fig. 5). For example, thyroid carcinoma (THCA) primary tumors shows correlation of 0.70 with THCA cell lines and 0.95 with SKCM cell lines. Put together, this result confirmed that a large panel of cell lines could capture the phosphorylation site mutation patterns of primary tumors with the similar origin.

Next, we used ANOVA to identify phosphorylation site mutations as markers associated with drug responses. As shown in Fig. 6A, we identified 70 statistically significant pairwise interactions between phosphorylation site mutations and drug responses by combining data in 1,001 cell lines [false discovery rate (FDR) < 30% and Cohen |d$| > 0.8]. As shown in Fig. 6A and B, high sensitivity (low IC₅₀ value) for four BRAF inhibitors: dabrafenib (P = 1.7 × 10⁻⁵⁷), PLX4720 (P = 1.0 × 10⁻²⁶), SB590885 (P = 5.2 × 10⁻²¹), and AZ628 (P = 9.4 × 10⁻⁸) was significantly associated with BRAF p.T599 mutations. In addition, five MEK inhibitors, bicalutamide (P = 1.7 × 10⁻⁵⁷), RDEA119 (P = 1.2 × 10⁻¹⁴), PD-0325901 (P = 1.2 × 10⁻¹²), trametinib (P = 8.8 × 10⁻⁹), and CI-1040 (P = 1.2 × 10⁻⁷), had a higher response to NRAS p.Y64–mutated cell lines in comparison with wild-type cell lines. Interestingly, one hotspot mutation site, TP53 p.S269 also lead to a higher response to bicalutamide (P = 2.2 × 10⁻¹⁵) in TP53 p.S269–mutated cell lines compared with wild-type ones. Studies of multiple-target receptor tyrosine kinase inhibitors have suggested high efficiency in targeted cancer therapies (41). For example, linifanib, which has high binding affinity with PDGF receptor β, KDR, and colony-stimulating factor 1 receptor, is effective orally in VEGF-induced uterine edema and corneal angiogenesis (42). Figure 6A and B indicated that a higher response of linifanib was significantly associated with several phosphorylation site mutations on energy homeostasis enzyme, creatine kinase, M-type (CKM; P = 7.4 × 10⁻⁸) and a cell metabolism gene, SLC25A1 (P = 3.1 × 10⁻⁷). Furthermore, quizartinib (AC220), a FDA-approved FLT3 inhibitor for the treatment of acute myeloid leukemia, was found significantly associated with a higher response in cancer cell lines harboring the rewired phosphorylation network by several proteins, such as MKNK1 p.S39 (P = 1.7 × 10⁻⁹), PKP2 p.Y161 (P = 1.6 × 10⁻⁸), CKM p.T313 (P = 2.5 × 10⁻⁷), and USP37 p.S170 (P = 2.2 × 10⁻⁷), in comparison with wild-type groups. This pattern of multiple targets is similar to the multiple-target RTK inhibitors like linifanib. Put together, we identified several promising actionable mutations with their mechanisms on altering phosphorylation sites. These markers may be useful for further characterizing drug sensitivity to MAPK signaling and ERK-signaling inhibitors (Fig. 7).

We next explored potential mechanisms involved in the resistance of PI3K signaling inhibitors. For example, a PDK1 inhibitor, OSU-03012, showed potential for treating various cancers, such as multiple myeloma (43), brain tumor (44), and thyroid cancer. Figure 6B revealed that somatic mutations altering ATP8A1 p.S1126 might correlate with the high risk of resistance of OSU-03012 (P = 7.6 × 10⁻⁸). In addition, QL-VIII-58, an mTOR inhibitor, presented high risk of resistance with somatic mutations altering ABLIM3 p.T276 (P = 1.1 × 10⁻⁷). FH535, a Tcf/beta-catenin signaling inhibitor, showed high potential for the treatment of several cancer types, such as lung, liver, and colon cancers. Figure 6B displayed that mutations around APBB1IP p.S652 were significantly correlated with the resistance of FH535 in pan-cancer analysis. A previous study has suggested that APBB1IP plays an important role in Rap1 signaling pathway that is closely involved with beta-catenin signaling (45). Collectively, rewired signaling network altered by phosphorylation site mutations is an important mechanism in mediating the resistance of cancer therapies for various signaling pathways, such as PI3K signaling and Wnt/β-catenin signaling (Fig. 7).

In this study, we developed an oncoproteomics-based framework to perform systematic investigation of the somatic mutations by altering phosphorylation sites and the related signaling networks in approximately 5,000 tumor samples across 16 major cancer types. After physically locating a total of 746,631 somatic mutations in protein sequences, we identified 8,057 unique missense mutations directly located at the direct phosphorylation sites and 97,277 mutations in the seven residues flanking regions around phosphorylation sites. In our comparison of the mutational load distribution of missense mutations at phosphorylation sites versus other sites of the protein sequences, we found a higher mutation rate at phosphorylation sites than that of the remaining protein sites in nine cancer types we examined. Furthermore, we found that missense mutations located around phosphorylation sites were more likely to be deleterious when compared to other residues. Based on the significance test implemented in our kinome-wide network module for cancer pharmacogenomics, KNMPx, we reported 47 proteins that harbored significantly mutated phosphorylation sites in pan-cancer analysis and 74 proteins in 16 individual cancer type analysis. Enrichment analysis indicated that these proteins were significantly enriched in two well-studied cancer gene data sets.

Recently, Creixell and colleagues identified six types of network-attacking mutations (NAM), including the changes in kinase and SH2 modulation, network rewiring, and the genesis and extinction of phosphorylation sites (19). They showed that signaling networks are both dynamically (e.g., mutant molecular logic gates) and structurally rewired (e.g., both constitutive activation and inactivation of kinase and SH1 domains) in cancer cells. They further demonstrated that the amino acids in the kinase domain encoding substrate specificity form sparse networks of nonconserved residues spanning distant regions (20). For example, they found that inter-residue allostery play important roles in specificity and an evolutionary decoupling of kinase activity and specificity. However, those two studies did not provide novel network algorithms and show the relationships of NAMs with clinical investigations (survival analysis) and drug responses. Compared with previous studies (19, 20), this study has two significant contributions. (i) Building on our observation of the clustering of phosphorylation site mutations in highly coexpressed kinase–substrate interactions, we proposed a new network algorithm to search for cancer-related modules consisting of a set of genes that are enriched with phosphorylation site mutations and tissue-specific kinase-substrate network. We demonstrated that network module genes were significantly enriched in well-known cancers for both individual cancer and pan-cancer analyses, as well as were associated with patient survival. (ii) We reported a kinome-wide landscape of cancer pharmacogenomic interactions through integrating somatic mutation-rewired signaling networks in approximately 4,997 primary tumors and 1,001 cancer cell lines using a statistical framework. Surprisingly, we found that mutation profiles in cell lines could highly reproduce the oncogenic phosphorylation site mutations identified in primary tumors (Fig. 5). We further demonstrated that many of these associations with drug sensitivity or resistance of various target therapeutic agents, including inhibitors targeting EGF, MAPK, PI3K, mTOR, and Wnt signaling pathways (Figs. 6 and 7). All dataset and predicted results in this study are free available at website: https://bioinfo.uth.edu/mutPhosphoProt. To our best knowledge, this is the first study of systematic oncoproteomics analysis for investigating phosphorylation site mutations mediating tumorigenesis and anticancer drug efficacy at the kinome-level. This study, although based on computational and statistical approaches, identified many promising mutation candidates for drug development in the rapidly evolving precision oncology.

Similar to many published computational works, our results might be influenced by potential data bias, data incompleteness, and confounding factors. For example, we only observed higher mutation rates at phosphorylation sites in nine of the 16 cancer types. This may be biased due to the incompleteness of the current phosphorylation data and the lack of sufficient number of samples. To overcome the potential bias from the phosphorylation sites, we further investigated the tissue-specific kinase–substrate network models, which may be potentially altered by phosphorylation site mutations. Although two tumors may not have any phosphorylation site mutations in common, they may share the biological networks affected by these mutations. By investigating phosphorylation site mutations in the context of network, we might be able to identify some rarely mutated proteins in the subnetworks linking to the well-characterized cancer-related proteins. This approach increased our power to detect potential biomarkers at the protein and network levels (46).

Another limitation is that this study only focused on missense mutations, which alter phosphorylation sites by amino acid substitution and excluded other types of mutations, such as nonsense mutations, insertions/deletions (indels) or gene fusion events. All these mutations may also play important roles in altering phosphorylation sites and leading to tumor initiation and progression. Moreover, our current study mainly focused on the distribution of phosphorylation site mutations. We did not deeply characterize the potential impact of mutations by gain-of-function or loss-of-function at the kinome-level. Integrating the functional impact of phosphorylation site mutations using high-throughput phenotyping assays with our present analysis framework will be one promising direction for our next-step work (47). In addition to phosphorylation signaling, dephosphorylation signaling by phosphatases also play critical roles in various diseases, such as cancer (48). In our future work, we will investigate both phosphorylation signaling by kinases and dephosphorylation signaling by phosphatases altered by somatic mutations, which may medicate tumorigenesis and anticancer drug responses.

No potential conflicts of interest were disclosed.

The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conception and design: J. Zhao, F. Cheng, Z. Zhao

Development of methodology: J. Zhao, F. Cheng, Z. Zhao

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Zhao, Z. Zhao

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Zhao, F. Cheng, Z. Zhao

Writing, review, and/or revision of the manuscript: J. Zhao, F. Cheng, Z. Zhao

Study supervision: F. Cheng

The authors would like to thank Dr. Pora Kim for developing the mutPhosphoProt website.

This work was partially supported by the NIH grant R01LM011177 to Z. Zhao.

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.