Construction and Validation of an N7-Methylguanosine–Related Prognostic Model for Acute Myeloid Leukemia

Introduction

Acute myeloid leukemia (AML) is a common hematological malignancy originating from hematopoietic stem and progenitor cells and is characterized by marked biological heterogeneity. It arises from the uncontrolled proliferation of clonal hematopoietic cells and is the most prevalent form of acute leukemia in adults, accounting for approximately 85% of all adult acute leukemia cases.¹ The National Cancer Institute’s SEER Program anticipated an estimated 22,010 new cases of AML in the United States in 2025, representing approximately 1.1% of all newly diagnosed cancers (https://seer.cancer.gov/statfacts/html/amyl.html). AML incidence increases with age, with a median age at diagnosis of approximately 68 years. Clinically, Moreover, AML commonly presents with anemia, bleeding, fever, infection, and extramedullary infiltration. Its pathogenesis is associated with various factors, including myeloproliferative neoplasms or myelodysplastic syndrome, environmental exposures, prior chemotherapy or radiotherapy, and inherited or acquired genetic predispositions, all of which contribute to AML development.² Beyond patient-related factors such as age and comorbidities, AML prognosis is predominantly based on the underlying biological characteristics of the disease.³ However, the intricate pathogenesis of AML remains unclear.

At present, prognosis assessment and treatment decision-making in AML rely primarily on genetics-based risk stratification. The European LeukemiaNet (ELN) (https://www.leukemia-net.org/European LeukemiaNet) updated its risk classification in 2022 for patients eligible for intensive chemotherapy, refining prognostic categorization by integrating the FLT3-ITD allelic ratio with NPM1 mutation status. More recently, a 2024 guideline introduced a distinct risk stratification framework for older or unfit patients receiving low-intensity therapy. Although these classification systems provide a crucial foundation for therapeutic decision-making, substantial prognostic heterogeneity persists within individual risk groups. This limitation highlights the need to identifying additional molecular determinants beyond conventional genetic markers to enable a more precise risk assessment. Therefore, elucidating the precise molecular mechanisms underlying AML development remains an urgent priority. Research efforts to improve AML prognosis must identify potential prognostic features and explore their underlying molecular mechanisms. Effective evaluation and guidance of treatment for patients with AML require the integrated consideration of clinical parameters and genetic prognostic markers.

Recent studies have emphasized the important roles of transcriptional dysregulation and epigenetic mechanisms in AML pathogenesis.^4,5 Among these mechanisms, RNA epigenetic modifications regulate gene expression at the post-transcriptional level and significantly influence RNA properties.⁶ To date, over 170 distinct RNA modifications have been identified, including N7-methylguanosine (m7G), one of the most prevalent post-transcriptional modifications. m7G is commonly present in the 5′ cap structure of mRNA as well as in tRNA and rRNA in eukaryotic cells.⁷ Accumulating evidence suggests that epigenetic modifications play a significant role in cancer initiation and progression. Although the role of m7G in cancer appears to be context dependent, most studies suggest that m7G modification promotes proto-oncogene expression, thus facilitating tumor development and progression.⁸ Additionally, m7G-related lncRNAs have been reported to hold potential value for cancer prognosis prediction and therapeutic guidance.⁶ In AML, emerging evidence demonstrates that epigenetic dysregulation contributes to disease initiation, progression, and therapeutic resistance. The m7G methyltransferase complex, composed primarily of methyltransferase-like 1 (METTL1) and WD repeat domain 4, catalyzes m7G modification in mammalian cells and has been shown to be significantly upregulated in patients with AML. Elevated expression of these enzymes is associated with poor prognosis, suggesting a close link between m7G modification and AML outcomes.⁶ Furthermore, recent studies indicate that m7G-related methyltransferases may contribute to drug resistance in AML cells.⁹ Despite these advances, the specific roles of m7G-related genes (m7G-RGs) in AML pathogenesis and prognosis remain incompletely defined.

In this study, we used bioinformatics methods to investigate the prognostic significance of m7G-RGs in AML and to explore their potential molecular mechanisms. Although interest in m7G modification in cancer has increased, a systematic evaluation of its prognostic relevance and associated immune landscape in AML remains limited. To address this knowledge gap, we integrated transcriptomic data from publicly available cohorts and applied a series of computational analyses—including differential expression analysis, weighted gene co-expression network analysis (WGCNA), and machine learning algorithms—to construct a robust m7G-based prognostic signature. Our analysis focused on elucidating the connections between the identified m7G-RGs and key disease mechanisms in AML, particularly immune cell infiltration and chemotherapy resistance–related pathways. Collectively, the aim of this study was to provide novel biomarkers and mechanistic insights that may facilitate improved risk stratification and more individualized therapeutic strategies for AML.

Materials and Methods

Data Extraction

Gene expression profiles and corresponding clinical data for AML were obtained from The Cancer Genome Atlas (TCGA) database through the University of California Santa Cruz Xena (https://xena.ucsc.edu/). The TCGA-LAML cohort served as training set 1 and included 151 AML samples, of which 132 samples with complete survival information were retained for downstream analyses. Two additional datasets—GSE114868 and GSE37642—were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). GSE114868, generated using the GPL17586 platform, was used as training set 2 and comprised 194 AML and 20 healthy control samples. The GSE37642 dataset served as the validation cohort and contained 984 samples; among these, 136 AML samples generated on the GPL570 platform were included in the analysis. In both GEO datasets, the samples comprised bone marrow mononuclear cells. Additionally, 29 m7G-RGs were curated from previously published literature (Supplementary Table 1).¹⁰

Differential Expression Analysis

Differential expression analysis was conducted on the GSE114868 dataset utilizing the limma package (version 3.52.4)¹¹ to identify differentially expressed genes (DEGs) between the control and AML samples. Genes with an absolute log₂ fold change (|log₂FC|) ≥ 1 and an adjusted p-value < 0.05 were considered statistically significant. Volcano plots and heat maps were generated using the ggplot2 package (version 3.3.6)¹¹ and the Complex Heatmap package (version 2.14.0),¹² respectively, to visualize DEG expression patterns. The heat map specifically displays the top 20 upregulated and top 20 downregulated genes.

Weighted Gene Co-Expression Network Analysis

WGCNA was conducted using the WGCNA package (version 1.71)¹³ to identify gene modules associated with m7G scores in the GSE114868 dataset. Using the curated m7G-RGs as the reference gene set, m7G scores for each sample were calculated using single-sample gene set enrichment analysis (ssGSEA) implemented in the GSVA package (version 1.44.5).¹⁴ Subsequently, differences in m7G scores between the AML and control groups were assessed using the Wilcoxon rank-sum test (p < 0.05). Prior to network construction, hierarchical clustering was performed to identify and exclude potential outliers or aberrant samples in the GSE114868 dataset. The optimal soft-thresholding power (β) was determined according to the scale-free topological criterion to ensure that gene interactions conformed to a near–scale-free distribution. A gene co-expression network was then constructed. Using the hybrid dynamic tree cut algorithm, a minimum module size of 100 genes was specified, and modules with high similarity were merged using a cut height of 0.3 to generate distinct gene modules. Subsequently, correlations between gene modules and clinical traits were calculated using the m7G score as the phenotypic trait (p < 0.05). Modules exhibiting the highest absolute correlation with the m7G score were defined as key modules, and the genes within these modules were designated as key module genes.

Identification and Enrichment Analysis of Intersection Genes

Intersection genes were identified by overlapping DEGs with the key module genes. To further investigate the signaling pathways and biological functions associated with these intersection genes, a functional enrichment analysis was conducted using the clusterProfiler package (version 4.4.4).¹⁵ This analysis encompassed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment (p < 0.05).

Identification of Genes Associated with AML Prognosis

In the TCGA-LAML cohort, the intersection genes were first evaluated using univariate Cox proportional hazards (PH) regression analysis implemented with the survival package (version 3.4–0)¹⁶ to identify genes significantly associated with AML (p < 0.05). Genes identified in the univariate analysis were subsequently assessed for compliance with the PH assumption, and only those satisfying the PH assumption (p > 0.05) were retained for further analysis. The retained genes were then subjected to least absolute shrinkage and selection operator (LASSO) regression using the glmnet package (version 4.1–4).¹⁷ Genes corresponding to the optimal penalty parameter (lambda.min) were selected and further included in a stepwise multivariate Cox regression analysis to recognize independent prognostic genes associated with AML survival.

Construction and Evaluation of the Prognostic Risk Model

A prognostic risk model was constructed based on the independent prognostic genes derived from multivariate Cox regression analysis using the TCGA-LAML dataset. The riskScore for each patient was calculated using the following formula:

where Expi denotes the expression level of gene (i), and Coefi indicates the corresponding regression coefficient derived from the multivariate Cox model. Based on the median riskScore, the 132 patients with available survival data in the TCGA-LAML cohort were split into high- and low-risk groups. Risk distribution curves were generated to evaluate survival outcomes across risk subgroups. Additionally, heat maps were generated to depict the expression patterns of prognostic genes between the two subgroups. Kaplan–Meier (KM) survival analysis was performed using the survminer package (version 0.4.9)¹⁶ to compare overall survival between the two risk subgroups. Meanwhile, a time-dependent receiver operating characteristic (ROC) was generated using the survivalROC package (version 1.0.3.1)¹⁸ to create receiver curves for evaluating the prognostic predictive capability of the risk model in patients with AML. Finally, the predictive performance of the risk model was independently validated in the GSE37642 dataset using the same analytical procedures.

Furthermore, to evaluate the clinical utility of the prognostic risk score within the ELN 2022 risk classification system, the corresponding risk stratification data (favorable, intermediate, and adverse) were obtained from the TCGA-LAML cohort. Differences in risk scores among the three ELN 2022 risk subgroups were evaluated using the Kruskal–Wallis test, with P-values adjusted by the Benjamini–Hochberg method (p < 0.05). Subsequently, KM survival curves were plotted to visualize the prognostic value of the risk score within each ELN risk subgroup.

Clinical Stratified Analysis

To estimate the correlation between the riskScore and clinical features, we compared the differences in riskScore across multiple clinical variables within the TCGA-LAML cohort, including age, sex, FAB classification, cytogenetic risk category, bone marrow blast percentage, and race (p < 0.05). The Wilcoxon rank-sum test was applied for comparisons between two groups, while the Kruskal–Wallis test was used for comparisons among multiple groups.

Independent Prognostic Analysis and Nomogram Construction

To identify independent prognostic markers, riskScore and other clinical parameters were sequentially subjected to univariate Cox regression analysis (p < 0.05), PH hypothesis assumption testing (p > 0.05), and multivariate Cox regression analysis in the TCGA-LAML dataset. Based on the independent prognostic indicators identified, we constructed a nomogram using the rms package (version 6.3–0)¹⁹ to forecast the 1-, 2-, and 3-year overall survival rates in patients with AML. The predictive performance and clinical utility of the nomogram were evaluated using decision curve analysis (DCA) and calibration curves.

Enrichment Analysis Between High- and Low-Risk Subgroups

To characterize differences in biological functions and signaling pathways between the two risk subgroups within the TCGA-LAML dataset, differential expression analysis was conducted. Subsequently, the DEGs were ranked according to log₂FC, and enrichment analysis was conducted using GSEA. The HALLMARK and KEGG gene sets obtained from the msigdbr package (version 7.5.1)²⁰ were used as reference gene sets, with an adjusted p < 0.05 considered statistically significant. Meanwhile, GSVA scores were calculated for each pathway in the TCGA-LAML dataset using the GSVA package (version 1.44.5). The C2: KEGG gene sets from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) served as the reference gene set, enabling comparison of pathway activity between the two risk subgroups (p < 0.05). In addition, ssGSEA was performed using the HALLMARK gene sets from MSigDB to quantify enrichment scores for 50 cancer-related pathways in each AML sample in the TCGA-LAML dataset. Differences in pathway enrichment scores between the two subgroups were analyzed using the Wilcoxon rank-sum test (p < 0.05).

Analysis of the Immune Microenvironment Between Risk Subgroups

In the TCGA-LAML cohort, we applied the ESTIMATE algorithm to calculate StromalScore, ImmuneScore, and ESTIMATEScore for each AML sample. Differences in these scores between the two risk subgroups were evaluated using the Wilcoxon rank-sum test (p < 0.05). Subsequently, immune cell infiltration was assessed using the ssGSEA algorithm in the TCGA-LAML dataset. An immune-related gene signature derived from the published literature served as the reference gene set to calculate enrichment scores for 28 immune cell types.²¹ Differences in immune cell infiltration between the two risk groups were analyzed using the Wilcoxon rank-sum test (p < 0.05). Furthermore, Spearman correlation analysis was conducted to assess the correlation among prognostic genes, the riskScore, and differentially infiltrating immune cell populations.

Drug Sensitivity Analysis and Regulatory Network Construction

Drug sensitivity analysis was conducted using the oncoPredict package (version 0.2)²² to predict the half-maximal inhibitory concentration (IC₅₀) values of antineoplastic agents in AML samples. Version 2 data from the Genomics of Drug Sensitivity in Cancer database (https://www.cancerrxgene.org), comprising 805 cancer cell lines, 17,419 genes, and 198 drugs, were used as the training dataset. Furthermore, comparisons were made between the IC₅₀ values of the two risk subgroups, and their association with riskScore were assessed. Drugs that met the screening criteria of p < 0.05 and |r| ≥ 0.3 were selected for further analysis. To explore potential regulatory mechanisms, the miRNet (https://www.mirnet.ca/) database was employed to forecast transcription factors (TFs) associated with the identified prognostic genes with ≥ 1 degree. Based on these interactions, an mRNA–TF regulatory network was constructed and visualized using Cytoscape software (version 3.9.1).²³

Expression Analysis of Prognostic Genes

To further investigate the expression patterns of prognostic genes, we evaluated their expression in the GSE114868 dataset. To validate the bioinformatics findings, reverse transcription quantitative polymerase chain reaction (RT-qPCR) was performed. Whole-blood samples were collected from five patients with AML and five healthy controls at the First Hospital of Lanzhou University, following approval by the institutional ethics committee. All participants signed an informed consent form. To ensure robustness and objectivity, each biological sample was analyzed in triplicate, sample processing was randomized, and the investigator performing qPCR data analysis was blinded to the group assignments until completion of the initial analyses. Total RNA was extracted from 10 samples using TRIzol reagent (Ambion, Austin, TX, USA) according to the manufacturer’s protocols. Complementary cDNA synthesis was subsequently performed using the SureScript First-strand cDNA Synthesis Kit (Servicebio, Wuhan, China) in accordance with the producer’s instructions. Subsequently, qPCR was performed utilizing the 2× Universal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China) following the manufacturer’s protocol. Primer sequences used for PCR are presented in Table 1. Gene expression levels were normalized to the internal reference gene GAPDH and calculated using the 2^−ΔΔCt method.²⁴

Table 1 The Primer Sequences of PCR

Statistical Analysis

All statistical analyses were conducted using R software (version 4.2.1). A p-value <0.05 was considered statistically significant.

Results

Identification of 2232 DEGs and 276 Key Module Genes in AML

In the GSE114868 dataset, 2232 DEGs were identified between AML and control groups (|log₂ FC| ≥ 1, adjusted p < 0.05). Among these DEGs, 997 genes were upregulated, and 1235 genes were downregulated in AML (Figure 1A and B). Subsequently, WGCNA was implemented to identify module genes associated with the m7G score in the GSE114868 dataset. Wilcoxon rank-sum testing indicated a notable disparity in m7G scores between the AML and control groups. The m7G score was markedly higher in the AML group, indicating the potential involvement of m7G-associated mechanisms in AML pathogenesis and progression (Figure 1C). Hierarchical clustering analysis revealed no outlier samples from the GSE114868 dataset (Figure 1D). The optimal soft-thresholding power (β) was determined to be 5, with the ordinate scale-free topology fit index (R²) approaching the 0.85 threshold, and mean connectivity approximated zero, indicating that the constructed network conformed to a scale-free distribution (Figure 1E). In total, 13 gene modules were identified (Figure 1F). Among them, the MEtan module exhibited the most significant correlation with the m7G score (r = −0.60, p = 1×10⁻²²) and was therefore defined as the key module. This module contained276 key module genes for subsequent analyses (Figure 1G).

Figure 1 The results of differential expression genes analyses in AML patients in GSE114868 dataset. (A and B) Volcano and heat maps of DEGs ((A) Volcano, (B) Heatmap, only log2FC sorted up-regulated. (C) m7G-Score in AML and control differences (****P < 0.0001). (D) Sample Clustering analysis. (E) Soft threshold filtering map. (F) Hierarchical dendrogram of genes and module clustering genes. (G) Relationship between modules and traits.

Enrichment Analyses of 116 Intersection Genes

A total of 116 intersection genes were identified by overlapping the 2232 DEGs with 276 key module genes (Figure 2A). Subsequently, functional enrichment analysis was implemented to explore the signaling pathways and biological processes (BPs) associated with these intersection genes. GO analysis revealed enrichment in 413 terms, including 54 cellular components, 58 molecular functions, and 301 BPs. These GO terms included neutrophil chemotaxis, neutrophil migration, and cytokine receptor activity (Figure 2B). Additionally, KEGG pathway analysis identified 15 significantly enriched pathways, including hematopoietic cell lineage, starch and sucrose metabolism, NOD-like receptor signaling, and IL-17 signaling pathways (Figure 2C). Collectively, these results provide preliminary insights into the biological functions and pathways associated with the intersection genes.

Figure 2 Functional enrichment results. (A) Venn diagram of DEGs and key module genes. (B) GO enrichment analysis of 116 intersection genes (the block size refers to the number of enriched genes; color refers to significance). (C) KEEG enrichment analysis of 116 intersection genes.

Prognostic Performance of the Risk Model in AML

Univariate Cox regression analysis revealed that 17 of the 116 intersection genes were significantly associated with overall survival in patients with AML (p < 0.05) (Figure 3A). All 17 genes satisfied the PH assumptions (p > 0.05). Subsequently, the genes were subjected to LASSO regression analysis, which identified 11 candidate genes: TM6SF1, CYP4F2, IL1R2, MTX1, MANSC1, SUSD3, NCF4, SLC22A4, TUBA4A, VSTM1, and RETN (Figure 3B). Based on these 11 genes, stepwise multivariate Cox regression analysis was performed, resulting in the identification of seven independent prognostic genes: TM6SF1, IL1R2, MTX1, SUSD3, SLC22A4, TUBA4A, and RETN (Figure 3C).

Figure 3 Acquisition of 7 prognostic genes. (A) The forest plot for univariate Cox analysis of 116 intersection genes. (B) LASSO analysis of 17 genes, (Bi) cross-validation curve, (Bii) path plot of regression coefficients. (C) The forest plot for multivariate Cox regression analysis of 11 genes.

Using these prognostic genes, we constructed a risk model and calculated the risk score as follows: riskScore = (TM6SF1 × −0.457) + (IL1R2 × 0.191) + (MTX1 × 0.661) + (SUSD3 × 0.15) + (SLC22A4 × −0.508) + (TUBA4A × 0.427) + (RETN × −0.107). According to the median value of riskScore, we categorized the AML samples into two risk subgroups. Based on the median risk score, patients with AML in both the TCGA-LAML and GSE37642 cohorts were stratified into high- and low-risk subgroups. In both datasets, risk distribution analyses demonstrated that patients in the high-risk group experienced higher mortality rates and shorter overall survival than those in the low-risk group. (Figure 4A and B). The expression patterns of the seven prognostic genes across the two risk subgroups are depicted in Figure 4C and D. KM survival analysis demonstrated notable variances in overall survival between the two risk subgroups in both cohorts. Specifically, patients categorized as high risk demonstrated significantly poorer survival outcomes than those classified as low risk (Figure 4E and F). Additionally, time-dependent ROC analysis showed that the AUC values for 1-, 2-, and 3-year survival exceeded 0.6 in both datasets, indicating satisfactory predictive performance of the risk model (Figure 4G and H). Collectively, these results confirm the robustness and predictive value of the proposed risk model.

Figure 4 The construction and validation of the risk model. (A and B) The risk curve analysis in both TCGA-LAML dataset (Training set) and GSE37642 dataset (Verification set). (C and D) The expression patterns of prognostic genes in both training set and verification set. (E and F) KM survival curve distribution in both training set and verification set. (G and H) ROC curve with AUC in both training set and verification set.

The performance of the risk score was further evaluated within the ELN 2022 risk stratification framework. Kruskal–Wallis testing revealed no significant differences in risk scores among the favorable, intermediate, and adverse ELN subgroups (p = 0.19). However, substantial heterogeneity in risk score distribution was observed within the intermediate-risk subgroup (Figure S1A). KM analysis demonstrated that within this subgroup, patients categorized as high risk had significantly poorer survival outcomes compared to their low-risk counterparts (P = 0.00088). Although survival differences between high- and low-risk groups did not reach statistical significance within the favorable and adverse subgroups, a consistent trend was observed wherein high-risk patients displayed lower survival probabilities (Figure S1B). Together, these findings indicate that the proposed risk model may provide complementary prognostic value by refining risk stratification within the heterogeneous intermediate-risk subgroup of the ELN 2022 classification, potentially enabling more precise risk stratification for these patients.

Construction of a Nomogram Based on riskScore and Age in AML

The association between the riskScore and clinical features was examined. Significant disparities in riskScore were observed across age groups, FAB classifications, and cytogenetic risk groups. Specifically, older patients exhibited higher riskScore values; riskScore differed significantly among FAB subtypes; and patients classified within the intermediate/normal and poor cytogenetic risk groups displayed higher riskScore values compared to those in the favorable group (p < 0.05) (Figure 5A). Subsequently, riskScore and relevant clinical characteristics were sequentially included in univariate Cox regression analysis (Figure 5B), PH assumption testing, and multivariate Cox regression analysis (Figure 5C) to identify independent prognostic factors. RiskScore (HR = 2.63 [1.93–3.60], p < 0.0001) and age (HR = 2.32 [1.46–3.68], p < 0.0001) were identified as independent predictors of overall survival. Based on incorporating riskScore and age variables, a prognostic nomogram was developed to estimate 1-, 2-, and 3-year survival probabilities for patients with AML (Figure 5D). The calibration curves demonstrated good agreement between predicted and observed outcomes, with slopes approaching 1 (Figure 5E), and decision curve analysis indicated that the nomogram provided greater net clinical benefit than either factor alone (Figure 5F). These findings indicate that the nomogram has acceptable predictive accuracy for AML prognosis.

Figure 5 Clinical features and nomogram. (A) The analysis of riskScore and different clinical features. (B) Univariate Cox regression analysis analysis. (C) PH hypothesis test and multivariate Cox regression analysis. (D) nomogram model. (E) The calibration curves of Nomogram predicting the OS of AML patients in 1-, 2-, 3- years (X-axis: predicted probability; Y-axis: actual probability). (F) DCA curves of nomogram in 1- (Fi), 2- (Fii), 3- (Fiii) years.

Signal Pathways Enriched by the Two Risk Subgroups Showed Significant Differences

GSEA revealed that several HALLMARK pathways were significantly enriched in the high-risk subgroup, encompassing heme metabolism, allograft rejection, interferon gamma response, and TNF-a signaling through NF-kb (Figure 6A). Additionally, KEGG pathway analysis showed significant enrichment of immune- and hematopoiesis-related pathways in the high-risk group, such as the chemokine signaling pathway, cytokine–cytokine receptor interaction, and hematopoietic cell lineage (Figure 6B). Consistent with these findings, GSVA identified 43 differentially enriched KEGG pathways between the two risk subgroups, including 37 upregulated and six downregulated pathways. Several pathways overlapped with those identified by GSEA, such as the chemokine signaling pathway and cytokine–cytokine receptor interaction (Figure 6C), further supporting the robustness of the enrichment results. Additionally, the enrichment scores for 19 cancer-related pathways differed significantly between the two risk subgroups, with consistently higher scores observed in the high-risk group compared to the low-risk group for pathways. Notably, enriched pathways included IL2 STAT5 signaling, IL6 JAK STAT3 signaling, and the p53 signaling pathway (Figure 6D). This collective enrichment pattern suggests that the high-risk AML subgroup is characterized by persistent immune-inflammatory activation and dysregulation of cancer-associated signaling pathways.

Figure 6 Investigating the pathological mechanisms of AML. (A) GSEA Enrichment Trend Chart of HALLMARK signaling pathways. (B) GSEA Enrichment Trend Chart of KEGG pathways. (C) Heatmap of GSVA. (D) Inter group differences in cancer-related pathways (*P < 0.05, **P < 0.01, ***P < 0.001).

Association Between riskScore and the Immune Microenvironment in AML

Significant differences were observed in ImmuneScore and ESTIMATEScore between the two risk subgroups, with higher scores in the high-risk group than in the low-risk group. These findings indicate that samples classified as high risk were associated with a greater overall immune cell component (Figure 7A). Immune cell infiltration analysis further revealed notable variances in the enrichment scores of 13 immune cell types, including activated dendritic cells, monocytes, and natural killer (NK) cells, between the two risk subgroups. Specifically, the high-risk subgroup exhibited higher immune cell enrichment scores than the low-risk subgroup (Figure 7B), consistent with the ESTIMATE-based results. Correlation analyses confirmed varying degrees of positive correlation between IL1R2 expression, TUBA4A expression, the riskScore, and differentially infiltrating immune cells. These findings suggest that IL1R2, TUBA4A, and riskScore may represent the level of immune cell infiltration within AML samples to some extent (Figure 7C). Collectively, these results indicate that the high-risk AML subgroup is characterized by a more active, yet potentially dysfunctional, immune microenvironment.

Figure 7 Correlation of the riskScore with the immune microenvironment in AML. (A) Inter group differences of tumor stromal score and immune score in two risk subgroups. (B) Differences in immune cell infiltration levels in two risk subgroups. (C) The correlations of prognostic genes, risk scores and immune cell. (*P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001).

Identification of Drug Sensitivity and TFs

Based on the drug sensitivity analysis, 25 drugs were obtained (Table 2). Among them, the IC50 values of GNE-317 and gemcitabine were lower and showed significant negative correlations with the riskScore (p < 0.05), suggesting enhanced sensitivity in high-risk patients. The results indicate that GNE-317 and gemcitabine may have potential therapeutic value for patients classified within the high-risk subgroup (Figure 8A). Additionally, 36 TFs associated with the prognostic genes were identified through regulatory prediction analysis. An mRNA–TF regulatory network was subsequently constructed. In this network, USF2 regulated SUSD3, TUBA4A, TM6SF1, and IL1R2 simultaneously (Figure 8B). These findings provide insights into potential regulatory mechanisms and may contribute to the development of personalized therapeutic strategies for AML.

Table 2 25 Drugs Were Obtained Based on Drug Sensitivity Analysis

Figure 8 Drug sensitivity analysis and TF-mRNA. (A) The correlation between IC50 values and risk scores of GNE-317 and Gemcitabine treatments (****P < 0.0001). (B) TFs-mRNA network (Red refers to prognostic genes; Yellow refers to transcription factors).

Prognostic Genes Were Downregulated in AML

Expression of the prognostic genes was performed using the GSE114868 dataset. All seven prognostic genes exhibited significant differences between AML and control groups, with significantly lower expression levels observed in the AML group (Figure 9A). Consistent with these findings, RT-qPCR validation confirmed the downregulation of the seven prognostic genes in the AML group compared with controls, with particularly pronounced reductions observed for TM6SF1, IL1R2, MTX1, and SLC22A4. (Figure 9B). The concordance between public transcriptomic data and clinical sample validation strongly supports the robustness and reliability of the identified dysregulated m7G-RGs in AML.

Figure 9 Expression of prognostic genes. (A) Expression of prognostic genes among AML and control groups. (B) The expression t of the 7 prognostic genes in the AML group and control groups by RT-qPCR. (ns P > 0.05, *P < 0.05, ***P < 0.001, and ****P < 0.0001).

Discussion

AML is an aggressive hematological malignancy characterized by high levels of genomic aberrations and marked molecular heterogeneity. Although a subset of patients achieves remission through intensive chemotherapy or allogeneic hematopoietic stem cell transplantation, treatment resistance and high relapse rate contribute to a low 5-year survival rate and poor prognosis.²⁵ m7G modification is one of the most prevalent RNA modifications and is present in various RNAs. It influences nearly all stages of mRNA metabolism, including splicing, polyadenylation, nuclear export, translation, and degradation.²⁶ With advances in the accurate identification of m7G modifications in RNA, their functions in regulating gene expression and various physiological processes have become increasingly clear. Abnormal m7G modification has been closely related to tumor initiation and progression through modulation of oncogenes and tumor suppressor gene expression.²⁷ At present, research on the m7G gene in AML remains limited. In this study, we established a prognostic risk model for AML based on m7G-RGs. Importantly, when evaluated against the current standard ELN 2022 genetic risk stratification system, our model demonstrated significant prognostic utility within the heterogeneous intermediate-risk category, effectively distinguishing patients with distinct survival outcomes. This finding suggests that the m7G-based signature provides complementary prognostic information beyond conventional genetic markers and may refine risk assessment for the clinical subgroup. Consequently, the proposed model represents a novel framework for exploring RNA modification-related mechanisms in AML and warrants further investigation for its potential role in supporting more personalized management strategies.

Prognostic Model Based on m7G-Related Genes and the Biological Significance of Core Genes

Through Cox regression and LASSO regression analyses, we identified seven prognostic genes—TM6SF1, IL1R2, MTX1, SUSD3, SLC22A4, TUBA4A, and RETN—associated with AML. These genes have been reported to play an important role in the development and progression of various malignancies, including lung adenocarcinoma,²⁸ gastric cancer,²⁹ hepatocellular carcinoma,³⁰ and breast cancer.^31,32 TM6SF1 encodes a lysosomal transmembrane protein critical for intracellular and extracellular protein transport. Previous studies have indicated reduced expression of TM6SF1 in AML,³³ and our research corroborated these findings, further demonstrating its downregulation in AML subgroups. This pattern suggests that TM6SF1 may serve as a potential biomarker for AML diagnosis and prognosis. The SLC22A4 gene is located on chromosome 5q23.3 and encodes a membrane-integrated protein belonging to the carnitine/organic cation transporter family. The activity and expression of drug transporter proteins substantially influence intracellular drug concentrations. This transporter family is involved in the cellular uptake of organic cations and metabolites, implicating SLC22A4 in metabolic regulation and cellular homeostasis. Furthermore, specific promoter genotypes of SLC22A4 have been shown to correlate significantly with the imatinib response in chronic myeloid leukemia (CML). Patients with the rs460089 GC genotype have a higher probability of achieving sustained major molecular response, and may serve as an independent predictor of treatment-free remission.³⁴ SLC22A4 is also expressed on primitive cells in AML and has shown potential prognostic value in patients receiving first-line cytarabine (Ara-C) therapy.³⁵ In addition, low methylation levels of the SLC22A4 gene are associated with improved clinical outcomes in both pediatric and adult patients with AML.³⁵ Metaxin-1 (MTX1) is predominantly localized to the outer mitochondrial membrane and is involved in necrosis induced by tumor necrosis factor (TNF). MTX1 is frequently overexpressed in various tumor tissues, and its expression levels are often positively correlated with tumor progression. In CML, MTX1 is closely associated with RAC2, a protein essential for the proliferation and survival of BCR-ABL1–transformed human hematopoietic stem/progenitor cells. RAC2 plays a critical role in maintaining mitochondrial stability.³⁶ RETN is predominantly expressed in immune cells, including monocytes, macrophages, and neutrophils, in several cellular signaling pathways. Elevated RETN expression has been reported in CML, suggesting its potential utility as a diagnostic biomarker. Additionally, increased RETN expression has been observed in relapsed pediatric acute lymphoblastic leukemia (ALL), suggesting its role in monitoring disease recurrence.³⁷ Consistently, RETN expression levels were significantly higher in relapsed pediatric ALL than those in newly diagnosed patients, supporting potential clinical value as a biomarker for detecting leukemia recurrence.³⁸ In AML, four types of mRNA, including TUBA4A, have been reported to enhance risk stratification and prognostic assessment, further supporting the relevance of this gene in AML outcome prediction.³⁹ As a hallmark gene of ferroptosis, IL1R2 has been reported to influence both the prognosis and tumor immune microenvironment in AML.²⁵ Additionally, SUSD3 has emerged as a potential risk factor associated with adverse AML outcomes.⁴⁰ This study further highlights that IL1R2 functions as a prognostic marker for m7G-related AML. Collectively, our analysis and the existing literature suggest that these seven identified m7G-RGs may serve as potential prognostic biomarkers. Their contributions to AML pathogenesis merit further investigation to assess their suitability as potential therapeutic targets.

Enrichment Pathways and Their Association with Therapy Resistance

We explored functional pathways related to riskScores using the HALLMARK and KEGG gene sets as reference databases. Several immune- and cancer-related pathways were significantly enriched in the high-risk subgroup, including the chemokine signaling pathway, cytokine–cytokine receptor interaction, IL6–JAK–STAT3 signaling, IL2-STAT5 signaling, and the P53 signaling pathway. Notably, the enrichment scores for these pathways were significantly higher in the high-risk group than in the low-risk group. Pathways such as IL6-JAK-STAT3 signaling and IL2-STAT5 signaling have been well established as key drivers of AML initiation and progression.⁴¹ Chemokines are small cellular signaling proteins that modulate cell migration, immune responses, and inflammatory responses and play essential roles in various pathological and physiological processes of tumors. Their effects are mainly mediated by binding to specific receptors within the tumor microenvironment (TME).⁴² AML cells release numerous CCL and CXCL chemokines and express several chemokine receptors, particularly CXCR4.⁴³ The CXCL12/CXCR4 signaling axis is critical for the homing of leukemia cells to protect bone marrow niches and strongly correlates with poor prognosis in patients with AML.⁴⁴ High CXCR4 expression is connected with elevated relapse rates and diminished overall survival rates in patients with AML.⁴⁵ Additionally, high CXCR4 expression remains an adverse prognostic factor in patients with AML who have a normal karyotype.⁴⁶ Furthermore, upregulation of CXCR4 expression during chemotherapy has been linked to the development of chemoresistance in leukemia.⁴⁷ In addition, interactions between IL-8 (CXCL8) and its receptors CXCR1/CXCR2 can activate downstream signaling pathways in both stromal and malignant cells, thereby contributing to chemotherapy resistance.⁴⁸ The IL-6/JAK2/STAT3 signaling pathway regulates a broad range of physiological and pathological processes—including inflammation, immune regulation, cell proliferation, and differentiation—and plays a crucial role in tumor development, metastasis, and invasion.⁴⁹ This pathway is crucial in the pathogenesis of leukemia. Emerging studies indicate that mesenchymal stem cells in the AML microenvironment contribute to chemotherapy resistance. This occurs through the activation of epithelial–mesenchymal transition programs mediated by IL-6/JAK2/STAT3 pathway, highlighting its potential as a therapeutic target.⁵⁰ Furthermore, transglutaminase 2 can exert its effects on T-LBL cell behavior by regulating the IL-6/JAK2/STAT3 signaling pathway.⁵¹ Collectively, the CXCL12/CXCR4 and IL-6/JAK2/STAT3 signaling pathways are critical in mediating chemotherapy resistance and tumor development in AML. These pathways reveal significant potential value as therapeutic targets and prognostic indicators.

Tumor Microenvironment Characteristics and Immune Escape

Immunotherapy has achieved meaningful therapeutic benefits in several hematological malignancies. As a crucial component of the bone marrow microenvironment, immune cells influence the occurrence and development of AML. Therefore, our analysis using ssGSEA showed that enrichment scores for immune cell–related gene signatures were significantly higher in the high-risk subgroup than in the low-risk subgroup. It should be noted that ssGSEA scores represent the relative activity of predefined gene sets rather than absolute cell proportions. Moreover, these inferred immune activities may be influenced by variations in leukemia blast burden across samples. Consequently, their precise cellular composition and functional state of immune cell populations within the AML microenvironment warrant further experimental validation through experimental and single-cell–based approaches. The immune components identified in this study included activated dendritic cells, central memory CD4 T cells, gamma delta T cells, B cells, macrophages, MDSCs, monocytes, NK cells, NKT cells, neutrophils, plasmacytoid dendritic cells, T follicular helper cells, and type 1 T helper cells. Lu et al⁵² evaluated immune cell infiltration within the bone marrow microenvironment of AML using computational algorithms and categorized the TME into immune-hot and immune-cold subtypes. Immune-hot subtypes are characterized by excessive activation of immune and inflammatory pathways, leading to increased secretion of inflammatory mediators and T-cell dysfunction, which may enhance responsiveness to immunotherapy. However, despite heightened immune activity, this subtype is often associated with poor prognosis, potentially due to the accumulation of inhibitory immune cell populations, such as regulatory T cells, MDSCs, and macrophages. These inhibitory immune cells contribute to an immunosuppressive microenvironment, inhibit normal immunity, and facilitate immune escape. NK cells play a critical role in antitumor immunity by directly recognizing and eliminating malignant cells, and by producing significant amounts of cytokines, including γ-interferon, TNF, granulocyte–macrophage colony-stimulating factor, and various chemokines that recruit other immune cells and promote robust secondary adaptive immune responses involving T and B cells.⁵³ Chen et al⁵⁴ speculated that tumor cells can actively remodel the TME through various mechanisms. In particular, tumor cells can suppress NK cell activity, thereby inducing immune dysfunction, promoting tumor progression, and remodeling the landscape of the TME.⁵³ Our research findings also suggest that immune cell infiltration in the high-risk groups may be related to immune escape mechanisms, contributing to the unfavorable prognosis observed in these patients. Correlation analysis between the expression of seven m7G-RGs and immune cell infiltration revealed significant associations involving IL1R2, TUBA4A, and the riskScore, as well as differential immune cell populations. These findings indicate that the IL1R2 and TUBA4A genes may play crucial roles in regulating immune infiltration in AML. IL1R2, situated on chromosome 2q12, functions as a decoy receptor for IL-1, and serves as an endogenous inhibitor that negatively regulates IL-1 signaling. This regulatory mechanism underscores the significance of IL1R2 in modulating immune responses and shaping the immune microenvironment in AML.⁵⁵

Potential Implications for Risk-Stratified Therapy

Despite advances in treatment, a substantial proportion of patients with AML fail to achieve remission following standard chemotherapy or experience early relapse, ultimately affecting their prognosis and survival. We utilized public databases to investigate distinctions in predicted drug sensitivity between risk subgroups and to examine the correlation between IC50 values and risk scores. Our analysis revealed that the IC50 values of GNE-317 and gemcitabine were lower, and significantly lower IC50 values were observed in the high-risk group compared with the low-risk group. This indicates that patients classified as high risk might exhibit higher sensitivity to these drugs in a computational context. Although these results indicate a hypothetical link between our risk stratification model and drug response, they remain hypothetical and warrant further experimental validation and prospective clinical investigation. Elucidating the mechanisms underlying differential drug sensitivity among risk groups may provide new insights into therapeutic optimization for AML. Gemcitabine and cytarabine are nucleoside analogs with similar structural features, both necessitating intracellular phosphorylation for activation and therapeutic efficacy. However, gemcitabine possesses unique mechanisms of action and self-potentiation that differentiate it from cytarabine.⁵⁶ Notably, gemcitabine demonstrates higher affinity for deoxycytidine kinase, resulting in a more efficient generation of its monophosphate and diphosphate metabolites compared with cytarabine and fludarabine.⁵⁷ When used in combination regimens, gemcitabine retains activity in relapsed and refractory AML, a clinical setting that remains particularly difficult to treat.⁵⁷ In addition, computational analyses have revealed that gemcitabine is effective against KMT2A-rearranged acute leukemia cell lines.⁵⁸ Furthermore, over 50% of patients with AML exhibit overactivation of the PI3K/AKT/mTOR signaling pathway.⁵⁹ The important components of this pathway represent attractive therapeutic targets in AML. Consistent with this notion, in vitro experiments have shown that the dual PI3K/mTOR inhibitor VS-5584 effectively inhibits AML cell proliferation.⁶⁰

Limitations and Future Perspectives

To the best of our knowledge, this study is the first to construct and validate an m7G–related prognostic risk model in AML. However, this study also has its limitations, which should be acknowledged. First, our model was developed and validated using retrospective public datasets; thus, prospective clinical studies are essential to confirm its reliability and clinical applicability in real-world settings. Second, the availability of comprehensive clinical data—including specific treatment regimens (eg, chemotherapy protocols or targeted therapies such as gilteritinib and venetoclax) and complete molecular mutation profiles—was limited in the analyzed cohorts. This restriction precluded more refined analyses that examined associations between the risk model, therapeutic responses, and genetic subtypes. Third, although standardized bioinformatics pipelines were applied, residual batch effects across datasets may have influenced the stability and generalizability of the findings. Finally, additional experimental studies are required to substantiate our findings and to elucidate the molecular mechanisms by which the seven m7G-RGs influence AML prognosis, immunotherapy, and drug resistance. Future research will focus on further investigating the expression patterns and functional roles of these genes in AML.

Conclusion

In this study, seven m7G-RG genes with prognostic value in AML were identified, and a novel m7G-based prognostic risk model was constructed. Our findings demonstrated that patients classified into the high-risk group exhibited significantly poorer survival outcomes. Importantly, when evaluated alongside the ELN 2022 genetic risk stratification system, our model provided significant prognostic refinement within the clinically heterogeneous intermediate-risk category, effectively identifying patients with distinct survival outcomes. Given that this model was derived from retrospective datasets, prospective validation is necessary to confirm its clinical utility. Ultimately, this m7G-based prognostic model may serve as a complementary tool to enhance risk stratification, particularly within the ELN 2022 intermediate-risk subgroup, potentially supporting the development of more individualized management strategies for patients with AML.