Identification of key ferroptosis genes and mechanisms associated with breast cancer using bioinformatics, machine learning, and experimental validation

Data download and preprocessing

Transcriptome data and clinical information related to breast cancer were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). A total of 1204 transcriptome data samples were included, consisting of 1091 breast cancer tumor tissue samples and 113 adjacent normal tissue samples. The clinical information contained 1097 breast cancer patients. After excluding breast cancer patients without transcriptome data and breast cancer patients with overall survival time less than 10 days, clinical data of 1060 breast cancer patients were selected for prognostic analysis. Table 1 shows the statistical description of the clinical data of 1060 breast cancer patients. The transcriptome data were processed as follows: the Ensemble ID in the transcriptome dataset was converted to gene names based on the gene annotation information in the TCGA database, and a new expression matrix was generated by extracting the protein-coding genes in the dataset. The expression matrix was then normalized and standardized. A total of 259 ferroptosis-related genes were obtained from the FerrDB database (http://www.zhounan.org/ferrdb/).

Materials and reagents

The human breast cancer cell line MCF-7 was purchased from the Cell Bank of the Chinese Academy of Sciences. DMEM medium and fetal bovine serum were purchased from Gibco (Grand Island, NY, USA); Op-timen medium was purchased from Procell Life Science and Technology (Wuhan, China); Lipofectamine 2000 transfection reagent was purchased from Invitrogen (Thermo Fisher Scientific, Carlsbad, CA, USA). Small interfering RNA targeting SLC2A1 (siRNA-SLC2A1), pcDNA-TXNIP overexpression plasmid (TXNIP-OE), and pcDNA-ATF3 overexpression plasmid (ATF3-OE) were designed and synthesized by Sangon Biotech (Shanghai, China); SLC2A1, TXNIP, and ATF3 primers were synthesized by Sangon Biotech (Shanghai, China). Primary antibodies against SLC2A1 (#4393), TXNIP (#13113), and ATF3 (#12691) were purchased from Abcam (Cambridge, UK), and anti-actin was purchased from Proteintech Group (Wuhan, China). HRP-conjugated secondary antibody was purchased from Elabscience Biotechnology (Wuhan, China). BCA protein assay kit was purchased from Elabscience Biotechnology (Wuhan, China). TRIzol reagent, PrimeScript RT reagent Kit, and SYBR Green Master Mix were purchased from Thermo Fisher Scientific (Wilmington, NC, USA). Reactive oxygen species (ROS) and Fe²⁺ detection assay kits were purchased from Yeasen Biotechnology (Shanghai, China), glutathione (GSH) detection assay kits were purchased from Beijing Solarbio Science and Technology (Beijing, China), malondialdehyde (MDA) assay kits were purchased from Nanjing Jiancheng Biological Engineering Research Institute (Nanjing, China). CCK-8 assay kit was purchased from New Cell and Molecular Biotech (Suzhou, China). GIEMS staining solution was purchased from Gibco (Grand Island, NY, USA).

Identification of differentially expressed genes associated with ferroptosis

Differential expression analysis of transcriptome data was performed using the R package “limma”. A fold change threshold of 2.5 and a P < 0.05 were set to identify differentially expressed genes (DEGs) between tumor (BC group, n = 1091) and adjacent tissues (NC group, n = 113). To further investigate the relationship between DEGs and ferroptosis, the identified DEGs were compared with ferroptosis-related genes downloaded from the FerrDB database (http://www.zhounan.org/ferrdb). By overlapping the two gene sets, a list of ferroptosis-related DEGs was identified, which may play critical roles in breast cancer.

Functional enrichment analysis

To further investigate the biological functions and pathways associated with ferroptosis, the differentially expressed genes (DEGs) related to ferroptosis were input into the DAVID database (https://david.ncifcrf.gov/) for functional enrichment analysis. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to identify GO terms and KEGG pathways that were significantly enriched (P < 0.05). The results of the functional enrichment analysis were visualized using R software. The identified enriched GO terms and KEGG pathways could provide insights into the molecular mechanism of ferroptosis in breast cancer.

Identification of key ferroptosis-related genes associated with breast cancer by machine learning

Four machine learning algorithms, namely Maximal Clique Centrality (MCC), Random Forest, Support Vector Machine Recursive Feature Elimination (SVM-RFE), and LASSO Regression, were employed to screen ferroptosis-related genes associated with the development, progression, and prognosis of breast cancer. The MCC algorithm was used to explore genes that potentially play a central regulatory role among these ferroptosis-related DEGs by calculating the topological parameters of the nodes in the PPI network. First, the STRING database (http://string-db.org) was utilized to perform Protein-Protein Interaction (PPI) analysis on the ferroptosis-related DEGs. After PPI network construction, the MCC module in the CytoHubba plugin in Cytoscape 3.7.1 software was used to calculate MCC value for each node. The nodes were ranked by MCC value and the top 20 were selected. Random Forest and SVM-RFE algorithms are commonly used effective algorithms for identifying potential disease biomarkers. In this study, R packages “randomForest” was utilized to construct the random forest model with the following key parameters: the number of trees (ntree) was set to 500, the number of variables randomly sampled as candidates at each split (mtry) was set to the default value of square root of the number of features. Based on the feature importance scores assigned by the model, all the ferroptosis-related genes were ranked in descending order of importance. The SVM-RFE approach was implemented using the “e1071” package in R. A linear kernel function was specified for the SVM classifier. Recursive feature elimination was carried out by iteratively removing the features with the smallest weights based on the fitted SVM model. After recursive elimination, the remaining ferroptosis-related genes were ranked in ascending order according to their weights assigned by the final SVM model. Moreover, by integrating clinical information of breast cancer patients and expression matrix of ferroptosis-related genes in tumor samples, LASSO regression was performed using the R package “rms” to identify ferroptosis-related genes closely associated with the prognosis of breast cancer patients. LASSO regression performs feature selection and model building by penalizing the model with L1 regularization to force irrelevant coefficients to zero. This process yields a feature subset with non-zero coefficients. Finally, the common ferroptosis-related genes identified under the four machine learning algorithms were considered as key ferroptosis-related genes associated with the development, progression, and prognosis of breast cancer.

Construction of prognostic risk model based on the key ferroptosis-related genes

To investigate the association between key ferroptosis-related genes and breast cancer prognosis, multivariate Cox regression analysis was employed to construct prognostic risk model. The risk score for each patient was calculated according to the following formula:

$$\text{Risk}\text{Score}={\displaystyle {\sum }_{i=0}^{n}Coe{f}_i\times ex{p}_i}$$

(n represents the number of genes composing the prognostic model, Coef is the risk coefficient of each prognostic gene, and exp represents the standardized expression level of each prognostic gene). Based on the risk scores of each patient, all tumor samples were classified into high-risk (n = 530) and low-risk (n = 530) groups. The R package “survival” was used to plot Kaplan-Meier survival curves, comparing the overall survival rates of high-risk and low-risk groups. Subsequently, Receiver Operating Characteristic (ROC) analysis was performed and ROC curves were plotted to evaluate the accuracy of this prognostic model for predicting the survival time of breast cancer patients.

Cell culture and transfection

The breast cancer cell line MCF-7 was cultured in high-glucose DMEM supplemented with 10% FBS and 1% penicillin-streptomycin at 37°C and 5% CO₂ in a cell culture incubator. Cells in logarithmic growth phase were digested with trypsin and seeded into 6-well plates at an appropriate density. After the cells grew to approximately 60% confluency, Lipofectamine 2000 transfection reagent was used to transfect siRNA-SLC2A1, TXNIP overexpression plasmid, and ATF3 overexpression plasmid mixed with Opti-MEM medium into MCF-7 cells. Empty plasmids (siRNA-NC and NC-OE) were also transfected into MCF-7 cells as negative controls. After incubating for 6 h in the cell culture incubator, the medium was replaced, and the transfection efficiency was detected by q-PCR and Western blotting 24 h later. The cells were divided into four groups according to the treatment: normal control group (Control), SLC2A1 knockdown group (siRNA-SLC2A1), TXNIP overexpression group (TXNIP-OE), and ATF3 overexpression group (ATF3-OE). The cells in each group were used for subsequent experiments 24 h after transfection.

Quantitative reverse transcription-PCR (q-PCR)

First, total RNA was extracted from the transfected cells and normal cells using the TRIzol reagent according to the manufacturer’s instructions. The RNA concentration and purity were determined using a NanoDrop spectrophotometer (Applied Thermo Fisher Scientific, Wilmington, NC, USA). Then, cDNA was synthesized from the extracted RNA using the PrimeScript RT reagent Kit according to the manufacturer’s instructions. Finally, q-PCR was performed using the SYBR Green PCR Master Mix on the ABI 7500 Real-Time PCR System (Applied Biosystems, Wilmington, USA). The primers utilized in q-PCR are detailed in Table 2.

Western blot

To validate the expression of ferroptosis-related genes in breast cancer cell line MCF-7, Western blot analysis was performed. Briefly, after transfection with siRNA-SLC2A1, TXNIP overexpression plasmid, and ATF3 overexpression plasmid, respectively, the cells were harvested and lysed using RIPA buffer containing protease inhibitor cocktail. The protein concentration was determined by the BCA protein assay kit. Equal amounts of protein were separated by SDS-PAGE and transferred onto PVDF membranes. After blocking with 5% non-fat milk, the membranes were incubated with primary antibodies against the selected ferroptosis-related genes and β-actin (loading control) at 4°C overnight. After washing with TBST, the membranes were incubated with HRP-conjugated secondary antibodies at room temperature for 1 h. The protein bands were visualized using ECL reagents and analyzed using ImageJ software. The relative expression levels of the ferroptosis-related genes were normalized to β-actin.

Cell proliferation assay

To assess cell proliferation, single-cell suspensions of the Control group and experimental groups (siRNA-SLC2A1, TXNIP-OE, ATF3-OE) were seeded into 96-well plates at a density of approximately 2000 cells per well and incubated at 37°C and 5% CO₂ in a cell culture incubator for 24, 48, and 72 h. At each time point, 10 μL of CCK8 working solution was added to each well, and the cells were incubated for 30 min at 37°C. The absorbance (A) was measured at 450 nm wavelength according to the instructions of the CCK8 kit.

Colony formation assay

To evaluate colony formation ability, cells from each group were digested with trypsin and made into single-cell suspensions. Approximately 1000 cells were seeded into each well of a 6-well plate and cultured at 37°C and 5% CO₂ in a cell culture incubator for 21 days, with the medium changed every 3 days. After incubation, the cells were stained with GIEMS solution for 15 min and photographed under a microscope.

Detection of cell ferroptosis

The levels of GSH, MDA, and ferrous iron ions (Fe²⁺) in cells were measured according to the instructions of the detection kit. First, the transfected cells were cultured for 24 h, and then collected and lysed using RIPA cell lysis buffer containing 1% protease and phosphatase inhibitor. The lysates were centrifuged at 13000 × g for 5 min, and the supernatant was collected. The cell lysates were mixed with the assay reagents, and the absorbance was measured using a microplate reader. The detection wavelengths for GSH, MDA, and Fe²⁺ were 412 nm, 532 nm, and 593 nm, respectively. Furthermore, the reactive oxygen species (ROS) levels in cells were measured according to the instructions of a ROS detection kit. Transfected cells were cultured for 24 h and collected to prepare a single-cell suspension. A suitable amount of 10 μmol/L DCFH-DA probe solution was added to the cells, and they were incubated at 37°C for 30 min. After washing the cells twice with serum-free cell culture medium, the fluorescence intensity was measured using a flow cytometer with an excitation wavelength of 488 nm and an emission wavelength of 525 nm.

Statistical analysis

Statistical analysis was performed using GraphPad Prism 8.3.0 software. All experimental data were repeated at least three times and presented as mean ± standard deviation (SD). Differences between groups were compared using the independent sample t-test and one-way analysis of variance (ANOVA). A p-value less than 0.05 was considered statistically significant.

Availability of data and materials

The data sets supporting the results of this article are included within the article. All data, models, and code generated or used during the study appear in the submitted article.

Screening of differentially expressed genes associated with ferroptosis

A total of 3031 differentially expressed genes (DEGs) were identified between tumor (BC group) and adjacent tissues (NC group), including 1275 significantly upregulated genes and 1756 significantly downregulated genes (Figure 1A). To further investigate the relationship between DEGs and ferroptosis, the identified DEGs were compared with ferroptosis-related genes downloaded from the Ferroptosis Database (FerrDB). By overlapping the two gene sets, a list of 52 ferroptosis-related DEGs was identified, which may play critical roles in breast cancer (Figure 1B). The expression levels of these ferroptosis-related DEGs are presented in a heatmap (Figure 1C).

Figure 1. Identification of ferroptosis-related DEGs in breast cancer. (A) Volcano plot of all DEGs between tumor (BC group, n = 1091) and adjacent tissues (NC group, n = 113). Red and blue dots represent significantly upregulated and downregulated genes, respectively. (B) Venn diagram showing the overlap between ferroptosis-related genes and DEGs. (C) Heatmap of differentially expressed ferroptosis-related genes in tumor (BC group, n = 1091) and adjacent tissues (NC group, n = 113). The color scale represents the relative expression level of each gene, with red indicating high expression and green indicating low expression.

Functional enrichment analysis

The GO analysis showed that the ferroptosis-related DEGs were mainly involved in biological process (BPs) such as positive regulation of transcription from RNA polymerase II promoter, negative regulation of transcription from RNA polymerase II promoter, and response to hypoxia, as well as cellular components (CCs) such as cytosol, plasma membrane, nucleus, and extracellular exosome. In addition, they were enriched in molecular functions (MFs) such as identical protein binding, protein kinase binding, ubiquitin protein ligase binding, and protein homodimerization activity. These results suggest that the ferroptosis-related DEGs may play important roles in the regulation of transcription, hypoxia response, and protein-protein interactions (Figure 2A).

Figure 2. Functional enrichment analysis of ferroptosis-related DEGs in breast cancer. (A) Gene Ontology (GO) analysis of ferroptosis-related DEGs. The top 10 enriched terms in each category (biological process, cellular component, and molecular function) are shown. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of ferroptosis-related DEGs. The top 10 enriched pathways are shown.

The KEGG pathway analysis revealed that the ferroptosis-related DEGs were significantly enriched in several cancer-related pathways, including Pathways in cancer, Central carbon metabolism in cancer, HIF-1 signaling pathway, and NOD-like receptor signaling pathway. These pathways are known to be associated with cancer development and progression, suggesting that the ferroptosis-related DEGs may play important roles in breast cancer (Figure 2B).

Screening of key ferroptosis-related genes

The 52 ferroptosis-related DEGs were analyzed by constructing a PPI network, which revealed that they interacted with each other (Figure 3A). Based on the maximal clique centrality (MCC) algorithm, we identified the top 20 genes with the best connectivity from the PPI network by calculating the topological parameters of the nodes (Figure 3B). SVM-REF analysis showed that a model containing 19 ferroptosis-related genes had the highest accuracy in differentiating breast cancer tissues from adjacent tissues (Figure 3C). The top 20 genes identified by random forest analysis were selected as the feature ferroptosis-related genes for differentiating breast cancer tissues from adjacent tissues (Figure 3D). In addition, LASSO regression analysis identified 7 ferroptosis-related genes that were closely associated with the prognosis of breast cancer patients (Figure 3E). By integrating the results of the 4 machine learning algorithms, we identified 3 common ferroptosis-related genes, namely TXNIP, SLC2A1, and ATF3 (Figure 3F). Among them, TXNIP and ATF3 were significantly downregulated in breast cancer tissues compared to adjacent tissues (P < 0.001), while SLC2A1 was significantly upregulated in breast cancer tissues compared to adjacent tissues (P < 0.001) (Figure 3G).

Figure 3. Screening of key ferroptosis-related genes in breast cancer. (A) Protein-protein interaction network of ferroptosis-related DEGs. (B) Top 20 genes with the best connectivity identified by the MCC algorithm. (C) SVM-REF analysis of the feature ferroptosis-related genes for differentiating breast cancer tissues from adjacent tissues. (D) Random forest analysis of the feature ferroptosis-related genes for differentiating breast cancer tissues from adjacent tissues. (E) LASSO regression analysis of the ferroptosis-related genes associated with the prognosis of breast cancer patients. (F) Venn diagram showing the common ferroptosis-related genes identified by the 4 machine learning algorithms. (G) Expression levels of TXNIP, SLC2A1, and ATF3 in breast cancer tissues (BC group, n = 1091) and adjacent tissues (NC group, n = 113).

Construction of the prognostic risk model

We used multivariate Cox regression analysis to construct a risk score model, which was calculated as follows: RiskScore = SLC2A1 × 0.813354160440223 − TXNIP × 0.473133422604816 − ATF3 × 0.21764547803389237. The survival status plot of breast cancer patients based on the prognostic model is shown in Figure 4A. Based on the risk score of each patient, all tumor samples were divided into high-risk and low-risk groups. The results showed that the overall survival rate of the high-risk group was significantly lower than that of the low-risk group (P < 0.001), as shown in Figure 4B. ROC curve analysis validated the good predictive ability of the prognostic model for survival. The AUC values of the prognostic model for predicting the survival of breast cancer patients at 365 days, 1095 days, and 1825 days were 0.79, 0.73, and 0.73, respectively, as shown in Figure 4C.

Figure 4. Construction and validation of the ferroptosis-related gene prognostic risk model in breast cancer. (A) Survival status plot of breast cancer patients based on the prognostic model. (B) Kaplan-Meier curve analysis of the overall survival rate of high-risk (n = 530) and low-risk (n = 530) groups. (C) ROC curve analysis of the prognostic model for predicting the survival of breast cancer patients at 365 days, 1095 days, and 1825 days. The AUC values were 0.79, 0.73, and 0.73, respectively.

Functional validation of key ferroptosis-related genes in breast cancer cells

We conducted in vitro transfections with siRNA-SLC2A1, TXNIP overexpression plasmids, and ATF3 overexpression plasmids to further investigate the effects of key ferroptosis-related genes on ferroptosis in breast cancer cells. Transfection results showed that compared to the normal control group (Control), SLC2A1 mRNA and protein levels were significantly downregulated in the SLC2A1 knockdown group (siRNA-SLC2A1) (P < 0.05), while TXNIP mRNA and protein levels were significantly upregulated in the TXNIP overexpression group (TXNIP-OE) (P < 0.05), and ATF3 mRNA and protein levels were significantly upregulated in the ATF3 overexpression group (ATF3-OE) (P < 0.05) (Figure 5A–5E). These results indicate that the transfection efficiency was good in all experimental groups.

Figure 5. Transfection efficiency detection. (A) q-PCR analysis of the transfection efficiency of siRNA-SLC2A1, TXNIP overexpression plasmids, and ATF3 overexpression plasmids. (B–E) Western blot analysis of the transfection efficiency of siRNA-SLC2A1, TXNIP overexpression plasmids, and ATF3 overexpression plasmids. The bar graphs show the mean ± SD of at least three independent experiments (n = 3). ^**P < 0.01, compared to the Control group, the difference is extremely significant.

Cell growth was assessed using CCK8 and colony formation assays. CCK8 results showed that at 24, 48, and 72 hours, cell viability in all experimental groups (siRNA-SLC2A1, TXNIP-OE, and ATF3-OE) was significantly decreased compared to the normal control group (Control) (P < 0.05) (Figure 6A). The results of the colony formation assay showed that the number of cell colonies in all experimental groups was significantly reduced compared to the normal control group (Control) (P < 0.05) (Figure 6B, 6C). These results indicate that downregulation of SLC2A1 or upregulation of TXNIP or ATF3 expression can significantly inhibit the growth of breast cancer cells.

Figure 6. Cell growth ability detection. (A) CCK8 assay was used to detect the effects of SLC2A1, TXNIP, and ATF3 on the cell viability of MCF-7 cells. (B, C) Colony formation assay was used to detect the effects of SLC2A1, TXNIP, and ATF3 on the colony formation ability of MCF-7 cells. The data presented are the mean ± SD of at least three independent experiments (n = 3). ^**P < 0.01, compared to the Control group, the difference is extremely significant.

To investigate the effects of key ferroptosis-related genes (TXNIP, SLC2A1, and ATF3) on ferroptosis in breast cancer cells, we measured the levels of GSH, MDA, and Fe²⁺ in the cells. The results showed that compared to the normal control group (Control), GSH, Fe²⁺, and ROS levels were significantly reduced (P < 0.05), and MDA levels were significantly increased (P < 0.05) in the experimental groups (siRNA-SLC2A1, TXNIP-OE, and ATF3-OE) (Figure 7A–7E). These results indicate that downregulation of SLC2A1 or upregulation of TXNIP or ATF3 expression can induce ferroptosis in breast cancer cells.

Figure 7. Ferroptosis detection in breast cancer cells. (A) GSH levels detection. (B) MDA levels detection. (C) Fe²⁺ levels detection. (D, E) ROS levels detection. The bar graphs show the mean ± SD of at least three independent experiments (n = 3). ^**P < 0.01, compared to the Control group, the difference is extremely significant.