Prediction of post-translational modification sites using multiple kernel support vector machine

Data collection and preparation

We adopted a dataset of experimentally identified PTMs used in our previous study (Wang, Jiang & Xu, 2015), which included 2,990 S sites and 1,961 T sites (phosphorylation, O-linked glycosylation, acetylation) collected from several major comprehensive PTM databases, including dbPTM (version 3.0) (Lee et al., 2006), PhosphoSitePlus (Hornbeck et al., 2012), Phospho. ELM (Diella et al., 2004), dbOGAP (Wang et al., 2011) and SysPTM (Li et al., 2009). The detailed information of this dataset was provided in Wang, Jiang & Xu (2015). For Y sites, we followed the procedure described in Wang, Jiang & Xu (2015) and collected 1,791 local sequences (containing phosphorylation, sulfation and nitration) from dbPTM database (Lee et al., 2006) and the supplementary material provided by Pan et al. (2014). In these two datasets, the negative data contained the target sites that are not experimentally to be modified for a specific PTM or kinase group. To further confirm the fairness of constructing the negative dataset, we further randomly selected local sequences on S/T/Y sites from the known PTM proteins as additional negative samples, according to the median values of the positive samples of all PTMs or kinase groups, respectively. It should be noted that these additional negative samples were not experimentally to be modified by any PTMs or kinase groups. Finally, we totally obtained 3,239 local sequences on S sites, 2,037 local sequences on T sites and 1,880 local sequences on Y sites for this study. The detailed positive/negative information about each PTM or kinase group was illustrated in Table S1.

Local sequence kernel similarity for target sites

After obtaining protein local sequences (10 upstream and 10 downstream residues and central residue has PTM) on S/T/Y sites, we computed the local sequence similarity for target sites using amino acid substitution matrix BLOSUM62, which has been proven to be efficient for calculating pairwise similarity (Gao et al., 2010). Then, the local sequence similarity between two samples t_i and t_j could be calculated as follows: (1)${\displaystyle S_{seq}\left(t_i,t_j\right)=\sum _{1\le x\le 21}BLOSUM62\left(t_i\left(x\right),t_j\left(x\right)\right)}$ where x is the window size of a local sequence and is set to 21 in this study. $t_i\left(x\right)$ (or $t_j\left(x\right)$) represents the amino acid located in the xth position of t_i (or t_j). Since the similarity between samples should be non-negative, we normalized S_seq using: (2)${\displaystyle K_{seq}\left(t_i,t_j\right)=\frac{S_{seq}\left(t_i,t_j\right)-min\left(S_{seq}\right)}{\max \left(S_{seq}\right)-min\left(S_{seq}\right)}}$ where $max\left(S_{seq}\right)∕min\left(S_{seq}\right)$ represents the largest/smallest number in the matrix, respectively. Thus, after applying this operation to the local sequence similarity, the matrix K_seq was obtained and could be considered as the local sequence kernel similarity, which was both symmetric and positive denite.

Gaussian interaction profile kernel similarity for target sites

As described in Fig. 1A, the actual relationships between target sites and PTMs could be represented as a bipartite network. Formally, given a set $T=\left\{t_1,t_2,\dots ,t_m\right\}$ of sites and $P=\left\{p_1,p_2,\dots ,p_n\right\}$ of PTMs, and an edge $E=\left\{e_{ij},t_i\in T,p_j\in P\right\}$ was drawn in the network if the target site t_i has been experimentally modified by this PTM p_j. For simplicity, we can further characteristic this bipartite network by an adjacency matrix A, in which each row denotes target sites and each column denotes different types of PTM. The entity $A\left\{i,j\right\}$ in row i and column j equals 1 if the target site t_i is modified by the PTM p_j, otherwise 0. Here, m is the number of target sites, and n is the category of the PTM types.

Due to the fact that the row i of adjacency matrix A indicates the interaction profile of a target site t_i (or t_j), which specifies the presence or absence of relationship with every PTM in the constructed bipartite network, we adopted a powerful kernel named Gaussian interaction profile kernel (GIP) that has been widely used in the area of drug-target interaction prediction (Van Laarhoven, Nabuurs & Marchiori, 2011). The definition of Gaussian kernel between target sites was using the follow equation: (3)${\displaystyle K_{GIP}\left(t_i,t_j\right)=exp\left(-\gamma {∥A_{ti}-A_{tj}∥}^2\right)}$ where A_ti (or A_tj) represents the interaction profile for the target site t_i (or t_j), namely the binary vector encoding relationship between sites t_i (or t_j) and each PTM. $∥\cdot ∥$ indicates the Euclidean distance between A_ti and A_tj andthe parameter γ is the kernel bandwidth. Generally, the kernel bandwidth can be obtained by cross-validation process, here in this study was set to 0.001. Finally, the Gaussian interaction profile kernel similarity for target sites, denoted by K_GIP, is an m by m symmetric matrix. It should be noted that K_GIP should be recalculated since the adjacency matrix A changed when performing cross-validation process.

Multiple kernel SVM

In this section first we briefly introduced the concepts of SVM for classification tasks. The detailed information were also provided in Huang & Wang (2006) and Vapnik (2000). Given a training dataset $T=\left\{\left(x_1,y_1\right),\left(x_2,y_2\right),\dots ,\left(x_n,y_n\right)\right\},x_i\in R^m$ and $y_i\in \left\{+1,-1\right\}$. For SVM with L1 soft margin formulation, we can use following equation to deal with the primal problem: (4)${\displaystyle minJ\left(\overrightarrow{\mathrm{w}},\overrightarrow{\xi }\right)=\frac{1}{2}\parallel \overrightarrow{\mathrm{w}}{\parallel }^2+C\sum _{i=1}^n{\xi }_i}$ (5)${\displaystyle \mathrm{s}.\mathrm{t}.y_i\left({\overrightarrow{\mathrm{w}}}^T\varphi \left(\overrightarrow{x_i}\right)+b\right)\ge 1-{\xi }_i,i=1\dots n}$ ${\displaystyle {\xi }_i\ge 0}$ where ξ_i ≥ 0 represent the non-negative slack variables and C is the regularization parameter. The aforementioned quadratic optimization problem could be solved by using the Lagrange function and differentiating with respect to $\overrightarrow{\mathrm{w}},b$ and ξ_i, then the primal problem would transform to the dual problem: (6)${\displaystyle max\sum _{j=1}^n{a}_i-\frac{1}{2}\sum _{i,j=1}^n{a}_i{a}_j{y}_i{y}_{j}K\left(\overrightarrow{x_i},\overrightarrow{x_j}\right)}$ (7)${\displaystyle \mathrm{s}.\mathrm{t}.\sum _{i=1}^n{y}_i{a}_i=0}$ ${\displaystyle 0\le a_i\le C,i=1\dots n.}$ In a classification task, the optimal ${\overrightarrow{a}}^{\mathrm{\ast }},{\overrightarrow{w}}^{\mathrm{\ast }},b^{\mathrm{\ast }}$ would be obtained and the final predictive model can be represented as: (8)${\displaystyle y_i\left(\sum _{j=1}^n{a}_j^{\ast }{y}_{j}K\left(\overrightarrow{x_i},\overrightarrow{x_j}\right)+b^{\ast }\right)=1,i=1\dots n}$ where $K\left(\overrightarrow{x_i},\overrightarrow{x_j}\right)=\varphi \left(\overrightarrow{x_i}\right)\varphi \left(\overrightarrow{x_j}\right)$ is the kernel function.

Here multiple kernels namely local sequence kernel and Gaussian kernel were integrated into the kernel function and was described as follows: (9)${\displaystyle K\left(\overrightarrow{x_i},\overrightarrow{x_j}\right)=\sum _{d=1}^m{\beta }_d{K}_d\left(\overrightarrow{x_i},\overrightarrow{x_j}\right),{\beta }_d\ge 0}$ where m = 2, K₁ and K₂ were local sequence kernel and Gaussian kernel, respectively. In this study, we defined the integrated kernel K as the custom kernel function, instead of using the default kernel of SVM. We used LIBSVM (v.3.17) (Chang & Lin, 2011) SVM implementation freely available for the MATLAB environment. In applying the SVM algorithm to our dataset, we used balanced penalization in the case of positive and negative training dataset of different sizes. In all experiments, we used the default C regularization parameter. The whole procedure of this work was illustrated in Fig. 1.

Performance assessment

In this study, ten-fold cross-validation as described in existing studies (Gao et al., 2010; Xu et al., 2014a; Xue et al., 2006) was applied to assess the predictive performance of the proposed method. For a given PTM, 9/10 randomly chosen samples were used as the training data while the remaining 1/10 were used as the test data. The ten-fold cross-validation tests were repeated 10 times. As a result, the original data set was covered successfully both in the training and in the test data. The final evaluation was based on the average of these ten performances. Receiver-operating characteristic (ROC) curve, which plots true positive rate (sensitivity, Sn) against false positive rate (1-specificity, 1 − Sp) by gradually changing different thresholds, was utilized to estimate the predictive ability of the method. Sn is defined as the proportion of true positives that are correctly observed by the classifier, whereas Sp is given by the proportion of true negatives that are correctly identified. The corresponding area under ROC curve namely AUC is also calculated, with AUC = 1 represents perfect performance and 0.5 means random performance. In addition, other conventional measurements such as precision (Pre), accuracy (Acc) and Matthews’s correlation coefficient (MCC) were also applied to assess the predictive performance, and the definitions were shown as below: (10)${\displaystyle Sn=\frac{TP}{TP+FN}}$ (11)${\displaystyle Sp=\frac{TN}{TN+FP}}$ (12)${\displaystyle Pre=\frac{TP}{TP+FP}}$ (13)${\displaystyle Acc=\frac{TP+TN}{TP+TN+FP+FN}}$ (14)${\displaystyle MCC=\frac{TP\times TN-FP\times FN}{\sqrt{\left(TP+FN\right)\times \left(TP+FP\right)\times \left(TN+FN\right)\times \left(TN+FP\right)}}}$ where TP, TN, FP and FN refer to true positives, true negatives, false positives and false negatives, respectively.

Comparison with existing methods for phosphorylation

To evaluate the power of the proposed method, first three common phosphorylation prediction methods, PPSP (Xue et al., 2006), GPS (version 3.0) (Xue et al., 2008) and NetPhos (version 3.1) (Blom et al., 2004) were used to make comparison. We took kinase groups CAMK, CMGC, CK1 and TKL as examples to illustrate the predictive performance. It should be stated that the ten-fold cross-validation process is not available for GPS and NetPhos, so the phosphorylation dataset was utilized as testing dataset to evaluate the predictive performance, which may lead to over-estimation of the predictive performance of these tools. However, our proposed method still obtained promising and competitive performance. The ROC curves were plotted for four methods to compare the predictive performance at each specificity level and displayed in Fig. 2. As shown in Fig. 2, the proposed method achieved significantly better overall performance for four kinase groups than all other prediction methods. Performance of other kinase groups on S∕T and Y sites were displayed in Figs. S1 and S2, respectively.

Besides the ROC curves, the corresponding AUC value for each phosphorylation kinase group on S/T/Y sites was also calculated for each method and displayed in Table 1. It was indicated that our proposed method was consistently better than GPS, PPSP and NetPhos. For example, the AUC achieved by the proposed method for kinase group CAMK on S sites was 14.7%, 24.2% and 18.2% higher than GPS, PPSP and NetPhos, respectively. For kinase group TKL on T sites, the corresponding AUC values were 88.6%, 71.2% and 81.4% for the proposed method, GPS and PPSP, respectively. Also, from Table 1 it can be seen that our proposed method achieved comparable or better performance than Wang, Jiang & Xu (2015) that also used both sequence information and site-modification relationships, demonstrating the feasibility and usefulness of the developed kernels in predicting PTM sites. In addition, in order to ensure that the redundancy between training and evaluation data was minimized, all protein sequences were grouped into ten sets using BLASTClust by following Dou, Yao & Zhang (2014). Then, the proposed method was compared with PPSP and Wang, Jiang & Xu (2015) by cross validation of these grouped ten sets. For GPS and NetPhos, cross validation process is not supported as they only provide the web servers to make prediction. The results listed in Table S2 indicated that our proposed method was also consistently better than PPSP and in general comparable to Wang, Jiang & Xu (2015). Taken together, the proposed method achieved comparable or better performance for the prediction of phosphorylation sites.

Additionally, by following the study of Fan et al. (2014), a threshold was set for each method such that the specicity of each method was equal to 95.0% (medium) or 99.0% (high). We took two kinase groups (CAMK and CMGC) as examples, and the corresponding measurements were computed and reported in Fig. 3. It suggested that the proposed method achieved comparable or better predictive performance than other prediction methods in all cases. For instance, with Sp = 95.0%, Sn, Acc, Pre and MCC values of kinase group CMGC on T site were increased by 29.1%, 3.51%, 10.3% and 21.7% compared with PPSP and had an improvement of 20.4%, 2.41%, 6.76% and 14.7% compared with GPS respectively. In addition, for kinase groups CAMK and CMGC, the precision values obtained by our proposed method were 59.8% and 85.1%, and the precision values of Wang, Jiang & Xu (2015) were 59.7% and 84.2%, respectively. Table S3 showed the detailed comparative results for kinase groups CAMK and CMGC with Sp = 99.0%. From this Table, we can see that our proposed method obtained better performance than other prediction methods, especially sensitivity. For example, for kinase group CMGC the proposed method obtained the Sn value of 50.4%, while the Sn values of GPS, PPSP, NetPhos and Wang, Jiang & Xu (2015) were 16.3%, 13.7%, 15.0% and 38.7%, respectively. Besides, PTMPred (Xu et al., 2014b) was also used to make comparison and the results were illustrated in Table S4, indicating that the proposed method compared favorably with it.

It is known that the control of false positive prediction results is usually critical in the field of computational bioinformatics (Xu & Wang, 2015). Hence, in addition to the aforementioned measurements, we used similar bar plot with those adopted in previous studies (Peng & Li, 2016; Xu & Wang, 2015) to indicate the number of true positives in top-ranked results. For each percentile p%, first we counted the number of true positives in the top ranked p%*total samples, then we calculated the fraction of true positives by dividing total positive samples. Here we took CAMK and CMGC for instance, results of five top 1, 2, 5, 10 and 20 percent of the total samples were compared (Fig. 4). It was observed that the proposed method gave most of the known sites higher ranks than other prediction methods investigated in this study. For example, for kinase group CAMK at the top20%, the fraction of true positives of the proposed method was 78.3% and the corresponding values of GPS, PPSP, NetPhos and Wang, Jiang & Xu (2015) were 55.8%, 40.2%, 36.9% and 74.3%, respectively. Also, Fig. S3 suggested that our method had comparable fraction of predicted sites with other prediction methods. In summary, the proposed method can obtain comparable or better performance for the prediction of phosphorylation sites.

Comparison with existing methods for other PTMs

In this section, we also made comparison with existing methods about other PTMs. For O-GlcNAc the proposed method was compared with several methods including YinOYang (Gupta & Brunak, 2001), O-GlcNAcPRED (Jia, Liu & Wang, 2013) and O-GlcNAcScan (Wang et al., 2011). The detailed ROC curves of different methods were illustrated in Fig. 5. The proposed method achieved an AUC value of 88.0%, and the corresponding AUC values of YinOYang, O-GlcNAcPRED and O-GlcNAcScan were 70.1%, 58.8% and 75.8% on S sites respectively (Fig. 5A). In addition, the AUC values achieved by the proposed method also had an improvement by 26.3%, 34.6% and 18.4% compared with YinOYang, O-GlcNAcPRED and O-GlcNAcScan, respectively on T sites (Fig. 5C). Therefore, the proposed method remarkably outperformed the predictive performance compared with YinOYang, O-GlcNAcPRED and O-GlcNAcScan on both S and T sites. Besides, we also studied the predictive performance of nitration and sulfation on Y sites compared with other existing methods. For sulfation, GPS-TSP (Pan et al., 2014), Sulfinator (Monigatti et al., 2002) and SulfoSite (Chang et al., 2009) were applied to compare the predictive performance, while GPS-YNO2 (Liu et al., 2011) and iNitro-Tyr (Xu et al., 2014c) were compared with the proposed method for nitration. As shown in Fig. 5D, for sulfation, the AUC values were increased by 15.7% compared with GPS-TSP. For nitration (Fig. 5B) the AUC value of proposed method was 7.0% and 30.1% higher than iNitro-Tyr and GPS-YNO2, respectively. Furthermore, the comparison of Sn, Acc, Pre and Sp with multiple types of PTM at the two stringency levels was listed in Table 2. Taking O-GlcNAc on S sites as an example, our proposed method obtained a precision of 45.1% at Sp = 95.0%, while the precision values of O-GlcNAcScan, O-GlcNAcPRED and YinOYang were 34.5%, 26.8% and 14.3%, respectively. For O-GlcNAc on T sites, the Sn value of our proposed method was 28.5% at Sp = 99.0%, while the corresponding values of O-GlcNAcScan, O-GlcNAcPRED and YinOYang were 10.9%, 4.85% and 3.03%, respectively. For sulfation on Y sites, with Sp = 95.0%, the precision value of the proposed method was 77.1% and the corresponding precision values of GPS-TSP, Sulfinator and SulfoSite were 69.4%, 54.3% and 60.9%, respectively. For nitration on Y sites, with Sp = 99.0%, the precision value was increased by 14.3% compared with GPS-YNO2, while was 1.7% lower than iNitro-Tyr, respectively. However with Sp = 95.0%, all measurements were higher than other prediction methods. In conclusion, aforementioned analyses suggested that proposed method outperformed other prediction methods in predicting multiple types of PTM on serine, threonine and tyrosine sites.

Analysis of the predicted potential PTM sites

Due to the difficulty of the experimental verification, the computational method is required to have the ability to detect unknown PTM sites (Xu & Wang, 2015). Hence, we extracted the top ten ranked candidate sites which were not experimentally modified by acetylation or O-GalNAc in our dataset according to the probability estimates of LIBSVM package, respectively. We manually checked these predicted results from UniProtKB database (Boutet et al., 2007) and literature. Table 3 showed the detailed top ten predicted results of acetylation, in which we found that some sites of proteins have been demonstrated to be modified by acetylation. The potential acetylation site with largest probability (0.771) was Thr2 of EBP. Interestingly, we found that this site can be modified by acetylation in the UniProtKB database (http://www.uniprot.org/uniprot/Q15125#ptm_processing). At the same time, in Table S5, we also listed the top ten ranked candidate sites for O-GalNAc. Interestingly, we found according to previous study (Carlsson, Lycksell & Fukuda, 1993) that the Ser207 of protein LAMP2 (probability: 0.702) could be modified by O-GalNAc. These results further demonstrated the proposed method had the ability to discover new target sites, which could be helpful for the subsequent experimental verification.