To predict heart-related diseases using abnormal m6A modification sites, we obtained methylation sites from five heart disease conditions and their corresponding control groups using Methylated RNA Immunoprecipitation Sequencing (MeRIP-seq) from the Gene Expression Omnibus (GEO) (Barre et al., 2013) dataset (Table 1). We also obtained single-base resolution m6A modification sites by m6A-SAC-seq (Hu et al., 2022) through the m6A-Atlas database (Liang et al., 2024) (Supplementary Table S1). In this study, we defined the overlapping sites obtained from MeRIP-seq and single-base resolution sites as true m6A methylation modification sites. We defined abnormal m6A modification as two cases: m6A modification upregulation and m6A modification downregulation. Methylation modification sites present only in the experimental group were defined as m6A upregulation sites, while methylation modification sites present in both the experimental and control groups were defined as their corresponding control group sites. We defined methylation modification sites present only in the control group as m6A downregulated sites and sites that were not methylated in either the experimental or control groups as their corresponding control group sites. In total, we obtained 20 datasets for analysis. In addition, we retrieved datasets for three additional cardiovascular diseases—dilated cardiomyopathy, ischemic cardiomyopathy, and aortic dissection—from the GEO database (Table 1) to serve as independent validation cohorts, thereby demonstrating the generalizability and robustness of our model.

In our predictions, we utilized the m6A-upregulated sites of each pathological state and their corresponding control sites, m6A downregulated sites with their control group sites to train two separate prediction models. For example, AC16EV_down was trained using evodiamine-induced cardiac toxicity samples, where sites with no methylation (positive sites in the experimental group but not in the control group) and sites with no methylation in both the experimental and control groups (negative sites) were used. AC16EV_up was trained using evodiamine-induced cardiac toxicity samples, where positive sites were those with methylation in the experimental group but not in the control group, and negative sites were those without methylation in both groups. When training the predictor, we first randomly selected 1,000 samples from each of the positive and negative site datasets, with 80% of the samples used for training and the remaining 20% for independent testing. We trained the model using full transcriptomic data, considering both unmodified and methylated sites from exons and introns.

Machine learning algorithms possess powerful data processing and pattern recognition capabilities, rendering them widely applicable in the field of bioinformatics, particularly in the prediction of DNA, RNA, and protein modifications. This study utilized the R program (version 4.4.2) and related packages to train models for the two abnormal m6A modifications in each pathological state, and optimized the models to obtain the final predictors. The sequence length is critical for prediction accuracy. To select the optimal sequence length, we tested 21, 31,41, 51, 61, 71, and 81 nucleotides (nt) centered around the m6A modification sites. In addition to the sequence length, selecting an appropriate feature-encoding method is equally important. We employed six encoding methods and combined all 63 single or combined features for training and testing, ultimately selecting the best feature combination for further optimization.

To determine the most suitable algorithm for model construction, we systematically compared four machine learning algorithms currently popular in bioinformatics research, including Support Vector Machines (SVM, e1071 package (Dimitriadou et al., 2008), Random Forests (RF, randomForest (Liaw and Wiener, 2002),generalized linear models (GLM, stats package (R Core Team R, 2020), and efficient gradient boosting (XGBoost, xgboost package (Chen and Guestrin, 2016). We evaluated the performance of the predictors through independent testing and explored the impact of the parameter selection. The primary performance metric was the area under the receiver operating characteristic curve (AUROC). Additionally, we calculated the accuracy (ACC), sensitivity (Sn), and specificity (Sp) to provide a comprehensive comparison of the algorithm performance (Equation 2). To rigorously evaluate model’s performance on potentially imbalanced datasets, we concurrently employed the Area Under the Precision-Recall Curve (AUPRC). This metric provides a more comprehensive assessment of a classifier’s ability to accurately identify positive samples (true m6A sites) with high precision, complementing AUROC analysis. S n = T P T P + F N S p = T N F P + T N A C C = T P + T N T P + F P + T N + F N (2)

To further elucidate the driving mechanisms behind the model’s predictions, we employed the Mean Decrease Gini from a random forest to assess feature importance (Breiman, 2001). This method measures how much each feature reduces node impurity across all splits in the forest. Specifically, for each split using a given feature, the decrease in impurity from the parent node to the child nodes is recorded and summed over all splits. The resulting value reflects the overall contribution of the feature to reducing the uncertainty of the model. A key advantage of this method is its efficiency, as it leverages the existing structure of the random forest without requiring additional computation.

For our predictive model, we used the Mean Decrease Gini to analyze feature importance across different feature extraction methods. By calculating the Mean Decrease Gini value for each feature, we identified the features that contributed most significantly to the model’s classification decisions. This not only enhanced the interpretability of the machine learning model, but also provided biological insights into the most distinctive features of aberrant m6A sites under cardiac pathological conditions.

To visualize cardiac morphology in vivo, the transgenic zebrafish line Tg(myl7:GFP) (Catalog ID: CZ56), obtained from the China Zebrafish Resource Center (CZRC), was used in this study. Zebrafish were maintained under standard conditions at 28.5 °C with a 14 h light/10 h dark cycle. Healthy embryos were collected and randomly divided into the control group and the Evodiamine treatment group. At 24 h post-fertilization (hpf), embryos in the treatment group were exposed to Evodiamine at a concentration of 50 ng/mL, while the control group was treated with 0.1% DMSO. The cardiac morphology of zebrafish larvae was observed at 3, 4, and 5 days post-fertilization (dpf) using a stereomicroscope. Pericardial cavity area and cardiac area were quantified to assess the extent of cardiotoxicity. ImageJ software was used for image processing and quantitative analysis.

Total RNA was extracted from 30 zebrafish larvae at 4 days post-fertilization (dpf) per group using Trizol reagent (Invitrogen, Carlsbad, CA, United States) method. Subsequently, total RNA was reverse-transcribed into cDNA templates using the Evo M-MLV RT Reaction Mix kit (Accurate Biology, Changsha, China). The primer sequences specific to the m6A modification-related genes, including METTL3, METTL5, METTL14, METTL16, RBM15B, RBM15, VIRMA, ZC3H13, and WTAP, are listed in Table 2. RT-qPCR was performed using the AriaMX real-time PCR system (Agilent, Santa Clara, CA, United States). The thermal cycling conditions comprised an initial denaturation step at 95 °C for 30 s, followed by 40 cycles at 95 °C for 5 s and 60 °C for 30 s. A dissociation curve analysis was performed according to the instrument’s default settings. The relative mRNA expression levels of each targeted gene were normalized to ACTB2 and calculated using the 2^−ΔΔCT method.

Recent studies have demonstrated that sequence-derived features are reliable and effective in capturing the inherent specificity of target sequences. Therefore, we explored six different encoding methods to compare their effectiveness in predicting m6A modification specificity in different cardiac pathological states. In addition, we performed motif analysis of abnormal m6A modification sequences from five cardiac pathological states using MEME Suite (Bailey and Elkan, 1994). The “RRAC” motif is frequently observed in sequences exhibiting abnormal m6A modification in five distinct cardiac pathological states (Figure 2A). Consequently, it can be deduced that cardiac diseases resulting from aberrant m6A modification may be driven by shared factors.

To explore the correlation between m6A modification sequences, m6A modification and biological functions, we performed Gene Ontology (GO) enrichment analysis on https://www.bioinformatics.com.cn (last accessed on 10 December 2024), an online platform for data analysis (Tang et al., 2023). This approach enabled us to identify the biological processes associated with cardiac pathological states and reveal their potential roles in cellular functions. As illustrated in Figure 2B, the top 10 GO biological process terms associated with the two m6A abnormalities in each cardiac pathological state are shown. We can observe that different biological processes are significantly enriched in the two types of m6A abnormal modifications under various cardiac pathological conditions. For example, evodiamine-induced cardiac toxicity caused by abnormal downregulation of m6A is associated with RNA splicing and DNA replication; whereas m6A upregulation makes evodiamine-induced cardiac toxicity more closely related to histone modification and peptidyl-lysine modification. Similarly, matrine-induced cardiac toxicity caused by abnormal downregulation of m6A is associated with histone modification and negative regulation of the cell cycle, while abnormal upregulation of m6A makes matrine-induced cardiac toxicity more closely related to RNA splicing and chromatin organization. For heart calcification, the biological processes associated with abnormal m6A downregulation are closely linked to tRNA metabolic processes and tRNA processing, while abnormal m6A upregulation is associated with RNA splicing and RNA splicing via transesterification reactions. In studies on cardiac hypertrophy, biological processes induced by abnormal downregulation of m6A primarily involve histone modification and peptidyl-lysine modification, while abnormal upregulation of m6A is associated with ncRNA metabolic processes and histone modification. Finally, TKI-induced cardiotoxicity caused by abnormal downregulation of m6A is associated with histone modification and peptidyl-lysine modification, while abnormal upregulation of m6A is more closely related to regulation of response to DNA damage stimulus and histone modification.

The amount of sequence information captured varies with the length of the input sequence window, and the selection of window size has a direct impact on the performance of the predictive models (Chen et al., 2021; Chen et al., 2015). To ensure optimal model performance, we systematically evaluated input sequences of varying lengths (21, 31, 41, 51, 61,71, and 81 nt), each centered on the m6A modification site (Figure 3). In the full transcriptome model, as the sequence length increased, the predictive performance for different cardiac pathologies also improved, as evidenced by the increase in the AUROC score. However, beyond a certain length, further increases in the sequence length lead to a decline or stabilization of the predictive performance. Based on these observations, we selected the optimal nucleotide sequence length for each prediction model (based on AUROC scores) to generate RNA sequence features for different pathological states. For AC16EV_down, AC16EV_up, AC16Mat_down, AC16Mat_up, Calcific_down, Calcific_up, hypertrophy_down, hypertrophy_up, TKIs_down, and TKIs_up, we selected sequence lengths of 31 nt, 41 nt, 51 nt, 41 nt, 41 nt, 71 nt, 21 nt, 41 nt, 81 nt, and 71 nt, respectively, to generate RNA sequence features with abnormal m6A modifications under different pathological conditions for the model training.

We trained the prediction models using m6A sites that were relatively upregulated or downregulated in each cardiac pathological condition and validated their performance on independent test sets. The optimal feature combinations were identified by calculating the AUROC values across the models. We considered all possible feature combinations derived from the six encoding schemes, including single-feature sets for each pathological state. The model performance (AUROC) for various feature combinations was evaluated, and the results showed that different models excelled when paired with different feature sets. Consequently, for each model, we selected the most suitable feature combination to establish the initial m6AHD (Figure 4).

Many machine learning algorithms have gained widespread recognition and application in the field of RNA modification prediction because of their excellent predictive accuracy and generalization capabilities (Chen et al., 2021; Chen et al., 2019a; Chen et al., 2019b; Liu et al., 2020b; Huang et al., 2018). To investigate which machine learning algorithm is more suitable for our project to predict different cardiac pathological states, we selected four algorithms—Support Vector Machine (SVM), Random Forest (RF), Generalized Linear Model (GLM), and efficient extreme gradient boosting (XGboost)—for a systematic comparison (Figure 5). We primarily evaluated the performance of the predictive models by calculating the area under the receiver operating characteristic curve (AUROC) on an independent test dataset. We also examined accuracy, sensitivity and specificity. Additionally, we assessed the area under the precision-recall curve (AUPRC) to validate model robustness. In line with the AUROC results, the Random Forest algorithm achieved the highest AUPRC scores for most pathological conditions, outperforming Support Vector Machines (SVMs), Generalised Linear Models (GLMs) and XGBoost. This suggests that the RF-based m6AHD framework strikes an excellent balance between precision and recall, effectively reducing the false positive rate. In summary, the Random Forest algorithm demonstrated the most stable and optimal performance when optimised sequence lengths and feature combinations were employed.

The parameter settings of the Random Forest model significantly influenced its predictive performance. Two key hyperparameters commonly adjusted are the number of trees (ntree) and the number of features considered at each split (mtry). The ntree parameter specifies the number of decision trees in the ensemble; more trees typically improve stability and reduce variance, although beyond a certain point, the gains diminish and computation becomes inefficient. In contrast, mtry controls the number of candidate predictors randomly selected at each node; lower values increase the diversity among trees and help prevent correlation (thus reducing variance), while higher values allow stronger predictors to be consistently chosen, potentially improving accuracy but risking overfitting. Thus, the careful tuning of ntree and mtry is essential to strike an optimal balance between bias, variance, and computational cost in Random Forest modeling.

In this study, we combined all parameter combinations for the ntree parameter ranging from 200 to 1,400 (in increments of 100) and the mtry parameter ranging from 50 to 100 (in increments of 1), and evaluated the predictive performance of the Random Forest models by calculating AUROC values (Figure 6). Based on the resulting AUROC values, the most appropriate ntree–mtry combination was determined. As a result, we finally determined the optimal predictive model for abnormal m6A modification for each cardiac pathological condition. The AUROC scores for the following model were determined: AC16EV_down, AC16EV_up, AC16Mat_down, AC16Mat_up, Calcific_down, Calcific_up, hypertrophy_down, hypertrophy_up, TKIs_down, and TKIs_up. The AUROC scores were found to be 0.795, 0.875, 0.745, 0.880, 0.728, 0.867, 0.749, 0.854, 0.773 and 0.854, respectively (Figure 7).

We evaluated the contribution of each feature to the overall classification task in the final model established based on the Mean Decrease Gini assessment inherent in the random forest model. As shown in the Figure 8, we present a bar chart of the Mean Decrease Gini rankings for the top 10 most important features in each model. It can be observed that the feature indicating whether a nucleotide has a double-ring structure consistently ranks highly in terms of feature importance across multiple models. This suggests that this structural feature plays a critical role in model predictions, indicating that purine nucleotides may frequently participate in m6A modification in RNA, potentially contributing to heart-related diseases.

Based on the optimized model described above, we estimated the crosstalk between m6A modification patterns in different cardiac pathological states (Figures 9A,B). It can be observed that, irrespective of the status of m6A modification (upregulated or downregulated), the AUROC scores obtained using disparate models on the same dataset are generally similar, indicating the presence of shared m6A modification sites across distinct cardiac pathological states. Consequently, we hypothesise that m6A-modified sequences exhibit analogous characteristics in heart-related diseases, and that the trained models are capable of more effectively capturing the effects of m6A modification on these diseases.

To further evaluate the reliability and generalizability of m6AHD, we assessed its performance on independent datasets. Specifically, we selected two datasets from the GEO database encompassing three cardiac pathological conditions—dilated cardiomyopathy, ischemic cardiomyopathy, and aortic dissection—and constructed the corresponding feature profiles. The finalized model was then applied to these independent cohorts for prediction. The results demonstrated that our model could effectively identify other cardiac diseases with favorable generalizability, further confirming the similarity of m6A modification patterns across cardiac disorders (Figures 10A,B). These findings highlight the robustness of our approach and underscore its significance for advancing research on the relationship between m6A modifications and cardiovascular diseases.

To verify the association between cardiac pathology and the m6A modification machinery, as inferred by our computational framework, we established a model of evodiamine-induced cardiotoxicity in zebrafish. As shown in Figures 11A–F, exposure to 50 ng/mL Evodiamine resulted in severe cardiac morphological defects across 3, 4, and 5 dpf. The treated group exhibited significantly enlarged pericardial cavity areas and cardiac areas compared to the control group (P < 0.05), indicating successful induction of heart failure-like phenotypes and edema.

To determine if this cardiotoxicity involves the dysregulation of m6A regulatory enzymes, we examined the expression of key m6A methyltransferases. RT-qPCR analysis (Figure 11G) revealed a widespread and significant downregulation of core component genes of the m6A methyltransferase complex and regulatory subunits in the 50 ng/mL group (P < 0.01).

It is worth noting that in our preliminary investigation, we exposed zebrafish to a lower concentration of Evodiamine. As demonstrated in Supplementary Figure S1, the low-dose group did not exhibit significant cardiac malformations or edema compared to controls. Correspondingly, RT-qPCR analysis demonstrated that the expression levels of m6A writers in this phenotypically normal group were preserved, showing no significant difference from the control group.

The collective in vivo results are consistent with the premise of the predictive model. The observation that significant suppression of m6A writers occurred specifically in the high-dose group with severe toxicity–but not in the low-dose group with preserved morphology–strongly suggests that the extent of aberrant m6A modification regulation is closely correlated with the severity of cardiac pathology. The data suggests that Evodiamine-induced cardiotoxicity may be driven by a threshold-dependent collapse of the m6A methylation machinery, leading to a global hypomethylation state that is essential for the disease phenotype.