M6AMRFS: Robust Prediction of N6-Methyladenosine Sites With Sequence-Based Features in Multiple Species

As one of the well-studied RNA methylation modifications, N6-methyladenosine (m6A) plays important roles in various biological progresses, such as RNA splicing and degradation, etc. Identification of m6A sites is fundamentally important for better understanding of their functional mechanisms. Recently, machine learning based prediction methods have emerged as an effective approach for fast and accurate identification of m6A sites. In this paper, we proposed “M6AMRFS”, a new machine learning based predictor for the identification of m6A sites. In this predictor, we exploited a new feature representation algorithm to encode RNA sequences with two feature descriptors (dinucleotide binary encoding and Local position-specific dinucleotide frequency), and used the F-score algorithm combined with SFS (Sequential Forward Search) to enhance the feature representation ability. To predict m6A sites, we employed the eXtreme Gradient Boosting (XGBoost) algorithm to build a predictive model. Benchmarking results showed that the proposed predictor is competitive with the state-of-the art predictors. Importantly, robust predictions for multiple species by our predictor demonstrate that our predictive models have strong generalization ability. To the best of our knowledge, M6AMRFS is the first tool that can be used for the identification of m6A sites in multiple species. To facilitate the use of our predictor, we have established a user-friendly webserver with the implementation of M6AMRFS, which is currently available in http://server.malab.cn/M6AMRFS/. We anticipate that it will be a useful tool for the relevant research of m6A sites.


INTRODUCTION
To date, more than 150 types of RNA modifications have been discovered (Maden, 1990;. Of these modifications, N6-methyladenosine (m 6 A) is the most common and abundant one and exists in various species. It is found to be closely associated with diverse biological processes, such as RNA localization and degradation , RNA structural dynamics (Roost et al., 2015), alternative splicing (Liu N. et al., 2015), primary microRNA processing (Alarcón et al., 2015), cell differentiation, and reprogramming , and regulation of circadian clock (Geula et al., 2015). Thus, identification of m 6 A sites is of great importance for better understanding of their functional mechanisms. In the past few years, high-throughput experimental methods, such as MERIP (Meyer et al., 2012) and m 6 A-seq (Dominissini et al., 2012), have been utilized to identify m 6 A modifications, and more and more m 6 A peaks have been characterized. However, they have the following limitations: (1) they cannot accurately locate the positions of m 6 A sites; (2) they are highly cost; and (3) they are not applicable for the large-scale identification of m 6 A sites. Hence, it is highly desirable to develop fast and accurate computational methods for the identification of m 6 A sites (Chen et al., 2015b. In recent years, machine learning based prediction methods have emerged as effective approach for predicting m 6 A sites. For example, Chen et al. (2015a) developed the first machine learning based predictor, called "iRNA-Methyl", for m 6 A site identification. They exploited physicochemical properties and sequence-order information embedded in PseDNC (pseudo dinucleotide composition) (Liu B. et al., 2015), and used support vector machine for model construction. Later, Liu Z. et al. (2016) proposed to incorporate more additional physicochemical properties coupled with a scalable transformation algorithm into their feature extraction model. To improve the predictive performance, Jia et al. proposed to fuse three types of feature descriptors, such as bi-profile Bayes, dinucleotide composition and KNN scores. Their results showed that this fusion strategy is able to achieve better performance than single one feature descriptor . Similarly, Xiang et al. (2016b) found that combining binary encoding scheme together with k-mer frequency could contribute to the improved performance. Recently, Zhou et al. (2016) developed "SRAMP", a powerful prediction tool using multiple types of feature descriptors, including positional binary encoding of nucleotide sequence, k-nearest neighbor encoding, nucleotide pair spectrum encoding, and secondary structure pattern, to train an ensemble predictive model with random forest for the identification of m 6 A sites. SRAMP is reported to achieve relatively good performance as compared to other predictors. More recently, Xiang et al. (2016a) proposed a new predictor called "RNAMethyPre", using compositional information and position-specific information to build predictive models for the prediction of m 6 A sites on both human and mouse. Additionally, in our previous study, we proposed to use deep learning algorithm to generate highlatent features to improve the predictive performance (Wei et al., 2018d). However, we found that most of existing predictors are species-specific. Currently, there is not any predictor that is capable of predicting m 6 A sites for multiple species.
For this purpose, we proposed a novel sequence-based predictor, namely "M6AMRFS" for detecting m 6 A sites in RNA sequences. For feature extraction (Mrozek et al., 2007(Mrozek et al., , 2013, we proposed a feature representation algorithm to encode sequences with dinucleotide binary encoding and local positionspecific dinucleotide frequency. To optimize the feature space, we combined the F-score algorithm with SFS (Sequential Forward Search) (Wei et al., 2018a,c,e) to improve the representation ability of our features. For model training, we trained the optimal feature representations under XGBoost algorithm. Our experimental results showed that the proposed M6AMRFS is able to achieve competitive and robust performance as compared to state-of-the-art predictors for four different species. To the best of our knowledge, this is the first predictor that is applicable for multiple species. Furthermore, we have established a user-friendly webserver that implements the proposed M6AMRFS, which is currently available in http:// server.malab.cn/M6AMRFS/. We anticipate that it will be a useful tool complementary for existing tools, facilitating to further reveal the functional mechanisms of m 6 A sites.

Benchmark Datasets
To predict the m 6 A sites in multiple species, we employed four benchmark datasets from four species, including Saccharomyces cerevisiae, Arabidopsis thaliana, Musculus, and Homo sapiens. The detail of the four benchmark datasets is listed in Table 1. For the four benchmark datasets, the positives are the sequences centered with true m 6 A sites, while the negatives are usually the sequences centered with adenines but without any m 6 A peaks detected. The datasets can be found in the following website: http://server.malab.cn/M6AMRFS/. Figure 1 illustrates the overall procedure of the proposed predictor. As we can see from Figure 1, there are two steps in the predictor. The first step is data pre-processing, including data clean and feature extraction. It filters out those irrelevant sequences from input sequences. Then, the resulting sequences are submitted into the feature representation algorithm, in which the sequences are encoded with feature vectors. The second step is feature optimization and model training. For feature space optimization, we used the F-score algorithm combined with SFS (Sequential Forward Search) to search for the optimal features. Afterward, the resulting optimal feature representations are fed into a well-trained XGBoost model to predict whether the sequences are true m 6 A sites or not. In our predictor, the predicted outcome for each sequence is 0 or 1, where 0 denotes non-m 6 A site and 1 denotes true m 6 A site.

Feature Representation
In this work, we present a new feature representation algorithm that combines two feature descriptors. One is named "Dinucleotide binary encoding" and the other is "Local positionspecific dinucleotide frequency", which are described as follows,

Dinucleotide Binary Encoding
The feature descriptor encapsulates the positional information of the dinucleotide at each position in the sequence. Obviously, there are a total of 16 possible dinucleotides. In this descriptor, each dinucleotide can be encoded into a 4-dimensional 0/1 vector. For example, AA is encoded as (0,0,0,0); AT is encoded as FIGURE 1 | Framework of algorithms proposed in this study. There are two main steps. In the first step, the input RNA sequences are filtered by removing those irrelevant sequences. Then, the remaining sequences are fed into the proposed feature extraction algorithm for feature representation. In the second step, the resulting feature representations are optimized by feature selection, and then, the optimal feature representations are predicted by a XGBoost model.

Local Position-Specific Dinucleotide Frequency
For a given sequence, the feature vector of this descriptor can be denoted as (f 2 , f 3 , . . ., f l ), where f i is calculated as follows, where l is the length of the given sequence, |N i | is the length of the i th prefix string {X 1 X 2... X i } in the sequence, and C (X i−1 X i ) is the occurrence number of the dinucleotide X i−1 X i in position i of the i th prefix string.

Feature Selection
Feature selection is an important process to improve the classification performance (Mrozek et al., 2009;Mrozek et al., 2014;Zeng et al., 2016;Zou et al., 2016a,b;Liu, 2017). Here, we used the F-score algorithm together with the SFS strategy to search the most discriminative features (Peng et al., 2005). Figure 2 illustrates the procedure of the feature selection strategy, which is described as follows. Firstly, the F-score algorithm is utilized to rank all the features from the highest scores to the lowest scores, generating a ranked feature list. Secondly, we added the features one by one from the ranked list, and respectively trained the predictive models. Lastly, the feature subset corresponding to the highest accuracy of the predictive model is used as the optimal features. The results of feature selection were discussed in section of "Results and Discussion".

XGBoost (eXtreme Gradient Boosting)
eXtreme Gradient Boosting, which was proposed by Chen and Guestrin (2016), has been shown to be a powerful classification algorithm. The general idea of XGBoost is to enumerate several candidates that may be the segmentation points according to the percentile method, and then to find the best segmentation point from the candidates for calculating the segmentation points.
The main advantage of XGBoost is to combine multithreading, data compression, and fragmentation methods to improve the efficiency of the algorithm as much as possible. Moreover, the regularization terms added by XGBoost in the loss function can be used to control the complexity of the model and avoid overfitting. Parameters, such as subsamples, max depth, and estimators, are utilized to optimize evaluation performance via parallelization program namely "Grid Search". For the implementation of XGBoost in our predictor, the range of max depth is set from 2 to 10; learning rate is ranged from 0.1 to 0.8; and estimators are ranged from 1 to 10.

Performance Evaluation
In this work, four commonly used performance metrics are used for performance evaluation, including Acc (accuracy), Sn (sensitivity), Sp (specificity), and MCC (Mathew's correlation coefficient), respectively (Zeng et al., 2015;Lai et al., 2017;Zhang et al., 2017;Cheng et al., 2018;Su et al., 2018;Tang et al., 2018;Wei et al., 2018b;Yang et al., 2018). They are formulated as follows where TP denotes true positive; TN denotes true negative; FP denotes false positive; and FN denotes false negative. Sn measures the predictive ability of a predictor for positive samples while Sp measures the predictive ability of a predictor for negative samples. Acc and MCC are two metrics measuring the overall performance of a predictor. Besides, we used Receiver Operating Characteristic (ROC) curve to intuitively evaluate the overall performance (Liu et al., 2013. It is plotted with true positive rate (TPR) against false positive rate (FPR) under different classification thresholds. The TPR is the same with sensitivity as described above, while FPR is calculated as 1-specificity. Area under ROC curve (AUC) is usually used as an evaluation metric . The value of AUC ranges from 0.5 to 1. If the AUC is close to 1, it indicates that the predictor has excellent performance. If the AUC approaches to 0.5, the predictor does not perform well for prediction.
Additionally, we used 10-fold cross validation method and jackknife test to evaluate the predictive performance (Wei et al., 2017a;Zeng et al., 2017a,b;Liao et al., 2018;Zou et al., 2018). The two evaluation methods were chosen since existing methods in the literature used them for performance evaluation.

Comparison of XGBoost and Other Classifiers
To evaluate the effectiveness of the XGBoost classifier, we compared it with five commonly used machine learning algorithms, including Random Forest (RF) (Liu B. et al., 2015;Li et al., 2016;Wei et al., 2017b), Naïve Bayes (NB), Logistic Regression (LR), K-Nearest Neighbors (KNN) (Huang and Li, 2018), Support Vector Machine (SVM) (Song et al., 2010(Song et al., , 2012(Song et al., , 2018Wang M. et al., 2014;Wei et al., 2017), and Gradient Boosting Decision Tree (GBDT) (Liao et al., 2018), respectively. For fair comparison, the machine learning algorithms were trained and evaluated with 10-fold cross validation on the benchmark datasets, respectively. The performance of different classifiers is illustrated in Figure 3. The detailed results are presented in Table 2.
As shown in Table 2 and Figure 3, XGBoost outperforms the other classifiers on three out of the four datasets, with the exception of Dataset-A101, for which the SVM classifier is slightly better than the XGBoost, which is the second best among the compared classifiers. For those datasets that the XGBoost outperforms other classifiers, the XGBoost is able to achieve higher Acc and MCC. To be specific, our Acc and MCC are 0.7314 and 0.4629 in the Dataset-S51, 0.6 and 1.1% higher

Impact of Feature Selection
In this study, we employed the F-score with the SFS for feature selection. The results of feature selection are summarized in Table 3 and illustrated in Figure 4 as well. As seen from Table 3,  terms of Acc, Sn, Sp, and MCC, respectively. After applying the feature selection, we observed that the performances in terms of all the metrics were improved. To be specific, the Acc and MCC were improved to 0.7425 and 0.4852, respectively. This indicates that the feature selection strategy to yield more informative features to distinguish true m 6 A sites from nonm 6 A sites. For the other datasets from different species, similar results were observed. We can see from Table 3 that almost all the performances were improved by using feature selection, demonstrating that feature selection is an effective way to enhance the predictive performance of the predictor. Moreover, Figure 4 illustrates the Acc of the features by varying the feature number when conducting feature selection. As seen in Figure 4, we pointed out the optimal feature number and their corresponding highest Acc for each dataset. The optimal feature number for the four datasets are 85, 57, 13, and 355, giving the highest Acc of 0.7425, 0.9102, 0.8924, and 0.8105, respectively.

Comparison With Other Feature Representation Algorithms
To examine the performance of the proposed feature algorithm, we evaluated and compared it with existing feature representation algorithms, including RFH, PseDNC, PCP (physical and chemical properties), KNN (K-Nearest Neighbors), and AthMethPre, respectively. These algorithms were reported to have relatively strong power for the identification of m 6 A sites. Thus, they were chosen for comparison. The results of the above algorithms were presented in Table 4. As we can see from Table 4, the proposed features are competitive with the best-performing AthMethPre other feature representation methods and remarkably outperform the other existing features in all the four datasets. Note that for the Dataset-S51 and the Dataset-A101, our method performs slightly worse than the best-performing AthMethPre; while for the other two datasets, our method is slightly better. As well known, for the genome-wide identification, the running time for a predictor is important as well. Therefore, we further compared the feature number of AthMethPre and our feature representation method. We found that the feature number of the AthMethPre method for each dataset are 540, 500, 500, and 740, while ours are 85, 57, 13, and 355, respectively. As can be seen, our feature numbers for all the four datasets are averagely much fewer than the AthMethPre method. This indicates that the computation time by our predictive models costs less. In general, it can be concluded that our features are at least effective for the representatives of m 6 A sites in multiple species with different sequence lengths.

Comparison With State-of-the-Art Predictors
To assess the effectiveness of our predictor, we compared it with existing predictors including pRNAm-PC (Liu Z. et al., 2016), MehtyRNA , and RFAthM6A (Wang and Yan, 2018), respectively. There were chosen since they were reported to have the best performance on the four benchmark datasets used in this work. The results were presented in Table 5. As shown in Table 5, M6AMRFS outperforms pRNAm-PC on the Dataset-S51. The Acc, Sn, Sp, and MCC by our predictor are 0.7425, 0.7521, 0.7339, and 0.4852, respectively. The performances are higher than that of the second best pRNAm-PC on this dataset. To be specific, our overall performances are 0.0451 and 0.0852 higher in terms of Acc and MCC, respectively. As for the other datasets (Dataset-H41 and Dataset-M41), we observed similar results that our overall performance outperforms the existing predictors. Only on Dataset-A101, our predictor performs slightly worse than RFAthM6A. To be concluded, our results demonstrate that the proposed predictor is better than existing predictors or at least competitive with existing predictors on multiple benchmark datasets from different species. Importantly, our predictor exhibits robust performance for multiple species, demonstrating that our predictor is able of capturing the characteristics of m 6 A sites in different species. This also implies that the m 6 A sites from different species might share the common patterns.

CONCLUSION
In this study, we have developed a machine learning based predictor, namely M6AMRFS, for the identification of m 6 A sites in multiple species. We have conducted a series of comparative study, and our experimental results indicate that our predictor is at least competitive as compared to previously published predictors. Importantly, we found that our predictor is able to achieve robust performance in several species. To the best of our knowledge, it is the first predictor that can provide predictions in multiple species. For further analysis, we found that the robust performance contributes to the following two possible reasons. One reason is the XGBoost classifier we used for model training. We have compared XGBoost with other machine learning algorithms. XGBoost is shown to perform better than other classification algorithms. The other reason is that our feature selection strategy helps to adaptively select the optimal features for specific species. We anticipate that the tool and webserver we have established will be useful for facilitating to reveal the functional mechanisms of m 6 A sites.

AUTHOR CONTRIBUTIONS
XQ and HC wrote the manuscript. HC developed the webserver and analyzed the results. XY analyzed the results. RS and LW designed the experiments. All authors read and approved the manuscript.

FUNDING
The work was supported by the National Natural Science Foundation of China (Nos. 61701340 and 61702361).