Uncovering Phosphorylation-Based Specificities through Functional Interaction Networks*

Protein kinases are an important class of enzymes involved in the phosphorylation of their targets, which regulate key cellular processes and are typically mediated by a specificity for certain residues around the target phospho-acceptor residue. While efforts have been made to identify such specificities, only ∼30% of human kinases have a significant number of known binding sites. We describe a computational method that utilizes functional interaction data and phosphorylation data to predict specificities of kinases. We applied this method to human kinases to predict substrate preferences for 57% of all known kinases and show that we are able to reconstruct well-known specificities. We used an in vitro mass spectrometry approach to validate four understudied kinases and show that predicted models closely resemble true specificities. We show that this method can be applied to different organisms and can be extended to other phospho-recognition domains. Applying this approach to different types of posttranslational modifications (PTMs) and binding domains could uncover specificities of understudied PTM recognition domains and provide significant insight into the mechanisms of signaling networks.

, Scansite (4), NetPhorest (3), GPS (9), KinasePhos (10), and many more. However, these methods depend on the availability of many known target sites for each modeled kinase. Here, we aim to tackle a more difficult challenge of predicting the specificity of kinases, without any direct knowledge of its experimentally determined target sites. The prediction tool Predikin (11) takes such an approach by trying to predict kinase specificity by examining 3D models of kinases bound to their substrate peptides. This analysis has identified residues in the kinases catalytic domain referred to as substrate determining residues (SDRs), which confer a preference for residues in the phosphosites flanking regions. Predikin uses these SDRs sites to match a new kinase sequence to a kinase with known specificity. In this way, Predikin also makes use of known kinase target site information. In addition, Predikin depends on the availability of protein structures and therefore cannot be easily scaled to kinase families without 3D structures nor to other PTM recognition domains.
We decided to take an alternative approach at predicting kinase specificity. Previous studies have shown that it is possible to use information regarding the interaction partners of a peptide-binding protein to identify potential motifs mediating these interactions (12). We reasoned that putative interaction partners of a kinase are more likely than random proteins to be phosphorylated by that kinase. Thus, phosphosites occurring on interaction partners of kinases should confer a bias in amino acid composition toward the kinase's specificity, which can be revealed by motif enrichment (Fig. 1). We tested this on human kinases to identify already known specificities, as well as other understudied kinases. We experimentally determined peptide targets of four understudied kinases and showed that predicted models closely resemble the experimentally identified sites. We extended our analysis to show that specificities can be predicted, not only for kinases, but also for other phospho-residue binding domains, such as 14-3-3 proteins and to an acetyl-lysine binding bromodomain. We also applied our method to mouse and showed that the predicted specificities of some kinases are conserved. We show here that it is possible to combine large-scale PTM data with protein network information to derive the specificity of PTM regulators and believe that this approach can be widely applicable to different PTM types.
Kinase Domain Prediction-Given a protein sequence, we used Kinomer (15), which uses multilevel hidden Markov models and HMMER (16) to identify protein kinases, and classify them into their appropriate kinase family. E-value cutoffs for each family were used as defined in (15). If a kinase was assigned more than one predicted family, the one with the highest E-value was used. These families were also used to determine if the kinase is Ser/Thr-specific or Tyr-specific. We assume that a kinase is either Ser/Thr-specific or Tyr-specific and do not account for dual specificity kinases.
Motif Enrichment-To identify motifs enriched within a set of phosphosites, compared with a background, we used the motif-x algorithm (17). Here, we used two background sets, defined as 10,000 15-mers centered on nonphosphorylated Ser/Thr or Tyr residues, depending on if the kinase is Ser/Thr-specific or Tyr-specific. All enrichments were carried out for Ϯ7 residues surrounding the central residue with occurrences Ն 10 and p Ͻ 10 Ϫ6 . Since the motif-x tool was only available via an online webserver, we reimplemented the tool for the R programming language, which can be found here: https://github.com/omarwagih/rmotifx.
Kinase Specificity Models-Specificity models were constructed as PWMs, which are commonly used to model specificities of linear motifs (18). PWMs can then be used to score peptides. We use an adapted version of the matrix similarity score (MSS), originally developed in the MATCH algorithm (19), as described in (20). The MSS ranges from 0 -1, where 0 represents no predicted binding, and 1 represents perfect predicted binding.
The performance of a given PWM was evaluated as the area under the ROC curve (AUC), which is the curve representing the relationship between the false positive rate and true positive rate as the MSS score cutoff is varied: Here, FP, TP, TN, FN represent the number of false positives, true positives, true negatives, and false negatives, respectively. The PWMs were used to score positive and negative sequences in order to generate these values. For a kinase of interest, we define the positive sequences as the set of phosphosites annotated to the kinase and the negative sequences as phosphosites annotated to any kinase not belonging to the same kinase family, as defined by Manning et al. (1).
In the case where the performance of experimental models were evaluated (i.e. using the gold-standard sequences), we performed 10-fold cross validation in which the kinase sequences are randomly split into 10 bins. Each bin is iteratively used as the test set, while the remaining nine are used to construct the PWM. This results in 10 AUCs, which are averaged to provide an unbiased proxy of the PWM's prediction power.
We supply a resource that contains all of the information used for the specificity predictions of each kinase. This can be accessed from http://evocellnet.github.io/kpred/. For each kinase, the user can find the specificity logo as well as the list of interacting partners and phosphosites used to generate the prediction. We also provide an R package to allow others to more easily use this method for their own PTM and species of interest. The package and a tutorial on how to use it are available via the help section of the prediction website.
Profiling In Vitro Kinase Substrates-Identification of in vitro kinase substrates was conducted as previously described (21). Briefly, lysate proteins were extracted from HeLa S3 cells at about 80% confluence in 15 cm dishes, and the total protein amount was measured by a BCA protein assay kit. Dephosphorylation was then carried out with TSAP (Promega, Madison, MI, USA) at 37°C for 1 h, and TSAP was inactivated by heating to 75°C for 30 min. For in vitro kinase reaction, each 100 g of dephosphorylated proteins (1 g/l) was reacted with 1 l of each recombinant kinase (0.5 g/l) or distilled water as a control at 37°C in kinase reaction buffer (40 mM Tris-HCl (pH 7.5), 20 mM MgCl 2 , 1 mM ATP) for ). The kinases were expressed as N-terminal GST-fusion protein using baculovirus expression system with SF9 cells and were purified using glutathione Sepharose chromatography. The reaction was stopped by heating to 95°C for 5 min. After protein reduction/alkylation, Lys-C/trypsin digestion (1/100 w/w) was performed and phosphopeptides were enriched by TiO2-based hydroxyl-acid-modified metal oxide chromatography (22).
Phosphopeptides were desalted by StageTips and analyzed by nanoLC-MS/MS using a self-pulled analytical column (150 mm length ϫ 100 m inner diameter) packed with ReproSil-Pur C18-AQ materials (3 m, Dr. Maisch, Ammerbuch, Germany). An Ultimate 3000 pump (Thermo Fisher Scientific, Germering, Germany) and an HTC-PAL autosampler (CTC Analytics, Zwingen, Switzerland) were used coupled to an LTQ-Orbitrap XL (Thermo Fisher Scientific). A spray voltage of 2,400 V was applied. The MS scan range was m/z 300 -1,500. The top 10 precursor ions were selected in MS scan by the Orbitrap with r ϭ 60,000 for MS/MS scans and the ion trap in the automated gain control (AGC) mode, where automated gain control values of 5.00 ϫ 105 and 1.00 ϫ 104 were set for full MS and MS/MS, respectively. To minimize repetitive MS/MS scanning, a dynamic exclusion time was set at 20 s with a repeat count of 1 and an exclusion list size of 500. The normalized CID was set to be 35.0. Mass Navigator v1.2 (Mitsui Knowledge Industry, Tokyo, Japan) with the default parameters for the LTQ-Orbitrap XL was used to create peak lists on the basis of the recorded fragmentation spectra. Peptides and proteins were identified by automated database searching using Mascot v2.3 (Matrix Science, London, UK) against SwissProt release 2010_11 (02/11/2010, 522,019 entries) with a precursor mass tolerance of 3 ppm, a fragment ion mass tolerance of 0.8 Da, and strict trypsin specificity allowing for up to two missed cleavages. Carbamidomethylation of cysteine was set as a fixed modification and oxidation of methionines; phosphorylation of serine, threonine, and tyrosine were allowed as variable modifications. Peptides were considered identified if the Mascot score was over the 95% confidence limit based on the "identity" score of each peptide and if at least three successive y-or b-ions with a further two or more y-, b-, and/or precursor-origin neutral loss ions were observed, based on the error-tolerant peptide sequence tag concept. After identification, phosphopeptides identified from the control samples were rejected. A randomized decoy database created by a Mascot Perl program gave a 1% false-discovery rate for identified peptides with these criteria. Phosphosite localization was evaluated using a site-determining ion combination method based on the presence of site-determining y-or b-ions in the peak lists of the fragment ions, which supported the phosphosites unambiguously.

Network-Based Prediction of Kinase-Substrate Specific-
ity-We hypothesized that the interaction network of a protein kinase should be enriched in its target proteins. This hypothesis was confirmed by the observation of a very significant enrichment of known kinase targets in the functional interaction or physical interaction partners of kinases (Supplemental Fig. 1). In order to predict kinase specificities, we then combined information on human protein interaction data and phosphorylation data derived from large-scale MS studies. Given that kinases bind their target proteins transiently, we used functional interactions derived from STRING (23) as a source for potential kinase interactors. We collected a total of 2,425,314 interactions in 22,523 proteins and compiled ex-perimentally determined phosphorylation sites from three public databases (PhosphoSitePlus (5), PhosphoELM (6) and HPRD (7)). Phosphosites were mapped back to proteins having information in STRING, resulting in 107,444 sites in 12,207 proteins. We identified 493 kinases in this reference proteome (81% serine/threonine and 19% tyrosine kinases) using the Kinomer prediction tool (15) (Methods). For a given kinase, all phosphosites occurring on the STRING partners of a kinase were collected and enriched for motifs using the motif-x algorithm (17) (Fig. 1, Methods). A random sample of 10,000 unphosphorylated Ser/Thr/Tyr sites were used as the background for enrichment. Phosphosites matching the most significant extracted motifs were then used to build a position weight matrix (PWM), which highlights the predicted specificity of the kinase and can be used to score novel phosphosites (Fig. 1). A survey of all known phosphorylation sites revealed a strong enrichment for prolines at position ϩ1 (Proϩ1) (Supplemental Fig. 2). This results in consistent enrichment of Proϩ1 motifs (Supplemental Fig. 3, Supplementary Results). To circumvent this, we require to know if a kinase is prolinedirected (i.e. prefers Proϩ1). We found that kinases of the CMGC family, including CDKs, MAPKs, GSKs, and CDK-like kinases have Proϩ1 motifs as shown in their experimental binding sites (Supplemental Fig. 4). Also, of all the non-CMGC kinases (with Ն20 known targets) only 1.57% (1/68) were found to be proline directed. Thus, if a kinase is not predicted as CMGC, phosphosites containing Proϩ1 are removed from foreground and background sets prior to motif enrichment.
There are two variable parameters in our method: the cutoff for the functional-interaction prediction score from STRING and the top k number of significant motifs extracted during the enrichment. To determine the best thresholds to use, we tested the predicted kinase specificity models against a set of 9,595 gold-standard kinase-substrate relationships. We carried out the benchmarking using a set of nine well-studied kinases from a diverse set of kinase families (ABL1, AKT1, ATR, AURKB, CDK2, CSNK2A1, GSK3B, MAPK1, and PRKACA) with well-recognized specificities in the literature. We varied the STRING cut-off, and the top k motifs extracted. The performance of the resulting PWM in each case was evaluated using the area under receiver operating characteristic (AUC) (Methods). We found that increasing the STRING score threshold, overall, resulted in higher AUCs (Supplemental Fig. 5). However, in most cases, we did not see a significant increase after a score cutoff of 400, and therefore, we used that cut-off throughout our analysis. We chose to select the top five motifs for two reasons. First, the AUCs among varied k values did not vary considerably. Second, overselecting motifs could mask the predicted specificity of the kinase (Supplemental Fig. 6). We also tested if difference sources of evidence (e.g. text mining, coexpression, interaction data) within the STRING database resulted in a different performance. Overall, using different evidence provided by the STRING database did not provide a significant increase in the AUCs, and a larger number of kinases can predicted using the combination of all evidence types (Supplemental Fig. 7).
Next, we checked to see how likely random models constructed without the network information performed in comparison to our predicted models. If a given kinase has n STRING interactions, and among those interactors there are m phosphosites (s 1 , s 2 , . . . , s m ), then m random phosphosites are selected from all known phosphosites in public databases (r 1 , r 2 , . . . , r m ). Specificity is then predicted, as previously described, using these sites. Random models, constructed using the top five motifs, were compared against the predicted models in their discriminative power against the gold-standard sequences, as measured by the AUC. We found that our predicted models for most of the nine kinases, with the exception of AKT1 and MAPK1, performed signifi-cantly better than random ( Fig. 2A and Supplemental Fig. 8). This does not mean that the AKT1 model is incorrect since it performs very well at predicting known AKT1 target sites (AUC ϭ 0.90, Fig. 2A). However, some kinases like AKT1 have specificities that are well modeled by the most represented motifs across all sites. Thus, in these cases, the network information appears to provide almost no gain compared with random sampling. In opposition to these kinases, ATR has a specificity that is very uncommon with a preference for glutamine at position ϩ1 that is very well recovered by this approach (Figs. 2A and 2F) but very unlikely to be observed in a random pool of phosphosites.
These results demonstrate the ability to integrate protein interaction information with large-scale data on protein phosphorylation to derive kinase specificity models.
Prediction of Kinase Specificity across All Human Kinases-Our method was applied to all kinases, resulting in predictions for 282/493 (57%) of kinases (Supplementary data). Kinases that did not result in a prediction either had a low number of partners or a scarcity of phosphosites on partners. We selected 85 kinases with Ն20 known target phosphosites as well as a prediction and compared how well the predicted models performed with respect to the kinase family (Fig. 2B). The average AUC across all kinases was 0.64 with 32% (27/85) of kinases having an AUC greater than 0.7 (Fig. 2B). We found that CMGC, PIKK, and AGC families performed best, whereas TKL, STE and TK kinases had a larger fraction of poorly performing models. Excluding the TKL, STE, and TK kinases, the average AUC increases to 0.68 with 44% of kinases (27/61) scoring higher than 0.7 (Fig. 2B). We speculated that the differences in performance across the different kinase families reflect different degrees of specificity in the kinase-substrate recognition. For example, many tyrosine kinases have additional targeting domains (i.e. PTB and SH2 domains). Also, several STE kinases are known to have an additional interaction surface for a "docking motif" (24,25). For these kinases, targeting is achieved by multiple interfaces or often aided by other mechanisms (see Discussion), and therefore, they might be less specific in the recognition of sequences around the phosphosite. In line with this reasoning, we found that kinases that had additional protein domains also had worse-performing models (p Ͻ 1.88 ϫ 10 Ϫ3 , Wilcoxon signed-rank test; Supplemental Fig. 9). We tested this notion more explicitly by comparing our predicted models with a proxy for kinase promiscuity. For the same set of kinases, we built PWMs using the gold-standard target sites and computed an experimental AUC using 10-fold crossvalidation (Methods), which reflects how well the gold-standard models perform at distinguishing their own sites. Interestingly, we found that AUCs of predicted models showed a strong correlation to that of the experimental (r ϭ 0.757, p Ͻ 2 ϫ 10 Ϫ16 ; Fig. 2C), suggesting that kinases with high true specificity are more likely to have high predictability.
We selected a few example kinases to demonstrate the specificity determinants captured by our approach (Figs. 3D-3G). The predicted specificities strongly resemble that of the experimental and, in many cases, are much more apparent. For example, the pronounced preference for glutamine at ϩ1 for ATR is recovered (Fig. 2F). The known Akt preference for an arginine at -3 and -5 is fully recovered in the predicted model. In addition, there is a more apparent preference for [ArgSerThr] at -2, and arginine at -1 that is not as apparent in the experimental specificity (Fig. 2G).
In attempt to identify features comprised by better performing models, we searched for relationships between the AUC of the predicted model and (1) the number of functional interacting partners, (2) the number of phosphosites on interacting partners, (3) the distribution of information content, and (4) the number of extracted motifs. We observed weak correlations (r Ͻ 0.361, Supplemental Fig. 10) for each of the individual features. However, we were able to achieve a higher correlation by combining a number of features using a linear regression model (r ϭ 0.542, p ϭ 8.37 ϫ 10 Ϫ8 ). This model can thus be used to assign a quantitative measure of confidence re-lated to the truth of predicted specificity, which we use to rank our predictions (Supplementary Data).
We tested several alterations of the method, such as using different background sets for motif enrichment as well as using only high confidence phosphosites, and overall did not observe a strong improvement in the AUC (Supplementary Results, Supplemental Figs. [11][12]. For example, to obtain a list of higher confidence sites, we tried excluding phosphosite positions that are supported only by one study or excluding from the analysis highly abundant proteins (Supplemental Fig.  11). To test the impact of removing the Proϩ1 sites, we used an alternative strategy whereby we did not filter Proϩ1 peptides but used as background all phosphosites instead of phospho-acceptor residues. Although we were still able to retrieve many correct predictions, the performance was lower in this alternative implementation.
Mass-Spectrometry-Based Validation of Kinase Specificity-We selected four kinases with few known target sites in the literature, across different kinases subfamilies for which we had a predicted model (CMGC kinases SRPK2 and HIPK2, AGC kinase AKT2, and PEK kinase EIF2AK4). For each of these kinases, we identified in vivo target sites using the phosphoproteomic approach described by Imamura et al. (21) (Fig. 3A). Briefly, HeLa cell extracts were treated with phosphatase to remove any existing phosphosites, and kinases were added in separate experiments. Phosphorylated extract were then subjected to trypsin digestion, phosphopeptide enrichment, and nanoLC-MS/MS (Fig. 3A). We identified a total of 483 novel phosphosites for these kinases (AKT2, n ϭ 248; EIF2AK4, n ϭ 91; HIPK2, n ϭ 106; SRPK2, n ϭ 38, available in Supplementary Data). The identified target sites were then compared with our predicted models for these kinases (Figs. 3B-3E). We found that all predicted models performed significantly better than random, and three of the four had an AUC Ն 0.7 at classifying the experimentally identified sites (Fig. 3F). These results are in line with the benchmarks performed and further support the validity of the approach described here. We note that the SRPK2 kinase was predicted to have a strong preference for serines and arginines at several positions. This motif was unusual given previously described models, though several elements of this motif are confirmed by the experimental sites (Fig. 3D). The kinase specificity of SRPK2 was very recently determined using a chemical genetic approach (26) that provides further validation for the predicted specificity model for this kinase.
Prediction of PTM Binding Specificities-To demonstrate the extendibility of our method to other types of linear motif specificities, we applied the same method to 14-3-3 proteins. 14-3-3 proteins are conserved single domain proteins capable of binding a phospho-serine or threonine and are responsible for tight regulation of several important pathways such as cell death, cell cycle control, and signal transduction (27). Previous studies have shown that binding sites these proteins have specificities toward their target phosphosites (28) (Fig.  4A). We applied the method to all seven human 14-3-3 proteins (Figs. 4A-4C, Supplemental Fig. 13) and similarly show that the recovered models are very good predictors of known binding sites (AUCϾ0.80) and perform significantly better than random (Fig. 4D). We recovered the well-known determinants such as arginine at position -3 and some preference for proline at position ϩ2. We found that there is little to no overlap between sites used to construct the models for each 14-3-3 protein, despite showing similar predicted specificity (Supplemental Fig. 14). This suggests that the same motif is recovered in each case from a different source of partner sites, adding to the confidence of the recovered models.
To highlight the broader usefulness of the method, we selected p300 that contains a bromodomain that binds acetylated lysines and its binding specificity has been well characterized (29). We then obtained a collection of human lysine acetylation sites and used the same network-based motif enrichment to predict the specificity of p300's bromodomain.
We found that the predicted specificity (Supplemental Fig. 13) is very similar to the known preference for KXXK or KXXXK (where X is any amino acid and both lysines are acetylated).
These results again show that PTM recognition specificity can be predicted by combining network information with PTM data.
Conservation of Kinase-Substrate Specificity-We next sought to predict kinase specificities in a different organism. We applied the same approach to mouse (Mus musculus), which contained 29,732 phosphosites and 2,425,424 STRING interactions. Using human kinases that had an AUC Ͼ 0.6, we identified one-to-one ortholog kinases in mouse, using the InParanoid resource (30). By applying our method to these kinases, we found a close resemblance of the specificity determinants of human and many of their corresponding mouse orthologs (Fig. 5). For 56 mouse kinases analyzed in this way, 19 (34%) showed a similar or better performance at predicting the known human kinase sites than the orthologous human model. This suggests that at least these 19 kinase pairs have very conserved kinase preferences. For the remaining cases, we cannot confidently say that there is a divergence in specificity since we cannot rule an incorrect prediction. DISCUSSION The advances in MS have expanded tremendously our knowledge of exact protein modifications sites for a number of different PTM types. However, there is almost no information regarding the regulatory interactions connecting regulators to target proteins. Determining the recognition preferences for PTM enzymes and binding domains in large scale is still an open problem and remains a limiting factor in achieving this goal. We used phospho-regulation as a model system and showed that it is possible to combine PTM information with interaction network data to derive accurate models of enzymes and binding domains. Predicted kinase motifs for 59% of human kinases are provided in supplementary material. In addition, a resource that contains all of the information used for the specificity predictions of each kinase can be accessed from http://evocellnet.github.io/kpred/. The code required to apply this approach can be found in the help page along with a tutorial.
We note that even though some models do not perform better than models created by random sampling of sites, this does not necessarily reflect the reliability of the predicted model. Some kinase specificities are well modeled by the most common motifs that are recovered from a random sample. For these cases, the added information from the network data does not result in a model that is more accurate than random. The power of our approach is therefore more obvious for regulators that have specificities that are less common such as the DNA damage kinase ATR. For this kinase, the recovered model is both accurate (AUC ϭ 0.94) and performs

FIG. 4. Prediction of 14-3-3 domain specificities. (A)
The specificity of 14-3-3 domains as constructed by experimentally verified substrates. (B and C) Prediction of specificities for two 14-3-3 proteins. Each example shows a logo representing the predicted specificity (left) and the top five extracted motifs and the number of phosphosites matching them (right). (C) The performance of each predicted model compared with models predicted using random phosphosites. All cases perform better than random sampling of phosphosites (.p Ͻ .01, *p Ͻ .05, **p Ͻ .01, *** p Ͻ .001, one-sided z-test). Error bars represent the median absolute deviation of 1,000 random models. much better than models produced by random sampling of sites.
In the current implementation of this approach, we assume that kinases that are not CMGC tend not to be proline directed and remove Proϩ1 phosphosites. This may result in mispredicting cases where a non-CMGC kinases is proline directed and also cases where CMGC kinases are proline directed. We tested an alternative approach that does not require the removal of Proϩ1 phosphosites, but this resulted in a lower overall performance. However, we note that, even when using Proϩ1 peptides, we can still obtain predictions that do not have Proϩ1. For example, the CK2a1 kinase is a casein kinase and therefore part of the CMGC group. For this group, we allow Proϩ1 peptides to be included in the predictions and we still obtain a strong bias for an acidic residue at the ϩ1 position. Additionally, we require to know the class of the kinase: serine/threonine or tyrosine to filter only phosphosites matching the class of the kinase. This is because phosphotyrosine is in many regards a different PTM from phosphoserine and phosphothreonine. In particular, it occurs at much lower frequency, so if we would not discriminate between these two types, the predicted specificities would be dominated by phospho-Ser/Thr.
It is important to take into account that most phosphosite information was retrieved from phosphoproteomics experi-ments that have used trypsin for protein digestion. Given that trypsin cleaves C-terminal to arginine and lysine residues, it is very possible to expect a bias for Arg/Lys residues in the phosphopeptides. However, we do not think this bias is a strong influence on the recovered motifs. If it was a strong influence, we would expect any bias to be equally possible at positions before or after the target site and also not specifically biased for Arg or Lys. Instead, arginine determinants are more frequent than Lys determinants and Arg determinants are not symmetrically distributed. Of the 202 Arg determined positions (defined as having Ͼ0.25 relative frequency at the position), 96% (194/202) are found before the phosphosite and 0.039% (8/202) are found after the phosphosite. There are only 19 positions where Lys is the major determinant, and these tend to be more distributed with 42.1% (8/19) occurring before and 57.89% (11/19) occurring after the phosphosite.
Different kinase families show different average performance in their predictions and that the degree of kinase specificity for the known target sites is correlated with the accuracy of the predicted models. These observations highlight the inherent limitation of the approach proposed here. PTM-interacting proteins that recognize their target sites mostly by residues flanking the target site will be more amenable to this approach than those that use multiple recognition mechanisms. These include docking motifs, colocalization, coex- pression, and scaffolding interactions [see (2) for a full review]. In addition this approach assumes that the recognition occurs in a linear epitope at the PTM position. However, it has recently been shown that kinase targeting can occur also in a 3D epitope instead of a linear motif (31). If a PTM enzyme or binding domain often recognizes the target site by a 3D epitope, then this linear motif enrichment strategy will not be appropriate. These observations should be taken into account for future use of this method for other PTM recognition domains. We show that this method can be applied to different modes of site-directed motif-binding domains, such as 14-3-3 domains and bromodomains, suggesting that the method could thus be extended further to analyze specificities of other PTM recognition domains. Finally, we applied this approach to study the conservation of kinase specificity between human and mouse kinases. For the cases that we analyzed, at least 34% appear to have conserved specificity. We suggest that this approach, in combination with an analysis of potential mutations in specificity-determining residues, could be used to identify PTM recognition domains with diverged specificities across species. Given that these regulators interact with many different target PTMs, it is expected that their specificity diverges slowly. This is in contrast to the fast changes in the PTMs targeted by these proteins that can diverge quickly (32,33). Conserved regulator specificity with diverged target sites is a scenario that is analogous to what is observed in transcriptional regulation (34). However, there have been cases described for divergence of transcriptionfactor specificity (35), so analogous cases of divergence of PTM recognition are likely to exist. In addition to studying the evolution of specificity, applying this method to different organisms could lend further confidence to the true specificity of a PTM recognition domain since models trained in different species could contribute complementary specificity determinants and ultimately be combined to provide better models.
In summary, we describe here a novel approach to predict PTM recognition motifs, and we believe it can be applicable to a wide range of recognition domains and contribute significantly to our understanding of these signaling systems.