On the existence of Hopf bifurcations in the sequential and distributive double phosphorylation cycle

1.   Introduction

Protein post-translational modification is an important chemical process, which plays a key role in regulating cell functions ^[1] and also changes the physical and chemical properties of protein. More than 400 post-translational modifications including methylation ^[2], acetylation ^[3], phosphorylation ^[4], and S-nitrosylation (SNO) ^[5] have been discovered so far. As SNO is a reversible post-translational modification of proteins, a large number of studies have shown that SNO plays an important role in multiple biological processes such as redox signal transduction ^[6], cell signal transduction ^[7], cell senescence ^[8], and transcription ^[9]. SNO is also related to many human diseases such as cancer ^[10], Alzheimer's disease ^[11], and chronic renal failure ^[12]. Therefore, a well-grounded understanding of SNO is of great significance for the study of basic biological processes ^[9,13] and the development of drugs ^[14]. In recent years, many SNO sites have been identified through molecular signals ^[15,16], but identification of SNO sites still faces some challenges including low accuracy, time-consuming and labor-intensive. With the continuous development of computer technology, a large number of computational models have been used to predict the specific sites of SNO modification.

Many post-translational modifications of proteins have been detected by a variety of computational models. Qiu identified phosphorylated ^[17] and acetylated ^[18] proteins with GO notations. GPS-SNO ^[19], SNOSite ^[20], iSNO-PseAAC ^[21], PreSNO ^[5] and RecSNO ^[22] have been applied to the prediction of SNO sites. The GPS-SNO, SNOSite and iSNO-PseAAC models use relatively small data sets. In addition, many negative samples in these data sets are now experimentally verified as positive samples. The data sets used by PreSNO and RecSNO are relatively large and new, but there is still room for improvement in the performance of the model.

On the basis of previous research, this work established two models for predicting SNO proteins and sites. In predicting SNO proteins, a bag-of-words model has been proposed on the basis of KNN scoring matrix obtained from proteins' GO annotation information ^[18], PseAAC ^[23,24] of amino acids sequence. Fusion of multiple features can more comprehensively reflect the information of the protein sequence and improve the prediction results. A combination of oversampling technique and random deletion method are applied to balance the training set since the issue is involved in imbalanced data sets. In predicting SNO sites, two feature extraction methods, TPC ^[25] and CKSAAP ^[26], are used to extract the features of protein sequence fragments. In order to eliminate the redundancy and noise information of the original feature space, elastic nets ^[27] are used to reduce the dimensionality of the feature space after the fusion strategy is performed on the original features. Random Forest severed as the classifiers and be verified with 5-fold cross-validation. The specific flow chart is shown in Figure 1.

2.   Materials

2.2.   Datasets for prediction SNO proteins

To obtain a scientifical prediction result, a strict benchmark data set is essential. The UniProKB has been accepted by most bioinformatics researchers. Here, the negative samples are extracted from the UniProKB and the positive sample are extracted from Xie's ^[28], which is a high-quality data set based on extensive literature research. The protein sequence can be expressed as:

where ${R}_{i}$ represents the $i$-th amino acid residue, and $L$ represents the length of the protein sequence.

In order to identify SNO proteins, we constructed a benchmark data set similar to dataset of Hasan et al. ^[5], which consists of 3113 SNO proteins. Every one of the positive samples, i.e., SNO proteins, has at least one SNO site. For negative samples, we randomly selected 18, 047 proteins without any SNO site from the UniProKB. In order to make the results more rigorous, the CD-HIT was used to remove 30% of the 3113 positive samples and 18, 047 negative samples. Finally, 2192 positive samples and 7809 negative samples are collected in the proposed benchmark data set.

2.2.   Datasets for predicting SNO Sites

The benchmark data set for predicting SNO sites are the same as Hasan et al. ^[5], which consists of 3383 positive samples and 3365 negative samples. A potential SNO(C) site-containing peptide sample can be generally expressed by

where the subscript $\mathrm{\xi }$ is an integer, ${R}_{-\mathrm{\xi }}$ represents the $\mathrm{\xi }$-th upstream amino acid residue from the center, the ${R}_{+\mathrm{\xi }}$ the $\mathrm{\xi }$-th downstream amino acid residue, and so forth. If the number of left or right residues of the center $C$ is less than $\mathrm{\xi }$, then the pseudo amino acid "$X$" would be used to supplement the sequence. The $(2\mathrm{\xi }+1)$-tuple peptide sample $\stackrel{-}{{\bf{P}}_{\mathrm{\xi }}}$ can be further classified into the following two categories:

where $\overline{\bf{P}}_{\mathrm{\xi }}^{+}$ denotes a true SNO segment $C$ with at its center, $\overline{\bf{P}}_\overline{\mathrm{\xi }}$ a corresponding false SNO segment, and the symbol $\in$ means "a member of" in the set theory.

3.   Feature extraction

3.1.   Feature extraction methods for predicting SNO proteins

3.1.1.   GO-KNN

GO-KNN ^[18] features were extracted on the basis of the GO annotations of proteins. In this work, we need to find out the GO terms of all protein sequences and calculate the distance between proteins. Take protein ${P}_{1}$ as an example, for anyone of other proteins, for example, ${P}_{2}$, then their GO terms can be listed as ${P}_{GO}^{1} = \left\{{GO}_{1}^{1}, {GO}_{2}^{1}, \cdots, {GO}_{M}^{1}\right\}$ and ${P}_{GO}^{2} = \{{GO}_{1}^{2}, {GO}_{2}^{2}, \cdots, {GO}_{N}^{2}\}$ are obtained. If there is no GO term for a protein, we will replace it with GO terms of its homologous protein. The distance between two proteins can be calculated with Eq (4):

where, ${GO}_{i}^{1}$ and ${GO}_{i}^{2}$ represent the $i$-th GO of ${P}_{1}$ and ${P}_{2}$, respectively. The $M$ and $N$ represent the numbers of GO, respectively, ∪and ∩ are the union and intersection in the set theory, and ⌊ ⌋ represents the number of elements in the set. Then, the GO-KNN features could be extracted according to the following steps: 1) Sorting the calculated distances in ascending order; 2) Selecting the first $k$ near neighbors of the test protein; 3) Calculating the percentage of positive samples in the $k$ neighbors. In this study, $k$ were selected as 2, 4, 8, 16, 62, 64, 128, 256, 512, 1024. In this way, a 10-dimensional feature vector $\left({x}_{1}{, x}_{2}, \cdots, {x}_{10}\right)$ could represents the protein ${P}_{1}$.

3.1.2.   BOW

A bag-of-words model ^[29] based on the physical and chemical properties of protein has been used in identifying GPCR-drug interaction. The main steps are listed as follows: 1) Encoding the protein sequence with its physical and chemical properties. Up to now, scientists have obtained various physical and chemical properties of 20 common amino acids ^[30]. After careful experimental comparison, hydrophilicity was selected as an indicator for the proposed model. 2) Designing wordbooks for protein. When the window sizes are 1, 2 and 3, and the step size of the moving window is 1, the coding sequence is divided into segments of different lengths. Segments of length 1 form wordbook ${WB}_{1}$, segments of length 2 form wordbook ${WB}_{2}$, and segments of length 3 form wordbook ${WB}_{3}$. When the window size is 2, the step size of moving the window is still 1. But the window at this time is different from the above, it is separated by an amino acid. At this time, the coding sequence is divided into fragments of length 2, and these fragments form the wordbook ${WB}_{4}$. 3) Clustering the word books. We divided the words in the wordbook ${WB}_{1}$ into 20 sub-groups according to the types of amino acids. Words in ${WB}_{2}, {WB}_{3}$ and ${WB}_{4}$ were clustered with K-means algorithm. The numbers of clusters were 16, 62 and 16. 4) Calculating the ratio of the number of each amino acid to the total number of words in the vocabulary with Eq (5).

here, $K$ is the number of clusters in the wordbook ${WB}_{j}$, ${X}_{i}^{{WB}_{j}}$ is the number of words in the $i$-th category of the wordbook ${WB}_{j}$, and $N$ is the total number of words in the wordbook ${WB}_{j}$. Then a 114-D feature vector was formed for a given protein sequence, i.e. $\left({X}_{1}^{{WB}_{1}}, \mathrm{ }\dots, \mathrm{ }{X}_{20}^{{WB}_{1}}, {X}_{1}^{{WB}_{2}}, \dots {, X}_{16}^{{WB}_{2}}, \mathrm{ }{X}_{1}^{{WB}_{3}}, ...{X}_{62}^{{WB}_{3}}, {X}_{1}^{{WB}_{4}}...{, X}_{16}^{{WB}_{4}}\right)$.

3.1.3.   PseAAC

PseAAC ^[23,24] is a very popular feature for bioinformatics. In this work, six physical and chemical properties of hydrophobicity, hydrophilicity, molecular side chain mass, PK1, PK2 and PI are used. We first used Eq (6) to transform the original physical and chemical properties of amino acids:

where, $a\in \{\mathrm{1, 2}, \cdots, 6\}$ and $i\in \{\mathrm{1, 2}, \cdots, 20\}$. ${W}_{a}^{0}\left(i\right)$ represents the value of the $a$th original physical and chemical properties of the $i$th amino acid. We substitute the values of the transformed physical and chemical properties with Eq (7):

where, ${W}_{1}\left({R}_{j}\right)$represents the value of hydrophobicity of ${R}_{j}$. By analogy, ${W}_{6}\left({R}_{i}\right)$represents the PI value of ${R}_{i}$. Then the correlation factor of each layer can be obtained by using the Eq (8):

where ${\theta }_{\lambda }$ represents the correlation factor of the $\lambda$-th layer of the protein sequence. Finally, the protein sequence is converted into a feature by Eq (9):

here, ${f}_{i}$ represents the frequency of the $i$-th amino acid, $\omega$ is 0.5, and $\lambda$ is 5. In this way, a 25-dimensional feature vector is formed.

3.2.   Data balancing methods for predicting SNO protein

In order to reduce the adverse effect of unbalanced data on the performance of the model, many methods for dealing with unbalanced data have been proposed, such as Synthetic Minority Oversampling Technique ^[31] (SMOTE) and Random Under Sampler ^[32] (RUS). SMOTE is a method proposed by Chawla et al. It has been used to predict protein sites ^[27] and improve the prognostic assessment of lung cancer ^[33]. RUS is a very simple and popular method of under-sampling. It can be used in pediatric pneumonia detection ^[34] and convolutional neural network performance improvement issues ^[35]. In this study, we combined these two methods to process the data. SMOTE is used to oversample the positive samples, and RUS is used to under-sample the negative samples. In the end, the number of processed positive samples is equal to negative samples. The specific process is shown in Figure 2.

3.3.   Feature extraction methods for predicting SNO sites

3.3.1.   CKSAAP

CSKAAP ^[25] has been widely used in protein site prediction ^[26] since it can effectively express internal laws for a given protein sequence. The protein fragment is composed of 20 common amino acids and a pseudo amino acid, which contains 441 residue pairs (AA, AC, ..., XX) for each $l$. Here $l$ represents the space between each residue pair. The following formula is used to calculate the characteristics of the fragment:

here ${N}_{AA}, {N}_{AC}, \cdots$ represent the number of times the corresponding amino acid pair appears in the fragment, $L$ is the length of the protein fragment. ${N}_{T} = L-l-1$. In this study, the values of $l$ are 0, 1, 2, 3, 4, and the corresponding ${N}_{T}$ are 40, 39, 38, 37 and 36, respectively. Then, a 2205-D feature vector is formed.

3.3.2.   TPC

Based on the structural properties of proteins, researchers have proposed the tripeptide composition (TPC). It has been used to predict protein subcellular localization ^[36] identification of plasmodium mitochondrial proteins ^[37]. TPC calculates the frequency of three consecutive amino acids, and then a protein fragment can be represented by a 9261-dimensional vector.

where ${N}_{i}$ represents the total number of $i$-th in 9261 tripeptides.

3.4.   Feature selection for predicting SNO sites

The elastic net proposed by Zou and Hastie ^[38] is an effective feature selection method. By introducing the ${L}_{1}, {L}_{2}$ norm into a simple linear regression model, the elastic net can not only perform continuous shrinkage and automatically select variables at the same time, but also predict related variables. At present, elastic nets have been widely used in protein site prediction ^[27,39] and achieved good results.

4.   Model evaluation metrics and operation engine

4.1.   Model evaluation metrics

In this study, four indicators were used to evaluate the performance of the models. They are accuracy (ACC), sensitivity (SN), specificity (SP) and Matthews correlation coefficient (MCC) ^[40], which are defined by Eq (12):

In predicting SNO proteins, $TP$ indicates the number of proteins that are predicted to have SNO sites and actually have SNO sites, and $TN$ is the number of proteins that are predicted to have no SNO sites that are actually not have SNO sites. $FP$ is the number of proteins without SNO sites but predicted to have SNO sites, $FN$ is the number of proteins with SNO sites but predicted to have no SNO sites. In addition, the area under the ROC curve AUC is also used to evaluate this model.

In predicting SNO sites, TP indicates the number of actual SNO sites predicted to be SNO sites, and TN indicates the number of non-SNO sites predicted to be not SNO sites. FP is the number of non-SNO sites predicted to be SNO sites, and FN is the number of actual SNO sites predicted to be non-SNO sites.

4.2.   Operation engine

4.2.1.   Random Forest

Random Forest ^[41] is an algorithm that integrates multiple trees through the idea of ensemble learning. Its basic unit is decision tree. As a highly flexible machine learning algorithm, Random Forest (RF) has been widely used in data analysis ^[42], biological information ^[43] and technological development ^[44].

4.2.2.   Naive Bayes

Naive Bayes (NB) ^[45] is a simple and effective classifier, which is widely used in software defect prediction ^[46], medical diagnosis ^[47] and biological information ^[48]. NB is based on the Bayes theorem and the assumption of the conditional independence of features, which greatly reduces the complexity of the classification algorithm.

4.2.3.   K Nearest Neighbor

K Nearest Neighbor (KNN) ^[49] is one of the supervised machine learning algorithms, which is widely used in face recognition ^[50], disease research ^[51] and engineering applications ^[52]. Its main idea is to judge the category of the predicted value based on the category of the $k$ points closest to the predicted value.

4.2.4.   eXtreme Gradient Boosting

XGBoost ^[53,54] is an improved algorithm for boosting based on GBDT ^[55]. XGBoost is an integrated lifting algorithm that integrates many basic models to form a strong model. Because of its advantages such as good prediction effect and high training efficiency, XGBoost has been widely used in the field of data analysis.

5.   Results and discussion

5.1.   Results and discussion for SNO proteins prediction

5.1.1.   Effect of features

In this research, GO-KNN, BOW and PseAAC three kinds of feature extraction methods were used to encode the protein sequence, and obtained 10-D, 114-D and 25-D feature vectors, respectively. These three kinds of features were fused into a 149-D feature vector ALL. Through the 5-fold cross-validation, the prediction results obtained by different feature extraction are shown in Table 1.

It can be seen from Table 1 that different features obtained varied prediction results. Among the three methods, GO-KNN has the highest ACC, Sn, MCC and AUC, of which are 82.72%, 48.36%, 0.4533 and 0.8521 respectively. The ACC, Sn, MCC and AUC of BOW are the lowest, of which are 78.83%, 13.04%, 0.1969 and 0.7359, respectively. But the Sp of BOW is 97.30%, which is the highest. After combining these three characteristics, ACC, Sn, Sp, MCC, AUC are 83.77%, 49.49%, 93.40%, 0.4840, 0.8593, respectively. Among them, ACC, Sn, MCC and AUC are all higher than those produced by GO-KNN. The results show that multi-feature fusion can improve a number of indicators. In order to better analyze the influence of different features on the prediction of SNO proteins, the prediction results obtained from the three features and their fusion features are shown in Figure 3.

It can be seen from Figure 3 that the three features and their fusion affect the five evaluation indicators to some extent. They are less effective on Sn and MCC, and better on ACC, Sp and AUC. Comparing these four feature codes, the ACC, Sn, MCC and AUC of the fusion feature ALL are improved. Multi-feature fusion can reflect sequence information more comprehensively, thereby improving prediction ability. Therefore, multi-feature fusion can be used to predict SNO proteins.

5.1.2.   Effect of the SR balancer

Here, SMOTE and RUS are denoted as SR balancer. We input the pre-balanced and post-balanced data sets into the model, and passed the 5-fold cross-validation to obtain the prediction results of ACC, Sn, Sp, MCC, AUC on the balanced and unbalanced data sets, as shown in the Table 2.

It can be seen from Table 2 that the balanced Sn and Sp are relatively balanced. In addition, Sn, MCC and AUC have improved. Therefore, in summary, it is very necessary to balance the dataset.

5.1.3.   Effect of classifiers

Classifiers play an important role in model prediction. This work used the above four classifiers to identify SNO proteins. After 5-fold cross-validation, the results of each classifier for ACC, Sn, Sp, MCC and AUC are shown in Table 3. It can be seen from Table 3 that the effect of random forest on various evaluation indicators is the best. In order to better compare the effects of different classifiers, the prediction results of the four classifiers are shown in Figure 4.

The area under the ROC curve can evaluate the predictive performance of the model. It can be seen from Figure 4 that when the random forest is used as a classifier, the area under the ROC curve is the largest. Therefore, random forest is the best choice for the proposed model.

5.2.   Results and discussion for SNO sites prediction

5.2.1.   Effect of features

In this study, two kinds of features, CKSAAP and TPC, were used, and the 2205-dimension and 9261-dimension feature vectors were obtained on the basis of above-mentioned algorithms. In order to better reflect the information of protein fragments, these features are fused into a 11, 466-dimension feature vector. Through 5-fold cross-validation, the prediction results obtained by different feature extraction are shown in Table 4.

It can be seen from Table 4 that the ACC, Sn and MCC of CKSAAP are higher than those of TPC. TPC performs better than CKSAAP on Sp and AUC. After feature fusion, Acc, Sn, MCC and AUC are all higher than single feature. Therefore, feature fusion is necessary for this issue.

5.2.2.   Effect of elastic net

Multi-information fusion can more comprehensively extract protein sequence information, but redundancy and noise information will also be generated. The dimensionality reduction method can not only retain important features, but the computational efficiency of the model will also be improved. In this paper, elastic net was used to reduce the dimensionality of the fused feature data set, and obtain the feature subset of 704. After 5-fold cross-validation, the prediction results of Random Forest are shown in Table 5.

The features after dimensionality reduction using elastic nets, except for Sn, all other evaluation indicators have been improved. In addition, because the feature dimension is greatly reduced after dimensionality reduction, the efficiency of the model is also significantly improved.

5.2.3.   Effect of different classifiers

Four kinds of classifiers, Random Forest, Naive Bayes, K-Nearest Neighbor and XGBoost, were tested in this work for predicting SNO sites. After 5-fold cross-validation, the results were shown in Table 6. From Table 6 we can get that Naive Bayes and K-Nearest Neighbors are relatively inferior. Except for Sp, all indicators of Random Forest were the best. In order to evaluate the performance of the classifier more comprehensively, the ROC curves of different classifiers are shown in Figure 5.

From Figure 5, we can clearly see that the area under the ROC curve of the random forest is the largest. Therefore, random forest has been selected as the classifier of the proposed model.

5.2.4.   Comparison with other methods

To further evaluate the performance of this model, and we compared it with the PreSNO and RecSNO models. The prediction results of three different methods for the same data set are shown in Table 7. From Table 7, we can see that the ACC, Sn and MCC models of this model are all the highest. In addition, the Sp and AUC of the model in this paper also have good results. Therefore, the performance of this model is better than PreSNO and RecSNO.

6.   Conclusions

In order to identify SNO proteins, we used GO-KNN, BOW and PseAAC to extract the sequence information. GO-KNN extracted KNN neighbor information based on protein GO information, and BOW and PseAAC extracted protein sequence information based on physical and chemical properties. In addition, we used the SR balancer to process the unbalanced data set, reduce the negative impact of the unbalance on the model. Finally, Random Forest was used to make predictions. For predicting SNO sites, CKSAAP and TPC were used to extract protein fragment information. In order to improve the computational efficiency and eliminate the redundancy and noise generated by the fusion features, we used elastic nets to reduce the dimensionality of the fusion features. These processes only need to require calculation models without any physical and chemical experiments, which can save experimental costs and improves work efficiency. We hope that this work will be helpful for solving biological problems with computational methods.

Acknowledgments

This work was supported by the grants from the National Natural Science Foundation of China (No. 31760315, 62162032, 61761023), Natural Science Foundation of Jiangxi Province, China (NO. 20202BAB202007).

Conflict of interest

The authors have declared that no competing interest exists.