Computational prediction and analysis of species-specific fungi phosphorylation via feature optimization strategy

Abstract

Protein phosphorylation is a reversible and ubiquitous post-translational modification that primarily occurs at serine, threonine and tyrosine residues and regulates a variety of biological processes. In this paper, we first briefly summarized the current progresses in computational prediction of eukaryotic protein phosphorylation sites, which mainly focused on animals and plants, especially on human, with a less extent on fungi. Since the number of identified fungi phosphorylation sites has greatly increased in a wide variety of organisms and their roles in pathological physiology still remain largely unknown, more attention has been paid on the identification of fungi-specific phosphorylation. Here, experimental fungi phosphorylation sites data were collected and most of the sites were classified into different types to be encoded with various features and trained via a two-step feature optimization method. A novel method for prediction of species-specific fungi phosphorylation-PreSSFP was developed, which can identify fungi phosphorylation in seven species for specific serine, threonine and tyrosine residues (http://computbiol.ncu.edu.cn/PreSSFP). Meanwhile, we critically evaluated the performance of PreSSFP and compared it with other existing tools. The satisfying results showed that PreSSFP is a robust predictor. Feature analyses exhibited that there have some significant differences among seven species. The species-specific prediction via two-step feature optimization method to mine important features for training could considerably improve the prediction performance. We anticipate that our study provides a new lead for future computational analysis of fungi phosphorylation.

Introduction

Reversible phosphorylation [1], one of the most important post-translational modifications (PTMs), plays a key role in diverse cellular processes including cell division [2], apoptosis [3], immune response [4] and signal transduction pathway [5]. The process of protein phosphorylation primarily occurs at amino acid residues such as serine (S), threonine (T) and tyrosine (Y) in substrate protein. However, not all of the three residues can be phosphorylated. Thus, the efficient identification of phosphorylation sites is very significant to the understanding of intracellular mechanism. To date, various experimental methods have been used for the identification of phosphorylation sites. Especially, the advent of high-throughput mass spectrometry (MS) [6] has rapidly accumulated fungi phosphorylation data. Furthermore, there is much evidence to indicate that fungi phosphorylation also occurs in more complex organisms [7–9], but its abundance and role remain largely unknown, primarily due to the lack of adequate analytical tools [10]. Indeed, since the discovery of phosphorylation in 1906 by Phoebus A. Leve [11], most researches have been made to analyse the pattern of phosphorylation sites in plants [12, 13], animals [14] and bacteria [15]; only limited advances have been made in the detection of fungi phosphorylation and the clarification of its biological function.

Compared with conventional experimental approaches, in silico prediction of phosphorylation sites has become a promising strategy, which could be conducted as preliminary analyses and screens out the potential targets for further confirmation. Although there are excellent methods to analyse protein phosphorylation at S/T/Y residues, few are suitable for the study of fungi phosphorylation. From the research, we also observed that understanding the molecular mechanisms of different species of fungi phosphorylation proteins was essential to regulate fungal pathogenicity of fungal biology. For example, for Cryptococcus neoformans (C. neoformans), one of the essential factors for its virulence is due to activation of cyclic adenosine 5′-monophosphate (cAMP)-dependent protein kinase A (PKA) pathway. PKA is known to phosphorylate various down-stream targets including transcription factors such as Nrg1 and Rim101, through which cAMP/PKA pathway modulates capsule size and melanin formation in C. neoformans. Identification of phosphoproteins in C. neoformans could serve as a discovery platform to establish a molecular network that confers virulence in C. neoformans [16].

Here we used the Fungi Phosphorylation Database (FPD) provided by Bai et al. that comprises 62 272 non-redundant phosphorylation sites in 11 222 proteins across 8 fungal species [17], including Aspergillus flavus (A. flavus) [18], Aspergillus nidulans (A. nidulans) [19], C. neoformans [16], Fusarium graminearum (F. graminearum) [20], Magnaporthe oryzae (M. oryzae) [21], Neurospora crassa (N. crassa) [22], Saccharomyces cerevisiae (S. cerevisiae) [23] and Schizosaccharomyces pombe (S. pombe) [23] to preprocess, and proposed a computational tool of PreSSFP via a two-step feature optimization method for the prediction of specific fungi phosphorylation sites. At first, five types of feature extraction strategies were used to formulate our protein peptide fragments based on the sequence information, evolutionary information and physicochemical properties. Then we performed a two-step feature selection to remove the redundant feature vectors and construct the optimization model on the basis of 10-fold cross-validation. By comparison, PreSSFP exhibited a competitive performance to other existing tools for predicting specific fungi phosphorylation sites. The fungi phosphorylation sites were classified into different species and PTMs sites for training, which could further improve the prediction performance. Taken all together, PreSSFP is provided as an open source tool and implemented in Matlab and JAVA at http://computbiol.ncu.edu.cn/PreSSFP.

Materials and methods

There are four main procedures to construct and evaluate PreSSFP system as follows (Figure 1): (i) construct valid benchmark dataset and independent dataset to train and test the PreSSFP for seven organisms separately; (ii) formulate the peptide samples in the datasets with an effective mathematical expression by extracting their various sequence, evolutionary and physicochemical properties feature; (iii) perform a two-step feature optimization algorithm based on the benchmark datasets to train the optimal feature subsets for each organism in a cross-validation manner; and (iv) objectively evaluate the optimal model for each organism by using the independent test datasets and develop the PreSSFP webserver.

Data collection and preprocessing

In this work, experimentally verified fungi phosphorylation data for eight model organisms, including A. flavus, A. nidulans, C. neoformans, F. graminearum, M. oryzae, N. crassa, S. cerevisiae and S. pombe datasets, were collected from publicly available FPD [16] and UniProt database [24]. Both A. flavus and A. nidulans are belonging to Aspergillus, so we consolidated them into a category of Aspergillus. Therefore, we obtained seven categories data in fungi. The homology protein sequences were removed with a 30% identity cutoff using CD-HIT [25]. For seven model organisms, independent datasets were constructed by randomly selecting 15% of all non-redundant protein entries and the remaining non-redundant proteins were used to construct training dataset. Experimentally verified phosphorylation serine (S) / threonine (T) / tyrosine (Y) residues were regarded as positive samples. All the remaining S/T/Y residues that had not been verified as phosphorylation sites in these proteins were considered as negative samples. Each site was represented as a peptide segment of the length of 15 with S/T/Y in the center. Based on that, a total of 94 proteins with 152 serine phosphorylation sites and 6762 non-phosphorylation sites were obtained as an independent dataset for Aspergillus in this study; the remaining 528 proteins containing 930 serine phosphorylation sites and 37420 non-phosphorylation sites were utilized as training samples (the other different species-specific phosphorylation serine, threonine and tyrosine datasets are shown in Table 1). From Table 1, we find that the number of tyrosine phosphorylation in all species except S. cerevisiae is too small and lacks statistical significance. Therefore, we regarded the tyrosine phosphorylation dataset in other six organisms as the second independent test dataset for cross-species prediction of tyrosine phosphorylation sites. Furthermore, the amount of serine phosphorylation data in S. cerevisiae is too big, so we divided it into 10 averages for training. Finally, the average value of all 10 training results is adopted. In addition, the number of negative samples is greater than positive samples, and a positive-to-negative sample ratio of 1:1 is randomly pooled as a training dataset and independent test dataset.

Features extraction and optimization

Sequence information feature

where |${f}_i$| represents the frequency of i-th amino acid pair in|$\Big\{ AA, AC, AD,\cdots, YY\Big\}$|⁠.

is used to encode peptide |$\mathrm{S}$|⁠, where |${p}_j\Big(j=1,2,\cdots, n\Big)$| denotes the posterior probability of each amino acid at the j-th position in the positive training samples, and |${p}_j\Big(j=n+1,n+2,\cdots, 2n\Big)$| denotes the posterior probability of each amino acid at the j-th position in the negative training samples.

Evolutionary information feature

where |$a$| and |$b$| are the two amino acids; |$M$| is the BLOSUM62 substitution matrix; |$Sim$| stems from the normalized amino acid substitution matrix; |$L$| denotes the number of up-stream or down-stream amino acids flanking each side of the target S/T/Y; and |$\mathit{\max}(M)$| and |$\mathit{\min}(M)$| represent the largest and the smallest number in the matrix, respectively.

Physicochemical properties feature

For biochemical reactions, physicochemical property (PSP) is the most instinctive feature and is successfully applied in the bioinformatics field [32, 33]. The properties of each of the 20 amino acids are omnifarious, which can be answerable for the diversity and specificity of protein function and structure. A large number of experimental and theoretical studies have been carried out to show different varieties of properties of individual amino acids and to represent them based on the numerical index. Amino Acid index (AAindex) database (Version 9.1) contains 544 amino acid indices and specifies the physicochemical properties of amino acids. Then, according to the values connected with each PSP, the amino acids around the phosphorylation sites can be encoded. All of the 544 physicochemical properties are examined with descending order on the basis of the value of information gain (IG) in seven species separately. The larger the value of IG is, the greater the impact of the corresponding amino acid residue on the phosphorylation site is. Thus, in our work, in terms of the IG values of all PSP, the top nine are selected and defined as informative features in the prediction model.

Two-step feature optimization

A two-step feature optimization was employed to select prominent feature vectors and then to construct the optimal feature subset. In the first step, we calculated the importance of each feature vector by using random forest (RF) classifier in Python [34–36], which ranked each input feature according to the mean accuracy of the given test data and labels. When we made use of the RF classifier, the parameter n_estimators was first optimized. For n_estimators, the larger the better, but also the longer it will take to compute. In addition, note that results will stop getting significantly better beyond a critical number of trees. Thus, we defined it in the range of 50 to 250 via grid search with cross-validation to obtain the best n_estimators on the basis of accuracy (the detailed information of n_estimators in different model is shown in Table S1). Then we used the best parameter to execute model classification and acquired the top 200 features as the optimal feature candidates after this step. The second step was step-wise feature elimination. First, the feature vectors ranking top 10 in the first step were merged as a training model to calculate the performance based on support vector machine (SVM). At each round, the next feature of the first-step feature list was added into the model and then assessed. If SVM accuracy of the model is improved by its addition as a feature or the SVM accuracy of the model is over 80%, the feature is accepted as the final feature subset to build the prediction model; otherwise, it is eliminated (the detail in Algorithm Two-step feature optimization).

Performance evaluation

where TP, TN, FP and FN denote the number of true positives, true negatives, false positives and false negatives, respectively. The receiver operating characteristic (ROC) curves and area under ROC (AUC) values were also carried out.

Results

Current progresses in the prediction of eukaryotic phosphorylation sites

It is believed that up to one-half of eukaryotic proteins undergo phosphorylation [39]. Thus, most of the previous studies were mainly focused on eukaryotic proteins phosphorylation, including animals, plants and fungi. Table 2 provides a simple overview of currently available eukaryotic protein phosphorylation site prediction methods, including the type of species, tool name, kinase or kinase family-specific prediction, feature information, prediction algorithm, reference and website. From Table 2, 47 phosphorylation site prediction tools have been developed. We have broadly divided them into two categories: tools for general phosphorylation site prediction and tools for species-specific phosphorylation site prediction. However, the tool for fungi-specific phosphorylation site prediction is not available due to lack of data. There have only tools for predicting yeast phosphorylation site, including NetPhosYeast [38] and PhosphoPICK [76]. Meanwhile, it can be seen that the biggest difference among these computational methods is the feature information and prediction algorithm. As shown in Table 2, most of the prediction algorithms for identifying protein phosphorylation sites were machine learning methods, including artificial neural networks (ANNs), SVM, RFs, logistic regression and Bayesian network, etc. For example, the first tool, NetPhos 3.1, was built to predict eukaryotic proteins phosphorylation sites based on ANN by Blom et al. in 1999 [40]. Kim et al. implemented PredPhospho 1.0 in SVMs algorithm to predict for four PK groups and four PK families. And its enhanced version-PredPhospho 2.0 could predict for 7 PK groups and 18 PK families [70, 71]. As we can see from Table 2, owing to the strong generalization ability and versatility of SVM, it is the most widely used prediction algorithm for identifying phosphorylation sites.

On the other hand, in Table 2, many different feature infor-mation and encoding methods have been used to formulate the protein or peptide samples, including sequence information, structural information, physicochemical properties, evolutionary information and functional information. Among them, sequence information is the most widely used. Some predictors, like NetPhos, used the whole protein sequence information to input training model, or other predictors, such as GANNPhos and AMS, used AAC or binary encoding to encode amino acid sequence fragment for model training. Phosphorylation sites have been observed to be preferentially located in disordered regions [90]. In 2004 and 2010, predicted disorder scores for phosphorylation sites were used as features in the phosphorylation predictor DISPHOS [43] and Musite [69]. It is true that disorder scores of the experimentally verified phosphorylation sites are rare, and most of them require predictive tools to predict disorder scores, but doing so will bias the results of the prediction model. Meanwhile, the use of physicochemical properties could play a role in the phosphorylation site prediction tools. For example, the predictor PHOSFER used 24 AAIndex indices to encode protein sequence [89]. It mainly adopted consensus fuzzy clustering method to cluster 544 amino acid properties into 8 clusters and then identified a set of 24 indices consisting of three individual AAIndex indices from each cluster. The results suggested it was necessary to extract important information from 544 AAIndex for encoding phosphorylation sequences. Evolutionary information, such as KNN, was very effective feature when used for predicting phosphorylation sites [69, 83]. Li et al. designed Quokka with KNN and sequence information feature in 2018, which provided a high-quality prediction [83]. Furthermore, in order to improve predictive performance of the model, some tools utilized multiple information integration to predict phosphorylation sites, including PostMod [59] and Musite [69]. The combination of features can lead to a certain improvement, while they have not discussed the importance and contribution of different features for the phosphorylation sites prediction. Many of the derived features may contain redundant information and not all of the features play an equally important role in the prediction. Thereby feature selection is very necessary to select the most informative features and their contribution to the phosphorylation site prediction should be discussed in detail.

Fortunately, lots of fungi phosphorylation sites have been identified in the past 6 years, which provides a great opportunity to design a tool for fungi phosphorylation analysis. In view of the above, by incorporating sequence information, physicochemical properties and evolutionary information that had been evidenced good performances for prediction PTMs sites in the previous tools with two-step feature optimization strategy, we proposed PreSSFP for predicting species-specific fungi phosphorylation sites.

Feature optimization results via two-step feature optimization

In the previous research [26, 91, 92], we found that the results of independent coding tests with single feature were less than the combination of all five types of features (AAC, DAAC, BPB, PSP and KNN). So we combined all five types of features as a model for access. However, the prediction performances of the model were not fully satisfactory. It may be owing to all features that are not equivalently essential for the model performance. Or some of feature vectors may be unwanted noise, which easily leads to bias in the prediction performance. Thus, it is generally indispensable to execute feature optimization so that we can reserve the important one. As we know, the principal components analysis (PCA) [93] and maximum relevance minimum redundancy (mRMR) [78] are usually used for feature selection. PCA is aimed at selecting the data with higher cumulative contribution rate after processing and transformation of the original data. The mRMR method is based on mutual information for quantifying relevance and redundancy. In Figure 2, taking serine phosphorylation in Aspergillus, tyrosine phosphorylation in S. cerevisiae and threonine phosphorylation in C. neoformans as an example, we used the PCA and mRMR to carry out feature selection, respectively, and compared their results with two-step feature optimization. From Figure 2, we found that prediction performance of the model by using two-step feature optimization is superior to that of PCA and mRMR based on the feature vectors of the same dimension. For example, the accuracy has improved by 5% and 4% for serine phosphorylation of Aspergillus in comparison with mRMR and PCA separately, and the accuracy has increased at a rate of 17% and 15% for threonine phosphorylation in C. neoformans, compared with mRMR and PCA, respectively.

Here we applied two-step feature optimization to dig the principal feature vectors when generating the serine phosphorylation sites predictive model of Aspergillus. To start with, we had a total of 580 dimension feature vectors by combining all five types of features. Then we set the parameter n_estimator in the range of 50 to 250 for acquiring the best n_estimator. By using gridsearchcv, we obtained that the best n_estimator was 190 for prediction Aspergillus serine phosphorylation. Next, we received the ranking of the feature importance via RF classifier on the basis of 190. The lower the ranking of the feature vector is, the greater the impact on the serine phosphorylation sites is. Accordingly, we obtained the top 200 feature vectors as the optimal feature candidates for the second step feature optimization. Finally, we selected the feature vectors ranking top 10 from the 200 feature vectors to construct an initial model and calculated the performance. Whereafter, the remaining 190 feature vectors were put into this model to train one by one. If SVM accuracy of the model has been improved, the feature vector is preserved to construct a new model for the next training. Otherwise, it is eliminated. From what has been described above, with the dimension decreasing from 580 to 146, the accuracy increased from 77.90% to 81.18% and the AUC value added from 77.41% to 80.65% in serine phosphorylation sites prediction for Aspergillus. Comparatively speaking, the chosen feature vectors have superior prediction performance to the combination of five types of features. Totally, 15 models are built, including 7 serine phosphorylation models, 7 threonine phosphorylation models and 1 tyrosine phosphorylation model. The prediction performance of combining five types of features and the selected feature vector by two-step feature optimization for predicting other species-specific fungi phosphorylation protein sites are shown in Table 3.

Analysis of feature importance and contribution

Sequentially, we further analysed which feature vectors are valuable to prediction model based on optimization features selected by two-step feature optimization. In Figure 3, taking serine phosphorylation in Aspergillus, tyrosine phosphorylation in S. cerevisiae and threonine phosphorylation in C. neoformans as an example, for phosphoserine in Aspergillus, we reconstituted 143-dimension new feature vectors from above five types of features, and the ratios of dimension of chosen feature vectors belonging to the five types of features are 0.04 (15/400), 1 (5/5), 0.72 (91/126), 0.75 (21/28) and 0.67 (14/21). For phosphothreonine in C. neoformans, 168-dimension new feature vectors were selected from the above five types of features, and 0.06 (27/400), 1 (5/5), 0.80 (101/126), 0.82 (23/28) and 0.57 (12/21) were the proportions of dimension of chosen feature vectors accounting for the five types of features. In view of these results, we observed that the ratio of KNN feature was notably higher than those of other four features, indicating that KNN features exerted a vital effect on the performance evaluation of the model and contributed to identifying phosphorylation S/T/Y sites. Meanwhile, from the RF classifier results, the importance ranking of KNN feature in seven species is also higher than those of other features. Correspondingly, KNN features exhibited the best performances in these models, which could be concerned with conserved residues and detected local sequence similarity. Phosphorylation-related clusters may exist in local sequences around phosphorylation sites [69], and KNN can capture the cluster information and distinguish them from the background. Therefore, KNN score is suitable to be used as a feature for predicting phosphorylation sites. Then, from Figure 3, the ratios of the PSP and BPB features were also relatively high, implying that the PSP and the frequency of each amino acid at each position were also comparatively important to predict phosphorylation S/T/Y sites. In fact, the physicochemical properties of amino acids surrounding phosphorylation sites had an impact on phosphorylation event [94]. In contrast, the ratios of DAAC feature vectors are relatively low. Meanwhile, from Figure 3, the constitution of final feature vectors is various in different species models. It shows that the impact of each feature on all species is distinct and then the performance evaluation of the prediction model is different as well. Furthermore, the feature optimization results of serine and threonine phosphorylation in different species also revealed that two-step feature optimization method, which fully considered the importance and contribution of each dimension feature vector, can achieve the higher prediction accuracy.

Feature analysis of phosphorylation of different species

After that, we further detected and visualized statistically significant differences about the compositional biases of sequences between phosphorylation and non-phosphorylation sites. Figure 4a displays the frequency information of amino acid occurrence in serine phosphorylation. For Aspergillus, aliphatic amino acids such as aspartic, proline, alanine, glutamic and glycine are enriched in the positive fragment, while polar amino acids including cysteine, threonine and serine are depleted. For S. pombe, lysine, arginine, aspartic and glutamic, which have charged, are enriched in the positive fragment, while aliphatic amino acids such as leucine, asparagine and tryptophane are depleted in the negative fragment. It shows that there have significant differences among different species. This result is certified by the two-sample logo results (Figure 4c–d) [95]. Similarly, from Figure S1, we also find that there have significant differences in threonine phosphorylation for seven species. Then, based on IG, we chose the rank top 9 from 544 PSP to construct our model (Figure 4b). From the analysis of feature importance, we find that the selected nine PSP features make a great contribution to our model indeed.

Figure 5a and b reflects the frequency information of amino acid pair occurrence. Taking serine phosphorylation in Aspergillus as an example, after the two-step feature optimization, we extracted the DAAC feature from the final optimized feature set to plot the radar map (Figure 5a). We can find that the SP and SS have the highest values in 15 amino acid pairs, which represent the two most important types of amino acid pair. Figure 5b shows the frequency information of 400 amino acid pairs for serine phosphorylation in Aspergillus. It can be seen that SP and SS also have higher frequencies than other amino acid pairs. Meanwhile, the motif scores of Figure S2 are in accordance with the above results [96]. Similarly, we make the same analysis for threonine phosphorylation in C. neoformans (Figure S3). From Figure S3, the frequency of TP is far higher than other 26 amino acid pairs and the heat maps of DAAC also display that the amino acid pair TP is significant, which are consistent with motif result. Subsequently, in Figure 5c and Figure S4a, we further analysed the posterior probabilities of each amino acid at each position in the positive and negative training datasets for serine phosphorylation in Aspergillus, respectively. From Figure 5c, we observe that positive samples tend to enrich at position −3, +1 and +3, which is identical with the two-sample logo results, while the negative samples have almost the same value at each position, indicating that there have significant differences between positive and negative samples. It can be seen from the previous analysis of feature importance that the proportion of the selected BPB feature to the total reaches 80%. Similarly, we observe that the average frequency value of positive samples at position −4, −3, −2, +2, +3 and +4 is relatively higher than negative samples in Figure S4a. Finally, we compared the average KNN score of phosphorylation sites with those of non-phosphorylation sites in Figure 5d and Figure S4b. Due to the gap of data, we selected different k values and comparison sets (the detail in Table S2). With the biggest difference between positive and negative, KNN score with different k values exhibited the best performance. For example, for Aspergillus, the average KNN score of serine phosphorylation is within 0.70−0.88, while the non-phosphorylation average KNN scores are within 0.46–0.53. From the above results of feature optimization, the Acc value of Aspergillus is 0.812, which have the higher accuracy than that of other modifications. Taken together, it shows that there have significant differences among different species in the same modification residue, which indicates the necessity of species-special identification of fungi phosphorylation sites.

Species-specific fungi phosphorylation site prediction of PreSSFP

In this work, 15 models for fungi phosphorylation sites prediction were constructed. On the basis of 10-fold cross-validation, we evaluated the performance of each model on each training dataset and the ROC curve is shown in Figure 6. For serine phosphorylation in seven species, the AUC values of Aspergillus, C. neoformans, F. graminearum, M. oryzae, N. crassa, S. cerevisiae and S. pombe were 0.807, 0.808, 0.802, 0.808, 0.807, 0.784 and 0.775, respectively. For threonine phosphorylation in seven species, the AUC values were 0.754, 0.811, 0.812, 0.824, 0.788, 0.751 and 0.734. For tyrosine phosphorylation in S. cerevisiae, the AUC value was 0.70. From the ROC curves, we knew that those models have relatively good confident predictions. In contrast, we used the tyrosine phosphorylation dataset of each species except S. cerevisiae as testing set to test tyrosine phosphorylation site model for S. cerevisiae, but the results of prediction was not fully satisfactory and achieved lower AUC scores and specificities (the detail in Table S3). In general, these justify the necessity of species-specific prediction and our prediction model PreSSFP is a robust predictor.

Comparison with other existing tools

For the purpose of evaluating the prediction performance of the PreSSFP objectively, general and species-specific predictors were selected to compare with PreSSFP. Taking into account the availability and representativeness of tools, general tool NetPhos 3.1, yeast-specific predictor NetPhosYeast, soybean-specific predictor PHOSFER and human-specific predictor iPhos-PseEn were employed. Because NetPhosYeast tools can only predict S/T phosphorylation in yeast proteins, we used it to predict S. cerevisiae, S. pombe and C. neoformans data. The tool NetPhos 3.1, PHOSFER and iPhos-PseEn were used to predict all seven species. We directly submitted testing data set of each species for the prediction, and the results of PreSSFP were used for a comparison. As shown in Table 4, for the serine phosphorylation prediction of S. pombe, the AUC value of PreSSFP is 0.784. Compared with NetPhos 3.1, PHOSFER and iPhos-PseEn, the AUC value of our model has improved by 19.3%, 13.69% and 14.27%, respectively. Meanwhile, the AUC value of NetPhosYeast is 0.676. Comparison of prediction performance of PreSSFP with other tools in threonine and tyrosine phosphorylation sites is listed in Table S4. For the threonine phosphorylation in C. neoformans, compared with NetPhos 3.1, PHOSFER and iPhos-PseEn, the AUC value of our model has increased at a rate of 34.4%, 16.33% and 10.38%, respectively; compared with the test results of NetPhosYeast, the AUC value of our model has improved by 24.18%. This suggests that PreSSFP is better than general tools, and there are some differences among serine, threonine and tyrosine phosphorylation sites in the same species.

The following reasons may account for such a big difference. First, with the large-scale identification of fungi phosphorylation sites using MS, the number of fungi phosphorylation sites has increased apparently. The data may cause a bias classification of the training model. For example, the NetPhosYeast only considered a total of 953 phosphoserine sites and 192 phosphothreonine sites from 675 yeast proteins. However, PreSSFP integrated all experimental fungi phosphorylation sites from FPD. Second, different organisms may significantly differ in sequences or structural patterns around the phosphorylation sites. For instance, comparative analysis of serine phosphorylation event between different species of eukaryotes performed by Frades et al. revealed that the animal-specific motifs are mainly basic amino acids and the plant-specific top discriminative n-grams contain many acidic amino acids [93]. The general tool NetPhos 3.1, yeast-specific predictor NetPhosYeast, soybean-specific predictor PHOSFER and human-specific predictor iPhos-PseEn are not suitable for predicting all species. The comparison results of PreSSFP with other tools just highlighted the necessity for developing systematic species-specific model to improve the prediction performance of phosphorylation sites. Third, the NetPhos 3.1 and NetPhosYeast only concentrated on the sequence information, and PHOSFER merely combined several amino acid indices to encode protein sequences and iPhos-PseEn just considered the evolutionary information, which could not extract fully features information. Nevertheless, PreSSFP integrated five feature strategies to ensure the complete extraction of sequence-derived information, evolutionary information and physical chemistry properties information.

Conclusions

In this study, on the basis of the primary protein sequences, an online tool named PreSSFP that used a two-step feature optimization method was developed for identifying potential fungi phosphorylation sites. The corresponding analyses and comparison with the existing tools demonstrated that PreSSFP was stabilized and satisfied in the prediction performance and considerably improved the prediction results of fungi phosphorylation sites. Feature analysis showed that the different fungi-species phosphorylation sites have some significant differences in sequenced-derived information, indicating that distinguishing different species was important for fungi phosphorylation sites prediction in order to improve prediction quality. Meanwhile, feature optimization exhibited that KNN feature was significant and exerted a great influence on predicting fungi phosphorylation sites. For the future prediction of fungi phosphorylation sites, newly identified fungi phosphorylation sites will be continuously collected and integrated into computational models, for a better prediction if available. In conclusion, we anticipate that PreSSFP can exert a complementary effect on existing approaches in fungi phosphorylation sites identification. Additionally, the combination of computational analyses with experimental verification could afford useful information for the understanding of the modified mechanism.

Funding

This work was supported by the National Natural Science Foundation of China (21665016 and 21305062) and the Natural Science Foundation of Jiangxi Province (20151BAB203022).

Cao Man is a graduate student at School of Sciences, Nanchang University. Her research focuses on integrative tyrosine post-translational modifications data for analysis and validation.

Guodong Chen is a graduate student at School of Sciences, Nanchang University. His current research focuses on developing novel data analysis algorithms and software of prokaryotes lysine acetylation.

Jialin Yu is a graduate student at School of Sciences, Nanchang University. His research focuses on deep learning and prediction of protein structure and function.

Shaoping Shi is an associate professor at School of Sciences, Nanchang University. Her research focuses on the development of novel data analysis algorithms and bioinformatics tools for prediction of protein structure and function.

References