The sequencing data of CML cohorts GSE13159, GSE144119, GSE4170, and GSE44589 were obtained from the Gene Expression Omnibus (GEO) database. The analysis cohort for this project was the GSE13159 cohort, which consisted of 76 CML samples and 74 normal samples. Raw sequencing data were downloaded and normalized for subsequent analysis. The validation cohort (GSE144119) included 48 newly diagnosed CML samples, 32 remission CML samples, and 17 normal samples that were converted to transcripts per kilobase million (TPM) values. For clinical validation purposes, transcriptome sequencing was performed on five chronic-phase CML samples, five blast crisis samples, and five normal control samples with written consent from patients approved by the Ethics Committee of the Second Affiliated Hospital of Nanchang University; these data were also transformed into TPM values for further validation. To differentiate between other types of leukemia such as acute lymphoblastic leukemia (750 cases), acute myeloid leukemia (542 cases), chronic lymphocytic leukemia (448 cases), and myelodysplastic syndromes (206 cases), a subset of the GSE13159 cohort was utilized. Furthermore, we used the imatinib-treated sample dataset from GSE44589 containing 198 sequenced samples to evaluate treatment response in CML patients. Additionally, single-cell RNA-seq data from the GES76312 cohort were employed to visualize clusters using the uniform manifold approximation and projection (UMAP) algorithm. Finally, we retrieved ferroptosis pathway genes from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).

The “limma” software package was employed for conducting differential expression analysis of FRG. Adjusted p values below 0.05 were considered significant, indicating the presence of differentially expressed FRGs (DEFRGs) between CML and normal samples. Subsequently, we performed Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on these genes using the “clusterProfiler” package (22). To quantify the activity of a biological pathway or gene set, we utilized the Gene Set Variation Analysis (GSVA) algorithm to calculate an enrichment score (23).

The Spearman method was employed for correlation analysis. The STRING database (https://string-db.org/) was used to analyze the PPI of DEFRG. Subsequently, the PPI network was visualized using Cytoscape software.

The estimation of immune cell infiltration was conducted by employing the deconvolution algorithm “CIBERSORT” to accurately quantify the proportions of 22 distinct immune cell types based on the gene expression profiles of individual samples (24).

Weighted correlation network analysis (WGCNA) was employed to identify the genes associated with ferroptosis scores in the GSE13159 cohort (25). Pearson correlation analysis was utilized to construct the adjacency matrix for all matched genes, and the scale-free topology of this matrix was established based on an optimal soft threshold power. Subsequently, the adjacency matrix was transformed into a topological overlap matrix (TOM). By employing the TOM dissimilarity measure, modules consisting of genes exhibiting similar expression patterns were identified through average linkage hierarchical clustering, with a minimum module size set at 30 and a cut height at 0.2. Finally, an evaluation of the correlation between module signature genes (MEs) and ferroptosis score was performed.

To identify diagnostic biomarkers for CML, three machine learning algorithms, namely support vector machine recursive feature elimination (SVM-RFE), least absolute shrinkage selection operator (LASSO), and random forest (RF) were employed to screen the diagnostic FRGs. Additionally, LASSO regression analysis was used to calculate regression coefficients for the diagnostic FRGs, and a CML risk score diagnostic model was constructed using the following formula:

where i represents the specific diagnostic FRG and “Coef” and “ExpGene” denote the regression coefficient and expression value of that particular FRG respectively. By constructing this risk score model, we can further assess the combined diagnostic value of FRGs.

To comprehensively assess inter-individual variations in CML patients, we employed the “ConsensusClusterplus” package to conduct a cluster analysis of CML samples based on the expression profiles of the diagnostic FRGs, aiming to identify distinct molecular subtypes within CML (26). The robustness and stability of the clustering results were confirmed through 1000 iterations. Additionally, principal component analysis (PCA) was utilized for classification validation.

The expression matrix and drug response data of blood cell lines from the Cancer Genome Project (CGP) database were utilized in this study to predict the half-maximal inhibitory concentrations (IC50) of CML samples to TKIs. This prediction was made using the “pRRophetic” package, a computational tool commonly used for such analyses (27). To further investigate the response of different risk score groups towards anti-PD-1 and anti-CTLA4 immune checkpoint inhibitors, we employed the “SubMap” algorithm available at a publicly accessible website called GenePattern. The SubMap algorithm is widely recognized for its ability to forecast treatment responses based on gene expression profiles. To assess the level of immune escape exhibited by tumor cells in CML samples, we computed the TIDE score using an established online resource known as Tumor Immune Dysfunction and Exclusion (TIDE).

We employed miRTarBase, miRDB, and TargetScan databases to predict the binding sites of miRNAs on CML diagnostic ARGs. Subsequently, we filtered out the miRNA-target pairs that were predicted by all three databases. The GSE90773 cohort was utilized to identify differentially expressed miRNAs between CML cells and normal cells, which served as the basis for constructing the miRNA regulatory network.

The CML cell line K562 was cultured in RPMI1640 medium supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin in a humidified incubator saturated with 5% CO2 at 37°C. The K562 cells were exposed to imatinib, and the concentration was gradually increased until the development of K562/IR cells capable of sustained growth in a medium containing 1μM of imatinib. This concentration is considered physiologically relevant and may simulate the peak plasma/serum level of imatinib (5μM). Transcriptome sequencing analysis was conducted on K562, K562/IR, K562/IR control, and erastin-treated K562/IR cells. The processing procedure employed in this study was based on our previous research (28). The concentration of imatinib was gradually increased until the induction of resistant cells was completed. Cell viability was assessed using the cell counting kit-8 (CCK-8) assay. For this assay, 5-e3 cells were seeded in 96-well plates, and each group was repeated three times. After the indicated culture time, 10 μL of CCK8 solution was added, followed by incubation at 37°C for 2 hours. The optical density (OD) value at 450 nm was measured using a microplate reader. Apoptosis detection involved staining cells with the Annexin V-PE/7-AAD apoptosis detection kit and subsequent examination in a flow cytometer. Additionally, reactive oxygen species (ROS) were detected using a fluorescent probe DCFH-DA in flow cytometry. The levels of GSH and GSSH were determined using Solarbio’s BC1175 and BC1185 kits, respectively. Bioss’ AK091 kit was used for GPX4 activity measurement. All reagents were employed following the manufacturer’s instructions. Cell homogenization was performed using lysate buffer to facilitate the reaction between REDOX substances in the sample and reagents, resulting in the formation of adducts that can be quantified through colorimetry.

All analyses were conducted using the R software and corresponding software packages. Differences between two or more groups were assessed using the Wilcoxon rank sum test and the Kruskal-Wallis test, respectively. The diagnostic value of biomarkers was determined through receiver operating characteristic (ROC) curve analysis. A bilateral P-value less than 0.05 indicates a statistically significant difference.

We conducted a comprehensive evaluation of ferroptosis activity and molecular characteristics in CML using transcriptomics analysis. The GSVA algorithm was utilized to calculate ferroptosis scores, revealing significantly lower ferroptosis scores in CML samples compared to normal samples ( Figure 1A ), while the ferroptosis score increased following treatment remission ( Figure 1B ). Patients in blast crisis (BC) exhibited even lower ferroptosis scores than those in the chronic phase (CP) ( Figures 1C, D ) (Due to the limited sample size and vulnerability to individual outliers, although Figure 1C does not exhibit a statistically significant difference, the overall trend persists that BC patients display lower ferroptosis scores compared to CP patients.), and individuals with major molecular responses displayed higher ferroptosis scores compared to non-responders ( Figure 1E ). Single-cell analysis consistently demonstrated a trend of decreased ferroptosis scores in CML patients, particularly those in BC, which subsequently increased after treatment with TKI ( Figures 1F-H ). Differential expression analysis indicated the down-regulation of numerous genes associated with ferroptosis in CML samples ( Figures 1I, J ), including those involved in iron ion homeostasis, mitochondrial outer membrane function, and ligase activity ( Figure 1K ). These differentially expressed genes were primarily enriched in signaling pathways related to ferroptosis, metabolic pathways, mineral absorption, and cysteine and methionine metabolism ( Figure 1L ). PPI network analysis identified STEAP3, TFRC NCQA4 TP53 IREB2 as hub genes within the network formed by these DEFRG ( Figure 1M ). Volcano plot analysis further revealed down-regulation of gene expression for various suppressors of ferroptosis in CML samples ( Figures 1N, O ). Therefore, we speculate that the observed lower ferroptosis score in CML may be attributed to an overall decrease in inhibition of this process within cancer cells indicating their heightened susceptibility towards undergoing cell death through the mechanism of the ferroptosis pathway. Ferroptosis is closely linked to lipid metabolism, and our findings reveal a significant increase in the activity of unsaturated fatty acids such as linoleic acid, arachidonic acid, and α-linolenic acid in CML ( Figure 1P ). Considering that the peroxidation of unsaturated fatty acids is a prerequisite for ferroptosis to occur, this result further supports the hypothesis that CML exhibits heightened susceptibility to ferroptosis. These results collectively indicate an aberrant regulation of ferroptosis in CML samples, which may have implications for the initiation and progression of the disease.

The relationship between the ferroptosis score and the immune microenvironment of CML as well as cancer pathways was further analyzed. It was observed that there were significant associations between the ferroptosis score and key tumor marker pathways, including xenobiotic metabolism, reactive oxygen species pathway, heme metabolism, and epithelial mesenchymal transition activities ( Figure 2A ). Moreover, positive correlations were found with glycolysis and hypoxia, while negative correlations were observed with Notch signaling and WNT beta-catenin signaling. These findings suggest that an increased activity in the ferroptosis pathway is accompanied by enhanced cancer cell metabolism. Immune infiltration analysis revealed a positive correlation between the ferroptosis score and eosinophil infiltration, M0 macrophage infiltration, as well as regulatory T cell (Treg) infiltration; meanwhile, a negative correlation was identified with naive CD4+ T cells ( Figure 2B ). Furthermore, a significant positive correlation was also found between the ferroptosis score and gene expression of immune checkpoints LAG3 and TNFRSF9 ( Figure 2C ), indicating potential immunosuppression among patients with high ferroptosis scores.

To gain a deeper understanding of the underlying mechanisms associated with ferroptosis in CML, we conducted WGCNA to explore the network of co-expressed genes significantly correlated with ferroptosis scores. The cluster dendrogram depicted the clustering characteristics of all CML samples ( Figure 2D ). Figures 2E, F illustrate the scale-free fit exponential and average connectivity analysis for various soft threshold powers. We set the cut height at 0.2 to include modules exhibiting a correlation coefficient greater than 0.8 ( Figure 2G ). Based on an optimal soft threshold power β=15 (unscaled R^2 = 0.9), WGCNA classified the top 5000 genes with the highest standard deviation into 23 independent co-expression modules ( Figure 2H ). The correlograms depicting module-trait relationships revealed that both yellow and sienna3 modules exhibited strong correlations with ferroptosis scores ( Figure 2I ). KEGG enrichment analysis demonstrated that these two modules were enriched in porphyrin and chlorophyll metabolism as well as metabolic pathways ( Figures 2J, K ). Additionally, yellow module genes were found to be associated with nitrogen metabolism, adipocytokine signaling pathway, mTOR signaling pathway, and mitophagy; while sienna3 module genes showed enrichment in hippo signaling pathway, glutathione metabolism, glycolysis/gluconeogenesis, and carbon metabolism. The findings suggest that metabolic reprogramming may contribute to the malignant proliferation of CML cells, while also enhancing the susceptibility of CML cells to ferroptosis by generating higher levels of ROS and unsaturated fatty acids (11, 29).

We conducted further analysis on the diagnostic value of FRG in CML. Three machine learning algorithms, namely LASSO, RF, and SVM-RFE, were employed for dimensionality reduction to select the most informative FRGs. From the DEFRGs, we identified 5, 6, and 6 variables that accurately distinguished CML samples from normal samples, respectively ( Figures 3A–E ). Among these variables, there were five overlapping diagnostic FRGs (ACSL6, SLC11A2, HMOX1, SLC38A1, and AKR1C3) included among them ( Figure 3F ). The expression levels of all five FRGs were significantly downregulated in CML samples compared to normal samples ( Figure 3G ). Using LASSO regression analysis, we developed a risk score model to assess the combined diagnostic value of FRG ( Figure 3H , Supplementary Table S1 ). The risk score levels were significantly elevated in the CML samples ( Figure 3I ). ROC curve analysis revealed high diagnostic AUC values for ACSL6 (0.818), SLC11A2 (0.864), HMOX1 (0.782), SLC38A1(0.783), AKR1C3(0.791), as well as for the risk score (0.920) ( Figure 3J ). The combination of these five FRGs further improved their diagnostic value.

We confirmed the diagnostic value of the five FRGS. In the GSE144119 cohort, we observed a significant decrease in expression levels of all five FRGS in CML samples, which showed partial restoration after treatment response ( Figure 4A ). Furthermore, the risk score levels were significantly increased in CML samples and exhibited a significant decrease after treatment remission ( Figure 4B ), thereby demonstrating the therapeutic evaluation value of FRG. ROC curve analysis revealed that ACSL6, SLC11A2, HMOX1, SLC38A1, AKR1C3, and the risk score model had AUC values of 0.949, 0.934, 0.868, 0.842, and 0.975 respectively ( Figures 4C-H ); thus confirming their diagnostic value in CML cases. In our clinically independent cohort study, we also observed a significant decrease in ACSL6, SLC11A2, HMOX1, and SLC38A1 expression in CML samples while AKR1C3 did not show a significant difference due to small sample size issues ( Figure 4I ). The risk score levels were also significantly increased in CML samples ( Figure 4J ). ROC curve analysis demonstrated an AUC value of 1 for the risk score model ( Figure 4K ). Clinical sample-based sequencing data further verified the high diagnostic value associated with these five FRGs in CML. In conclusion, we have identified highly reliable FRGs which could potentially serve as a novel adjunctive tool for clinical diagnosis and treatment decision-making in patients with CML.

We conducted a comprehensive analysis to evaluate the differential diagnostic value of the five FRGs. The GSE13159 cohort included sequencing data from 750 ALL samples, 542 AML samples, 448 CLL samples, and 206 MDS samples. Interestingly, the expression levels of most FRGs, including SLC38A1, SLC11A2, and HMOX1, were found to be lower in CML samples compared to other types of hematologic tumors. Conversely, ACSL6 exhibited higher expression levels ( Figure 5A ). Furthermore, subsequent calculations revealed that CML samples displayed the highest risk score ( Figure 5B ). ROC curve analysis demonstrated that the risk score effectively distinguished CML from other hematological malignancies with high accuracy (AUC=0.844) ( Figure 5C ). The diagnostic value of FRG has been systematically evaluated, and we have also endeavored to investigate the regulatory mechanisms governing FRG expression. In this study, our focus lies on miRNA, as we aim to construct a miRNA regulatory network to identify potential miRNAs that could inhibit FRG expression by binding to FRG in CML cells ( Figure 5D ).

To comprehensively analyze the biological significance of FRGs in CML, we utilized the expression profiles of the five diagnostic FRGs in CML samples to identify two distinct molecular subtypes, namely Cluster C1 and Cluster C2, employing a consensus clustering algorithm ( Figure 6A , Supplementary Table S2 ). The distribution characteristics of these two molecular subtypes were further confirmed by PCA, revealing significant and discernible differences ( Figure 6B ). Subsequently, through heatmap visualization, it was observed that ACSL6, SLC11A2, HMOX1, and AKR1C3 exhibited up-regulation in subtype C1 while SLC38A1 displayed higher expression levels in subtype C2 ( Figure 6C ). To explore additional distinctions between these subtypes at a biological level, immune infiltration analysis demonstrated that subtype C1 had an increased proportion of CD8+ T cells, follicular helper T cells, activated dendritic T cells, and eosinophils compared to subtype C2 ( Figure 6D ). Furthermore, there were notable variations in the expression levels of immune checkpoint genes; specifically within subtype C1 where PD-L1, CTLA-4, HAVCR2, PD-1, and CD80 showed elevated expressions ( Figure 6E ). This suggests that subtype C1 may exhibit certain immunosuppressive tendencies leading to potential exhaustion of CD8+ T cells. These findings were corroborated by higher TIDE scores for subtype C1 ( Figure 6F ). Conversely, C2subtype appeared more likely to benefit from immunotherapy ( Figure 6G ).

Additionally, our GSVA analysis revealed that the C1 subtype demonstrates heightened activation of signal transduction pathways such as hedgehog signaling and TNFA signaling via NFKB ( Figure 6H ). Moreover, we observed increased activity in cancer-promoting pathways including hypoxia and reactive oxygen species pathway. In contrast, the C2 subtype exhibited elevated activity in proliferation-related pathways such as G2M checkpoint, E2F targets, and MYC targets V1. Notably, C1 displayed a higher ferroptosis score while C2 had a higher risk score ( Figures 6I, J ). Drug prediction analysis indicated that imatinib, nilotinib, dasatinib, and bosutinib demonstrated greater efficacy against subtype C1 compared to subtype C2 ( Figures 6K-N ). These findings will significantly contribute to the development of personalized treatment strategies for patients with CML.

The expression of five FRGs was detected in CML cell lines K562 and imatinib-resistant cell lines K562/IR. In comparison to K562, SLC38A1 expression showed a slight up-regulation in K562/IR, whereas ACSL6, SLC11A2, and AKR1C3 expressions were down-regulated (HMOX1 gene expression was not detected and therefore not shown) ( Figure 7A ). In our study above, our preliminary analysis indicated that CML cells may exhibit sensitivity to ferroptosis, while CML cells in blast crisis demonstrate resistance towards TKI treatment and potentially higher sensitivity. To validate these findings, we conducted in vitro experiments. However, it was observed that the CML cell line K562 did not display sensitivity to erastin-induced ferroptosis ( Figure 7B ); nevertheless, erastin exhibited a certain cytotoxic effect on imatinib-resistant K562 cells (K562/IR) with an IC50 of 5.099 μM ( Figure 7C ). Furthermore, treatment of K562/IR cells with the ferroptosis inhibitor Fer-1 significantly restored cellular viability ( Figure 7D ). Compared to K562 cells, there was a significant increase in ROS levels within K562/IR cells which further escalated after erastin treatment-indicating ROS as a crucial factor for inducing ferroptosis ( Figure 7E ). Additionally, it was discovered that low-dose erastin enhanced the therapeutic sensitivity of imatinib towards K562/IR cells by reducing the IC50 from 3.184 μM to 1.886 μM ( Figure 7F ). Moreover, low-dose erastin promoted apoptosis levels in K562/IR cells treated with imatinib ( Figures 7G, H ). GSH and GPX4 are important indicators of ferroptosis. We found that after erastin treatment, GSH content, GSH/GSSH ratio, and GPX4 enzyme activity of K562/IR cells were significantly decreased, and GPX4 mRNA expression level was slightly increased ( Figures 7I-L ), indicating that erastin inhibited GSH production. In turn, the GPX4 enzyme activity is reduced, which can not inhibit the production of excess ROS, resulting in ferroptosis of K562/IR cells. Finally, we also detected GPX4 expression in K562 and K562/IR cells and CML samples, and the results showed that GPX4 expression in K562/IR cells was lower than that in K562 cells, and there was no significant difference in GPX4 expression between BC-CML and CP-CML samples and normal samples ( Figures 7M, N ). These results suggest that important mechanisms of ferroptosis resistance in CML-resistant cells may not be regulated by GPX4.