The gene expression data of the transcriptome (such as, mRNA and LncRNA) and the clinical information of BC patients were downloaded from the TCGA dataset (19), of which the clinical information included the survival time, survival status, age, gender, grade, stage, T stage, N stage, and M stage. Specifically, we preliminarily screened 1,178 cases of transcriptome data, such as 112 cases of the normal and 1,066 cases of the tumor, and 1,053 cases of clinical data that include 911 cases of survival and 142 cases of death. A list of ferroptosis-related genes was extracted from FerrDb (http://www.zhounan.org/ferrdb/operations/download.html), and the expression of ferroptosis-related genes was extracted.

For picking out the target ferroptosis-related LncRNAs (FRlncRNAs), a co-expression network analysis was adopted to show the relationship between lncRNAs and ferroptosis-related genes. The LncRNAs with |r| ≥ 0.5 and P < 0.001 were confirmed as FRlncRNAs. Then, the univariate Cox regression analysis was used to screen the prognosis-related FRlncRNAs, and the results were presented in the form of forest plot. Further, we drew a heatmap and compared the differential expression of these FRlncRNAs in normal tissues and tumor tissues using the rank sum test.

According to the FRlncRNAs related to prognosis, the hierarchical consensus clustering was used to perform the classification of TCGA cohort (20). To obtain robust classification, we adopted an unsupervised consensus approach implemented in the R package “Consensus Cluster Plus” (21). Moreover, the relative change in area under the cumulative distribution function (CDF) curve was employed to determine the optimal number of clusters, k, which was further verified by the total within sum of squares (WSS) and the gap statistics. The difference of survival probability and clinical information (age, stage, T stage, N stage, and M stage) between clusters were investigated.

Inspired by the success of immunotherapy in the patients with BC in recent years (22), we further explored whether there was an immunological explanation for the survival differences between the clusters. Set PD-L1 gene as the strongest example, we carried out Pearson's correlation coefficient to test the co-expression and correlation between the hub gene and the prognostic FRlncRNAs in tumor tissues. Besides, ESTIMATE algorithm was performed to calculate the immune and stromal scores to quantify the presence of stromal cells and the infiltration of immune cells in tumor samples. To observe the differences of immune cells among clusters in the tumor microenvironment in a more detailed way, CIBERSORT, a gene expression-based deconvolution algorithm to describe the cell constitution of tissues (23), was performed to intuitively display the distribution of immune cells, which was showed in the violin plot generated by the vioplot package.

Gene Set Enrichment Analysis determines whether the predetermined gene sets have statistically significant differences between the two biological states in a computational method (24). In view of the consensus clustering, we conducted GSEA analysis in clusters with the aim of mining survival differences. By adjusting the p-value, the enrichment pathway for each phenotype was classified via the normalized enrichment score (NES). An NES >1 and false discovery rate (FDR) < 0.05 denoted statistical significance.

Based on the results of univariate Cox regression analysis, the least absolute shrinkage and selection operator (LASSO) regression was applied to select the optimal survival-related FRlncRNAs, which were involved in the final modeling. According to the coefficients derived from LASSO regression and expression levels of FRlncRNAs included in final models, the prognostic risk score formula was constructed as follows:

where β is the regression coefficient.

We randomly divided the patients with BC from TCGA into two groups according to the ratio of 1:1, one group as the training set and the other group as the validation set. In the training set, we calculated the risk scores of the patients with BC and classified these patients into the high-risk group and the low-risk group on the basis of the median risk score as the threshold. The Kaplan–Meier (K–M) survival analysis was used to prove whether there was a survival difference between the two groups. The receiver operator characteristic (ROC) curve was built by using the survival ROC package to assess the efficiency of the prognostic model. In the test group, the same processes were performed to validate the prognostic model of this group. Independent cohort validation is important for prognostic signatures. In the current study, the GSE69031 cohort was used to validate our OS signature (25). The expression data of the genes included in the final signature were obtained and substituted into the equation for risk score calculation. All patients in this cohort were stratified into low- or high-risk groups. The prediction accuracy of signature in the independent validation cohort was evaluated by ROC curve and K–M survival analyses.

The clinical information, such as age and stage, was integrated to testify the independence of the prognostic model with the combination of risk score. For that, the univariate and multivariate Cox regression analyses were performed to verify the independence of the prognostic model in the training cohort and test cohort, respectively.

Through stratification of clinical data, K–M curves were drawn, respectively to study whether the genetic risk scores were applicable to patients in different groups. Specifically, we divided patients into the two categories based on age >65 and age <65, T1-2 and T3-4, N0 and N1–3, M0 and M1, Stage I–II and Stage III–IV, respectively, and calculate the difference in survival curves between the high-risk and low-risk patients in each category, so as to expand the applicability of risk scores.

We enrolled all eligible patients and divided them into the high-risk and low-risk groups according to the calculation above. The clinical information (age, stage, T stage, N stage, and M stage), immune score, and clusters were integrated to exhibit the differences between high-risk and low-risk groups. More importantly, we also considered the expression differences of the survival-related FRlncRNAs involved in the final modeling, hoping to provide a theoretical basis for finding new therapeutic targets for BC.

Through the co-expression analysis, we sorted out the 63 FRLncRNAs, and visualized the relationship by co-expression network (Figure 1). In the Figure 1, the red nodes represent the ferroptosis-related genes, and the blue nodes represent LncRNAs co-expressed with ferroptosis-related genes. These LncRNAs were used as FRlncRNAs for the subsequent analyses.

To explore the influence of FRlncRNAs on the prognosis of patients with BC, the univariate Cox regression analysis was used to preliminarily determine 31 FRlncRNAs related to the prognosis, which were visualized by forest map (Figure 2). If the hazard ratio (HR) >1, the higher the expression level, the higher the risk of patients. On the contrary, HR < 1 means that the higher the expression level, the lower the risk of patients. Further, the differential analysis was employed to exhibit the expression level of 31 prognosis-related FRlncRNAs between the normal tissues and tumor tissues. LINC00702, AC121247.2, HSD11B1-AS1, NR4A1AS, AC011472.4, LIPE-AS1, etc., were lower expressed in tumor tissues, conversely, AP001434.1, LINC02257, LINC01655, LINC01614, AF015262.1, LMNTD2-AS1, etc., were highly expressed in tumor tissues (Figure 3).

Our study confirmed that some FRlncRNAs were related to the prognosis of patients, and were expressed in a heterogeneous manner among tumor patients. In order to better understand the intertumoral heterogeneity of BC, we conducted an unsupervised consensus analysis to explore the influence of FRlncRNAs on the occurrence and development of BC from multiple perspectives. Using the similarities in the expression of prognosis-related FRlncRNAs, we choose the value of k = 2 as optimal selection (Figures 4A,B). Consequently, the two clusters of samples were determined as follow: C1 (n = 726, 80.8%) and C2 (n = 172, 19.2%). We then applied K–M curves to compare the differences in survival between different clusters and found that C1 tended to carry a good prognosis (Figure 4C). Next, we researched whether there was a correlation between the clinical data and clusters (Figure 4D). It was obvious that the expressions of FRlncRNAs in the upper right corner of heatmap were upregulated in C2, with a more red color, such as AP001434.1, LINC02257, LINC01655, LINC01614, and AF015262.1. More interestingly, we noticed that most of these FRlncRNAs upregulated in C2 with poor prognosis were also upregulated in tumor tissue, suggesting that these FRlncRNAs may act as tumor promoters to accelerate tumor migration and progression. At the same time, we could also compare whether there were differences in clinical traits in different clusters. Figure 4D shows no significance in the clinical parameters indicating no differences between these characteristics in the various clusters.

To investigate the correlation of FRlncRNAs with immune features, we conducted further analysis from the perspective of immunology. The correlation analysis was performed to explore the relationship between prognosis-related FRlncRNAs and PD-L1 gene (Figure 5). The result implies that PD-L1 was significantly associated with some FRlncRNAs, such as LINC02328, LINC01871, LMNTD2-AS1, LINC00702, LINC02384, AP001434.1, TRG-AS1, LIPE-AS1, HSD11B1-AS1, USP30-AS1, LINC02446, AC099329.2, PRKAR1B-AS1, and LINC02084. Furthermore, we investigated the differences of immune microenvironment between the clusters. Intriguingly, we found that there was no difference in immune cells between C1 and C2 (Figure 6A), whereas C2 was associated with higher stromal scores compared with C1 (Figure 6B). Furthermore, violin plot analysis showed that the levels of cell infiltration (Figure 6C), such as B cells naïve, plasma cells, T cells CD8, T cells CD4 memory activated, T cells follicular helper, natural killer (NK) cells activated, monocytes, macrophages M1, dendritic cells resting, and neutrophils were higher in C1 than in C2, whereas the levels of macrophages M0 and macrophages M2 were lower in C1 than in C2. Collectively, the turbulent changes of immune cells in the tumor immune microenvironment may support the conclusion that C2 had a poor prognosis to some extent.

The explanation from the immune point of view made us understand the impact of target FRlncRNAs on tumor microenvironment more thoroughly. We further conducted GSEA analysis to explore the biological signal pathways with obvious differences between the clusters. Some metabolism-related pathways, in relation to “adipocytokine signaling pathway,” “linolenic acid metabolism,” “arachidonic acid metabolism,” “fatty acid metabolism,” “ether lipid metabolism,” and “glutathione metabolism” were significantly enriched in C1, among which the pathways related to lipid metabolism accounted for the majority (Figures 7A–F). In addition, we found that immune-associated pathway was enriched, such as “T cell receptor signaling pathway” and “TGF-beta signaling pathway” (Figures 7G,H). Besides, the enrichment of some carcinogenic pathways was also exposed, such as “pathway in cancer,” “apoptosis,” “cytosolic DNA sensing pathway,” and “base excision repair” (Figures 7I–L).

First, according to the ratio of 1:1, the patients with BC were randomly divided into a training set and a validation set, both of which included 448 patients with complete information, as described in the previous literature (26). Second, according to 31 prognostic FRlncRNAs calculated by univariate Cox regression analysis, we then performed LASSO regression analysis to pick out the optimal prognosis-related FRlncRNAs with nonzero coefficients (Figures 8A,B). Consequently, 12-FRLncRNA signature was determined. The formula of the final model was listed as follows: risk score = −0.250^*expression of LINC01871+0.091^*expression of LINC00393+1.058^* expression of AC121247.2-0.330^* expression of LINC02384-0.336^* expression of LIPE-AS1-0.072^* expression of HSD11B1-AS1+0.361^* expression of AC010655.2+0.018^* expression of LINC01419+0.051^* expression of PTPRD-AS1+0.137^* expression of AC099329.2+0.024^* expression of OTUD6B-AS1+0.157^* expression of LINC02266. We then calculated the risk scores and divided patients into high- and low-risk groups by the median risk score in the training set (Figure 9A). The relationships between the survival status and survival times of patients with BC ranked by risk scores were depicted in Figure 9B. In addition, a heatmap was plotted to show the expression profiles of 12 FRLncRNAs (Figure 9C). To study the relationship between the risk score and survival probability, K–M curves were carried out in Figure 9D. Patients in the high-risk group had a worse survival probability, whereas those in the low-risk group had a better survival probability. Furthermore, the ROC curve was plotted to verify the predictive ability of the models, whose AUC of the ROC curves was 0.693, revealing a positive prognostic ability (Figure 9E).

To verify the results of the training set, we employed the same models on the patients in the testing set (Figure 10). The results showed that the high-risk patients showed significantly worse survival probability than the low-risk groups (Figure 10D), whose AUC of ROC was 0.655 (Figure 10E), which suggested that the prognostic model could satisfactorily predict the prognosis of patients with BC. In addition, the prognostic values of the risk score of patients with BC in the GSE69031 cohort were calculated. The results indicated that patients in the high-risk group showed a worse prognosis than the low-risk patients (Figure 10F). The AUC values of signature to predict the OS was 0.681 (Figure 10G). Generally, our signature showed satisfactory performance in the independent cohorts, which indicated that these signatures are robust prognostic biomarkers.

In view of the complexity of a variety of clinical factors, further univariate and multivariate Cox regression analyses were performed to explore the independence of risk scores (Figures 11A,B). Collectively, the risk score was an independent prognostic factor, independent of other clinical factors, in both the training and validation sets (P < 0.05). Generally, our research indicated that the patients with higher risk scores were a worse prognosis.

First of all, we divided patients into the two categories based on age >65 and age <65, T1-2 and T3-4, N0 and N1-3, M0 and M1, Stage I-II, and Stage III-IV. By employing the same models, we plotted K–M survival curves for each subgroup (Figure 12). For example, in Figures 12C,D, patients with T stage at T1 and T2 showed statistically significant differences in their survival curves, and patients with T3 and T4 also showed statistically significant differences in their survival curves, indicating that the risk score was applicable to the patients with different T stages (P < 0.05). Similarly, in other grouping variables, such as age <65 and age> 65 (Figures 12A,B), in N0 and N1–3 patients (Figures 12E,F), in early and late-stage patients (Figures 12G,H), value of P < 0.05, indicating that the risk score of this gene was applicable to the different groups of patients. Nevertheless, among M0 and M1 stage patients, the difference was not statistically significant due to the small number of M1 stage patients, which required further data to verify. Overall, in the subgroup analyses, we confirmed that the high-risk patients were a worse prognosis than the low-risk patients in all subgroups (all p < 0.05).

Through a comprehensive analysis of the heatmap, we can easily distinguish the high-risk and low-risk FRLncRNAs (Figure 13). Specifically, LINC01871, LINC02384, LIPE-AS1, and HSD11B1-AS1 were highly expressed in low-risk groups, demonstrating that these FRLncRNAs were low-risk LncRNAs. On the contrary, LINC00393, AC121247.2, AC010655.2, LINC01419, PTPRD-AS1, AC099329.2, OTUD6B-AS1, and LINC02266 were highly expressed in the high-risk groups, demonstrating that these FRLncRNAs were high-risk LncRNAs. At the same time, we could also compare whether the clinical traits were different between the high- and low-risk groups. It could be seen that there were differences among the N stages, immune scores, and clusters (P < 0.05). More specifically, patients in the low-risk groups were more likely to be assigned to C1 and had a higher immune score, which were consistent with a better prognosis.