Protein phosphorylation is one of the most studied posttranslational modifications (PTMs) that plays essential roles in many vital biological processes like cell metabolism, cell motility, apoptosis, replication, transcription, environmental stress responses, DNA repair, immunological responsiveness, and cell cycle control. [1][2][3][4][5][6][7] As much as 30% of all eukaryotic proteins are phosphorylated and disruption in the pathway of phosphorylation is associated with the pathological progression of diseases such as Parkinson's, Alzheimer's, cancer, and heart disease, 6,[8][9][10] phosphorylation is commonly found on serine (S), threonine (T), tyrosine (Y), and histidine (H) residues of proteins. 11 Experimental techniques like low throughput 32 P-labeling 12,13 and high throughput mass spectrometry 14,15 are used to detect phosphorylation sites in proteins. However, these techniques are time-consuming and labor-intensive. In this regard, several computational approaches for the prediction of the general phosphorylation sites have been developed. NetPhos 16 is an artificial neural network (ANN)-based approach for prediction of phosphorylation sites. RF-Phos 17 is a random forest-based approach for prediction of general phosphorylation sites.
Similarly, PhosPred-RF 18 is another phosphorylation site predictor based on random forest.
Recently, deep learning (DL) architectures have been used to predict various PTMs in proteins. Unlike machine learning (ML)-based models, the DL architecture does not require manual feature extraction. For example, MusiteDeep 19 is a DLbased approach that utilizes one-hot encoding and convolutional neural network (CNN) with an attention layer to predict phosphorylation sites. Likewise, DeepPhos 20 uses densely connected CNN (DCCNN) blocks with different filter sizes and windows to learn multiple representations of sequences to predict phosphorylation sites. Furthermore, Wang and Xu devised a capsule network (CapsNet) 21 -based architecture for prediction of protein phosphorylation sites.
Furthermore, DeepPSP 22 merges local window (51) and global window (2000) information using squeeze-and-excitation blocks and LSTM blocks to further improve the phosphorylation PTM site prediction. Among the existing approaches, DeepPSP is quite unique in the sense that this approach integrates both local and global information for predicting the phosphorylation sites. Note that these deep learning methods utilize a one-hot encoding scheme for representation of amino acids. Chlamy-EnPhosSite 23 is an organism-specific phosphorylation site predictor forChlamydomonas reinhardtii based on an ensemble approach that combines long short-term memory and convolutional neural network models in which a supervised word embedding scheme is used to encode the protein sequences. Similarly, PhosIDN 24 uses an integrated deep learning architecture to improve phosphorylation prediction by combining local sequence (using one-hot encoding) and protein-protein interaction (from the STRING database) 25 information. Admittedly, considerable progress has been made in the development of general phosphorylation site prediction methods. However, for all current predictors of phosphorylation sites, the input features are either one-hot encoding (each amino acid is represented with a binary vector), embedding encoding (a learned representation where amino acids that are similar have similar representation), or handcrafted features (physiochemical features extracted from protein sequence).
Recently, the transformer-based 26 language model has shown huge potential to unfold the meaningful latent representation of the sentences/language by implementing a multihead selfattention-based mechanism with masking. By considering the protein sequences as sentences, Elnaggar et al. (ProtTrans) 27 developed a pretrained protein language model (pLM), namely, ProtT5-XL-UniRef50 trained on 2.5 billion protein sequences from the UniRef50 database. The representation learned by this model has been used in various prediction tasks, and the results demonstrate that the information on the evolutionary context of a sequence, contact map, taxonomy, long-range dependencies, protein structure, subcellular localization, physicochemical properties, and function is encoded in their distributed representation. [28][29][30][31][32][33][34] Thus, a more effective model of phosphorylation site prediction may be established by using the knowledge distilled by this language model.
In this work, we present a novel general phosphorylation site prediction tool named LMPhosSite that integrates per-residue contextualized embedding (aka embeddings that are obtained based on its context) from the pretrained protein language model (ProtT5) with local word embedding. LMPhosSite utilizes stacked generalization to improve the prediction of phosphorylation sites by combining two different models trained using embeddings obtained from the supervised embedding layer and embeddings from the ProtT5 language model. Initially, these two modules are passed to their corresponding deep learning architectures, and finally, the stacked generalization of these two models is performed using a meta-model. When compared with the other existing phosphorylation site prediction tools, LMPhosSite exhibits an improved prediction performance.
To train and evaluate our models, we used the data set from DeepPSP. 22 The data set in DeepPSP were experimentally identified phosphorylation sites collected from the SWIS-SPROT, 35 dbPTM, 36 phosphoELM, 37 and PhosphoSite-PLUS. 38 Subsequently, the CD-HIT 39 tool was used to remove the homologous sequences using a cutoff threshold of 0.5. Finally, the sequences were randomly divided into the training and test set in a ratio of 9:1. The annotated phosphorylation sites from these sequences were defined as positive sites, and any remaining S, T, and Y sites that were not annotated as phosphorylated within the same protein sequence were defined as negative sites. To extract the local information, positive windows were generated with the annotated phosphorylated sites in the middle and an equal number of amino acids on both sides flanking the phosphorylated sites. Negative windows were generated in a similar manner. When the site of interest was located near the N-or C-termini of the protein, pseudoresidues "-" were added to make the window sizes to be of the same length. Duplicates were removed from both the positive and negative data sets. Additionally, we also experimented with a CD-HIT cutoff of 0.3 on the homologous sequences and the corresponding data set and the results are presented in the Supporting Information (Tables S5-S9). However, for comparison purposes with the existing approaches, we use the data set obtained using a CD-HIT cutoff of 0.5.
As S and T residues can be phosphorylated by the same specific kinase (hydroxyl groups in their side chains), the S and T data sets were combined. However, as Y residues undergo different enzymatic processes because of pi-electron and the conjugated electrons 40 (phenol ring on the side chain), the Y data set was kept separate. Like many existing approaches, we developed one model for the combined S and T residues and a separate model for the Y residues. Table 1 shows the total number of positive and negative sites in the training and independent test data sets. To avoid overestimation, we made sure that no protein sequences from the independent test data set are present in the training data set as ProtT5 can learn representation for other sites from the same protein. Also, the number of positive ST sites is slightly different than the original DeepPSP data set as some protein sequences have been updated in the UniProt.
In developing a statistical model for the discrimination of protein phosphorylation sites, one pivotal step is numerical encoding through an encoding scheme that assigns a numerical representation to each amino acid that can accurately reflect the intrinsic correlations with the desired targets. 41 In this study, similar in spirit to DeepPSP, we employed both global and local sequence information for the prediction of protein phosphorylation sites. For extracting global information, we used the entire protein sequence that includes the site of interest as input to a pretrained protein language model to extract per-residue contextualized word embedding. Additionally, for encoding the local information, we utilized a supervised embedding layer to obtain word embedding for the window sequence centered around the site of interest in a supervised manner. We described in detail these two types of embeddings below.
The post-translational classification approaches such as Deep-IPs, 42 LMSuccSite, 43 DeepRMethylSite, 44 and pLMSNOSite 45 showed improved prediction performance using supervised word embedding encoding 46 obtained using a supervised embedding layer (SEL). The SEL learns the representation from integer-encoded sequences through back-propagation in a supervised manner. We used the SEL to encode the local window sequence. First, the 20 canonical amino acids and one pseudoresidue "-" were converted into specific integers ranging from 0 to 20. These make the inputs for the embedding layer that lies at the beginning of our DL architecture. Initially, the weights in the embedding layer were randomly initialized, and these weights were updated during the training process. The key arguments in the embedding layer are vocabulary size, output_dim (size of vector space), and input_length (size of input windows). The output from the embedding layer has the dimension of input_length × output_dim. Based on the 10-fold cross-validation results, for our case, we used a vocabulary size of 23, output_dim of 21, and size of input windows to be 51. This embedding was fed into one of the base-models.
Recently, embeddings from protein language models (pLMs) were utilized to predict binding residues, 47 signal peptides, 48 subcellular localization, 49 and protein structural features. 50 In this regard, here, we employ a protein language model (pLM) called ProtT5-XL-UniRef50 (ProtT5) 27 to obtain contextualized per-residue embedding for the site of interest. ProtT5 is based on the T5 architecture 51 and trained solely on unlabeled protein sequences from BFD (Big Fantastic Database; 2.5 billion sequences including meta-genomic sequences) 52 and Uni-Ref50. 53 To obtain the contextualized per-residue embedding for the site of interest from ProtT5, the overall protein sequence containing the site of interest was fed into a pretrained ProtT5 27 model and the fixed-length per-residue features were extracted from the last encoder layer. This contextualized (1024 length feature vector) embedding of the site of interest produced by ProtT5 was then fed into the other base-model.
LMPhosSite integrates two base-models: a model for learning encoding from local sequence information obtained from the window sequence and another model for learning contextualized embedding from the global sequence information obtained using the overall protein sequence, using stacked generalization. Initially, these base-models are trained independently. Finally, score-level fusion of these two models is performed by using a meta-classifier based on a fully connected network. The proposed architecture of LMPhosSite is shown in Figure 1.
The input to this model is the (local) window sequence with the site of interest in the middle. For the embedding layer, a window sequence of size 51 (as shown in Figure 2) was taken around the site of interest, and features were extracted subsequently for classification based on this window. The architecture of this model starts with a supervised embedding layer followed by two layers of a one-dimensional convolution neural network (1D-CNN) and a fully connected layer. The first and second layer 1D-CNN were configured to 128 and 64 filters, respectively, and the size of the filter was set to 3. Additionally, the 1D-CNN layer uses the ReLU activation function and a dropout probability of 0.3. All generated feature maps from 1D-CNN were then sent into the MaxPooling layer with a pool size of 2. The MaxPooling layer samples the output feature map of the convolution layers and extracts the most obvious features. Subsequently, the output feature maps from 1D-CNN were flattened and sent to a fully connected layer with 256 neurons that has a ReLU activation function and a dropout probability rate of 0.4. The output of this model is a probability (of being phosphorylated) score between 0 and 1. Note here that the architecture and the associated hyperparameters for this base-model were obtained using 10-fold cross-validation with grid-search. Essentially, several DL architectures such as 1D-CNN, 2D-CNN, LSTM, and BiLSTM were implemented for this model, and 1D-CNN was selected based on the 10-fold cross-validation results.
The input to this model is the contextualized embedding of the site of interest (size = 1024) obtained from ProtT5 extracted using the overall sequence. Initially, we extracted a static embedding of size n × 1024, where n is the size of the overall sequence, from the last hidden layer of the ProtT5 architecture. However, we only used embedding of the "site of interest" for our purpose. For instance, given a protein sequence of length of 100 and the site of interest located at the 25 th position in the sequence, the ProtT5 will output a contextualized feature vector of size of 1 × 1024 for every 100 amino acid positions (totaling 100 × 1024 size feature matrix). We only utilized the feature vector (dimension: 1 × 1024) corresponding to the site of interest, i.e., 25 th amino acid. In our case, the sites of interest are S, T, and Y where they can either belong to the positive set or a negative set. We presented a schematic representation of the overall architecture and how features are extracted for the site of interest using ProtT5 in the Figure 1. The input sequence of length n with the site of interest ("T", denoted in red) was fed to ProtT5. From these embeddings, we only used the 1024 length feature vector of the site of interest ("T") for classification purposes.
The architecture of this model consists of an ANN with two hidden layers and one output layer. The first hidden layer has 128 neurons, and the second hidden layer has 32 neurons. The output of this model is also a probability (of being phosphorylated) score between 0 and 1. Similarly, various architectures such as RF, SVM, XGBoost, logistic regression, and ANN were implemented for this model, and an ANN-based architecture was selected based on 10-fold cross-validation results.
Next, stacked generalization of the two base-models was performed to obtain a meta-classifier using score-level fusion. The input to this classifier is the two probability outputs from the base-models: the supervised embedding layer model and the ProtT5-based model. The architecture for this model used a basic ANN with one hidden layer having two neurons. All the hyperparameters of individual base-models and the meta classifier used in LMPhosSite are presented Table 2.
Note here that the architectures and the associated hyperparameters for the meta-classifier were also obtained using 10fold cross-validation with grid-search.
In this study, phosphorylated sites are considered as positives sites, and nonphosphorylated sites are considered as negatives sites. Additionally, phosphorylated sites and nonphosphorylated sites predicted correctly by the model are true positive (TP) and true negative (TN), respectively. The negative sites misclassified as positive sites are false positive (FP) and positive sites misclassified as negative sites are false negative (FN). Ten-fold cross-validation was utilized to evaluate the performance of the models. For the results of 10-fold cross-validation, unless otherwise noted, all performance metrics are reported as the mean value ± one standard deviation. Four metrics including precision (Pre), recall (Rec), F1-score, and Matthew's correlation coefficient (MCC) were used to evaluate the performance of the models (with a probability decision threshold = 0.5). Furthermore, the area under the receiver operating characteristics (ROC) curve and area under the precision-recall curve (PrAUC) were also used as performance metrics. Equation 1 describes precision, recall, F1-score, and MCC where TP, TN, FP, and FN represent true positive, true negative, false positive, and false negative, respectively. a Dense layers represent the fully connected layers in TensorFlow.
b The size of convolution layers means the kernel sizes, and the size of dense layers denotes the number of neurons in hidden states.
■
As discussed in the Materials and Methods section, LMPhosSite integrates two types of embeddings using stacked generalization by performing score-level fusion. Using 10-fold cross-validation, we first find the best architectures and hyperparameters for the individual models and subsequently find the best hyperparameters for the meta-classifier using the 10-fold crossvalidation strategy. Essentially, we analyzed the comparative performance of various architectures for base-models and the meta-classifier using 10-fold cross-validation techniques. Finally, we compare the performance of LMPhosSite against existing phosphorylation tools using independent testing. Below, we discuss the results in detail.
The recent studies [54][55][56] unravel that neighboring residues can influence the phosphorylation status of the site of interest (in this case serine, threonine, and tyrosine residues). In that regard, to capture local sequence information for the site of interest, a local window sequence centered around the site of interest (S, T, or Y) surrounded by an equal number of flanking residues on both sides is generally taken as input. However, the window size is also a parameter. To determine the window size, we performed a 10-fold cross-validation for the supervised embedding layer model. Essentially, we experimented with various window sizes ranging from 15 to 61 in an increment of 2 for combined ST and Y phosphorylation training data sets (Figure 2) and 3 showed the mean Matthew's correlation coefficient (MCC) produced using different window sizes on the combined ST and Y phosphorylation training data set using 10-fold cross-validation. The detailed results of the analysis are presented in Tables S1 andS2. Further window sizes were not analyzed due to the sheer size of the windows and the corresponding increase in the number of pseudoresidues "-" that were required at higher window sizes (for residues near Nand C-termini).
From Figure 2, it can be observed that the best performing window size was 51 (highest MCC, 0.434 for ST and highest MCC, 0.3041 for Y) for the combined ST and Y phosphorylation training data set. Hence, 51 was selected as the value of window size for the combined ST and Y phosphorylation residue for subsequent analysis. Interestingly, the DeepPSP 22 phosphorylation prediction method also uses the same window size (= 51).
Subsequently, the model architectures for the two basemodels and the meta-model are selected based on 10-fold crossvalidation results.
Here, we discuss the comparative performance of the basemodels and the meta-classifier using 10-fold cross-validation on the ST phosphorylation training data set. The results are presented in Table 3. The performance metrics used are precision, recall, MCC, and F1-score.
As mentioned above, the training set consists of 165,483 phosphorylated sites and 165,483 (under sampled from 878,133) nonphosphorylated sites. Consequently, the 10-fold cross-validation was performed on the combined ST phosphorylation training data set (330,966 training examples) to choose the best performing meta-classifier (PrT5 + Em). We also performed experiments using cost-sensitive learning by using all the negative sets, and the result of this analysis is shown in Table S4. The mean MCC of the supervised word embedding basemodel (Em) is 0.434 ± 0.004, and the mean MCC of the ProtT5 base-model (ProtT5) is 0.461 ± 0.005. Additionally, it can be observed from Table 3 that the meta-classifier produces a mean MCC, mean precision, mean recall, and mean F1-score of 0.502 ± 0.004, 0.766 ± 0.006, 0.721 ± 0.007, and 0.743 ± 0.002, respectively. Thus, the MCC of the meta-classifier is higher than the MCC of both base-models for the prediction of phosphorylation site prediction, indicating that meta-classifier performs better than the base-models for phosphorylation prediction. Hence, we select this meta-model that combines the two base-models as the final model for the prediction of phosphorylation sites and we call this model as LMPhosSite. As described earlier, the architecture of the meta-classifier is a basic ANN with one hidden layer having two neurons. In conclusion, the meta-model (PrT5 + Em) that combines the local sequence information from the supervised embedding layer (embedding (Em)) and the global information from the ProtT5 model (ProtT5 (PrT5)) produces the best MCC. Note that the architecture for the base-models was independently optimized using 10-fold cross-validation.
Subsequently, we assessed the performance of LMPhosSite and the base-models on the independent test set for ST and Y sites. The results of LMPhosSite on independent ST sites are presented in Table 4. As seen from the table, LMPhosSite produces a precision, recall, Matthew's correlation coefficient (MCC), and F1-score of 38.78%, 67.12%, 0.390, and 49.15%, respectively. The LMPhosSite was able to classify 82,275 samples as true negative, 12,452 samples as true positive, 19,650 samples as false positive, and 6,099 as false negative for the ST independent data set. Similarly, when we tested the LMPhosSite on Y phosphorylation sites and tested on the Y phosphorylation test data set, it achieved a precision, recall, Matthew's correlation coefficient (MCC), and F1-score of 34.90%, 62.03%, 0.298, and 44.67%, respectively. The LMPhosSite was able to classify 10,735 samples as true negative, 2,010 samples as true positive, 3,748 samples as false positive, and 1,230 as false negative for the Y independent data set. The independent test results for basemodels are shown for informative purposes only, and the final model (meta-classifier) is selected based on 10-fold crossvalidation.
Figure 4 shows the receiver operating characteristic (ROC) curve for base-models and LMPhosSite on the combined ST independent test data set. It can be observed that LMPhosSite (PrT5 + Em) has the highest area under the curve (0.823) compared to the base-models. Additionally, we also plotted the precision-recall curve for the base-models and LMPhosSite on the combined ST independent test data set, and the results are shown in Figure 5. It can be observed that LMPhosSite also has the highest precision-recall area under the curve (PrAUC = 0.532).
Moreover, we also plotted ROC and precision-recall curve of the phosphorylation Y independent data set, and the corresponding ROC and precision-recall curves are shown in Figures 6 and7, respectively. The AUC and PrAUC of LMPhosSite are the highest compared to the base-models for the Y phosphorylation independent data set too.
We further compared the performance of LMPhosSite with some publicly available predictors including DeepPSP, 22 MusiteDeep, 19 Musite, 57 and CapsNet, 21 5 that LMPhosSite trained on a combined ST phosphorylation data set produced an MCC, precision, recall, and F1-score of 0.3905, 38.78%, 67.12% and 49.15%, respectively.Moreover, LMPhosSite trained on the Y phosphorylation data set produced an MCC, precision, recall, and F1-score of 0.2984, 34.90%, 62.03%, and 44.66%, respectively. LMPhosSite produces the highest MCC (0.3905 for ST and 0.2984 for Y) among the compared approaches for both ST and Y phosphorylation independent testing. Moreover, the confusion matrices of all the predictors are listed in Table S3.
These results demonstrate that LMPhosSite performs better than the compared approaches across all these performance metrics except for recall for combined ST and produces the best results across all these performance metrics except for recall for the Y phosphorylation data set.
Additionally, to discern the classification effectiveness of the features from the base-models as well as the features from the stacked generalization model, we used t-distributed stochastic neighbor embedding (t-SNE) 58 to project these features into the two-dimensional space (Figure 8). The main purpose of t-SNE visualization in this work is to visually observe the separation boundaries between classes in two-dimensional space. The t-SNE plot was generated using 50 as the value for the learning rate to visualize the learned ProtT5 features and features obtained from the penultimate dense layer of the trained model. Note that these plots were generated using 8,274 (4,137 positive and 4,137 negative) randomly chosen data points.
For the features extracted from the language model ProtT5 of phosphorylated or nonphosphorylated token "S/T", there are some clusters for positive and negative samples; overall, the positive and negative samples are still not distinct (Figure 8A). Figure 8D represents the t-SNE plot of the feature vectors generated from the concatenation layer of the trained deeplearning architecture, which shows that negative samples (blue points) are concentrated at the third quadrant, while positive samples (orange points) are concentrated at the first quadrant, which indicates that the ensembled DL architecture that processes features from two information sources with MLP as a meta-classifier largely clusters positive and negative samples in 2 space. For further elaboration, we also provide the t-SNE plot generated from the penultimate layer of trained base-models ProtT5 and embedding encoding in Figure 8B and Figure 8C, respectively.
Protein phosphorylation is one of the most important posttranslational modifications. In this work, we developed a stacked generalization approach called LMPhosSite to predict protein phosphorylation sites. LMPhosSite uses a supervised embedding layer to encode local sequence information and a pretrained protein language model (ProtT5) to encode global sequence information and then integrates these two types of information using a stacked generalization approach. The novelty of the approach is the use of global information extracted from contextualized embedding obtained from a pretrained protein language model (ProtT5) for the site in conjunction with supervised word embedding obtained from the local sequence window. The independent test results show that LMPhosSite achieves better performance than the compared approaches. Hence, it can be concluded that LMPhosSite is a promising general phosphorylation site prediction tool.
The improved performance of the LMPhosSite can likely be attributed to two things: the use of contextual protein language models to extract features from the overall protein sequence for the site of interest and the novel architecture that uses stacked generalization. To our knowledge, LMPhosSite is the first approach that utilizes the distilled information from large pretrained protein language models for the prediction of phoshorylation sites. Similar in spirit to DeepPSP, our model combines the global information (the ProtT5-based embeddings obtained for the site from the whole protein sequence) with the local information (the supervised word embedding obtained from the window sequence using the supervised embedding layer). Since these pLM embeddings can be easily extracted for any protein sequence, LMPhosSite provides robust and fast predictions for phosphorylation sites.
The improvement in the performance of these approaches that utilize language models may be improved by (i) fine-tuning a pretrained protein language model using the proteins of the same characteristics (in this case, phosphorylated proteins), (ii) training newer protein language models using newer language models like GPT-4, and (iii) combining language model-based features with other physio-chemical features. Finally, with the development of newer and more powerful protein language models, the prediction performance of the approaches that make use of distilled information from these language models is likely to improve.
https://doi.org/10.1021/acs.jproteome.2c00667 J. Proteome Res. 2023, 22, 2548-2557