Risk stratification of acute myeloid leukemia: Assessment using a novel prediction model based on ferroptosis-immune related genes

1.   Introduction

Traditional treatment strategies for acute myeloid leukemia (AML) include intensive chemotherapy, hypomethylating medications, and hematopoietic stem cell transplantation. The 5-year overall survival (OS) rate remains poor, at no more than 30% ^[1]. Strategies for AML treatment have recently evolved away from traditional therapies toward immunotherapy and targeted therapy. However, more than half of the patients still experience relapse or become treatment refractory. The median OS for these patients is only 0.5 years ^[2]. It is therefore critical to find a more accurate prognostic risk stratification that can assist in accurately identifying patients who will respond to specific treatments.

Ferroptosis is a modulated form of cell death involving iron overload and the overproduction of reactive oxygen species ^[3]. Activation of ferroptosis triggers the death of oncocytes, especially in cancers resistant to traditional treatment ^[4,5,6]. Although the precise effects and mechanisms of ferroptosis in the pathogenesis of AML remain unclear, some studies have shown that patients with AML frequently exhibit abnormal expression of ferroptosis-associated genes. Promoting ferroptosis in AML cells through treatment with the multikinase inhibitor sorafenib could help in reversing the drug resistance of tumor cells ^[7,8]. This suggests that ferroptosis could be used as a new indicator for AML risk stratification.

Immunotherapy, including antibody-based treatment, tumor vaccines, immune checkpoint inhibitors, and adoptive cellular immunotherapy, has shown great clinical application prospects for AML ^[9]. However, research on immunotherapies for AML has achieved fewer breakthroughs, compared to those for plasma cell neoplasms and lymphoid malignancies. Studies have demonstrated interactions between ferroptosis and immune regulation in the AML microenvironment ^[10,11]. Preclinical studies have shown that some ferroptosis inducers may strengthen the clinical therapeutic effect of immune checkpoint inhibitors in tumors that tend to respond poorly to immunotherapy ^[12]. Immune checkpoint inhibitors and ferroptosis activators synergistically inhibit hepatocellular carcinoma cell proliferation. However, the mechanism underlying this has not yet been fully elucidated ^[13,14]. As a result, screening and validating the involvement of ferroptosis–immune related genes (FIRGs) in AML prognosis and treatment sensitivity is important.

The present study aimed to identify differential FIRGs in patients with AML. A risk score model on the foundation of five FIRGs was constructed to forecast the AML prognosis. Finally, the potential therapeutic value of the risk assessment model was determined. Our findings could have clinical implications for the prognosis and treatment of AML.

2.   Materials and methods

2.1.   Data collection from public databases

The transcriptomic data of 129 samples from patients with AML with both clinical and information on survival were acquired from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). TCGA-LAML cohort was used as the training dataset. Two validation datasets, both generated on the GPL96 platform, were acquired from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The transcriptome data and survival information for 162 and 417 AML samples were obtained from GSE12417 and GSE37642, respectively. The GSE164894 data set was used to compare the risk score between AML patients who did and did not achieve complete remission (CR). The GSE60926 data set was further used to compare the risk scores between acute lymphoblastic leukemia patients with and without central nervous system (CNS) involvement. Principal component analysis (PCA) analysis was employed during exploratory data analysis (EDA) ^[15].

2.2.   Assessment of ferroptosis-related enrichment score

FIRGs were retrieved from the MsigDB data base (http://www.gsea-msigdb.org/gsea/msigdb/) and FerrDb data base (http://www.zhounan.org/ferrdb/). After merging and de-duplication, 275 FIRGs were identified. The ferroptosis-related enrichment score (ferroptosis score) for each AML sample in the TCGA-LAML dataset was calculated using single sample gene set enrichment analysis (ssGSEA) ^[16,17]. The optimal cutpoint was obtained with the surv_cutpoint function from the R package 'survminer' (v0.4.9) based on the ferroptosis score of each sample ^[18]. The 129 AML patients from TCGA database were then separated into high- and low-ferroptosis score groups based on the cutpoint. To compare the differences in survival of the high- and low-ferroptosis scoring groups, Kaplan-Meier (K-M) analysis and log-rank testing were used ^[19]. Furthermore, to validate the relationship between the ferroptosis enrichment score and ferroptosis, differential expression of ferroptosis signal pathway marker genes was analyzed in the two groups. Genes involved in the ferroptosis pathway were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/; map04216).

2.3.   Identification of FIRGs

The 'limma' R package was used to identify differentially expressed genes (DEGs) between two ferroptosis score groups in the TCGA dataset ^[20]. The cutpoint for DEGs was P-value < 0.05 and log2 fold change |(FC)| > 1. To visualize the variations in genetic expression levels between the two groups, a Volcano plot was created. The expression of the top 50 DEGs was depicted in a heat map. The ImmPort database (https://www.immport.org/home) yielded a total of 1793 genes related to immunity. A Venn diagram between DEGs and 1793 immune related genes were used to identify FIRGs.

2.4.   Functional and pathway enrichment analyses of FIRGs

The 'clusterProfiler' program was adopted to perform gene ontology (GO) enrichment as well as KEGG pathway analyses on the FIRGs to explore the common functions and pathways ^[21]. Biological processes (BP), cellular components (CC), and molecular functions (MF) linked to FIRGs were mainly defined using GO enrichment analysis. The biological pathways associated with FIRGs were discovered using KEGG pathway analysis. The P-value threshold for enrichment significance was adjusted to 0.05.

2.5.   Construction and validation of the prognostic risk score model

To find FIRGs that were substantially linked with the OS of AML patients in the training sample, univariable Cox regression analyses were used. To detect FIRG signatures and determine their respective coefficients value, LASSO regressive analyses of candidate FIRGs were performed ^[22,23]. Using the 'glmnet' R package ^[24], a tenfold cross verification was utilized to find optimal values for the alpha and lambda parameters. The following equation was introduced to calculate the risk score: Risk score = β1X1 + β2X2+…+ βnXn, where β refers to the regressive coefficient and X1 denotes the genetic expression level. The median risk score was computed via the 'survminer' program in R program to separate AML patients into high-risk and low-risk groups. The K-M analysis and log-rank test were then employed to contrast the OS of the two groups. The K-M 'survival ROC' program in R software was employed to plot receiver operating characteristic (ROC) curves and measure the area under the ROC curve (AUC) values to predictive efficiency of the risk model. Better efficiency was indicated by larger AUC value ^[25]. In addition, a risk plot was created with R software's 'pheatmap' package ^[26]. The effectiveness of the risk scoring model in distinguishing the high-risk group from the low-risk group was further assessed via t-distributed stochastic neighbor embedding (t-SNE) analyses ^[27]. This approach was used to validate both GSE12417 and GSE37642 datasets. Moreover, integrated discrimination improvement (IDI) analysis was performed to compare the overall improvement of the risk prediction models ^[28,29,30].

2.6.   Independent prognostic analysis and construction of a nomogram

Univariable regression and multivariable Cox regression analyses were performed to find independent prediction factors of OS in the TCGA-LAML dataset. Using the R packages 'rms' and 'nomogramEx' ^[11], a nomogram including the risk scoring and clinicopathological parameters was created. Corresponding and correction curves of 1-, 3-, and 5-year OS were constructed with the 'regplot' R package to further verify the accuracy of this nomogram ^[31]. Decision curve analysis (DCA) curves at 3 and 5 years were also drawn to evaluate the prognostic model ^[32].

2.7.   Immune infiltration analysis by ssGSEA

The ssGSEA was introduced to determine the enrichment score of 24 immunocytes ^[33]. A box plot was drawn to visualize the differences in infiltration of 24 immune cell types between the low- and high-risk groups. Associations between the prognostic model genes and differential immune cells were analyzed using Spearman's correlation analysis. Correlation coefficient |cor| > 0.3 and P-value <  0.05 were considered statistically significant.

2.8.   Assessment of the immunotherapeutic and chemotherapeutic responses

Based on the half-maximal inhibitory concentration (IC50) reported for every patient on the Genomics of Drug Sensitivity in Cancer web site, the 'pRRophetic' R package was employed to study the chemotherapeutic response ^[34]. In addition, an unsupervised subclass mapping method (SubMap) was applied to predict response to anti-programmed cell death ligand 1 (PD1) and anti-cytotoxic T-lymphocyte associated protein 4 (CTLA4) immunotherapy. A data file containing 47 cutaneous melanoma samples and expression data of 791 genes was obtained from a prior study ^[35]. The R package 'pheatmap' was utilized to visualize these findings.

2.9.   Statistical analyses

The R programming language was used for all studies. The function and references for R packages are listed in Table S1. Data from multiple groups were compared employing the t-test, Wilcoxon rank-sum test, Kruskal-Wallis H test, and Chi-square test. If not specified above, a P-value less than 0.05 was considered statistically significant.

3.   Results

3.1.   Establishment and validation of ferroptosis score

The transcriptome data and clinical features of the patients with AML were acquired from the TCGA and GEO databases. The TCGA-LAML cohort was chosen as the training dataset. GSE12417 and GSE37642 were the validation datasets. Table S2 provides details on clinical characteristics of the training dataset and two validation datasets. Exploratory data analysis of the TCGA-LAML cohort revealed no major changes in the global gene expression of AML patients with various survival statuses or durations (Figure S1).

The workflow for the study is shown in Figure 1. The study included 129 AML patients identified from the TCGA database. The ferroptosis score of every sample was obtained by ssGSEA. Using the median ferroptosis score as the threshold, all AML cases were separated into a high-ferroptosis score group and a low-ferroptosis score group. The OS of the low-ferroptosis group was superior to that of the high-ferroptosis group (median OS, 66 vs. 25.3 months, P < 0.05; Figure 2A). Furthermore, the differences in gene expression of biomarkers involved in ferroptosis signaling between two groups were assessed. The expression of 11 genes (ATG3, ATG5, GCLC, GCLM, GPX4, LPCAT3, PCBP1, PRNP, SAT1, SLC3A2, and SLC39A8) was substantially different between the two groups. These findings suggested a significant association between the predicted ferroptosis score and ferroptosis in patients with AML (Figure S2).

3.2.   Identification of FIRGs and functional enrichment analyses

Limma software was adopted to examine genetic expression changes between two groups. The 1061 DEGs in the groups (Figure 2B, Figure S3, Table S3) included 104 downregulated DEGs and 957 upregulated DEGs. Using the intersection of DEGs and immunity-related genes, 157 FIRGs were discovered (Figure 2C, Table S4).

The FIRGs were examined by GO and KEGG enrichment analysis. GO analysis of BP revealed that FIRGs were most enriched in immune response, leukocyte chemotaxis, and phagocytosis. In terms of CC, the FIRGs were highly enhanced on the outer surface of the plasma membrane. Regarding MF, FIRGs were enriched in signaling receptor activator and cytokine and chemokine activity, in addition to receptor binding. The top ten enrichment results for each GO category are illustrated in Figure 2D. KEGG pathway enrichment results revealed that the FIRGs were remarkably enriched in cytokine-cytokine receptor interplay, chemotactic factor signal path, neutrophil extracellular trap formation, phagosome, and signaling pathways for tumor necrosis factor (TNF), interleukin (IL)-17, mitogen-activated protein kinase, nuclear factor-kappa B (NF-κB), and Toll-like receptor (TLR) (Figure 2E).

3.3.   Establishment and verification of five-gene prognosis risk score model based on FIRGs

The prognostic significance of the FIRGs was evaluated using univariable Cox regression analyses on the training set. Thirty-three FIRGs were significantly associated with the OS of AML patients (P < 0.05; Figure 3A). To construct a prognosis model, these 33 FIRGs were further included in a LASSO regression analysis. Five FIRGs (BMP2, CCL3, EBI3, ELANE, and S100A6) were screened (Figure 3B, C). The risk scoring was calculated according to the expression levels and regression coefficients as: risk score = 0.01394 × BMP2 + 0.07846 × CCL3 + 0.11465 × EBI3 + (−0.0783) × ELANE + 0.16425 × S100A6. Based on the average risk score, AML patients were divided into high- and low-risk groups. Based on the five-gene prognostic model, Figure 4J shows risk score distribution, OS results, and model gene expressions in the two risk groups. Patients in the high-risk group had lower survival durations than patients in the low-risk group. Reduced expression of ELANE and increased expression of the other four genes were linked to higher risk. As indicated in Figure 4A, high-risk patients had a much lower OS in contrast to low-risk patients (median OS, 27.9 vs. 116.8 months, P < 0.001). Time-dependent ROC curves were constructed to illustrate predictive performance. In the training set, the 1-, 3-, and 5-year AUC was 0.698, 0.763, and 0.776, respectively, indicating good efficiency (Figure 4D). Meanwhile, t-SNE analysis confirmed that the patients in the two risk groups had been allocated (Figure 4G). The collective findings demonstrated that the constructed model performed well in predicting OS of patients with AML.

We computed the risk scores of individuals in the GSE12417 and GSE37642 datasets to verify the dependability of the constructed model. The outcomes of the two validation sets' analyses were consistent with those of the TCGA cohort. The high-risk group in GSE12417 also exhibited a low survival rate (Figure 4B, K), with an AUC > 0.6 in the ROC curve at 1, 3, and 5 years (Figure 4E). The t-SNE analysis illustrated that the patient distributions in these two groups were remarkably diverse (Figure 4H). The risk model validation results from dataset GSE37642 were consistent with those from GSE12417 (Figure 4C, F, I, and L). The results of the analyses of the two validation sets matched the results from the TCGA cohort. Furthermore, we compared our risk score model to risk models utilized in prior AML studies (PMID: 35111195; PMID: 35359394), as identified by the IDI. The discrimination of our model was higher than that of two recently reported models (the model in PMID 35111195: IDI 10.2%; 95% CI: 3.7 to 16.5%; P = 0.018; the model in PMID 35359394: IDI 10.2%; 95% CI: 3.2 to 16.1%; P = 0.016). The results indicated that our risk model significantly improved predictive accuracy (Figure S4). The five-gene prognostic model could successfully forecast the OS of patients with AML.

3.4.   Association between risk score model and clinicopathological features

We explored the associated clinical data of patients with AML in TCGA to determine the association between the model and clinicopathology traits. The analysis revealed that the risk scoring differed significantly between age groups, FAB subtypes, and OS (Table 1). Stratification analysis revealed the prognosis signature as a promising predictor of survival of the AML clinical subtypes. In addition, the CR rate in the high- and low-risk patients was explored. Patients who achieved CR appeared to have lower risk score compared to the patients without CR, but there was no statistical significance between the two groups (Figure S5A). Low-risk patients seemed to achieve non-significantly higher CR rates than high-risk patients (Figure S5B). We then investigated the relationship between CNS involvement and risk score in acute lymphoblastic leukemia patients. As shown in Figure S6A, these patients displayed significantly higher risk score than patients without CNS leukemia. Patients with high-risk had a higher proportion of CNS leukemia (Figure S6B). The result has a potential reference value for exploring the rate of CNS leukemia in different risk score groups of AML patients.

3.5.   Risk factor as an independent prognostic factor and nomogram construction

Univariable and multivariable Cox regression analyses were performed to explore if risk scoring could be considered a prognostic factor for individuals with AML independent of other clinical variables. Age, ferroptosis group, and risk score were all linked to patient prognosis in the univariate analysis (Figure 5A). The risk score might be utilized to predict patient prognosis independently, according to the multivariate Cox regression analysis (Figure 5B). In multivariate analysis, age was also a significant prognostic factor.

The independent prognostic factors of age and risk score were used to create a nomogram to forecast 1-, 3-, and 5-year OS (Figure 5C). The nomogram's C-index was 0.680. Calibration plots revealed that the nomogram succeeded in predicting the three survival probabilities in AML patients (Figure 5D). DCA of the risk model and nomogram was then performed to assess the prediction effectiveness (Figure S7A, B).

3.6.   Comparison of immune infiltration between two risk groups

Next, ssGSEA was used to calculate the percentage of 24 immune cells in the two groups. Low-risk patients had more T helper and Th2 cells than high-risk patients (Figure 6A). In comparison to the low-risk group, the high-risk group had considerably more activated dendritic cells (aDCs), CD8⁺ T cells, cytotoxic T cells, induced DCs (iDCs), neutrophils, macrophages, natural killer (NK) CD56dim cells, NK CD56bright cells, Th17 cells, Th1 cells, T follicular helper cells (TFHs), and regulatory T cells (Treg). These results suggested that the immunological milieu of the two risk groups was significantly different. The Pearson correlation between different immune cells and prognostic genes was then examined (Figure 6B). S100A6 was found to have a strong relationship with iDCs, macrophages, and neutrophils. T helper cell levels of chemokine (C-C motif) ligand 3 (CCL3) and S100A6 exhibited a significant negative correlation. TFH, Th1 cells, macrophages, Tgd, neutrophils, NK CD56bright cells, and iDCs were also correlated with bone morphogenetic protein 2 (BMP2). CCL3 was linked to macrophages, iDC, neutrophils, NK CD56bright cells, NK CD56dim cells, Th17 cells, Th1 cells, and TFH cells. A substantial positive connection was found between Epstein-Barr virus-induced gene 3 (EBI3) and aDCs, NK CD56dim cells, TFHs, and macrophages.

3.7.   The five-gene prognostic risk model predicts the sensitivity to chemotherapy and immunotherapy

The SubMap algorithm was utilized to predict the clinical activity of immune checkpoint inhibitors in the two risk groups (Figure 6C). The high-risk group was more responsive to anti-PD-1 immunotherapy (P < 0.05). In addition, using drug sensitivity data obtained from the Genomics of Drug Sensitivity in Cancer website, we estimated the chemosensitivity of the two groups to 133 drugs. High-risk patients had larger IC50 values for midostaurin (a multitargeted kinase inhibitor), NSC 87877 (inhibitor of Src homology-2 domain-containing protein tyrosine phosphatase-2), and pyrimethamine (inhibitor of dihydrofolate reductase) (Figure 6D). The findings indicate the value of the five-gene prognostic risk model in predicting treatment sensitivity to immunotherapy and targeted therapy.

4.   Discussion

An accurate prognostic signature in patients with AML might help guide individualized treatment plans and improve therapeutic efficacy ^[36]. Ferroptosis is an iron-dependent form of cell death driven by the increase of lipid peroxidation to fatal levels. There is growing interest regarding the interaction between ferroptosis and AML treatments ^[37]. Furthermore, ferroptosis affects the tumor immune microenvironment and potentially influences the outcome of immunotherapy ^[38,39]. Comprehensive research of the role of FIRGs in AML prognosis may reveal new avenues in the prediction of patient prognosis.

We used ssGSEA to determine the ferroptosis-related enrichment scores for each TCGA-LAML sample with a collection of known ferroptosis-related genes. The samples were subsequently separated into high- and low-ferroptosis score groups, with significant DEGs identified in the two groups. This method may provide theoretical support for the discovery of novel ferroptosis-related genes. Intersecting the DEGs with immune related genes identified 157 FIRGs. Functional enrichment analysis revealed that 157 FIRGs were enriched in the TNF, phosphoinositide 3-kinase (PI3K)-Akt, and IL-17 signaling pathways. These three pathways are associated with AML pathogenesis. In AML patients, abnormally elevated TNF-α levels could promote Treg cell development through the TNF-α–TNFR2 pathway, thus aggravating the progression of AML ^[40]. Furthermore, abnormal activation of the PI3K-Akt or IL-17 pathway is also linked with discouraging prognosis in AML ^[41,42].

We analyzed prognosis-related FIRG signatures in patients with AML. Through univariate Cox and LASSO regression analyses, five signature genes (BMP2, CCL3, EBI3, ELANE, and S100A6) were screened for the prognostic risk score model. Among them, BMP2, CCL3, and EBI3 are related to the progression and recurrence of AML. For example, previous studies implicated BMP2 in the promotion of stemness maintenance and self-renewal in leukemia stem cells, causing drug resistance and recurrence of AML ^[43,44]. CCL3 is a poor prognostic factor for AML, and elevated CCL3 in the leukemic environment suppresses erythropoiesis ^[45]. EBI3 participates in the formation of IL-35, which can contribute to an immunosuppressive AML microenvironment and facilitate immune escape of leukemic cells ^[46,47]. In addition, ELANE and S100A6 may also be associated with other hematopoietic malignancies ^[48]. Their roles in AML deserves further study.

The present study discovered that individuals with high-risk scores have poor prognosis, which is consistent with the prognostic effect of models based only on ferroptosis-related genes ^[29]. External validation confirmed the efficacy of this risk model in predicting the survival of AML patients. Importantly, the model was verified as an independent factor for OS in individuals with AML in both univariable and multivariable Cox regression analyses. To forecast the survival of individuals with AML, we developed a nomogram and assessed the validity of independent prognostic factors. Despite these accomplishments, further studies are required to verify the prognostic importance of the ferroptosis-immune related risk score.

We also investigated the association between the risk score model and immunocyte infiltration as well as therapeutic sensitivity to understand more about the implications of the model for AML therapy. The high-risk score group showed higher infiltration of Treg, Th17, and CD8⁺ T cells. Tregs facilitate immune escape of AML, further promoting disease progression and relapse ^[49,50]. Treg cells inhibit effector T cells through direct and indirect mechanisms, including the secretion of immunosuppressive cytokines, increased ATP hydrolysis, and increased adenosine production ^[51,52,53]. The levels of Th17 cells and their related cytokines (IL-17, IL-1 β, and IL-6) are increased in newly diagnosed and incomplete remission patients with AML, suggesting that Th17 cells may be involved in the pathogenesis of AML ^[54]. This was confirmed by the significant enrichment of FIRGs in the IL-17 pathway, as discussed previously. Notably, we presently observed that the high-risk group had a higher CD8⁺ T cells infiltration than the low-risk group. Although CD8⁺ T cells are crucial effector cells involved in the killing of cancerous cells, their function may be suppressed by AML blasts ^[55]. The sensitivity of different risk score groups to chemotherapy and immunotherapy was also investigated. Anti-PD-1 treatment produced a higher response rate in the high-risk group. PD-1/PD-L1 inhibitors have been previously demonstrated to fully restore the functional activity of pre-existing intratumoral tumor-infiltrating lymphocytes ^[55,56]. Therefore, we speculated that PD-1 therapy may be important in the recovery of immune function of immunosuppressed CD8⁺ T cells in patients with AML with a high-risk score. The low-risk group appeared to be more sensitive to midostaurin, NSC 87877, and pyrimethamine treatments. Our findings suggest that a risk score model based FIRGs is able to forecast therapeutic sensitivity and could help clinicians make more personalized treatment decisions.

5.   Conclusions

A prognostic risk score model was constructed and validated based on five ferroptosis–immune associated genes to accurately predict OS of AML patients. This risk score model could provide a basis for optimizing individualized clinical treatments.

Acknowledgments

We thank Dr. Yu-chen Li (Nanjing Agricultural University) and Dr. Sen-lin Ji (Affiliated Drum Tower Hospital, Nanjing University) for the help in our data analysis.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.