Figure 2 shows the procedure of dataset construction and preprocessing. 7,520 human SADPr sites with high confidence (i.e., ADPr peptides with Andromeda scores >40 and localization probability >0.75) were collected from the literature (Larsen et al., 2018; Hendriks et al., 2019; Buch-Larsen et al., 2020; Nowak et al., 2020) (Figure 2A). 151,227 human pS sites were obtained from the database PhosphositePlus (Hornbeck et al., 2012) and the literature (Luo et al., 2019) (Figure 2A). We compared both datasets and found 3,250 pSADPr peptides, 147,977 pS peptides, and 4,270 SADPr peptides. We also collected 80,096 unmodified serine (UM) sites after removing modified serine sites (i.e., pSADPr, SADPr and pS) from the reported dataset (Luo et al., 2019).

Each serine site of the above datasets was represented by a 41-residue-long sequence segment with the serine at the center (Sha et al., 2021). CD-HIT (Li and GodzikCd-hit, 2006; Huang et al., 2010) was applied to eliminate the homologous peptides by setting the threshold to 60% sequence identity, which is valuable for avoiding overestimation. Specifically, we combined the pSADPr peptides with SADPr peptides, pS peptides, and UM peptides, respectively, and clustered them using CD-HIT. Accordingly, we obtained 4,959 clusters, 30,106 clusters and 66,526 clusters. We selected one sequence randomly from each cluster according to the criterion: One pSADPr peptide was selected if it was included in the cluster; otherwise, one of the other peptides was selected. After that, 2,378 pSADPr, 2,581 SADPr, 27,728 pS and 64,148 UM peptides were collected (Figures 2B–D). Furthermore, each of the three datasets was divided into 11 groups, where ten groups were used as a cross-validation dataset, and the rest group was considered an independent test dataset (Figures 2B–D). It should be noted that if the central serine residue is located near the N or C terminus of the protein sequence, the complement symbol “_” was added to the input sequences at the affected terminus to ensure the length was maintained. All these data are available at http://edeepsadpr.bioinfogo.org/.

We selected five encoding features representing the input peptides for the model construction. They included the One-Hot encoding (OH) (Wang D. et al., 2020), the Enhanced Amino Acid Composition Encoding (EAAC) (Chen et al., 2018), the Enhanced Grouped Amino Acids Content encoding (EGAAC) (Chen et al., 2018), the ZSCALE Encoding (ZSCALE) and the Word Embedding (WE).

We constructed five classifiers based on Convolutional Neural Network (CNN). They included the model combined with the One-Hot Encoding (CNN_OH), the model with the Word Embedding Encoding (CNN_WE), the model with the ZSCALE Encoding (CNN_ZSCALE), the model with the EAAC encoding (CNN_EAAC) and the model with the EGAAC encoding (CNN_EGAAC). We took the CNN Model with the One-Hot encoding (CNN_OH) as an example to demonstrate the architecture (Figure 3).

(1) Input layer. Each sequence is converted into a feature vector with One-Hot encoding.

Several statistical measures were used to evaluate prediction performance, including sensitivity (SN), specificity (SP), overall accuracy (ACC), Matthew correlation coefficient (MCC) and the area under the receiver operating characteristic (ROC) curve (AUC). The definitions of SN, SP, ACC, and MCC are given as follows: S N = T P T P + F N S P = T N T N + F P A C C = T P + T N T P + F P + T N + F N M C C = T P × T N − T N × F P T P + F N × T N + F P × T P + F P × T N + F N

In the above formulas, TP, TN, FP, and FN are the number of true positives, true negatives, false positives, and true negatives, respectively.

We created three datasets for constructing classifiers to predict pSADPr sites (Figure 2). The first dataset was the pSADPr-SADPr dataset, containing pSADPr and SADPr peptides. The related model was used to recognize pSADPr sites from known SADPr sites (Figure 2B). The second was the pSADPr-pS dataset, including pSADPr and pS peptides (Figure 2C). The third was the pSADPr-UM dataset, containing pSADPr and UM peptides (Figure 2D). Because the vast majority of serine residues are unmodified in the human proteome, the model based on the third dataset was expected to recognize pSADPr sites from the human proteome (Figure 2D). Each of the three datasets contained two parts: cross-validation and independent test datasets (Figures 2B–D).

We explored the characteristics of the pSADPr crosstalks by comparing pSADPr-containing and other peptides in the three datasets through the Two-Sample-Logo program (Vacic et al., 2006). For the pSADPr-SADPr dataset, the amino acid R was significantly enriched at positions −2 and −3 (i.e., P-2 and P-3), whereas K was depleted at P-1 (Figure 4A). For the rest datasets, the pSADPr crosstalks showed similar characteristics (Figures 4B, C). Specifically, K was enriched entirely except P+1 and G was enriched at P1 and P2; D and E were depleted at P-3 to P+5 and L was depleted entirely. The maximum enriched/depleted value (29.3%) for the pSADPr-pS dataset was similar to that (32.0%) for the pSADPr-UM dataset, and both were more than twice as large as that (13.2%) for the pSADPr-SADPr dataset (Figure 4). It indicates that the differences between pSADPr and SADPr sites are smaller than those between pSADPr and pS/UM sites. In other words, it is easy to distinguish pSADPr sites from pS/UM sites, compared to recognizing pSADPr sites from SADPr sites.

The human serine kinase family contains a few subfamilies, each with its characteristics. We explored which subfamilies preferred phosphorylating the pSADPr sites. To perform this analysis, we used the human pS sites as the background and the pSADPr sites as the test dataset. We employed the GPS program (Wang C. et al., 2020) to predict pS sites for each subfamily from both datasets (Figure 5). We found that four subfamilies (i.e., AGC, CAMK, STE and TKL) tended to phosphorylate pSADPr sites (p < 5.0 × 10⁻²⁶, hyper-geometric test). In comparison, two subfamilies (i.e., CK1 and CMGC) prefer not to phosphorylate pSADPr sites (p < 5.1 × 10⁻²⁹, hyper-geometric test). For example, 68% of pSADPr sites could be phosphorylated by the AGC subfamily, whereas only 44% of pS sites are modified by this subfamily (p = 2.3 × 10⁻¹⁷⁴, hyper-geometric test). This observation suggests that the pSADPr sites may be related to specific subfamilies of serine kinases.

In the three datasets, the pSADPr-pS and pSADPr-UM datasets were imbalanced because the numbers of pS and UM peptides were far more than the number of pSADPr peptides (Figures 2C, D). To explore the effect of the imbalanced dataset on the predictor’s performance, we built the related balanced cross-validation dataset where the number (2,162) of randomly selected pS or UM peptides was the same as that of pSADPr peptides. We constructed the CNN_OH models related to the imbalanced and balanced datasets and evaluated their prediction performances in terms of the independent test. The CNN_OH model based on the imbalanced dataset had better performance than the counterpart constructed using the balanced dataset (p = 0.002 for both pSADPr-pS and pSADPr-UM datasets, Wilcoxon rank sum test; Figure 6). Therefore, we chose the imbalanced dataset for model construction.

We constructed five CNN classifiers (i.e., CNN_OH, CNN_WE, CNN_EAAC, CNN_EGAAC and CNN_ZSCALE) to recognize pSADPr sites from the three datasets and compared their prediction performances. Here, we used the pSADPr-SADPr dataset to demonstrate the process. Three out of the five classifiers (i.e., CNN_OH, CNN_WE and CNN_ZSCALE) showed similar performances and superiority over the rest two (i.e., CNN_EAAC and CNN_EGAAC) in ten-fold cross-validation and independent test (Table 1; Figure 7 and Supplementary Figure S1). For instance, the CNN_OH model had an AUC value of 0.712, larger than that (0.659) of the CNN_EAAC model in the cross-validation. We repeated this analysis for the pSADPr-pS and pSADPr-UM datasets and made similar observations that the three classifiers had the best performances (Supplementary Tables S1, S2; Supplementary Figures S1–S5). Furthermore, we compared the classifiers’ performances for the three datasets. We found that the AUC values (0.921 and 0.953) of the CNN_OH classifiers for pSADPr-pS and pSADPr-UM datasets were significantly larger than that (0.712) for the pSADPr-SADPr dataset. These results were consistent with our observation that the differences between pSADPr and SADPr sites are smaller than those between pSADPr and pS/UM sites (Figure 4). Since the One-Hot feature is the simplest compared to the WE and ZSCALE features, we chose the CNN classifier with the One-Hot scheme as the representative of the three classifiers.

A stacking-based ensemble learning architecture is one of the ensemble techniques in which multiple learning models are integrated to produce one optimal predictive model, which performs better than the base models taken alone. In the stacking ensemble architecture, a meta-learner is trained to output a prediction based on the different base learner’s predictions. The stacking ensemble architecture has been used to improve the prediction performance in various bioinformatics applications (e.g., lysine acetylation site prediction) (Mishra et al., 2019; Zhang et al., 2021; Basith et al., 2022). Here, we introduced the two-stage stacking ensemble approach to improve the performance of the pSADPr site prediction (Figure 8). In the first stage, different CNN algorithms (e.g., CNN_OH, CNN_WE and CNN_ZSCALE) were selected to construct base classifiers. Specifically, ten base classifiers for each CNN algorithm were built and validated using the ten-fold cross-validation dataset. The base classifiers were then used for prediction in the independent test dataset, and their prediction results were averaged. Therefore, each CNN algorithm corresponds to the validation result and the averaged result for the independent test dataset. In the second stage, the validation and the averaged results were merged as a meta cross-validation dataset and a meta-independent test dataset, respectively (Figure 8). The former dataset was used to train and validate a meta-classifier, whereas the latter was employed to evaluate the meta-classifier’s performance. Here, we constructed the meta-classifier using the random forest algorithm (RF), which was optimized using the GridSearchCV package. The optimized parameters included max_depth as 8, max_features as “sqrt,” min_samples_leaf as 20, min_samples_split as 300 and n_estimators as 100.

According to the above analysis, the three classifiers (i.e., CNN_OH, CNN_WE and CNN_ZSCALE) had better performances than two other classifiers (i.e., CNN_EAAC and CNN_EGAAC) for all three datasets. Based on the observation, we fused them as base classifiers to build the two-stage stacking ensemble approach with a good performance. We started with the fusion of the three best classifiers until we fused all the classifiers. The related stacking models included Stacking_O+Z+W, Stacking_O+Z+W+E and Stacking_O+Z+W+E+EG, where O stands for OH, Z for ZSCALE, W for WE, E for EAAC and, EG for EGAAC. For the pSADPr-SADPr dataset, the three stacking models showed similar performances in meta ten-fold cross-validation and independent test (Table 2; Figure 9 and Supplementary Figure S6). For instance, their average AUC/MCC values were around 0.719/0.313 in cross-validation (Table 2). The stacking models for the two other datasets (pSADPr-pS and pSADPr-UM) also performed similarly (Supplementary Figures S7–S10).

We compared the performances of the CNN-based models and the stacking ensemble models for each of the three datasets. We found no statistical difference between the CNN_OH model and these stacking ensemble models for each dataset (Figure 9 and Supplementary Figures S9, S10). The observation that the meta-classifiers perform similarly to the base classifier is consistent with the previous report for predicting bacterial Type IV secreted effectors, in which the meta-classifier and base classifier performed similarly (Xiong et al., 2018). It suggests that the base classifiers may have sufficient predictive ability, and the stacking ensemble architecture does not constantly improve prediction accuracy.

We developed an online prediction tool for predicting human pSADPr sites extensively from different conditions, dubbed EdeepSADPr. This tool consists of three models, each corresponding to the prediction from the SADPr dataset, the serine phosphorylation dataset or the human proteome. As the CNN_OH classifier had no less predictive performance than other methods, we selected this classifier to construct EdeepSADPr. The usage of this tool was described as follows. After the model selection, the input sequence with the fasta format would be uploaded. The prediction results were output in tabular form with five columns: sequence header, position, sequence, prediction score, and prediction category. The predicted results can also be downloaded as a data file. EdeepSADPr is accessible via http://edeepsadpr.bioinfogo.org/.