The role of ferroptosis in breast cancer patients: a comprehensive analysis

Breast cancer (BC) affects the breast tissue and is the second most common cause of mortalities among women. Ferroptosis is an iron-dependent cell death mode that is characterized by intracellular accumulation of reactive oxygen species (ROS). We constructed a prognostic multigene signature based on ferroptosis-associated differentially expressed genes (DEGs). Moreover, we comprehensively analyzed the role of ferroptosis-associated miRNAs, lncRNAs, and immune responses. A total of 259 ferroptosis-related genes were extracted. KEGG function analysis of these genes revealed that they were mainly enriched in the HIF-1 signaling pathway, NOD-like receptor signaling pathway, central carbon metabolism in cancer, and PPAR signaling pathway. Fifteen differentially expressed genes (ALOX15, ALOX15B, ANO6, BRD4, CISD1, DRD5, FLT3, G6PD, IFNG, NGB, NOS2, PROM2, SLC1A4, SLC38A1, and TP63) were selected as independent prognostic factors for BC patients. Moreover, T cell functions, including the CCR score, immune checkpoint, cytolytic activity, HLA, inflammation promotion, para-inflammation, T cell co-stimulation, T cell co-inhibition, and type II INF responses were significantly different between the low-risk and high-risk groups of the TCGA cohort. Immune checkpoints between the two groups revealed that the expressions of PDCD-1 (PD-1), CTLA4, LAG3, TNFSF4/14, TNFRSF4/8/9/14/18/25, and IDO1/2 among others were significantly different. A total of 1185 ferroptosis-related lncRNAs and 219 ferroptosis-related miRNAs were also included in this study. From the online database, we identified novel ferroptosis-related biomarkers for breast cancer prognosis. The findings of this study provide new insights into the development of new reliable and accurate cancer treatment options.


Background
Breast cancer (BC) affects the breast tissue and easily metastasizes to the bones and lungs. Globally, BC has very high incidence rates 1 . It is the second most common cause of mortalities among women 2 , accounting for an estimated 24% of diagnosed cases and 15% of mortality cases among them 3 . However, due to the asymptomatic early stages, BC is often diagnosed in the advanced stages, leading to serious outcomes. Early diagnosis and intervention are important for favorable outcomes 4 . Based on genomic and transcriptomic sequences, BC can be allocated into five inherent molecular subtypes, including Luminal A/B, HER-2 enriched, Basal-like and Claudinlow 5 . Despite of the advances in diagnosis, surgery, chemotherapy and radiotherapy, BC is still an extremely malignant tumor with poor survival outcomes 6 . It has been found that BC has become tolerant to some antitumor drugs with apoptotic effects 7 . Thus, it is important to evaluate other potential forms of cell death to overcome chemotherapeutic resistance and to uncover new biomarkers and therapeutic targets for improving BC prognosis.
Studies on ferroptosis in tumors are on the ascendancy. Ferroptosis is an iron-dependent cell death mode that is characterized by intracellular accumulation of reactive oxygen species (ROS), which is distinct from apoptosis and autophagy in morphology 8 . The main mechanism of ferroptosis is that, under the action of divalent iron or ester oxygenase, it catalyzes the unsaturated fatty acids highly expressed on the cell membrane to cause liposome peroxidation, thereby inducing cell death. Dysregulation of iron metabolism can enhance the risk of cancer and facilitate tumor growth. Compared to normal cells, cancer cells are more dependent on iron, a phenomenon known as iron addiction 9 . Activation of ferroptosis pathways can overcome the drug resistance mechanisms associated with the existing chemotherapeutic agents, and ultimately provide new therapeutic options for treating cancer. Ferroptosis enhances metabolic dysfunctions that lead to the production of cytoplasmic and lipid ROS, independent of mitochondria while dependent on NADPH oxidases in certain cellular environments 10 . A recent study devised an Fe 2+ -based metal-organic skeleton, which transfers Fe 2+ to cancer cells, thereby triggering Fenton reaction, producing excess ROS, and inducing ferroptosis in BC 11 . Moreover, ferroptosis resistance can be mediated by a prominin2-MVB-exosome-ferritin pathway in BC cells. Prominin2 is a pentaspanin protein involved in regulation of lipid dynamics 12 . However, very few studies have systematically evaluated ferroptosis-related molecular signatures to predict overall survival (OS) in BC patients. We established a prognostic multigene signature based on ferroptosis-related differentially expressed genes (DEGs) obtained from The Cancer Genome Atlas (TCGA). The signatures were validated in the Gene Expression Omnibus (GEO) cohort. In addition, the role of ferroptosisrelated miRNAs, lncRNAs, and immune responses in BC were also evaluated.

Results
Enrichment analysis of ferroptosis-related genes BP of target genes were enriched in response to oxidative stress, multicellular organismal homeostasis, cofactor metabolic processes, and response to metal ion activities among others; MF were mainly enriched in ferric iron binding, aldo-keto reductase (NADP) activity, and oxidoreductase activity among others; CC were mainly enriched in the regulation of epithelial cell proliferation, lipid droplet, apical plasma membrane, and astrocyte projection among others; KEGG were mainly enriched in ferroptosis, hypoxia-inducible factor (HIF)−1 signaling pathway, NOD-like receptor signaling pathway, central carbon metabolism in cancer, Kaposi sarcoma-associated herpesvirus infection, and PPAR signaling pathway (Fig. 1). A total of 65 (28 upregulated and 37 downregulated) DEGs were identified.

Building the prognostic ferroptosis-based signature
To screen for prognostic ferroptosis, differentially expressed ferroptosis genes were evaluated in the univariate COX analysis. Thirty-three ferroptosis genes that were found to be significant in the univariate COX analysis were included in the multivariate COX analysis ( Fig. 2A). As a result, 15 differently expressed genes (ALOX15, ALOX15B, ANO6, BRD4, CISD1, DRD5, FLT3, G6PD, IFNG, NGB, NOS2, PROM2, SLC1A4, SLC38A1, and TP63) were selected as independent prognostic factors for BC patients (Table S1). Thus, the formula for our model was: Risk Score = (0.177 × expression ALOX15 ) + (0.535 × expression ANO6 ) + (0.473 × expression CISD1 ) + (3.530 × expression DRD5 ) + (0.189 × expression G6PD ) + Fig. 1 Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) outcomes of ferroptosis-related genes through the R language "ggplot2" and "org.Hs.eg.db" package (adjust p value <0.05). The org.Hs.eg.db package implements the conversion between gene symbols and entrezIDs based on mget function. The X-axis represents the gene ratio, that is, the number of genes enriched in each GO/the total number of genes and the Y-axis represents the specific function or pathway that is enriched. The size of the circle represents the number of genes enriched on each GO, and the redder the color represents the higher the degree of enrichment on GO.  ). In addition, the risk score of each patient in our study cohort was calculated.

Survival outcomes and multivariate examination
The effect of the expression of 15 ferroptosis genes on patient survival was determined through Kaplan-Meier curve. It was found that the expression of ALOX15, ANO6, BRD4, CISD1, FLT3, IFNG, NOS2, PROM2, and TP63 significantly affected patient's OS (Table S1). In addition, as shown in Fig. 2B, the Kaplan-Meier curves of the two cohorts established the predictive value of the molecular signature. The high-risk group exhibited a poorer survival outcome when compared to the low-risk group (p = 1.645e−13). We also established patient's risk survival status plot, and as the patient's risk score increased, the mortality rate also increased (Fig. 2C). Using our molecular signature, we predicted the AUC values of the patients' 1, 3, 5-year survival rate as 0.74, 0.71, and 0.75, respectively (Fig. 2D). The ROC curves were also used to determine if the expression patterns of survival-related ferroptosis could be used as early markers for the prediction of the occurrence of BC. We found an AUC of 0.719, implying that the sensitivity and specificity of this prognostic model is moderate (Fig. 2E). We then identified risk factors to aid in the establishment of a 15-ferroptosis based prognostic model. The risk score and age were confirmed as independent prognostic factors for OS (Fig. 3). A combination of clinical pathology and prognostic models were used to construct a nomogram (C-index: 0.764). The combination enhanced the predictive value of OS at 1, 3, and 5 years and, therefore, could be used to inform clinical management (Fig. 4). Furthermore, we used GEO to verify the accuracy of our prognostic signature with an AUC of 0.754 (Fig. S1).

Gene set enrichment analyses
The samples were allocated into high-and low-risk groups to distinguish between their potential functions and elucidate the significant survival differences using GSEA. Most of the metabolic-related pathways such as citrate cycle tca cycle, glycolysis gluconeogenesis, pentose phosphte pathway, n glycan biosynthesis, propanoate metabolism, o glycan biosynthesis, steroid biosynthesis, starch and sucrose metabolism, biosynthesis of unsaturated fatty acids ( Fig. S2 and Table 1). Fig. 2 Evaluate the accuracy of the ferroptosis-related gene signature based on TCGA. A Univariate COX analysis identified 33 ferroptosisrelated genes; red means upregulated and greed means downregulated genes; B Kaplan-Meier indicated that high-risk group (red) exhibited a poorer survival outcome when compared to the low-risk group (blue); C The top graph represents the classification of patients into high (red) and low (blue) risk groups based on risk scores, the middle graph represents as the patient's risk score increased, the mortality rate also increased, the bottom heatmap represents the change of 15 ferroptosis-related genes as the risk score increases; D The AUC values of the patients' 1, 3, 5-year survival rate; E the AUC of the signature was 0.719, exhibiting superior performance than the traditional clinicopathological features in predicting the prognosis.

Immunity and gene expression
The heatmap of immune responses based on five algorithms is shown in (Fig. 5). The correlation between immune cell subpopulations and related functions were evaluated using ssGSEA. T cell functions, including CCR score, immune checkpoint, cytolytic activity, HLA, inflammation promotion, para inflammation, T cell costimulation, T cell co-inhibition, and type II INF response were significantly different between the low-risk and highrisk group of the TCGA cohort (Fig. 6A). Studies have characterized the importance of checkpoint inhibitorbased immunotherapies 13 . Therefore, we determined the differences in the expression levels of immune checkpoints between the two groups and found that the expression levels of PDCD-1 (PD-1), CTLA4, LAG3, TNFSF4/14, TNFRSF4/8/9/14/18/25, IDO1/2 et al were significantly different (Fig. 6B).

miRNA and lncRNA analysis
We identified 1185 ferroptosis-related lncRNAs. The relationship between ferroptosis-related genes and lncRNAs is shown in Fig. 7A. Twelve differentially expressed lncRNAs , and OTUD6B-AS1) were selected as independent prognostic factors (Table S2). The multivariate COX analysis results demonstrated that the risk score based on ferroptosis-related lncRNAs and age were confirmed as independent prognostic factors for OS (Fig. 7B). The Kaplan-Meier curves, risk survival status plot, and AUC value are presented in Fig. S3. Targetscan 14 and  miranda 15 database through the ENCORI platform were used to predict the potential binding miRNA of mRNA and only satisfying both two databases been considered potential miRNA, ultimately, 219 ferroptosis-related miRNAs were identified. In addition, Cytoscape was used to establish mRNAs and miRNA networks. The plug-in Molecular Complex Detection (MCODE) of Cytoscape was used to detect densely connected regions in the networks. Four hub miRNAs (hsa-miR-181d-5p, hsa-miR-152-3p, hsa-miR-497-5p, and hsa-miR-155-5p) were identified in the network (Fig. 8).

Discussion
Ferroptosis can eliminate the adaptive characteristics of malignant cells, remove cells that cannot obtain key nutritional factors, are infected or damaged by environmental changes. It plays a key role in inhibiting tumorigenesis and is a potential strategy for tumor treatment. In this study, we identified a novel and useful ferroptosisrelated prognostic gene signature based on the TCGA dataset. The signature was validated in the GSE7390 dataset of BC patients. Meanwhile, we also evaluated the role of ferroptosis-related miRNAs and constructed a ferroptosis-related lncRNAs signature. Potential immune responses involving the tumor microenvironment, tumorinfiltrating immune cells, immune function, and immune checkpoints were also evaluated. Therefore, this study informs on the potential biomarkers and therapeutic targets of the ferroptosis signaling pathways.
Analysis of the 259 ferroptosis-related genes and their potential functions revealed that KEGG were mainly enriched in the HIF-1 signaling pathway, NOD-like receptor signaling pathway, central carbon metabolism in cancer, and PPAR signaling pathway. It has been reported that the efficacy of FG-4592 (an inhibitor of prolyl hydroxylase of HIF) pretreatment is attained by suppressing ferroptosis in the advanced stages of kidney injury through Akt/GSK-3β-mediated Nrf2 activation 16 . In addition, Yu et al. 17 reported the correlation between central carbon catabolic pathways and ferroptosis in HNSCC cisplatin resistance. PPARα activity is essential for the promotion of ferroptosis in cancer by MDM2 and MDMX (tumor suppressor p53 negative regulators) 18 . In this study, 15 ferroptosis genes (ALOX15, ALOX15B, ANO6, BRD4, CISD1, DRD5, FLT3, G6PD, IFNG, NGB, NOS2, PROM2, SLC1A4, SLC38A1, and TP63) were selected as independent prognostic factors for BC. By upregulating GPX4 and down-regulating arachidonic 15lipoxygenase (ALOX15), the silenced GSK-3 protein blocked erastin-induced ferroptosis and secreted few ROS molecules in BC. The role for ALOX15B (ALOX15 type B), and to a lesser extent ALOX15, in cholesterol homeostasis in human macrophages have been documented 19 . ANO6 (TMEM16F) serves as a Ca 2+ activated chloride and scrambles membrane phospholipids to disclose phosphatidylserine at the cell surface. Activation of ANO6 is essential during ferroptosis mediated cancer cell death 20 . A high expression of BRD4 in cancer tissues has been associated with poor prognosis. Ferroptosis was found to be induced by BRD4 knockdown 21 . CISD1 (mitoNEET), an iron-containing outer mitochondrial membrane protein, negatively regulates ferroptosis cancer cell death, and promotes the proliferation of cancer cells, adjuvant tumor growth and metastasis 22,23 . Dopamine receptor D1-like receptors (DRD5) regulates tumor behavior 24 . FLT3 inhibitors, as potent protectors, were found to inhibit glutamate toxicity and ROS generation 25 . G6PD is associated with mitochondrial dysfunction 26 and human neuroglobin. It binds proteins of the cellular iron metabolism under ferroptosis stress 27 . Prominin2  (PROM2) can suppress ferroptosis sensitivity 28 . SLC38A1 which is an important regulator of lipid peroxidation 29 and TP63 gene polymorphism enhances the risk of BC development 30 . Five immune algorithms were then used to assess the relative immune cell infiltrations of each sample, and therefore, establish the efficacy of gene signatures for predicting the OS of BC patients. Our prognostic signature was found to be a reliable tool for predicting the prognosis of BC patients. The induction of ferroptosis combined with immune checkpoint inhibitors (ICIs) showed a synergistically enhanced antitumor activity, even in ICI-resistant tumors 31 . A limited number of studies have evaluated the relationship between ICI and ferroptosis. This study provides important avenues for future research, especially PD-1 and CTLA4 targets. microRNA (miRNA) and long non-coding RNA (lncRNA) are key mediators for regulating ferroptosis.
Twelve ferroptosis-related lncRNAs were used to build the lncRNA signature. Four hub miRNAs that could be associated with ferroptosis were identified. LINC00618 promotes ferroptosis by elevating lipid ROS and iron levels, and by suppressing the expression levels of solute carrier family 7 member 11 (SLC7A11) in human leukemia 32 . lncRNA ZFAS1 knockdown attenuates ferroptosis by sponging miR-150-5p to downregulate SLC38A1 expression in pulmonary fibrosis 29 . Since Nrf2 inhibits iron absorption, limits ROS production and up-regulates the expression of SLC7A11, miRNA can affect ferroptosis by regulating the Nrf2 expression [33][34][35] . miRNA is also involved in the regulation of iron transport, storage, utilization and absorption. It has been documented that miR-20a regulates the expression of iron exporter ferroportin (FPN) and decreases iron export, leading to intracellular iron retention, which in turn favors cell proliferation in lung cancer 36 . Fig. 6 The correlation between immune cell subpopulations and related functions based on ssGSEA in high (red) and low (blue) risk groups, and the differences in the expression levels of immune checkpoints between the high (red) and low (blue) risk groups. A Immune function; B Immune checkpoints. Ferroptosis, as a form of cell death provides a new therapeutic strategy for tumor treatment. However, various issues such as the interconnection between ferroptosis and other cell deaths, and host immunogenicity have not been established. In this study, we integrated ferroptosis biomarkers to assess treatment outcomes and to reveal precise therapeutic targets for BC. However, more studies using independent cohorts and predictive ferroptosis functional experiments should be performed to validate our signature. The limitations of our study are that; our outcomes were not clinically validated and the insufficient data from the small sample size limits the statistical power of our findings.
In summary, we identified novel ferroptosis-related biomarkers markers for BC prognosis. This study provides novel insights for the development of reliable and accurate cancer treatment tools.

Data collection
RNA-seq (113 normal and 1109 tumor) data and miRNAs-seq (104 normal and 1103 tumor) data for BC patients was obtained from TCGA. The RNA-seq data of the GSE7390 (198 tumor tissues) dataset obtained from the GEO database was validated and its clinical characterizes presented in Table 2. Ferroptosis-related genes were obtained from FerrDb, a web-based resource. FerrDb provides a comprehensive and up-to-date list of ferroptosis markers, regulatory molecules and ferroptosis-disease associations 37 . We obtained 259 (Driver: 108; suppressor: 69; marker: 111) ferroptosis-related genes. Pearson correlation analysis was used to determine the correlations among ferroptosis-related genes and lncRNAs. A ferroptosis-related lncRNA association was considered positive with a correlation coefficient |R 2 | > 0.3 and p < 0.001.
The collected clinicopathological data included gender, age, stage, TMN classification, survival status and the number of survival days. FDR < 0.05 and |log 2 FC| ≥ 1 were determined as differentially expressed ferroptosis-related genes. Moreover, we evaluated the functions of the identified ferroptosis-related genes using Gene ontology (GO) function analysis (biological processes (BP), molecular functions (MF), and cellular components (CC)) and Kyoto Encyclopedia of Genes and Genomes (KEGG) through the R language ggplot2 package.

Development of the ferroptosis-related prognostic gene signature
Lasso-penalized Cox regression and Univariate Cox regression analyses were used to build the ferroptosisrelated prognostic gene signature. The regression coefficient β value (coef) of each gene is obtained in the Cox model, and the positive sign indicates that the risk (mortality) is higher, so the prognosis is worse. The signature was defined as risk score = (CoefficientmRNA1 × expression of mRNA1) + (CoefficientmRNA2 × expression of mRNA2) + ⋯ + (CoefficientmRNAn × expression mRNAn). miRNA-based signature and lncRNA-based signature were also used the same formula to calculate. The associated clinical data for BC patients were also downloaded and evaluated. Based on the median, we denoted the data as either low-risk (<median number) or high-risk (≥median number) group. The Kaplan-Meier survival analysis was used to determine the survival rate.

Building a predictive nomogram
We developed a nomogram for predicting the occurrence of cancer events, such as recurrence or death. The time-dependent receiver operating characteristic (ROC) curve was used to evaluate the predictive accuracy of the developed prognostic signatures for BC patients. Univariate/multivariate Cox regression analyses were used to evaluate the connection between ferroptosis-related genes and clinicopathological manifestations. To understand the mechanisms underlying the gene signatures in the KEGG pathways, we used Gene set enrichment analyses (GSEA) to search for rich terms in C2 in the TCGA-BC cohort. P < 0.05 and false discovery rate (FDR) q < 0.05 were considered to be statistically significant.

Immune analysis and gene expression"
We compared CIBERSORT 38,39 , ESTIMATE 40 , MCPcounter 41 , single-sample gene set enrichment analysis (ssGSEA) 42 , and TIMER 43 algorithms to assess the CC or cell types of immune responses in heterogeneous samples among high-risk and low-risk groups based on ferroptosis-related gene signatures. ssGSEA was then used to quantify tumor-infiltrating immune cell subgroups during immune responses and immune functions between the two groups. Potential immune checkpoints were extracted from previous literature.

Statistical analysis
Statistical analyses were performed using R version 3.5.3 and Bioconductor. The unpaired Student's t test and the Wilcoxon test were used to evaluate the normal distribution variables and the non-normal distribution variables, respectively. To identify differential genes in DEGs assessments, the Benjamini-Hochberg method was used to convert the p value to FDR. ssGSEA-normalized BC DEGs data were compared to a genome using "GSVA" (Rpackage). For each statistical analysis, a p ≤ 0.05 was considered to be statistically significant. The remaining value after removing the sample with missing data from the total number of people.