Genetic and transcriptional alterations of m6A regulators in AML
We investigated a total of 23 m6A regulators including 8 writers, 2 erasers, and 13 readers in this study (Table S1). The differences in copy number variation (CNV) of all regulators were explored. The CNV of m6A regulators was prevalent in AML. The CNV was found to be focused on the deletion predominantly, while YTHDF2, YTHDF3, FTO, YTHDF1, and ELAVL1 had more prevalent amplification frequency (Figure 1A, Table S2). The locations of CNV alteration in 23 m6A regulators on their respective chromosomes were displayed in Figure 1B (Table S3, Table S4). Moreover, we compared AML samples to normal samples using PCA analysis based on paired tumor-normal specimens and found that AML samples were clearly distinguished from normal samples by the 23 m6A regulators (Figure S1A). After evaluating the incidence of somatic mutations of 23 m6A regulators in AML, we concluded that the alteration frequency was only 1.04% among 193 samples with 2 missense mutations identified in METTL3 and RBMX (Figure 1C). Notably, m6A regulators are widely known for their collaboration in cancer progression. The relevance of the co-expression of regulators was therefore analyzed, and LRPPRC and VIRMA owned a significant correlation with the highest correlation coefficient (R = 0.66) (Figure 1D, Table S5 and Table S6). Intrigued by these findings, we further compared the mRNA expression levels between AML and normal tissues. The results displayed all of the 23 m6A regulators showed significant differences between AML and normal tissues except for YTHDF2 and METTL3 (Figure 1E). Beyond that, m6A highlighted distinct patterns in different immune cell types, myeloid malignancies along with lymphoid malignancies (Figure S1B, Figure S1C). The observed clustering suggested that the cell-of-origin and different disease groups influenced the repertoire of m6A regulators. YTHDF3 presented relatively higher expression compared with the other 22 m6A regulators in M2-macrophage. At the same time, AML highly expressed IGF2BP2.
Identification of m6A subtypes in AML
151 patients from the TCGA dataset with available OS data and clinical information were enrolled in our study for further analysis (Table S7). The results of univariate Cox regression analysis revealed the prognostic values of 23 m6A regulators in patients with AML (Table S8). Figure 2A showed the comprehensive landscapes of 23 m6A regulators concerning their interactions, connections, and prognostic values in AML in the form of a m6A regulator network. Significant correlations existed among writers, erasers, and readers. We further explored the expression pattern of m6A regulators in AML and conducted the R package of ConsensusClusterPlus to classify the patients based on the expression level of 23 m6A regulators. It was clear that the clustering model was optimal when K = 3 during the process of K ranging from2 to 9 (Figure S2A- S2C). The entire cohort was eventually allocated into three distinct modification patterns consisting of 111 cases in cluster A, 33 cases in cluster B, and 7 cases in cluster C (Figure 2B). The Kaplan-Meier curves for the three main m6A modification subtypes revealed the worst overall survival in patients with m6Acluster A (P = 0.003, Figure 2B). The prognosis of m6Acluster-B and m6Acluster C was similar in the Kaplan-Meier analysis, so we combined m6Acluster B with m6Acluster C and m6Acluster A was compared with it using univariate analyses. m6Acluster A had significantly poorer prognosis (HR = 0.3998, P < 0.01, Table S9). Moreover, we validated the biological characteristics of three m6A subtypes. The m6Acluster B was preferentially related to FAB subtypes M3, and the m6Acluster C subtype was markedly correlated with intermediate cytogenetic risk and FTO (Figure 2C).
Associations of TME with three m6A subtypes
The association of 23 TILs with three m6A clusters was investigated by using the CIBERSORT algorithm in the TME of AML. We observed profound differences in immune cell infiltration including activated B cell, CD56dim natural killer cell, immature B cell, mast cell, regulatory T cell, type 1 T helper cell, type 17 T helper cell, and type 2 T helper cell among the three subtypes (Figure 3A). In addition, we performed GSVA enrichment analysis to detect the biological characteristics of three distinct m6A modification patterns. GSVA enrichment analysis showed that m6Acluster A was significantly enriched in metabolic pathways such as biosynthesis of unsaturated fatty acids, propanoate metabolism, pyruvate metabolism, citrate cycle TCA cycle, glyoxylate and dicarboxylate metabolism, pyrimidine metabolism, purine metabolism, cysteine and methionine metabolism (Figure 3B and Table S10). m6Acluster B presented enrichment pathways highly related to metabolic pathways including linoleic acid metabolism (Figure 3B and Table S10), and carcinogenic activation pathways such as ECM receptor interaction (Figure 3C and Table S11). Likewise, m6Acluster C was prominently associated with metabolic pathways which encompassed cysteine and methionine metabolism, glyoxylate and dicarboxylate metabolism (Figure 3C and Table S11). What’s more, m6Acluster A correlated with immune pathways and immune-related diseases simultaneously, including antigen processing and presentation, chemokine signaling pathway activation, asthma, and systemic lupus erythematosus (Figure 3D and Table S12). Figure 3E illustrated cell type assignment for the scRNA expressed genes in AML.
m6A phenotype-related DEGs in AML
To further investigate the underlying m6A-related transcriptional expression differences within three m6A modification patterns, we did some further analysis. Actually, we applied a Venn diagram to illustrate the DEGs among the three m6A modification patterns and 70 m6A phenotype-related overlap genes were presented (Figure S3A, Table S13-S16). To further elaborate on the function of DEGs, we carried out a GO enrichment analysis. The result elucidated that enrichment of the biological processes related to embryonic organ morphogenesis/development (Figure S3B). Table S17 showed the expression level of 70 overlap DEGs. Univariate Cox regression study of 70 genes was analyzed and 25 genes with P-value less than 0.05 were screened out (Table S18). Then we performed an unsupervised clustering algorithm to classify the entire cohort into three m6A gene signature subtypes, named as m6A gene cluster A-C, respectively (Figure 4A). Further survival analysis indicated three m6A genomic phenotypes showed significant prognostic differences in AML samples, of which 45 cases were in gene cluster A, 47 cases in gene cluster B, and 59 cases in gene cluster C (P = 0.002, Log-rank test, Figure 4B). Patients in gene cluster C were proved to be related to a better prognosis, while patients in gene cluster A and gene cluster B were associated with the outcome of poorer prognosis. We also observed that m6A gene cluster C patterns were associated with good cytogenetic risk, while m6A gene cluster A presented poor cytogenetic risk and m6A gene cluster B showed intermediate cytogenetic risk. Patients with alive status were largely concentrated in the m6A gene cluster A patterns. Interestingly, VSTM1 and M1AP were markedly decreased in the m6A gene cluster B and increased in the m6A gene cluster C (Figure 4C). The three m6A gene clusters indicated significant differences in the expression of VIRMA, IGF2BP2, FMR1, IGF2BP3, ELAVL1, and YTHDC2 (Figure 4D).
In addition, the immune cell infiltration of activated B cell, activated CD4 T cell, activated dendritic cell, macrophage, NK cell, plasmacytoid dendritic cell and type 2 T helper cell showed statistically significant differences across the three m6A gene clusters (Figure 5A). Meanwhile, m6A gene cluster A was prominently associated with metabolic pathways including mucin type O glycan biosynthesis, nitrogen metabolism, heme biosynthesis, lysine degradation, epinephrine biosynthesis, glycogen degradation, pantothenate and CoA biosynthesis, vitamin B6 metabolism, retinoic acid metabolism, phenylalanine tyrosine and tryptophan biosynthesis (Figure 5B, Figure 5C and Table S19, Table S20). m6A gene cluster B was enriched in the pathway associated with retinoid metabolism, glyoxylate and dicarboxylate metabolism, steroid biosynthesis, nicotinamide adenine dinucleotide biosynthesis, ubiquinone, and other terpenoid quinone biosynthesis, folate biosynthesis, oxidative phosphorylation, glycogen degradation, cysteine and methionine metabolism, phenylalanine tyrosine and tryptophan biosynthesis, pantothenate and CoA biosynthesis, pentose phosphate, primary bile acid biosynthesis, and retinoic acid metabolism (Figure 5C, Figure 5D, and Table S20, Table S21). m6A gene cluster C presented enrichment pathways prominently associated with retinol metabolism, ascorbate and aldrate metabolism, steroid hormone biosynthesis and metabolism, porphyrin and chlorophyll metabolism, metabolism of xenobiotics by cytochrome P450, drug metabolism by cytochrome P450, epinephrine biosynthesis, lysine degradation, and N glycan biosynthesis (Figure 5C, Figure 5D, and Table S20, Table S21).
Construction of m6A score and the relevance of clinical features.
Considering the individual heterogeneity and complexity of m6A modification, a scoring system based on the signature genes related to m6A, named m6Ascore, was established to quantify the m6A modification pattern of single AML patient (Table S22). Figure 6A illustrated an alluvial diagram of three m6A clusters, three m6A gene clusters, m6Ascore, and OS status. Remarkably, m6A cluster B exhibited the highest m6A score, followed by m6A cluster A and m6A cluster C (Figure S4A). The results suggested that the m6Ascore of m6A gene cluster B was the lowest, while that of m6A gene cluster C was the highest (Figure S4B). Spearman analysis was applied to better understand the correlation between immune cells and the m6Ascore (Figure S4C and Table S23). m6Ascore was found to be positively associated with regulatory T cell. To further shed light on the impact of m6Ascore in AML, we next explored its potential value for predicting clinical outcomes. Overall survival curves of Kaplan-Meier analyses revealed that patients with high m6Ascore had a significantly prominent prognosis (P = 0.001, Log-rank test, Figure 6B). Similarly, by univariate analyses of the TCGA data set, the high m6Ascore was significantly related to better OS (HR = 2.442, P < 0.01, Table S24). To explore the distribution variations of the somatic mutations between low and high m6Ascore, we evaluated samples from the TCGA dataset using the maftools package. The top five mutated genes in the high m6Ascore group were CEBPA, KIT, FLT3, TET2, and NRAS, in parallel, the top five mutated genes in the low m6Ascore group were FLT3, DNMT3A, NPM1, IDH2, and IDH1 (Figure 6C, Figure 6D). The Forest plot revealed the alteration frequencies of the top five mutated genes in two groups using the chi-square test (Figure S4D). Delightfully, we found m6Ascore correlated strongly with tryptophan metabolism, pyruvate metabolism, glycogen degradation, cysteine and methionine metabolism, arachidonic acid metabolism, and ADP ribosylation (Figure 6E).