MRM-BERT: a novel deep neural network predictor of multiple RNA modifications by fusing BERT representation and sequence features

ABSTRACT

RNA modifications play crucial roles in various biological processes and diseases. Accurate prediction of RNA modification sites is essential for understanding their functions. In this study, we propose a hybrid approach that fuses a pre-trained sequence representation with various sequence features to predict multiple types of RNA modifications in one combined prediction framework. We developed MRM-BERT, a deep learning method that combined the pre-trained DNABERT deep sequence representation module and the convolutional neural network (CNN) exploiting four traditional sequence feature encodings to improve the prediction performance. MRM-BERT was evaluated on multiple datasets of 12 commonly occurring RNA modifications, including m⁶A, m⁵C, m¹A and so on. The results demonstrate that our hybrid model outperforms other models in terms of area under receiver operating characteristic curve (AUC) for all 12 types of RNA modifications. MRM-BERT is available as an online tool (http://117.122.208.21:8501) or source code (https://github.com/abhhba999/MRM-BERT), which allows users to predict RNA modification sites and visualize the results. Overall, our study provides an effective and efficient approach to predict multiple RNA modifications, contributing to the understanding of RNA biology and the development of therapeutic strategies.

KEYWORDS:

• RNA modifications
• deep learning
• sequence features
• Bidirectional encoder representations from transformers (BERT)
• convolutional neural network
• web server

Background

Chemical modification is a highly specific and effective way to regulate the function of biomolecules. Almost all types of RNA undergo chemical modification. At least 160 types of chemical modifications have been found in RNA [1], which play important roles in biological processes such as transcription, splicing, translation and gene silencing, affecting the structure, function and interaction of RNA [2]. Therefore, RNA modifications are crucial for almost every aspect of the RNA life cycle and various cellular and developmental processes. RNA modifications also play important roles in diseases, for example, the m⁶A methylation modification writer (the catalysing enzyme of RNA modification) methyltransferase 3 (METTL3) can promote the proliferation of acute myeloid leukaemia (AML) cells [3] and the deficiency of pseudouridine synthase 10 (PUS10), a writer of pseudouridine modification, can protect prostate cancer cells from tumour necrosis factor-related apoptosis-inducing ligand (TRAIL)-induced apoptosis [4].

The target and position of RNA modifications are of the utmost importance in determining the effects of RNA modifications, so it is particularly important to accurately and efficiently identify the RNA modification sites [5]. Recently developed high-throughput sequencing techniques for detecting RNA modification sites have greatly accelerated the functional study of RNA modifications. Taking m⁶A methylation as an example, there are various sequencing approaches for m⁶A identification including antibody-dependent methods like MeRIP-seq [6] and miCLIP [7] and antibody-independent methods like m6A-REF-seq [8], MAZTER-seq [9] and DART-seq [10]. However, these experimental methods are expensive, labour-intensive, and highly sensitive to the experimental conditions. For example, MeRIP-seq requires a large amount of RNA input and cannot map m⁶A sites at single-base resolution, miCLIP requires special antibodies and sophisticated enzyme crosslinking procedures, m6A-REF-seq and MAZTER-seq can cover only a subset of modification sites, while DART-seq is also complex to conduct. Therefore, there is an urgent need for computational methods that can predict RNA modification sites based on multimodal sequencing data generated by various conventional sequencing methods.

Several algorithms have been proposed to predict RNA modification sites based on RNA sequence features, such as SRAMP [11], DL-M6A [12], WHISTLE [13], Gene2vec [14], MTDeepM6A-2S [15], IRNA-m5U [16], m5C-HPromoter [17], and pRNAm-PC [18] These predictors extract various sequence features and utilize different machine learning or deep learning framework to predict RNA modification sites with improving accuracy. However, these methods also have some limitations, such as most of the previous predictors focus only on one specific type of modification, their datasets are mostly from one experimental method, and the sequence representation is not fully comprehensive. To this end, MulitRM [19] assembled datasets of 12 common modifications from multiple experiments and used the word2vec [20] embedding method combined with a long short-term memory (LSTM) attention network to predict RNA modification sites, achieving simultaneous prediction of multiple types of modifications with competitive performance. On the other hand, we note that the sequence representation could be enriched by diverse sequence embeddings and encodings, while the predictor could also be further calibrated with more sizable test samples.

In this study, we proposed a novel method for predicting multiple types of RNA modification sites by combining pre-trained models and various feature encodings. The term pre-trained models here is specifically referred to the deep neural networks that are pre-trained on large-scale unlabelled data. Such models can take full advantage of the large amount of unlabelled data, learn general features and representations of sequences, and thus relieve the limited availability and quality of labelled data, which is also often the case for RNA modifications. Pre-trained models have achieved great success in the field of natural language processing (NLP) [21]. For example, bidirectional encoder representations from transformers (BERT) [22], GPT-3 [23] and other models have reached the state-of-the-art performance in various NLP tasks. Similarly, DNA or RNA sequences can also be considered as a kind of language, with certain grammatical and semantic rules. Therefore, deep features and representations of DNA or RNA sequences can be learned through pre-trained models, and then transferred to learn downstream prediction tasks to improve the accuracy and efficiency of prediction. A typical representative of pre-trained models is DNABERT [24]. As its name implies, DNABERT treats DNA sequences as a kind of text that is similar to human language, and uses two tasks, i.e. masked language model (MLM) and next sentence prediction (NSP) to learn the semantic representation of DNA sequences. DNABERT achieves better results than traditional deep learning models in downstream tasks such as genome fragment classification, genome fragment similarity calculation and genome fragment clustering. Nonetheless, RNA sequences have implicated structural and functional information that are not fully represented by their DNA counterpart. As the result, DNABERT, although highly comprehensive, is not well suited for RNA modification site prediction. To overcome this limitation, we devised a deep learning framework, MRM-BERT, that can fine-tune the pre-trained model DNABERT and enhance the performance of the prediction model by integrating multiple sequence features through a convolutional neural network. We used this method to predict 12 common RNA modification sites from multiple datasets. In the following sections, we will first introduce the establishment of MRM-BERT, and then demonstrate the performance and usage of this predictor.

Methods and materials

Dataset

The RNA modification site data were derived from the MultiRM dataset [19], which consisted of datasets generated by 20 different base-resolution techniques covering 12 different types of RNA modifications (m⁶A, m¹A, m⁵C, m⁵U, m⁶Am, m⁷G, Ψ, AtoI, Am, Cm, Gm and Um). It is noteworthy that the MultiRM dataset had already discarded modification sites identified by ambiguous techniques such as ordinary MeRIP-seq methylation peak calling combined with motif search. As for the negative samples (non-modification sites), they were randomly sampled from the non-modified bases of the same transcript containing the positive sites [19]. In total, the whole dataset contains more than 300,000 sites.

We adopted two strategies for dataset partition. First, to facilitate direct comparison with MultiRM, we used their data partition strategy as ‘MultiRM’s fixed validation & testing’. We divided the dataset of each RNA modification type into three groups, i.e. training set, validation set and test set. We required that both the positive and negative test samples used here were exactly the same as those used by MultiRM. Second, we noted that the proportions of validation and test sets in the original partition of MultiRM were small, and the model could be better calibrated and evaluated with more testing samples. Therefore, we also used the ‘K-fold cross-validation & testing’ data partition strategy with training-validation-test setting, i.e. randomly taking 10% of all samples as the test set, and then dividing the remaining 90% of the data evenly into five parts. Each time, different four parts are used for training and the remaining 1 part for validation, repeated five times.

To address the issue about overlapping/redundant sequence windows, we performed CD-HIT-EST on each modification dataset before model training and testing, removing sequences with a similarity of no less than 0.9. More specifically, for the K-fold cross-validation and testing, the intra-dataset similarities of both training and test sets are <0.9, and the inter-dataset similarity between training and test sets is also <0.9. As for the MultiRM’s fixed validation and testing, to ensure a fair comparison with MultiRM, we did not remove redundant sequences within the test set; however, the inter-dataset similarity between training and test sets and the intra-dataset similarity of training set was still reduced to <0.9 via CD-HIT-EST de-redundancy procedure.

In addition, to ensure a completely independent test set, we collected m⁶A, m⁵C, Psi, and m¹A sites from mouse data downloaded from m6A-Atlas [25], m5C-Atlas [26] and RMBase 3.0 [27,28]. Other modification types had insufficient data to constitute a testing set. CD-HIT-EST-2D was applied to eliminate all mouse testing sequences similar to human training sequences as what has been described above. Data partitions and their respective purposes are shown in Supplementary Figure S1.

Selection of sequence length

To determine the optimal sequence length, we tested sequence lengths of 21, 31, 41, … , 121,131,141. Results indicate a steady increase in the average AUC from length 21 to 101. However, at lengths 121 to 141, there is no significant improvement in AUC compared to length 101 (Supplementary Figure S2). Therefore, we chose 101 as the final length of the sequence windows.

Traditional sequence feature encodings

For the RNA sequence of each sample, we retained 101 bases with the modification/non-modification site at the centre position. To establish an informative combination of traditional sequence feature encodings, different encodings were evaluated and selected by assessing the performance of XGBoost classifiers based on the sequence encodings. More specifically, the following sequence encodings were evaluated:

1. PseKNC (Pseudo K-tuple Nucleotide Composition): This method extracts features by considering the local sequence patterns and the global sequence patterns of consecutive k nucleotides [29].

2. PseDNC (Pseudo Dinucleotide Composition): This method extracts features by considering the local sequence patterns and global sequence patterns of consecutive dinucleotides [29].

3. CKSNAP (k-spaced nucleotide pairs): This method extracts features by considering k-spaced nucleotide pairs [11,30].

4. KMER: This method extracts features by considering all possible combinations of k nucleotides [31].

5. DPCP (Dinucleotide Physicochemical Properties): This method extracts features by considering the physico-chemical properties of dinucleotides [32].

6. NCP (Nucleotide Chemical Property): This method extracts features by considering the chemical properties of nucleotides [33].

7. BINARY: This method extracts features by converting RNA sequences into binary codes [34].

8. ENAC (Enhanced Nucleic Acid Composition): This method extracts features by considering the composition and arrangement of nucleic acids [30].

9. EIIP (Electron–Ion Interaction Pseudopotential): This method extracts features by considering the pseudopotential of electron–ion interaction [35].

10. ANF (Atomic Number Frequencies): This method extracts features by considering the frequencies of different atoms in RNA sequences [33]

Then, we performed prediction of each modification type one by one using the above sequence encodings, using XGBoost [36] as the backend classifier, additionally, we explored all possible combinations, examining the relevance of these features. Finally, NCP, EIIP, ENAC, Binary, CKSNAP and Kmer sequence encodings were chosen as the representative features, as shown in the Supplementary Table S1.

Deep learning framework

As shown in Figure 1, MRM-BERT was composed by integrating two deep neural network modules: one was the CNN module that took the selected traditional sequence encodings as the input; the other was the BERT module that fine-tuned the DNABERT pre-trained model to fit the RNA modification prediction task.

Figure 1. The deep learning framework of MRM-BERT.

MRM-BERT integrated two deep neural network modules to achieve the prediction of multiple RNA modifications. First, the RNA sequences around the modification sites were represented by the selected traditional feature encodings, and then a CNN module based on these traditional feature encodings was trained. Second, the sequences around the modification sites were transformed into 3-mer vectors, and a BERT module was established by fine-tuning the DNABERT representation. Finally, the fully connected (FC) layers from the two deep neural network modules were fused to obtain the final classification prediction.

1. CNN module: The four feature matrices (NCP, EIIP, ENAC, Binary) were input into the CNN module. Considering that bidirectional LSTM would focus on the two ends of the input while the sequences near the centre position were of the utmost importance for RNA modification, CNN would be a better choice to integrate various feature encodings. The inputs were firstly passed through two one-dimensional convolutional layers and max-pooled to reduce the parameter size, resulting in four feature maps, respectively. We concatenated the four feature maps with other two feature matrices (CKSNAP and Kmer) and then passed the combined feature map through two convolutional layers, for each feature matrix, a bundle of k 1D CNN filters, rather than a single 2D CNN filter, was applied along the sequence length (n) dimension. The results of these 1D CNN filters could be finally integrated by the fully connected layer to obtain the CNN-based representation of RNA modification sites.

2. Fine-tuned BERT module: For each sequence, we firstly performed segmentation and masking, dividing the base sequence of length 101 into 99 k-mers of length 3 (i.e. 3-mers) using a sliding window of width 3. For the training set, we masked 15% of the k-mers and replaced them with corresponding tokens. After preprocessing, we converted the obtained k-mers and masks into index embedding in the dictionary, and added position embedding to obtain a fused embedding. The fused embedding was input into the pre-trained DNABERT model that was composed of 12 transformer blocks and then passed through a fully connected layer to fine-tune the DNABERT model and obtain the BERT-based representation of RNA modification sites.

Model training procedure

The CNN and BERT modules were trained or fine-tuned using the same training and validation sets. Subsequently, the parameters of these two networks were fixed and the integrated network was trained. For the CNN module, a grid search algorithm was used, the hyper-parameters of the Adam optimizer were set to a learning rate of 1e-5, weight decay of 0.001, and the learning rate scheduler was used to determine an appropriate learning rate increase and decrease curve. For the BERT module, a pre-trained model with a k-mer of 3 was used, and through grid search, the hyper-parameters were optimized as learning rate of 3e-5, warmup percent of 0.1. These hyper-parameters were determined to gradually increase the learning rate at the beginning of training while avoiding unstable training results when the learning rate is too high. The hidden dropout probability of 0.1 could reduce overfitting and improve generalization performance, while the weight decay of 0.01 could control model complexity to avoid overfitting. Half-precision floating-point arithmetic was also used to increase training speed. All other hyper-parameters were set to default.

Evaluation metrics

We used widely applied evaluation metrics to assess the performance of the model and compare it with other methods. Specifically, the following six evaluation metrics were considered, where TP represents true positive samples, TN represents true negative samples, FP represents false positive samples, and FN represents false negative samples:

• Accuracy refers to the proportion of correctly classified samples in the total number of samples. Its calculation formula is:

  (1) $\mathit{ACC}=\frac{\mathit{TP}+\mathit{TN}}{\mathit{TP}+\mathit{FP}+\mathit{FN}+\mathit{TN}}$(1)

• Sensitivity, also known as Recall, refers to the proportion of positive cases correctly classified. Its calculation formula is:

  (2) $\mathit{SN}=\frac{\mathit{TP}}{\mathit{TP}+\mathit{FN}}$(2)

• Specificity, which refers to the proportion of negative cases correctly classified. Its calculation formula is:

  (3) $\mathit{SP}=\frac{\mathit{TN}}{\mathit{TP}+\mathit{FP}}$(3)

• F1 Score, which is the harmonic mean of Precision and Recall. It is suitable for unbalanced classification problems and can be used as an evaluation indicator when the ratio of positive to negative samples is imbalanced. The calculation formula for the F1 score is:

  (4) $\mathit{F}1=\frac{2\ast \mathit{Precision}\ast \mathit{SN}}{\mathit{Precision}+\mathit{SN}}$(4)
  (5) $\mathit{Precision}=\frac{\mathit{TP}}{\mathit{TP}+\mathit{FP}}$(5)

• MCC is the Matthews Correlation Coefficient, a value between −1 and 1, which indicates the correlation between the classifier and the true label. Its calculation formula is:

  (6) $\mathit{MCC}=\frac{\mathit{TP}\ast \mathit{TN}-\mathit{FP}\ast \mathit{FN}}{\sqrt{\left(\mathit{TP}+\mathit{FP}\right)\ast \left(\mathit{TP}+\mathit{FN}\right)\ast \left(\mathit{TN}+\mathit{FP}\right)\ast \left(\mathit{TN}+\mathit{FN}\right)}}$(6)

• AUC, the Area Under the Receiver Operating Characteristic Curve (ROC) Curve, is the metric to assess the threshold-independent overall performance of the predictors. The ROC curve is plotted with sensitivity as the y-axis and 1-specificity as the x-axis, reflecting the performance of the classifier at different thresholds. The closer the ROC curve is to the upper left corner, the higher the AUC value, and the better the predictor.

Results

Selected traditional feature encodings for the CNN module of MRM-BERT

Our framework fuses CNN deep neural network module and fine-tuned pre-trained BERT module to predict 12 common types of RNA modifications (Figure 1). Given a set of modifiable sites, the network will learn to extract RNA sequence features from the nucleotide sequence and integrate them with the semantic feature encoding of the sequences, to predict the possible modification map from the RNA sequences.

The CNN module combines various traditional sequence encodings. We first evaluated the performance of different encodings and all possible combinations, then selected the best-performing ones as the favourable input of the CNN module. For this purpose, we applied K-fold cross-validation & testing partitioning setting for the dataset (see Methods section for details) and evaluated the performance on the validation set. Notably, because the performance of a CNN model is usually more sensitive to hyperparameter selection and model training strategy, we instead exploited XGBoost, an effective machine learning framework that is less sensitive to training configuration to evaluate the performance of different sequence encodings.

Indeed, not all traditional sequence encodings are informative for RNA modification site prediction. Taking AtoI, Psi, m¹A, m⁵C and m⁶A, five modification types with extensive modification site annotations as examples, we noted that individual encoding methods alone often could not achieve satisfactory performance. On the other hand, among the combinations of six models, the most outstanding one has essentially reached the highest accuracy and AUC among all combinations (Figure 2). Therefore, these six sequence encodings were selected as the input encodings of the CNN module.

Figure 2. Comparison of prediction performance of different sequence encodings combinations. the performance shown in this figure was evaluated on the validation set of the K-fold cross-validation & testing partitioning strategy. A. Boxplots of AUC and accuracy comparing XGBoost models combining different number of encodings for AtoI modification prediction. B. Boxplots of AUC and accuracy comparing XGBoost models combining different number of encodings for psi modification prediction. C. Boxplots of AUC and accuracy comparing XGBoost models combining different number of encodings for m [1]A modification prediction. D. Boxplots of AUC and accuracy comparing XGBoost models combining different number of encodings for m [6]A modification prediction. E. Boxplots of AUC and accuracy comparing XGBoost models combining different number of encodings for m [5]C modification prediction. F. The best-performing feature combinations and their worst AUC rankings (at worst within 15, or at worst within 30) across all modification types.

Integration of fine-tuned BERT module further improves the overall accuracy

To demonstrate the effectiveness of the fusion model that integrates the CNN module with the fine-tuned pre-trained BERT module, we used the CNN and BERT modules separately to build the predictor. The performance evaluation is summarized in Table 1. In general, the performance of the BERT module is better than the CNN model, indicating that the pre-trained sequence representation could more sophisticatedly depict the sequence pattern around modification sites. This advantage is more pronounced for modification types with less known modification site annotation (e.g. m⁵U and Nm) and those with less prominent sequence motif (e.g. AtoI). Moreover, in most cases except m⁶Am, the model integrating CNN and BERT modules can achieve improved performances in comparison with at least one of the two modules, indicating the effectiveness of the integration.

Table 1. The AUC of the MRM-BERT model versus models using individual deep neural network module.

Model	m^6A	m^5C	Am	Cm	Gm	Um	m^5U	m^1A	m^6Am	m^7G	AtoI	Psi
CNN	0.8705	0.9268	0.9166	0.9421	0.9181	0.9208	0.9496	0.9364	0.9154	0.7742	0.7422	0.7652
BERT	0.8993	0.8775	0.9388	0.9197	0.9418	0.9045	0.9518	0.9505	0.9264	0.6786	0.8409	0.8407
CNN+BERT	0.9004	0.9018	0.9404	0.9321	0.9450	0.9187	0.9549	0.9510	0.9267	0.7430	0.8426	0.8452

The performance shown in this table was evaluated on the validation set of the K-fold cross-validation & testing partitioning strategy.

Performance of MRM-BERT on independent test datasets

On the validation set of the K-fold cross-validation & testing partitioning strategy, the CNN-BERT has shown competitive performance in various aspects (Supplementary Figure S3 and Table S2). To avoid the possible overestimation of performance information leakage, we further evaluated the performance of MRM-BERT on independent test samples from this data partitioning strategy.

By this dataset partitioning, the proportion of test sets is fixed at 10%, and the remaining 90% samples are evenly divided into five parts for cross-validation. There is no overlap between the test samples and the training or validation samples. Using the procedure, a sizable and independent test set was prepared for the performance evaluation. The evaluation results are summarized in Table 2, and the performance of MRM-BERT on these test samples is comparable to that on the validation samples, indicating the robustness of MRM-BERT’s performance.

Table 2. Performance metrics of MRM-BERT on the test set from the K-fold cross-validation & testing partitioning strategy.

Mertrics	m^6A	m^5C	Am	Cm	Gm	Um	m^5U	m^1A	m^6Am	m^7G	AtoI	Psi
Acc	0.8505 ±0.0100	0.9440 ±0.0331	0.9566 ±0.0294	0.9567 ±0.0129	0.9784 ±0.0102	0.9429 ±0.0260	0.9606 ±0.0175	0.9336 ±0.0099	0.8923 ±0.0266	0.8859 ±0.0579	0.9230 ±0.0297	0.7522 ±0.0331
Sn	0.8737 ±0.0104	0.9245 ±0.0416	0.9324 ±0.0497	0.9517 ±0.0145	0.9669 ±0.0155	0.9340 ±0.0272	0.9460 ±0.0236	0.9288 ±0.0072	0.8339 ±0.0526	0.8589 ±0.0607	0.8693 ±0.0608	0.6498 ±0.1207
Sp	0.8274 ±0.0197	0.9632 ±0.0253	0.9789 ±0.0121	0.9602 ±0.0166	0.9873 ±0.0101	0.9511 ±0.0316	0.9709 ±0.0139	0.9377 ±0.0124	0.9550 ±0.0139	0.9107 ±0.0684	0.9631 ±0.0142	0.8369 ±0.0940
MCC	0.7020 ±0.0195	0.8887 ±0.0657	0.9141 ±0.0575	0.9107 ±0.0264	0.9562 ±0.0207	0.8859 ±0.0520	0.9187 ±0.0362	0.8664 ±0.0199	0.7927 ±0.0475	0.7729 ±0.1167	0.8437 ±0.0597	0.5071 ±0.0604
AUC	0.9246 ±0.0066	0.9827 ±0.0122	0.9815 ±0.0152	0.9794 ±0.0072	0.9933 ±0.0076	0.9789 ±0.0119	0.9873 ±0.0054	0.9731 ±0.0050	0.9544 ±0.0054	0.9373 ±0.0398	0.9715 ±0.0163	0.8492 ±0.0238
F1	0.8535 ±0.0089	0.9422 ±0.0346	0.9532 ±0.0329	0.9473 ±0.0153	0.9750 ±0.0118	0.9404 ±0.0268	0.9520 ±0.0214	0.9282 ±0.0105	0.8884 ±0.0304	0.8786 ±0.0614	0.9053 ±0.0394	0.6992 ±0.0655

The performance shown in this table was evaluated on the test set of the K-fold cross-validation & testing partitioning strategy.

We further compared MRM-BERT with MultiRM, the state-of-the-art predictor of multiple RNA modifications. Notably, MultiRM adopted a different dataset partitioning strategy, where a fixed number of samples (50 positive samples for each modification type) were assigned to the test set. To ensure fair comparison with MultiRM, we changed the dataset partitioning to ensure that the same positive and negative test samples were used as in MultiRM. We re-trained MRM-BERT using a dataset that excluded these test samples and compared the performance of MRM-BERT on the original test samples of MultiRM. It is noteworthy that because of the relatively small size of these test sets, the performance would be different than those observed in more sizable datasets. Nonetheless, as shown in Table 3 and Supplementary Table S4, MRM-BERT show improved overall performances than MultiRM. The most prominent improvements were observed in modification types with less known modification site data or less clear sequence motif, including Am, mA [1] and AtoI, with increases in AUC greater than 0.1. Considering that such modification types can be better handled by the fine-tuned BERT module, these improvements would be largely attributed to the integration of BERT pre-trained deep neural network that enables a more comprehensive representation of sequence features of RNA modifications.

Table 3. Comparison of the AUC results on the original test set of MultiRM’s fixed validation and testing partitioning strategy.

AUC	m^6A	m^5C	Am	Cm	Gm	Um	m^5U	m^1A	m^6Am	m^7G	Atol	Psi
CNN+LSTM+Attention	0.7272	0.7732	0.6088	0.7872	0.9304	0.8004	0.8592	0.7548	0.8324	0.4968	0.5116	0.6452
HMM+LSTM+Attention	0.7472	0.8892	0.6984	0.8352	0.8942	0.8916	0.9404	0.8404	0.9276	0.5284	0.5908	0.7216
MultiRM	0.8558	0.9056	0.7892	0.8604	0.9256	0.8780	0.9480	0.7792	0.8912	0.6772	0.6698	0.8528
MRM-BERT	0.8912	0.9072	0.9248	0.9344	0.9164	0.9524	0.9776	0.9460	0.9520	0.6852	0.8744	0.8861

The performance shown in this table was evaluated on the original test set of MultiRM’s fixed validation & testing partitioning strategy. For MultiRM and its variants, the performance metrics are from its original literature.

Additionally, we evaluated the cross-modification prediction performance for the 12 modification types, as shown in Supplementary Figure S4. Generally, the models performed best on the same modification type, while they also showed good performance among similar modification types (for example, the model for Am also performed well among Cm, Gm and Um).

In order to compare MRM-BERT and MultiRM on a completely independent test set, we introduced mouse data of four common RNA modification types (m⁶A, m⁵C, Psi, and m¹A sites). Due to the cross-species nature of the independent test set from mouse, the performances of MRM-BERT and MultiRM on the mouse test set are not as impressive as those on the human test sets. Nonetheless, MRM-BERT could still outperform MultiRM on the mouse data, indicating its better robustness (Supplementary Table S3).

MRM-BERT online web server

To facilitate the community, we have made MRM-BERT available as an online web server. MRM-BERT web server can be accessed at http://117.122.208.21:8501. It is free to all users and does not require registration. MRM-BERT is compatible with common web browsers. It integrates the trained model parameters proposed in this paper, provides prediction and display of prediction results for RNA sequences, and provides advanced visualization tools to explore and analyse predicted modification sites in a user-friendly interactive manner. For users who would have trouble in accessing this web server, we also provided the source code in GitHub (https://github.com/abhhba999/MRM-BERT), together with the runnable models and datasets in FigShare (doi: 10.6084/m9.figshare.24873195, reference number 24,873,195).

A sample prediction result is shown in Figure 3. Users simply need to paste the sequence into the corresponding window on the server and select the required modification-type prediction option. MRM-BERT will automatically perform sequence embedding, feature calculation, deep neural network inference, predict RNA modification sites, and generate functionally associated interactive graphic representations. This visualization of results is helpful for screening target sites for subsequent experimental assays.

Figure 3. A sample prediction result from the online MRM-BERT webserver.

Discussion

In this work, we have developed a multi-label deep neural network model that integrates traditional sequence feature encodings and implicated sequence representation from the fine-tuned pre-trained model. This model is capable of simultaneously predicting multiple several common RNA modifications. To fully exploit the multi-level features of RNA sequences, we employed the pre-trained model DNABERT, which was trained on a large amount of DNA data, and fine-tuned it for RNA tasks, thereby obtaining high-dimensional features applicable to RNA tasks at various levels. On this basis, we integrated a deep neural network utilizing multiple RNA sequence features and achieved the improved performances, validating the effectiveness of our model design. These results suggest the enormous potential of pre-trained models in RNA modification prediction, i.e. pre-trained models can be adapted to different types of RNA modifications, enhancing the model’s generalization capability and adaptability, and can be fine-tuned or used for transfer learning on specific tasks, improving the flexibility and scalability of the model.

The BERT module plays a greater role more often, while in a few modifications, such as m⁵C, m⁷G, Cm, CNN contributes more prediction effects. Our CNN+BERT achieved the best performance for approximately two-third of the cases. Indeed, for particular modification types, the CNN would be a little better than CNN+BERT. However, since the objective of this study is to simultaneously predict multiple modification types under the same computational framework, we chose the CNN+BERT framework, which showed superior performances for most cases. On the other hand, for ad hoc model of a specific modification type, sequence encodings, deep learning framework and inclusion of pre-trained models could be further changed and optimized, but this is not the objective of this study.

To facilitate the potential users, we have constructed a web server with a user-friendly interface, and made the dataset, code and model publicly available on the website. Our model is already capable of predicting multiple different types of RNA modifications, but in reality, its capability depends on the dataset and can be applied to other types of RNA modifications with appropriate training. Nonetheless, since the BERT representation has been fine-tuned, the deep learning framework of MRM-BERT can be easily transferred to the prediction of other types of RNA modifications, once a substantial amount of new modification data become available.

On the other hand, MRM-BERT has certain limitations. Although our model can accurately predict possible modification sites of different modification types, this does not ensure reasonable prediction of the modification states of these sites on a transcript. In other words, the in vivo modification state of each predicted site is determined by the epitranscriptome regulatory dynamics, therefore not all of the predicted sites are modified in the same transcript at the same time. Interestingly, recently emerging nanopore-based third-generation sequencing and single-cell sequencing would provide useful indications for the modification state dynamics of each modification site. Therefore, further integration with these new generation of sequencing data will be a promising direction for the context-dependent prediction of modification sites and states for dynamic epitranscriptome analysis.

In addition, directly integrating more updated feature encoding methods into the pre-trained model may further enhance the performance of the pre-trained model in extracting features, which is beneficial to subsequent predictions. However, equipment limitations make it difficult for us to spend months using GPUs for customized large-scale pre-training. As the performance of laboratory equipment improves, we hope to have the opportunity to train large pre-trained models that integrate more feature encodings ourselves in the future.

In conclusion, our results demonstrate that the use of pre-trained models to predict RNA modification sites from RNA sequences is an effective and cutting-edge computational method, bringing new opportunities and challenges to RNA modification research. We anticipate that with the advancement of deep learning methods and sequencing methods, we will be able to achieve better computational mapping of the epitranscriptome landscape in the future.

Authors’ contributions

YZ and LW conceptualized the study; LW and YZ designed the methodology; LW performed the analysis; LW and YZ established the online software; LW drafted the manuscript; YZ supervised the study and revised the manuscript. All authors read and approved the final manuscript.

Acknowledgments

We thank Leibo Liu at Peking University for his technical assistance in configuring the server environment.

Disclosure statement

There is no relevant conflict of financial or non-financial interest.

Data availability statement

The data that support the findings of this study are openly available in FigShare at doi: 10.6084/m9.figshare.24873195, reference number 24,873,195. The source code of this study is openly available in Github at https://github.com/abhhba999/MRM-BERT

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15476286.2024.2315384

Additional information

Funding

This study was supported by the National Key Research and Development Program of China (2021YFF1201201 to YZ) and National Natural Science Foundation of China (32070658 to YZ).