PhosVarDeep: deep-learning based prediction of phospho-variants using sequence information

Dataset and pre-process

Here, in order to train and evaluate our method, we collect over 2,000 potential phospho-variants in PhosSNP (Ren et al., 2010) as the positive set and adopt three negative sets generated by Patrick et al. (2017), which contain information of protein UniProt index, amino acid variation and phosphorylation sites. The positive dataset in PhosSNP we adopt is based on exact string matching reference/variant sequence with experimentally identified human phosphorylation sites. Accordingly, either identical hit in reference or variant sequence is considered as potential phospho-variants (Ren et al., 2010). As for negative sets, as described in Ralph et al., the three high confidence negative sets adopted in this study are generated by different criteria under the fact that phosphorylation sites are less likely to occur (1) in solvent inaccessible/buried regions of a protein; (2) in transmembrane domains; (3) in proteins that do not interact either directly or through mediators with a kinase (Patrick et al., 2017). Next, we map the collected phospho-variant data with protein sequences by its corresponding UniProt index and amino acid variation so that each phospho-variant corresponds to a pair of reference and variant protein sequences (Bateman et al., 2015). Table 1 shows the sizes of the positive set and three negative sets used for phospho-variant prediction, with respect to different phosphorylation site types, i.e., serine (S)/threonine (T) or tyrosine (Y). Given a phospho-variant, we intercept the protein fragments centered at a phosphorylation site on its corresponding pair of reference and variant sequences. Then each protein fragment is encoded by the one-hot encoding strategy into a L × N two-dimension matrix, where L represents the window size of the protein fragment, and N is set to 21 according to the total number of common amino acids (Min et al., 2017). Besides, we apply CD-HIT with similarity threshold of 40% to all the collected data to reduce the sequence redundancy by following previous studies (Pan et al., 2014; Zhao et al., 2012).

After completing the data pre-processing, the positive set is combined with each negative set to form three phospho-variant datasets, and for each phospho-variant dataset, we train two deep learning models using phospho-variant data on S/T and Y sites respectively. Meanwhile, we use a performance evaluation strategy commonly adopted in deep learning methods for sequential analysis (Zhou & Troyanskaya, 2015; Khurana et al., 2018), that is, each dataset is randomly divided into strictly non-overlapping training, validation and test sets and the ratio is set to 6:2:2 in our study. In this way, we adopt the training data to tune the model weights and utilize the validation data to prevent overfitting (Zhou & Troyanskaya, 2015). The test set is adopted to evaluate the performance of PhosVarDeep and to implement the comparison of PhosVarDeep and other phospho-variant prediction approaches.

Siamese-like architecture of PhosVarDeep

In this study, we design a Siamese-like deep-learning framework for phospho-variant prediction. As a classic metric learning method, Siamese neural network is first introduced in Bromley et al. (1993) and has been successfully adopted in tasks such as signature verification (Chopra, Hadsell & LeCun, 2005), face verification (Cao, Ying & Li, 2013) and object tracking (Bertinetto et al., 2016). To measure similarity between pairs of inputs, Siamese neural network learns a function mapping input patterns to the representation or target space where similar pairs will get close and non-similar pairs will get away from each other (Zagoruyko & Komodakis, 2015). Generally, Siamese neural network contains two subnetworks sharing same configuration, weights, and parameters, which ensures that two similar inputs are transformed into similar feature representations. Here, we adopt a deep-learning framework in Siamese style to learn a complex non-liner mapping for separating the positive and negative samples of the phospho-variants. Figure 1 shows the proposed architecture of PhosVarDeep, which consists of two identical subnetworks with shared weights and a prediction module. Specifically, within each subnetwork, the PhosFEN and the CNN module extract a high dimensional feature representation of an input sequence. Next, the prediction module is introduced to integrate the outputs of the two subnetworks as combined features to generate prediction results of phospho-variants. The details of those networks are described as follows.

PhosFEN

In each subnetwork, PhosFEN presents a sequence feature encoding network to extract the respective features from the reference and variant local sequences of a given phospho-variant. Specifically, the input of PhosFEN is one-hot local sequence e_s ∈ R^L×21 (s = s₁, s₂), with s₁ and s₂being the reference and variant sequences, respectively. Then the sequence features related to phosphorylation are extracted by PhosFEN as follows: (1)${\displaystyle f_s=\text{PhosFEN}\left(e_s;W^F\right)}$ where W^F represents all the parameter matrices and bias items in PhosFEN. In this way, the feature map f_s is exported from PhosFEN as general phosphorylation-related local features, and the pair of f_s1 and f_s2 is used as input to the following CNN module.

In this work, since the amount of training data of phospho-variants is much smaller than the size of parameters in deep learning models, overfitting is likely to occur in the training process. In order to solve the problem of insufficient training data, we have studied the feature learning capabilities of the existing deep-learning models for phosphorylation site prediction and adopt DeepPhos (Luo et al., 2019) in our study through transfer learning (Yosinski et al., 2014). DeepPhos is a pre-trained deep learning model consisting of densely connected CNN blocks, which can capture the complex deep sequence representations relevant to phosphorylation better than other deep-learning models. We transfer the whole layers of the model to PhosFEN except its final three layers including the flatten layer, the fully connected layer, and the output layer in our study.

CNN module

After the general phosphorylation sequence features are extracted from PhosFEN, they are fed into the CNN module to generate variant-aware phosphorylation sequence features. As shown in Fig. 1, each CNN module is comprised of multiple convolutional layers and max pooling layers, and a flatten layer. The convolutional layer generates a feature map by convoluting the input with a set of convolution kernels. Mathematically, for input phosphorylation-related features f_s (s = s₁, s₂), let h_i,s be the ith feature map in output, and the convolutional layers can be described as follows:

(2)${\displaystyle h_{1,s}=\alpha \left(W_1^k{f}_s+b_1^k\right)}$ (3)${\displaystyle h_{i,s}=\alpha \left(W_i^k{h}_{i-1,s}+b_i^k\right)2\le i\le M}$ where $W_i^k$ and $b_i^k$ refer to the parameter matrices and bias items in the ith convolutional layer, with k being the number of convolutional kernels. α represents ReLU activation function that can realize the nonlinear transformation, M refers to the number of convolutional layers and here is set to 3. In this way, the variant-aware phosphorylation sequence features can be generated, and then by the flatten layer, they are transformed into a pair of one-dimensional tensors (h_s1, h_s2) ∈ R^d. In addition, to relieve the risk of overfitting during the training process, a dropout layer is added after each convolutional layer to discard some neurons randomly (Srivastava et al., 2014). The details of the CNN Module are listed in Table 2.

Prediction module

Taking a pair of sequence features h_s1 and h_s2 as inputs, the prediction module utilizes a multi-layer DNN consisting of three fully connected layers to integrate the feature pair and compute a prediction score for phospho-variant. In detail, the two input features are concatenated and then fed into the multi-layer DNN to capture the abstract combined features c_v ∈ R^u, here u refers to the number of neurons in the final fully connected layer. Next, the output layer with softmax as the activation function is used to generate prediction scores of the positive and negative phospho-variant, which can be calculated as follows: (4)${\displaystyle P\left(y=1|\left(s_1,s_2\right)\right)=\frac{1}{1+exp\left(-c_v{W}_v\right)}}$ (5)${\displaystyle P\left(y=0|\left(s_1,s_2\right)\right)=1-P\left(y=1|\left(s_1,s_2\right)\right)}$ where W_v ∈ R^u×2 represents the weight matrix of softmax function. The details of the multi-layer DNN and output layer are listed in Table 2.

Training

PhosVarDeep is a unified deep learning framework for phospho-variant prediction and is trained to classify the phospho-variants into two classes: positive phospho-variants and negative phospho-variants. Accordingly, for the aim of minimizing training error, here we adopt a binary cross-entropy as loss function: (6)${\displaystyle L_c=-\frac{1}{N}\sum _{j=1}^n{y}^{j}lnP\left(y^j=1|{\left(s_1,s_2\right)}^j\right)+\left(1-y^j\right)lnP\left(y^j=0|{\left(s_1,s_2\right)}^j\right)}$ where N refers to the size of training data, ${\left(s_1,s_2\right)}^j$ represents the pair of reference and variant sequences for the j th input phospho-variant and y^j refers to the corresponding class label. We freeze all the layers of PhosFEN and train the CNN module and prediction module jointly on phospho-variant training data, in which the weights and biases of conventional layers and fully connected layers are the parameters to be estimated. Besides, as a widely used stochastic gradient descent algorithm, Adam optimizer (Kingma & Ba, 2015). is adopted in the training process. At the same time, we use mini-batch training strategy in this study to randomly divide small proportions of the training samples in each iteration into optimizer loops. In addition, to deal with the problem of data imbalance, we follow the previous study (Wang et al., 2017) to apply a bootstrapping strategy in our deep learning method.

Performance assessment

In this study, several widely used measurements are leveraged to assess the performance of our proposed PhosVarDeep, which include specificity (Sp), sensitivity (Sn), precision (Pre), overall accuracy (Acc), F1 scores and Matthew’s correlation coefficient (MCC). Their detailed definitions are as below:

(7)${\displaystyle Sp=\frac{TN}{TN+FP}}$ (8)${\displaystyle Sn=\frac{TP}{TP+FN}}$ (9)${\displaystyle Pre=\frac{TP}{TP+FP}}$ (10)${\displaystyle Acc=\frac{TP+TN}{TP+TN+FP+FN}}$ (11)${\displaystyle F1=\frac{2\times Pre\times Sn}{Pre+Sn}}$ (12)${\displaystyle MCC=\frac{TP\times TN-FP\times FN}{\sqrt{\left(\left(TP+FN\right)\times \left(TP+FP\right)\times \left(TN+FN\right)\times \left(TN+FP\right)\right)}}}$ where TP, FP, TN, and FN denote true positives, false positives, true negatives, and false negatives, respectively. The other measurements are calculated based on TP, FP, TN, and FN. Moreover, we plot the receiver operating characteristic curve (ROC) and calculate the area under ROC curve (AUC) to evaluate the overall performance.

Determining the PhosFEN model

To begin, in order to extract general phosphorylation sequence features, we transfer phosphorylation prediction models to PhosFEN in our method. There are several pre-trained models with good phosphorylation-related sequence feature learning capabilities. Among them, Musitedeep (Wang et al., 2017) and DeepPhos (Luo et al., 2019) are both CNN models achieving significantly better performance than previous methods. To determine which model provides more useful information for phospho-variant prediction, we compare the two phosphorylation site predictors using an existing comprehensive phosphorylation site dataset (Luo et al., 2019) collected from Phospho.ELM, PhosphositePlus, HPRD, dbPTM and SysPTM. The Roc curves and AUC values of the two models are shown in Fig. 2. It is obvious that DeepPhos consistently achieved higher performance on both S/T sites and Y sites than Musitedeep. For S/T sites, the AUC value of DeepPhos is 3.0% higher than that of Musitedeep. For Y sites, compared with Musitedeep, the AUC value obtained by DeepPhos is increased by 5.7%. The above results suggest that DeepPhos has better phosphorylation-related feature learning capability on phospho-variant data. Accordingly, for accurate phospho-variant prediction, we employ pre-trained DeepPhos in PhosFEN to capture general phosphorylation sequence features in phospho-variant data.

Performance evaluation of PhosVarDeep

In this part, to evaluate the performance of PhosVarDeep in extracting and integrating reference and variant sequence information for phospho-variant prediction, we conduct an ablation study by comparing three different model configurations: (1) PhosFEN*: in this case, the paired outputs of PhosFEN are directly combined and fed into a fully connected layer to get the prediction results; (2) PM*: in this case, we employ the prediction module to integrate the outputs of PhosFEN to predict phospho-variant; (3) PhosVarDeep.: in this case, for the paired outputs of PhosFEN, we utilize the CNN module to further learn variant-aware phosphorylation sequence features and employ the prediction module to generate combined features for final prediction.

We apply all above methods on three phospho-variant test sets and Table 3 lists their AUC values on S/T and Y sites. It is obvious that PM* produces higher performance than PhosFEN*, which indicates the prediction module effectively contributes to the performance of phospho-variant prediction. For example, on test set2 the AUC values of PM* are 0.930 and 0.917 on S/T and Y sites respectively, which yield 1.5% and 1.9% improvement over PhosFEN*. More importantly, it can also be clearly observed that PhosVarDeep achieves better performance than PM* on all the test sets. For example, on test set3 the AUC values obtained by PhosVarDeep on S/T and Y sites reach 0.906 and 0.874 respectively, which are 5.8% and 6.2% higher than the corresponding AUC values of PM*. These results demonstrate that the CNN module can effectively learn variant-aware phosphorylation sequence features that are useful for predicting phospho-variant. At the same time, by integrating the CNN module and prediction module, PhosVarDeep consistently outperforms PhosFEN* with remarkable improvements on AUC value. For example, on test set1, the AUC value is enhanced from 0.845 (PhosFEN*) to 0.946 (PhosVarDeep) on S/T sites, and similarly, on Y sites the AUC value is increased from 0.827 to 0.919.

In addition to AUC values, we also utilize Sp, Sn, Acc, Pre, F1 and MCC to verify the effectiveness of the proposed method. Similar to previous studies (Luo et al., 2019; Wang et al., 2021), we calculate the values of other measurements when Sp is equal to high stringency level (95.0%) or medium stringency level (90.0%), respectively. The detailed measurement values on S/T and Y sites are separately displayed in Tables 4 and 5. From these results, we can see that the performance of PM* is better than PhosFEN* across all the test sets. Take test set2 on S/T sites as an example, when Sp is set at 90.0%, the Sn, Acc, MCC, Pre and F1 values of PM* are 7.2%, 3.6%, 7.0%, 0.8% and 4.2% higher than PhosFEN*, respectively. Also, Tables 4 and 5 clearly show the superior performance of PhosVarDeep. For example, on test set1 PhosVarDeep manages to obtain F1 scores of 0.912 at middle stringency level on S/T sites, with an improvement of 8.0% and 4.8% compared with PhosFEN* and PM*, respectively. On test set3 the Acc values of PhosVarDeep at high stringency level are increased by 13.1% and 7.4% on Y sites, respectively. The above results further verify that our proposed PhosVarDeep can effectively extract and combine the information of reference and variant sequences, which leads to a significant improvement in the performance of phospho-variant prediction.

Comparison with existing methods

To further assess the performance of PhosVarDeep for phospho-variant prediction, we compare it against MIMP (Wagih, Reimand & Bader, 2015), a sequence-based method and PhosphoPICK-SNP (Patrick et al., 2017), a motif analysis and context-based method. It is of note that we follow (Patrick et al., 2017) to perform the comparison using those phospho-variants for which we can obtain prediction results from MIMP.

As shown in Table 6, PhosVarDeep obtains the best performance on AUC values in comparison with other methods. Take test set2 as an example, compared with MIMP and PhosphoPICK-SNP, on S/T sites PhosVarDeep obtains 12.6% and 10.6% improvement for AUC values, respectively. Also, on Y sites, the AUC value achieved by PhosVarDeep is 0.922, which has an improvement of 9.5% over the next-best method. Furthermore, we compute the values of Sn, Acc, Pre, MCC and F1 for all the methods on the test sets on S/T sites, and the results are shown in Fig. 3. It can be observed that PhosVarDeep consistently performs better on all the metrics than MIMP and PhosphoPICK-SNP. For example, on test set3 at high stringency level, PhosVarDeep obtains Sn of 0.627, Acc of 0.790, MCC of 0.612, Pre of 0.927 and F1 of 0.748 on S/T sites, while the corresponding values of the next-best method are 0.310, 0.633, 0.341, 0.863 and 0.456, respectively. In addition to performing comparison on the three phospho-variant test sets, we adopt a list of experimentally confirmed phospho-variants compiled by Patrick et al. (2017) and there are 20 examples after deleting the overlapping data in the training set. The scores predicted by PhosVarDeep and MIMP are shown in Table 7. The closer the prediction score is to 1.0, the more likely it is to change the phosphorylation state. From Table 7, we can see that PhosVarDeep remarkably outperforms MIMP for predicting the experimentally confirmed positive examples of phospho-variants. To sum up, these results indicate that PhosVarDeep is a highly competitive and efficient method in predicting phospho-variant based on sequence information.

Visualization of learned features

The combined features extracted by PhosVarDeep and the original combined one-hot encoding features of reference and variant sequences are visualized in this section to intuitively shows the ability of our proposed deep learning method in phospho-variant prediction. Here, we use a popular visualization algorithm t-SNE (Van Der Maaten & Hinton, 2008) and observe the difference between positive and negative examples of phospho-variants. Take test set3 as an example, the original combined sequence features and the abstract combined features extracted by our model are shown in Fig. 4. It suggests that we can hardly separate the positive examples of phospho-variant from negatives by original combined features, while as the extracted combined features show separate trends, we can distinguish these two classes more clearly by the deep representations of PhosVarDeep. These results demonstrate that original reference and variant sequences can be combined and transformed into meaningful representations with stronger discriminant power by PhosVarDeep, which can be helpful for further analysis of phospho-variant prediction.