Genome-wide analysis to predict protein sequence variations that change phosphorylation sites or their corresponding kinases Open Access

Abstract

We define phosphovariants as genetic variations that change phosphorylation sites or their interacting kinases. Considering the essential role of phosphorylation in protein functions, it is highly likely that phosphovariants change protein functions. Therefore, a comparison of phosphovariants between individuals or between species can give clues about phenotypic differences. We categorized phosphovariants into three subtypes and developed a system that predicts them. Our method can be used to screen important polymorphisms and help to identify the mechanisms of genetic diseases.

INTRODUCTION

Protein phosphorylation is involved in various important processes: development and learning at the organism level, and the cell cycle, differentiation and apoptosis at the cellular level ( 1 , 2 ). Phosphorylation can change the subcellular localization of a protein, its life span and its affinity for other proteins or DNA ( 3 ). Therefore, the addition or deletion of phosphorylation sites through phosphovariants can lead to functional variations in proteins that can result in phenotypic variations or genetic diseases. By our definition, phosphovariants are variations that change phosphorylation sites or their interacting kinases. We propose three subtypes of phosphovariants. First, some variations occur directly at phosphorylation sites, and these sites will be removed if the phosphoreceptors are replaced with amino acids other than serine, threonine or tyrosine. Conversely, replacement of another amino acid with a serine, threonine or tyrosine may add a new phosphorylation site. Second, variations adjacent to phosphorylation sites can result in the removal or addition of phosphorylation sites. Third, variations may change the kinases that recognize phosphorylation sites, without changing the phosphorylation site itself. We divided phosphovariants into type I, II and III, respectively, according to the above descriptions (Figure 1 ).

We developed PredPhospho (version 2), a web-based computer program that predicts phosphorylation sites, and PhosphoVariant, a database for human phosphovariants. Even the advanced laboratory techniques used to analyze phosphorylation sites, such as mass spectrometry (MS), cannot analyze all types of proteins ( 4 , 5 ). For example, peptides that are either too small or too large in mass can be easily missed. Moreover, membrane proteins cannot be obtained in sufficient quantities for analysis ( 5 ). Even when proteins can be analyzed with MS, it is very time consuming and expensive to make thousands of variant proteins and select the phosphovariants. PredPhospho can predict the phosphorylation sites in kinase-specific ways, using the support vector machines (SVMs) derived from statistical learning theory proposed by Vapnik and Chervonenkis in 1995 ( 6 ). In our study, we searched for known phosphovariants and tried to predict other possible phosphovariants among human variations.

METHODS

PredPhospho

We created classifiers of various kinases by training SVMs with phosphorylation site sequences and nonphosphorylated site sequences. In other words, our classifiers determine whether serine, threonine or tyrosine residues within a sequence can be phosphorylated or not. ‘Phosphorylated site sequences’ refers to peptide sequences with a serine, threonine or tyrosine residue located at the center, and which are phosphorylated. Conversely, ‘nonphosphorylated site sequences’ are sequences with a serine, threonine or tyrosine residue located at the center, which have not been found, yet, to be phosphorylated. We obtained phosphorylated site sequences from public databases: the Swiss-Prot (release 54.8) and the Human Protein Resource Database (HPRD, release 7). Nonphosphorylated site sequences were taken from laboratory data confirmed by MS (see Supplementary Data ).

Manning et al. ( 7 ) found 518 human protein kinase genes in the human genome sequence, using the hidden Markov model (HMM) profile, and confirmed the identities of more than 90% of the identified kinase genes using cDNA cloning. They also classified the protein kinase superfamily into nine broad groups, and subdivided the groups into 134 families and 204 subfamilies, using sequence comparisons of the kinase catalytic domains. We classified the phosphorylated site sequences according to their kinases and created the classifiers in a kinase-specific manner. Because of the limitations of the phosphorylated sequence data presently available in public databases, we can make classifiers for only six kinase groups: AGC, CAMK, CK1, CMGC, STE and TK; and 18 kinase families: AKT, CAMK2, CAMKL, CDK, CK1, CK2, GSK, IKK, JakA, MAPK, PDGFR, PIKK, PKA, PKC, RSK, Src, STE20 and Syk (all abbreviations are shown in the footnote to Supplementary Table S1 ). The detailed algorithms and methods were described in the Supplementary Data .

Evaluation of the system

The performance status of the prediction for each kinase group and family is shown in Supplementary Table S2 . The performance of the prediction with combinations of all the kinase group models or all the family group models is not the numerical multiplication of the performance of each model. Therefore, to evaluate the performances of the predictions at the kinase group level or at the family level, we tested two proven real data sets, which were compiled with MS experiments. Data set I was created by Olsen et al. ( 5 ), who identified phosphorylation sites in proteins from HeLa cells and classified the phosphorylation sites according to their definition. Four classes are based on their localization probabilities: <0.25, 0.25 ≤ <0.75 without kinase motifs, 0.25 ≤ <0.75 with kinase motifs and ≥0.75. We used only monophosphopeptides that had localization probabilities of at least 0.75. Data set II was derived from the paper of Beausoleil et al. ( 4 ) and we used those of their phosphopeptides with a localization certainty of >99.4% (see Supplementary Figure S2 ). We selected phosphorylation sites and nonphosphorylated sites from these two data sets. To avoid overestimating the performance of PredPhospho, we discarded sequences that were more than 70% identical to sequences used for training PredPhospho to sequences used for training PredPhospho ( Table 4 ). We also avoided using false nonphosphorylated sites by omitting those sites that are listed as phosphorylation sites in Swiss-Prot or HPRD. We tested the two kinds of data sets not only with our PredPhospho, but also with Scansite.

Prediction of phosphovariants

We extracted information about human genetic variations from SwissVariant of the Swiss-Prot database. SwissVariant includes single amino acid polymorphisms and missense mutations ( 8 ). The number of variations listed in SwissVariant was 33 651. We consulted the Swiss-Prot database and HPRD about their effects, and the references to these variations and phosphorylation sites. With PredPhospho and Scansite, we predicted the phosphorylation sites and related kinases for the original sequences and the variant sequences. The phosphovariants could be identified when the phosphorylation sites or interacting kinases were altered between the original sequence and the variant sequence. If the phosphorylation site is in the same location as the variation, it is type I. In type II phosphovariants, the variation is not in the same location as the phosphorylation site. We added the symbol (+) to types I and II when the phosphovariants added new phosphorylation sites, and (–) when the phosphovariants removed phosphorylation sites [e.g. type I (+) or type I (–)]. Type III phosphovariants are caused by changes in the types of kinases involved, rather than in the phosphorylation site itself, regardless of the locations of the variations. One variation can include more than one class of phosphovariant, because one variation can affect two or more phosphorylation sites. We predicted phosphovariants at the kinase group level and the family level. The predictions at the family level are more sensitive, but less specific, than those at the group level. To minimize false negatives, we varied the specificity options (95%, 97%, 98% or 99%) according to the specificity of each model ( Supplementary Table S3 ). The specificity options are described in the Supplementary Data .

Sequence logos

We obtained 562 phosphorylation site sequences recognized by the CMGC kinase group from Swiss-Prot and HPRD. We trimmed the sequences as six symmetric residues centered phospohrylation sites. We aligned the sequences and obtained a sequence logo using the web program ( http://weblogo.berkeley.edu/logo.cgi ).

WWW programs

The PredPhospho (version 2) and PhosphoVariant were implemented using the PERL (version 5.8.8) programming language and MySQL (version 5.0.18). The PhosphoVariant is a database for the definite and possible human variants changing phosphorylation sites and their interacting kinases. They are available at: http://phosphovariant.ngri.go.kr/seq_input_predphospho2.htm and http://phosphovariant.ngri.go.kr , respectively.

Bug fixing and minor changes of database will be done, whenever needed. The new version of Swiss-Prot and HPRD will be incorporated into our database, annually at least.

RESULTS

Type I phosphovariants

The substitution of phosphoreceptor amino acids with amino acids other than serine, threonine or tyrosine causes the elimination of phosphorylation sites and can be classified as type I (–) phosphovariants according to our classification. We found 50 type I phosphovariants by matching the locations of the variations and those of phosphorylation sites registered in the Swiss-Prot database and the Human Protein Resource Database (HPRB, Table 1 ) ( 9 , 10 ). Of these phosphovariants, 19 are known to cause Mendelian-inherited diseases and 18 are associated with cancers. Another 13 phosphovariants are polymorphisms.

Conversely, new phosphorylation sites can be created by variations and we defined these as type I (+) phosphovariants. We found an example in the study of Nousiainen et al. ( 11 ) The cell line that they used had a Gly766Ser mutation in the probable ATP-dependent RNA helicase DDX27 (Swiss-Prot ID, Q96GQ7), which was identified as phosphorylated. Similarly, the polymorphisms in Tables 1–3 are good examples of the addition of phosphorylation sites by sequence variations. For example, if we postulate that the isoleucine at amino acid 571 of CTP synthase 1 (Swiss-Prot ID, P17 812) is changed to serine, then it also represents a type I (+) phosphovariant, rather than a type I (–) phosphovariant ( Table 1 ).

Type II phosphovariants

It is more difficult to find examples of type II phosphovariants than to find type I phosphovariants, because we cannot definitely say that a phosphorylation site is changed by a substitution near a phosphorylation site. However, when some kinases (although not all kinases) recognize phosphorylation sites, the specific amino acids near the phosphoreceptor are important. In Figure 2 , we present sequence logos of phosphorylation site sequences for the CMGC kinase group. The proline residues at position +1 relative to the phosphorylation sites are important in the phosphorylation site sequences of the CMGC kinase group, especially the CDK kinase family and the MAPK kinase family. Most (84%) of the phosphorylation sites of the CMGC group registered in the Swiss-Prot database and the Human Protein Resource Database (91% of the CDK and 87% of the MAPK kinase families) have +1 proline residues. If the proline is substituted with another amino acid, it is highly probable that the adjacent phosphorylation site will be abolished. The phosphorylation site at Ser112 of peroxisome proliferator-activated receptor gamma protein (Swiss-Prot ID, P37231) is eliminated by the Pro113Gln substitution ( 12 ). We found three other polymorphisms that abolish phosphorylation sites of the CMGC kinase group, but the removal of these phosphorylation sites has not yet been confirmed ( Table 2 ). The presence of specific amino acids does not directly affect phosphorylation by kinases other than those of the CMGC kinase group, but sequences near the phosphorylation site must be considered. Kinases recognize the residues surrounding the target phosphorylation site, and the amino acids bordering phosphorylation sites are, in turn, affected by other nearby residues ( 13 ). Hence, when the relevant kinases are not members of the CMGC kinase group, it is difficult to predict type II phosphovariants simply by database matching, without specific programs that predict phosphorylation sites.

Type III phosphovariants

Type III phosphovariants are those variations that change only the type of kinase involved, without affecting the phosphorylation site itself. For example, Ser386 of the tyrosine protein phosphatase, nonreceptor type 1 (PTPN1, Swiss-Prot ID, P18031) is phosphorylated by cell division cycle 2 (CDC2) kinase, a member of the CMGC kinase group, and by casein kinase 2 (CK2), a member of the ‘Other’ kinase group ( 14 , 15 ). The Pro387Leu substitution reduces 75% of the phosphorylation by CDC2 in vitro ( 16 ). However, it has not been confirmed that Pro387Leu inhibits the recognition of Ser386 by CK2. Only about 5% of the sites phosphorylated by CK2 that are registered in the Swiss-Prot database or HPRD have a proline residue at position +1 relative to the phosphorylation site. There is also no known consensus sequence for CK2 that contains proline at that location ( 17 ). Therefore, we infer that the proline residue is not essential for phosphorylation by CK2 and that Pro387Leu will have little effect on phosphorylation by CK2. Therefore, we consider Pro387Leu of PTPN1 a type III phosphovariant because it inhibits the recognition of Ser386 by CDC2 kinase but has little effect on its phosphorylation by CK2 ( 16 , 18 ).

Kinases that recognize serine and threonine differ from the kinases that recognize tyrosine. The substitution of phosphorylated serine or threonine for tyrosine, or vice versa , can remove a phosphorylation site or change the type of kinase that recognizes it. Changes between serine and threonine can also cause changes in the phosphorylation site and the responsive kinase. Therefore, the phosphovariants in Table 3 (Panel b) are either type I or type III phosphovariants.

Performance of PredPhospho and Scansite

We developed prediction models for six kinase groups and 18 kinase families. Their accuracy ranged from 70.80 to 94.67% at the kinase family level and from 71.77 to 91.18% at the kinase group level ( Supplementary Table S2 ). We tested our prediction models using two real laboratory data sets compiled with MS. For the six kinase group models, the sensitivities were 79.40% with data set I and 75.47% with data set II, but the specificities were as low as 60.62% for data set I and 61.04% for data set II because of the accumulation of false negatives by all six kinase group models. When we modified the specificity to >95% (see Supplementary Data for the modification of the specificity for each model), the specificities increased to 72.09% and 72.39%, respectively, for each data set, whereas the sensitivities decreased to 73.24% and 65.76%, respectively ( Table 5 ). At a specificity of >99%, the specificities changed to 95.79% and 96.62%, respectively, and the sensitivities to 23.39% and 20.05%, respectively. Scansite is a widely used web-based prediction software for phosphorylation sites ( 19 ). We also tested Scansite with the same data sets. When we applied Scansite with the low-stringency option to the data sets, the specificities were 52.60% and 57.06%, respectively, and the sensitivities were 84.47% and 83.92%, respectively, whereas with the high-stringency option, the specificities were 96.77% and 95.71%, respectively, and the sensitivities were 16.39% and 13.60%, respectively ( Table 6 ). The performance of PredPhospho with no modification to the specificity is similar to that of Scansite used with the low-stringency option. The performance of PredPhospho with a specificity of >99% was similar with that of Scansite used with the high-stringency option. The data sets were analyzed in the same way with the 18 kinase family models ( Tables 5 and 6 ). The family-wise prediction generally had greater sensitivity and lower specificity than the group-wise prediction.

Phosphovariants predicted with PredPhospho and Scansite

The numbers of phosphovariants predicted with PredPhospho and Scansite are shown in Table 7 . In the Supplementary Data , we present the results for phosphovariants that were predicted with PredPhospho with the >99% specificity option. The sensitivity and specificity of the prediction of type I phosphovariants will be the same as the result shown in Table 5 , because we can predict type I phosphovariants simply by the location of phosphorylation sites, with no knowledge of the kinds of kinases involved. The notions of type II and type III phosphovariants include the kinds of kinases that recognize the phosphorylation sites. Therefore, the prediction of type II and III phosphovariants differs from that of type I. Not only the phosphorylation site but also the type of kinase must be identified, because important amino acid residues flanking the phosphorylation sites, which guide kinases to the site, may differ according to the kinase involved. Therefore, the general performance of our prediction of type II and type III phosphovariants will be somewhat different from those shown in Table 5 . Instead, the performance of the each kinase-specific prediction can be judged from the performances shown in Supplementary Table S2 .

The proportion of phosphovariants predicted by PredPhospho is shown in Table 8 . These data were selected with the >99% specificity option at the kinase family level, and are related to the confirmed phosphorylation sites in human or orthologous proteins, for which kinase information is not yet available. As described in the Supplementary Data , only phosphorylation sites with available kinase information were used as training models for PredPhospho. Hence, the phosphorylation sites in Table 8 are predicted ones with PredPhospho, because the phosphorylation sites were not included in the training data for PredPhospho. If a specific site is shown to be a phosphorylation site in human or orthologous proteins, then the site is definitely located on the surface of the protein, and therefore, is accessible to a kinase. The predicted phosphovariants related to a proven phosphorylation site are more likely to be true than are those related to an unproven phosphorylation site. Numerous phosphorylation sites in humans or other species are constantly being identified. The priority for further research among predicted phosphovariants can be decided based on the confirmation of specific phosphorylation sites.

DISCUSSION

Changes in phosphorylation sites cause various diseases by numerous mechanisms. Some proven mechanisms of the phosphovariants shown in Tables 1–3 are related to changes in the protein's affinity for DNA, inducing hyperphosphorylation and the inhibition of ubiquitinization. For example, microphthalmia-associated transcription factor (MITF) activates the transcription of the tyrosinase gene. The Ser405Pro change in MITF eliminates the phosphorylation site at Ser405 and inhibits the binding of MITF to DNA. As a result, the mutation causes Waardenburg syndrome type IIa, the symptoms of which include depigmentation and sensorineural hearing loss ( 20 ). Abnormal phosphorylation can also cause disease by increasing phosphorylation at other sites. The hyperphosphorylation of tau protein induces neurofibrillary tangles and the accumulation of these tangles can result in Alzheimer's disease and frontotemporal dementia (FTD). Paradoxically, serine threonine protein kinase N (PKN) interrupts the phosphorylation of other sites by phosphorylating Ser637 and Ser669 of tau protein ( 21 ). The mutations Ser637Phe and Ser669Leu of tau protein eliminate the recognition sites for PKN and induce the hyperphosphorylation of tau protein. FTD and the respiratory failure with dementia are known to be related to Ser637Phe and Ser669Leu, respectively ( 22 , 23 ). Phosphorylation at Ser32 of NF-κB inhibitor α (Swiss-Prot ID, P25963) causes ubiquitinization and results in the activation of NF-κB ( 24 ). The Ser32Ile substitution of NF-κB inhibitor α (Swiss-Prot ID, P25963) eliminates the phosphorylation site at Ser32. Consequently, the ubiquitinization of NF-κB inhibitor α is inhibited by the Ser32Ile variant of NF-κB inhibitor α and NF-κB cannot be activated. The Ser32Ile variant of NF-κB inhibitor α causes autosomal dominant anhydrotic ectodermal dysplasia with immunodeficiency ( 25 ). These are a few examples of phosphovariants. Considering the numerous functional roles played by phosphorylation in vivo , there must be many mechanisms by which phosphovariants can cause specific diseases, and these must be identified.

Some phosphovariants do not cause Mendelian-inherited diseases, but change an individual's susceptibility to disease. The Pro113Gln substitution (dbSNP id, rs1800571) of peroxisome proliferator-activated receptor gamma (Swiss-Prot ID, P37231), which eliminates the phosphorylation site at Ser112, is known to cause obesity ( 12 ). The Pro387Leu substitution (dbSNP ID, rs16995309) of tyrosine-protein phosphatase nonreceptor, type 1 (Swiss-Prot ID, P18031) is associated with type II diabetes mellitus ( 16 , 18 ). We only found these two polymorphic phosphovariants that are related to disease susceptibility. As shown in Tables 1–3 , 21 phosphovariants are polymorphisms, and the biological significance of 19 of these polymorphisms is not yet known. In addition, we do not know the biological significance of most of the polymorphic phosphovariants predicted in our study. Considering the importance of phosphorylation in protein function, polymorphic phosphovariants may well be involved in specific diseases or phenotypes.

Apart from the characteristics of the three types of phosphovariants already suggested, there are other fundamental differences between the type I and other phosphovariants. Type I phosphovariants completely add or remove a phosphorylation site because kinases can only donate a phosphor moiety to an amino acid with a hydroxyl group. However, type II and III phosphovariants can significantly affect kinase kinetics without completely changing the kinase's recognition site. For example, the Pro387Leu substitution of PTPN1 removes only 75% of the phosphorylation of Ser386 by CDC2 kinase in vitro , rather than 100% ( Table 3 ) ( 16 ).

We can explain phenotypic variations and diseases in terms of phosphovariants in more cases than we have anticipated. Of the human proteins registered in the Swiss-Prot database, 25.5% are phosphoproteins and 60.9% of these phosphoproteins have multiple phosphorylation sites. The protein with the greatest number of confirmed phosphorylation sites is the serine/arginine repetitive matrix protein 2 (Swiss-Prot ID, Q9UQ35), with 195 phosphorylation sites. Although we could only determine 62 phosphovariants with a database search, many more phosphovariants must exist. Furthermore, if we count the mutations that add phosphorylation sites and the type II and III phosphovariants that cannot be found without specific programs, the number of mutations associated with changes in phosphorylation sites is much greater. We did not consider haplotypes in this study, for simplicity. However, if two or more nearby variations are frequently linked, a nearby phosphorylation site can be altered, although each variation individually does not affect the phosphorylation site. Therefore, phosphovariants may account for a much greater proportion of human variation than we have anticipated.

Several issues must be resolved in future studies. First, we regarded nonannotated serine, threonine and tyrosine residues as nonphosphorylated sites. Although we tried not to select false nonphosphorylated sites by using real MS data sets, some of the nonphosphorylated sites might be determined to be phosphorylated sites, in the long run. Second, type II and type III phosphovariants can be interchanged. For example, removed phosphorylation sites associated with type II (–) phosphovariants predicted with PredPhospho may be recognized by kinases that are not included in our prediction models. In such cases, these variations are type III phosphovariants, not type II (–). Conversely, if a phosphorylation site is falsely classified as a site recognized by multiple kinases, instead of by one true kinase, and if a variation affects only some of these kinases, this is a type II phosphovariant incorrectly predicted to be a type III phosphovariant. Other points that must be improved are the low sensitivity achieved with the high-specificity option and the low specificity achieved with the no-specificity option of PredPhospho ( Table 5 ). Moreover, the number of types of kinases that we can predict must be increased and haplotypes should be considered in future studies.

Our method can be used in pathophysiological studies of mutations and in the selection of polymorphisms of clinical and phenotypical importance. Many of the papers that have described the variations, shown in Tables 1–3 , did not mention that the variations could be related to changes in phosphorylation sites. This could be attributable to the lack of a specific database that connects mutations with phosphorylation sites, or the lack of a general understanding of the association between phosphorylation and mutation. The type I (+), II and III phosphovariants we have defined cannot be identified simply by database analyses. Specific programs are required to identify these phosphovariants. Accordingly, many nonsense point mutations whose functional mechanisms are unknown can be reconsidered in terms of phosphovariant. Furthermore, if some mutations are predicted to be phosphovariants with our system, further research will clarify the cause of the associated disease or protein function. Our system can be used to select meaningful variations among endless numbers of newly identified polymorphisms. As sequencing techniques advance, a large number of genetic variations are emerging. At present, comparison of whole genomes of individuals is possible, because the human genome can be sequenced in two months ( 26 ). A comparison of phosphovariants between individuals or between species can be undertaken before amino acid variations or nucleic acid variations are compared in whole genomes. A reverse genetic approach for unknown protein functions or phenotypic variations is possible with proven phosphovariants. The screening and prediction of phosphovariants can be a starting point for further research.

FUNDING

This work was supported by Health Fellowship Foundation. Funding for open access charge: intramural grants from Korea National Institute of Health, the Korea Center for Disease Control of the Republic of Korea (Project No.: 091-4845-301-210).

Conflict of interest statement . None declared.

REFERENCES