iCar-PseCp: identify carbonylation sites in proteins by Monte Carlo sampling and incorporating sequence coupled effects into general PseAAC

1Cancer and many other major diseases are often caused by varieties of subtle modifications in biological sequences, typically by various types of post-translational modification (PTM or PTLM) in protein, post-replication modification (PTRM) in DNA and post-transcription modification (PTCM) in RNA. In order to reveal the pathological mechanisms of these diseases and find new and revolutionary strategies to treat them, considerable efforts have been made in order for identifying the possible modified sites in proteins (see, e.g.,), DNA, and RNA sequences. For a systematic introduction about this, see two recent review articles.

2In vivo, PTM is one of the most efficient biological mechanisms for regulating physiology as well as for expanding the genetic code. But when body's well-designed proteolysis or other repair systems are overwhelmed by excess reactive oxygen species (ROS), the oxidative stress may occur, weakening the damage-repairing ability. This may also bring about varieties of PTMs on proteins, including nitration, carbonylation, sulfhydration and glutathionylation. Among these PTMs, the protein carbonylation has been used as a biomarker for severe oxidative protein damage due to its relative early formation, stability, and irreversibility. Actually, protein carbonylation is an early stage of diseases induced by external oxidative stress, aging and obesity. It may cause numerous major human diseases, including Alzheimer's disease, diabetes, Parkinson's disease, chronic renal failure, chronic lung disease, sepsis and so forth. Therefore, the information of carbonylation sites in proteins is indispensable not only for in-depth understanding many important biological processes but also for precisely aiming targets in developing effective drugs against the aforementioned diseases.

3Mass spectrometry is one of the most common techniques to analyze the carbonyl level of a protein and determine its carbonylation sites. So far four types of amino acid residues have been found more prone to carbonylation; they are lysine (K), proline (P), arginine (R), and threonine (T). But it would take much longer time and need more labors to utilize the conventional experimental techniques alone to determine the carbonylation sites in proteins. Facing the rapid growth of biological sequences, we are challenged to develop automated methods as a complimentary approach to experimental methods.

4Actually, some investigators have made efforts to do so. Maisonneuve et al., based on their spectrometry analysis, proposed some empirical rules to identify the hot spots of carbonylation. Recently, Lv et al. and Xu et al. developed two different bioinformatical tools to predict the protein carbonylation sites. These methods did have contribution in stimulating the development of this area. Since the topic's importance as well as the urgency of demanding more powerful high throughput tools in this area, further efforts aiming at prediction of protein carbonylation sites are definitely needed.

5Here, we are to develop a new and more powerful predictor by (1) using the Monte Carlo sampling approach to optimize the training dataset, (2) incorporating the vectorized sequence-coupling model into the general PseAAC, and (3) installing the random forest (RF) algorithm to operate the prediction system.

6As shown in many recent relevant papers, to establish a biological sequence-based statistical predictor that not only can be easily used by most experimental scientists to get their desired results but also can inspirely stimulate theoretical scientists to create various other prediction methods, we should observe the Chou's 5-step rules or guidelines: (1) benchmark dataset preparation; (2) mathematical representation of biological sequence samples; (3) calculation algorithm; (4) cross-validation; (5) web-server establishment. Below, let us to address the five guidelines one-by-one. To match the rubric style of the Oncotarget journal, however, the order in addressing them may be changed.

7A new and more powerful predictor, called iCar-PseCp, has been established for predicting the protein carbonylation sites. Moreover, to maximize users' convenience, the point-to-point instructions are given below.

8Click the web-server at http://www.jci-bioinfo.cn/iCar-PseCp, the top page of the iCar-PseCp will be prompted on your computer screen (Figure 1).

9In the input box (Figure 1), enter your query protein sequences, which can be done by either typing or copying/pasting manner. The entered query protein sequences should be in the FASTA format. Not familiar with FASTA? Just click the button of Example.

10You can see the prediction results by clicking the Submit button. If you use the Sequence_K in the Example window as the input and check on the K button, after 15 seconds or so since your submitting, you will see the following on your screen: Sequence_K contains 9 K residues, of which 5 are predicted to be of carbonylation site and they are at the sequence positions 2, 14, 41, 68 and 95. If you use the Sequence_P as the input and check on the P button, you will see: Sequence_P contains 10 P residues, of which 5 are of carbonylation site and at positions 95, 122, 142, 145, and 149. If you use the Sequence_R as the input and check on the R button, you will see: Sequence_R contains 8 R residues, of which 3 are of carbonylation site and at the positions 14, 41, and 75. If you use the Sequence T as the input and check on the T button, you will see: Sequence_T contains 7 T residues, of which 1 is of carbonylation site and at the positions 14. Compared with experimental observations, the above (9 + 10 + 8 + 7) = 34 predicted results contain no false positive result but 5 false positive results, which are the 2nd and 13th K residues in sequence_K, the 142th and 145th P residues in sequence_P, and the 75th R residue in sequence_R. In other words, the total number of carbonylation sites involved in the above predictions is N+ = 3 + 3 + 1 + 1 = 8, while the total number of non-carbonylation sites investigated is N− = 6 + 7 + 6 + 6 = 25. Substituting these data into Eq.9, we have Sn = 100%, Sp = 80.00% and Acc = 84.80%, and MCC = 0.7018, quite consistent with the rates reported in Table 1 via the rigorous cross validation on the 250 benchmark proteins.

11If you have a lot of query protein sequences and need a lot of computational time, you can choose to use the batch prediction. To do so, just use the Browse button to select the desired file (in FASTA format of course) and follow the online instruction.

12The benchmark dataset used in this study is available by clicking the button of Supporting Information on the top of Figure 1.

13To see the key papers used to develop this server, just click on the button of Citation.

14The success scores achieved by the iCar-PseCp predictor via the 10-fold target cross validation for K-, P-, R-, and T-type carbonylation are shown in Table 1. Meanwhile, the corresponding rates by PTMPred and CarSpred are also listed there. As we can see from Table 1, compared with its counterparts, although the Acc values obtained by the iCar-PseCp are within the ± 4%, its Sn and Sp values are more than 20% and 5–9% higher than those by PTMPred and CarSpred, indicating that the results predicted by the previous methods contain much more false negative and positive events. Particularly, the MCC values achieved by iCar-PseCp are about 2 or 3 times higher than those of its counterparts, indicating that the new proposed predictor is significantly more stable.

15Graphical approach is a useful vehicle for analyzing complicated biological systems as demonstrated by a series of previous studies (see, e.g.,. Here, to provide an intuitive comparison, the graph of Receiver Operating Characteristic (ROC) was utilized to show the advantage of iCar-PseCp over the PTMPred and CarSpred. In Figure 2 the red and green graphic lines are the ROC curves for the PTMPred and CarSpred, respectively; while the blue graphic line for the proposed predictor iCar-PseCp. The greater the area under the AUC is, the better the predictor will be. As we can see from Figure 2, the area under the blue curve is remarkably greater than that under the red or green line, once again indicating that the proposed predictor is indeed much better than PTMPred and CarSpred predictors. Therefore, iCAR-PseCp will become a very useful bioinformatics tool for relevant basic research and drug development as well.

16Why can the proposed method enhance the prediction quality so significantly? First, the coupling effects among the amino acids around the carbonylation sites are taken into account via the conditional probability approach, which has been proved to be indeed very useful in a series of previous studies. Second, the predictor is trained by a balanced benchmark dataset via Monte Carlo sampling, and hence many false prediction events as occurring in the cases of PTMPred and CarSpred trained by very imbalanced and skewed datasets can be completely avoided.

17The data used in the current study was derived from the 230 carbonylated protein sequences from human and 20 carbonylated protein sequences from Photobacterium and Escherichia coli.

18For facilitating description later, the Chou's peptide formulation was adopted. It was used for studying enzyme specificity, signal peptide cleavage sites, hydroxyproline and hydroxylysine sites, methylation sites, nitrotyrosine sites, protein-protein interaction, and protein-protein binding sites. According to Chou's scheme, a potential carbonylation site-containing peptide sample can be generally expressed by

19where the symbol ⊛ denotes the single amino acid code K, P, R, or T, the subscript ξ is an integer, R−ξ represents the ξ-th upstream amino acid residue from the center, the R+ξ the ξ-th downstream amino acid residue, and so forth. The (2ξ + 1) -tuple peptide sample Pξ(⊛) can be further classified into the following two categories:

20where denotes a true carbonylation segment with K, P, R, or T at its center, a false segment with K, P, R, or T at its center, and the symbol ∈ means “a member of” in the set theory.

21In literature the benchmark dataset usually consists of a training dataset and a testing dataset: the former is used for training a model, while the latter for testing the model. But as pointed out in a comprehensive review, there is no need to artificially separate a benchmark dataset into the two parts if the prediction model is examined by the jackknife test or subsampling (K-fold) cross-validation since the outcome thus obtained is actually from a combination of many different independent dataset tests. Thus, the benchmark dataset for the current study can be formulated as

22where the positive subset only contains the samples of true carbonylation segments, and the negative subset only contains the samples of false carbonylation segments (see Eq.2); while ∪ represents the symbol for “union” in the set theory.

23The detailed procedures to construct the benchmark dataset are as follows. (1) As done in, slide the (2ξ + 1) -tuple peptide window along each of the aforementioned 230 + 20 = 250 protein sequences used by, and collected were only those peptide segments that have K, P, R, and T at the center. (2) If the upstream or downstream in a protein sequence was less than ξ or greater than L−ξ where L is the length of the protein sequence concerned, the lacking amino acid was filled with a dummy residue X. (3) The peptide segment samples thus obtained were put into the positive subset if their centers have been experimentally annotated as the carbonylation sites; otherwise, into the negative subset. (4) Using the CD-HIT software, the aforementioned samples were further subject to a screening procedure to winnow those that had ≥ 30% pairwise sequence identity to any other in a same subset.

24Note that the length of peptide samples and their number thus generated would depend on the ξ value. But preliminary tests had indicated that it would be most promising when ξ = 7 or the sample's length was 2ξ + 1 = 15. Accordingly, hereafter we only consider the case of ξ = 7; i.e., the samples with 15 amino acid residues. Thus, the benchmark datasets thus obtained for are given in Supporting Information S1, S2, S3, and S4, respectively. Listed in Table 2 is a summary of their sizes.

25With the avalanche of biological sequence generated in the post-genomic age, one of the most important problems in computational biology is how to formulate a biological sequence with a discrete model or a vector, yet still considerably keep its sequence order information or essential feature. This is because all the existing machine-learning algorithms can only handle vector but not sequence samples, as elaborated in.

26To address this problem, the pseudo amino acid composition or PseAAC was proposed. Ever since the concept of pseudo amino acid composition or Chou's PseAAC was proposed, it has rapidly penetrated into many biomedicine and drug development areas and nearly all the areas of computational proteomics (see, e.g., as well as a long list of references cited in).

27Because it has been widely and increasingly used, recently three powerful open access soft-wares, called ‘PseAAC-Builder’, ‘propy’, and ‘PseAAC-General’, were established: the former two are for generating various modes of Chou's special PseAAC; while the 3rd one for those of Chou's general PseAAC, including not only all the special modes of feature vectors for proteins but also the higher level feature vectors such as “Functional Domain” mode (see Eqs.9–10 of), “Gene Ontology” mode (see Eqs.11–12 of), and “Sequential Evolution” or “PSSM” mode (see Eqs.13–14 of). Inspired by the successes of using PseAAC to deal with protein/peptide sequences, three web-servers were developed for generating various feature vectors for DNA/RNA sequences. Particularly, recently a powerful web-server called Pse-in-One has been developed that can be used to generate any desired feature vectors for protein/peptide and DNA/RNA sequences according to the need of users' studies.

28According to the general PseAAC, the peptide sequence of Eq.1 can be formulated as

29where

30and

31In Eq.5 is the conditional probability of amino acid R−7 occurring at the left 1st position (see Eq.1) given that its closest right neighbor is is the conditional probability of amino acid R−6 occurring at the left 2nd position given that its closest right neighbor is R−5, and so forth. Note that in Eq.5, only and are of non-conditional probability since the right neighbor of R−1 and the left neighbor of R+1 are always ⊛ (namely Lys, Pro, Arg, or Thr, respectively). All these probability values can be easily derived from the positive training subsets taken from Supporting Information S1, S2, S3, and S4, respectively as done in. Likewise, the components in Eq.6 are the same as those in Eq.5 except for that they are derived from the corresponding negative training subsets therein.

32As we can see from the Supporting Information S1, S2, S3, and S4, the negative subset in each of them is much larger than its corresponding positive one in number of samples. Although this might reflect the real world in which the non-carbonylation sites are always the majority compared with the carbonylation ones, a predictor trained by such a highly skewed benchmark dataset would inevitably have the bias consequence that many carbonylation sites might be mispredicted as non-carbonylation ones. Therefore, it is important to find an effective approach to minimize this kind of bias consequence. To realize this, we adopted the Monte Carlo simulation to expand the samples of positive subset. The concrete procedures are as follows.

33Step 1. Suppose (i = −7, −6, …, −1, +1, …, +6, +7; i ≠ 0) is the probability of the 20 native amino acids occurring at the i-th position of the carbonylation samples that can be derived from a training dataset in the positive subsets of Supporting Information S1, S2, S3, or S4, respectively.

34Step 2. For simplicity, let us formulate the probability thus obtained according to the alphabetical order of the single-letter code of the 20 native amino acids (note that the dummy amino acid X introduced in the Benchmark Dataset section was treated as the 21st amino acid); i.e.,

35Step 3. Generate a random number between 0 and 1; if

36then the k-th amino acid is drawn for an expanded positive sample at its i-th subsite. For example, if k = 2 and i = −7, then the amino acid thus drawn should be C for the left 1st sequence position (cf. Eq.1); if k = 19 and i = −6, then the amino acid drawn should be W for the left 2nd sequence position; if k = 20 and i = +7, then the amino acid drawn should be Y for the right last sequence position; and so forth.

37Step 4. Repeat the above steps until the number of positive (the original plus the expanded) samples is the same as the negative samples.

38At first glance, the rationale of the above Monte Carlo sampling procedure seems like a circular argument. But it is correct as elucidated in. Particularly, these expanded positive samples were used only for training a model but not used for testing it, as well be further discussed later.

39The random forests (RF) algorithm is a powerful algorithm and has been used in many areas of computational biology (see, e.g. [,,,, −]). The detailed procedures of RF and its formulation have been very clearly described in, and hence there is no need to repeat here.

40For the current study, all the involved peptide samples were converted into a 14-D (dimensional) vector according to Eq.4, and then entered into the RF operation engine as the input. And the output would indicate whether the center residue ⊛ of the query peptide is a “carbonylation site” or “non- carbonylation site”. Note that, in using the current prediction method, one must observe the self-consistency principle: if the center residue of a query peptide is ⊛ = K then the corresponding training data must be taken from if the center residue of a query peptide is ⊛= P, then the training data must be taken from and; and so forth (see Eq.3).

41The predictor established via the above procedures is called “pCar-PseCp”, where “i” stands for identify”, “Car” for “carbonylation site”, “Pse” for “general PseAAC”, and “Cp” for “sequence coupled effect”.

42As pointed out in the Introduction section, one of the keys in establishing a useful predictor is how to properly evaluate its anticipated success rates. To realize this, we need to consider the following two things: one is what metrics or scales should be adopted to quantitatively measure its prediction quality; the other is what validation method should be utilized to calculate or derive the metrics values. Below, we are to address the two problems.

43The following four metrics are usually used in literature to measure the quality of binary classification: (1) overall accuracy or Acc; (2) Mathew's correlation coefficient or MCC; (3) sensitivity or Sn; and (4) specificity or Sp (see, e.g.,). Unfortunately, the conventional formulations for the four are not intuitive and that most experimental scientists feel difficult to understand them, particularly for the one of MCC. Interestingly, by using the Chou's symbols and derivation in studying signal peptides, the aforementioned four metrics can be easily converted into a set of following equations:

44where N+ represents the total number of carbonylation sites investigated whereas the number of true carbonylation sites incorrectly predicted to be of non-carbonylation site; N− the total number of the non-carbonylation sites investigated whereas the number of non-carbonylation sites incorrectly predicted to be of carbonylation site.

45According to Eq.9, it is crystal clear to see the following. When meaning none of the true carbonylation sites are incorrectly predicted to be of non-carbonylation site, we have the sensitivity Sn = 1. When meaning that all the carbonylation sites are incorrectly predicted to be of non-carbonylation site, we have the sensitivity Sn = 0. Likewise, when meaning none of the non-carbonylation sites are incorrectly predicted to be of carbonylation site, we have the specificity Sp = 1; whereas meaning that all the non-carbonylation sites are incorrectly predicted to be of carbonylation sites, we have the specificity Sp = 0. When meaning that none of carbonylation sites in the positive dataset and none of the non-carbonylation sites in the negative dataset are incorrectly predicted, we have the overall accuracy Acc = 1 and MCC = 1; when and meaning that all the carbonylation sites in the positive dataset and all the non-carbonylation sites in the negative dataset are incorrectly predicted, we have the overall accuracy Acc = 0 and MCC = −1; whereas when and we have Acc = 0.5 and MCC = 0 meaning no better than random guess. Therefore, using Eq.9 has made the meanings of sensitivity, specificity, overall accuracy, and Mathew's correlation coefficient much more intuitive and easier-to-understand, particularly for the meaning of MCC, as concurred recently by many investigators (see, e.g.,).

46Note that, however, the set of equations defined in Eq.9 is valid only for the single-label systems. For the multi-label systems whose emergence has become more frequent in system biology and system medicine, a completely different set of metrics are needed as elaborated in.

47With a good set of metrics to measure the predictor's quality, the next thing to consider is what kind of validation method should be adopted to calculate the metrics values.

48The following three cross-validation methods are often used in statistics to derive the metrics values for a predictor: independent dataset test, subsampling (or K-fold cross-validation) test, and jackknife test. Among these three, however, the jackknife test is deemed the least arbitrary that can always yield a unique outcome for a given benchmark dataset as elucidated in and demonstrated by Eqs.28–32 therein. Accordingly, the jackknife test has been widely recognized and increasingly used by investigators to examine the quality of various predictors (see, e.g.,). However, to reduce the computational time, in this study we adopted the K-fold cross-validation, as done by most investigators with SVM and random forests algorithms as the prediction engine.

49When conducting the K-fold cross-validation for the current predictor iCAR-PseCp, however, some special consideration is needed. This is because a dataset, after expanding by Monte Carlo sampling, may contain many hypothetical positive samples. It would be fine to use such an expanded dataset to train a prediction model, but certainly not for validation. This is because the validation should be made on a testing dataset that only contains experiment-confirmed samples without any added hypothetical samples. To ensure this, a special cross-validation, the so-called target cross-validation, has been introduced here. During the target cross-validation process, only the experiment-confirmed samples are picked out from the testing dataset for validating and scoring. The detailed procedures of the target K-fold cross-validation (without losing the generality, let us consider K = 10) can be described as follows.

50Step 1. Before expanding the positive samples, both the original positive and negative subsets were randomly divided into 10 parts with about the same size. For example, for in Supporting Information S1, after such evenly division we have and where the symbol ≜ means that the divided 10 datasets are about the same in size, and so are their subsets.

51Step 2. One of the 10 sets, say was singled out as the testing dataset and the remaining nine sets as the training dataset.

52Step 3. Based on the training dataset, use Eqs.4–6 to derive the sequence-coupled information. Also, based on the same training dataset, use Monte Carlo sampling to expand its positive subset making it have the same size as the negative subset.

53Step 4. Use the sequence-coupled information and the expanded training dataset obtained in Step 3 to train the model and perform the prediction for each of the samples in the testing dataset.

54Step 5. Repeat Steps 2–4 until all the 10 divided sets had been singled out one-by-one for testing validation.

55Step 6. Substituting the average scores obtained from the above 10-round tests into Eq.9 to calculate Sn, Sp, Acc, and MCC.

56It is crystal clear to see from the above steps that the validation was made only for experiment-confirmed samples, and that none of information from the testing datasets was ever used to train the predictor.

57The iCar-PseCp predictor is a new bioinformatics tool for identifying the carbonylation sites in proteins. Compared with the existing predictors in this area, its prediction quality is much better, with remarkably more stability and less false predictions. For the convenience of most experimental scientists, we have provided its web-server and a step-by-step guide, by which users can easily obtain their desired results without the need to go through the detailed mathematics. The reason of including them in this paper is for the integrity of the new prediction method, and that these techniques, such as sequence-coupled approach and Monte Carlo sampling, may be of use as well in developing other tools in computational biology.

58We anticipate that iCar-PseCp will become a very useful high throughput tool, or at the very least, a complementary tool to the existing methods for predicting the protein carbonylation sites.