Injecting structure-aware insights for the learning of RNA sequence representations to identify m6A modification sites

Problem description and datasets

The task of identifying m6A modification sites is conventionally formulated as a binary classification problem. Each sample in a dataset represents a segment of the RNA sequence encompassing the modified central adenine and its adjacent nucleotides. Formally, the objective is to construct and train a model to map RNA sequences to a label set Y = {y}, where y ∈ {0, 1}. The fundamental objective is to effectively encode RNA sequences. This model takes the intelligible encoding of RNA sequences as input and strives to accurately assign labels to these representations, optimizing the binary cross-entropy loss function.

In this study, we leverage RNA sequence samples obtained from datasets of three distinct species as benchmark datasets for this classification task. Specifically, we choose the RNA m6A modification site dataset of Arabidopsis thaliana (A101 dataset) (Wan et al., 2015) for our investigation, comprising 2,100 samples with a balanced ratio of positive and negative samples set. To mitigate the impact of sample imbalance, we used a 1:1 balance of positive and negative samples, ensuring that the model learns equally from both classes, thereby enhancing its generalization ability. Negative samples were selected from the same gene regions as the positive samples to ensure biological consistency in the training data. In addition, negative samples were randomly selected from the same gene regions as the positive samples and are not required to conform to the DRACH motif. The A101 dataset consists of RNA sequences with a length of 101 nucleotides. This length provides sufficient coverage of the upstream and downstream sequence information around the m6A modification site while considering computational complexity. Additionally, we introduced two independent datasets, human-liver and mouse-heart, to further evaluate our model. The human-liver dataset contains 2,634 positive samples and 2,634 negative samples, while the mouse-heart dataset includes 2,201 positive samples and 2,201 negative samples. In both datasets, the sequence length is 41 nucleotides. The sequence lengths of the datasets are consistent with those used in previous studies (Dao et al., 2020; Rehman, Tayara & Chong, 2022), and all negative samples across the datasets are constructed from a central adenine and flanking nucleotides of equal length. This length is selected based on preliminary analysis and proven effective in capturing local contextual information while minimizing noise and redundancy from longer sequences. This approach generally helps to improve dataset quality by: 1) enhancing the effectiveness and reliability of computational recognition, validating the model’s ability to learn latent features for modification site identification; 2) strengthening the model’s generalizability, improving its applicability and robustness across different species, tissues, and experimental conditions.

Model description

M6A-SAI consists of three integral components that work in tandem to enhance m6A modification site prediction. The process begins with sequence representation initialization, where RNA sequences are encoded into embeddings using a Transformer encoder, capturing nucleotide types and their positional relationships. These embeddings are then used to construct a similarity graph, which is refined by the structure-aware optimization block to produce an optimized awareness graph, addressing limitations in the initial graph by focusing on key local regions and incorporating structural insights. Finally, the self-correlation fusion graph convolution framework processes both the similarity and awareness graphs through distinct graph convolution blocks to generate separate low-dimensional representations, which are then integrated by a fusion block to combine local and global perspectives, resulting in the final, comprehensive RNA sequence representations. The overall architecture and its components are illustrated in Fig. 1.

Sequence representation initialization

To capture more nuanced associations between RNA sequences, we integrate Transformer-based encoding with Manhattan distance metrics to construct a similarity graph that reflects the semantic relationships inherent in RNA sequences. Within the Transformer encoder framework, two core components are utilized: the positional encoding block, which incorporates positional information to differentiate nucleotide locations, and the self-attention mechanism, which enables the model to weigh and integrate dependencies between nucleotides across the sequence. Additionally, to delve deeper into the structural relationships within RNA sequences, we incorporate the calculation of Manhattan distance between vector representations of the output from the Transformer encoding block. Moreover, we rigorously constrain the values of the structure matrix within the set {0, 1}. The details of these components are outlined below.

Position encoding of Transformer is a functional encoder, that is, the position vector $\mathbf{p}=\left[p_i\right]$ can be computed as:

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

(2) $$w_i={\displaystyle \frac{1}{m^{2i/10\text{,}000}}}$$where $m$ is the number of sequences. Thus, the RNA sequences are represented by K-mer and position encoding. We obtain the sequence representation matrix ${\mathbf{R}}_s=\left[{\mathbf{r}}_{\mathbf{j}}\right]\in {\mathbb{R}}^{n\times d}$, where ${\mathbf{r}}_j$ can be computed as:

(3) $${\mathbf{r}}_j={\mathbf{r}}_{k-mer}+\mathbf{p}$$where, ${\mathbf{r}}_{k-mer}$ denotes the K-mer representations. After computing sequence representations, the encoder further optimizes the representation matrix ${\mathbf{R}}_s$through the self-attention mechanism:

(4) $$\{\begin{array}{c}\mathbf{Q}={\mathbf{W}}_q{\mathbf{R}}_s\hfill \\ \mathbf{K}={\mathbf{W}}_k{\mathbf{R}}_s\hfill \\ \mathbf{V}={\mathbf{W}}_v{\mathbf{R}}_s\hfill \\ {\mathbf{R}}_s=\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{x}\left({\displaystyle \frac{\mathbf{Q}{\mathbf{K}}^T}{\sqrt{d}}}\right)\mathbf{V}\end{array}$$where, $\mathbf{Q}$, $\mathbf{K}$ and $\mathbf{V}$ are the query matrix, key matrix and value matrix respectively. $\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{x}(\cdot )$ 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.

Structure-aware optimization block

We believe that it is a valid way to describe the associations of RNA sequences through internal relationships captured by the self-attention mechanism. To this end, we also measure this structural association by the Manhattan distance and construct the similarity graph $G_{s\mathrm{m}}=\left(V,E\right)$ for input sequences, where $V$ is the node set and E is the edge set. For an arbitrary edge $e_{ij}$ of node $v_i$ and $v_j$, we determine it as:

(5) $$\begin{array}{c}{e}_{ij}=\sum _{k=1}^n\left|{\mathbf{r}}_i^k-{\mathbf{r}}_j^k\right|\end{array}.$$

In addition, we ensure that the values of elements in $G_{sm}$ are strictly constrained in the set of $\left\{0,1\right\}$ by threshold filtering, as shown below:

(6) $$e_{ij}=\{\begin{array}{c}1if{e}_{ij}\le 0.5\hfill \\ 0otherwise.\hfill \end{array}$$

The similarity graph constructed using Manhattan distance may suffer from excessive smoothing, even with the application of threshold filtering to alleviate this issue. To address this limitation, we propose the integration of structure-aware blocks that are specifically designed to construct a perception graph. This approach enhances the similarity graph by incorporating constraints that emphasize local topological features, thereby improving the representation of structural nuances within the graph. The structure-aware block comprises three components: shortest path distance encoding and perceptual graph construction. Shortest path distance encoding is employed to learn potential global topology information of nodes within the similarity graph. The construction of the perception graph takes into account the importance of nodes to ascertain potential relationships between nodes in the perception graph. The specific descriptions of these two components are detailed as follows.

For an undirected graph, the shortest path distance (SPD) is one of the main ways to measure its topological properties. To this end, for a given similarity graph, we randomly select K nodes to construct the anchor set $\left\{{\mathcal{V}}_1,{\mathcal{V}}_2,\dots ,{\mathcal{V}}_K\right\}$. We further compute the SPD vector between all nodes in the similarity graph and this anchor set. For a given node ${\mathbf{v}}_i$, its SPD vector can be computed as:

(7) $$\begin{array}{c}{\mathbf{p}}_{v_i}=\left(\psi \left(v_i,{\mathcal{V}}_1\right),\psi \left(v_i,{\mathcal{V}}_2\right),\dots ,\psi \left(v_i,{\mathcal{V}}_K\right)\right)\end{array}.$$

Here, $\psi (\cdot )$ represents the computation process of a given node with the corresponding anchor. Taking $\psi \left(v_i,{\mathcal{V}}_K\right)$ as an example, it is calculated as follows:

(8) $$\begin{array}{c}\psi \left(v_i,{\mathcal{V}}_K\right)=\mathrm{m}\mathrm{i}{\mathrm{n}}_{{\mathcal{V}}_j\in {\mathcal{N}}_{v_i}}\left({\mathcal{P}}_{G_{sm}}\left(v_i,{\mathcal{V}}_j\right)\right)\end{array}$$where ${\mathcal{V}}_j\in {\mathcal{N}}_{v_i}$ denotes the existence of at least one path between node $v_i$ and anchor ${\mathcal{V}}_j$ in ${\mathbf{G}}_{sm}$. $\mathrm{m}\mathrm{i}\mathrm{n}(\cdot )$ denotes min-take operation.

After obtaining the SPD vectors of each node, combined with the output of the Transformer encoder, we reconstruct the awareness graph ${\mathbf{G}}_a$. This graph is used to re-measure the probability of a node establishing a link relationship. Specifically, take the node $v_i$ and $v_j$ as an example, we quantify this likelihood based on Eq. (9):

(9) $$\begin{array}{c}{G}_a\left(v_i,v_j\right)={\displaystyle \frac{1}{2}(z_i{z_j}^T+p_{v_j}{p_{v_j}}^T)\ast Rank}\end{array}$$where $W_h$ is a learnable weight matrix, $||$ is a concatenation operation, and $z_i$ and $z_j$ are the output representations of Transformer encoder of node $v_i$ and $v_j$, respectively. $Rank$ is a PageRank matrix.

Self-correlation fusion graph convolution framework

In this section, our objective is to acquire low-dimensional representations for nodes in both the similarity graph and perception graph, utilizing these representations for the identification of m6A modification sites. To achieve this, we propose the self-correlation fusion graph convolution framework, which comprises two main components: an self-correlation graph convolutional block for representation learning and a fusion block for integrating node representations from different graphs. The self-correlation graph convolutional block is designed to thoroughly capture the topological properties inherent in both the similarity and awareness graphs, independent of pre-existing node features. In contrast, the fusion block is utilized to integrate these low-dimensional node representations, which encapsulate topological information, thereby enhancing the quality of node representations for the subsequent classification task. Detailed descriptions of the framework are provided below.

We define the adjacent matrix $\mathbf{A}\in {\mathbb{R}}^{m\times m}$ for the similarity graph $G_{s\mathrm{m}}$, and define the node initial representation matrix as $\mathbf{E}\mathbf{r}\in {\mathbb{R}}^{m\times d}$. Each value in the matrix is determined and iteratively optimized by the network. $\mathbf{A}$ and $\mathbf{E}\mathbf{r}$ are input into a three-layer self-correlation graph neural networks to get the node embedding matrix ${\mathbf{R}}_{sm}$ of $G_{s\mathrm{m}}$:

(10) $$\begin{array}{r}\{\begin{array}{c}{\mathbf{R}}_{sm}^{\left(i+1\right)}={\mathbf{R}}_{sm}^{\left(i\right)}+{\zeta }^{\left(i+1\right)}Gcov{\left(\mathbf{A},\mathbf{E}\mathbf{r}\right)}^{\left(i+1\right)}\hfill \\ {\zeta }^{\left(i+1\right)}={\displaystyle \frac{1}{I+1}\left(i+1\right)}\hfill \end{array}\end{array}$$where, $\zeta$ 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:

(11) $$H^{\left(l+1\right)}=\delta {\left(\mathbf{A}\right)}^{-\frac{1}{2}}\left(\delta \left(\mathbf{A}\right)-{\displaystyle \frac{1}{2}\left(\mathbf{A}+{\mathbf{A}}^T\right)}\right)\delta {\left(\mathbf{A}\right)}^{-\frac{1}{2}}{H}^{\left(l\right)}{\mathrm{W}}^{\left(l\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)}=\mathbf{E}\mathbf{r}$. Similarly, we can compute the node embedding matrix for the awareness graph.

To attain a thorough understanding of RNA sequences, we employ a hybrid approach to integrate learned embedding representations. This methodology processes two distinct embeddings of the same sequence, addressing both local and global perspectives within the embedding space. We then establish weight relationships between the embeddings derived from similarity matrices and those obtained from sequence structure graphs. This is achieved through the implementation of specialized components: a local information extraction block and a global information extraction block. Each component is specifically designed to compute the weights associated with its respective embedding representation, thereby facilitating a more nuanced integration of local and global information.

Specifically, for a given similarity embedding matrix ${\mathbf{M}}^{sm}\in {\mathbb{R}}^{m\times d}$ and the corresponding awareness embedding matrix ${\mathbf{M}}^{aw}\in {\mathbb{R}}^{m\times d}$, we first obtain the weighted embedding matrix ${\mathbf{M}}^{we}\in {\mathbb{R}}^{m\times d}$ as:

(12) $${\mathbf{M}}^{we}={\mathbf{M}}^{sm}+{\mathbf{M}}^{aw}.$$

We then enter ${\mathbf{M}}^{we}$ into the local and global information extraction blocks, respectively. For the local and global information extraction block, the extraction process is described as follows:

(13) $$\{\begin{array}{c}{\mathbf{M}}_{lo}^{we}=f_{1D}\left(ReLU\left(f_{1D}^{\mathrm{\prime }}\left({\mathbf{M}}^{we}\right)\right)\right)\hfill \\ {\mathbf{M}}_{gl}^{we}=f_{1D}\left(\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left(f_{1D}^{\mathrm{\prime }}\delta \left({\mathbf{M}}^{we}\right)\right)\right)\hfill \end{array}$$where ${\mathbf{M}}_{lo}^{we}$ and ${\mathbf{M}}_{\mathrm{g}\mathrm{l}}^{\mathrm{w}\mathrm{e}}$ denote the output of the local extraction block and that of the global extraction block respectively. $f_{1D}$ and $f_{1D}^{\mathrm{\prime }}$ represent a one-dimensional convolution layer containing normalized functions respectively. 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 obtain the fusion representation matrix as:

(14) $$\{\begin{array}{c}{\mathbf{M}}_{final}^{we}={\mathrm{W}}_p{\mathbf{M}}^{sm}+\left(I-{\mathrm{W}}_p\right){\mathbf{M}}^{aw}\hfill \\ {\mathrm{W}}_p={\psi }_{so}\left({\mathbf{M}}_{lo}^{we}+{\mathbf{M}}_{gl}^{we}\right).\hfill \end{array}$$

Modification site identification

In this article, we use a support vector machine (SVM) to classify the final embedding matrix, a process represented by $p_i=\mathrm{S}\mathrm{V}\mathrm{M}\left({\mathbf{M}}_{final}^{we}\right)$, where $p_i$ is a identification score. In classification problems, SVM identifies a hyperplane to divide data points of different classes, aiming to maximize the margin between the hyperplane and the support vectors for enhanced generalization performance. SVM’s adaptability to high-dimensional spaces is a notable strength, achieved through the introduction of kernel functions that map data into higher-dimensional spaces, allowing the algorithm to address nonlinear classification problems. The algorithm excels in handling small sample sizes, high-dimensional datasets, and cases where data is not linearly separable.

Metrics for experiments

For a binary classification problem, the prediction outcomes can be categorized into four distinct groups: true positive (TP), false positive (FP), true negative (TN), and false negative (FN). To comprehensively assess the predictive performance, we employ several key evaluation metrics, including accuracy (Acc), F1-score, precision (Prec), sensitivity (Sn), specificity (Sp), and Matthews correlation coefficient (MCC). Leveraging these metrics, we construct both the receiver-operating characteristic (ROC) curve and the precision-recall (PR) curve to visualize the prediction results. Additionally, we compute the area under the PR curve (AUPR) and the area under the ROC curve (AUC) to provide quantitative assessments of our model’s performance. We evaluate our model using ten-fold cross-validation. Specifically, the dataset is divided into 10 approximately equal subsets, with each subset used as a validation set in turn, while the remaining nine subsets are used for training. This process repeats 10 times. This method is particularly important for small datasets, as it allows each data point to be used for both training and validation, maximizing data utilization and reducing evaluation bias due to limited data.

Performance on A101 datasets

In this article, we conduct a comprehensive evaluation of our model through 10-fold cross-validation. Specifically, we partition the dataset into 10 folds, employing nine folds for training and reserving one fold for validation in each iteration. This process is repeated ten times to ensure the validation of each fold, with a crucial consideration that validation set data do not overlap with the training set data. Six metrics, including AUPR, AUC, Acc, F1 score, Prec, and Sen, are employed to quantify the model’s performance. The results, detailed in Fig. 2 and Table 1, reveal promising outcomes. PR and ROC curves exhibit consistently good performance across all folds, with an average AUC of 91.96% and AUPR of 90.03%. This underscores the model’s strong overall performance, further supported by an MCC value of 73.0%. Additionally, other indicators such as Acc, F1, Prec, and Sen consistently exceed 80%. Notably, the Acc value of 86.1% signifies superior recognition performance in positive samples, highlighting the model’s efficacy in identifying potential positive instances. To emphasize the impact of structural perception, we visualize the scatter distribution of the validation set based on tSNE (Van der Maaten & Hinton, 2008), as depicted in Fig. 3. These experimental results demonstrate that the integration of structural information enhances the model’s ability to recognize m6A modification sites.

Performance comparison with other predictors

The comparative evaluation of three methods, namely M6AMRFS, BERMP, and RFAthM6A, is conducted through 10-fold cross-validation on the A101 dataset. Specifically, M6AMRFS utilizes two feature descriptors, namely dinucleotide binary coding and local site-specific dinucleotide frequency, to encode RNA sequences. It employs the F-score algorithm in conjunction with sequence forward search (SFS) to enhance feature representation and utilizes XGBoost (Chen & Guestrin, 2016) as a downstream classifier. BERMP utilizes gated recurrent units (GRU) (Dey & Salem, 2017) to represent RNA sequences and adopts an end-to-end training process for site identification. RFAthM6A attempts to classify various features from RNA sequences using machine learning methods. The comparison results, presented in Table 2 indicate that the proposed model outperforms the other three models across most indicators. The evaluation values for Acc and MCC reach 86.12% and 73.00%, respectively. This suggests the superiority of the proposed model. Notably, compared with the second method, the proposed model achieves a 5.07% improvement in Acc and a 0.45% improvement in MCC. M6AMRFS exhibits the lowest value in the MCC, while BERMP and RFAthM6A demonstrate similar values across all indicators. The experimental results affirm the effectiveness of the proposed model.

Ablation analysis

In this study, to assess the efficacy of the employed components, we introduce two additional variants and conduct validation using our model on the A101 dataset. The variants are labeled as M6A-SAI without Transformer (w/o Tr) and M6A-SAI without graph optimization (w/o Go). Specifically, w/o Tr employs K-mer representation in lieu of Transformer, while w/o Go omits the structure-aware graph optimization component. The comparison results, detailed in Fig. 4A, reveal a substantial decrease in all five indicators (AUC, AUPR, Acc, F1, and Prec) for both M6A-SAI and its variants. The degree of decrease is nearly identical for both variants, with the maximum decrease in AUC reaching 41.99% and the maximum decrease in Acc being 35.41%. This underscores the pivotal roles played by both the Transformer and the graph optimization process in the recognition task. To further underscore this point, we generated a scatter plot comparing our model with the two variants based on tSNE (Van der Maaten & Hinton, 2008). The results, displayed in Figs. 4B–4D, indicate that the scatter patterns of the two variants are nearly indistinguishable. This suggests the presence of a considerable number of false positive and false negative samples in the recognition results for both variants.

Furthermore, we discuss the individual contributions of each component. Specifically, we compare M6A-SAI with the model that solely utilizes fully connected and convolutional layers. We denote the model employing only fully connected layers as M6A-SAI-MLP, while the one utilizing graph convolutional layers is referred to as M6A-SAI-GCN. We conducted ten independent runs for M6A-SAI variants under different random seeds. As depicted in Table 3, Tables S1 and S2, we find that M6A-SAI-MLP exhibits inferior performance compared to both M6A-SAI and M6A-SAI-GCN on three datasets. It is probably because M6A-SAI-MLP solely extracts information from the feature matrix without capturing neighbor information adequately. Additionally, we observe a significantly superior performance of M6A-SAI over M6A-SAI-GCN. Although M6A-SAI-GCN can learn graph structure information, it is susceptible to noise present in the initial graph which hampers effective information aggregation. Conversely, by incorporating a graph optimization mechanism, M6A-SAI filters out this noise thereby enhancing node representation and improving prediction performance.

Analysis of the number of layers in graph neural networks

In this study, we establish the number of network layers for the graph convolution as 4. To elucidate the rationale behind this decision, we construct variants of the model with different numbers of network layers in the graph convolution as the backbone network. The comparative results of these variants with the model are depicted in Fig. 5. It is noteworthy that by considering the variant with four network layers as the critical point, we observe that when the number of layers is less than four, the model’s performance metrics increase with the growing number of network layers. However, when the number of layers exceeds four, the model’s metrics begin to decline with an increasing number of layers. Notably, with six layers, the model registers a reduction in AUC. We posit that this phenomenon is attributable to the initial stages, where augmenting the number of network layers facilitates the model in learning more intricate and abstract features, thereby enhancing its ability to represent data. However, as the number of network layers continues to increase, the model’s parameters also escalate, amplifying the risk of overfitting. In scenarios with limited training data and an excessively complex model, this may lead to favorable performance on the training set but poorer generalization to unseen data. Consequently, selecting an optimal number of network layers is crucial to strike a balance between feature representation capacity and the risk of overfitting.

Analysis of the type of graph neural networks

In our model, we adopt the simplest graph convolutional network (GCN) as the backbone network, primarily for the sake of model simplification. To emphasize this choice, we further explored three variants of graph neural networks (GNNs), including ChebNet (Tang, Li & Yu, 2019), graph attention network (GAT) (Velickovic et al., 2017), and GraphSAGE (Hamilton, Ying & Leskovec, 2017), as potential replacements for GCN in the backbone network. The comparison results, illustrated in Fig. 6 and detailed in Table 4, reveal that while different variants exhibit certain effects within our model, none outperforms GCN across all metrics. It is noteworthy that the other variants generally employ more intricate network structures, such as GAT, which introduces an attention mechanism, thereby increasing computational costs without yielding significant improvements in results. Despite their complexity, these variants do not surpass the performance of the simpler GCN in the context of m6A modification site identification. In summary, the results indicate that GCN stands out as the most suitable GNN type for our model in the task of m6A modification site identification, striking a balance between effectiveness and model simplicity.

Performance on independent data

To evaluate the generalization M6A-SAI, we adopt the datasets with different data size from different species. Specifically, we compare M6A-SAI with M6AMRFS, BERMP, RFAthM6A, AdaRM (Song et al., 2023), iRNA-m6A (Dao et al., 2020) and TS-m6A-DL (Abbas et al., 2021). Table 5 demonstrates that M6A-SAI outperforms BERMP (the second) on human-liver dataset, achieving an Acc that is 0.9% higher, MCC that is 3% higher, Sn that is 1.2% higher and Sp that is 3.9% higher. Table 6 demonstrates that M6A-SAI outperforms TS-m6A-DL (the second) on mouse-heart dataset, achieving an Acc that is 3.7% higher, Mcc that is 7.4% higher, Sn that is 2.9% higher, and Sp that is 4.6% higher. Besides, we notice that BERMP achieves satisfactory prediction performance on mouse data but worse performance on human-liver dataset, but M6A-SAI exhibits remarkable performance on all datasets, which illustrate M6A-SAI has good generalization compared with other methods.

Analysis of the overfitting of M6A-SAI

To illustrate the effectiveness of M6A-SAI in an intuitive fashion, we visualize learning curve to validate that there is overfitting in training process. As shown in Fig. 7, we find that with the increase in sample size, the performance of M6A-SAI on the test set exhibits gradual improvement and eventually stabilizes. Furthermore, there is no significant disparity in performance between the M6A-SAI trained on the training set and tested on the test set, indicating that despite operating on a small dataset, there is not overfitting in training process.

Performance on single-cell datasets

To further evaluate the generalization performance of M6A-SAI, we apply M6A-SAI to single-cell datasets. Despite the original design of M6A-SAI not being intended for such data, we believe that M6A-SAI can still be utilized for predicting m6A modification sites at the single-cell level. Specifically, we employ the dataset provided by Scm6A and followed the data processing approach of Scm6A to obtain the model input. Subsequently, M6A-SAI is employed for predicting m6A modification sites. As shown in Fig. 8, although M6A-SAI does not surpass Scm6A (RF) in terms of AUC, it outperforms Scm6A (KNN), Scm6A (LinearSVR), and Scm6A (LR) significantly. It is probably because that M6A-SAI is not originally designed for single-cell data, which often contains numerous sparse values, impacting the performance of M6A-SAI.