HD-6mAPred: a hybrid deep learning approach for accurate prediction of N6-methyladenine sites in plant species

Datasets

High-quality datasets are essential for model construction and evaluation. For the convenience of the comparison, we utilized the same datasets which are previously employed in 6mA prediction (6mA-StackingCV (Huang, Huang & Luo, 2023) and i6mA-vote (Teng et al., 2022). The datasets include: the Rosaceae training dataset, the Rosaceae testing dataset, the Rice dataset, and the Arabidopsis dataset. In which, the Rosaceae training dataset is employed to train the model and debug parameters, and the other three datasets, also named independent testing sets, are utilized to test the performance of the model. The Rice dataset is built by Lv et al. (2019), whereas the Rosaceae and Arabidopsis datasets are curated by Hasan et al. (2021). In summary, the Rosaceae training dataset contains 29,237 positive samples (6mA site-containing sequences) and 29,433 negative samples (non 6mA site-containing sequences), and the Rosaceae testing dataset consists of 7,298 positive samples and 7,300 negative samples. The numbers of positive samples and negative samples are 153,635 and 153,629 in the Rice dataset, respectively; and they are 31,414 and 31,843 in the Arabidopsis dataset. Each sequence is 41 bp long, with nucleobase A located at position 21. The information of these datasets has been given in 6mA-StackingCV (Huang, Huang & Luo, 2023), and the specific data can be downloaded from https://github.com/Xiaohong-source/6mA-stackingCV provided in the article.

Feature extraction

Feature selection is vital to the effectiveness of a predictive model. We explored four different feature encoding ways including one-hot encoding (Alipanahi et al., 2015; Kelley, Snoek & Rinn, 2016; Zhou & Troyanskaya, 2015), EIIP (Alakuş, 2023), NCP (Chen et al., 2015) and ENAC (Huang et al., 2018), which are widely used to represent DNA sequences for predicting 6mA sites.

In one-hot encoding, each nucleotide is encoded as a 4-dimensional binary vector, i.e., nucleotides A, C, G, and T, respectively denoted as A = [1, 0, 0, 0], C = [0, 1, 0, 0], G = [0, 0, 1, 0], and T = [0, 0, 0, 1] (Alipanahi et al., 2015; Kelley, Snoek & Rinn, 2016; Zhou & Troyanskaya, 2015). EIIP expresses nucleotides according to their electron–ion energy distribution in a given DNA sequence, and the four nucleotides are respectively represented as A = 0.1260, C = 0.1340, G = 0.0806 and T = 0.1335 (Alakuş, 2023). NCP encodes nucleotides based on their chemical properties, and the four nucleotides are expressed as A = [1, 1, 1], C = [0, 1, 0], G = [1, 0, 0] and T = [0, 0, 1], respectively (Chen et al., 2015). ENAC is the frequency of occurrence of a nucleotide in a sequence starting from the 5′ end to the 3′ end, determined through

$$\begin{array}{rl}\mathrm{V}=& [{\displaystyle \frac{{\mathrm{N}}_{\mathrm{A},\mathrm{w}\mathrm{i}\mathrm{n}1}}{\mathrm{S}},{\displaystyle \frac{{\mathrm{N}}_{\mathrm{C},\mathrm{w}\mathrm{i}\mathrm{n}1}}{\mathrm{S}},{\displaystyle \frac{{\mathrm{N}}_{\mathrm{G},\mathrm{w}\mathrm{i}\mathrm{n}1}}{\mathrm{S}},{\displaystyle \frac{{\mathrm{N}}_{\mathrm{T},\mathrm{w}\mathrm{i}\mathrm{n}1}}{\mathrm{S}},...,{\displaystyle \frac{{\mathrm{N}}_{\mathrm{A},\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{L}-\mathrm{S}+1}}{\mathrm{S}},{\displaystyle \frac{{\mathrm{N}}_{\mathrm{C},\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{L}-\mathrm{S}+1}}{\mathrm{S}},}}}}}}\\ & {\displaystyle \frac{{\mathrm{N}}_{\mathrm{G},\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{L}-\mathrm{S}+1}}{\mathrm{S}},{\displaystyle \frac{{\mathrm{N}}_{\mathrm{T},\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{L}-\mathrm{S}+1}}{\mathrm{S}}}}]\end{array}$$where S denotes the width of the sliding window (Huang et al., 2018) and it equals five in this study, ${\mathrm{N}}_{\mathrm{M},\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{j}}$ represents the occurrence number of nucleotide M in the window $\mathrm{j}$, $\mathrm{M}\in \{$A, C, G, T}, j = 1, 2, …., L-S+1.

To obtain the optimal combination of features, we implemented a hold-out search strategy similar to that used in 6mA-StackingCV (Huang, Huang & Luo, 2023). This strategy can test the performance of various feature combinations in distinguishing 6mA sites from non-6mA sites. Unlike 6mA-StackingCV (Huang, Huang & Luo, 2023), which aims to identify better features within a conventional ML framework, we employed BiGRU and CNN as our learning algorithms to achieve improved performance with the deep learning paradigm. To illustrate the generalization of selected features, we only carried out a series of experiments on the Rosaceae training dataset. The Rosaceae training dataset was split into the hold out-training dataset and hold out-testing dataset as a ratio of 8:2. First, we computed the ACC of the BiGRU using a single feature encoding and selected the feature with the largest ACC value as the first baseline feature, named as the first-order encoding. Second, we combined the first-order encoding with the remaining three single encoding ways. The increasing ACC values mean that the added feature encoding contributed positively to the performance. If there exist ACC values that are more than that of the first-order encoding, then the combination having the highest ACC score is chosen as the updated encoding, denoted as the second-order encoding. The process continued, adding remaining single encoding ways to the current combination, until no further ACC improvements were observed. The feature combination with the highest ACC value is used to construct the 6mA sites predictor. The strategy for selecting the optimal encoding ways under the CNN is the same as above.

Deep learning approach

As illustrated in Fig. 1, in constructing the model, we mainly employed three deep learning modules: BiGRU, CNN and attention mechanism. As an improved form of recurrent neural network (RNN), GRU can address the vanishing gradient problem often encountered in the processing of long sequences and has demonstrated excellent performance in addressing sequential challenges. Conversely, CNN excels at achieving parallel computation and capturing local dependencies in data (Wang et al., 2024). To effectively integrate the outputs from both networks, we employed an attention mechanism that fuses features extracted by the BiGRU and CNN. Firstly, the sequences were fed to the BiGRU and CNN modules according to their corresponding optimal encoding ways. Subsequently, the attention mechanism assigns weights to the features extracted by both the BiGRU and CNN, highlighting important features and filtering out less relevant information. Finally, the features refined by the attention mechanism were processed through a fully connected layer (Dense) and a sigmoid activation function, resulting in a score between 0 and 1 in identifying the DNA 6mA sites. The details of the methods involved in the model were listed follows.

BiGRU: In DNA sequence data, the BiGRU module is designed to capture both forward and backward dependencies (Nguyen-Vo et al., 2023; Jia et al., 2024). It is composed of forward and backward GRUs, which read the sequences from two opposite directions. Special gating mechanisms regulate the retention and forgetting of information from previous states. These mechanisms generate candidate hidden states, which are then used to calculate the current hidden state. Finally, the outputs of the forward and backward GRUs are concatenated along the feature dimension to obtain a comprehensive feature representation.

CNN: CNN utilizes convolution and pooling operations to extract essential local features. Through the introduction of non-linearity using the ReLU activation function, the network can efficiently reduce data dimensions via pooling. In this study, the computation can be expressed as follows Eqs. (1)–(3):

(1) $$\mathrm{Y}\left[\mathrm{i}\right]=\sum _{\mathrm{j}}\mathrm{W}\left[\mathrm{j}\right].\mathrm{X}\left[\mathrm{i}+\mathrm{j}\right]$$ (2) $$\mathrm{R}\mathrm{E}\mathrm{L}\mathrm{U}\left(\mathrm{x}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left(0,\mathrm{x}\right)$$ (3) $${\mathrm{Y}}_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}}\left[\mathrm{i}\right]=\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{X}\left[\mathrm{i}:\mathrm{i}+\mathrm{k}\right]\right)$$where $\mathrm{X}$ is the input sequence, $\mathrm{W}$ denotes the convolution kernel, $\mathrm{Y}\left[\mathrm{i}\right]$ represents the convolution output, $\mathrm{W}\left[\mathrm{j}\right]$ is the convolution weight, $\mathrm{X}\left[\mathrm{i}+\mathrm{j}\right]$ is the local region of the input sequence. k is the size of the pooling window, and ${\mathrm{Y}}_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}}\left[\mathrm{i}\right]$ is the maximum pooling. While the max function selects the maximum value from the defined window, effectively reducing feature dimensions.

Attention mechanism: The essence of the attention mechanism is to calculate the values representing the significance of different components. It comprises three components: Query (Q), Key (K) and Value (V). The similarity of Q and K is usually quantified through dot products:

(4) $$\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(Q,K\right)=\mathrm{Q}\cdot {\mathrm{K}}^{\mathrm{T}}$$ (5) $$\mathrm{A}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}=\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{x}\left({\displaystyle \frac{\mathrm{Q}\cdot {\mathrm{K}}^{\mathrm{T}}}{\sqrt{{\mathrm{d}}_{\mathrm{k}}}}}\right)\mathrm{V}$$ (6) $$\mathrm{Q}={\mathrm{W}}_{\mathrm{Q}}\cdot {\mathrm{h}}_{\mathrm{G}\mathrm{R}\mathrm{U}}$$ (7) $$\mathrm{K}={\mathrm{W}}_{\mathrm{K}}\cdot {\mathrm{h}}_{\mathrm{C}\mathrm{N}\mathrm{N}}$$ (8) $$\mathrm{V}={\mathrm{W}}_{\mathrm{V}}\cdot {\mathrm{h}}_{\mathrm{C}\mathrm{N}\mathrm{N}}$$where ${\mathrm{d}}_{\mathrm{k}}$ denotes the size of the K vector, $\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{x}$ is used to normalize the weights. ${\mathrm{W}}_{\mathrm{Q}}$, ${\mathrm{W}}_{\mathrm{K}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{W}}_{\mathrm{V}}$ are the learned weight matrices. In our model, the Query originates from the output of BiGRU, while the Key and Value are derived from the output of CNN.

Output layer: The ultimate result is generated by a dense layer that maps the outcomes of the attention mechanism to the final predictions as Eq. (9).

(9) $${\mathrm{y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=\mathrm{\sigma }\left({\mathrm{W}}_{\mathrm{o}}\cdot \mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}+{\mathrm{b}}_{\mathrm{o}}\right)$$where ${\mathrm{W}}_{\mathrm{o}}$ represents the weight matrix in the output layer, ${\mathrm{b}}_{\mathrm{o}}$ means the bias term, and $\sigma$ denotes the sigmoid function.

Loss function: The model employs binary cross entropy (Ma, Liu & Qian, 2004) as its loss function, appropriate for binary classification tasks, and the formula is as follows:

(10) $$\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}=-{\displaystyle \frac{1}{\mathrm{N}}\sum _{\mathrm{i}=1}^{\mathrm{N}}\left[{\mathrm{y}}_{\mathrm{i}}\log \left({\mathrm{y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d},\mathrm{i}}\right)+\left(1-{\mathrm{y}}_{\mathrm{i}}\right)\log \left(1-{\mathrm{y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d},\mathrm{i}}\right)\right]}$$where ${\mathrm{y}}_{\mathrm{i}}$ is the true label and ${\mathrm{y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d},\mathrm{i}}$ is the predicted probability.

Detail procedure and hyperparameters setting

In this section, we give the detail procedure and hyperaramters about HD-6mAPred.

For the BiGRU module, the number of BiGRU layers is set to 2, with 64 units per layer. The activation function is tanh, and the gated activation function is sigmoid. The dropout rate is 0.5.

For the CNN module, the number of convolutional layers is 2, with 128 filters per layer. The kernel size is 3, and the activation function is the rectified linear unit (ReLU). After each convolutional layer, the max pooling layer has a pooling size of 2. After the two convolutional operations, a flatten layer converts the multi-dimensional output into a one-dimensional vector, with an additional dropout rate of 0.5.

The outputs of the two subnetworks are connected through the concatenate layer, and the attention mechanism is used for feature weighting. The dimension of query, key, and value in the attention mechanism is 16. The output of attention is processed through layer normalization, and the final output is obtained from a fully connected layer (Dense), whose activation function is sigmoid.

For model optimization, the optimizer is Adam, and the initial learning rate is 0.0001. We utilized a 5-fold cross-validation method for training and testing the model. During each fold’s training phase, we set the maximum number of epochs to 100 and the batch size to 64. An early stopping strategy was implemented to monitor the validation set loss. Training ceases if there is no improvement after five consecutive epochs, and the model reverts to the optimal weights obtained.

Evaluation metrics

To test the performance of HD-6mAPred, we implemented a 5-fold cross-validation. Specifically, all samples were partitioned into five groups randomly, the four of which were used to train the model, while the remaining one was designated for testing the model. This procedure is repeated five times to ensure that each sample serves as a testing set at least once. The final testing results were obtained by averaging the metrics from the five iterations. Four standard evaluation metrics consisting of accuracy (ACC), Matthews correlation coefficient (MCC), sensitivity (SN) and specificity (SP) are employed for evaluating the performance of HD-6mAPred. These metrics can be formulated by Eqs. (11)–(14).

(11) $$\mathrm{A}\mathrm{C}\mathrm{C}={\displaystyle \frac{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}+\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}}$$ (12) $$\mathrm{M}\mathrm{C}\mathrm{C}={\displaystyle \frac{\mathrm{T}\mathrm{P}\times \mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}\times \mathrm{F}\mathrm{N}}{\sqrt{\left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}\right)\left(\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{N}\right)\left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}\right)\left(\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}\right)}}}$$ (13) $$\mathrm{S}\mathrm{N}={\displaystyle \frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}}$$ (14) $$\mathrm{S}\mathrm{P}={\displaystyle \frac{\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}}}$$where, TP (true positives) and TN (true negatives), respectively denote the number of the correctly predicted 6mA and non-6mA samples; FP (false positives) and FN (false negatives) are the number of the falsely predicted non-6mA and 6mA samples, respectively. In addition, we utilized the receiver operating characteristic (ROC) curves to visually test the prediction performance of HD-6mAPred. AUC is computed as the area beneath the ROC curve.

Selection of the optimal feature combination

We conducted a series of five-fold cross-validations on the Rosaceae training set using both the BiGRU and CNN models to identify the optimal feature combinations. ACC and AUC values of various features and their corresponding combinations were compared.

As presented in Table 1, the BiGRU model demonstrates that ENAC encoding yields the highest ACC of 0.903 and AUC of 0.957 for a single feature. Consequently, we selected ENAC encoding as our first-order feature. Subsequently, we examined the combination of ENAC encoding with the other three feature encoding ways: one-hot, EIIP and NCP. Among the combinations, ENAC combined with one-hot (denoted as ENAC+one-hot) achieved the highest ACC of 0.914 and AUC of 0.831. Accordingly, ENAC+one-hot is selected as the optimal second-order feature. Further investigations into the combination of ENAC+one-hot with other single features illustrates that ENAC+one-hot+NCP has the highest ACC of 0.920 and AUC of 0.860, but those of them in ENAC+one-hot+NCP+EIIP are descreased. Therefore, the third-order feature ENAC+one-hot+NCP is selected as the optimal feature for the BiGRU model.

For the CNN model, as shown in Table 2, NCP encoding provided the highest ACC (0.940) and AUC (0.979) when considered as a single feature. Therefore, we chose NCP encoding as our first-order feature. We then explored the combination of NCP with the other three encodings: one-hot, EIIP and ENAC. Among the combinations, NCP combined with ENAC (denoted as NCP+ENAC) achieved the highest ACC and AUC values. Accordingly, NCP+ENAC is selected as the optimal second-order feature. Further investigations into the combination of NCP+ENAC with other single features revealed that when combining all three encoding ways, the highest ACC value achieved is 0.942, which is lower than the 0.943 obtained from NCP+ENAC. Therefore, NCP+ENAC is chosen as the optimal feature combination for the CNN model.

Motif analysis

To analyze position-specific differences in 6mA and non-6mA-containing sequences for Rosaceae, rice and Arabidopsis, the two-sample logo software (Vacic, Iakoucheva & Radivojac, 2006) at a level of p < 0.05 was used. As shown in Fig. 2, the ‘A’ nucleotide is at the 21st position in 41-nt DNA sequences.

Figure 2 reveals distinct over-representation (enrichment) and under-representation (depletion) of nucleotides at specific positions. For Rosaceae, the ‘A’ base is enriched at positions 12, 17–19, 25 and 28 but depleted at 22, 24, 2 and 30. The ‘G’ base is enriched at positions 22, 23 and 26, while depleted at 27 and 28. The ‘C’ base is over-represented at positions 24, 27 and 30 but under-represented at 20, 25, 28 and 29. In rice, the ‘A’ base is enriched at positions 15–18, 25, and 28. The ‘C’ base is more abundant at positions 19, 23, 27 and 30. The ‘G’ base is over-represented at 13, 20, 26 and 29 and under-represented at 27. For Arabidopsis, ‘A’ is enriched at positions 15–18, 20, 25, 29 and 32; ‘C’ at 19 and 27; ‘G’ at 20, 22–24, and 28; and ‘T’ is depleted at 14–20, 25–26, and 37–41. Given these sequences containing 6mA are enriched and depleted with some nucleotides at some positions, it is speculated that the sequence information contribute in discriminating the 6mA from non-6mA.

Model analysis

To verify the effectiveness of combining BiGRU, CNN, and an attention mechanism in improving predictive performance, we compared six different models using their corresponding optimal encoding ways on the independent testing sets of Rosaceae, rice and Arabidopsis. Among these models, one model is our proposed HD-6mAPred, the other five models are BiGRU, CNN, the combination of BiGRU and CNN (denoted as BiGRU+CNN), and the two models that incorporate attention mechanisms, referred to as BiGRU+Attention and CNN+Attention, respectively.

Table 3 presents a comprehensive performance comparison of these models. Notably, the HD-6mAPred exhibited superior performance across all testing species. Specifically, for Rosaceae, HD-6mAPred achieved an ACC value of 0.996, outperforming BiGRU by 0.481, CNN by 0.050, BiGRU+Attention by 0.414, CNN+Attention by 0.064, and BiGRU+CNN by 0.050. For rice, the ACC value for HD-6mAPred is 0.952, which is 0.014 higher than BiGRU, 0.016 higher than CNN, 0.015 higher than BiGRU+Attention, 0.024 higher than CNN+Attention, and 0.016 higher than BiGRU+CNN. In Arabidopsis, HD-6mAPred achieved an ACC of 0.937, exceeding BiGRU by 0.100, CNN by 0.074, BiGRU+Attention by 0.082, CNN+Attention by 0.085, and BiGRU+CNN by 0.074. Furthermore, the HD-6mAPred demonstrated significant improvements in MCC, SN and SP, indicating its robust capability to accurately identify both true positive and true negative 6mA sites. The ROC curves, illustrated in Fig. 3, further support these findings, with AUC values for HD-6mAPred being 1.000 for Rosaceae, 0.985 for rice, and 0.980 for Arabidopsis. These results demonstrate the effectiveness of employing a hybrid model that integrates BiGRU, CNN and the attention mechanism to enhance predictive performance for 6mA site identification.

Performance comparison with existing methods

In recent years, several excellent computational methods have been developed to identify 6mA sites. We conducted a comparison of our HD-6mAPred approach with these existing methods using 5-fold cross-validation across three independent testing sets: the Rosaceae testing set, the Rice testing set and the Arabidopsis testing set. Table 4 lists the comparison methods and results. We complied the code of CNN6mA (Tsukiyama, Hasan & Kurata, 2022) on our test datasets, and the values of evaluation metrics for other existing methods were directly cited from the results reported by i6mA-Caps (Rehman et al., 2022) or 6mA-StackingCV (Huang, Huang & Luo, 2023). As shown in Table 4, among all the comparison methods, 6mA-StackingCV and i6mA-Caps are the latest and the best methods, respectively, for predicting 6mA sites in the three plant species. Therefore, we mainly reported the comparison results with these two methods.

On the Rosaceae dataset, HD-6mAPred demonstrates the most excellent results, surpassing all other approaches. Its values of ACC, MCC, SN and SP are 0.996, 0.993, 0.995 and 0.998, respectively. Compared with 6mA-StackingCV, they are 0.036, 0.073, 0.036 and 0.037 higher, respectively; compared with i6mA-Caps, they are 0.029, 0.055, 0.029 and 0.03 higher, respectively.

On the Rice dataset, although the SN value (0.955) of HD-6mAPred is slightly lower than that of several methods, including Meta-i6mA, iDNA6mA-Rice, MM-6mAPred, 6mA-vote, 6mA-StackingCV and CNN6mA, and its MCC value (0.905) is lower than that of CNN6mA (0.995), HD-6mAPred performed best in the other two evaluation metrics. Its ACC and SP are 0.952 and 0.949, respectively. Compared with 6mA-StackingCV, these two metrics are 0.107 and 0.223 higher, respectively. In addition, compared with i6mA-Caps, HD-6mAPred increased the ACC by 0.012, the MCC by 0.025, the SN by 0.004, and the SP by 0.02.

On the Arabidopsis dataset, HD-6mAPred performed best in three of the four evaluation metrics, and its SP value (0.948) is almost the same as that of i6mA-Fuse_FV (0.949). The values of ACC, MCC, SN and SP of HD-6mAPred are 0.937, 0.875, 0.927 and 0.948, respectively. Compared with 6mA-StackingCV, these values are increased by 0.155, 0.299, 0.250 and 0.082, respectively; compared with i6mA-Caps, they are increased by 0.069, 0.135, 0.105 and 0.033, respectively.

To visually highlight the superiority of HD-6mAPred, we presented a performance comparison between HD-6mAPred and the latest methods 6mA-StackingCV, as well as the currently best methods i6mA-Caps by bar charts. The results are illustrated in Fig. 4. Overall, our method achieved the best performance in almost all evaluation metrics. These results not only demonstrate the excellent predictive capabilities of HD-6mAPred, but also evidence its strong generalization ability across different species.

In addtion, to estimate the performance of the model on imbalanced datasets, we randomly extracted 10% of the positive samples from each species while keeping the number of negative samples unchanged, constructing imbalanced datasets with a ratio of 1:10 for the three species, and then tested the performance of our model on these datasets through 5-fold cross-validation. And we obtained some intresting results. For example, for Rosaceae, ACC, MCC, SN, SP and AUC were 0.983, 0.966, 0.984, 0.982, and 0.997, respectively; for rice, the corresponding values were 0.954, 0.908, 0.958, 0.949, and 0.986, respectively. For Arabidopsis, the corresponding values were 0.928, 0.856, 0.915, 0.940, and 0.672, respectively. These results further demonstrate the robust generalization ability of the HD-6mAPred model.