Integrated analysis of single-cell RNA-seq and bulk RNA-seq reveals RNA N6-methyladenosine modification associated with prognosis and drug resistance in acute myeloid leukemia

Introduction: Acute myeloid leukemia (AML) is a type of blood cancer that is identified by the unrestricted growth of immature myeloid cells within the bone marrow. Despite therapeutic advances, AML prognosis remains highly variable, and there is a lack of biomarkers for customizing treatment. RNA N6-methyladenosine (m⁶A) modification is a reversible and dynamic process that plays a critical role in cancer progression and drug resistance.

Methods: To investigate the m⁶A modification patterns in AML and their potential clinical significance, we used the AUCell method to describe the m⁶A modification activity of cells in AML patients based on 23 m⁶A modification enzymes and further integrated with bulk RNA-seq data.

Results: We found that m⁶A modification was more effective in leukemic cells than in immune cells and induced significant changes in gene expression in leukemic cells rather than immune cells. Furthermore, network analysis revealed a correlation between transcription factor activation and the m⁶A modification status in leukemia cells, while active m⁶A-modified immune cells exhibited a higher interaction density in their gene regulatory networks. Hierarchical clustering based on m⁶A-related genes identified three distinct AML subtypes. The immune dysregulation subtype, characterized by RUNX1 mutation and KMT2A copy number variation, was associated with a worse prognosis and exhibited a specific gene expression pattern with high expression level of IGF2BP3 and FMR1, and low expression level of ELAVL1 and YTHDF2. Notably, patients with the immune dysregulation subtype were sensitive to immunotherapy and chemotherapy.

Discussion: Collectively, our findings suggest that m⁶A modification could be a potential therapeutic target for AML, and the identified subtypes could guide personalized therapy.

1 Introduction

Acute myeloid leukemia (AML) represents an infrequent hematologic malignancy distinguished by the aberrant expansion of precursor cells in the myeloid lineage, resulting in a perturbed differentiation process. In recent decades, allogeneic hematopoietic stem cell transplantation (allo-HSCT) has emerged as a pivotal therapeutic intervention, substantially enhancing the overall survival of eligible AML patients (1). However, the effectiveness of AML treatment and the issue of relapse continue to pose significant clinical challenges. Until recently, there has been a considerable gap in meeting the treatment needs of AML patients.

RNA methylation, particularly N6-methyladenosine (m⁶A) methylation, constitutes the prevalent intramolecular modification observed within eukaryotic mRNA molecules. Such modification is involved in regulating various aspects of RNA, including splicing, stability, localization, and translation (2, 3). In recent times, m⁶A modifications have not solely surfaced as a novel stratum of epigenetic control within the realm of cancer; they have also exhibited substantial therapeutic promise for addressing diverse forms of malignancies (3–6). Notably, current studies have highlighted the potential of m⁶A modifications in elucidating the pathogenesis and therapeutic aspects of AML. Specifically, the m⁶A methyltransferases METTL3 and METTL14 can control and/or maintain myeloid leukemia cells (7–9). The m⁶A readers IGF2BP2 and IGF2BP3 promote AML development in an m⁶A-dependent behavior by controlling the expression level of critical genes in the glutamine metabolism pathways (10, 11). The RNA-binding protein YBX1 plays a pivotal role in sustaining myeloid leukemia cells through an m⁶A-dependent mechanism, wherein it regulates the stability of BCL (12). FTO plays an oncogenic role in AML as an m⁶A RNA demethylase. Small molecule inhibitors that target FTO have the potential to be used in the treatment of AML (13, 14). These findings provide valuable insights into the crucial and intricate involvement of m⁶A modification and its regulators in AML, highlighting a promising avenue for AML treatment.

In this study, we identified two distinct patterns of m⁶A modification intensity based on the m⁶A activity score. We found that high-intensity modes significantly impact the prognosis of AML. Additionally, we conducted a comprehensive analysis of the role of m⁶A modification in AML subtypes, including their genomic alterations, tumor immune microenvironment, and immunotherapy implications. This examination notably enriches our comprehension of the molecular processes that contribute to m⁶A modification’s role in the development of AML and its impact on the effectiveness of drug treatments. Overall, our findings provide a solid foundation for the development of m⁶A modification-targeting therapies for AML and suggest that this approach could be an effective and specific strategy for cancer treatment.

2 Materials and methods

2.1 Data download and processing

We got patient data for acute myeloid leukemia (AML) from the TCGA (The Cancer Genome Atlas) database, which was downloaded from the following source: https://xena.ucsc.edu/. RNAseq data from a total of 151 samples were re-analyzed. Clinical survival information was obtained for 132 patients and somatic mutation information was obtained for 143 patients. To verify, we conducted searches in the GEO database using the keywords ‘acute myelogenous leukemia’ and ‘LAML’ to identify relevant datasets. Specifically, we included datasets in our analysis if they met the following criteria: (1) Adequate sample size (n > 200). (2) Using RNA-seq instead of microarrays. (3) Count matrix is provided for further analysis. GSE106291 (n = 250) and GSE146173 (n = 246) were selected to further analysis (15, 16). We utilized the dataset GSE178926, which comprises targeted immune gene expression profiles derived from pre-treatment bone marrow samples of acute myeloid leukemia patients undergoing treatment with pembrolizumab and the hypomethylating agent azacytidine (ClinicalTrials.gov Identifier: NCT02845297) as described by Rutella et al. in 2022 (17), for the purpose of immunotherapy prediction. Single-cell RNA sequencing (scRNA-seq) data from 12 patients diagnosed with acute myeloid leukemia (AML) were obtained from the dataset provided by van Galen et al. (18).

2.2 Re-analysis of scRNA-seq

The single-cell RNA sequencing (scRNA-seq) count data underwent normalization using the ‘Seurat’ R package (19), followed by log-transformation with an offset value of 1 and subsequent scaling. We discerned genes displaying significant variability through the utilization of the ‘FindVariableFeatures’ function, wherein the ‘vst.method’ parameter was configured as ‘vst’. We conducted Principal Component Analysis (PCA) on the highly variable genes, utilizing the top 30 principal components for subsequent analyses. Cell type annotation was obtained from Zeng et al. (20). In order to represent the outcomes visually, we employed t-distributed stochastic neighbor embedding (t-SNE) to reduce the complexity of the dataset. The RunTSNE function was used to generate a 2-dimensional t-SNE plot based on the top 30 principal components. The resulting t-SNE plot was then used to visualize the clustering results. Subsequently, we employed the ‘FindAllMarkers’ function to identify cluster-specific markers, with the ‘method’ parameter set to ‘MAST’ (21).

2.3 m⁶A modification activity score

We obtained 23 m⁶A modification enzymes (FTO, CBLL1, FMR1, HNRNPC, HNRNPA2B, IGF2BP1, IGF2BP2, ELAVL1, IGF2BP3, ALKBH5, LRPPRC, METTL3, METTL14, RBM15, RBM15B, VIRMA, WTAP, YTHDF1, YTHDF2, YTHDC1, YTHDC2, YTHDF3, ZC3H13) and scored their m⁶A modification activity using the AUCell R package (22). The enzymes were used as a gene set to calculate the area under the curve (AUC) value. This value was used to rank gene expression for each cell, reflecting the proportion of highly expressed genes in the gene set for that cell. To identify the active gene set, we used the “AUCell_exploreThresholds” function to calculate the threshold value (0.048). Using the AUC scores of each cell, we colored the cell clustering UMAP embedding to show which cells were active.

2.4 Regulon analysis and signature enrichment

To assess the activity of transcription factor (TF) regulons in the context of single-cell RNA sequencing (scRNA-seq) data, we conducted a regulon analysis employing the SCENIC framework. Following the guidelines outlined by Van de Sande et al., we leveraged the Docker image of pySCENIC and utilized logarithmically transformed expression counts obtained from AML cells for input data (23). The identification of putative transcription factors (TFs) was carried out with default parameters and a compiled list of human TFs from Lambert et al. (24). To refine potential TF-target interactions within each regulon, we incorporated CisTarget, which entailed the utilization of established human TF motifs databases annotated at intervals of 500 base pairs, 5 kilobases (kb), and 10 kb from transcriptional start sites. Additionally, a drop-out masking strategy was applied during this refinement process. Subsequently, the enrichment of these refined TF regulons was assessed using AUCell, and enrichment scores were scaled for visualization purposes.

2.5 GRN construction

We generated gene regulatory networks (GRNs) for both m⁶A-active and m⁶A-inactive cell populations using the R package bigSCale2 (25, 26). To do so, we first separated the m⁶A-active and m⁶A-inactive cells and generated a cell count matrix for each m⁶A modification state using Seurat’s GetAssayData function. We then filtered the resulting matrices to remove genes with ensemble identifiers and passed them to bigSCale2 to construct the networks. The networks were generated under the “normal” clustering parameter, with an edge cutoff set to the top 0.9 quantile for correlation coefficient. We visualized the networks using the R package igraph, and each network’s layout was derived from the Fruchterman-Reingold algorithm.

2.6 Construction of cell-specificity m⁶A related signature

The single-cell approach enables a detailed exploration of the m⁶A modification model. We employed a likelihood-ratio test to establish a m⁶A-related gene signature for each cell-type-specific cluster. This entailed identifying differentially expressed genes between cells exhibiting m⁶A activity and those that are m⁶A inactive within each cluster. We focused on identifying m⁶A related genes based on the size of effect value (log2FC). The m⁶A-related genes within the various cell types were identified using the following criteria: |log2FC| > 0.25, percentage of cells (pct) < 0.1, and p-value < 0.05. We prefer to obtain m⁶A-related genes specific to each cell type, we therefore used and evaluated the cell specificity of m⁶A-related genes. Cell-specific genes within the different cell types were identified using the following criteria: log2FC > 0.25, percentage of cells (pct) < 0.1, and p-value < 0.01. Cell-specificity m⁶A related signature emphasizes strict cell specificity and comprehensive m⁶A modification profile, we thus relaxed the screening criteria for m⁶A-related genes and limited the screening criteria for cell markers.

2.7 Reanalysis of bulk RNA-seq

DEseq2 software was used to detect Differentially expressed genes (DEGs). Count matrix was scaled by variance-stabilizing transformation. We identified differentially expressed genes (DEGs) among the three subtypes using a threshold of |log2FC| ≥ 1 and a false discovery rate (FDR) ≤ 0.01. Subsequently, we conducted Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, visualizing the results through Metascape. Furthermore, Molecular Signatures Database V7.4 of hallmark gene sets were used for GSEA by the R package ClusterProfiler.

2.8 Clustering, ESTIMATE, and ssGSEA

For TCGA-LAML, we integrated m⁶A related genes from five types of blood cancer cells and further identified the genes that determine survival status by Cox regression. We conducted hierarchical clustering analysis using the ‘ConsensusClusterPlus’ package (R implementation, K = 3) based on the 40-genes (Cox regression p-value < 0.05) from cell-specificity m⁶A Related Signature. We utilized ESTIMATE to assess the stromal score, immune score, ESTIMATE score, and tumor purity for each LAML sample (27). The ssGSEA was conducted to determine the enrichment levels of immune-related pathways in each sample.

2.9 Comparing tumor immune microenvironment

We conducted a systematic ssGSEA analysis to assess the activation levels of 17 immune pathways, as defined by ImmPort (28). Additionally, we quantified the expression levels of 78 immunomodulators based on a prior study, with the aim of exploring the immune microenvironment in LAML. To discern differences in immune signatures, including immune pathways and immunomodulators, among LAML subtypes, we employed both the Kruskal-Wallis test and the Wilcoxon rank-sum test.

2.10 Survival analysis

We defined overall survival (OS) as the duration between diagnosis and either death or last follow-up. We evaluated the differences in OS within subtypes in each cohort using Mantel-Cox log-rank tests, implemented using the R package ‘survival’. Kaplan-Meier plots were generated using the R package ‘survminer’ to visualize the survival curves for each cluster. We also derived univariate and pairwise hazard ratios (HRs) for each cluster using Cox proportional hazards regression.

2.11 Construction of composite machine learning model

The genes that satisfy the two conditions are used to build the machine learning (ML) model: 1. m⁶A related genes; 2. FDR<0.01 in the difference analysis of the three subtypes. We employed a systematic machine learning (ML)-based framework to construct subtype prediction models. The process involved several key steps: 1. Data Preprocessing: This step encompassed imputing missing values and scaling continuous features. Continuous features were standardized to have a zero mean and unit variance. 2. Data Splitting: Following preprocessing, we divided the data into two sets, a training dataset for model development, and a validation dataset for assessing model accuracy. The data split was random and maintained a ratio of 7:3. 3. Model Development: For each of the five ML algorithms (random forests, support vector machines, naïve Bayes, K-nearest neighbors, and neural networks), we developed optimal models using the training dataset. 4. Validation: To ensure robustness, we required consistent results from at least three models for each sample. All of the aforementioned modeling steps were implemented using R (version: 4.1.2).

2.12 Genomic analysis

We conducted Copy Number Alteration (CNA) analysis using GISTIC2.0 on the TCGA LAML cohort. We examined variations in amplification or deletion events at the gene level across the three subtypes. To visualize the Copy Number Variation (CNV) data, we utilized the ‘ComplexHeatmap’ package in R, creating a waterfall plot. Additionally, we calculated the Tumor Mutation Burden (TMB) by determining the number of mutations per patient.

2.13 Immunotherapy response score generation

Immunotherapy signature feature selection was performed using differentially expressed genes between CR (Complete Remission) and NR (Non-Responders) samples in the GSE178926 dataset. Sample identity information was established based on Rutella et al. study (17). We selected genes based on their association with CR compared to NR samples, utilizing Elastic Net penalized logistic regression (glmft algorithm, glmnet package) with alpha = 0.6 and lambda = 0.06 to accommodate the high degree of correlation observed among some genes. We opted for Leave-One-Out Cross-Validation (LOOCV) to maximize our training set for validation, especially given our small dataset size (n = 33), as LOOCV exhibits low bias.

2.14 Chemotherapeutic response prediction

Chemotherapeutic response predictive model is based on the Cancer Genome Project (CGP) data from the Genomics of Drug Sensitivity in Cancer (GDSC) project, using gene expression and drug sensitivity data from cancer cell lines. The R package “pRRophetic” utilized ridge regression to estimate the half-maximal inhibitory concentration (IC50) for each sample (29). Model training was performed using blood-type cell lines through 10-fold cross-validation based on the GDSC’s training set. A total of 45 cell lines derived from hematological cancers were utilized. To mitigate batch effects, we employed the ‘combat’ algorithm with default settings and considered only tissue types labeled as ‘blood’ Duplicate gene expression data were averaged.

2.15 Statistical analyses

The statistical analyses were conducted utilizing R software, version 4.1.2. The comparison of two groups was established via Wilcoxon-rank sum test, while analysis of more than three groups was accomplished through execution of Kruskal-Wallis test. The benchmark for significance was established with p-value, wherein a level of 0.05 or lower was deemed significant, while a level of 0.01 or lower was deemed extremely significant.

3 Results

3.1 Reanalysis of scRNA-seq data identifies active and inactive m⁶A methylation cells in AML

we conducted a comprehensive reanalysis of scRNA-seq data obtained from 13,653 cells from 12 patients diagnosed with AML (GSE116256) (18, 20). Cell identity information was established based on prior studies encompassing five leukemic cell types (LSPC, ProMono-like, Mono-like, GMP-like and cDC-like) and seven non-leukemic immune cell types (B, T, natural killer, plasma, monocytes, CTL and cDCs) (Supplementary Figure 1A). Our study examined 23 key genes (see method) involved in the m⁶A modification process in AML cells. There is no evidence of cell-specific expression for these genes (Supplementary Figure 1B). We used AUCell to score m6A sets within individual cells to further understand m⁶A activity. The AUC values across all cells showed two peaks, with 8465 cells showing relatively higher AUC values when the AUC threshold was set to 0.048 (Figure 1A). Next, we categorized AML cells into two distinct groups based on their m⁶A modification: m⁶A active (AUC value > 0.048) vs. m⁶A inactive (AUC value ≤ 0.048). We performed differential gene analysis on cells with varying m⁶A activity and observed 59 significantly upregulated genes and 46 significantly downregulated genes (|log2FC| > 0.25 and p-value < 0.01, Figure 1B). Furthermore, we performed pathway enrichment analysis on these genes, revealing their impact on myeloid leukocyte cytokine production, mRNA metabolism, and immune effector functions (Figure 1C). We found that the Hallmark gene set, which is comprised of a broad range of biological processes, provided a more comprehensive representation of differences between m⁶A-active and m⁶A-inactive cells (Figure 1D). Finally, we found that a higher proportion of leukemic cells (χ² = 88.81, p-value = 4.34e−21, Supplementary Figure 1C), such as LSPCs, Mono-like, ProMono-like, GMP-like and cDC-like blasts were identified as m⁶A-active cells than immune cells (χ² = 422.95, p-value = 8.17e−84, Figure 1E). This suggests that m⁶A modification activity may play a critical role in regulating a wide range of biological processes and underscores the importance of considering the overall functional landscape of cellular processes when studying the effects of m⁶A modification.

FIGURE 1

Figure 1 Identification of active and inactive m⁶A methylation cells in AML. (A) Score of 23 m⁶A modification activity. The threshold was chosen as 0.048 and the m⁶A modification score of 8,465 cells exceeded the threshold value (The dotted lines overlaying each distribution represent Gaussian fits to the distribution data.). (B) The volcano plot illustrates the differential expression of genes between m⁶A modification-active cells and m⁶A modification-inactive cells. (C) GO gene set enrichment analysis results were shown as lollipop plot where on x-axis, -log of FDR adjusted p-value for GO terms were shown. (D) Barplot of significantly activated pathways in the m⁶A-active (blue) and m⁶A-inactive (green) cells. (E) Stacked barplots showing the frequencies of m⁶A-active and m⁶A-inactive cells in 14 cell types (p-value > 0.05: ns, p-value < 0.05: *, 0.05 < p-value < 0.01:**, 0.01 < p-value < 0.001:***, 0.001 < p-value < 0.0001:****).

3.2 Effects of m⁶A modification status on gene regulatory networks and biological processes in leukemic cells

We identified leukemic cells, including LSPCs, cDC-like, GMP-like, Mono-like and ProMono-like blasts, which exhibited the most significant impact on m⁶A modification, as evidenced by the observation of more up- and down-regulated genes in their active and inactive states compared to other cell types (Figure 2A). To fully appreciate the complexity of m⁶A modification, it is essential to complement gene expression analysis with an understanding of the underlying gene regulatory network. In this regard, we have constructed regulatory networks for leukemic and immune cells with distinct m⁶A modification patterns using SCENIC, a transcription factor-based gene regulatory network. Compared with immune cells, leukemia cells showed greater differences in transcription factor activity between m⁶A-active and m⁶A-inactive cells (Figure 2B). Some transcription factors have higher activity in leukemic cells with active m⁶A modification, such as BATF, GABPA, E2F8, E2F2, E2F3, E2F7, ELF1, YY1, HOXA9 (30–35) (Supplementary Figure 1D, Supplementary Table 1). The transcription factors ELF1, YY1, and E2F are key regulatory elements with a ubiquitous impact across diverse cell types, exerting a significant influence on the modulation of gene expression changes associated with M⁶A modifications (Figure 2C; Supplementary Figures 1D, E). Tumor necrosis factor α-inducible protein 8 (TNFAIP8) is a novel anti-apoptotic molecule that plays a role in AML chemoresistance (31, 36). ELF1 was reported to be responsible for the upregulation of TNFAIP8 expression in human AML patients (31). YY1 has been reported to bind to the promoter region of METTL3 and promote its expression, resulting in increased AML cell proliferation (34). E2F transcription factor 1 (E2F1) was reported to be involved in AML cell differentiation recently (35).

FIGURE 2

Figure 2 Effect of m⁶A modification status on transcriptional regulatory networks and pathway enrichment in four types of leukemia cells (A) Bubble plot showing the genes that were differentially expressed between m⁶A-active and m⁶A-inactive in each cell type. Node sizes correspond to -log10 of FDR adjusted p-value. (B) The volcano plots of differential transcription factor activity analysis of single-cell RNA sequencing data performed on each cell type comparing m⁶A-active cells vs. m⁶A-inactive cells. (C) Regulatory networks of transcription factors in four major types of leukemia cells (cDC-like, GMP-like, LSPCs and ProMono-like). (D) Pathway enrichment network of four major types of leukemia cells.

To identify the unique m⁶A-related genes in each of these cell types, we developed a three-step workflow (Supplementary Figure 1F). Firstly, we obtained the marker genes for each cell type, which enabled us to accurately classify the cells. Secondly, we performed differential expression analysis between m⁶A active and inactive cells in each cell type. Finally, we obtained a list of unique m⁶A-related genes for each cell type by overlapping the two sets of genes. We performed pathway enrichment analysis of the identified m⁶A-related genes for each cell type, which revealed their biological significance. In cDC-like blasts, m⁶A modification affected pathways related to antigen processing and presentation and 3’-UTR-mediated mRNA stabilization (Figure 2D). In GMP-like blasts, m⁶A modification affected pathways related to negative regulation of mRNA metabolic process and maintenance of DNA methylation (Figure 2D). In LSPCs, m⁶A modification affected pathways related to DNA topological change and mRNA methylation (Figure 2D). In ProMono-like blasts, m⁶A modification impacted pathways related to negative regulation of mRNA metabolic process and mRNA modification (Figure 2D). In Mono-like blasts, m⁶A modification affected pathways related to RNA metabolic process and transport (Supplementary Figure 1G). These findings provide insights into the distinct roles of m⁶A modification in different cell types and highlight the importance of this epigenetic mechanism in regulating cellular function and its relationship to AML progression.

3.3 The impact of m⁶A modification on gene regulatory networks in AML

To better understand how leukemia cell or immune cell global regulatory networks are altered in m⁶A modification pattern, a gene regulatory network (GRN) reasoning methodology called bigSCale(25, 26) was implemented to m⁶A-active or m⁶A-inactive cells (Supplementary Table 2). In leukemic cells, the m⁶A-active and m⁶A-inactive GRN exhibited the similar densities (Figure 3A). To assess the biological relevance of the differences between m⁶A-active and m⁶A-inactive networks in leukemic cells, the major distinguishing variables PageRank delineating their topology were investigated. We distinguished the top 5 influencer genes ranked by PageRank, as a proxy for a gene’s influence on the network. Many of these genes have been previously linked to AML pathology, such as CDK1 and CKS1B (37, 38) (Figure 3A). As an illustration, two genes, CDK1 and CKS1B, known for their involvement in G1/S and G2/M phase transitions of the eukaryotic cell cycle, exhibited no significant differential expression (p-value > 0.05) but concurrently showed increased centrality measures (as depicted in Figure 3A). This finding is particularly intriguing, given that leukemic cells, despite being the most deregulated cell type in the initial AML analysis, did not originally highlight the significance of CDK1 and CKS1B in this context. In each LPSCs, cell cycle distribution was determined by Zeng et al. (20). Cycling LSPCs exhibited greater m⁶A modification activity (Supplementary Figure 1H). However, in immune cells, the topology of the m⁶A-active GRN showed a relatively higher density in active cells compared with inactive cells, as reflected by gene-gene relationships and modularity (Figure 3B). Overall, our findings suggest that m⁶A modification pattern in leukemia cells or immune cells may have different effects on GRN topology, which could provide important insights into the pathological mechanisms of diseases.

FIGURE 3

Figure 3 Gene regulatory network of m⁶A-inactive and m⁶A-active cells in leukemic cells and immune cells. (A) Gene regulatory network of m⁶A-inactive and m⁶A-active cells in leukemic cells (inactive: Cutoff = 0.9, Node = 284, Edge = 883; active: Cutoff = 0.9, Node = 256, Edge = 873). (B) Gene regulatory network of m⁶A-inactive and m⁶A-active cells in immune cells (inactive: Cutoff = 0.9, Node = 146, Edge = 260; active: Cutoff = 0.9, Node = 2036, Edge = 11078).

3.4 Classification of AML subtypes using the m⁶A-related genes in leukemia cells

We integrated the unique m⁶A-related genes of cDC-like blasts, GMP-like blasts, LSPCs, Mono-like blasts and ProMono-like blasts, and performed a consensus clustering analysis on TCGA LAML samples (Supplementary Table 3, see method). Our clustering analysis revealed three distinct subtypes with high consistency and robustness, which were identified based on the cumulative distribution function (CDF) of the consensus matrix (Figure 4A). Subsequently, we conducted survival analysis on these subtypes using Kaplan-Meier curves and found significant differences in OS between the subtypes (Figure 4B). We performed differential gene expression analysis on each subtype and identified pathways that were significantly enriched (Figure 4C). The first subtype, characterized by the poorest survival, demonstrated significant enrichment in immune system processes such as negative regulation of immune system processes, innate immune response, and adaptive immune response, in addition to pathways related to cell activation, cytokine stimulus, and inflammatory response. Based on these findings, we defined this subtype as the immune dysregulation (ID) subtype (Figure 4C). The second subtype, with intermediate survival, exhibited significant enrichment in pathways related to ovulation, response to xenobiotic stimulus, and neutrophil-mediated immunity, along with pathways related to chemotaxis, blood circulation, and hormone processing. We defined this subtype as the hormone regulation (HR) subtype (Figure 4C). The third subtype, which had the best survival, demonstrated significant enrichment in pathways related to cell morphogenesis, extracellular matrix organization, and lipid transport, as well as pathways related to cell junction organization and signaling by VEGF. Based on these findings, we have defined a new subtype, called the cellular adaptation (CA) subtype (Figure 4C; Supplementary Table 4).

FIGURE 4

Figure 4 Classification of AML subtypes using the m⁶A-related genes in leukemia cells by K-means analysis. (A) K = 3 was identified as the optimal value for consensus clustering. (B) Kaplan-Meier curves of overall survival (OS) among the three subtypes in the TCGA LAML cohort. (C) The enrichment statistics of ssGSEA signaling pathways of immune dysregulation subtype, hormone regulation subtype and cellular adaptation subtype. (D) The expression levels of m⁶A-associated genes of immune dysregulation subtype, hormone regulation subtype and cellular adaptation subtype. (E) Sankey plot showing the correlation between the classification of AML subtypes using the m⁶A-related genes and French-American-British classification of acute myeloid leukemia. (F) The survival curve demonstrates the prognostic outcomes of FAB classification from TCGA-LAML. (G) Kaplan-Meier curves of OS among the three subtypes in two other GEO datasets (GSE146173 and GSE106291).

To comprehensively characterize the landscape of m⁶A modifications across subtypes, we have examined the expression levels of m⁶A-associated genes. Our analysis revealed that each subtype exhibits a distinctive m⁶A expression profile, suggestive of subtype-specific regulation of m⁶A modification (Figure 4D). For example, the immune dysregulation subtype was characterized by high expression of IGF2BP3, IGF2BP2, and FMR1, and low expression of ELAVL1, FTO, and METTL14 (Figure 4D). We further investigated the correlation between AML subtype classification using m⁶A-related genes and the French-American-British (FAB) classification AML (39). The classification obtained from m⁶A-related genes, although not specifically correlated with FAB classifications, demonstrates a robust prognostic diagnostic capability (Figures 4E, F). The three subtypes are characterized from different aspects, including m⁶A gene expression, m⁶A modification activity of leukemic cells, and subtype distribution. To examine the robustness of our subtype identification, we have developed a composite machine-learning classifier based on m⁶A-related genes and the DEGs among the three subtypes (see method, Supplementary Table 5). We further validated the stability and reproducibility of the three subtypes in two additional GEO datasets, providing strong evidence of their statistical robustness (Figure 4G; Supplementary Figure 2; Supplementary Tables 6, 7).

3.5 Comparison of genomic alterations of the three AML subtypes

We utilized oncoplot to assess the genomic alterations in all subtypes of the TCGA LAML cohort (Figures 5A, B). Our analysis revealed that RUNX1, DNMT3A, and IDH2 mutations were most frequently observed in the immune dysregulation subtype, while being least prevalent in the other two subtypes (Figure 5B). Additionally, TTN, MUC16, KIT, ZFHX4, and FCGBP mutations were mainly observed in the cellular adaptation subtype (Figure 5B). Although the TMB of the immune dysregulation subtype was found to be higher than the remaining two subtypes, the difference was not statistically significant (Figure 5C). Similarly, the Burden of Copy Number gain was higher in the cellular adaptation subtype compared to the other two subtypes, whereas the Burden of Copy Number Loss was higher in the immune dysregulation subtype than in the other two subtypes; however, these differences failed to achieve statistical significance (Figure 5D). Notably, the copy number variation in KMT2A was specifically present in the immune dysregulation subtype (Figure 5E), highlighting the potential role of this gene in driving immune dysregulation in AML.

FIGURE 5

Figure 5 Comparison of genomic alterations of the three AML subtypes. (A) The genomic alterations among the three subtypes of the TCGA LAML cohort. (B) Mutation percentage of mostly mutated genes. (C) Tumor mutant burden difference among the three subtypes in the TCGA LAML cohort. (D) The Burden of Copy Number gain and the Burden of Copy Number Loss among the three subtypes in the TCGA LAML cohort. (E) Distinct Copy number alterations (CNA) profile among the three subtypes.

3.6 The three AML subtypes exhibited different immune statuses

The immune differences in the three subtypes of AML were explored through a comprehensive analysis of the immune microenvironment. We employed the ESTIMATE algorithm to calculate immune scores, revealing that patients in the immune dysregulation subtype exhibited significantly higher ESTIMATE scores compared to those in the other two subtypes (Figure 6A). However, the prognosis of the cellular adaptation subtype was found to be better than that of the other two subtypes, despite its higher immune score. Further, the activation states of immune-related pathways also were measured (Figure 6B). Several pathways crucial for immune function activation, including ‘Antigen processing and presentation,’ ‘NK cell cytotoxicity,’ and the ‘TCR signaling pathway,’ consistently exhibited upregulation in both the immune dysregulation subtype and the cellular adaptation subtype. Consistent with the ESTIMATE algorithm, all immune pathways within the tumor microenvironment get the highest activation level in the immune dysregulation subtype. We conducted a more detailed examination of the expression levels of 78 immunomodulators in each subtype (as shown in Figure 6C). Notably, the majority of differentially expressed immunomodulators exhibited the highest levels in both the immune dysregulation subtype and the cellular adaptation subtype. This includes key genes related to antigen presentation (HLA-B, HLA-C, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DQB1, HLA-DQB2, HLA-DRA, HLA-DRB1, HLA-DRB5, MICA and MICB). The total expressions of 78 immunomodulators were also the highest in the immune dysregulation subtype. However, it must be emphasized that unlike the cellular adaptation subtype, the suppressor proteins, such as CD274, CD276, PDCD1LG2, BTN3A1, BTN3A2, SLAMF7 and C10orf54 in the immune dysregulation subtype are also activated (Figures 6C, D; Supplementary Figure 3A). This is particularly remarkable because there was a significant difference in survival risk between the two subtypes (Figures 4B, F), but the immune dysregulation subtype and the cell adaptation subtype exhibits highly active immune functions. This compelled us to further analyze the key differences that determine the survival risk in these two subtypes and conducted pathway enrichment analysis on the top 100 differentially expressed genes in each subgroup. As expected, the immune dysregulation subtype not only exhibited highly activated immune responses, but also highly enriched pathways involved in negative regulation of immune responses (Figure 6E). Further GSEA confirmed our findings: in the immune dysregulation subtype, immune response and immune response processes were negatively regulated (Figure 6F), while the cell adaptation subtype exhibited enrichment in pathways related to DNA synthesis and protein synthesis (Figure 6G).

FIGURE 6

Figure 6 The three AML subtypes exhibited different immune statuses. (A) The ESTIMATE scores among the three subtypes in the TCGA LAML cohort. (B) The activation degree of 17 immune pathways in each tumor sample among the three subtypes. (C) The expression levels of 78 immunomodulators among the three subtypes. (D) Boxplot showing the activation status of suppressor proteins in different subtypes. (E) The pathway enrichment analysis on the top 100 differentially expressed genes in CA and ID subgroup. (F, G) GSEA results showing the activated signaling pathways in ID and CA subgroup. p-value > 0.05: ns, p-value < 0.05:*, 0.05 < p-value < 0.01: **, 0.01<p-value<0.001:***.

As a result, the clinical prognosis and biological characteristics of each subtype in two GEO cohorts closely resembled those of the corresponding subtypes identified in the TCGA LAML cohort (see Supplementary Figures 3A, B). Furthermore, akin to the TCGA LAML cohort, the levels of immune pathway activation and expression of immunomodulators were highest in the immune dysregulation subtype, reinforcing the robust and reliable nature of our identification of immune features among the three LAML subtypes (see Supplementary Figures 3A, B).

3.7 The three AML subtypes exhibited different drug resistance

To further investigate whether the immune dysregulation subtype had a better response rate to immune checkpoint blockade, we conducted a reprofiling of primary bone marrow (BM) samples obtained from 33 adult patients diagnosed with either newly diagnosed or relapsed/refractory AML who underwent treatment with AZA+Pembro (Azacitidine in Combination with Pembrolizumab, ClinicalTrials.gov NCT02845297). We investigated differentially expressed genes (DEGs) at the baseline between patients who eventually achieved Complete Remission (CR) and those who did not respond (NR). To precisely delineate immunotherapy responses distinguishing CR and NR samples in an unbiased manner, we conducted further refinement of the gene list between CR and NR samples. We employed Elastic Net penalized logistic regression (as shown in Figure 7A) to ensure the inclusion of only the most significant genes in the model, while also accounting for the high degree of correlation observed among certain immunotherapy responses. Collapsing our immunotherapy signature into a gene-weighted response score, we confirmed its ability to distinguish CR from NR samples through ROC analysis in the cohort (Figures 7B, C). Fortunately, the immune dysregulation subtype is more sensitive to AZA+Pembro therapy compared to other subtypes. (Figure 7D). This conclusion was validated in the other two publicly available GEO datasets, confirming that the immune dysregulation subtype is more sensitive to AZA+Pembro therapy compared to other subtypes (Figures 7E, F). Further, we employed molecular subtyping to stratify patients for precision medicine. To assess whether the three AML subtypes exhibited varying sensitivities to different drug profiles, we utilized the Genomics of Drug Sensitivity in Cancer (GDSC) database. Our findings revealed that different AML subtypes displayed sensitivity to distinct compounds, reinforcing the significance of molecular subtyping. (Figure 7G). We found that the Immune dysregulation patients were more sensitive to Pazopanib and Dasatinib (Figures 7H, I), and more tolerance to Vorinostat, Metformin, AZD6244 and Nutlin-3a etc. (Figure 7G). To investigate whether the drug sensitive result of the three subtypes of AML were repeatable in other datasets, we did a similar analysis on two GEO cohorts. As a result, the chemotherapy response of each subtype in two GEO cohorts were the same as the corresponding subtype identified in the TCGA LAML cohort (Supplementary Figure 4, Supplementary Tables 8, 9).

FIGURE 7

Figure 7 The three AML subtypes exhibited different drug resistance. (A) Features selected by Elastic-Net regression to differentiate between CR and NR samples. (B) ROC curves for the performance of GSE178926 cohort in predicting immunotherapy response. (C) Differentially expressed genes between CR and NR samples in GSE178926 cohort (Response score is the immunotherapy response score calculated by the elastic network model). (D) Boxplot showing the AZA+Pembro response score of the TCGA LAML cohort and (E, F) two GEO cohorts. (G) The heatmap showing the sensitivity of the three AML subtypes to different compounds. (H, I) Sensitivity of the three AML subtypes to Pazopanib, Dasatinib.

4 Discussion

Numerous investigations have firmly established the critical role of m⁶A modification in the etiology of AML and have identified singular molecular subtypes of AML based on changes in the m⁶A modifying enzymes (40, 41). Despite these advances, most of these studies focused on the bulk tissues, while largely neglecting the relevance of the cell hierarchy to AML initiation and progression(42). Neglecting the heterogeneity of m⁶A modification at the cell level may impair our understanding of fundamental mechanisms underlying AML, given that cellular properties significantly shape the microenvironment and ultimately determine the disease outcome (18, 20). In order to bridge this gap in knowledge, we undertook this investigation to incorporate information on m⁶A modification status and cell-level features in AML with prudence.

Several methods, including principal component analysis (PCA), single-sample gene set enrichment analysis (ssGSEA), and gene set variation analysis (43) (GSVA) are commonly used in bulk RNA-seq analysis to define the activation level of biological processes. These methods have also been applied extensively to investigate m⁶A modification activity in AML patients (44, 45). However, due to the significant sparsity issues in single-cell sequencing data, these conventional methods may not be well-suited for single-cell data analysis (46). Additionally, the 23 m⁶A modifying enzymes differ from conventional biological process gene sets as they consist of three types of enzymes: writers, readers and erasers, each performing a distinct function (47). These enzymes can form numerous combinations to participate in the m⁶A modification process in cells (48). Therefore, we used AUCell (22) to determine the m⁶A modification activity of each cell, taking into account the characteristics of single-cell data and the properties of m⁶A modifying enzymes. We categorized each cell into an active or inactive group based on its m⁶A modification activity. Notably, our results showed that m⁶A modification was more active in AML cells than in non-leukemic immune cells, highlighting the potential role of m⁶A modification in AML pathogenesis.

Transcription factor regulatory networks within cells are frequently constructed, offering an exciting opportunity for high-resolution identification of distinct transcriptional states and transitions, such as the m⁶A modification status (22). In a recent study, Zeng et al. employed the SCENIC network to delineate a subset of cells characterized by activation of CDK6 and E2F3, subsequently categorizing them as leukemia cells associated with cell cycle initiation (20). Subsequently, Guo et al. conducted VIPER to elucidate the transcription factor (TF) activity in AML progenitor cells (44). Our investigation aimed to elucidate the impact of differential m⁶A activity on regulatory networks. Notably, our findings demonstrated that m⁶A modifications exert a more pronounced effect on the transcription factor regulatory network within leukemia cells, with transcription factors in m⁶A-active leukemia cells exhibiting more conspicuous transcriptional regulatory activity. However, it is noteworthy that there is no significant correlation between transcription factor regulatory activity and m⁶A modification activity in immune cells. Specifically, TFs such as YY1, E2F1, and E2F8 were found to be closely associated with m⁶A activity within leukemia cells. Extensive evidence supports the role of YY1 and E2F1 in the regulation of cell differentiation and proliferation (34, 35). Additionally, we considered the topological changes in gene regulatory networks, which may reflect the complexity of intracellular gene expression(26). Importantly, AML cells with active m⁶A modifications and those with inactive m⁶A modifications exhibited similar gene-gene interaction densities. In contrast, gene interaction networks in non-AML immune cells were significantly influenced by m⁶A modifications. This suggests that m⁶A modifications in AML cells may impair their ability to coordinate gene networks effectively.

Numerous previous studies have established the significant influence of m⁶A modification on the prognosis of patients with AML(48). Nevertheless, the precise mechanism underlying this effect has not been thoroughly investigated. In this study, we present an innovative hypothesis suggesting that m⁶A modification primarily affects leukemic cells, thereby altering the prognostic status of AML patients. Our analysis revealed the existence of a high-risk subtype of AML termed “immune dysregulated subtype,” which exhibits RUNX1 mutations and KMT2A copy number variations (49, 50). Notably, our findings shed light on AML molecular pathogenesis and may offer crucial implications for personalized therapeutic strategies. To this end, we aim to determine the altered immune status and sensitivity to immunotherapy of this high-risk subtype. After identifying this high-risk subtype, our goal is to determine its altered immune status and sensitivity to immunotherapy. It is important to note that this subtype is active in the immune system and has shown sensitivity to immunotherapy. Our analysis revealed that Pazopanib and Dasatinib are highly sensitive to this subtype, and potential therapies for treating AML using these drugs have been reported in the literature (49, 51). This information on immunotherapy can be used to support clinical treatment and diagnosis.

Our study was the first to demonstrate the m⁶A modification information in AML at the single-cell level. We found that the m⁶A modification status of leukemic cells are strongly correlated with patient prognosis. By combining m⁶A modification information and cellular information, we identified a high-risk subtype and elucidated its clinical features. Furthermore, our study developed a high-precision machine learning composite model that integrates multiple mainstream machine learning algorithms. The robustness of our findings has been validated through one-by-one confirmation in a mainstream cohort. Nevertheless, we do acknowledge certain limitations of our study. In our study, we employed AUCell to assess the activation status of m⁶A-related genes based on their expression rank. This approach allowed us to form sets of genes that could potentially act in coordination to mediate m⁶A modifications. However, it is important to note that m⁶A enzymes often function in a collaborative and interdependent manner, and pinpointing precise coordination mechanisms can be challenging (48). Additionally, distinguishing between immune cells and leukemia cells in bulk RNA-seq data can be challenging because these cell types may share common biological processes or differential gene expression profiles (20).

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

ZL: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Writing – original draft, Writing – review & editing. XL: Data curation, Formal Analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. LW: Funding acquisition, Supervision, Writing – review & editing. HZ: Data curation, Investigation, Visualization, Writing – original draft. SW: Visualization, Writing – original draft. GY: Funding acquisition, Project administration, Supervision, Writing – review & editing. DW: Funding acquisition, Project administration, Supervision, Writing – review & editing. JC: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing. JH: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. The National Natural Science Foundation of China (81770216, 81900180). The Ministry of Science and Technology, PR China, (2019YFE0119500). Henan Province Science and Technology Project (212102310894).

Acknowledgments

The present study acknowledges with heartfelt appreciation the profound gratitude to Dr. Baofa Sun from Nankai University for his sagacious counsel and valuable suggestions that have greatly enriched the manuscript. The authors also extend their gratitude to Wei Geng from Huazhong University of Science and Technology for his invaluable contribution in data analysis and visualization. The remarkable expertise and support of these individuals have undoubtedly contributed to the sublimity of this research project, and the authors remain sincerely grateful. We are also grateful for the assistance of The High-Performance Computing Center of Henan Normal University.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1281687/full#supplementary-material

Supplementary Figure 1 | Effect of m⁶A modification status on active and inactive m⁶A methylated cells in AML. (A) t-SNE visualization of 13,653 single-cell transcriptomes (points), with similar cells positioned closer together. Cell type annotation were obtained from Zeng et al. (20). (B) Heatmap showing the expression of 23 selected cell-type-specific genes (rows) across 13,653 single cells ordered by cell-type annotations as shown in A (columns). (C) The stacked bar chart displays the proportions of m⁶A modification-active and inactive cells in leukemia and immune cells. (D) Line plots show the number of cell types of transcription factors with significantly different activity between m⁶A-active cells and m⁶A-inactive cells in immune cells (blue) and leukemic cells (orange). (E) The network plot shows the regulatory relationship between transcription factors and differentially expressed genes (m⁶A-active cells vs. m⁶A-inactive cells) in Mono-like blasts. (F) The workflow shows the screening process of cell-specific m⁶A-related genes. (G) The network plot shows the pathways the pathways enriched by Mono-like blasts-specific m⁶A-related genes. (H) Boxplot shows the m⁶A modification activity of LSPCs in the three states (Quiescent, Primed, Cycle).

Supplementary Figure 2 | The Robustness of composite machine learning classification models for the identification of subtypes in the GEO dataset. Kaplan-Meier curves of overall survival (OS) among the three subtypes classified by different machine learning models in the GSE106291 cohort (up) and the GSE146173 cohort (down).

Supplementary Figure 3 | Identification of immune features of the three AML subtypes in the GEO dataset. (A) The ESTIMATE scores and the activation degree of 17 immune pathways among the three subtypes in the GSE106291 cohort and GSE146173 cohort. (B) The expression levels of 78 immunomodulators among the three subtypes in the GSE106291 cohort and GSE146173 cohort.

Supplementary Figure 4 | The sensitivity of the three AML subtypes to different compounds in the GEO dataset. (A) The heatmap showing the sensitivity of the three AML subtypes to different compounds (GSE106291). (B) Sensitivity of the three AML subtypes to Pazopanib, Dasatinib (GSE106291). (C) The heatmap showing the sensitivity of the three AML subtypes to different compounds (GSE146173). (D) Sensitivity of the three AML subtypes to Pazopanib, Dasatinib (GSE146173).

Supplementary Table 1 | Results of differential transcription factor activity analysis.

Supplementary Table 2 | Results of BigScale2-based gene regulatory network analysis.

Supplementary Table 3 | Cell-specific m⁶A-associated gene set.

Supplementary Table 4 | The classification information of the TCGA.

Supplementary Table 5 | Gene set used to construct models that identify molecular subtypes.

Supplementary Table 6 | Accuracy statistics for multiple machine learning models.

Supplementary Table 7 | The classification information of the validation set.

Supplementary Table 8 | Results of drug susceptibility prediction.

Supplementary Table 9 | IC50 range table for each cell line corresponding to the drug.

Supplementary File 1 | Count matrix of AML patients from TCGA database.

Supplementary File 2 | Count matrix of AML patients from GSE106291 database.

Supplementary File 3 | Count matrix of AML patients from GSE146173 database.

References

1. Xuan L, Liu Q. Maintenance therapy in acute myeloid leukemia after allogeneic hematopoietic stem cell transplantation. J Hematol Oncol (2021) 14:4. doi: 10.1186/s13045-020-01017-7

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell (2017) 169:1187–200. doi: 10.1016/j.cell.2017.05.045

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Huang H, Weng H, Chen J. m(6)A modification in coding and non-coding RNAs: roles and therapeutic implications in cancer. Cancer Cell (2020) 37:270–88. doi: 10.1016/j.ccell.2020.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Deng X, Su R, Weng H, Huang H, Li Z, Chen J. RNA N(6)-methyladenosine modification in cancers: current status and perspectives. Cell Res (2018) 28:507–17. doi: 10.1038/s41422-018-0034-6

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Ma S, Chen C, Ji X, Liu J, Zhou Q, Wang G, et al. The interplay between m6A RNA methylation and noncoding RNA in cancer. J Hematol Oncol (2019) 12:121. doi: 10.1186/s13045-019-0805-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Li X, Ma S, Deng Y, Yi P, Yu J. Targeting the RNA m(6)A modification for cancer immunotherapy. Mol Cancer (2022) 21:76. doi: 10.1186/s12943-022-01558-0

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Barbieri I, Tzelepis K, Pandolfini L, Shi J, Millan-Zambrano G, Robson SC, et al. Promoter-bound METTL3 maintains myeloid leukemia by m(6)A-dependent translation control. Nature (2017) 552:126–31. doi: 10.1038/nature24678

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Vu LP, Pickering BF, Cheng Y, Zaccara S, Nguyen D, Minuesa G, et al. The N(6)-methyladenosine (m(6)A)-forming enzyme METTL3 controls myeloid differentiation of normal hematopoietic and leukemia cells. Nat Med (2017) 23:1369–76. doi: 10.1038/nm.4416

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Weng H, Huang H, Wu H, Qin X, Zhao BS, Dong L, et al. METTL14 Inhibits Hematopoietic Stem/Progenitor Differentiation and Promotes Leukemogenesis via mRNA m(6)A Modification. Cell Stem Cell (2018) 22:191–205.e199. doi: 10.1016/j.stem.2017.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Weng H, Huang F, Yu Z, Chen Z, Prince E, Kang Y, et al. The m(6)A reader IGF2BP2 regulates glutamine metabolism and represents a therapeutic target in acute myeloid leukemia. Cancer Cell (2022) 40:1566–1582.e1510. doi: 10.1016/j.ccell.2022.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhang N, Shen Y, Li H, Chen Y, Zhang P, Lou S, et al. The m6A reader IGF2BP3 promotes acute myeloid leukemia progression by enhancing RCC2 stability. Exp Mol Med (2022) 54:194–205. doi: 10.1038/s12276-022-00735-x

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Feng M, Xie X, Han G, Zhang T, Li Y, Li Y, et al. YBX1 is required for maintaining myeloid leukemia cell survival by regulating BCL2 stability in an m6A-dependent manner. Blood (2021) 138:71–85. doi: 10.1182/blood.2020009676

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yoon KJ, Ringeling FR, Vissers C, Jacob F, Pokrass M, Jimenez-Cyrus D, et al. Temporal control of mammalian cortical neurogenesis by m(6)A methylation. Cell (2017) 171:877–889 e817. doi: 10.1016/j.cell.2017.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Huang Y, Su R, Sheng Y, Dong L, Dong Z, Xu H, et al. Small-molecule targeting of oncogenic FTO demethylase in acute myeloid leukemia. Cancer Cell (2019) 35:677–691.e610. doi: 10.1016/j.ccell.2019.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Herold T, Jurinovic V, Batcha AMN, Bamopoulos SA, Rothenberg-Thurley M, Ksienzyk B, et al. A 29-gene and cytogenetic score for the prediction of resistance to induction treatment in acute myeloid leukemia. Haematologica (2018) 103:456–65. doi: 10.3324/haematol.2017.178442

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Bamopoulos SA, Batcha AMN, Jurinovic V, Rothenberg-Thurley M, Janke H, Ksienzyk B, et al. Clinical presentation and differential splicing of SRSF2, U2AF1 and SF3B1 mutations in patients with acute myeloid leukemia. Leukemia (2020) 34:2621–34. doi: 10.1038/s41375-020-0839-4

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Rutella S, Vadakekolathu J, Mazziotta F, Reeder S, Yau TO, Mukhopadhyay R, et al. Immune dysfunction signatures predict outcomes and define checkpoint blockade-unresponsive microenvironments in acute myeloid leukemia. J Clin Invest (2022) 132. doi: 10.1172/JCI159579

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Van Galen P, Hovestadt V, Wadsworth Ii MH, Hughes TK, Griffin GK, Battaglia S, et al. Single-cell RNA-seq reveals AML hierarchies relevant to disease progression and immunity. Cell (2019) 176:1265–1281.e1224. doi: 10.1016/j.cell.2019.01.031

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell (2021) 184:3573–3587.e3529. doi: 10.1016/j.cell.2021.04.048

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zeng AGX, Bansal S, Jin L, Mitchell A, Chen WC, Abbas HA, et al. A cellular hierarchy framework for understanding heterogeneity and predicting drug response in acute myeloid leukemia. Nat Med (2022) 28:1212–23. doi: 10.1038/s41591-022-01819-x

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Finak G, Mcdavid A, Yajima M, Deng J, Gersuk V, Shalek AK, et al. MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol (2015) 16:278. doi: 10.1186/s13059-015-0844-5

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods (2017) 14:1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Kumar N, Mishra B, Athar M, Mukhtar S. Inference of gene regulatory network from single-cell transcriptomic data using pySCENIC. Methods Mol Biol (2021) 2328:171–82. doi: 10.1007/978-1-0716-1534-8_10

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Lambert SA, Jolma A, Campitelli LF, Das PK, Yin Y, Albu M, et al. The human transcription factors. Cell (2018) 172:650–65. doi: 10.1016/j.cell.2018.01.029

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Iacono G, Mereu E, Guillaumet-Adkins A, Corominas R, Cuscó I, Rodríguez-Esteban G, et al. bigSCale: an analytical framework for big-scale single-cell data. Genome Res (2018) 28:878–90. doi: 10.1101/gr.230771.117

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Iacono G, Massoni-Badosa R, Heyn H. Single-cell transcriptomics unveils gene regulatory network plasticity. Genome Biol (2019) 20:110. doi: 10.1186/s13059-019-1713-4

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumor purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, et al. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci Data (2018) 5:180015. doi: 10.1038/sdata.2018.15

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9:e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Pulikkan JA, Dengler V, Peramangalam PS, Peer Zada AA, Müller-Tidow C, Bohlander SK, et al. Cell-cycle regulator E2F1 and microRNA-223 comprise an autoregulatory negative feedback loop in acute myeloid leukemia. Blood (2010) 115:1768–78. doi: 10.1182/blood-2009-08-240101

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Pang Y, Zhao Y, Wang Y, Wang X, Wang R, Liu N, et al. TNFAIP8 promotes AML chemoresistance by activating ERK signaling pathway through interaction with Rac1. J Exp Clin Cancer Res (2020) 39:158. doi: 10.1186/s13046-020-01658-z

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Jiang F, Wang XY, Wang MY, Mao Y, Miao XL, Wu CY, et al. An immune checkpoint-related gene signature for predicting survival of pediatric acute myeloid leukemia. J Oncol (2021) 2021:5550116. doi: 10.1155/2021/5550116

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Fiskus W, Boettcher S, Daver N, Mill CP, Sasaki K, Birdwell CE, et al. Effective Menin inhibitor-based combinations against AML with MLL rearrangement or NPM1 mutation (NPM1c). Blood Cancer J (2022) 12:5. doi: 10.1038/s41408-021-00603-3

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Li M, Li M, Xia Y, Li G, Su X, Wang D, et al. HDAC1/3-dependent moderate liquid-liquid phase separation of YY1 promotes METTL3 expression and AML cell proliferation. Cell Death Dis (2022) 13:992. doi: 10.1038/s41419-022-05435-y

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Ou ZY, Wang K, Shen WW, Deng G, Xu YY, Wang LF, et al. Oncogenic FLT3 internal tandem duplication activates E2F1 to regulate purine metabolism in acute myeloid leukaemia. Biochem Pharmacol (2023) 210:115458. doi: 10.1016/j.bcp.2023.115458

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Kumar D, Whiteside TL, Kasid U. Identification of a novel tumor necrosis factor-alpha-inducible gene, SCC-S2, containing the consensus sequence of a death effector domain of fas-associated death domain-like interleukin- 1beta-converting enzyme-inhibitory protein. J Biol Chem (2000) 275:2973–8. doi: 10.1074/jbc.275.4.2973

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Ibrado AM, Kim CN, Bhalla K. Temporal relationship of CDK1 activation and mitotic arrest to cytosolic accumulation of cytochrome C and caspase-3 activity during Taxol-induced apoptosis of human AML HL-60 cells. Leukemia (1998) 12:1930–6. doi: 10.1038/sj.leu.2401218

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Grey W, Ivey A, Milne TA, Haferlach T, Grimwade D, Uhlmann F, et al. The Cks1/Cks2 axis fine-tunes Mll1 expression and is crucial for MLL-rearranged leukaemia cell viability. Biochim Biophys Acta Mol Cell Res (2018) 1865:105–16. doi: 10.1016/j.bbamcr.2017.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Bennett JM, Catovsky D, Daniel MT, Flandrin G, Galton DA, Gralnick HR, et al. Proposals for the classification of the acute leukemias. French-American-British (FAB) co-operative group. Br J Haematol (1976) 33:451–8. doi: 10.1111/j.1365-2141.1976.tb03563.x

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Xu ZJ, Wen XM, Zhang YC, Jin Y, Ma JC, Gu Y, et al. m6A regulator-based methylation modification patterns and characterization of tumor microenvironment in acute myeloid leukemia. Front Genet (2022) 13:948079. doi: 10.3389/fgene.2022.948079

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Yuan S, Cong Z, Ji J, Zhu L, Jiang Q, Zhou Y, et al. Analysis of m6A-related signatures associated with the tumor immune microenvironment and predict survival in acute myeloid leukemia. Ann Transl Med (2022) 10:902. doi: 10.21037/atm-22-3858

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Kuksin M, Morel D, Aglave M, Danlos FX, Marabelle A, Zinovyev A, et al. Applications of single-cell and bulk RNA sequencing in onco-immunology. Eur J Cancer (2021) 149:193–210. doi: 10.1016/j.ejca.2021.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

44. Wu J, Xiao Y, Sun J, Sun H, Chen H, Zhu Y, et al. A single-cell survey of cellular hierarchy in acute myeloid leukemia. J Hematol Oncol (2020) 13:128. doi: 10.1186/s13045-020-00941-y

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Du A, Wu X, Gao Y, Jiang B, Wang J, Zhang P, et al. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in acute myeloid leukemia. Front Immunol (2021) 12:789914. doi: 10.3389/fimmu.2021.789914

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Holland CH, Tanevski J, Perales-Patón J, Gleixner J, Kumar MP, Mereu E, et al. Robustness and applicability of transcription factor and pathway analysis tools on single-cell RNA-seq data. Genome Biol (2020) 21:36. doi: 10.1186/s13059-020-1949-z

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. BioMed Pharmacother (2019) 112:108613. doi: 10.1016/j.biopha.2019.108613

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Deng LJ, Deng WQ, Fan SR, Chen MF, Qi M, Lyu WY, et al. m6A modification: recent advances, anticancer targeted drug discovery and beyond. Mol Cancer (2022) 21:52. doi: 10.1186/s12943-022-01510-2

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Al-Harbi S, Aljurf M, Mohty M, Almohareb F, Ahmed SOA. An update on the molecular pathogenesis and potential therapeutic targeting of AML with t(8;21)(q22;q22.1);RUNX1-RUNX1T1. Blood Adv (2020) 4:229–38. doi: 10.1182/bloodadvances.2019000168

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Du W, Xu A, Huang Y, Cao J, Zhu H, Yang B, et al. The role of autophagy in targeted therapy for acute myeloid leukemia. Autophagy (2021) 17:2665–79. doi: 10.1080/15548627.2020.1822628

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Kessler T, Koschmieder S, Schliemann C, Crysandt M, Mikesch JH, Von Stillfried S, et al. Phase II clinical trial of pazopanib in patients with acute myeloid leukemia (AML), relapsed or refractory or at initial diagnosis without an intensive treatment option (PazoAML). Ann Hematol (2019) 98:1393–401. doi: 10.1007/s00277-019-03651-9

PubMed Abstract | CrossRef Full Text | Google Scholar