Public transcriptome and clinical data were collected from the Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database (24) up to May 28, 2022. Three eligible AML cohorts (TCGA-LAML, GSE12714 (25), and GSE37642 (26)) were selected for further investigation. Patients without survival data were excluded. For the TCGA dataset, the RNA transcriptome data were acquired from the Genomic Data Commons (GDC) database via the TCGAbiolinks package (27) and normalized to fragments per kilobase of exon model per million mapped fragments (FPKM) values, while the clinical data were acquired via the UCSC Xena database (28). For GEO datasets, GEOquery package was used to obtain the normalized matrix files and clinical data (29). The baseline information of the included patients is summarized in Table S1 .

MMRGs and immune-related genes were obtained from GeneCards database (30). A total of 50 MMRGs were obtained with a relevance score >5 ( Table S2 ), and 27 immune-related genes were obtained with a relevance score >20. Twenty DNA repair genes were also obtained from published literature (31) ( Table S3 ).

Somatic mutation profile of AML patients was obtained from the TCGA database. Mutation signatures of MMRGs were summarized and presented using the R package Maftools (32). The chromosomal locations of MMRGs were presented using the R package RCircos (33).

The relative abundance of a gene set enriched in a sample is widely detected using the single sample gene set enrichment analysis (ssGSEA) algorithm. In this study, MMs of patients in the three AML datasets were calculated using GSVA package (34) through the ssGSEA algorithm based on the expression data of MMRGs. The patients in the three datasets were separated into high- and low-MMs groups based on the cut-off values. Cut-off values with the lowest significant (P-value< 0.05) log-rank P-value were chosen as the optimal cut-off values (35). Two-group comparisons were analyzed with Wilcoxon rank sum test. One-way ANOVA and Tukey’s Honestly Significance Difference test were used to compare MMs among different French–American–British (FAB) subgroups.

First, the limma package was used to determine the differentially expressed genes (DEGs) between high- and low-MMs groups in the TCGA-LAML dataset (36). DEGs with the adjusted P-value< 0.05 and |logFC|> 0 were included for analysis. The DEGs were used for WGCNA via the WGCNA package (37) to determine the gene modules related to MMs with highly correlated expression profiles across samples. WGCNA was run with networkType = “signed”, minModuleSize = 30, and softpower = 21.

The patients were classified in TCGA-LAML datasets using the R package ConsensusClusterPlus (38) based on the expression of MMRGs via the consensus clustering algorithm (39). The cluster size was set between 2 and 8. About 80% of the samples were randomly selected 50 times to ensure classification stability. Wilcoxon rank-sum test was used to determine the expression differences of MMRGs between different subtypes of AML patients. P-value< 0.05 was considered statistically significant.

The MMRGs with prognostic value (P-value< 0.1) were selected among the 17 MMRGs via univariate Cox analysis for the least absolute and selection operator (LASSO) regression. LASSO regression was performed via the glmnet package (40) based on prognostic MMRGs with parameter family = “cox” and ten-fold cross-validation, and repeated 1,000 times to prevent overfitting. A prognosis model was built via the multivariate Cox regression based on the MMRGs selected through LASSO regression. Per-patient risk score was calculated based on the prognostic model formula.

Patients were categorized into high- and low-risk groups based on the optimal cut-off value. An internal verification nomogram for overall survival (OS) prediction was developed using the rms package to verify the prognostic accuracy of the model. Additionally, calibration curves and decision curve analysis (DCA) were used to verify the model’s performance (41). Further external validation was performed using the GSE12417 cohort.

Bone marrow specimens were fixed in paraformaldehyde, decalcified, and embedded in paraffin. Tissue sections were deparaffinized and rehydrated, and then stained with H&E, or incubated with antibodies, specifically ECHS1 (Cat No. 11305-1-AP; proteintech, China), and NUDFS2 (Cat No. R27071; Zenbioscience, China). Antigens were detected using diaminobenzidine. Images were captured using a microscope equipped with a SPOT camera under high-power fields (400x). The mean intensity optical density (IOD) was used to quantify the immunohistochemical expression measured by ImageJ. This study was performed according to the Declaration of Helsinki, and received approval from the Ethics Committee of Zhongnan Hospital of Wuhan University.

The limma package was used to determine the DEGs between high- and low-risk groups of the TCGA-LAML cohort (36) for further analysis. Adjusted P-value<0.05 and |logFC| >1 were set as the inclusion criteria.

Gene Ontology (GO) analysis (42), Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis (43), and Disease Ontology (DO) analysis (44) were performed using the clusterProfiler package to reveal the characteristics of the DEGs (45). Besides, Gene Set Enrichment Analysis (GSEA) (46) was performed using the clusterProfiler package (45). The c2.all.v7.5.1.entrez.gmt gene set was acquired via the Molecular signatures database 3.0 (MSigDB) (47) as the reference gene set. Significantly enriched pathways were determined at P-value< 0.05 and FDR< 0.25.

PPI networks were constructed via the STRING database based on the DEGs (48) (confidence level 0.400) and visualized using Cytoscape software. MCC (Matthews Correlation Coefficient metric) (49), MNC (the maximal neighbourhood coefficient), and Degree algorithms were used to identify the hub genes.

Tumor immune activity in AML patients was assessed using the ESTIMATE package based on TCGA-LAML expression profiles (50). The Tumor Immune Dysfunction and Exclusion (TIDE) immune scores were estimated through the TIDE database (51) to predict patients’ immunotherapy responses. Also, the infiltration abundance of 22 immune cells was determined using the CIBERSORTx database (52) based on patients’ transcriptome data. Besides, enrichment scores of 28 tumor-infiltration-associated immune cells were calculated using the ssGSEA algorithm in the GSVA package to represent the immune infiltration levels (53) in each sample. The correlation of different immune cells within two groups was calculated using the spearman algorithm. The expression matrix of the TCGA-LAML dataset was combined to calculate correlations between immune cells and hub genes in different groups. The correlation heat map was developed through the pheatmap package.

RBPs interacting with hub genes were predicted using ENCORI database (54). mRNA-RBPs interaction pairs were then screened with clipExpNum > 2 and clipIDnum > 2 as screening criteria and plotted mRNA-RBPs interaction network. In addition, potential drugs or small molecule compounds interacting with hub genes were predicted using the drug-gene interaction database (DGIdb) (55). Cytoscape software was used to visualize mRNA-RBPs and mRNA-drugs interaction networks.

The relationships between hub genes and drug sensitivity were explored based on their transcriptome profiles and drug-sensitive profiles downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) database (56), the Cancer Cell Line Encyclopedia (CCLE) database (57), and the CellMiner database (58).

R software (Version 4.1.2) was used for data processing and statistical analyses. Students t-test and Wilcoxon rank sum test were used for two-group comparisons, while Kruskal-Wallis tests were used for three or more group comparisons. Chi-square or Fisher’s exact tests were used to compare categorical variables. The survival package was used to perform survival analysis. Kaplan-Meier (KM) curves were used to express the survival difference between the two groups, and its significance was measured using a log-rank test. Spearman correlation analysis was used to calculate the correlation of distinct variables if not specified. The P-value was adjusted by the Benjamini and Hochberg method.

A general overview of the study is presented in Figure 1 . The three datasets (TCGA-LAML, GSE12417, and GSE37642) were first normalized using the limma package, with the TCGA-LAML dataset containing 151 samples ( Figures S1A, B ), the GSE12417 dataset containing 163 samples ( Figures S1C, D ), and the GSE37642 dataset containing 422 samples ( Figures S1E, F ). The expression profiles of the three datasets were consistent among samples after normalization ( Figures S1A, F ). The mutation profile of the 31 MMRGs in AML samples was analyzed to evaluate the mutation characteristic of the 31 MMRGs (BCS1L, COX10, COX15, DGUOK, ECHS1, ETHE1, FBXL4, LIAS, MPV17, NDUFA1, NDUFA13, NDUFAF5, NDUFS1, NDUFS2, NDUFS4, NDUFS7, NDUFV1, NFU1, POLG, POLG2, SDHA, SDHD, SLC25A4, SUCLA2, SUCLG1, SURF1, TACO1, TIMM8A, TMEM70, TRIT1, TTC19). The MMRGs had few single nucleotide polymorphisms (SNPs) in patients in the TCGA-LAML cohorts, only three genes (SUCLA2, SURF1, and POLG) with somatic mutations among 31 MMRGs in three AML patiens ( Figure 2A , Figure S2A ), mainly missense mutation. Single nucleotide variant (SNV) with C>T was the most prevalent, followed by T>A ( Figure 2A ).

RCircos package was used to annotate the chromosome location of 31 MMRGs to analyze the location of the 31 MMRGs on human chromosomes ( Figure 2B ). The MMRGs were mainly distributed in the 2, 5, 10, 11, and 17 chromosomes, of which most were distributed on the second chromosome (6 MMRGs) ( Figure 2B ). MMRGs located on the same chromosome were closely related at the genome level.

To investigate the overall expression of MMRGs in AML, MMs of each patient were calculated to represent the mitochondrial metabolism levels of AML patients. The patients were divided into two groups after combining the prognostic information and MMs of AML patients in the three datasets (excluding the AML patients lacking OS and OS. time) based on the MMs levels. Patient survival outcomes were significantly different per grouping (P-value< 0.05) depending on the best cut-off values for three datasets; TCGA-LAML (P-value = 0.038, Figure 2C ), GSE12417(P-value = 0.025, Figure 2D ), GSE37642 (P-value = 0.0017, Figure 2E ). Moreover, the levels of MMs were significantly different between the two groups in the three datasets (P-value< 0.001) ( Figures 2F–H ). What’s more, we found that MMs constructed based on the expression of 31 MMRGs in the TCGA-LAML dataset could accurately diagnose the M0-M5 stages of AML patients (P-value< 0.05, Figure S2B , Table S4 ).

A total of 3,798 DEGs were identified (|logFC| > 0, adjusted P-value< 0.05) (3,036 upregulated and 762 downregulated genes) in the high-MMs group ( Figure 3A ). Furthermore, the expressions of the MMRGs were significantly different in the TCGA-LAML dataset ( Figure 3B ). To obtain genes closely associated with mitochondrial metabolism, WGCNA were conducted to identified modules in which genes were highly correlated. The modules with similarity lower than 0.25 ( Figure 3D ) were merged to obtain 11 modules at the scale independence of 0.9 ( Figure 3C ) (MEbrown, MEgrey60, MElightcyan, MEmagenta, MEpurple, MEcyan, MEmidnightblue, MEgreenyellow, MElightgreen, MElightyellow, and MEgrey). The correlations between the 11 module eigengenes and different MMs groups were then investigated ( Figure 3E ). The genes of the MEbrown (|r|=0.52, P-value =5e^-11), MEpurple (|r|=0.51, P-value =1e^-10), and MEcyan (|r|=0.55, P-value =2e^-12) modules with correlation coefficients greater than 0.5 were selected for the analysis. The 31 MMRGs were intersected with genes in MEbrown ( Figure 3F ), MEpurple ( Figure 3G ) and MEcyan ( Figure 3H ) modules to obtain the 17 MMRGs modules (BCS1L, COX10, DGUOK, ECHS1, ETHE1, MPV17, NDUFA1, NDUFA13, NDUFS2, NDUFS7, NDUFV1, POLG, SDHA, SUCLG1, SURF1, TACO1, and TIMM8A).

Additionally, the correlations between the 17 MMRGs modules in the TCGA-LAML cohorts were explored. The 17 MMRGs modules were significantly correlated according to the TCGA-LAML dataset ( Figure S3A ). BCS1L and NDUFV1 (r = 0.701, P-value< 0.001, Figure S3B ), MPV17 and NDUFV1 (r = 0.645, P-value< 0.001, Figure S3C ), DGUOK and NDUFA1 (r = 0.653, P-value< 0.001, Figure S3D ), NDUFA13 and NDUFS7 (r = 0.681, P-value< 0.001, Figure S3E ), NDUFS7 and NDUFV1 (r = 0.643, P-value< 0.001, Figure S3F ) were the top five pairs of genes with the highest correlations ( Figures S3B–F ).

Two AML subtypes (cluster1 and cluster2) were identified based on the expression of the MMRGs modules via an unsupervised clustering method in the TCGA-LAML to investigate differential expression of the MMRGs modules in AML patients ( Figure 4A ). Cluster1 and cluster2 contained 71 samples and 71 samples, respectively. Principal component analysis (PCA) showed that the two AML subtypes had significant differences ( Figure 4B ). The consensus clustering had the best stability when k = 2 ( Figures 4C, D ). Subsequent survival analysis for the two AML clusters indicated that the clinical outcomes of two AML subtypes were significantly different (P-value = 0.018, Figure 4E ). Furthermore, the expression levels of 17 MMRGs modules were significantly different between the two AML subtypes ( Figure 4F ).

Furthermore, the MMRGs modules were significantly upregulated in cluster1 in the TCGA-LAML dataset ( Figure S4A ). The expression of COX10 and TIMM8A was significantly different between the two AML subtypes (P-value< 0.01). However, the expression of the other 15 MMRGs was significantly different between the two AML subtypes (P-value< 0.001). The ROC curves showed that expression levels of COX10 (AUC = 0.650, Figure S4C ), NDUFA1 (AUC = 0.673, Figure S4H ), and TIMM8A (AUC = 0.630, Figure S4R ) had lower accuracy in distinguishing two AML subtypes, while the expression levels of BCS1L ( Figure S4B ), DGUOK ( Figure S4D ), ECHS1 ( Figure S4E ), ETHE1 ( Figure S4F ), MPV17 ( Figure S4G ), NDUFA13 ( Figure S4I ), NDUFS2 ( Figure S4J ), NDUFS7 ( Figure S4K ), NDUFV1 ( Figure S4 L), POLG ( Figure S4M ), SDHA ( Figure S4N ), SUCLG1 ( Figure S4O ), SURF1 ( Figure S4P ), and TACO1 ( Figure S4Q ) could accurately distinguish the two AML subtypes ( Figures S4B–R ).

First, nine MMRGs modules correlated with the survival outcomes with P-value< 0.1 were selected via univariate Cox regression analysis ( Figure 5A , Table S5 ) for further analysis. Finally, only five MMRGs modules (ECHS1, NDUFA1, NDUFS2, SDHA, SUCLG1) were selected through the LASSO regression analysis with minimized lambda ( Figures 5B, C ) for the construction of the prognosis model ( Figure 5D ).

Further validation of the prognostic value of the five MMRGs was conducted by analyzing the survival information of AML patients in TCGA-LAML dataset ( Table 1 ). A prognostic model was then constructed based on the expression of five MMRGs using multivariate Cox regression ( Table 2 ).

The prognosis model of the MMRGs was evaluated using nomogram analysis ( Figure 5E ). The expression of NDUFA1 and SDHA performed better significantly in the Cox regression model than other variables. In addition, the nomogram performed best at predicting 1-year OS based on calibration plots ( Figures 5F–H ). DCA at 1-year ( Figure 6A ), 2-year ( Figure 6B ), and 3-year ( Figure 6C ) were performed to assess the clinical applicability of the nomogram. The results indicated that the nomogram had a higher net clinical benefit at 1-year. The KM curve showed that low-risk patients had markedly longer survival times than other patients (P-value< 0.001, Figure 6D ). In addition, risk scores of two groups were significantly different (P-value< 0.001, Figure 6E ).

Further assessment was conducted to assess the predictive capability of prognosis model in the GSE12417 dataset ( Figures 6F, G ). Similarly, high-risk patients had worse outcomes (P-value = 0.045, Figure 6F ) and higher risk scores (P-value< 0.001, Figure 6G ) in the GSE12417 dataset, indicating that the prognosis model developed by MMRGs could accurately predict the survival outcomes of AML patients.

DEGs were obtained via differential analyses on the TCGA-LAML dataset. A total of 38 DEGs were identified (|logFC| > 1, adjusted P-value< 0.05) (33 upregulated and 5 downregulated genes) in the high-risk group displayed in the volcano plot and heatmap ( Figures 6H–I ).

Expressions of ECHS1 and NDUFS2 were analyzed by IHC staining in bone marrow (BM) biopsies of 10 AML patients, and 10 non-neoplastic patients. Control BM samples showed weak ECHS1 ( Figure 7A ) and NDUFS2 expression ( Figure 7B ), and most cells were negative. However, most cells showed strongly positive signals of ECHS1 and NDUFS2 in AML samples ( Figures 7A, B ). According to the statistics of mean IOD from IHC images, AML samples had significantly increased expression levels of ECHS1 ( Figure 7C ) and NDUFS2 ( Figure 7D ) when compared to normal BM samples.

To explore the biological process between two risk groups, GO and KEGG analyses were conducted based on the 38 DEGs. GO analysis showed that the DEGs were mainly involved in immune responses, including positive regulation of cytokine, neutrophil mediated immunity, T cell activation and antigen processing and presentation ( Figure 8 , Table S6 ). The GO enrichment pathways of the 38 DEGs were mainly concentrated in the BP pathway ( Figure 8E ). KEGG analysis showed that the DEGs were mainly enriched in four pathways, hematopoietic cell lineage, phagosome, antigen processing and presentation, and Th1 and Th2 cell differentiation ( Figure 8A , Figure S5 , Table S7 ). These results indicated that 38 DEGs may play a crucial role in regulating AML TME.

To further explore the correlation between the 38 MMRGs’ prognosis-correlated DEGs and disease, DO analysis was performed. The DEGs were significantly (P-value< 0.05) enriched in parasitic protozoa infectious disease, multiple sclerosis, hypersensitivity reaction type IV disease, demyelinating disease, vasculitis, psoriatic arthritis, lung disease, sarcoidosis, parasitic infectious disease, and malaria ( Figure S6A , Table S8 ). However, HMOX1 and HLA-DRB1 had the most correlations with diseases, and they were significantly correlated with seven diseases among the top ten diseases ( Figure S6B ).

To reveal the influence of MMRGs on AML occurrence, TCGA-LAML dataset was analyzed via GSEA to identify biological processes involved in two different groups. Surprisingly, these genes were remarkably associated with mitochondrial metabolism, drug resistance, and immunity, indicating that mitochondrial metabolism is crucial in immunoregulatory in TME and tumor drug resistance. Furthermore, genes highly in the high-risk group were significantly correlated with LSCs ( Figure S6D ), inflame pathway ( Figure S6E ), Ebola virus infection in host ( Figure S6F ), leukemia with MLL fusions ( Figure S6G ), mitochondrial translation ( Figure S6H ), glycolysis ( Figure S6I ), and multiple drug resistance ( Figure S6J ) ( Figures S6C–J , Table S9 ) in the TCGA-LAML cohort.

To further explore the molecular mechanisms underlying influence of mitochondrial metabolism on AML, A PPI network was built based on the 38 DEGs ( Figure 9A ). The top 10 DEGs were determined using Cytoscape’s cytoHubba plugin via MCC, MNC, and Degree algorithms ( Figures 9B–D ). Eight hub genes (CD14, CD74, HK3, HLA-DRB1, HLA-DRB5, LILRB2, S100A8, and S100A9) were identified with high connections ( Figure 9E ). Additionally, HLA-DRB1, S100A8 and S100A9 had the highest functional similarity with other hub genes ( Figures 9F, G ). The eight hub genes might play a pivotal role in immunoregulatory and treatment resistance.

To further explore the influence of mitochondrial metabolism on TME, the ESTIMATE, CIBERSORT, and ssGSEA algorithms were employed to characterize the TME. Results showed that the stromal scores were comparable between the two groups (P-value > 0.05, Figure 10A ). It was also observed that the high-risk group had higher immune scores (P-value = 0.007, Figure 10B ), and higher ESTIMATE scores (P-value = 0.018, Figure 10C ), but lower tumor purity compared with the low-risk group (P-value = 0.018, Figure 10G ), suggesting that high-risk patients had high immune cell infiltration. Interestingly, we found that the risk score showed a positive nonsignificant association with the stromal score (R = 0.091, P-value = 0.29, Figure 10D ), however, the risk score was significantly related to the immune score (R = 0.28, P-value = 0.00099, Figure 10E ), ESTIMATE score (R = 0.21, P-value = 0.012, Figure 11F ), and tumor purity (R = -0.2, P-value = 0.017, Figure 10J ).

Although the high-risk group showed higher immune cell infiltration, M2 macrophages, gamma delta T cells, myeloid-derived suppressor cells (MDSCs) and regulatory T cells (Tregs) were significantly higher in the high-risk group ( Figures S7 , S8 ). Additionally, inhibitory receptors (IRs), such as CTLA4 and IL2RA, were significantly overexpressed in high-risk group patients ( Figure S9 ). These indicated that patients in the high-risk group might have a strong immunosuppressive TME.

Given that immunotherapy has been considered as an important treatment for tumors, we also evaluate the immunotherapy response in AML patients in the TCGA-LAML dataset using the TIDE algorithm. The TIDE immunotherapy scores in high-risk patients were significantly higher compared with low-risk patients (P-value = 0.008, Figure 10H ), suggesting that responses to immunotherapy might be higher in low-risk patients. Further analysis revealed that the risk scores were significantly related to TIDE immunotherapy scores (R = 0.280, P-value< 0.001, Figure 10K ). High-risk patients showed significantly higher MMs (P-value< 0.001, Figure 10I ), and risk scores were significantly positively related to MMs (R = 0.556, P-value< 0.001, Figure 10L ).

To explore the role of hub genes in immune regulation, the correlations between immune infiltrating cells and hub genes were analyzed. By CIBERSORTx analysis, most immune cells were negatively correlated, with the highest positive correlation obtained between NK cells resting and Tregs in both groups, and the highest negative relationship was detected between monocytes and plasma cells ( Figures S7C-D ). Meanwhile, in both groups, 8 hub genes and 22 immune cells were found to be significantly correlated (P-value< 0.05) in both groups, and monocytes showed a significantly positive correlation with hub genes ( Figures S7D, E ). Through ssGSEA analysis, positive relationships were observed among the 15 immune cells in both groups ( Figures S8B, C ). Meanwhile, in the low-risk group of TCGA-LAML dataset, there were significant correlations (P< 0.05) between 15 immune cells and 8 hub genes, and most of them showed a significant positive correlation ( Figure S8D ). Specifically, CD14, HK3, LILRB2, S100A8, and S100A9 were significantly correlated with 15 immune cells. Similarly, significant positive correlations were found between 15 immune cells and 8 hub genes in high-risk group ( Figure S8E ), among which S100A8, S100A9, CD14, and LILRB2 were significantly correlated with 15 immune cells. Through correlation analysis, we found that the 8 hub genes were significantly positively correlated with the 11 immune-related genes ( Figure S9C ). These findings indicated that these hub genes were closely correlated with tumor immunity.

Given that tumor cells are characterized by genomic instability and altered metabolism, to determine the correlation between mitochondrial metabolism and DNA repair, we extracted 19 DNA repair genes from published literature. The analysis revealed that ERCC1, FANCC, FEN1, MGMT, MLH1, RAD23A, and XPC were upregulated in high-risk patients whereas MBD4 and WRN were downregulated in low-risk patients ( Figure S9B ). As displayed in Figure S9D , 8 hub genes were significantly related to 19 DNA repair genes, with WRN showing a negative correlation with the 8 hub genes.

To better illustrate the characteristic differences between the two groups, we compared clinical traits of the two groups. The proportion of patients under 60 years old in the low-risk group was significantly higher than that of high-risk group ( Figure S9E ), while there was no difference in sex between the high- and low-risk groups ( Figure S9F ). The proportion of OS events in low-risk group was markedly lower than that in the high-risk group ( Figure S9G ).

To investigate the regulatory mechanism of hub genes in AML and whether all of the hub genes are druggable targets, mRNA-RBP interaction analysis and mRNA-drug interaction analysis were conducted. An mRNA-RBP network with eight hub genes and 45 RBPs ( Figure 11A ) was constructed (69 pairs of mRNA-RBP interaction relationships) ( Table S10 ). Meanwhile, 39 potential drugs targeting six hub genes (CD14, CD74, HLA-DRB1, HLA-DRB5, S100A8, S100A9) were identified ( Figure 11B , Table S11 ). Besides, 31 drugs or molecular compounds targeting the HLA-DRB1 gene were identified.

What’s more, the drug sensitivity was predicted based on the mRNA expression profile and drug activity data of the eight hub genes in the GDSC, CCLE, and CellMiner database via pRRophetic algorithm. A total of 50 and 13 drugs interacted with eight hub genes in the GDSC and CCLE databases, respectively ( Figures 11C, D ). Fourteen drugs interacted with the other seven hub genes (except for the CD14 gene) in the CellMiner database ( Figure 11E ). These results indicate that the eight hub genes were all druggable targets. The majority of the hub genes are related to immune regulation, however small molecular inhibitors, such as BCL2 inhibitors, might target these genes. This suggested that treatment targeting mitochondrial metabolism might improve the TME of AML patients.

To evaluate the performance of the MMRGs prognosis model, univariate and multivariate Cox regression were conducted based on risk scores combined with 2 independent clinical variables (age, and gender). The results are shown in the forest plot ( Figure 12A ). Results indicated that the risk score remained an independent predictor of survival in the multivariate Cox regression analysis (P-value< 0.001, Table 3 ). A nomogram was developed which revealed that the risk score contributed the most risk points compared with other clinical variables ( Figure 12B ). Additionally, we conducted calibration analysis to evaluate the correctness of the model. The results indicated that the model could accurately predict AML patients’ OS at 1-, 2-, and 3-year ( Figures 12C–E ). The clinical benefit of this prognosis model at 1-, 2- and 3-year was also determined using the DCA ( Figures 12F–H ). This model displayed many clinical benefits over time.

The AML patients in the TCGA-LAML cohorts were grouped according to the OS to explore differences in risk scores between different OS event groups. Patients in the death group had higher risk scores ( Figures 12I ). Patients with lower risk scores had longer survival as determined from the KM curves ( Figure 12J ). Data shown in Figures 12K revealed that the prognosis model based on MMRGs risk score had high predictive accuracy (AUC1 = 0.808, AUC2 = 0.813, AUC3 = 0.806).