INTRODUCTION
Protein phosphorylation is a kind of post-translational modification and has been shown to be one of the most essential regulatory and signaling mechanisms in the cell (
Zolnierowicz and Bollen, 2000). The process is catalyzed by protein kinases, in which the γ phosphate on ATP or GTP is transferred to the substrates. In eukaryotic cells, phosphorylation usually takes place on Serine (S), Threonine (T) or Tyrosine (Y) of the substrate protein. The phosphate on substrates can be removed by phosphatases, so the phosphorylation process is reversible: it is determined by the balance between the protein kinases and phosphatases. This reversible character allows the phosphorylation process to work like a switch in a living cell. Specific substrates could be activated by protein kinases under the simulation of an external signal. After the signal wanes, the activated substrates could be inactivated by phosphatases and wait for the next signal. Phosphorylation can regulate a variety of important protein functions, including subcellular localization, protein degradation and stabilization, as well as biochemical activities (
Cohen, 2000;
Ficarro et al., 2002;
Manning et al., 2002a;
Zannini et al., 2012). There are usually a series of phosphorylation processes involved in a normal biological function
in vivo (
Ubersax and Ferrell, 2007;
Cai et al., 2012). It was also implicated in various pathological processes, such as cancer (
Finn and Lu, 2008;
Ollila and Makela, 2011, insulin resistance (
Tanti and Jager, 2009), polyglutamine disease (
Zhou et al., 2008) and Alzheimer's disease (
Chung, 2009).
In an eukaryotic cell, about 30%–50% of the proteins can be phosphorylated (
Pinna and Ruzzene, 1996). There are also hundreds of kinases within an eukaryotic genome, for instances: 518 protein kinases in humans (
Manning et al., 2002b), 540 kinases in mice (
Caenepeel et al., 2004) and 251 kinases in Drosophila (
Morrison et al., 2000). The enzymes must be specific and act only on a defined subset of cellular targets to ensure signal fidelity. Cells must have a mechanism to control the phosphorylation process involving so many kinases and protein substrates simultaneously and precisely. This mechanism is mainly realized by the specific recognization of protein kinases to substrates, which determines the exact time and place for phosphorylation to occur. Thus, the identification of the involved kinases and their phosphorylation sites is the first step to understand the mechanism.
Currently there are a number of computational methods for phosphorylation site prediction. Generally, these methods can be divided into two categories: non-kinase-specific and kinase-specific phosphorylation site prediction. For non- kinase-specific phosphorylation site prediction, there exist NetPhos (
Blom et al., 1999), DISPHOS (
Iakoucheva et al., 2004), PHOSIDA (
Gnad et al., 2007), etc.. For kinase-specific phosphorylation site prediction, there exist GPS (
Xue et al., 2010), NetPhosK (
Blom et al., 2004), KinasePhos (
Wong et al., 2007), PPSP (
Xue et al., 2006), etc.. Although these methods could predict whether S/Y/T sites could be phosphorylated or not, they cannot predict whether a phosphorylation site is functional or not, which is a major issue for further experimental researches.
For the phosphorylation sites with known functions, it had been demonstrated that they are under strong functional constraints and are evolutionarily more conserved than those with no characterized functions (
Landry et al., 2009;
Ba and Moses, 2010). So the evolutionary information of S/Y/T sites can be incorporated to identify the most likely functional phosphorylation sites. Since most phosphosites occur in disordered regions and the conservation of phosphorylation sites is also influenced by the region in which the residue is located (
Landry et al., 2009), it is necessary to consider the relative conservation of an S/Y/T site against its flanking region. Based on these considerations, we developed a prediction method that incorporated both absolute and relative conservation information of S/Y/T sites to facilitate the identification of the most likely functional phosphorylation sites.
RESULTS
Web server development
We developed a web server, which incorporated NetPhos (
Blom et al., 1999) and NetPhosK (
Blom et al., 2004) to predict general and kinase-specific phosphorylation sites. Users can also upload prediction results from other phosphorylation site prediction tools as a formatted table (the table template can be downloaded from the web server). PhosphoSitePlus (
Hornbeck et al., 2012) was incorporated to mark whether a predicted phosphorylation site is experimentally validated or not. Then, the server calculated the absolute and relative conservation score for each possible phosphorylation site with Rate4site. Both scores were normalized to the range of 0-1, and the larger the relative and absolute conservation score is, the more likely the phosphorylation site is functional. The access of the web server is “http://lifecenter.sgst.cn/ppps/en/home.do”. Using the web server, users could predict most likely functional phosphorylation sites using UniProt ID, the query box is shown as Fig. 1, and the query result is shown as Fig. 2. The web server could also search the predicted results in the constructed human, rat and mouse database and the query box is shown as Fig. 3. Users could also search the database by a specific kinase and the query box is shown as Fig. 4.
Human, rat and mouse database construction
We selected the phosphorylation sites predicted by both NetPhos and GPS (
Xue et al., 2010), and the conservation scores could be calculated to construct the human, rat and mouse database. For human, rat and mouse proteome, the number of the phosphorylation sites predicted by NetPhos, GPS or both, and the statistics of the final database are shown in Table 1. The “Database sites” column indicates the selected predicted phosphorylation sites in the database. The “Protein sequences” column indicates the number of protein sequences containing the selected predicted phosphorylation sites in the database. The “Experimentally validated sites” column indicates the number of experimentally validated phosphorylation sites in the database.
Analysis of relative and absolute conservation scores of human, rat and mouse database
We compared the density distribution of the relative and absolute conservation score of human, rat and mouse database using R stats package (Fig. 5). The density distribution of the relative conservation score of human, rat and mouse are consistent, the density distribution of the absolute conservation score are also consistent. For both relative and absolute conservation score density distribution, there are two peaks: one at about 0.025 and one at about 0.80. So our method not only can predict which phosphorylation sites are most likely to be functional, but also can give clues to which phosphorylation sites are least likely to be functional, thus can help relevant researchers to select more conserved and important phosphorylation sites to perform further studies. The density distribution of the relative and absolute conservation score intersect at about 0.85. For conservation score larger than 0.85, the distribution density of the relative conservation score is larger than the absolute conservation score. It may be explained that some phosphorylation sites are more conserved against its flanking region than against the overall protein sequence. For the majority of conservation score less than 0.85, the distribution density of the absolute conservation score is larger than the relative conservation score, indicating that some phosphorylation sites are less conserved against their flanking region than against the overall protein sequence. Previous studies showed that the conservation of phosphorylation sites is also influenced by the region in which the residue is located (
Landry et al., 2009), the density distribution of the absolute and relative conservation score also demonstrated that it is necessary to consider both the relative and absolute conservation of an S/Y/T site against its flanking region and the overall protein sequence, respectively.
General prediction results
We used the selected proteins as mentioned in the Materials and methods section to do the KEGG pathway enrichment analysis. The results for human, rat and mouse are shown in Table 2, Table 3 and Table 4, respectively. We found that the functions of protein phosphorylation in all the enriched KEGG pathways are supported by previous studies.
Enrichment analysis for human
For the human proteome, we selected a total of 1755 proteins containing 2834 predicted phosphorylation sites (the selection criteria of top 10% is 0.918232). We matched the 1755 selected UniProt protein IDs to their gene ids using the R package org.Hs.eg.db and used all the gene ids in the org.Hs.egUNIPROT table within this R package as the background. The cutoff of the P value in the enrichment analysis was set to 0.001. There are 25 enriched KEGG pathways (Table 2), in all of which protein phosphorylation has been demonstrated to play important roles, as shown in the References column of Table 2.
Enrichment analysis for rat
For the rat proteome, we selected a total of 85 protein sequences containing 107 predicted phosphorylation sites (the selection criteria of top 10% is 0.923121). We matched the 85 selected UniProt protein IDs to their gene ids using the R package org.Rn.eg.db and used all the gene ids in the org.Rn.egUNIPROT table within this R package as the background. The cutoff of the P value in the analysis was set to 0.001. There are 8 enriched KEGG pathways (listed in Table 3) for the selected rat proteins. The important roles of protein phosphorylation in all of these KEGG pathways had been supported by previous studies.
Enrichment analysis for mouse
For the mouse proteome, we selected a total of 555 protein sequences containing 794 predicted phosphorylation sites (the selection criteria of top 10% is 0.917137). We matched the 555 selected UniProt protein IDs to their gene ids using the R package org.Mm.eg.db and used all the gene ids in the org.Mm.egUNIPROT table within this R package as the background. The cutoff of the P value in the analysis was set to 0.001. There are 8 enriched KEGG pathways (Table 4) for the selected mouse proteins. The important roles of protein phosphorylation in all of these enriched KEGG pathways have been supported by previous studies.
Protein specific prediction results
We used two well-studied proteins, p53 and Cyclin-dependent kinase inhibitor 1B, in which phosphorylation plays an important role, to demonstrate the usefulness of our method for individual protein phosphorylation studies.
p53
Our method totally predicted 23 phosphorylation sites in p53 (Table 5). We ranked the predicted phosphorylation sites by their relative conservation scores. Within these 23 sites, 14 sites have been experimentally validated to be phosphorylated. And according to the annotation of UniProt (Version 196), 9 phosphorylation sites have been supported to be functional. Our method predicted 5 of these 9 functional phosphorylation sites (site 15, 46, 392, 315 and 9). The ranks of the relative conservation score of these 5 sites were 2, 3, 5, 6 and 19, respectively. p53 serine 15 phosphorylation could direct its interaction with B56γ and the tumor suppressor activity of B56γ-specific protein phosphatase 2A (
Shouse et al., 2008). p53 serine 46 could be phosphorylated by HIPK2 upon UV irradiation, which could regulate p53 apoptotic activity and is required for acetylation by CREBBP (
D'Orazi et al., 2002;
Hofmann et al., 2002;
Chang et al., 2005;
Lee et al., 2009). Phosphorylation at serine 9 by HIPK4 could increase the repression activity of p53 at p53 repressive promoters (
Arai et al., 2007). Phosphorylation of serine 392 stabilizes the tetramer formation of tumor suppressor protein p53 and could stimulate the DNA-binding ability of p53 (
Sakaguchi et al., 1997;
Kapoor et al., 2000). Phosphorylation of p53 at serine 315 after irradiation damage could stimulate p53-dependent transcription (
Blaydes et al., 2001).
We can see that 4 of these 5 sites were within the top 6 sites ranked by the relative conservation score. For site 9, it may be explained that the function of site 9 phosphorylation is relatively less important for biological activities and a previous study has demonstrated that the specific recognition of Ser9 appears to be dependent upon additional determinants of p53 beyond the N-terminal 65 amino acids (
Soubeyrand et al., 2004). But for site 9, we can also find that the relative conservation score (0.3185200) is much larger than the absolute conservation score (0.1658500), indicating it is more conserved against its flanking region than against the overall protein sequence.
Cyclin-dependent kinase inhibitor 1B
For cyclin-dependent kinase inhibitor 1B, our method predicted totally 19 phosphorylation sites (Table 6), within which 10 sites have been experimentally validated. According to the annotation of UniProt (Version 138), a total of 5 phosphorylation sites have been supported to be functional. Our method predicted 3 of these 5 functional phosphorylation sites, i.e. site 187, 10 and 198. The rank of the relative conservation score of these 3 sites were 2, 6 and 12, respectively. Phosphorylation of threonine 187 leads to protein ubiquitination and proteasomal degradation (
Boehm et al., 2002;
Fujita et al., 2002;
Motti et al., 2004;
Hao et al., 2005;
Sabile et al., 2006). Phosphorylation of serine 10 is the major site of phosphorylation in resting cells, which takes place at the G(0)-G1 phase and leads to protein stability (
Boehm et al., 2002;
Fujita et al., 2002;
Motti et al., 2004). Phosphorylation of threonine 198 is required for interaction with 14-3-3 proteins (
Fujita et al., 2002,
2003;
Motti et al., 2004). The relative conservation scores of all these three sites were larger than 0.7. The high relative conservation score of other sites may be explained by the possibility that the function of these sites may not have been studies or these sites may not work by phosphorylation directly. However there were also 6 predicted phosphorylation sites with relative conservation scores less than 0.2. It may give further studies a clue that researchers could pay less attention to these sites than those having higher conservation scores.
DISCUSSION
In this work, a prediction web server was developed to facilitate the identification of the most likely functional protein phosphorylation sites by incorporating both the absolute and relative evolutionary conservation scores. The larger the relative and absolute conservation score is, the more likely the phosphorylation sites is functional. To facilitate the usage of our method, we also selected and integrated two existing computational methods: NetPhos and NetPhosK for general and kinase-specific phosphorylation site prediction, respectively. Using our method, we predicted the most likely functional sites of the human, rat and mouse proteomes and built a database for the predicted phosphorylation sites. By the analysis of overall prediction results, we demonstrated that protein phosphorylation plays an important role in all the enriched KEGG pathways. By the analysis of protein-specific prediction results, we also demonstrated the usefulness of our method for individual protein studies. Our method would help to characterize the most likely functional phosphorylation sites for further studies in this research area.
MATERIALS AND METHODS
Web server development
In our pipeline, we first predicted all possible phosphorylation sites for a protein with NetPhos (
Blom et al., 1999) and NetPhosK (
Blom et al., 2004), which are two existing computational methods for non-kinase-specific and kinase-specific prediction, respectively. The prediction results from other phosphorylation site prediction methods can also be provided as a formatted table (the table template can be downloaded from the web server). We incorporated PhosphoSitePlus (
Hornbeck et al., 2012), which is a database containing experimentally validated phosphorylation sites, to mark whether a predicted phosphorylation site is experimentally validated or not.
For conservation score calculation, we used customized Rate4Site with default parameters (
Pupko et al., 2002;
Mayrose et al., 2004), which can compute the evolutionary rate for each site in a multiple sequence alignment. The alignments of ortholog families were downloaded from the NCBI HomoloGene database (
Sayers et al., 2012). For the calculation of the absolute conservation score, we normalized the evolutionary rate r of a phosphorylated site according to the rates of all the residues in the protein, i.e.
Where μ(rall) and σ(rall) are the mean and standard deviation of the evolutionary rates of all residues. For the calculation of relative conservation score, we normalized r according to the rates of the flanking residues around the phosphorylated site, 5 to the left and 5 to the right, i.e.
We transformed the zabs and zrel scores to [0, 1] by using the probability function of the standard normal distribution.
Human, rat and mouse database construction
We downloaded the proteome sequences of human, rat and mouse from UniProt (
Consortium, 2012). Using our method, we predicted the phosphorylation sites of these proteomes. To guarantee the prediction accuracy of our method, we took the phosphorylation sites predicted by both general and kinase-specific methods and having the conservation scores to construct the human, rat and mouse database. We then incorporated PhosphoSitePlus (
Hornbeck et al., 2012) to mark whether the predicted sites have been experimentally validated or not.
General and protein-specific prediction
Since the function information of phosphorylation sites is limited, it is difficult to construct a benchmark dataset to test the overall prediction performance of our method. We ranked the predicted phosphorylation sites in the human, rat and mouse database according to the relative conservation score and selected the top 10% of the experimentally validated phosphorylation sites. Then we selected the protein sequences containing these top 10% sites and did KEGG enrichment of such proteins to find whether protein phosphorylation plays an important role in the enriched KEGG pathways.
We used two well-studied proteins, p53 and Cyclin-dependent kinase inhibitor 1B in which phosphorylation plays an important role, to demonstrate the usage of our method for individual protein studies.
Higher Education Press and Springer-Verlag Berlin Heidelberg 2012