The expression landscape of ICDRGs

The results of the differential analysis showed that most ICDRGs such as ATG5, BAX, and CALR were up-regulated in AML samples (Fig. 1A), and some ICDRGs such as CD8A, CD8B, and IFNGR1 were down-regulated. There was no significant difference in the expression of ICDRGs such as ENTPD1, IFNA1, and IFNG between AML and normal samples. Univariate Cox regression analysis showed that most of the ICDRGs were risk factors for AML and were associated with poor prognosis of patients (Fig. 1B). Moreover, CALR, CASP1, CD4, CXCR3, FOXP3, IFNGR1, IL10, P2RX7, and PRF1 were significantly correlated with the prognosis of AML patients (p < 0.05) (Fig. 1C). We also found that almost all ICDRGs were positively correlated with each other except CALR, which was negatively correlated with multiple ICDRGs (Fig. 1B). Copy number variation analysis showed that the frequency of copy number gain of IL10, NLRP3, and CASP1 was significantly increased, while the frequency of copy number loss of MYD88, IL6, and IFNG was significantly increased (Fig. 1D). Somatic mutation analysis showed that the overall mutation rate of ICDRGs was low (Fig. 1E).

Fig. 1

Genetic characteristics of ICDRGs. A The interaction of ICDRGs in AML patients and its relationship with prognosis. B Differential expression analysis of ICDRGs between AML samples and normal samples. C Univariate Cox regression analysis was used to identify ICDRGs significantly associated with prognosis. D CNV frequency of ICDRGs in the TCGA-LAML cohort. E Somatic mutations of ICDRGs in the TCGA-LAML cohorts

Identification of ICD-related classifications in AML

To further explore the expression characteristics of ICDRGs in AML, consensus clustering algorithm was used to classify the meta-cohort. The results showed that when the number of clusters was 2, the classification effect was the best. 581 patients with AML were divided into two kinds of ICD-related molecular subtypes: Cluster A (n = 375) and Cluster B (n = 206) (Fig. 2A). The t-SNE algorithm showed significant differences between these two ICD-related molecular subtypes, which confirmed the reliability of the classification (Fig. 2B). The heatmap showed that most ICDRGs were upregulated in Cluster B compared with Cluster A (Fig. 2C). Kaplan–Meier curve analysis also showed that there was a significant difference in overall survival between the two subtypes, and the prognosis of patients in Cluster A was better (Fig. 2D). Thus, the two ICD-related molecular subtypes identified based on ICDRG expression are significantly different.

Fig. 2

Identification of ICD-related molecular subtypes and analysis of differences in TME characteristics between subtypes. A Two molecular subtypes were identified by consensus clustering. B The TSNE algorithm was used to verify the accuracy of molecular subtypes. C Heatmap shows the expression characteristics of ICDRGs between subtypes. D Survival analysis between subtypes. E–H Differences in TME scores (E), immune cell infiltration (F), immune checkpoint expression (G), and TIDE score (H) between subtypes

Analysis of TME differences between classifications

We further explored TME differences between the two ICD-related subtypes. TME analysis showed that Cluster B had higher TME scores, including stromal, immune, and ESTIMATE scores (Fig. 2E). The results of immune infiltration analysis showed that adaptive immune cells such as naive and memory B cells, plasma cells, CD8+ T cells, resting CD4+ T cells, regulatory T cells (Tregs), and innate immune cells such as resting NK cells and mast cells in Cluster A exhibited higher infiltration levels than those in Cluster B. Monocytes and M2 macrophages were mainly enriched in Cluster B (Fig. 2F). Therefore, we hypothesized that the better prognosis of Cluster A patients might be related to their higher proportion of anti-tumor immune cell infiltration. The expression levels of key immune checkpoints such as PD-1, PD-L1, and CTLA4 were not significantly different between the two groups, while HAVCR2, CD86, and TNFRSF9 were up-regulated in Cluster B (Fig. 2G). We further calculated the TIDE score to evaluate the immune escape ability of tumor cells. Strikingly, the TIDE score of Cluster A was significantly higher than that of Cluster B (Fig. 2H), suggesting that there may be some suppression of anti-tumor immune effects in Cluster A.

Constructing an ICD scoring system to explore the potential correlation between ICD and TME

Subsequently, we compared the differences in enrichment scores between tumor marker pathways between the subtypes, and the activity of these pathways was essentially higher in Cluster B, indicating that the subtype had more active signaling (Fig. 3A). In addition, we used the ssGSEA algorithm to determine the ICD score to evaluate ICD activity and explain the interaction between ICD and TME. First of all, the ICD score quantified ICD subtypes well, and the ICD score of Cluster B was significantly higher than that of Cluster A (Fig. 3B). The ROC curve also confirmed that the ICD score could effectively distinguish ICD subtypes (AUC = 0.900) (Fig. 3C). Second, high ICD scores were associated with lower infiltration of anti-tumor immune cells (Fig. 3D). Interestingly, the infiltrating proportion of Tregs was significantly negatively correlated with the ICD score and positively correlated with the TIDE score (Fig. 3E, F). Therefore, we hypothesized that infiltration of Tregs might be involved in the immunosuppression of Cluster A. Finally, the ICD score was also significantly positively correlated with the expression of immune checkpoints such as HAVCR2, CD86, TNFRSF9, TIGIT, IDO1, and LAG3 (Fig. 3G). TIDE also had a strong correlation with TME scores (Fig. 3H). These findings indicated that ICD could participate in regulating the TME in AML.

Fig. 3

Construction of ICD scoring system and correlation analysis with TME. A Differences in enrichment scores of tumor marker gene sets between ICD molecular subtypes. B Differences in ICD score between subtypes. C ROC curve analysis was used to verify the ability of the ICD score to discriminate subtypes. D Correlation between ICD score and immune cell infiltration. E, F Correlation of TIED score with Treg infiltration and ICD score. G Correlations of the immune checkpoint genes with ICD score. H Correlations of the ICD score with immune score, stromal score, and ESTIMATE score

Analysis of biological differences between classifications

To fully understand the potential biological functions of ICD-related subtypes in AML, we identified 444 DEGs associated with ICD-related subtypes through differential analysis. Volcano map shows that most DEGs are up-regulated in Cluster B (Fig. 4A). Then, we performed a functional enrichment analysis of DEGs. According to KEGG analysis, These ICD subtype-related DEGs are enriched in immune-related pathways such as cytokine-cytokine receptor interaction, NOD-like receptor signaling pathway, and cell adhesion molecules (CAMs) and antigen processing and presentation (Fig. 4B). The GO annotation showed that the main functions of these genes include positive regulation of cytokine production, immune response-regulating signaling pathway, and immune receptor activity (Fig. 4C). GSEA analysis showed that the activities of metabolism-related pathways such as fructose and mannose metabolism, galactose metabolism, glycolysis/gluconeogenesis, oxidative phosphorylation, and tryphan metabolism in cluster B were significantly higher (Fig. 4D). These results again suggest that ICD is significantly correlated with immune signaling, and that abnormal metabolic signaling may also contribute to poor prognosis in Cluster B patients by promoting malignant proliferation of tumor cells.

Fig. 4

Functional analysis of DEGs between ICD molecular subtypes. A The volcano plot shows the differential expression characteristics of genes between molecular subtypes. Red and blue dots indicate genes with significantly upregulated expression in Cluster B and Cluster A, respectively. B, C KEGG (B) analysis and GO annotation (C) of DEGs between ICD molecular subtypes. D GSEA analysis revealed significantly enriched signaling pathways in Cluster B

Prognostic predictive value of ICD subtype-related DEGs

Through univariate Cox regression analysis, we found that 34 ICD subtype-related DEGs were significantly associated with AML prognosis (p < 0.05). LASSO regression analysis was used to further reduce dimensionality and construct a risk score model with FGR, TFEB, KDM5B, SH3TC1, VNN1, TRIB1, HIP1, HTATIP2, AHR, CRIP1, THBS1, and IL1R2 as model factors (Fig. 5A, B) (Additional file 1: Table S2). AML patients were divided into high-risk and low-risk groups based on the best cut-off value (Fig. 5C). Compared with the low-risk group, the high-risk group had significantly more patients who died (Fig. 5D). Except for KDM5B, the expression of model genes was significantly higher in the high-risk group (Fig. 5E). Survival analysis showed that the overall survival of the high-risk score group was significantly shorter than that of the low-risk score group (Fig. 5F). ROC curve analysis confirmed the prognostic prediction efficacy of the risk score model, and the AUC values for predicting 1-, 3-, and 5-year overall survival were 0.698, 0.685, and 0.702, respectively (Fig. 5G). In the nine AML cohorts including the meta-cohort, it was confirmed that the prognosis of patients in the high-risk group was significantly worse (p < 0.05) (Fig. 6A, B). The multi-cohort data confirmed the prognostic predictive value of the risk score model, which was confirmed by ROC curve analysis (Fig. 6C). The TCGA-LAML cohort contained more clinical data. Univariate and multivariate Cox regression analysis confirmed that the risk score was an independent factor in predicting the prognosis of AML (p < 0.001) (Fig. 6D, E).

Fig. 5

Construction of risk scoring model. A The penalty coefficient of the minimum tenfold cross-validation error point was calculated to determine the corresponding model gene. B Determination of model gene coefficients. C–E Based on the optimal cut-off value, patients in the meta-cohort were divided into high- and low-risk score groups (C); the survival status distribution (D), and model gene expression (b) in high- and low-risk score groups. F Survival analysis between high- and low-risk score groups. G time-dependent ROC curve analysis of risk scores in the meta cohort

Fig. 6

Validation of risk score model. A, B Survival analysis between high- and low-risk groups in the validation cohorts (A) and the constituent cohorts of the meta-cohort (B). C time-dependent ROC curve analysis of risk scores in the meta-cohort and the validation cohorts. D, E Univariate and multivariate Cox regression analysis of clinicopathological factors and risk score in the TCGA-LAML cohort

Potential molecular mechanisms affecting the prognosis of patients in different risk groups

To better reveal the factors affecting the prognosis of patients in different risk groups, we systematically evaluated the differences in TME characteristics and clinicopathological factors between the two groups. The alluvial diagram shows that almost all patients in the low-risk group belong to cluster A, while patients in cluster B belong to the high-risk group, indicating that the risk score model further groups AML patients according to clinical outcomes based on ICD related molecular subtypes (Fig. 7A). The ICD score of patients in the high-risk group was significantly higher than that in the low-risk group, and the risk score was also significantly positively correlated with the ICD score, indicating that the risk score can also accurately reflect the ICD characteristics (Fig. 7B, C). The analysis of differences in immune cell infiltration and checkpoint expression showed that the difference between high-risk and low-risk groups was consistent with the difference between ICD-related subtypes, that is, the low-risk group showed increased infiltration of innate and adaptive immune cells, including Tregs, while the high-risk group showed up-regulated expression of immune checkpoints such as HAVCR2, LAG3, CD86 and TNFRSF (Fig. 7D, E). The low-risk group similarly showed a higher TIDE score (Fig. 7F). In addition, the activity of most tumor marker gene sets in Cluster A was higher than that in Cluster B (Fig. 7G). Combining these analysis results, we speculate that although the prognosis of low-risk patients is better, the infiltration of Tregs may hinder the immune effect of anti-tumor immune cells, and the high expression of immune checkpoints and active cancer-promoting signaling pathways may be the reasons for the poor prognosis of high-risk group.

Fig. 7

Differences in TME characteristics between high- and low-risk score groups. A Alluvial plots show the distribution of molecular subtypes, risk-score groups, and vital status of patients. B, C Differences in ICD scores between high- and low-risk groups and the correlation between risk scores and ICD scores. D–G Differences in immune cell infiltration (D), immune checkpoint expression (E), TIDE score (F), and tumor marker pathway enrichment score (G) between high- and low-risk groups

Differences in genomic traits and clinicopathological features between low- and high-risk groups

We further analyzed the differences in clinicopathological characteristics between risk groups in the TCGA-LAML cohort. First, compared with the low-risk group, patients in the high-risk group were older (aged > 60 years), had poor cytogenetic risk, had more FAB classifications of M0, M4, and M5, and had a higher proportion of platelets and white blood cells (WBC) (Fig. 8A). Correspondingly, these patients also had higher risk scores (Fig. 8B). In addition, there was no significant difference in the proportion of patients with common somatic mutations between high-risk and low-risk groups (Fig. 8C), but the risk score of patients with negative FLT3 and NPM1 mutations was significantly lower than that of patients with positive (Fig. 8D). According to somatic mutation analysis, in the two cohorts (TCGA-LAML and Beat AML), FLT3, DNMT3A, NPM1, TP53, and RUNX1 mutations were most likely to occur in the high-risk group, and the mutated genes in the low-risk group mainly included FLT3, NPM1, CEBPA, TET2, IDH2, and WT1 (Fig. 8E). These results suggest that the differences in clinicopathological characteristics and somatic mutations may be important reasons for the differences in AML prognosis.

Fig. 8

Differences in clinical characteristics and somatic mutations between high- and low-risk groups. A, B The differences in the proportion of conventional clinicopathological factors between high- and low-risk groups and the differences in risk scores between different clinicopathological factors were compared. C, D The difference in the proportion of somatic mutation positive and negative patients between high- and low-risk groups and the difference in risk scores between somatic mutation positive and negative patients. E Differences in overall somatic mutation frequency between high- and low-risk groups. FAB: French–American–British; WBC: white blood cell

Differences in chemosensitivity and response to immunotherapy between risk groups

We predicted sensitivity to commonly used drugs in AML, with the lower-risk group having a lower IC50 to doxorubicin, cytarabine, and midostaurin and being more sensitive to them (Fig. 9A, C). Analysis of drug sensitivity data from ex vivo leukemia cells in the Beat AML cohort showed that the low-risk group was more sensitive to treatment with GSK-1838705A and Venetoclax. While the high-risk group was more sensitive to Elesclomol, Flavopiridol, MK-2206, Nilotinib, Panobinostat, and Selumetinib (AZD6244) (Fig. 9D). Moreover, in both Beat AML and GSE14468 cohorts, we observed that the proportion of patients who responded to induction chemotherapy was higher in the low-risk group than in the high-risk group, and that patients who did not respond to induction chemotherapy had significantly higher risk scores than those who did (Fig. 9E, F). We used the TIDE algorithm to predict the possibility of patients in different risk groups responding to immunotherapy, and the results showed that the high-risk group (40%, 158/396) was more likely to respond to immunotherapy than the low-risk group (19%, 35/185) (P = 4.426e−07) (Fig. 9G). We also used subclass mapping to compare the expression profiles of high- and low-risk groups with another dataset containing 47 melanoma patients who responded to immunotherapy [19]. Excitingly, in the meta-cohort, we observed that the high-risk group was more likely to respond to anti-PD-1 (P = 0.001) and anti-CTLA4 (P = 0.013) treatments (Fig. 9H). This may be related to the higher expression level of immune checkpoints in the high-risk group.

Fig. 9

Differences in chemotherapy sensitivity and immunotherapy response between high- and low-risk groups. A–C Sensitivity prediction of cytarabine, doxorubicin, and midostaurin for AML in high- and low-risk groups. D Chemotherapeutic agents with differential sensitivity between risk groups in the Beat AML cohort. E, F The differences in risk scores of patients with or without response to induction chemotherapy and the proportion differences between risk groups. G Prediction of the proportion of patients with or without response to immunotherapy in different risk score groups. H Prediction of response to anti-PD-1 and anti-CTAL4 immunotherapy in different risk score groups

The expression of risk score model genes was validated in a clinical real-world cohort

We collected five normal samples, five myeloid leukemia chronic-phase samples, and five myeloid leukemia acute-phase samples for transcriptome sequencing. Compared with the normal samples, the expression of KDM5B and IL1R2 was up-regulated in both chronic and acute phase samples, and the expression of TRIB1, HIP1, VNN1, and IL1R2 was only up-regulated in the chronic samples, The expression of FGR was up-regulated in chronic phase samples and down-regulated in acute phase samples, and the expression of TFEB, SH3TC1, and AHR was down-regulated in both chronic and acute phase samples. Moreover, there was no significant change in the expression of THBS1, CRIP1, and HTATIP1 (Fig. 10).

Fig. 10

The clinical cohort was used to validate the expression of risk score model genes. Differences in the expression of risk-scoring model genes among normal samples, myeloid leukemia chronic-phase samples, and myeloid leukemia acute-phase samples