{'fraction_non_alphanumeric': 0.044991078256436404, 'fraction_numerical': 0.025830571841277933, 'mean_word_length': 4.763712047012732, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CPPred: coding potential prediction based on the global description of RNA sequence

Abstract The rapid and accurate approach to distinguish between coding RNAs and ncRNAs has been playing a critical role in analyzing thousands of novel transcripts, which have been generated in recent years by next-generation sequencing technology. Previously developed methods CPAT, CPC2 and PLEK can distinguish coding RNAs and ncRNAs very well, but poorly distinguish between small coding RNAs and small ncRNAs. Herein, we report an approach, CPPred (coding potential prediction), which is based on SVM classifier and multiple sequence features including novel RNA features encoded by the global description. The CPPred can better distinguish not only between coding RNAs and ncRNAs, but also between small coding RNAs and small ncRNAs than the state-of-the-art methods due to the addition of the novel RNA features. A recent study proposes 1335 novel human coding RNAs from a large number of RNA-seq datasets. However, only 119 transcripts are predicted as coding RNAs by the CPPred. In fact, almost all proposed novel coding RNAs are ncRNAs (91.1%), which is consistent with previous reports. Remarkably, we also reveal that the global description of encoding features (T2, C0 and GC) plays an important role in the prediction of coding potential.

Nevertheless, the biological significance and function of most ncRNAs remain unclear comparing with coding RNAs. Rapid and accurate coding potential prediction of transcripts is critical for analyzing these data. From a computational perspective, distinguishing between coding RNAs and ncRNAs is a binary classification task and various tools have been developed (34)(35)(36)(37)(38)(39)(40)(41)(42)(43)(44)(45). In 2006, Liu et al. present an SVM-based tool (43), namely COCN, which could predict coding RNAs from ncRNAs on the basis of a hybrid feature set. However, COCN is slow in calculating abundant datasets. CPC (Coding Potential Calculator) (38) uses support vector machine to differentiate coding RNAs from ncRNAs, which is developed by Kong et al. in 2007. Six biologically meaningful features are extracted, such as ORF quality, ORF coverage, ORF integrity, sequence similarity with known proteins. However, the performance of CPC depends on the quality of multiple sequence alignment. The coding potential calculator CPC1 to CPC2 is updated in 2017 (39). The CPC2 (a predictor of coding potential based on ORF length, Fickett score, ORF integrity and isoelectric point) is much faster and more accurate than CPC1, in particular for lncRNAs. Moreover, the model of CPC2 is species-nonspecific. In 2013, Wang et al. present a logistic regression model CPAT (42) for differentiating ncR-NAs from coding RNAs, which uses four features (ORF coding RNAs on the basis of sequence features, such as ORF length, ORF coverage, ORF integrity, Fickett score, Hexamer score, Isoelectric point (pI) of a predicted peptide, Grand average of hydropathicity (Gravy) of a predicted peptide, estimation of the stability (Instability) of a predicted peptide and global descriptor (CTD) features. The CTD (composition (C), transition (T) and distribution (D)) is originally proposed for predicting protein folding class, which is global protein sequence descriptors established by Dubchak's work (48). In this work, CTD is used to denote the global transcript sequence descriptors. The CTD features include nucleotide composition, nucleotide transition and nucleotide distribution. It should be noted that the CTD features are firstly proposed by us to distinguish between coding RNAs and ncRNAs in eukaryotes. Distinguishing coding RNAs from ncRNAs based on the CTD features in prokaryote have been explored in 2009 (49); however, overlooked in recent literature. Although the nucleotide composition is as well already represented in the Fickett score, the nucleotide transition and nucleotide distribution are novel features. From Figure 4, the nucleotide transition and nucleotide distribution features (T2, C0 and GC) are important features to classify coding RNAs and ncRNAs. Besides, the mRNA secondary structure of around the start and stop codons has an important potential impact on ribosome pausing (50,51). Thus, RNA secondary structural feature contributes to predicting protein coding potential. Moreover, CTD is built to predict protein folding class primitively (48), so in this work, the CTD features are connected with RNA structural features. The CTD features not only are beneficial to prediction of coding potential (49), but also prove to be important features in RNA-binding protein prediction (52), RNAs functional identity (53) and promoter recognition (54). Subsequently, we trained the model on human dataset, and then tested it on human, mouse, zebrafish, fruit fly and Saccharomyces cerevisiae. The data of several popular species are integrated to avoid species specificity. An integrated model is built. The testing results show that CPPred with higher Matthews correlation coefficient (MCC) (55) is particularly more effective on sORF data when compared with other programs. Moreover, to highlight the performance of CP-Pred on the sORF data, we compared it with sORF finder (56). The sORF finder is developed by Hanada et al. and using nucleotide composition to identify sORF. It is proved by Cheng et al. to be superior to other tools for predicting sORF (57). Additionally, CPPred is convenient because it only needs FASTA format sequence files as inputs.

Dataset
Two models are built for distinguishing between coding RNAs and ncRNAs. The first one (Human-Model) is trained by human data, and then tested on human, mouse, fruit fly, zebrafish and S. cerevisiae. The second one (Integrated-Model) is built for integrated species, including human (Homo sapiens), mouse (Mus musculus), zebrafish (Danio rerio), fruit fly (Drosophila melanogaster), S. cerevisiae, nematode (Caenorhabditis elegans) and thale cress (Arabidopsis thaliana). As the former, we downloaded  For ncRNAs, the data that have no source comments and are not annotated with Havana in the corresponding of gff3 file are removed. After that, the number of coding RNAs and ncRNAs is 50 040 and 36 244, respectively. We randomly selected two-thirds of the data as training set, a collection of 33 360 coding RNAs and 24 163 ncRNAs, which is called Human-Training. Then, the rest of the data are stored as a testing set. At the same time, we reduced redundancy between the testing and training set using CD-hit with sequence identity cutoff ≥80%. Finally, 8557 coding RNAs and 8241 ncRNAs are kept as Human-Testing. Then, the sequences with ORF shorter than 303 nucleotides in length are extracted from coding RNA in Human-Testing. Meanwhile, the same amount of considerable length ncRNAs from Human-Testing are selected randomly.  Figure 2). The transcript status of coding RNAs is 'KNOWN'. However, for ncRNAs, we removed 1053 ncRNA sequences without source comments, and these sequences are not annotated with Havana in the corresponding gff3 file. Eventually, the total numbers of annotated coding RNAs and ncRNAs are 50 040 and 36 244, respectively. We randomly selected two-thirds (41) of data as a training set, that is, the set including 33 360 coding RNAs and 24 163 ncRNAs, which is called Human-Training. Then, the remaining data are stored as a testing set. To improve the robustness of the assessment of accuracy, we reduced redundancy between the testing and training sets with a threshold of 80% (62)(63)(64) by the open-source program cd-hit-est-2d in CD-hit (65), which uses a short word filter to avoid unnecessary alignments and has been widely used as a clustering algorithm, resulting in 8557 coding RNAs and 8241 ncRNAs as Human-testing. Besides, the basic command of cd-hit-est-2d is cd-hit-est-2d -i training human -i2 testing human -o testing human redundancy -c 0.8 -n 5. Meanwhile, we examined the different threshold (<80%) for removing the redundancy much more stringently between the testing and training sets. The redundancy between the testing and training sets is removed by program cd-hit-2d in CD-hit with threshold of 75%, 70% and 60%, respectively. The performance of CPPred is shown in Supplementary Table S10 on the testing sets. The last row in Supplementary Table S10 is the performance of CPPred in the manuscript. The result shows that removing redundant sequences with ≥80% similarity is enough to obtain independence between test and training sets. Besides, the BLASTCLUST is also tried and used to exclude redundant RNA sequences with the similarity threshold of 75%, 70%, 60%, 50%, 40%, 30%, 20% and 10%, respectively. The result shows the same conclusion. That is because RNA sequences only have four letters, it is easy to get high identity between unrelated sequences. When we run a regular blast search, we often find top hits at very high identities. In general, comparing RNAs at low identity <75% may not be very effective. The threshold of 80% is a stringent cutoff. For example, the sequence redundancy in the datasets of CPC2 (39), lncRScan-SVM (45) and lncADeep (66) is removed by using CD-hit with thresholds of 0.9, 0.8 and 0.95, respectively. The sequence redundancy of CONC's datasets (43) is removed by using NCBI BLASTCLUST with the '-L 0.7' option, and the sequence similarity of CPC's  (37) and mRNN (71), no redundancy is removed. For the reason that the prediction of sORF is difficult, the sequences with ORF fragments <303 nucleotides in length from coding RNAs in Human-Testing are selected to build an interesting and challenging testing set. Meanwhile, the same amount of comparable length ncRNAs from Human-Testing were filtered out randomly (34,39). As a result, 641 coding RNAs and 641 ncRNAs are kept as Human-sORF-Testing. Based on the same building method of the human testing sets, we constructed Mouse-Testing, Mouse-sORF-Testing, Zebrafish-Testing, Zebrafish-sORF-Testing, S. cerevisiae-Testing, S. cerevisiae-sORF-Testing, Fruit-fly-Testing and Fruit-fly-sORF-Testing. These datasets are downloaded from Ensembl, as shown in Table 1. Afterward, several popular species are integrated for the sake of eliminating the problems caused by the specificity of species and the differences between the databases. The data for human, mouse, zebrafish, fruit fly, S. cerevisiae, nematode and thale cress are downloaded from NCBI Ref-Seq including 525 316 coding RNAs and 55 198 ncR-NAs. Furthermore, to reduce the computational effort and balance the proportion of human negative-positive data, we randomly selected 52 530 coding RNAs and 27 600 ncRNAs as training set with the same percentage of each species, which is called Integrated-Training. The redundancy of the remaining is removed with sequence identity cutoff ≥80% (62-64). The Integrated-Testing (balance data) is constructed with 13 903 coding RNAs and 13 903 ncR-NAs (34,39). As the same building procedure with Human-sORF-Testing, the Integrated-sORF-Testing is obtained, which contains 11 634 small coding RNAs and 11 634 small ncRNAs.

CPPred features
To predict the coding potential of RNA sequences, we extracted features from recently published scientific literature (39,42), and novel CTD features are added.
We used four features proposed by CPAT (42), including ORF length, ORF coverage, Fickett score and Hexamer Score. Similarly, ORF integrity and isoelectric point were derived from CPC2 (39). Next, the Gravy and Instability mentioned by CPC2 are also added. In addition, the algorithm of Hexamer score and Fickett score were discussed in detail (42,72,73). The Fickett score is calculated by con-sidering eight properties of coding sequences. Four of them are composition values with the frequencies of the four nucleotides from RNA sequence. The other four parameters are position value, which reflect the degree to codon preference. Furthermore, the feature of ORF length is used to predict the coding potential (38,39,42,43), but it relies on the full-length transcript (41). In this study, although the ORF length requires annotation of the full-length transcript, we still selected it because of its identifiable power and ease of calculation. The ORF length is the length of the maximum open reading frame, which starts with a start codon and ends with a stop codon (UGA, UAA or UAG). Here, the AUG is selected as the start codon. Although the use of non-AUG has been constantly described (26,29,30,33,(74)(75)(76), there is no clear consensus on how to choose translation start sites (77).
In particular, we added 30 new features, which were CTD features. In this study, CTD is used to denote the global transcript sequence descriptors. The transcript is a sequence containing four types of nucleotides A, T, G and C. The nucleotide composition (first index C) describes the percent composition of each nucleotide in a transcript sequence, which is contained in Fickett score. The nucleotide transition (second descriptor T) describes the percent frequency with conversion of four nucleotides between adjacent positions. Subsequently, we calculated five relative positions along the transcript sequence of each nucleotide, with the 0 (first one), 25%, 50%, 75% and 100% (last one), to describe the nucleotide distribution (last descriptor D).
For example, the RNA sequence is ACTTGCAGCC CCCCGCCTGTCCCGAG CCGCGCGGGCGCCAGC TCAGTTTGTCCGCGGCGG, which contains 5 adenines (As), 9 thymines (Ts), 20 guanines (Gs) and 26 cytidines (Cs). The features of first descriptor C are 5/60 = 0.083, 9/60 = 0.15, 20/60 = 0.33 and 26/60 = 0.43, respectively. We use A, T, G and C to represent the four features. For the second descriptor T, there is zero transition between A and T, five transitions between A and G, four transitions between A and C, five transitions between T and G, six transitions between T and C and twenty transitions between G and C. Therefore

Feature selection
The 38 features as mentioned may include redundant features, so a feature selection process is used to filter out redundant features for coding potential prediction. In this work, mRMR-IFS method is used to select the best subset of features (78)(79)(80)(81). The mRMR program is developed by Peng et al. (82), which selects good features based on PAGE 5 OF 10 Nucleic Acids Research, 2019, Vol. 47,No. 8 e43 mutual information with the minimal redundancy, maximal relevance criteria. First, the mRMR program is used to rank the 38 features in the training set. Then, Incremental Feature Selection (IFS) (78) is used to increase the features one by one in descending based on the mRMR ordering. For each additional feature, a new subset of feature is generated. Therefore, a total of 38 feature subsets are generated for 38 sorting features.
According to the obtained 38 feature subsets, we select the corresponding feature sets from the whole training set. Through 10-fold cross-validation on the training set, the best subset of features is selected and served as the final model.

SVM classifier
We used the Libsvm (83)   The ORF coverage is defined as the ratio of ORF to the length of a transcript. The ORF length, Hexamer score and Fickett score are discussed in the 'CPPred features' section. The integrity of the ORF is defined as whether the ORF starts with a start codon (AUG) and ends with a stop codon (UGA, UAA or UAG). PI, Gravy and Instability are calculated by the ProtParam. After that, using mRMR-IFS, the best feature subset is selected and used as input to the SVM classifier. Eventually, we got the final model, which is tested and evaluated by the testing sets.

Pipeline of CPPred
In CPPred, an SVM model is applied to calculating coding potential of a transcript by using features derived from RNA and protein sequences (see Figure 3), which is designed for distinguishing between coding RNAs and ncR-NAs. First, we constructed a training dataset, which contained coding RNAs and ncRNAs. Then, 38 features are calculated for each RNA or protein sequence. The uninformative features are reduced by using mRMR-IFS and the best feature subset is picked out. Based on the feature subset, the SVM classifier is used to obtain a model on the training set. Finally, the CPPred is analyzed and evaluated on testing sets.

Feature selection by the mRMR-IFS method
Here, mRMR-IFS (78) method is chosen for feature selection. For each feature subset, the corresponding features are selected and 10-fold cross-validation is performed on the training set. In Figure 4, for Human-Training, the best predictive performance is obtained by using top 37   In addition, the value of MCC dramatically increased to 0.920 on the Human-Training in Figure 4 after adding the seventh feature and then the MCC tended to be stable. The top seven features include ORF integrity, ORF Coverage, Instability, T2 (CTD feature), C0 (CTD feature), isoelectric point and ORF length (Supplementary Table S7). In particular, the seventh feature is the length of ORF. At the same time, for the Integrated-Training, the MCC value (MCC = 0.925) increased suddenly after the fifth feature was added. The top five features contain ORF coverage, ORF integrity, GC (CTD feature), Instability and ORF length (Supplementary Table S8). Moreover, the fifth feature is also the length of ORF. In conclusion, the length of the ORF is a more crucial feature for distinguishing between coding RNAs and ncRNAs, which is consistent with the findings by Wang et al. (42).
The top seven features of Human-Training and the top five features of Integrated-Training are important to differentiate coding RNAs from ncRNAs. Among these important features, ORF coverage and ORF length are from CPAT (42). The isoelectric point, length and integrity of the ORF are from CPC2 (39). Of the remaining important features, T2, C0 and GC are derived from CTD features, which indicate the importance of CTD features for distinguishing between coding RNAs and ncRNA. Besides, in Figure 4, it also reveals that the ORF length, ORF coverage, ORF integrity and Instability are shared among species, while pI, T2 and C0 are human-specific.

Performance of CPPred (Human-Model)
To evaluate our method CPPred, for Human-Model, we compared CPPred with CPAT (42), CPC2 (39) and PLEK (34) on testing sets of human, mouse, zebrafish, S. cerevisiae and fruit fly. Moreover, for sORF testing sets of human, mouse, zebrafish, S. cerevisiae and fruit fly, the method sORF finder (56) has also been added for comparison. The results of Human-Testing are shown in Table 2 Table 3. Significantly, CPPred is much better than the other three methods (CPAT, CPC2 and sORF finder) with MCC of 0.654 versus 0.472, 0.125 and 0.262. Overall, from the two testing sets, our method performs better than CPAT and CPC2, however, slightly worse than PLEK, probably due to the high level of redundancy between the two human testing sets and PLEK's training set (84). Herein, we downloaded the human training set of PLEK from https://sourceforge.net/ projects/plek/. Subsequently, we used the CD-HIT tool with a value of 0.8 as the threshold of sequence identity to analyze the redundancy of between human testing set of CP-Pred and human training set of PLEK. The results show that 3944 (3944/8557 = 46.1%) coding RNAs and 5411 (5411/8241 = 65.7%) ncRNAs are redundant sequences between the human testing set of CPPred and human training set of PLEK, respectively. Besides, the predictive performance of PLEK dropped behind the other methods (CP-Pred, CPAT, CPC2 and sORF finder) when testing with mouse, zebrafish, S. cerevisiae and fruit fly (Tables 4 and  5; Supplementary Tables S1-6). Moreover, from Table 4 Table 5.
In addition, for the testing sets of zebrafish and S. cerevisiae, CPPred also achieved the best performance. The MCCs of CPPred in the Zebrafish-Testing, Zebrafish-sORF-Testing, S. cerevisiae-Testing and S. cerevisiae-sORF-Testing are 0.9, 0.387, 0.9 and 0.440, respectively (Supplementary Tables S1-4 It is noteworthy that, in most cases, the CPPred in the testing of sORF is much better than CPAT, CPC2, PLEK and sORF finder.

Performance of CPPred (Integrated-Model)
Due to the specific differences among the species, we tested Integrated-Model in the Integrated-Testing and Integrated-sORF-Testing, and the results are presented in Tables 6 and  7, respectively. CPPred performed better than CPAT (42), CPC2 (39), PLEK (34) and sORF finder (56) (see MCC, AUC and ACC in Tables 6 and 7). Noteworthy, the model of CPC2 is species-neutral (39), and the performance of CPC2

OCTD features
To highlight the importance of CTD features, we built a model (OCTD-Model) by using only 30 CTD features (OCTD). The OCTD-Model is trained by Human-Training and then tested on Human-Testing and Human-sORF-Testing. As compared with OCTD-Model, we obtained NCTD-Model on Human-Training by using non-CTD features (i.e. ORF length, ORF coverage, Fickett score, Hexamer Score, Gravy, pI, ORF integrity and Instability) and then also tested on Human-Testing and Human-sORF-Testing. The results are shown as Supplementary  sORF-Testing (Table 5), which may be the main reason for a higher accuracy rate predicted by PLEK than CPPred. Moreover, CPPred performed worse than sORF finder. The reason may be that the sORF finder has a higher false positive rate due to base on a single feature (Tables 3, 5 and 7;  Supplementary Tables S2 and S6).

CONCLUSION
In this work, based on SVM classifier algorithm, we developed a tool CPPred to predict coding potential using multiple features, which are extracted from CPAT (42) and CPC2 (39), and CTD features are added particularly. Here, we used CTD features to predict coding potential for the first time in eukaryotes. Moreover, we found that the features of T2, C0 and GC (CTD features) play a key role in predicting coding potential. Our method CPPred is trained on Human-Training and Integrated-Training to obtain the Human-Model and the Integrated-Model. As the former, CPPred is tested on the dataset of human, mouse, zebrafish, S. cerevisiae and fruit fly, obtaining AUC from 0.72 to 0.99. Besides, our CP-Pred is compared with other methods CPAT, CPC2, PLEK and sORF finder. The CPPred outperforms CPAT, CPC2, PLEK and sORF finder on the testing sets of mouse, zebrafish and S. cerevisiae. However, CPPred does not perform as well as PLEK on human testing sets, which may be due to the high level of redundancy between the testing set of human and PLEK's training set (84) (Tables 4 and 5). We analyzed the testing set of human and PLEK's training set, and found 46.1% coding RNAs and 65.7% ncR-NAs redundant sequences between the human testing set of CPPred and human training set of PLEK, respectively. For fruit fly, CPPred performs worse than CPAT, which may be due to the fact that CPAT is trained on fruit fly while CP-Pred is trained on human dataset. Thus, the second model is built by Integrated-Training. From Tables 6 and 7, we compared CPPred with other tools (CPAT, CPC2, PLEK and sORF finder) and found some improvement in MCCs by >5% and >11% on Integrated-Testing and Integrates-sORF-Testing, respectively. Moreover, the CTD features are particularly important for predicting the coding potential of sORF datasets (Supplementary Table S9). Overall, the results demonstrate that CPPred performs well on long RNA datasets and much better than other tools on sORF datasets.

DATA AVAILABILITY
Source code was implemented in Python and is freely available at http://www.rnabinding.com/CPPred.