Porpoise: a new approach for accurate prediction of RNA pseudouridine sites

Overall framework of Porpoise

The design and performance evaluation process of Porpoise is summarized in Figure 1. As can be seen, there exist four major steps, which include dataset cleaning/preparation, feature engineering, stacked model training and optimization and performance evaluation. At the first step, we collect the benchmark training and independent test datasets from the related literature and online databases [11]. At the second step, we comprehensively extract 18 types of RNA sequence-based features and evaluate each type of features coupled with nine commonly used machine learning algorithms. The scatter plot in the “feature engineering” stage in Figure 1 is the t-SNE plot of the pseudouridine sites of three different species based on the PSTNPss features. At the third step, we construct a series of stacked ensemble learning models through different combinations of base classifiers and optimize the model for each of the three species. In the fourth step, we comprehensively evaluate the optimized stacked models by performing cross-validation and independent tests against several existing state-of-the-art approaches. Finally, we implement and make publicly available an online webserver of Porpoise.

Data collection

Current predictors for RNA pseudouridine sites, including iRNA-PseU [11], PseUI [12], iPseU-CNN [13], XG-PseU [14] and EnsemPseU [15], are all trained and evaluated using the same datasets, previously constructed in the iRNA-PseU work [11]. These datasets were collected from the RMBase database [21], which included three benchmark datasets, named H_990 (H. sapiens), S_628 (S. cerevisiae) and M_944 (M. musculus), for model training; and two other independent test datasets, named H_200 and S_200, for performance validation and comparison between different methods. In this study, we use the same datasets to construct the models of Porpoise and benchmarked the performance of Porpoise against the state-of-the-art approaches. All the samples in both H. sapiens datasets (H_990 and H_200) and the M. musculus dataset (M_944) have the RNA sequences with 21 nucleotides with uridine at the center. In comparison, the samples in S. cerevisiae (S_628 and S_200) have the RNA sequences with 31 nucleotides with uridine at the center. Both training and independent test datasets are balanced with the ratio of positive and negative samples set as 1:1. For the training datasets, H_990, S_628 and M_944 contained 495, 314 and 472 positive samples and the same number of negative samples. For independent test datasets, both H_200 and S_200 encompassed 100 positive and 100 negative samples, respectively. In addition, to further evaluate the performance of Porpoise and benchmark it against other state-of-the-art methods, we collect another independent test dataset from the m6A-Atlas database [22], a comprehensive public database for multiple RNA modifications. This independent test dataset includes 3,137 H. sapiens, 2,702 M. musculus and 733 S. cerevisiae pseudouridine sites and is named as Psi_test.

Feature engineering

In this study, we comprehensively test 18 types of RNA sequence encoding schemes to identify the best feature encoding combinations for training the stacked models, including binary feature, autocorrelation (AC), cross-covariance (CC), trinucleotide-based auto-cross covariance (TACC), accumulated nucleotide frequency (ANF), nucleic acid composition (NAC), pseudo nucleic acid composition (PseDNC), pseudo k-tuple composition (PseKNC), the composition of k-spaced nucleic acid pairs (CKSNAP), di-nucleotide composition (DNC), tri-nucleotide composition (TNC), electron-ion interaction pseudopotentials of trinucleotide (EIIP), enhanced nucleic acid composition (ENAC), nucleotide chemical property (NCP), electron-ion interaction pseudopotentials of trinucleotide (PseEIIP), position-specific trinucleotide propensity based on double-strand (PSTNPds), position-specific trinucleotide propensity based on single-strand (PSTNPss) and reverse complement Kmer (RCKmer). All these sequence encoding schemes can be calculated using our developed open-source iLearn and iLearnPlus software packages [23, 24]. Finally, four specific encoding schemes, including binary features, PseKNC, NCP and PSTNPss, are selected according to their predictive performance and used for training the stacked models. These four encoding schemes are described in the following sections.

Pseudo k-tuple composition (PseKNC)

The Pseudo k-tuple composition (PseKNC) feature is a type of pseudo nucleic acid composition feature that considers both the local and long-range sequence information [23]. PseKNC combines the k-tuple nucleotide composition, which is defined as:

(1)

where

(2)

where is the frequency of oligonucleotide, which is normalized to , is the factor and is defined as:

(3)

in Equation (3) is the correlation function, which is defined as follows:

(4)

where is the number of physicochemical indices. The PseKNC encoding used six indices, including ‘Rise (RNA)’, ‘Roll (RNA)’, ‘Shift (RNA)’, ‘Slide (RNA)’, ‘Tilt (RNA)’ and ‘Twist (RNA)’ for RNA sequences. is the numerical value of the -th () physicochemical index of the dinucleotide at the position .

Position-specific trinucleotide propensity based on single strand

The Position-specific trinucleotide propensity based on single strand (PSTNPss) encoding scheme employs a statistical rule to encode the RNA sequences. In general, there are 64 (i.e. ) types of trinucleotides (e.g., ‘AAA’, ‘AAC’, ‘AAG’, …, ‘UUU’). Therefore, for a given L-bp-long RNA sequence, the trinucleotide position specificity can be formulated as a 64(L-2) matrix:

(5)

where . and represent the frequencies of the -th trinucleotide () at the -th position appearing in the positive () and negative () datasets, respectively. equals ‘AAA’, equals ‘AAC’, …, equals ‘UUU’. Accordingly, an L-bp-long sequence can be represented as the following feature vector :

(6)

where denotes the transpose operator and is defined as follows:

(7)

Therefore, in this study, the samples in H_990 and M_944 are encoded as 21 − 2 = 19 PSTNPss features, while the samples in S_628 are represented as 31 − 2 = 29 PSTNPss features.

Stacking ensemble learning framework of Porpoise

Porpoise is a meta-learning framework designed for improving the prediction of RNA pseudouridine sites. Stacking is an effective ensemble learning strategy that integrates the information of various classifiers to enable the construction of a robust prediction model. Such strategy has been successfully applied in a number of recent bioinformatics and computational biology studies [25–29]. The stacking strategy contains two major steps, and the corresponding classifiers at each step are termed as base-classifier and meta-classifier. A set of base-classifiers are applied and built at the first step, and the outputs of base classifiers are then used as the input to train the meta-classifier at the second step.

In this study, we comprehensively assess the performance of nine commonly used machine learning algorithms for RNA pseudouridine site prediction, including Adaptive Boosting (AdaBoost) [30], extremely randomized trees (ERT), extreme gradient boosting (XGBoost) [31], Gradient Boosting Decision Tree (GBDT) [32], logistic regression (LR), k-nearest neighbors (KNN), support vector machine (SVM), random forest (RF) and Gaussian Naive Bayes (GaussianNB). Firstly, we evaluate the predictive performance of all the 18 types of features. For each machine learning algorithm, we train 18 individual classifiers based on each type of features and select the best performing one according to the Matthew's Correlation Coefficient (MCC) as the candidate base classifier. All the classifiers are built and optimized using the scikit-learn package in Python with ten times 10-fold cross-validation tests. Using this strategy, we obtained nine candidate base classifiers in accordance with nine different machine learning algorithms. In order to determine the optimal combination of base-classifiers for the stacked models, we evaluate eight different combinations of base classifiers. Firstly, we rank these nine base classifiers in terms of their classification performance. We set C as the ranked pool of candidate base classifiers, , where achieved the best MCC. Eight base-classifier combinations are generated by taking base classifiers from , where Ensemble 1 includes , Ensemble 2 includes , Ensemble 3 includes , Ensemble 4 includes , Ensemble 5 includes , Ensemble 6 includes , Ensemble 7 includes and Ensemble 8 includes , respectively. Logistic regression is used as the meta-classifier to train the stacked model. We comprehensively evaluate these eight base-classifier combinations for each of the three species and select the combination that achieved the best performance as the final model, then optimize the hyperparameters of the final model using the Bayesian optimization algorithm [33]. We implement our stacking strategy using the Stacking Cross-Validation algorithm provided in the 'mlxtend' package [34] in Python.

Performance evaluation

The predictive performance of Porpoise and other existing methods is evaluated and compared in terms of several commonly used performance measures [35–37], including sensitivity (Sn), specificity (Sp), precision, accuracy (ACC), F1 score, Matthew's Correlation Coefficient (MCC) and area under the receiver-operating curves (AUC). Sn, Sp, Precision, ACC, F1 and MCC are, respectively, defined as:

(8)

where , , and denote the numbers of true positives, true negatives, false positives and false negatives, respectively.

Performance evaluation of different feature encoding schemes

In this section, we comprehensively investigated and evaluated the performance of variant models trained using 18 different feature encodings in combination with nine commonly used machine learning algorithms (as detailed in Stacking ensemble learning framework of Porpoise section) by performing ten times 10-fold cross-validation tests. We built the corresponding 18 classifiers with the default parameters based on different feature encoding schemes for each machine learning algorithm and evaluate the model performance based on the average performance metrics of ten times 10-fold cross-validation tests. The performance evaluation results are provided in Supplementary Tables S1–S9. As mentioned in Stacking ensemble learning framework of Porpoise section, we selected the best-performing one in terms of MCC as the candidate base classifier for each algorithm. As a result, nine candidate base classifiers for each species are summarized in Supplementary Table S10.

Overall, we do not observe any single type of features that was consistently outperforming other features for any species, and the predictive performance of the models trained using single types of features was generally limited. Despite this, we make several important observations: First, two gradient boosting algorithms, XGBoost and GBDT, achieved the best performance among the nine compared machine learning algorithms. In particular, XGBoost achieved the best performance on H_990 and M_944. In contrast, GBDT achieved the best performance on S_628. Here, it should be noted that we ranked the classifiers according to the MCC values, as we want to construct the best classifier to achieve a better balance between sensitivity and specificity. Second, we find that the PSTNPss features contributed more to the model performance as opposed to other types of features, particularly for these two best-performing algorithms XGBoost and GBDT. For example, the XGBoost model trained using the PSTNPss features outperformed across all three species compared with the XGBoost models trained using 18 different types of features in accuracy, sensitivity, MCC and AUC (Supplementary Table S3). Figure 2 provides a comparison of the AUC values of 18 XGBoost models trained using these different types of features. XGBoost models trained using the PSTNPss features perform the best amongst 18 compared models across all three species. In addition, the GBDT model trained with PSTNPss features also achieved the best performance on H_990 and M_944 datasets in terms of all performance metrics, with the only exception of the specificity on H_990, where the RCKmer-based model performed the best (Supplementary Table S4). Third, the prediction performance of all base classifiers was not satisfying. This indicates that the single model trained with a single type of features could not effectively address this prediction task. As such, we were highly motivated to develop an effective ensemble learning strategy to enhance the model performance.

The optimal base-classifier combinations

Porpoise employed the stacking strategy to build the ensemble learning model. In order to identify the optimal base classifiers for the stacked model from the nine candidate base classifiers, we evaluated eight different combinations of the base classifiers as detailed in Stacking ensemble learning framework of Porpoise section. The performance evaluation results of three species are summarized in Supplementary Tables S11–S13 and Figure 3. As can be seen, it is clear that, for H_990, Ensemble3 achieved the best performance in terms of almost all performance metrics, with the only exception of Specificity. For S_628 and M_944, Ensemble1 secured the best-prediction performance with respect to almost all performance evaluation metrics. This means that the stacking model could achieve a competitive performance based on a few base classifiers on all these three datasets. Specifically, Ensemble3 of H_990 stacked an XGBoost and a GBDT classifier trained with PSTNPss features, a GaNB model trained with PseKNC and an SVM model trained with binary features. The optimal combination for S_628 ensembled a GBDT model and an XGBoost model trained using the PSTNPss features. As a comparison, Ensemble1 of M_944 stacked an XGBoost classifier trained using the PSTNPss features and an AdaBoost classifier trained with the NCP features.

Based on the optimal combinations of base-classifiers, we next conducted feature selection and hyperparameter optimization for each base classifier. Although each base classifier is only trained using one single type of features, some of the features are of high dimension and particularly sparse, which might have a negative impact on model training. Therefore, it might be helpful to further apply feature selection to identify more relevant and contributing features for model training. In this study, we employed a two-step feature selection strategy that combines the mRMR (minimum redundancy maximum relevance) [38] and incremental feature selection (IFS) algorithms [39, 40] to identify the optimal feature subsets. At the first step, we ranked all the initial features according to the information gain. We set F as the whole set of all ranked features, , where n is the number of features. Then, the IFS was applied to the training dataset by performing ten times 10-fold cross-validation tests to evaluate the relative importance of all the features included in F. At each iteration, IFS constructed n features subsets by adding one feature from to the candidate feature subset . The i-th feature subset can be defined as . As a result, the feature subset that achieved the highest MCC value is regarded as the optimal feature subset.

Following feature selection, we employed the Bayesian optimization algorithm [33] to optimize the hyper-parameters of the stacked models. The best-performing parameter values and the corresponding numbers of features after feature selection on the training dataset are provided in Supplementary Table S14. In addition, IFS feature selection results are shown in Supplementary Figures S1–S3. These results demonstrate that, for H_990, all 19 PSTNPss features were selected as optimal features for XGBoost and GBDT classifiers; for S_628, 22 out of 29 PSTNPss features were selected as the optimal features for GBDT and XGBoost classifiers, respectively; while for M_944, 15 out of 19 PSTNPss features were chosen as the optimal features for the XGBoost classifier. Furthermore, after feature selection, the total number of features fed into the stacked model was relatively small. There were only 109 features for H_990, 44 features for S_628 and 50 features for M_944, respectively. To some extent, feature selection helps reduce the stacked models’ complexity, alleviate the resources required for model training and improve the model interpretability.

The base classifiers of Porpoise were all trained using a single type of features. For example, on H_990, XGBoost only applied the PSTNPss features for model training, feature selection and parameter optimization. In comparison, a variety of diverse features could provide more useful information for improving the performance of machine learning algorithms. Therefore, a question will naturally arise—Will the models perform better if an optimal subset of features could be selected from a more comprehensive and larger feature set? To address this question, we subsequently compared and evaluated two contrasting feature selection strategies. We named the feature selection strategy used in Porpoise as Strategy 1, while Strategy 2 is also used for making a comparison. In Strategy 2, we first ranked all the 18 types of features using the mRMR algorithm. As the dimensionality of all 18 types of features is very high, which requires intensive computing resources to execute the IFS strategy. To deal with this, we used a percentage-based IFS strategy to select the optimal features and reduce the computing time. This percentage-based IFS strategy was conducted ten times in the following manner: in the first time, the top 10% of the ranked features were used to train and evaluate the model performance; in the second time, the top 20% of the ranked features were employed to train and evaluate the models, so on and as so forth; each time 10% more features were added, and the model performance was accordingly evaluated; in the 10-th time, all the features are fed into the classifier. Finally, the feature subset that attained the highest MCC value was regarded as the optimal feature subset. Both Strategy 1 and Strategy 2 were repeated ten times to calculate the average performance.

Moreover, we conducted a student’s t-test to assess the statistical significance of the MCC values between Strategy 1 and Strategy 2. Figure 4A plots the distributions of the MCC scores for Strategy 1 and Strategy 2 evaluated on the training datasets, while Figure 4B shows the distributions of the MCC scores on the independent test datasets. The results indicate that Strategy 1 indeed significantly outperformed Strategy 2 in terms of MCC on both benchmark training and independent test datasets, with an only exception on the independent test dataset of S. cerevisiae (S_200), for which the P-value is not significant. Nevertheless, it is apparent that the average MCC values of Strategy 1 are higher than those of Strategy 2; on the other hand, the standard deviation of Strategy 1 is lower than that of Strategy 2. This suggests that Strategy 1 provided better performance than Strategy 2 and achieved a more stable performance than Strategy 2. Therefore, we adopt Strategy 1 to optimize Porpoise models in the following studies.

Performance comparison between Porpoise and state-of-the-art methods

In this section, we employed the same benchmark training and independent test datasets previously used by several state-of-the-art methods to rigorously evaluate and compare the predictive performance of Porpoise against the other methods. Tables 3 and 4, respectively, summarize the performance comparison results of Porpoise and state-of-the-art predictors including iRNA-PseU [11], PseUI [12], iPseU-CNN [13], XG-PseU [14], EnsemPseU [15], RF-PseU [16] and MU-PseUDeep [17] on the same benchmark training and independent test datasets. As can be seen from Table 3, Porpoise achieved the overall best performance compared with the other seven tools across all three species H. sapiens, S. cerevisiae and M. musculus on the same training datasets in terms of accuracy and MCC. On H_990, Porpoise achieved the best performance in accuracy, sensitivity and MCC, while MU-PseUDeep achieved the best specificity. Porpoise’s accuracy, sensitivity and MCC were 5.93, 18.21 and 6.05% higher than the second-best method MU-PseUDeep, respectively. On S_628, Porpoise secures the best performance in terms of all the four measures i.e. accuracy, sensitivity, specificity and MCC, which were 4.89, 3.01, 2.37, and 8.78%, respectively, higher than the second-best performance value. The performance improvement of Porpoise on M_944 was relatively small, on which it achieved the best accuracy, specificity and MCC with a marginal increase of 1.75, 4.67 and 1.85%, respectively.

The superior performance of Porpoise on the benchmark training datasets suggests that the stacking ensemble learning strategy employed can enhance the effectiveness of model training. However, to examine if models trained using this strategy is subjected to overfitting, we further performed an independent test to evaluate and validate the trained models using the independent test datasets. Using the same independent test datasets, we evaluated and compared the predictive performance of Porpoise with that of other existing methods. The results are summarized in Table 4. As can be seen, the Porpoise models also achieved a very competitive performance on both H_200 and S_200. In particular, on H_200, Porpoise achieved the best accuracy, sensitivity and MCC, which were 2.35, 4.3 and 7.13% higher than the second-best performance values, respectively. On S_200, Porpoise achieved the best performance in terms of three performance measures, with an increase of 6.5, 3.0 and 13.27% in accuracy, sensitivity and MCC, respectively. In addition, the performance on the independent test was also very close to that on the cross-validation test on the benchmark training datasets: for example, Porpoise achieved an accuracy of 77.35% and MCC of 55.13% on H_200, which was very close to the performance on the training dataset, on which Porpoise attained an accuracy of 78.53% and MCC of 58.45%, respectively. Besides, Porpoise achieved an accuracy of 83.5% and MCC of 67.27% on S_200, which was slightly higher than that on the training dataset (81.69 and 63.38%, respectively). Taken together, these results demonstrated that Porpoise achieved a competitive performance and compared favorably with the existing state-of-the-art methods, and there was no overfitting for the Porpoise models.

Model interpretation based on the kernel SHAP

In this section, we employed the Shapley Additive explanation (SHAP) algorithm [41] to help interpret the Porpoise models. Although we applied the mRMR feature selection algorithm to rank the features and selected the optimal features using an incremental feature selection. However, the directionality of these optimal features was still unknown. The directionality here means when a feature takes a higher value, we do not know whether this feature indeed corresponds to larger or smaller feature importance for the performance of the classifier. To address this, we used the SHAP algorithm, which can assign each feature an importance value for the model prediction, thereby providing an interpretation of the stacked model of Porpoise [28, 42]. Specifically, in the current work, we used the kernel SHAP instead of the SHAP algorithm to interpret the stacked model of Porpoise, as the stacked model is too complicated to interpret using the latter, which is designed for explaining the single-tree ensemble models. The kernel SHAP applies a specially weighted local linear regression to estimate the SHAP values for any model. Figure 5 shows the top 20 features ranked based on the SHAP value for the three different species, where each row represents the SHAP value distributions of a feature. The color of data points in the figure indicates the corresponding feature value, and redder points mean higher feature values, while bluer points indicate lower feature values. The color (high or low) and the SHAP (positive or negative) value indicate the features’ directionality.

For each panel in Figure 5, the top 20 ranked features sorted according to the sum of SHAP values incorporating all training samples. The points with different colors illustrate the distribution of each feature’s impact on the output of the Porpoise model, where a SHAP value of >0 indicates that the prediction favors the positive class (i.e. pseudouridine site). In contrast, a SHAP value of <0 means that the prediction tends to be of the negative class (i.e. non-pseudouridine site). The results revealed that the PSTNPss features play a predominant role for Porpoise, which was reflected by the presence of all the top 20 important features for H. sapiens and S. cerevisiae models, and the top 15 critical features for the M. musculus model were PSTNPss features. Altogether, these results once again confirm and highlight the importance and contribution of the PSTNPss features to the performance of Porpoise.

Performance evaluation on the Psi_test dataset

In this section, we further evaluated the predictive performance of Porpoise and compared it with other state-of-the-art methods on the Psi_test dataset collected from the m6A-Atlas database. We submitted the FASTA sequences in the Psi_test dataset to the webservers of the existing approaches and accordingly calculated the percentage of corrected predicted samples. For Porpoise, we used the model trained on the training datasets to perform the prediction. The performance comparison results are provided in Table 5. As can be seen, Porpoise achieved the best predictive performance and secured 80.5, 83 and 88.1% in H. sapiens, M. musculus and S. cerevisiae, which was 2.8, 5.2, and 9.9% higher than the second-best performed values. In addition, we applied the webserver of Porpoise to make a time-usage comparison with other existing webservers and provided the comparison results in Supplementary Table S15. The results shown that Porpoise was the second-ranked time-saving webserver to generate the prediction results among the four compared webservers. PseUI was the most time-saving webserver for H. sapiens and M. musculus, while PseUI required considerable time for S. cerevisiae, 30 times longer than that of XG-PseU, which was ranked first for S. cerevisiae. It should be noted that this time-usage analysis was performed only for a rough time comparison purpose between different webservers, which may also be influenced by the status of the internet.

Webserver and software development

In this section, to facilitate community-wide efforts in performing high-throughput analysis and prediction of novel potential RNA pseudouridine sites from RNA sequence data, we implemented an online webserver of Porpoise, which is publicly accessible at http://web.unimelb-bioinfortools.cloud.edu.au/Porpoise/. The web page of Porpoise was developed based on PHP, while the webserver is managed by Apache HTTP Server and configured in a 16-core Linux server machine with 64GB RAM and 500GB hard disk. To use the webserver, users only need to input their sequences of interest or upload an input sequence file in the FASTA format and specify the species type(s). Users’ submitted tasks will be processed at the server-side to complete the prediction. Upon completion, the prediction results will be returned to the webpage or sent via the user’s optionally provided email addresses. A detailed user manual for using the Porpoise web server can be found on the help page of the webserver. The server has been tested using several popular web browsers and works well, including Internet Explorer, Mozilla Firefox, Google Chrome and Safari. In addition, we also make available a stand-alone version of Porpoise, which can be downloaded from the download webpage. Users can download and run the trained models of Porpoise in cases where they need to process large-scale data sets locally.