M6ATMR: identifying N6-methyladenosine sites through RNA sequence similarity matrix reconstruction guided by Transformer

Problem description and datasets

One of the focuses of RNA modification research is site recognition. For this task, the computational methods are usually to convert the problem into a binary classification problem, which takes the RNA sequence information as the initial input to the classification model and gets the probability value of modification sites. The method in our article follows this paradigm. For a given sequence $\mathrm{X}=\left\{x_1,x_2\dots x_n\right\}$, we summarize the model objective as $\mathrm{Y}=f_{\sigma }\left(X\right)\in \left\{0,1\right\}$, where $f_{\sigma }(\cdot )$ is the optimal mapping relationship, and $\mathrm{Y}$ is the prediction label. If the prediction result is a positive sample, $\mathrm{Y}=1$, otherwise $\mathrm{Y}=0$. It is worth noting that, as in most studies, the inclusion of m6A modification sites is used as the classification criteria for positive and negative samples. That is, for a given sequence $\mathrm{X}$, if its central position is the modified site, sequence $\mathrm{X}$ is regarded as a positive sample. Negative samples do not contain modification sites. To this end, we select the RNA m6A modification site dataset of Arabidopsis thaliana (A101 dataset) (Wan et al., 2015) for study. The dataset contains 2,100 samples, in which the ratio of positive and negative samples is 1:1.

Model description

M6ATMR implements the identification task of m6A modification sites based on several procedures. First, the RNA sequences were processed into readable coding representation based on k-mer algorithm, and the similarity between RNA sequences was measured by Manhattan distance to further construct the similarity matrix. After that, M6ATMR learned the structural information of RNA sequences through the Transformer encoder, and converts this structural information into the structural graph by Manhattan distance. It then inputted the similarity matrix and structure graph into the self-correlation graph neural network to iteratively update the similarity matrix and obtain the low-dimensional representations of each RNA sequence. These low-dimensional representations were passed through a local-global fusion block to generate the final fusion representation, which was fed into the RF for identification. For the convenience of description, section 2 describes M6ATMR in detail from four parts: similarity matrix calculation, structure graph learning, similarity matrix optimization, and local-global representation fusion. The framework of M6ATMR is shown in Fig. 1, and the details of each part are as follows.

Similarity matrix calculation

For similarity matrix calculation, the key step is to transform RNA sequences into the readable representations using k-mer frequency statistics. Statistical k-mer frequency information can reveal the distribution law of seed sequences and is an important tool to study sequence similarity. In sequence $\mathrm{X}$, a substring of length $K$ refers to $K$ monometric units starting from any position in $\mathrm{X}$, which is called k-mer. In this article, considering that each of the three adjacent nucleotides in an mRNA molecule is organized into a group that represents the pattern of a particular amino acid in protein synthesis, therefore, to make the model biological interpretation, we set $K$ to 3. The algorithm requires that the starting position of the sequence set $\mathrm{\Omega }=\left\{X_1,X_2,\dots ,X_m\right\}$ is aligned. For k-mer sequences with fixed $K$ value, at the offset position $l\left(0\le l\le n-k\right)$, we count the occurrence frequency of different substrings in k-mer sequence with length $K$ starting from the offset position $l$:

(1) $$R_X=\mathit{\text{Nom}}\left({\oplus }_{i=l}^{n-k}\mathit{\text{fre}}\left(M_i\right)\right)$$where $R_X\in {\mathbb{R}}^{n\times 1}$ denotes the k-mer representation of sequence $X$, and $fre(\cdot )$ represents the frequency of the k-mer substring $M_i$, and $\oplus$ represents the operation of concatenating all substring frequencies. $Nom(\cdot )$ represents a normalized operation.

After representing all RNA sequences using statistical k-mer frequency, we attempt to calculate the similarity between each sequence and construct an initial similarity matrix of RNA sequences. We construct matrices between k-mer representations for the following reasons: Firstly, earlier studies on the prediction of biomedical entity associations (Shao et al., 2020) have shown that similarity matrices can reflect the potential association between different entities, which provides strong support for the establishment of the association between different sequences of nucleotides. Secondly, constructing a graph data structure based on the similarity matrix allows us to capture the relationships between multiple sequences, which creates conditions for further sequence extraction. Therefore, it is a reasonable choice to transform the sequence representing the problem into the optimization task of the similarity matrix. In this article, we choose to employ the Manhattan distance to further measure the degree of similarity between nucleotides on different sequences:

(2) $$S=\left(\begin{array}{ccc}{d}_1^m& \cdots & d_1^m\\ ⋮& \ddots & ⋮\\ d_m^1& \cdots & d_m^m\end{array}\right)$$where, in the similar matrix $\text{S}\in {\mathbb{R}}^{m\times m}$, each element represents the Manhattan distance between k-mer sequences of corresponding positions. Taking $d_1^m$ as an example, its value is the Manhattan distance between sequence representation $R_1$ and $R_m$:

(3) $$d_m^1={\sum }_{j=1}^n\left|R_1^j-R_m^j\right|$$where, the $R_1^j$, $R_m^j$ represent the $jth$ k-mer value in sequence state $R_1$ and $R_m$ respectively, and $|\cdot |$ denotes the operation of taking the absolute value.

Structure graph learning

When obtaining the initial similarity matrices of the RNA sequences, we use a Transformer encoder to capture the structural information and learn the structural graph of RNA sequences. To this end, we apply the encoder part of the Transformer to further process the k-mer sequence. For a given Transformer encoder, there are two basic components: the position encoding block, and the self-attention mechanism. In addition, in order to further explore the structural relationship between RNA sequences, we also calculate the Manhattan distance between vector representations of the output of the Transformer encoding block, and strictly constrain the value of the structure matrix within the set of $\left\{0,1\right\}$. The details are as follows.

Position encoding of Transformer is a functional encoder, that is, position vectors are calculated for each element in the sequence:

(4) $$\begin{array}{c}{\overrightarrow{p_t}}^{\left(i\right)}=f{\left(t\right)}^{\left(i\right)}:=\{\begin{array}{c}\sin \left(w_i\cdot t\right)i=2k\\ \cos \left(w_i\cdot t\right)i=2k+1\end{array}\end{array}$$where $w_i$ denotes the frequency, which is calculated as follows:

(5) $$w_i={\displaystyle \frac{1}{{10000}^{2i/d}}}$$where $d$ is the output dimension of the neural network. It is worth noting that the length of this position vector is equal to the length of k-mer representations. Thus, the RNA sequences are represented by k-mer and position encoding:

(6) $$R_t=R+\overrightarrow{p_t}$$where, $\overrightarrow{p_t}$ denotes the position vectors of all elements in k-mer representations. After calculating the position vectors, the encoder further optimizes the representation $R_t$ through the self-attention mechanism:

(7) $$\{\begin{array}{c}{Q}^i=w^q{X}^i\\ K^i=w^k{X}^i\\ V^i=w^v{X}^i\\ g_X^i=\gamma \left({\displaystyle \frac{Q^i{(K^i)}^T}{\sqrt{d}}}\right)V^i\end{array}$$where, $Q^i$, $K^i$ and $V^i$ are the query matrix, key matrix and value matrix respectively. $\mathrm{\gamma }$ denotes the $softmax$ activation function. Through the self-attention mechanism, the Encoder can reveal the potential associations within the sequence and further excavate the structural associations of nucleotides.

In this regard, we believe that it is a valid way to describe the structural association of RNA sequences through internal relationships captured by the self-attention mechanism. To this end, we also measure this structural associations by the same method in section 2.3 and construct the structure graph $G_{st}$ for sequences. Since the model samples the same data and distance formulas during the calculation process, the structure graph has a potential correlation with the similarity matrix, which indicates that it is reasonable to optimize the similarity matrix further through $G_{st}$. In addition, we ensure that the values of elements in $G_{st}$ are strictly constrained in the set of $\left\{0,1\right\}$ by threshold filtering, as shown below:

(8) $$G_{st}^{\left(i,j\right)}=\{\begin{array}{c}1ifG\left(i,j\right)>0.5\\ 0otherwise\end{array}$$

For each element $G\left(i,j\right)$ in $G_{st}^{\left(i,j\right)}$, updating the value to 0 if $G\left(i,j\right)$ is less than 0.5, otherwise updating the value to 1.

Similarity matrix optimization

To optimize the sequence similarity matrix and update the sequence structure graph, we design a self-correlation graph neural network that does not depend on prior node representations. In traditional graph neural networks, the embedded learning process relies on the existing representations of the nodes or edges in the graph. These prior representations serve as the starting point for the learning algorithm to update and refine the embeddings based on the graph structure. For instance, in some studies related to drug-drug association prediction, the SMILES of drugs are often applied as the prior information to serve as the initial input to the graph neural network. However, for similar matrix and structure graphs of sequences, the initial information of nodes is difficult to be obtained. In addition, both the similarity matrix and the structure graph describe the self-correlation property of the sequence, which makes the introduction of additional prior information may lead to misleading information. Therefore, inspired by the self-attention mechanism, we generate the learnable initial node representations based on the input matrix information. The initial representations are constantly updated and optimized during the process of graph convolution and participate in the optimization of the final embeddings. Taking the processing of similarity matrix as an example, we describe in detail the learning process of representation of similarity graph.

For a given similarity matrix $S\in {\mathbb{R}}^{m\times m}$, we define the node initial representation matrix as $Er\in {\mathbb{R}}^{m\times m}$. Each value in the matrix is determined and iteratively optimized by the network. $S$ and $Er$ are input into the three-layer self-correlation graph neural networks to get the embeddings related to the similarity matrix $S$:

(9) $$\begin{array}{c}\{\begin{array}{c}{R}_s^{\left(i+1\right)}=R_s^{\left(i\right)}+{\alpha }^{\left(i+1\right)}Gcov{\left(Er,S\right)}^{\left(i+1\right)}\\ {\alpha }^{\left(i+1\right)}={\displaystyle \frac{1}{I+1}\left(i+1\right)}\end{array}\end{array}$$where, $\alpha$ is the scaling superparameter to prevent the elements in the similarity matrix representation from becoming infinitesimal during graph convolution, and $I=3$ denotes the number of convolution layers. $Gcov(\cdot )$ represents the convolution process. The hidden layer representation of layer $\left(i+1\right)$ and the representation of layer $\left(i\right)$ satisfies the following equation:

(10) $$H^{\left(l+1\right)}=\sigma (S^{-{\displaystyle \frac{1}{2}}}\left(diag\left(D{e}_s\right)-{\displaystyle \frac{1}{2}\left(S+S^T\right))S^{-{\displaystyle \frac{1}{2}}}{H}^{\left(l\right)}{W}^{\left(l\right)}}\right)$$where $D{e}_s$ denotes the degree matrix of the similar matrix $S$, and $diag(\cdot )$ denotes the diagonalization operation. $W$ is the learnable weight. The representation of the hidden layer is initialized to $Er$, that is, $H^{\left(0\right)}=Er$. Similarly, sequence structure graphs are fed into the self-correlation neural networks in the same way. The network further employs learnable initial representations to mine the self-correlation of sequences, which also ensures consistency between representations learned from similarity matrices and that learned from structure graphs.

In addition, another task of the networks is to get better RNA representations by optimizing the similarity matrix. Hence, we apply the reconstructed similarity matrix and sequence structure graph to calculate the loss of the networks:

(11) $$\mathcal{L}\left(\hat{S},S,\widehat{G_{st}},G_{st}\right)=BCELoss(\hat{S},\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{m}\left(\mathrm{S}\ast S^T\right)+BCELoss(\widehat{G_{st}},\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{m}\left(G_{st}\ast {G_{st}}^T\right)$$where $BCELoss(\cdot )$ denotes the binary cross-entropy loss. In the optimization process, we employ the difference between the reconstructed element values and the original element values to measure the performance of the representations.

Local-global representation fusion

In order to obtain comprehensive knowledge of RNA sequences, we propose a local-global strategy for fusing learned embedded representations. We process two kinds of representations of the same sequence respectively from the local and global perspectives of embedding representations, and further determine the weight relationship between embedding representations from similarity matrices and embedding representations from sequence structure graphs. Specifically, we design a local information extraction block and a global information extraction block respectively to calculate the weights of the two types of embedded representations.

For a given similarity matrix embedding representation $e^{sm}$ and structure graph embedding representation $e^{st}$, we first treat them as residue sequences and weight them to obtain an overall representation $e^a$:

(12) $$e^a=e^{sm}+e^{st}$$

We then enter $e^a$ into the local and global information extraction blocks. For the local and global information extraction block, the extraction process is described as follows:

(13) $$\{\begin{array}{c}{e}_{lo}^a=f\left(\vartheta \left(f^{\prime }\left(e^a\right)\right)\right)\\ e_{gl}^a=f\left(\vartheta \left(f^{\prime }\left(\delta \left(e^a\right)\right)\right)\right)\end{array}\text{}$$where $e_{lo}^a$ and $e_{gl}^a$ denote the output of the local extraction block and that of the global extraction block respectively. $f$ and $f^{\mathrm{\prime }}$ represent a one-dimensional convolution layer containing normalized functions respectively, and $\vartheta$ represents the $ReLu$ function. For the global information extraction block, we add the global average pooling layer $\delta$ on the basis of the local information extraction block. After that, we further calculate the weight difference $w_f$ between the two representations:

(14) $$w_f=sigm\left(e_{lo}^a+e_{gl}^a\right)$$

We utilize this weight difference to further integrate the two types of embedded representations to obtain the similar-structural representation $e^{ss}$:

(15) $$e^{ss}=w_f\ast e^{sm}+\left(1-w_f\right)\ast e^{st}$$

In addition, we consider the indispensable role of the k-mer representations of the sequences for the recognition of m6A modification sites, and further integrate these representations with the similar-structural representations to obtain the final embedding representations $e^{fi}$:

(16) $$e^{fi}=w_{f^{\mathrm{\prime }}}\ast e^{ss}+\left(1-w_{f^{\mathrm{\prime }}}\right)\ast e^{km}$$where, $e^{km}$ denotes the K-mer representations of the sequences, and $w_{f^{\mathrm{\prime }}}$ is the weight difference between $e^{km}$ and $e^{ss}$.

Experiments setting

For a binary classification problem, its prediction states can be divided into the four categories: true positive (TP), false positive (FP), true negative (TN), false negative (FN). Thus, we select some predictive indicators to evaluate the prediction effect, and the calculation processes of these indicators are as below. Moreover, we obtain the area under the precision-recall curve (AUPR) and area under the receiver-operating characteristic curve (AUC) for evaluating our model.

(17) $$Accuracy\left(Acc\right)={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(18) $$F1={\displaystyle \frac{prec\times Sn}{prec+Sn}}$$

(19) $$Precision\left(Prec\right)={\displaystyle \frac{TP}{TP+FN}}$$

(20) $$Sensitivity\left(Sn\right)={\displaystyle \frac{TP}{TP+FN}}$$

(21) $$Specificity\left(Sp\right)={\displaystyle \frac{TN}{TN+FP}}$$

(22) $$\begin{array}{c}MatthewsCorrelationCoefficient\left(MCC\right)\\ ={\displaystyle \frac{TP\times TN-FP\times FN}{\sqrt{\left(TP+FP\right)\left(TP+FN\right)\left(TN+FP\right)\left(TN+FN\right)}}}\end{array}$$

In our study, we set the learning rate to 0.0001 and set the number of layers as three in the self-correlation graph neural network. Considering the time and complexity of training, we reduced the number of heads from eight to six in the Transformer encoder. We set the embedding size of the encoder to 32, the hidden dimension of the feed-forward layer to 128, and the number of encoder blocks to six.

Performance on A101 datasets

We conduct a rigorous 10-fold cross validation to evaluate the performance of our proposed model on the A101 dataset. The dataset is systematically partitioned into ten subsets of equal size, ensuring non-overlapping test sets in each fold. For each fold, we utilize nine subsets for training and one for testing. During evaluation, we consider six essential indicators: AUPR, AUC, Acc, F1 score, Prec, and Sen. The results of the 10-fold cross validation are presented through PR curves and ROC curves in Fig. 2 and summarized in Table 1. Across the 10 folds, the AUC curve demonstrates remarkable stability, with the maximum AUC reaching 93.87% and the minimum at 89.79%. Similarly, the PR curve displays consistent performance, with the highest AUPR at 93.17% and the lowest at 87.66%. Table 1 reveals outstanding performance in various evaluation metrics. The Acc achieves an impressive 84.43%, signifying a high correct identification rate for both TN and TP samples. Additionally, the MCC attains a value of 83.72%, reflecting the overall strength of our model. Furthermore, the values of other metrics, including F1 score, Prec, and Sen, surpass 80%, indicating the reliability of our model. These experimental findings substantiate the robustness and efficacy of our proposed model in identifying m6A modification sites.

Performance comparison of Classifiers

In some studies, the downstream classifiers have shown to significantly influence the classification performance of RNA sequence representations generated by models, potentially leading to model instability. To demonstrate the stability and efficacy of our proposed model, we conduct a comparison experiments using five additional classifiers: logistic regression, eXtreme Gradient Boosting (XGBoost), Light Gradient Boosting Machine (lightgbm), CatBoost, and support vector machines (SVM). Logistic regression employs maximum likelihood estimation to predict model parameters, yielding binary results by minimizing cross-entropy loss during data training. XGBoost, an integrated classification algorithm, employs multiple simple base learners to iteratively train input data, continually reducing the discrepancy between model and input values. In contrast, lightgbm stands out with its advantage of low memory usage and faster training speed. CatBoost is designed to extract the most information from given data and is particularly effective for small machine-learning datasets. SVM, as a binary classification model, aims to find an optimal hyperplane for sample segmentation. For evaluation, we utilize 10-fold cross validation on the A101 dataset, as mentioned in “Performance on A101 datasets”. The same set of indicators, AUPR, AUC, Acc, F1 score, Prec, and Sen are employed for assessment. The experimental results, presented in Fig. 3, Fig. 4, and Table 2, indicate that when RF is used as the downstream classifier, the model achieves the highest performance, with AUC at 91.27% and AUPR at 90.40%. While XGBoost exhibits relatively inferior performance compared to logistic regression and RF, the overall classification effect of all three classifiers remains relatively favorable. The difference in AUC values between XGBoost and RF is 7.84%. These findings support the notion that the RNA sequence representations learned by our model are not easily influenced by the choice of downstream classifiers, indicating the model’s stability and effectiveness in feature extraction and m6A modification site identification.

Features selection

In this study, we integrate three types of features for comprehensive analysis, including feature representations from similarity matrices, feature representations from sequence structure diagrams, and k-mer sequence representations. K-mer representations are pre-coded representations based on the frequency count of sequence k-mer substrings, while the other two features are dynamically learned through neural networks. To achieve a holistic understanding of the sequences, we perform local-global integration of these features. To illustrate the importance of combining these three features, we construct two additional features, drop-raw and drop-trans, to validate our model’s performance. Drop-raw represents features without similarity matrix information, and drop-trans represents features without sequence structure graph information. The comparison results of the three types of features are presented in Fig. 5. Our model demonstrates the best recognition performance when all features are utilized. Drop-raw performs slightly better than drop-trans in terms of AUC and AUPR, suggesting that similarity matrices exert a stronger influence on the model compared to sequence structure graphs. Overall, all three types of features are essential and significantly enhance the effectiveness of site recognition. The local-global fusion block is a crucial component of our model, which integrates multiple features from both local and global perspectives by combining learned similarity matrix features with sequence-structure graph features. To demonstrate the necessity of this fusion block, we design another strategy using only weighted fusion, and the results are also presented in Fig. 5. The outcomes reveal that the recognition performance of the model is inferior when only weighted fusion is employed compared to the local-global fusion block. This underscores the significance of the local-global fusion block in achieving superior performance and highlighting its role in effectively capturing comprehensive sequence information.

Performance comparison of predictors

In this study, we conduct a rigorous 10-fold cross validation to compare our model with five other existing methods on the A101 dataset. The compared methods include M6AMRFS, BERMP, RFAthM6A, BERT-m7G (Zhang et al., 2021), and the model designed by Le & Ho (2022). M6AMRFS encodes RNA sequences using two feature descriptors, dinucleotide binary coding, and local site-specific dinucleotide frequency. It enhances feature representation through the F-score algorithm combined with sequence forward search (SFS) and employs XGBoost as the downstream classifier. BERMP utilizes GRU to represent RNA sequences and adopts an end-to-end training process for site recognition. RFAthM6A attempts to classify various types of features derived from RNA sequences using machine learning methods. Our model, M6ATMR, adopts the transformer encoder to extract sequence representations and uses a stacking ensemble classifier for predicting m6A sites. We also consider two other transformer-based models, BERT-m7G, and the model designed by Le & Ho (2022). In their approaches, BERT-m7G uses bidirectional encoder representations from transformers (BERT) to extract sequence representations, while Le & Ho (2022) use a pre-trained transformer to explore features and a convolutional neural network for further feature extraction. The comparison results are presented in Table 3. Our model achieves an Acc value of 84.42% and a MCC value of 83.72%. These indicators demonstrate that our model outperforms the other five methods in most aspects. Compared to the other models, our approach exhibits a remarkable improvement, with a maximum of 11.09% higher accuracy and 16.13% higher MCC value. The model designed by Le & Ho (2022) has the lowest MCC value among all methods, while BERMP and RFAthM6A show similar performance across various indicators. However, we note that the specificity (Sp) value of our method is slightly lower than that of BERMP and RFAthM6A, indicating a minor deficiency in predicting true negative samples. Nevertheless, overall, the experimental results clearly demonstrate the effectiveness of our model, which stands out as a superior approach for m6A site prediction on RNA sequences.