Identifying glycolysis-related LncRNAs for predicting prognosis in breast cancer patients.

PURPOSE: Functions associated with glycolysis could serve as targets or biomarkers for therapy cancer. Our purpose was to establish a prognostic model that could evaluate the importance of Glycolysis-related lncRNAs in breast cancer.
METHODS: Gene expressions were evaluated for breast cancer through The Cancer Genome Atlas (TCGA) database, and we calculated Pearson correlations to discover potential related lncRNAs. Differentially expressed genes were identified via criteria of FDR < 0.05 and |FC|> 2. Total samples were separated into training and validating sets randomly. Univariate Cox regression identified 14 prognostic lncRNAs in training set. A prognostic model was constructed to evaluate the accuracy in predicting prognosis. The univariate and multivariate Cox analysis were performed to verify whether lncRNA signature could be an independent prognostic factor The signature was validated in validating set. Immune infiltration levels were assessed.
RESULTS: Eighty-nine differentially expressed lncRNAs were identified from 420 Glycolysis-related lncRNAs. 14 lncRNAs were correlated with prognosis in training set and were selected to establish the prognostic model. Low risk group had better prognosis in both training (p= 9.025 e -10) and validating (p= 4.272 e -3) sets. The univariate and multivariate Cox analysis revealed that risk score of glycolysis-related lncRNAs (P< 0.001) was an independent prognostic factor in both training and validating sets. The neutrophils (p= 4.214 e -13, r=-0.223), CD4+ T cells (p= 1.833 e -20, r=-0.283), CD8+ T cells (p= 7.641 e -12, r=-0.211), B cells (p= 2.502 e -10, r=-0.195) and dendritic cells (p= 5.14 e -18, r=-0.265) were negatively correlated with risk score of prognostic model. The Macrophage (p= 0.016, r= 0.0755) was positively correlated with the risk score.
CONCLUSION: Our study indicated that glycolysis-related lncRNAs had a significant role to facilitate the individualized survival prediction in breast cancer patients, which would be a potential therapeutic target.

Introduction

Glycolysis was firstly described as a phenomenon that cancer cells metabolize glucose in a manner by Warburg[1, 2]. Glucose metabolism alteration is an important hallmark of tumor cells, which has been applied to the cancer diagnosis and assessment of cancer response. Interestingly, glycolysis is also a potential therapeutic strategy for cancers[3, 4]. It is increasingly evident that glycolysis is able to facilitate tumor growth and promote chemotherapy resistance[5, 6, 7].

Breast cancer is a highly heterogenous disease that is the most commonly diagnosed malignancies in women. Clinical tumor staging, lymph node status, histological grade and molecular characteristics are considered as prognostic factors by current guideline of AJCC (American Joint Committee on Cancer, https://cancerstaging.org). Investigation of prognostic biomarkers may provide new insight to develop effective therapeutic strategy.

Long noncoding RNAs (lncRNAs) were defined as a type of RNA with more than 200 nucleotides in length and barely encode proteins. LncRNAs are capable of regulating expressions of mRNA at various levels and take part in many critical biological processes. They have a correlation with the diagnosis, treatment and prognosis of breast cancer[8, 9]. However, the mechanism by which lncRNAs regulate gene transcription remains largely unknown. The prognosis of breast cancer patients has a correlation with the status of immune infiltrates [10]. This study aims to identify glycolysis-related lncRNAs and explore the status of immune infiltrates in breast cancer patients.

Materials and methods

Datasets and study cohort

Data of RNA expressions and clinical information of tumor tissues and normal breast tissues (1053 tumor tissues and 111 normal tissues) were obtained from the TCGA database (http://www.cancergenome.nih.gov/). Patients with incomplete survival data or being followed up for less than 30 days were excluded. The genes of Glycolysis were identified from GESA Gene Set (http://www.gsea-msigdb.org/gsea/msigdb/cards/ KEGG_GLYCOLYSIS_GLUCONEOGENESIS).

Bioinformatic analysis and statistical analysis

Genome Reference Consortium (GRCh38) was applied to identified lncRNAs Pearson’s correlation was calculated to screen out glycolysis-related lncRNAs identified as square of correlation coefficient R 0.3 and 0.001. Differentially expressed genes (DEGs) were identified via limma package in R v4.02. (http://www.r-project.org/) with the cut-off criteria of FDR 0.05 and FC 2. The heatmap and volcano map were drawn. Univariate and multivariate Cox regression analyses were performed through R package Survival. Finally, 15 glycolysis-related lncRNAs were identified among total samples.

Entire samples were randomly separated into training and validating sets. In training set, univariate Cox regression analysis were performed again to identify prognostic lncRNAs among training set. We identified 14 prognostic lncRNAs in training set and each gene was compressed into a single score. Sankey diagram and co-expression network were depicted via R v4.02 and Cytoscape software 3.8.0. The risk score was calculated using the following formula to construct the prognostic model.

Kaplan-Meier curves were drawn to compare the survival differences. Clinicopathological variables (e.g., age, stage and TNM staging) and risk scores were compared through Cox univariate and multivariate analyses. The ROC curves were drawn according to clinicopathological variables and risk score. To test the feasibility of our model, the value was validated in validating set.

Data of tumor-infiltrating immune cells were downloaded from the Tumor Immune Estimation Resource (TIMER) database. The infiltrating status of neutrophils dendritic, macrophage T cells and B cells were evaluated to identify their correlation with risk score.

Results

Identification of DEGs and prognostic glycolysis-related lncRNAs

There were 14142 lncRNAs in the breast cancer dataset downloaded in the TCGA database, and among them, 420 were glycolysis-related ones identified according to selection criteria (correlation coefficient R 0.3 and 0.001). Finally, 87 differentially expressed lncRNAs were selected as DEGs (FDR 0.05 and FC 2). The expressions of 87 DEGs were shown in the Additional file 1. The heatmap and volcano map were drawn in Fig. 1A–B. Univariate Cox regression analysis was performed to select the prognostic lncRNAs with 0.05 as the threshold among entire samples. At last, a total of 15 prognostic lncRNAs were identified (Additional file 2). The correlation coefficient and p value were shown in the Additional file 3.

Figure 1.

A: Heatmap of DEGs Blue and red indicate lower expression and higher expression. B: Volcano map of expression of Glycolysis-related lncRNAs. C: The Sankey diagram listed the relationship of 14 lncRNAs and 18 mRNAs. D: The co-expression network listed the relationship of 14 lncRNAs and 18 mRNAs.

Construction of the prognostic model

We randomly divided the entire samples into a training set (480 patient samples) and a validating set (478 patient samples) at a 1:1 ratio. The characteristics of two sets were presented in Table 1. There were no differences in clinicopathological variables between the two sets. We identified 14 prognostic lncRNAs via univariate Cox regression analysis in the training set (Table 2). Sankey diagram (Fig. 1C) and co-expression network (Fig. 1D) listed the relationship of 14 lncRNAs and 18 mRNAs. Among 14 prognostic lncRNAs, 8 lncRNAs have a positive (HRs 1 and 0.05) influence with the prognosis and 6 lncRNAs have a negative (HRs 1 and 0.05) influence with the prognosis (Fig. 1C).

Table 1

The characteristics of training and validation set

Covariates	Type	Total	Test	Train	p value
Age	< 60	520(54.28%)	256(53.56%)	264(55%)	0.7013
Age	≥ 60	438(45.72%)	222(46.44%)	216(45%)
Stage	Stage I–II	725(75.68%)	367(76.78%)	358(74.58%)	0.4737
Stage	Stage III–IV	233(24.32%)	111(23.22%)	122(25.42%)
T	T1-2	811(84.66%)	404(84.52%)	407(84.79%)	0.9781
T	T3-4	147(15.34%)	74(15.48%)	73(15.21%)
M	M0	943(98.43%)	470(98.33%)	473(98.54%)	0.9935
M	M1	15(1.57%)	8(1.67%)	7(1.46%)
N	N0	451(47.08%)	235(49.16%)	216(45%)	0.2202
N	N1-3	507(52.92%)	243(50.84%)	264(55%)

Table 2

Univariate Cox analysis of 14 prognostic glycolysis-related lncRNAs in training set

Id	HR	HR.95L	HR.95H	p value
LINC01614	1.3964652	1.1282007	1.7285178	0.0021532
U62317.1	0.7247198	0.5649086	0.9297412	0.0113057
AL031316.1	0.3779058	0.1869022	0.7641044	0.0067497
AC092142.1	1.7194702	1.1797157	2.5061782	0.0048056
AL136084.3	0.5839348	0.3437030	0.9920769	0.0466593
USP30-AS1	0.4388315	0.2645420	0.7279490	0.0014247
AC092718.4	1.3708613	1.0326176	1.8199001	0.0291104
LINC01235	1.3284102	1.0960503	1.6100298	0.0037931
AC010503.4	1.4792999	1.0081596	2.1706169	0.0453383
AC002546.1	0.2185639	0.0555115	0.8605451	0.0296486
LINC02446	0.5825388	0.3873295	0.8761311	0.0094582
MIAT	0.4990005	0.2745908	0.9068095	0.0225518
LINC01929	1.7198242	1.2426473	2.3802372	0.0010748
LINC01857	0.5763394	0.3767619	0.8816366	0.0110597

Each prognostic lncRNA obtained a risk score, and thus a prognostic model was constructed according to the formula. Finally, breast cancer patients were separated into high risk group and low risk group according to the median risk score for survival analysis. Kaplan-Meier curves indicated that low-risk patients were associated with better prognosis ( 9.025 e -10, Fig. 2A). Survival status and prognostic index distribution were shown in Fig. 2B and C. We applied univariate and multivariate Cox regression analyses to verify whether the risk score of glycolysis-related lncRNAs was an independent prognostic factor for breast cancer. We found that risk score of glycolysis-related lncRNAs ( 0.001) and age ( 0.001) were independent prognostic factors (Fig. 2D and E). The AUC for age, stage, tumor size, lymph node statue, metastasis statue and the risk score was 0.749, 0.795, 0.804, 0.689, 0.554 and 0.819, respectively. Multi-parameter ROC curves revealed that AUC of risk score of glycolysis-related lncRNAs was larger than that of other clinical traits (Fig. 2F). The model was validated in the validating set. Patients in the low-risk group had a better overall survival than those in the high-risk group ( 4.272 e -3, Fig. 3A). The similar survival status and prognostic index distribution were observed in validating set (Fig. 3B and C). The univariate and multivariate Cox analysis revealed that risk score of glycolysis-related lncRNAs ( 0.001) and age ( 0.001) were independent prognostic factors (Fig. 3D and E). The AUC for age, stage, tumor size, lymph node statue, metastasis statue and the risk score was 0.841, 0.56, 0.54, 0.534, 0.58 and 0.642 in the validating set (Fig. 3F). It indicates that a prognostic risk scoring system involving glycolysis-related lncRNAs was more accurate in predicting prognosis of breast cancer patients.

Figure 2.

A: Kaplan-Meier Curve of low-risk group and high-risk groups according to the risk score in training set. B: Distribution of prognostic index in training set. C: Survival status of patients in training set D: Forest plot of Cox univariate analysis in training set E: Forest plot of Cox multivariate analysis in training set. F: Multi-parameter ROC curves for risk score, age, stage, lymph node statue metastasis statue and tumor size in training set.

Figure 3.

A: Kaplan-Meier Curve of low-risk group and high-risk groups according to the risk score in validating set. B: Distribution of prognostic index in validating set. C: Survival status of patients in validating set D: Forest plot of Cox univariate analysis in validating set E: Forest plot of Cox multivariate analysis in validating set. F: Multi-parameter ROC curves for risk score, age, stage, lymph node statue metastasis statue and tumor size in validating set.

The infiltrating status of immune cells

The neutrophils ( 4.214 e -13, 0.223), CD4+ T cells ( 1.833 e -20, 0.283), CD8+ T cells ( 7.641 e -12, 0.211), B cells ( 2.502 e -10, 0.195) and dendritic cells ( 5.14 e -18, 0.265) had a negative correlation with risk score of prognostic model. The Macrophage ( 0.016, 0.0755) had a positive correlation with the risk score (Fig. 4).

Figure 4.

The correlation between neutrophils, CD4+ T cells CD8+ T cells B cells macrophage, dendritic cells and risk score.

Discussion

Breast cancer is a heterogeneous disease that is commonly detected in clinical practice. Although multiple therapeutic strategies for breast cancer are optional, there were half a million deaths of breast cancer in 2012 [11]. In the future, cancer treatment based on tumor biology and early therapy response is the trend to enhance clinical outcomes of cancer patients [12]. Numerous studies have found that lncRNAs are critical regulators of tumorigenesis, growth, tumorigenesis and drug resistance of breast cancer[13, 14, 15, 16]. In our analysis, we investigated expression levels of 420 lncRNAs in breast cancer, narrowing down to 87 DEGs. We divided the entire samples into a training set (480 patient samples) and a validating set randomly. In training set, we identified 6 oncogenic lncRNAs and 8 tumor-suppressor lncRNAs. We constructed a prognostic model according to prognostic genes for breast cancer. We found that glycolysis-related lncRNAs were independent prognostic factors for breast cancer via Cox regression analyses. This result was validated in validating set Lactic acid is a key production induced by glycolysis. The accumulation of lactic acid may promote biological processes via tumor microenvironment [17]. On the other hand, the differentiation of CD4+ and CD8+ T cells relies on glycolysis[18]. The status of immune infiltrates was evaluated to identify their correlation with risk score of glycolysis-related LncRNAs. We observed that risk score was negatively correlated with neutrophils, CD4+ T cells, CD8+ T cells, B cells and dendritic cells. Previous studies found that CD8+ T cells was predictive factor for pathological complete remission to primary systemic therapy and was associated with better clinical outcomes[19, 20, 21]. Although CD4+ T cell may serve as an antitumor factor [22]. Recent study found that high CD4+ T cells had a non-favorable effect on breast cancer patient[21]. The main mechanism of anti-tumor immune response of CD4+ T cells is still being explored.

So far, among 14 prognostic lncRNAs, the function AC002546.1, AC010503.4, AC092142.1, AC092718.4, AL031316.1, AL136084.3 and U62317.1 have not been previously explored in breast cancer or other cancers. LINC01929 has been found positively correlated with tumor progression in liver cancer, oral squamous cell carcinoma and lung cancer[23, 24, 25]. Che et al. indicated that LINC01929 promotes tumor progression through miR-137-3p/FOXC1 axis in oral squamous cell carcinoma [23]. So far, LINC01929 has not been reported in breast cancer yet. We found that LINC01929 acted as a carcinogenic lncRNA in breast cancer, and its specific oncogenic mechanism remains unclear and is worth exploring. In the present study, LINC02446 was consistently identified as a prognostic gene and might serve as a promising therapeutic biomarker. LINC01614 was identified as a non-favorable prognostic biomarker in breast cancer[26, 27], which is consistent with our analysis. Previous studies have indicated that the oncogenic LINC01614 stimulated the development of osteosarcoma, gastric cancer cells, non-small cell lung cancer and lung adenocarcinoma, which provide a new insight into the development of cancer treatment[28, 29, 31]. MIAT was firstly discovered as a susceptible locus for myocardial infarction[32]. In recent years, a growing number of evidence has indicated the vital effect of MIAT in regulating the development of human cancers[33]. Previous study indicated that silence of MIAT inhibits proliferation of MCF-7 cells, and overexpression of it yields the opposite trend[34]. In addition, MIAT has a correlation with 5-fluorouracil (5-FU) resistance of MCF-7 cell. Yao et al. found that the expression level of GRP78, MIAT, OCT4 and AKT were dysregulated in 5-FU resistant MCF cells. These genes contributed to drug resistance [35]. However, we found that upregulation of MIAT was associated with a better overall survival. It is reported that the MIAT is lower-expressed in hormone receptor-negative tumor tissues than that in hormone receptor-positive ones. Hormone receptor-positive subtypes account for 60% of breast cancers and their prognosis is better than that of hormone receptor-negative subtypes[11, 36]. In addition, expression level of MIAT increases in early-stage breast cancer, rather than stage III and IV[34, 37]. These may attribute to the inconsistent findings.

Our findings suggested that glycolysis-related lncRNAs may serve as therapeutic targets against breast cancer progression. Several limitations existed in our analysis Our data was downloaded from the TCGA database, and most of breast cancer patients were primarily Americans. Breast cancer patients from other countries should be explored. On the other hand, all samples only record the breast cancer related events and we are limited in discussing the event of other organ sites. Experiments of in vitro and in vivo were lacked to validate our data mining results, which are needed in the future.

In summary, we have identified transcripts corresponding to 14 genes that may serve as prognostic markers in training set and constructed a prognostic model of breast cancer. This analysis provides a promising method in predicting individualized survival of breast cancer. In the prognostic model, risk score were correlated with infiltrating status of neutrophils dendritic, macrophage T cells and B cells. Our results might provide potential therapeutic targets for breast cancer.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.