RNA-seq data and clinical information of 256 P-AML samples were collected at the first presentation/diagnosis from the Therapeutically Applicable Research to Generate Effective Treatment (TARGET) AML program and were available from NCI’s data portal (collected on 18/September/2021). Out of these 256 samples, 145 have been harmonized by the NCI’s Genomic Data Commons (GDC) (Zhang et al., 2021). Raw expression counts of these 145 AML patients were extracted TARGET-AML program of GDC (collected on 18/August/2021), as the TARGET-discovery cohort. The remaining 111 cases were used as a validation cohort and raw expression counts were gathered from the (NCI)’s data portal (TARGET-validation cohort). For model validation in the adult AML cohort, RNA-Sequencing data of 151 adult patients with AML and the corresponding clinical information were extracted from the The Cancer Genome Atlas (TCGA)-LAML program of GDC. The “DESeq2” package based on the negative binomial distribution was used to normalize the raw count expression data (Love et al., 2014). Detailed clinical information about the data cohorts is provided in Table 1, and the workflow is briefly depicted in Figure 1.

We obtained 259 ferroptosis-related genes from the FerrDb database (http://www.zhounan.org/ferrdb, collected on 16/November/2021), the first manually curated resource for regulators and markers of ferroptosis, which was released in January 2020 (Zhou and Bao, 2020) (Supplementary Table S1). LncRNAs were extracted by the GENCODE v20 annotation (http://www.gencodegenes.org) (Derrien et al., 2012). Ferroptosis-related lncRNAs co-expressed with ferroptosis-related genes were identified according to Pearson’s correlation analysis by the correlation test function of R (correlation coefficient Cor >0.8, p < 0.001).

Univariate Cox regression was applied to screen prognosis-related mRNAs and lncRNAs in the discovery cohort, using p < 0.05 as the cutoff. Then, model-making procedures were conducted for three lists of prognostic ferroptosis-related signatures separately, ferroptosis-related mRNA genes (FG), ferroptosis-related lncRNAs (FL), and ferroptosis-related mRNA genes combined with ferroptosis-related lncRNAs (FGL), and consisted of the following steps. Step 1, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was conducted on the prognosis-related signatures, and the optimal penalty parameter “λ” was selected by the 10-fold cross-validation method. Step 2, the prognostic risk score formula was established as follows: risk score = expression of gene 1×Coef1+expression of gene 2 × Coef 2+ ⋯ +expression of gene n × Coefn. Coef indicates the regression coefficients of each signature selected from the LASSO regression analysis. According to the risk model, samples in the discovery cohort were given a risk score and then divided into high- and low-risk groups using the median score as the cutoff. The survival curves of the patients in the high- and low-risk groups were drawn with the R-package “Survival,” and the survival time of the two groups was compared by log-rank test. The receiver operating characteristic (ROC) curves were drawn using the R package “survival ROC” for validation of the risk model and the AUC values of 1-, 3-, and 5-year survival were calculated. The best prognostic model identified with prognostic significance was selected for further analysis and the same algorithm was performed in the TARGET-validation cohort, TARGET-combined cohort, and TCGA-LAML cohort (151 adult patients), with the same coefficients derived for the discovery dataset.

Tumor microenvironment, immune, and stroma scores were calculated based on the gene expression data using ESTIMATE (Yoshihara et al., 2013) and xCell (https://xcell.ucsf.edu) tool (Aran et al., 2017). Estimating the proportion of immune and cancer cells (EPIC) was applied to estimate the infiltration ratio of eight types of immune cells (Racle et al., 2017). Furthermore, based on the review of relevant literature (Gong et al., 2018; Rowshanravan et al., 2018; Qin et al., 2019; Wang et al., 2020), we explored the difference in the expression of 10 potential immune checkpoint genes between groups. Detailed information is provided in Supplementary Note S1.

Gene set enrichment analysis (GSEA) on the 50 hallmark gene sets was used to identify the potential molecular mechanisms or potential functional pathways that involve the prognostic model (Liberzon et al., 2015). Significant gene sets were based on the following parameters: normalized enrichment score (NES) | >1, nominal p-value < 0.05, and false discovery rate (FDR) q-value < 0.05. Furthermore, 19 m6A-related genes and regulators based on Juan Xu’s research (Li et al., 2019) were further screened and investigated in high- and low-risk groups.

To evaluate the potential of FGL models in clinical practice for P-AML treatment, the half-maximal inhibitory concentration (IC₅₀) of commonly administered chemotherapeutic drugs in the TARGET-combined cohort was calculated using the algorithm R package “oncoPredict” (Maeser et al., 2021). The algorithm predicts the clinical drug response value (IC₅₀) in patients based on gene expression data in tumor, which is derived from the ridge regression model based on drug sensitivity data from the Genomics of Drug Sensitivity in Cancer (GDSC) database.

Kaplan–Meier curves were plotted to estimate overall survival (OS) time and event-free survival (EFS) time, and the data were statistically compared with the log-rank test. Univariate and multivariate analyses were performed by the Cox proportional hazard model. The ROC curve was established, and the area under the curve (AUC) was calculated to determine the predictive value of the prognostic model. Five-fold cross-validation was applied for model evaluation (Supplementary Note S2). All statistical analyses were carried out using R (3.5.2) software and p < 0.05 was taken as being statistically significant. The Wilcoxon test was used to compare the immune scores, immune infiltrate, and expression of genes between groups.

Normalized gene expression data of the TARGET-discovery cohort included 244 ferroptosis-related genes. Univariate Cox regression analysis revealed 39 ferroptosis-related genes with significant prognostic value for P-AML (p < 0.05) (Supplementary Table S2). By ferroptosis-related lncRNA co-expression analysis and univariate Cox regression analysis, we identified 8 prognostic ferroptosis-related lncRNAs out of 2,734 co-expressed lncRNAs (Supplementary Table S2). These 39 ferroptosis-related genes and 8 lncRNAs and merged lists were served as candidate lists for LASSO Cox regression analysis, respectively. After optimal parameter (lambda) selection in the LASSO regression, three prognostic ferroptosis-related signature models with 20 FG components, 6 FL components, and 22 FGL components were built, respectively (Supplementary Figure S1; Supplementary Table S3). The expression level of each gene and LASSO regression coefficient (β) were integrated to calculate the FG, FL, and FGL risk score for each patient (Supplementary Table S4).

In the TARGET-discovery cohort, patients were stratified into high- and low-risk groups using the median risk score as the cutoff value. Statistical differences in overall-survival probability have been identified for the patient groups stratified by the cutoff point of FG, FL, and FGL risk scores, and the survival rate of P-AML patients in the low-risk group of the three models was significantly higher than that in the high-risk group (all p < 0.05) (Figures 2A–C). However, only the FGL model was validated in both the TARGET-validation cohort (p = 0.04, Figures 2D‐F) and the TARGET-combined cohort (p < 0.001, Figure 3A), with significant differences in the survival rate between the two risk groups. The same results have been found for the event-free survival rate (EFS) (Supplementary Figure S2). Thus, the FGL risk score was selected for all subsequent analyses. The distribution and status of OS of the TARGET-combined cohort were then analyzed by ranking the risk scores (Figure 3B). The results showed that patients with higher FGL risk scores had a worse prognosis. Expression profiles of the 22 ferroptosis-related signatures are listed in the heatmap of Figure 3B. The AUC corresponding to 1-, 3- and 5-year OS in the TARGET-combined cohort were 0.70, 0.68, and 0.70, respectively, indicating that the predictive efficiency of the model was good (Figure 3C).

FLT3_ITD_positive, WT1 mutation, 2017 ELN classification system, and FGL risk score were considered significant risk parameters in the univariate analysis (p < 0.05) (Figure 4A), and further multivariate Cox analysis indicated that the FGL risk score was the only independent risk parameter in the discovery cohort (HR = 3.772, 95% CI = 2.529–5.625) (Figure 4B). In the TARGET-combined cohort, the results indicated that the higher FGL risk score was also the only independent poor prognosticator for OS (HR = 1.515, 95% CI = 1.344–1.708) (Figures 4C,D). Several prognostic models have been established or validated in the TARGET cohort recently: LSC17 (Duployez et al., 2019), LSC6 (Elsayed et al., 2020), yang_10_genes (Yang et al., 2020), docking_16_genes (Docking et al., 2021), and cai_3_genes (Cai et al., 2021). Risk scores were further calculated based on the coefficient defined in these studies and used further for correlation analysis with FGL risk score. FGL risk score was negatively correlated with LSC17 risk score and positively correlated with the docking_16_genes model and cai_3_genes model (Figures 5A,B). FGL risk model demonstrated the best predictive performance compared with previous prognostic models, 2017 ELN classification system, and other prognostic molecular characteristics in TARGET-discovery and -combined cohort, respectively (Figures 5C,D). The five-fold cross-validation method was applied to give a robust estimation of the performance of the FGL model. As shown in Supplementary Table S5, in the testing stage, the AUC of the FGL model ranged from 0.693 to 0.741 in the TARGET-discovery cohort and 0.693 to 0.741 in the combined cohort. Collectively, the risk score calculated according to the 22 ferroptosis-related signatures could serve as an independent and stable prognostic parameter for P-AML patients. However, we further explored the predictive ability of the FGL prognostic model for prognosis in adult AML patients. OS was not significantly different in the high-risk and low-risk groups of the TCGA-LAML cohort (p = 0.37, Supplementary Figure S3).

We next investigated the clinical, molecular, and immune features of the low- and high-risk groups in the TARGET-combined cohort, considering the sample size. The results identified a significant difference between the two groups with respect to the distribution of inv (16) mutation, WT1 mutation, 2017 ELN classification system, and FAB subtype (all p < 0.05, Table 2). Of the 42 patients with inv (16), 7.14% (3/42) were in the high-risk group and 92.86% (39/42) were in the low-risk group (p < 0.001) (Figure 6A); WT1 mutation occurs in 6.85% (17/256) of AML patients and 13 of them were in the high-risk group (p = 0.04) (Figure 6A). Consistent with expectations, the number of people identified as adverse-risk by the 2017 ELN classification system was significantly higher in the high-risk group than that of the low-risk group (14.06 vs. 7.81%), while the number of people identified as favorable-risk was on the contrary (27.34 vs. 60.16%). The morphological subtypes in order of frequency among P-AML cases were M2 (66/256, 25.78%), M4 (62/256, 24.22%), M5 (49/256, 19.14%), M1 (33/256, 12.89%), M7 (7/256, 2.73%), and M0 and M6 (3/256 and 1.17%). No M3 case was identified in the present study. M4 subtype was more common in the low-risk group than in the high-risk group (35.16 vs. 13.28%) (Figure 6B). The mutation landscapes in the FGL high- and low-risk groups showed top mutated genes with a frequency above 5%, and NARS mutation was the predominant alteration in both the high- and low-risk groups. The second highest mutation rate (30.49%) was found for KIT in the high-risk group, while the mutation rate is 5.19% in the low-risk group (p < 0.001) (Figures 6C,D).

We identified 112 (43.75%), 103 (40.23%), and 28 (10.94%) patients in the TARGET-combined cohort classified as favorable-, intermediate-, and adverse-risk groups, respectively, according to the 2017 ELN classification system. Our results validated the prognostic significance of the revised 2017 ELN classification system in the TARGET-combined cohort (p < 0.001) (Figures 7A,B). Individuals in the favorable group defined by the 2017 ELN classification system had significantly better OS. However, no significant prognostic difference between the intermediate and adverse groups was found (p = 0.76 for OS; p = 0.96 for EFS). Then, in the 2017 ELN intermediate-risk subgroup of the TARGET-combined cohort (N = 103), we found 70 patients grouped with a high FGL risk score and 33 patients grouped with a low FGL risk score, which could be well risk-stratified by the FGL scoring system (p = 0.003 for OS; p = 0.0047 for EFS, Figures 7C,D).

Examining the molecular trends of divergence across two risk groups using GSEA revealed 8 hallmark gene sets significantly perturbed (Figure 8A, and detailed results for 50 hallmark gene sets are shown in Supplementary Table S6). Six tumor-related hallmarks, apoptosis, hypoxia, TNFA signaling via NFKB, reactive oxygen species pathway, oxidative phosphorylation, and p53 pathway, were significantly enriched in the low-risk group, while bile acid (BA) metabolism pathway and spermatogenesis were found significantly enriched in the high-risk group.

To determine whether the FGL risk score was related to tumor immunity, we next evaluated the correlation between FGL risk and immune score and immune cell infiltration. We observed significantly higher immune score, microenvironment score, stromal score, and ESTIMATE score in the low-risk group (Figure 8B and Figure 8C; all p < 0.05). Strong relationship between infiltration levels of several immune cells and FGL risk score has also been established. As shown in Figure 8D, infiltration levels of cancer-associated fibroblasts (CAFs), CD4_T cells, and endothelial were significantly upregulated in the high-risk group, while infiltration levels of B cells, macrophages, and other cells were significantly upregulated in the low-risk group (all p < 0.05). Notably, we observed a statistically significant difference between the two groups in terms of the expression of several important immune checkpoint genes. The expression of PDL1, CTLA4, TIGIT, and PDL2 were more highly expressed in the high-risk group (all p < 0.05). However, a substantial increase in the expression of TIM3 was found in the low-risk group (p < 0.0001) (Figure 8E).

In addition, we investigated the expression of m6A-related genes between the two risk groups, and the results showed that most of them were more highly expressed in the high-risk group (all p < 0.05) (Figure 8F). Significant correlations were also identified for the expression levels of the 22 ferroptosis-related signatures and m6A-related genes (Supplementary Figure S4). According to m6A2Target, a comprehensive database for the target gene of writers, erasers, and readers (WERs) of m6A modification in a cancer cell line (Deng et al., 2021), 13 out of the 22 ferroptosis-related signatures in the FGL model were potential target genes of WERs of m6A modification in the leukemia cell line (Supplementary Table S7).

Sensitivity to 198 anticancer drugs was compared between the high- and low-risk groups to provide potential treatment guidance for P-AML patients. The results for 153 drugs were not considered for differences analysis because more than 20% of the samples were missing from the predicted IC₅₀ values. Since the predicted values of the samples differ significantly, differences analysis was applied after the removal of potential outliers (extremely high IC₅₀ values). ROUT method was used for outlier identification (setting Q to 5%) (Motulsky and Brown, 2006). The results demonstrated that the IC₅₀ values of SB505124_1194 and JAK_8517_1739 were significantly lower in patients within the FGL high-risk group, which implies that patients within the FGL high-risk group might benefit from treatment using SB505124_1194 and JAK_8517_1739 (Figure 8G).