Metabolism-related long non-coding RNAs (lncRNAs) as potential biomarkers for predicting risk of recurrence in breast cancer patients.

Metabolism affects the development, progression, and prognosis of various cancers, including breast cancer (BC). Our aim was to develop a metabolism-related long non-coding RNA (lncRNA) signature to assess the prognosis of BC patients in order to optimize treatment. Metabolism-related genes between breast tumors and normal tissues were screened out, and Pearson correlation analysis was used to investigate metabolism-related lncRNAs. In total, five metabolism-related lncRNAs were enrolled to establish prognostic signatures. Kaplan-Meier plots and the receiver operating characteristic (ROC) curves demonstrated good performance in both training and validation groups. Further analysis demonstrated that the signature was an independent prognostic factor for BC. A nomogram incorporating risk score and tumor stage was then constructed to evaluate the 3 - and 5-year recurrence-free survival (RFS) in patients with BC. In conclusion, this study identified a metabolism-related lncRNA signature that can predict RFS of BC patients and established a prognostic nomogram that helps guide the individualized treatment of patients at different risks.

Introduction

Breast cancer (BC) is the most common malignancy in women and leads to cancer death among females worldwide [1]. It is a heterogeneous disease which is mainly manifested in biological behavior, gene expression profiles, clinical prognosis, and therapeutic response, even at the same tumor stage [2]. Although the overall survival rate of BC patients has greatly improved over the past decades, frequent drug resistance, recurrence, and metastasis are still the leading obstacles to current clinical treatment [3]. However, currently used biomarkers have limited clinical value in predicting prognosis and metastatic recurrence. These features highlight the urgent need for the development and validation of novel prognostic signatures to accurately predict clinical outcomes of BC patients and provide guidance for personalized treatment.

Long non-coding RNAs (lncRNAs) represent a class of new and potent tumor regulator RNAs, with > 200 nucleotides and low protein-coding potential [4]. Emerging evidence has suggested that lncRNAs can directly combine with DNA, RNA, or proteins to regulate gene expression in the form of RNA at various levels, thus resulting in the alteration of physiological and/or pathological processes, including cell cycle, differentiation, proliferation, apoptosis, transcriptional regulation, and carcinogenesis by modulating gene expression at the post-transcriptional level [5,6]. The competing endogenous RNA (ceRNA) hypothesis is based on a large-scale regulatory network system representing the complex cross-talk between coding and non-coding RNAs, and their interactions can be predicted using several tools [7,8]. lncRNAs have been proven to be oncogenic factors or tumor suppressors in the development of BC. For example, lncRNA BCRT1 promotes BC progression by targeting the miR-1303/PTBP3 axis [9]. In contrast, lncRNA FGF14-AS2 has been identified as a tumor suppressor that suppresses BC metastasis by regulating the miR-370-3p/FGF14 axis [10].

Energy metabolism reprogramming, which can promote rapid cell growth and proliferation, is an emerging hallmark of cancer [11]. Several studies have reported that oncogenes and tumor suppressor genes such as PI3K [12], MYC [13], K-Ras [14], BRCA1 [15], and p53 [16] can regulate metabolic pathways at multiple levels in different cellular compartments. For example, rewiring of metabolism induced by oncogenic K-Ras in cancer cells involves both glucose and glutamine utilization and sustains enhanced, unrestricted growth. Combination therapy with the glutaminase inhibitor CB-839 and PI3K/aldolase inhibitor NVP-BKM120 effectively reduce cell growth in tumor xenografts [14]. Metabolic alterations are a feature of both treatment-naive and treatment-resistant BC. Therefore, understanding the mechanisms involved will allow us to identify potential tumor vulnerabilities for new therapeutic approaches. Accumulating evidence has confirmed the existence of crosstalk between lncRNAs and cellular metabolism in several types of cancer. Liu et al. [17] revealed that lncRNA AGPG can regulate PFKFB3-mediated tumor glycolytic reprogramming. Another recent study [18] indicated that lncRNA GLS-AS impairs GLS-mediated metabolism and inhibits malignant progression in pancreatic cancer. However, the actions of lncRNA-mediated metabolism in BC are not fully understood. Given the potential influence of lncRNAs and metabolism on the occurrence and development of BC, we sought to identify novel metabolism-related lncRNAs and construct prognostic signatures that are related to recurrence-free survival (RFS) in BC patients. This risk model will assist in understanding the molecular mechanisms underpinning BC and provide new ideas for metabolism-targeting therapies. First, metabolism-related lncRNAs were identified using Pearson correlation analysis. We established a metabolism-related lncRNA signature and genomic-clinicopathologic nomogram which can accurately evaluate individual prognostic markers in patients with BC. In summary, the current investigation may enable great progress in the exploration of prognostic metabolism-related lncRNAs and shed new insights into the possible mechanisms involved in BC development.

Materials and methods

Data sources and preprocessing

Data from The Cancer Genome Atlas (TCGA) and GTEx of breast cancer were acquired from the UCSC Xena Database (http://xena.ucsc.edu/), including 1,104 tumor samples and 292 normal samples. Samples with incomplete clinical data or a follow-up time of < 30 days were excluded. After the removal of patients lacking recurrence data, only 318 BC patients remained in TCGA cohort. Gene expression and clinical data of 241 BC patient samples from the Gene Expression Omnibus (GEO) database (GSE21653, https://www.ncbi.nlm.nih.gov/geo/) were used as a validation set. Metabolism-related gene sets were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways database, which included 186 gene sets from the GSEA website (http://software.broadinstitute.org/gsea/downloads.jsp#msigdb). According to the criteria described in previously published articles (|Pearson R| > 0.4, p-value < 0.05), Pearson correlation analysis was used to identify metabolism-related lncRNAs [19,20].

Construction and validation of lncRNAs prognostic signature

Univariate Cox proportional hazards regression analysis was performed to screen metabolism-related lncRNAs significantly associated with patient RFS in the training set. A p-value < 0.05 was set as the cutoff criterion. To minimize the risk of over-fitting and to remove highly related genes, LASSO Cox regression was performed using the glmnet and survival package in R software [21]. Multivariate Cox analysis was subsequently used to construct a prognostic signature. Risk scores were obtained based on a linear combination of the expression levels and weighted regression coefficients obtained by multivariate Cox analysis. Risk score = expression of gene1 × β1 + expression of gene2 × β2 +⋯+ expression of genen × βn. where β represents the coefficient value. All samples were grouped into low – and high-risk groups using the Survminer package in R software. Log-rank tests were performed to compare the prognostic outcomes between the two subgroups. Furthermore, receiver operating characteristic (ROC) analysis was applied to evaluate the predictive capability of the above signature using the SurvivalROC R package [22]. Then, we performed principal component analysis (PCA) using the ‘prcomp’ function of the ‘stats’ R package.

Validation of metabolic signature

To validate the homogeneity of the established metabolism-related lncRNA signature, the GSE84437 dataset were used as a validation set. BC samples were assigned to low – and low-risk groups using the formula described above. Survival analysis and log-rank tests were performed to compare RFS between the two risk groups. The predictive ability of the prognostic signature was assessed using the time-dependent ROC curve.

Development and validation of a prognostic nomogram

To quantitatively estimate BC prognosis in clinical practice, a prognostic nomogram that integrated both the metabolism-related lncRNA signature and clinicopathological characteristics was generated based on the multivariate Cox regression model. The predictive ability of the nomogram was evaluated using calibration plots and the ROC curve for 3 – and 5-year RFS.

Gene Set Enrichment Analysis (GSEA)

To explore prognostic biomarkers involved in potential biological processes, GSEA was performed on the metabolism-related prognostic signature in the low – and high-risk groups. The classical gene set in the Molecular Signatures Database (MSigDB), ‘c2.cp.kegg.v7.0.symbols.gmt (curated)’ was considered. FDR < 0.25, and p-values < 0.05 were regarded as reference values.

Statistical analyses

All statistical analyses were performed using R v.4.0.4. Cox proportional hazard regression analysis was used to determine whether the metabolism-related lncRNA signature could be an independent prognostic indicator. Stratification analysis was conducted to further validate the prognostic performance of the metabolism-related lncRNA signature. Statistical significance was set at p < 0.05, unless otherwise noted. The nomogram created by the ‘Survival’ and ‘RMS’ packages of R v/3.6.3 was used to provide a visual risk prediction.

Results

Following in-depth studies of metabolic reprogramming, researchers have gradually realized the importance of fatty acid metabolism in BC. However, understanding of the association between lncRNAs and metabolism in BC is limited. In this study, we developed and validated a novel metabolism-related lncRNA signature and a nomogram integrating the signature and TNM stage for predicting individual survival among BC patients via the integrated analysis of TCGA and GEO datasets. Furthermore, using this signature and TNM staging, a nomogram was constructed to predict the 3 – and 5-year RFS of patients with BC. Hence, our study developed a metabolism-related lncRNA prognostic signature that could facilitate prognostic stratification in patients with BC and provide new ideas for metabolism-targeting therapies.

Identification of metabolism-related lncRNAs in BC

Figure 1 presents a flowchart of the data analysis steps involved in this study. First, we identified 13,361 lncRNAs and 19,318 mRNAs from the TCGA dataset, and obtained 945 metabolism-related genes from the gene set database of KEGG pathways. Pearson correlation analysis was performed on the expression values of these lncRNAs and metabolism-related genes to screen-out highly metabolism-related lncRNAs with |Pearson R| > 0.4, and p-values < 0.05. As a result, 1,887 metabolism-related lncRNAs were identified in the TCGA dataset. Consequently, 241 metabolism-related lncRNAs were selected by taking the intersection of the TCGA and GEO datasets (Supplementary Table 1).

Figure 1.

Flow chart illustrating development of the metabolism-related lncRNA signature identified in this study

Table 1.

Coefficients and multivariable Cox model results in breast cancer

LncRNA	Coef	HR	95% CI	P value
FOXD2.AS1	0.402925894	1.496	1.013–2.209	0.042
A1BG.AS1	−0.454896016	0.635	0.390–1.031	0.066
C9orf163	1.002781719	2.726	1.392–5.336	0.003
GSN.AS1	0.744058391	2.104	1.302–3.403	0.002
LINC00893	−1.386910573	0.250	0.097–0.641	0.004

Development and validation of lncRNAs prognostic signature

We performed univariate Cox regression, LASSO regression, and multivariate Cox regression analyses of 241 metabolism-related lncRNAs. A total of 12 prognostic genes related to RFS were selected from the metabolism-related lncRNAs using univariate Cox regression analyses (Figure 2a). Next, we conducted LASSO regression analysis, which identified 10 crucial lncRNAs (Figure 2b and c). Then, multivariate Cox regression analysis was conducted, and five lncRNAs were finally selected to establish a prognostic signature (Table 1). Among the five lncRNAs, FOXD2-AS1, C9orf163, and GSN-AS1 were risk prognostic genes with HR > 1, while A1BG-AS1 and LINC00893 were protective genes with HR < 1 (Figure 2d). Then, a 5-lncRNA metabolic signature was built using the adjusted regression coefficients of each gene and associated risk scores. Risk scores (RSs) were calculated using the following formula: RS = (0.40292 × level of FOXD2-AS1) + (1.00278 × level of C9orf163) + (0.74405 × level of GSN-AS1) + (− −0.45489 × level of A1BG-AS1) + (−1.38691 × level of LINC00893). After scoring each patient’s risk through the signature, we divided the BC patients into low – and high-risk groups. The risk score distribution and survival status of the signature in the training set showed that the signature had a significant value for evaluating the prognosis of patients with BC (Figure 2e and f). PCA analysis indicated discernible dimensions between the low – and high-risk patients (Figure 2g). Kaplan-Meier analysis demonstrated that patients in the high-risk group exhibited significantly poorer RFS than those in the low-risk cohort (Figure 2h). To further assess the model performance, the time-dependent ROC was mapped. The AUC values were approximately 0.735 for the training set (Figure 2i).

Figure 2.

Development and evaluation of metabolism-related lncRNA signature related to RFS in TCGA cohort. (a) Forest plot of 12 candidate metabolism-related lncRNAs selected by univariate Cox regression analysis. (b) LASSO coefficient profiles of the 10 candidates. (c) Tenfold cross-validation for tuning parameter selection in the LASSO model. (d) Forest plot of 5 candidate metabolism-related lncRNAs selected by multivariate Cox regression analysis. (e) Patient survival status distribution by RS. (f) Patient survival status distribution of the low-risk group and the high-risk group. (g) PCA plot. (h) Kaplan-Meier method was used to plot the RFS curve for the high RS and low RS groups. (i) ROC curve of 5 lncRNA signatures

Validation of the metabolism-related signature

The GSE21653 dataset was separately used to determine the validity and robustness of the signature in an independent validation cohort (Figure 3a-e). We calculated the RS of the metabolism-related lncRNA prognostic signature for each patient in the validation set using the same formula. As expected, PCA analysis demonstrated that patients in the different risk groups were distributed in two directions (Figure 3c). Survival analysis indicated that the results of the Kaplan-Meier survival curve analysis were statistically significant (p < 0.05, Figure 3d). Additionally, the AUC for RFS was 0.695 in the validation cohort, suggesting a favorable performance of the signature (Figure 3e).

Figure 3.

Validation of metabolism-related lncRNA signature related to RFS in GEO cohort. (a) Patient survival status distribution by RS. (b) Patient survival status distribution of the low RS and the high RS. (c) PCA plot. (d) Kaplan-Meier method was used to plot the RFS curve for the high RS and low RS groups. (e) ROC curve of 5 lncRNAs signatures

Correlation between signature and clinicopathological characteristics

All patients were grouped according to their RS values to analyze any association between the metabolism-related signature and clinicopathological characteristics. The results revealed that there were distinct differences in the RSs between patients involving tumor stage, surgery type, and tumor status (p = 0.009, 0.04, < 0.001, respectively; Supplementary Figure 1). However, age (p = 0.43) and margin status (p = 0.33) were not significantly correlated with the RS.

Development of a nomogram for prognostic prediction

To determine whether the prognostic significance of the signature is dependent on clinicopathological parameters, univariate and multivariate Cox analyses were performed to analyze the following variables: RS, age, surgery type, margin status, and tumor stage (Figure 4a and b). As shown in Figure 4a, of the univariate analysis, tumor stage and RS were significantly correlated with RFS of BC patients as well as the results of multivariate analysis (Figure 4b). To evaluate whether the metabolism-related lncRNA prognostic model could serve as a prognostic indicator for RFS in subgroups of patients with different clinical characteristics, we stratified subgroups by age (age ≤ 60 years and age > 60 years), pathological stage (stage I–II and stage III–IV), surgery type (lumpectomy, modified radical mastectomy, and others), and tumor status (tumor-free and tumor-positive). The results indicated that the RFS in low-risk patients based on age (p-value = 0.001 in ≤ 60 and p-value = 0.004 in > 60), pathological stage (p = 0.003 in stage I–II and p-value = 0.005 in stage III–IV), surgery type (p-value = 0.035 in lumpectomy, p-value = 0.041 in modified radical mastectomy, and p-value = 0.002 in other), and tumor status (p-value = 0.020 in tumor-free and p-value = 0.032 in with tumor) was significantly higher than those of high-risk patients (Supplementary Figure 2).

Figure 4.

Identification of RFs and development of prognostic nomogram. (a) Univariate Cox regression analysis. (b) Multivariate Cox regression analysis. (c) Prognostic nomogram incorporating RFs for predicting the probability of 3 – and 5-year RFS in BC patients. (d) Time-dependent ROC for 3 – and 5-year RFS predictions of the nomogram. (e, f) Calibration curves for predicting 3 – and 5-year RFS. The nomogram-predicted probability of survival is plotted on the x-axis; actual survival is plotted on the y-axis

To provide surgeons with a quantitative tool for prognostic prediction in BC patients, we constructed a nomogram that integrated the RS based on the expression levels of five metabolism-related lncRNAs and tumor staging (Figure 4c). With the help of the nomogram, prognosis can be effectively predicted based on the individual characteristics of the patient. The nomogram was validated using ROC and calibration plots. The 3 – and 5-year AUC values of the nomogram for RFS were 0.774 and 0.745, respectively (Figure 4d). In addition, the calibration plots for the 3 – and 5-year RFS showed strong consistency between nomogram predictions and actual observations (Figure 4e and f).

GSEA

According to the criteria of FDR < 0.25, and p < 0.05, 18 significantly altered pathways were enriched in the high-risk group of the signature, including the P53 signaling pathway, DNA replication, galactose metabolism, fructose and mannose metabolism, amino sugar and nucleotide sugar metabolism, pyrimidine metabolism, and other pathways (Figure 5).

Figure 5.

GSEA analysis of differentially expressed genes in high – and low-risk groups. GSEA implied remarkable enrichment of metabolism-related phenotypes in the high-risk group

Discussion

Accurate prognostic prediction and individualized clinical treatment strategies are the basis of precision medicine [23]. Most of the established clinical markers for treatment response and prognosis in BC are based on clinical features, and their accuracy and specificity are limited. Traditionally, AJCC TNM staging is mainly based on anatomical information and cannot adequately assess recurrence in patients with BC. Thus, it is of great urgency to explore the molecular mechanisms and to screen specific cancer-related biomarkers for risk evaluation that can be employed to predict the survival of BC patients and promote the development of effective therapies for BC. Following in-depth studies of metabolic reprogramming, researchers have gradually realized the importance of metabolism in BC [24,25]. However, metabolism-associated lncRNAs for predicting personalized survival are lacking. Given the potential influence of lncRNAs and metabolism on the occurrence and development of BC, the identification of novel metabolism-associated biomarkers might provide important prognostic information and novel therapeutic targets.

A recent study based on public databases established a tumor microenvironment-related signature model to predict the prognostic features of BC [26]. However, metabolism-associated lncRNA signatures are lacking. In this study, we constructed and validated a novel metabolism-related lncRNA signature and a prognostic nomogram incorporating RSs and TNM staging for predicting individual survival among BC patients via the integrated analysis of TCGA and GEO datasets. First, we identified metabolism-related lncRNAs using Pearson’s correlation analysis. Univariate, LASSO, and multivariate Cox regression analyses were used to build a metabolism-related prognostic signature, which could be an independent prognostic factor for BC patients. This metabolism-related lncRNA prognostic model could serve as a prognostic indicator for RFS in subgroups of patients with different clinical characteristics. Patients with a low RS had superior RFS than those with a high RS. The predictability of the signature was confirmed by the tdROC curve, PCA analysis, and an independent GEO dataset. The GSEA results revealed that these metabolism-related lncRNAs are involved in many key biological functions related to tumor metabolic pathways. Furthermore, using this signature and TNM stage, we established a nomogram for predicting the 3 – and 5-year RFS of BC patients. The predictive ability of the model was validated by the discrimination and calibration curves. Hence, this study developed a metabolism-related lncRNA prognostic signature that could facilitate prognosis stratification in patients with BC and provide new ideas for metabolism-targeting therapies.

Accumulating evidence shows that lncRNAs also play irreplaceable roles in regulating metabolism in BC through different mechanisms. Ma et al. [27] found that lncRNA FGF13-AS1 can inhibit glycolysis and the ‘stemness’ properties of BC cells. Mechanistically, FGF13-AS1 can reduce the half-life of MYC mRNA by binding to IGF2BPS and disrupting the interaction between IGF2BPS and MYC mRNA. In addition, MYC transcriptionally inhibits FGF13-AS1, forming an FGF13-AS1/IGF2BPs/Myc feedback loop in this signaling pathway. Du et al. [28] revealed that lncRNA MIR210HG can promote the Warburg effect and tumor growth by enhancing HIF-1α translation in triple-negative BC. Li et al. [29] demonstrated that lncRNA SNHG3 actively regulates PKM expression by acting as a miR-330-5p sponge, inhibiting mitochondrial oxidative phosphorylation, increasing glycolytic carboxylation, and enhancing the proliferation of BC cells. Surprisingly, in addition to FGF13-AS1, MIR210HG, and SNHG3 mentioned above, many of these lncRNAs have been shown to play important roles in metabolism [30,31].

Among the 5 lncRNA signatures identified in this study, 4 lncRNAs have been previously reported. A1BG-AS1 and C9orf163 have been shown to be immune-related lncRNAs and are associated with prognosis in patients with BC and pancreatic cancer, respectively, [32,33]. Robust evidence has shown that there is an intimate relationship between immunity and metabolism in cancer [34,35]. We speculate that A1BG-AS1 and C9orf163 may affect immunity in BC by regulating metabolism, although further experiments are required to verify such speculation. In addition, Bai et al. [36] revealed that A1BG-AS1 functions as a tumor suppressor in hepatocellular carcinoma (HCC) and may promote tumorigenesis and progression through ceRNA-related mechanisms in HCC cells. Metabolic pathways such as glycolysis have distinct, essential roles in tumorigenicity and cancer progression [14]. Wang et al. [37] indicated by Pearson’s correlation analysis that FOXD2-AS1 is associated with glycolysis-related lncRNAs. However, the authors did not explore how FOXD2-AS1 affects the development and progression of BC by regulating metabolism. FOXD2-AS1 is also reported to be upregulated in tumor tissues and associated with poor prognosis in BC cases [38,39]. Knockdown of FOXD2-AS1 can inhibit the proliferation, migration, and invasiveness of BC cells and decrease the growth of transplanted tumors in vivo. FOXD2-AS1 promotes malignancy and tumorigenesis through different mechanisms, including the S100 calcium binding protein A1/Hippo signaling pathway and FOXD2-AS1/miR-150-5p/PFN2 axis [38,39]. Similarly, LINC00893 has been demonstrated to be expressed at low levels in thyroid cancer cells, which is consistent with our results. Overexpression of LINC00893 can inhibit the proliferation and migration of thyroid cancer cells by blocking the AKT pathway through PTEN upregulation [40].

Study limitations

Certain limitations of the present study should be considered. First, due to the lack of standardized treatment strategies in the TCGA and GEO databases, the predictive accuracy of the model needs to be further confirmed using prospective multicenter randomized controlled trials. Second, the information from the TCGA database is limited and incomplete, which may reduce predictive accuracy of the model. Furthermore, the prognostic model in this study was constructed using bioinformatics analysis. Therefore, validation using cellular experiments and animal and tissue models warrants further investigation.

Conclusions

We identified a metabolism-related lncRNA signature that might act as an independent prognostic variable in patients with BC. Moreover, we established a prognostic nomogram incorporating the gene signature and tumor staging to predict the RFS of BC patients. This novel model may serve as a reliable and reproducible tool for prognostic prediction in individual BC cases.

Click here for additional data file.