A novel ferroptosis-related long noncoding RNA signature for relapse free survival prediction in patients with breast cancer.

Ferroptosis is the process of cell death dependent on iron. Growing evidence suggests that ferroptosis plays vital roles in the biological process of many cancers. However, just a small number of ferroptosis-related lncRNAs have been explored in depth. Ferroptosis-related lncRNAs in breast cancer (BC) were identified by co-expression analysis based on The Cancer Genome Atlas database (TCGA). The whole set was divided into a training set and a test set with a 1:1 ratio. Univariate Cox regression and LASSO analyses were performed to establish a signature in the 3 sets. Kaplan-Meier analysis and receiver operating characteristic (ROC) for the 3 sets validated the effectiveness and robustness of the signature. Besides, we also explore the relationship between this and clinical characteristics, immune cell infiltration and tumor microenvironment. Meanwhile, the nomogram was drawn by screening indicators of independent recurrent prediction. Finally, we evaluated the relationships between the signature and tumor microenvironment. We identified 391 ferroptosis-related lncRNAs and constructed a 5 lncRNAs-based signature in the training, test, and whole sets, stratifying patients into high-risk and low-risk groups. According to survival analysis, patients in the high-risk groups had worse relapse free survival (RFS) compared to the low risk-groups. The ROC curves indicated that the recurrent signature had a promising predictive capability for BC patients. Moreover, an independent factors-based nomogram model could offer the quantitative prediction and net benefit for the recurrence of BC patients. Finally, the microenvironment, including tumor mutational burden (TMB), immune cell functions and immune checkpoints, showed big differences between the 2 groups. The 5 ferroptosis-related lncRNAs and their signature might be novel promising biomarkers and immunotherapy targets for patients with BC.

1. Introduction

Breast cancer (BC) remains one of the most common malignant diseases worldwide and one of the leading causes of cancer death in women.[ Annually, over 279,100 women are newly diagnosed with BC and 42,690 individuals died from BC in 2020, and its morbidity for BC exhibits a younger trend in recent years.[ Currently, the main treatments for BC are mainly divided into local interventions (surgery therapy and radiation therapy) and systemic therapy (chemotherapy, endocrine therapy, and targeted therapy). With the progress made in treatment methods for BC, the 5-year survival rate and 10-year survival rate of BC patients have been significantly improved.[ However, it was estimated that more than 30% of early-stage BC would develop distant metastasis, and many BC patients were already in the late stage at diagnosis, whose average 5-year survival rate remains poor.[ More importantly, even with aggressive treatments, nearly 30% of BC patients can still occur recurrence due to tumor heterogeneity. Reducing the morbidity and mortality of BC continues to be a challenge for worldwide public health particularly for metastatic BC. Therefore, exploration of the molecular mechanisms underlying BC occurrence and development is critical to find effective therapeutic targets and molecular biomarkers of recurrence.

Ferroptosis, a recently discovered form of type of cell death with autophagic and caspase-independent cell death, is biochemically and genetically distinct from apoptosis and necrosis, and autophagy.[ This process is characterized by iron-dependent and reactive oxygen species with a variety of cytological changes, including smaller volume and increased membrane, a vanished mitochondria crista, and an increased membrane density.[ Ferroptosis is considered an adaptive function to eradicate tumor cells. It exerts antitumor effects by clearing cells that are undernourished in the environment or damaged through infection and ambient stress.[ Further investigations indicated that patients treated with the combined therapy of ferroptosis inducer erastin and chemotherapeutic drugs often gained remarkable anticancer efficacy and better outcomes than traditional chemotherapy alone.[ It has been proved that ferroptosis level was associated with prognosis and therapy resistance in hepatocellular carcinoma and colorectal cancer.[ Likewise, ferroptosis plays an important role in BC. Prominin has been demonstrated to reduce the ferroptosis of BC cells and promote tumor progression by inhibiting iron exchange.[ Additionally, another study has reported that a new treatment option for refractory BC may be developing from siramesine and lapatinib, which appear to activate the ferroptosis process in BC cells.[ These researches indicate that ferroptosis plays a significant role in BC progression and provides a new a novel method of treating BC. Long noncoding RNAs (lncRNAs) are a type of noncoding RNAs with a length of > 200 nucleotides and are involved in various biological activities.[ Accumulating researches indicated that lncRNAs contributed to various types of cancer, and were involved in cancer cell processes such as proliferation, migration, apoptosis.[ In addition, lncRNAs also could regulate ferroptosis via multiple mechanisms. One recent study revealed that lncRNA LINC00336 inhibited ferroptosis in lung cancer by functioning as a competing endogenous RNA via binding to ELAVL1.[ In a related study, Mao et al found that lncRNA P53RRA promotes ferroptosis in lung cancer via higher retention of p53 in the nucleus.[ For BC, Chen et al found that metformin may induce ferroptosis by inhibiting autophagy through lncRNA H19, and this discovery may potentially facilitate the development of new cell-based therapies for BC.[ This discovery may facilitate the development of novel therapies for the treatment of BC. Given the lncRNAs and ferroptosis play important roles in BC, more biomarker of ferroptosis-related lncRNAs for BC need to be identified, which not only helps understand specific mechanism of ferroptosis in BC but also improves the prediction of BC patient’s prognosis. However, only a minority of lncRNAs in ferroptosis-related lncRNA have been examined in depth, and there are few reports about systematically evaluating ferroptosis-related lncRNAs and their recurrence prediction value in BC.

In this study, we investigated the expression of ferroptosis-related lncRNAs across BC in The Cancer Genome Atlas (TCGA) and ferroptosis datasets. On this basis, we established a ferroptosis-based lncRNAs recurrent signature to predict the relapse free survival (RFS) of individual BC patients. We then explored the roles of tumor mutational burden (TMB), immune cell functions and immune checkpoints in BC recurrence. Thus, we hope to highlight the potential functions of ferroptosis-related lncRNAs, as well as provide a promising recurrent tool for BC.

2. Materials and Methods

2.1. Data collection and processing

The gene expression data of BC patients, somatic mutation data and the corresponding clinical information were downloaded from TCGA-BRCA program (https://portal.gdc.cancer.gov/). Ferroptosis-related genes were extracted from FerrDb database.[ Ultimately, 259 ferroptosis-related genes were obtained, which were presented in Supplemental Digital Content (Table S1, http://links.lww.com/MD/G924). Then, we screened ferroptosis-related differentially expressed genes (DEGs) between BC samples and cancer-adjacent normal samples with |log2foldchange (FC)| >1 and an adjusted P value < 0.05 as the thresholds using the R software (4.1.0) “limma” package. We conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the R package “clusterProfiler” to evaluate the functions associated with the ferroptosis-related DEGs. The correlation between ferroptosis-related DEGs and lncRNAs was calculated using Pearson correlation analysis. A lncRNA with correlation coefficients (|R|) more than 0.4 and a P value <0.001 was considered as ferroptosis-related lncRNA. The corresponding clinicopathological information of the patients are summarized in Supplemental Digital Content (Table S2, http://links.lww.com/MD/G925). All data used in this study were extracted from public databases, so there was no requirement for the approval from the local ethics committee.

2.2. Establishment and validation of the ferroptosis-related lncRNAs-based signature

The whole set was randomly split into a training set and a test set in the ratio of 1:1. In the train set, univariate Cox regression analysis was performed to evaluate recurrent associated ferroptosis-related lncRNAs. The lncRNAs with a P value < 0.01 were included in the least absolute shrinkage and selection operator (LASSO) regression cox analysis with 10-fold cross-validation to establish a ferroptosis-related lncRNAs signature. L1-norm was applied to penalize the weight of the feature. The minimum lambda was defined as the optimal value for lncRNAs number. The risk score for each patient was calculated using the following formula: risk score = β1(lncRNA1) *expr(lncRNA1) + β2(lncRNA2)*expr (lncRNA2) + …… + βn (lncRNAn)*expr (lncRNAn), where βn means the coefficient of lncRNAs, expr is the expression of lncRNA. According to the median risk score in the training set, the patients in the 3 sets were classified into high-risk and low-risk groups, respectively. Kaplan-Meier log-rank test was used for evaluation of RFS difference between the high-risk and low-risk groups. Time-dependent receiver operating characteristic (ROC) curves were used to display the performance of the signature.

2.3. Exploration of signature in functional enrichment and microenvironment

Gene set enrichment analysis (GSEA) 4.1.0 (https://www.broadinstitute.org/gsea/index.jsp)[ was utilized to interpret the different functional phenotypes of gene expression between the high-risk and low-risk groups based on the Hallmark gene sets. TMB and mutation counts were computed from somatic mutation data. R package “maftools” was used to display the mutation information of patients. At the same time, we assessed the infiltrating immune cell composition and function between high-risk and low-risk groups by multiple algorithms including TIMER,[ CIBERSORT,[ QUANTISEQ,[ MCPCOUNT,[ XCELL,[ EPIC.[ Expression of potential immune checkpoints in the 2 groups was also compared.

2.4. Building a predictive nomogram

Univariate and multivariate Cox analyses were conducted to explore the relationship between RFS and age, tumor location, tumor (T), lymph node (N) metastasis, distant metastasis (M), stage, estrogen receptor (ER), progesterone receptor (PR), human epidermal growth factor receptor 2 (HER2) and risk score. Variables which had a P value < 0.05 in the univariate and multivariate analyses were included for nomogram construction. The concordance index (C-index), calibration curves, decision curve analysis (DCA) were applied to evaluate the accuracy and benefit of the model.

3. Results

3.1. Differential expression and enrichment analyses for ferroptosis-related genes

According to the criteria of |log2FC|>1 and adjusted P value < 0.05, 77 ferroptosis-related DEGs (Fig. 1A) were obtained (38 downregulated and 39 upregulated). Then, GO and KEGG analysis for ferroptosis-related DEGs was conducted to uncover the function and biological pathways. From the results of GO enrichment analysis (Fig. 1B), ferroptosis-related DEGs were mainly concentrated in response to oxidative and chemical stress (biological processes terms), lipid droplet, apical plasma membrane and apical part of cell (cellular components terms), antioxidant activity, oxidoreductase activity and oxygen binding (molecular functions terms). KEGG analysis showed that the ferroptosis-related DEGs were assigned into hypoxia-inducible factor (HIF)-1 signaling pathway and Amebiasis (Fig. 1C).

Figure 1.

Related analyses for ferroptosis-related DEGs. (A) Volcano plot. (B,C) GO and KEGG enrichment analysis. GO = Gene Ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes.

3.2. Identification of a ferroptosis-related LNCRNAS-based signature in patients with BC

A total of 391 ferroptosis-related lncRNAs were obtained through co-expression analysis. Meanwhile, we constructed a co-expression network to reflect the relationship among ferroptosis-related genes and lncRNAs (Fig. 1D). The BC patients were allocated randomly to the training (531 patients) and test (530 patients) sets according to a 1:1 ratio. Next, we carried out univariate Cox regression analysis for the 391 ferroptosis-related lncRNAs. Among them, 21 lncRNAs were significantly associated with the RFS of BC patients from the training set (P value < 0.05). LASSO regression with 10-fold cross-validation was adopted to screen hub from the 21 ferroptosis-related lncRNAs, and 5 lncRNAs (LINC01235, LINC02166, AL133467.1, TGFB2-AS1, LINC02266) were filtered with optimal lambda value. Meanwhile, a signature was established based on the 5 lncRNAs using LASSO cox regression analysis (Fig. 2). The coefficients of the 5 lncRNAs were showed in Supplemental Digital Content (Table S3, http://links.lww.com/MD/G926). BC patients were ranked based on risk score value and categorized into high-risk and low-risk groups using the intermediate risk score (Fig. 3A). Kaplan-Meier curve analysis showed that BC patients in the high-risk group had worse RFS compared to the low risk-group (Fig. 4A). In addition, time-dependent ROC curves demonstrated an area under the curve (AUC) of 5 lncRNAs-based signature was 0.743 for 5-year RFS (Fig. 4D). To validate the robustness of the signature, we built the same signatures in the test and whole sets with the same algorithm cutoff (Fig. 3B,C). Similarly, high-risk groups had a worse RFS than low-risk groups in the test and whole sets (Fig. 4B,C). The AUCs of the signatures for 5-year RFS were 0.692 and 0.708 in test and whole sets (Fig. 4E,F). To further explore whether the signature is applicable to predict the recurrence of patients with different clinicopathological features, we performed stratified analyses for the whole set with Kaplan–Meier curves. As a result, except for stage HER2 positive and left location patients, the RFS of the low-risk group continued to be superior to that of the high-risk group in all the other subgroups (Fig. 5).

Figure 2.

Construction of A signature for BC patients based on ferroptosis-related lncRNAs in the training, test and whole sets. (A) Selection of the optimal parameter (lambda) via 10 times cross-validation. (B) LASSO coefficient profiles of 20 ferroptosis-related lncRNAs. BC = breast cancer, LASSO = least absolute shrinkage and selection operator, lncRNA = long noncoding RNAs.

Figure 3.

Ferroptosis-related lncRNA signature risk score analysis. (A–C) Risk score distribution, survival status scatter plots, and heat map of the signature in the training, test, and whole sets. lncRNA = long noncoding RNAs.

Figure 4.

Prediction performances of the signature for BC patients. (A–C) Kaplan-Meier analysis of high and low patients stratified by the median risk score in the training, test and whole sets. (D–F) The time-dependent ROC curves for 5-year RFS predictions by the signature in the training, test, and whole sets. BC = breast cancer, RFS = relapse free survival, ROC = receiver operating characteristic.

Figure 5.

Kaplan–Meier curves were performed for patients stratified by clinicopathological features in the whole set.

3.3. Functional enrichment and microenvironment for signature

GSEA analysis was conducted to find the differences between the high-risk and low-risk groups in the Hallmark pathways. The results revealed that the genes in the high-risk group were enriched in androgen response, MTORC signaling, MYC targets, protein secretion and unfolded protein response. While the pathways in the low-risk group were allograft rejection, apoptosis, P53 pathway, TNFA signaling and WNT beta catenis signaling (Fig. 6). We investigated the relationship between clinicopathological features and the signature of BC patients. As shown in Fig. 7A, 2 groups existed significant differences in 4 characteristics (stage, Her2, PR, and status). Additionally, we obtained mutation distribution and the total TMB for each patient based on the mutation data of the whole set. Interestingly, the names of first 20 genes of highest somatic mutation were mostly similar in the high-risk and low-risk groups, but the high-risk group had higher total mutation frequency than the low-risk group (Fig. 7B,C). In addition, patients in the high-risk group presented higher TMB than the low-risk group (Fig. 7D). We then identified whether TMB was an independent biomarker for BC patients. BC cohorts are split into high-TMB and low-TMB groups as per the median TMB value. According to the findings, the TMB can be used alone for predicting RFS of patients (Fig. 7E). Besides, when the risk score and TMB were combined they also could effectively predict RFS of patients (Fig. 7F). Immune cells are an important constituent of the tumor microenvironment that regulates tumor growth and metastasis.[ Thus, we carried out the immune infiltration and functions for the whole set. The heatmap of the infiltration level of the immune cells estimated by multiple algorithms was shown in Fig. 8A. Immune cell functions including cytolytic activity, human leukocyte antigen (HLA), T cell co-stimulation differed significantly among the 2 groups (Fig. 8B). The expression levels of 27 genes of immune checkpoints, which are crucial for immune immunotherapies, were also explored. The expression of ADORA2A, CTLA4, IDO2, LAG3 among others were considerably lower in the high-risk group compared to the low-risk group (Fig. 8C).

Figure 6.

GSEA analysis based on the median value of the risk scores in high-risk and low-risk groups. GSEA = Gene set enrichment analysis.

Figure 7.

Estimation of the mutation information using ferroptosis-related lncRNAs recurrent signature in the whole set. (A) Heatmap for the signature and clinicopathological characteristics. (B, C) Waterfall plot of mutation profile for BC patients in low-risk and high-risk groups. (D) TMB difference between the high-risk and low-risk patients. (E) Survival curves for the RFS of the high-TMB and low-TMB groups. (F) Survival curves for patients stratified by both TMB and signature. BC = breast cancer, lncRNA = long noncoding RNAs, TMB = tumor mutational burden.

Figure 8.

Immune microenvironment analysis in the whole set. (A) Heatmap for immune cell infiltration based on multiple algorithms between the high-risk and low-risk groups. (B) Differences in immune cell-related function between the high-risk and low-risk groups. (C) Expression of immune checkpoints in the high-risk and low-risk groups.

3.4. Establishment of a nomogram for BC patients

Age, tumor location, T, N, M, stage, ER, PR, and HER2 are important and common clinical pathological factors for BC patients, as well as were also recorded for the great majority of the patients. Therefore, we performed univariate and multivariate Cox regression to identify the factors affecting the RFS based on signature and these clinical variables in the whole set. Univariate Cox regression analysis indicated that age, M, N, T, stage, and risk score were associated with RFS (Fig. 9A). Also, the multivariate Cox analysis indicated that age, M, N, T, and risk score could serve as the independent predictors of RFS (Fig. 9B). Then, a nomogram was established based on these independent predictors for predicting the probability of 1-, 3-, and 5-year RFS. As shown in Fig. 9C, the risk score and stage were the largest contributions to the effect of patients’ RFS among the factors. The range of the C index of the nomogram was (0.63, 082) in different years (Fig. 9D). The calibration curve of the nomogram demonstrated a good agreement between nomogram prediction and actual probability (Fig. 9E). More importantly, The DCA curve showed that the clinical benefit of the nomogram is higher than the extreme curve (Fig. 9F).

Figure 9.

The nomogram model for the prediction of BC prognosis for 1-, 3-, 5-year RFS in the whole set. (A) Nomogram to anticipate recurrent probabilities. (B) The graph showing concordance index changes over time. (C) Calibration curve. (D) The decision curve analysis. BC = breast cancer, RFS = relapse free survival.

4. Discussion

BC is one of the most life-threatening diseases with high morbidity and mortality rate among women worldwide. Recently, significant efforts have been devoted to finding novel treatment targets and strategies through elucidating mechanisms and etiology in order to improve the prognosis of BC.[ Tumor-selective cell death is the therapeutic goal of malignant tumors, and ferroptosis occurring in a variety of tumors which may offer therapeutic advantages.[ Accumulating studies have reported that lots of genes play vital roles in regulating the ferroptosis process in many ways.[ However, there is little to research on the ferroptosis-related lncRNAs for patients with BC. Thus, it is necessary to explore their roles in tumor progress and clinical value for BC based on the large-scale dataset.

In the current study, we firstly screened 77 ferroptosis-related DEGs. Then enrichment analysis was implemented for them. KEGG enrichment analysis showed that the DEGs were involved in response to oxidative and chemical stress, apical plasma membrane, antioxidant activity, HIF-1 signaling pathway and Amebiasis. Disruption of the glutathione-dependent lipid peroxide repair systems will cause ferroptosis through the accumulation of lipid-based reactive oxygen species.[ Heme Oxygenase is an antioxidant that may mediate ferroptosis as new neoadjuvant chemotherapy against tumors.[ We identified 391 ferroptosis-related lncRNAs, 21 lncRNAs of which were found to be associated with the RFS of BC patients. Next, a novel ferroptosis-based lncRNAs signature was established for BC recurrence prediction. Among the lncRNAs in signature, LINC01235 is highly expressed in BC and represents a more aggressive phenotype.[ LINC02166 is also identified as an autophagy-related lncRNA and can help increase the prognostic value of BC.[ Additionally, the contribution of other lncRNAs in BC has not been reported. Subsequently, BC patients were segregated into the high-risk and low-risk groups based on the median risk score, and the high-risk group exhibited a poorer RFS compared with the low-risk group. The results of ROC analysis suggest that the signature had moderate accuracy in predicting RFS. Both univariate and multivariate Cox analyses revealed 5 independent predictor factors for constructing a nomogram, which may quantitatively predict the RFS rate of individual patients in clinic.[ The C-index and calibration curves for different years indicated that nomogram could accurately judge the recurrence of BC patients. Furthermore, DCA cure proved that nomogram could bring a better net benefit to BC patients than either the “treat all” or “treat none” strategies.

Successful treatment of BC is a significant challenge for modern medicine. Decades of basic and clinical research suggest that the immune system plays a vital role in tumor prevention and antitumor.[ In recent years, tumor immunotherapy has established itself as an important adjacent therapy to traditional tumor therapies. Strategies that combine conventional therapies and immunotherapy may provide more optimistic outcomes for patients.[ TMB refers to the total number of somatic mutations in specific regions of a tumor genome,[ and it is helpful for predicting sensitivity to immunotherapy and prognosis in BC.[ This study reported that the TMB in the high-risk group exceeded that in the low-risk group. Besides, patients stratified by both TMB and signature has different outcomes. Likewise, immune checkpoint and tumor infiltrating immune cells are both the factors significantly associated with the immunotherapy effect.[ We found that immune cell functions were significantly stronger in the low-risk group than those in the high-risk group, especially in cytolytic activity, HLA. Out of 27 immune checkpoints, the great majority of genes were significantly highly expressed in the low-risk group. The above results suggest that the signature-based grouping may provide inferences for doctors when selecting immunotherapy options.

In the clinical, pathological stage is considered as the most important factor for tumor prognostic prediction. However, many patients at the same stage have different clinical processes and outcomes, indicating that tumor heterogeneity is influenced by multi-faceted aspects.[ Hence, it is necessary to explore novel predictive and therapeutic biomarkers and their functions. Constructing a ferroptosis-related lncRNAs signature may offer a new method for personalized recurrent prediction. The strength of this study is that the signature is only composed of 5 lncRNAs; so it is more probably to be applied in clinical practice. The results about function enrichment and microenvironment also provide a new direction for future studies on the roles of lncRNAs in ferroptosis and tumor immunotherapy. Meanwhile, some drawbacks should not be ignored in the present study. Firstly, some clinical information downloaded from the TCGA database is missing, such as operation methods, which may cause bias and errors. Secondly, only internal validation was performed from a single cohort in this study but lacked independent data validation and prospective studies to verify its clinical utility. Lastly, more experiments in vitro and vivo are needed to further elucidate the regulatory mechanism of the ferroptosis-related lncRNAs in BC. For the above disadvantages, our research group will focus on the issues and refine the mechanism of lncRNAs involvement in ferroptosis and the clinical application value of the signature.

5. Conclusions

In summary, our study identified a 5 ferroptosis-related lncRNAs signature for predicting RFS of patients with BC based on the public database. This signature may assist us with discovering novel biomarkers for clinical management. In addition, the signature also may predict the treatment response to immunotherapy and provide therapeutic targets for BC patients.

Acknowledgment

We thank all the researchers involved in the consolidation and submission of the data from the TCGA database, which may provide convenience and possibility of tumors studies in a large cohort.

Author contributions

All authors read and approved the final manuscript.

Conceptualization: Yuzhi Wang, Yi Zhang.

Data curation: Yi Zhang, Yunfei Xu.

Formal analysis: Yuzhi Wang, Yi Zhang.

Funding acquisition: Yi Zhang.

Investigation: Yi Zhang, Yuzhi Wang.

Methodology: Yuzhi Wang, Yunfei Xu.

Project administration: Yuzhi Wang, Yi Zhang.

Resources: Yuzhi Wang.

Software: Yi Zhang, Yunfei Xu.

Supervision: Yuzhi Wang, Yunfei Xu.

Validation: Yuzhi Wang, Yunfei Xu.

Visualization: Yi Zhang.

Writing – original draft: Yuzhi Wang.

Writing – review & editing: Yi Zhang, Yunfei Xu.