iRNAm5C-PseDNC: identifying RNA 5-methylcytosine sites by incorporating physical-chemical properties into pseudo dinucleotide composition

1Post-transcriptional modifications (PTCM) of RNA plays a paramount role for the metabolism processes of RNAs, such as for their splicing export, immune tolerance, and transcription. So far, more than 100 distinct PTCMs have been identified in tRNAs, rRNAs, Mt-tRNAs, miRNAs, lincRNAs, miscRNAs, protein-coding genes, pseudogenes, etc.. Among these modifications, the methylation of the 5-methylcytosine (m5C) is an epigenetic one formed by the action of RNA methyltransferases (Figure 1). The m5C modification is well investigated in DNA, but the corresponding studies in cellular RNA were mainly confined to tRNA and rRNA.

2Actually, the m5C modification site in RNA has various biological functions, including the one that can regulate RNA metabolism in both eukaryotes and prokaryotes. Furthermore, it plays a key role in yeast cell fate decision. It is also significant for animal (such as mouse) and human embryonic development.

3Although many efforts have been made by using biological experiments to determine the m5C sites in RNA (see, e.g.), it is time-consuming and expensive to completely rely on the experimental approaches alone. Facing today's explosive growth of uncharacterized RNA sequences, it is highly demanded to develop computational approach to help getting the information.

4Very recently, in a pioneering study, Feng et al. proposed an interesting method to identify RNA m5C sites via the powerful PseKNC approach. But no web-server has been provided for their method, and hence its practical application value is quite limited. In view of this, the present study was initiated to fill such an empty area.

5A predictor called “iRNAm5C-PseDNC” has been established. The success rates achieved by it on the benchmark dataset constructed based on experimental servations (Supplementary Information 1) are

6where the definitions for the metrics Sn, Sp, Acc, and MCC are given in Eq.13 of the MATERIALS AND DISCUSSION section later.

7Since it is the first web-server predictor ever developed for identifying the m5C sites in RNA sequences, it is not possible to demonstrate its power by comparing with its counterparts for exactly the same purpose. Nevertheless, we can indirectly show its power via a cohort of the anticipated success rates (Table 1) reported from the five powerful web-server predictors in genome and proteome analyses. As we can see from the table, the iRNAm5C-PseDNC is with the highest score for Acc (see column 3), and the same is true for MCC (column 4), indicating the proposed predictor is not only high in overall accuracy but also quite stable.

8Also, it is instructive to point out that, among the four metrics in Eq.13, the most important are the Acc and MCC. The metrics Sn and Sp are used to measure a predictor from two completely opposite angles, and they are actually constrained with each other. Therefore, it is meaningless to use only one of the two for comparison. When, and only when, both Sn and Sp of the predictor A are higher than those of the predictor B, can we say A is better than B. In other words, a meaningful comparison in this regard should count the scores of both Sn and Sp, or even better, the rate of their combination that is none but the score of MCC.

9Now, let us use graphic analysis to further show the proposed predictor's quality. Graphs are a useful vehicle for studying complicated biological systems because they can provide intuitive insights, as demonstrated by a series of previous studies (see, e.g.,). To provide an intuitive illustration, the graph of Receiver Operating Characteristic (ROC) was adopted as given in Figure 2, where the green line is the ROC for iRNAm5C-PseDNC. The area under the ROC curve is called the AUC (area under the curve). Being within the region of 0 and 1, the greater the AUC is, the better the predictor would be. For the current predictor, the AUC is 0.9626, which is very close to 1, the value for a perfect predictor.

10As shown in a series of recent publications (see, e.g.,), papers with a user-friendly and publicly accessible web-server will significantly enhance their impacts; this is particularly true for those papers that were aimed at developing various prediction methods. In view of this, the web-server for iRNAm5C-PseDNC has been established at http://www.jci-bioinfo.cn/iRNAm5C-PseDNC. Furthermore, to maximize users’ convenience, a step-to-step guide of how to use it is given in Supplementary Information 2.

11As practiced in a series of recent studies in complying with the 5-step rule proposed in, to establish a really useful sequence-based predictor for a biological system, one should make the following five steps very clear: (1) how to construct or select a valid benchmark dataset to train and test the predictor; (2) how to formulate the biological sequence samples with an effective mathematical expression that can truly reflect their essential correlation with the target concerned; (3) how to introduce or develop a powerful algorithm (or engine) to run the prediction; (4) how to properly conduct cross-validation tests to objectively evaluate the anticipated accuracy; (5) how to provide a web-server and user guide to make people very easily to get their desired results. In the rest of this paper, we are to address these point-by-point.

12To make the description simpler and cleaner, the Chou's sequential scheme, which had been used by many previous investigators for analyzing the enzyme specificity, signal peptide cleavage sites, nitrotyrosine sites, hydroxyproline or hydroxylysine sites, methylation sites, protein-protein interaction, protein-protein binding sites, carbonylation sites, and phosphorylation sites, was adopted in this study. According to Chou's scheme, a potential RNA m5C modificationsite sample can be generally expressed by

13where the center symbol denotes the single nucleic acid code cytosine (C), the subscript is an integer, represents the -th upstream nucleotide from the center, the denotes the -th downstream nucleotide, and so forth. The -tuple RNA sample can be further classified into the following two categories:

14where represents a true sample with C at its center, a false one with C at its center, and the symbol means “a member of” in the set theory.

15In literature the benchmark dataset usually consists of a training dataset and a testing dataset: the former is for training a model, while the latter for testing it. But as elucidated in a comprehensive review, there is no need to artificially separate a benchmark dataset into the two parts if the prediction model is tested by the jackknife or subsampling (K-fold) cross-validation because 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

16where the positive and negative subsets, and, only contain the true and false m5C samples, and, respectively (see Eq.3); while denotes the symbol of “union” in the set theory.

17The benchmark dataset used in this study was derived from RMBase, which is a resource for decoding the landscape of RNA modifications from high-throughput sequencing data. The detailed procedures are as follows. (1) The genomic sequences downloaded from RMBase are in the form of DNA; to make the entire description of this paper in a coherent manner, we first change the code T to U for all the genomic sequences taken from RMBase and make them become RNA sequences. (2) As done in, by sliding the -tuple nucleotide window along each of the RNA sequences thus obtained, collected were only those RNA segments with at the center. (3) If the upstream or downstream in a RNA sequence was less than ξ or greater than where L is the length of the RNA sequence concerned, the lacking code was filled with the same code of its nearest neighbor. (4) The RNA segment samples thus obtained were put into the positive subset if their centers were experimentally annotated as the sites; otherwise, into the negative subset. (5) To reduce redundancy and bias, none of included RNA segments had pairwise sequence identity with any other in a same subset. By strictly following the above procedures, we obtained an array of benchmark datasets with different values, and hence different lengths of RNA samples as well (see Eq.2), as illustrated below

18But it was observed via preliminary tests that when (i.e., the RNA samples formed by 41 nucleotides), the corresponding results were most promising. In other words, we observed a turning point for the success rates at. After this point, the success rates would become going down with the increase of such parameter. Accordingly, hereafter we only consider the 41-tuple nucleotide samples without explicitly mentioning the parameter any more.

19The benchmark dataset thus obtained is given in Supplementary Information 1, which can also be downloaded at http://www.jci-bioinfo.cn/iRNAm5C/Supp-S1.pdf. It contains 1,900 RNA samples, of which 475 belong to the positive subset and 1,425 to the negative subset.

20An RNA samples in the aforementioned benchmark dataset can be generally expressed as

21where represents the 1st nucleotide of the RNA sample at its sequence position 1, the 2nd nucleotide at its position 2, and so forth. Except for, they can be any of the four nucleotides; i.e.,

22Based on the sequential model of Eq.6, one could directly utilize BLAST to perform statistical analysis. Unfortunately, this kind of straightforward and intuitive approach failed to work when a query RNA sample did not have significant similarity to any of the character-known RNA sequences. To overcome this problem, investigators have shifted their focus to the discrete or vector model. The reason of doing so is also due to the fact that nearly all the existing machine-learning aorithms can be directly used to handle vector models but not sequences, as elaborated in.

23One of the well-known vector models for DNA/RNA sequences is the k-tuple nucleotide (or k-mers) composition; i.e.,

24where represents the normalized occurrence frequency of the i-th k-mer, and the symbol T is the transpose operator.

25When, Eq.8 reduces to

26where,, and are the normalized occurrence frequencies of adenine, cytosine, thymine, and uracil in the RNA sequence, respectively.

27When, Eq.8 reduces to

28where is the normalized occurrence frequencies of AA in the RNA sequence, is that AC, is that AG, is that AU, and so forth.

29As we can see from above, the vector's dimension will rapidly increase with the k value, causing the so-called “high-dimension disaster” or overfitting problem. This will significantly reduce the deviation tolerance or cluster-tolerant capacity, and make the prediction model contain a lot of noise and very unstable.

30Therefore, the k-mers approach is useful only when the value of k is very small. In other words, it can only be used to incorporate the local or short-range or local sequence-order information, but certainly not the long-range or global sequence-order information. To approximately cover the long-range sequence-order effects, one popular and well-known method is to use the pseudo components that were originally introduced in dealing with protein/peptide sequences and recently extended to deal with DNA/RNA sequences.

31According to the concept of pseudo components, the RNA sequence can be generally formulated by

32where the subscript is integer and its value as well as the components will depend on how to extract the desired information from the RNA sequence of Eq.6.

33In this study, we used the approach called “physical-chemical property matrix combined with auto/cross-covariance” proposed by Liu et al. to define the components in Eq.7. According to that approach, the vector components in Eq.7 are given by

34where is an integer within the range from 0 to 39. Using exactly the same calculation approach as elaborated in, we found that was optimal choice for the current study. As for how to calculate the concrete values in Eq.12, see ref. where a crystal clear description had been given and hence there is no need to repeat here.

35Being a powerful algorithm, the random forest (RF) has been increasingly used to analyze various different problems in computational biology (see, e.g.). The essence of RF is to compare each individual classifier as a tree, and the combination of many such classifiers as a forest. In this study, 100 trees were used for the forest, and dimension of the random subspace was 22. Each tree in the forest is trained with different part of the benchmark dataset, and hence may yield a different result. The final outcome is determined via a vote from all the trees. For more information about RF, see where a very detailed description has been given, and hence there is no need to repeat here.

36The final predictor obtained via the aforementioned procedures is called as iRNAm5C-PseDNC, where “i” stands for “identify”, and “RNAm5C” for “RNA 5-methylcytosine modification sites”, and “PseDNC” for “pseudo dinucleotide composition”.

37One of the important procedures in developing a new prediction method is how to objectively evaluate its anticipated success rate. To address this, we need to consider two issues. (1) What metrics should be used to quantitatively reflect the predictor's quality? (2) What kind of test approach should be utilized to score the metrics?

38The following metrics are generally used to measure the prediction quality from four different angles: (1) Acc for measuring the overall accuracy of a predictor, (2) MCC for its stability, (3) Sn for its sensitivity, and (4) Sp for its specificity. Unfortunately, their conventional formulations as given in lack intuitiveness and most experimental scientists feel difficult to understand them, particularly for the MCC. Interestingly, using the Chou's symbols introduced in studying signal peptides, Xu et al. and Chen et al. converted them into a set of four intuitive equations, as given by

39where represents the total number of the true m5C sites investigated, while is the number of the true m5C sites incorrectly predicted to be of false m5C site; is the total number of the false m5C sites investigated, while is the number of the false m5C sites incorrectly predicted to be of true m5C site.

40According to Eq.13, it is crystal clear to see the following. When meaning none of the true m5C sites are incorrectly predicted to be of false m5C site, we have the sensitivity. When meaning that all the true m5C sites are incorrectly predicted to be of false m5C site, we have the sensitivity. Likewise, when meaning none of the false m5C sites are incorrectly predicted to be of m5C site, we have the specificity; whereas meaning that all the false m5C sites are incorrectly predicted to be of true m5C sites, we have the specificity. When meaning that none of true m5C sites in the positive dataset and none of the false m5C sites in the negative dataset are incorrectly predicted, we have the overall accuracy and; when and meaning that all the true m5C sites in the positive dataset and all the false m5C sites in the negative dataset are incorrectly predicted, we have the overall accuracy and; whereas when and we have and meaning no better than random guess. Therefore, Eq.13 has made the meanings of sensitivity, specificity, overall accuracy, and stability much more intuitive and easier-to-understand, particularly for the meaning of MCC, as concurred recently by many investigators (see, e.g.,).

41Note that, however, the set of equations defined in Eq.13 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 or biomedicine, a completely different set of metrics are needed as elaborated in.

42Now let us discuss what kind of test method should be used to score the four metrics in Eq.13. In statistical analysis, the following three methods are often used to test a predictor: (1) independent dataset test, (2) subsampling (or K-fold cross-validation) test, and (3) jackknife test. Of 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. Accordingly, the jackknife test has been widely recognized and increasingly used by investigators to examine the quality of various predictors (see, e.g.,).

43Accordingly, here we also used the jackknife test to check the quality of iRNAm5C-PseDNC predictor. During the jackknifing process, both the training dataset and testing dataset are actually open, and each sample will be in turn moved between the two. The jackknife test can exclude the “memory” effect. Also, the arbitrariness problem rooted in the independent dataset and subsampling tests can be completely avoided because the outcome obtained by the jackknife cross-validation is always unique for a given benchmark dataset.