We retrieved the AML datasets [GSE9476 (Stirewalt et al., 2008) and GSE114868 (Huang et al., 2019)] from the Gene Expression Omnibus (GEO) database. Samples were grouped based on the information provided by authors. GSE9476 comprises 20 normal controls and 26 AML patient samples. GSE114868 comprises 20 normal controls and 194 AML patient samples. Both cohorts have been normalized.

R version 4.3.3 was used to complete this research. Differentially expressed genes (DEGs) of AML were obtained using the package “Limma” (McCarthy and Smyth, 2009). The package “ClusterProfiler” was applied to perform enrichment analysis (Yu et al., 2012). Finally, the package “ggplot2” was used to generate images for visualization.

The package “ESTIMATE” was used to assess stromal scores and immune scores (Yoshihara et al., 2013). “CIBERSORT” was applied to examine the extent of immune cell infiltration within each sample (Newman et al., 2015), based on the known leukocyte expression matrix LM22, and the permutations (PERM) were set to 1,000 to obtain reliable results.

Correlation analysis was performed using the “Corrplot” package.

Weighted Gene Co-expression Network Analysis (WGCNA) was applied to explore the crucial module of AML (Langfelder and Horvath, 2008). We used the genes with coefficients of variation in the top 25% as input data. A soft-threshold value of 14 was selected to ensure the construction of a stable and reliable co-expression network. The dynamic tree-cutting algorithm was used to aggregate genes with similar biological characteristics into the same module. Based on the association index between the module and AML, the key modules of AML were identified.

Seven overlapping genes were obtained by intersecting the three candidate gene sets. The “pROC” package was applied to construct the ROC curve (Robin et al., 2011). The candidate genes were ranked by their AUC values, with the top three designated as hub genes.

Using the threshold determined by the minimum p-value from the log-rank test, AML patients were split into groups with high and low expression levels. The Kaplan–Meier survival curve was then generated to evaluate survival differences between these two groups (Gyorffy, 2024). The GSE76008 cohort was used to clarify the differences in hub genes between leukemia stem cell (LSC)-positive and LSC-negative cells (Ng et al., 2016). The GSE83533 cohort was used to investigate the differences between diagnostic and relapsed AML samples (Li et al., 2016).

The acquisition and utilization of the samples were performed in accordance with the principles of the Declaration of Helsinki. Bone marrow was collected from AML patients and healthy donors of allogeneic hematopoietic stem cell transplantation patients at the Central People’s Hospital of Zhanjiang, and we used human lymphocyte separation medium to isolate the bone marrow mononuclear cells (BMNCs) according to the manual. BMNCs were placed in TRIzol for preservation.

Real‐time quantitative polymerase chain reaction (qRT-PCR) was applied to evaluate the difference in hub genes between healthy controls and AML patients. Total RNA was isolated using TRIzol reagent, and 1 μg of RNA was reverse-transcribed into stable complementary DNA (CDNA) using reverse transcriptase. A volume of 20 µL of qPCR reaction mixture was prepared, containing SYBR, primers, nuclease-free water, and template cDNA. After obtaining the CT values, relative gene expression was assessed using the comparative CT method [2^(-ΔΔCT)]. The primer sequences used were as follows:

Statistical analyses in R were performed using the “ggpubr” package, and images were generated using “ggplot2.” SPSS 19.0 was used for statistical analysis. An independent-samples t-test or a Mann–Whitney U test was selected based on whether the sample conformed to a normal distribution and whether variances between groups were equal. All data were presented as mean ± standard error of the mean (SEM). The images were edited using Adobe Photoshop (PS) 2022 software.

Figure 1 illustrates the overall workflow of this study. GSE9476 and GSE114868 datasets were used to investigate biomarkers of AML. Through differential analysis, WGCNA, and ROC curve analysis, we identified three immune-related hub genes. Furthermore, through prognostic analysis, immune-infiltration profiling, and validation in clinical samples, we clarified the critical roles of these hub genes in the initiation and progression of AML.

The GSE9476 dataset, comprising 20 normal controls and 26 AML samples, was utilized to explore DEGs between the two groups. In total, 1,389 DEGs were obtained, characterized by |log2FC|>1 and adj. p-value<0.05, comprising 450 upregulated genes and 939 downregulated genes (Figure 2A). A heatmap was used to demonstrate the top 10 upregulated genes, namely, CD34, SPINK2, SMYD3, DEPTOR, ATF3, HOXA5, HOXA10, ATP8B4, CLEC11A, and FLT3, whereas the expressions of IL18RAP, CYP4F3, FPR2, CD14, PLBD1, C5AR1, TGFB1, LEF1, IL7R, and HBB were downregulated in AML (Figure 2B). Gene Ontology (GO) Enrichment Analysis clarified that these DEGs were prominently involved in immune response pathways (Figure 2C). KEGG analysis confirmed that hematopoietic cell lineage and Th1 and Th2 cell differentiation were highly associated with AML (Figure 2D). Moreover, Gene Set Enrichment Analysis (GSEA) also suggested that the immune processes may be markedly involved in AML (Figure 2E). Genes with the criteria of |log2FC|>2 and adj. p-value <0.05 in GSE9476 were considered candidate gene set A. Furthermore, to identify feature genes pivotal for the initiation and progression of AML, the GSE114868 dataset was additionally incorporated into the analysis. Differential expression analysis of this dataset yielded a total of 2,850 DEGs (|log2FC|>1 and adj. p-value<0.05), comprising 1,304 upregulated genes and 1,546 downregulated genes (Figure 2F). Finally, genes with the criteria of |log2FC|>3.5 and adj. p-value <0.05 of GSE114868 were considered candidate gene set B.

Compared with the normal controls, AML was characterized by the enrichment of immune-related pathways. Therefore, the CIBERSORT algorithm was applied to explore the immune landscape in the training dataset GSE9476. A stacked bar plot and a heatmap displayed the proportions of immune cells and correlations among 22 immune cell types (Figures 3A,B). We next analyzed the immune cell infiltration in each sample. It was suggested that AML exhibited increased infiltration by immunosuppressive cells (Tregs, M2-like macrophages, plasma cells, and resting mast cells), while the fraction of pro-inflammatory cells (M1-like macrophages and activated CD4⁺ memory T cells) and naïve T/B cells were significantly reduced (Figure 3C), which were characterized by the immunosuppressive microenvironment. These results suggested an imbalance between immunosuppressive and pro-inflammatory cells in AML, which further shaped the dysfunctional immune microenvironment.

WGCNA was applied to identify the key module of AML in GSE9476. A robust co-expression network was established using the power of 14 (Figure 4A). Based on hierarchical clustering and the principle of dynamic tree cutting, a clustering dendrogram was constructed (Figure 4B). Genes with resembling expression models were clustered into a gene module, resulting in a total of nine gene modules, and the red module was the most significant (R ² = 0.94, P < 0.001) related to AML (Figures 4C,D). We performed enrichment analyses on the genes in the red module to explore their potential biological functions. These genes were also notably abundant in immune pathways. The GO enrichment results mainly involved lymphocyte differentiation and immune receptor activity (Figure 4E). KEGG pathway analysis indicated that the DEGs potentially participated in hematopoietic cell lineage and T-cell receptor signaling (Figure 4F). Finally, genes with |GS|>0.75 and |MM|>0.8 in the red module were verified as candidate gene set C.

Seven key genes were obtained by intersecting the three candidate gene sets (Figure 5A), which were strongly correlated (Figure 5B). In the GSE9476 dataset, the expressions of CCR7, SLC16A6, MS4A1, CD79A, IL-7R, and ARG1 were markedly downregulated, while that of FLT3 was upregulated in AML samples (Figure 5C). The consistent findings were subsequently validated in the GSE114868 dataset (Figure 5D). The ROC curve can quantitatively evaluate the disease diagnostic ability of indicators; a higher AUC value indicates stronger diagnostic performance. In some studies, it had been applied to screen for biomarkers (Gong et al., 2025; Li et al., 2025). We applied the ROC curve to assess the diagnostic ability of the seven key genes (Figure 5E). Ranked by the AUC value, the highest AUC was FLT3, followed by CCR7, SLC16A6, MS4A1, IL7R, CD79A, and ARG1. Given that the relationship between FLT3 and AML was well established, we selected the top three genes CCR7, SLC16A6, and MS4A1 as the final hub genes for further study.

AML was characterized by an imbalance in the immune microenvironment and abnormality of the immune process (Figures 2, 3). We further investigated the role of hub genes in AML immunity. The ESTIMATE algorithm was applied to assess the correlation between three genes and the TIME. The results confirmed that CCR7, SLC16A6, and MS4A1 showed a strong positive association with the immune, stromal, and estimate scores, suggesting that these factors could be crucial in remodeling the TIME (Figures 6A–C).

The immune score merely indicated the overall quantity of infiltrating immune cells, but not the actual immune state within the TIME. Therefore, the CIBERSORT algorithm was used to analyze the associations between three genes and 22 types of immune cells in GSE114868. Further analysis revealed strong correlations between key biomarkers (SLC16A6, CCR7, and MS4A1) and specific immune cells. They were positively correlated with some of the pro-inflammatory cell types, such as activated dendritic cells (DCs), monocytes, and activated mast cells (Figure 6D). On the contrary, hub genes were negatively correlated with anti-inflammatory cell types like M2-like macrophages, resting mast cells, and plasma cells (Figure 6D), while these cells were significantly increased in AML (Figure 3). These findings indicated that hub genes might exert an anti-tumor effect by shaping the pro-inflammatory phenotype. In AML, imbalances in the intricate interplay between pro- and anti-inflammatory cytokines can create a tumor-promoting microenvironment that impacts the proliferation and survival of leukemia cells (Binder et al., 2018). A variety of immune-regulatory factors were selected to explore their associations with hub genes, which have been reported to have a definite relationship with AML (Luciano et al., 2022). As a result, hub genes were positively associated with anti-oncogenic cytokines (TNFSF10, TGF-β, IL4, IL1RN, IL10, and IFN-γ) in AML (Figure 6E). SPP1 and KITLG, as oncogenic factors in AML, were negatively correlated with hub genes (Figure 6E). The findings demonstrated that three genes may influence TIME by modulating immune cell infiltration and contributing to the regulation of cytokines.

As previously described, hub genes may affect the development of AML by reshaping the TIME. We next explored the association among three genes and prognosis from three aspects: survival curves, relevance to LSC, and recurrence. Patients with high expression of three genes consistently displayed better prognosis, whereas those with low expression typically faced shorter survival times (Figure 7A). Although most patients can achieve remission through initial treatment, most relapses lead to a poor overall survival. Therefore, recurrence is a key factor affecting the prognosis of AML. We compared the expression of three genes at the time of diagnosis and relapse from paired samples. The results confirmed that three genes were further decreased at relapse, and the trends of CCR7 and SLC16A6 were statistically significant (Figure 7B). Recurrence was usually driven by a rare subpopulation of LSC. Moreover, the result verified that three genes were significantly reduced in LSC⁺ cells (Figure 7C). These results suggested that hub genes may intervene in survival outcomes by regulating the activity of LSC subpopulations and affecting patient recurrence. Furthermore, we used the GSE37642 dataset to examine the associations between the three hub genes and the molecular characteristics of AML patients. As generally recognized, AML patients with RUNX1–RUNX1T1 fusion have a favorable prognosis (Sun et al., 2024); the expressions of CCR7 and MS4A1 were upregulated in this subgroup. Conversely, RUNX1 mutation denotes poor prognosis (Tang et al., 2009), and the expressions of both CCR7 and MS4A1 were downregulated in patients with RUNX1 mutations. Differential expression of the hub genes was also observed across distinct FAB subtypes of AML (Supplementary Figure 1).

It is considered that the BM cell niche promotes leukemogenesis. We collected BMNCs from 5 healthy controls and 13 AML patients to detect the expression of three genes between two groups. The clinical and molecular characteristics of the 13 AML patients can be found in Supplementary Table 1. The results confirmed that, relative to healthy controls, these genes were markedly reduced in AML patients (Figures 8A–C). These results further indicated that three genes, as tumor-suppressor genes, were implicated in the onset and progression of AML.