A novel glycolysis-related gene signature for predicting prognosis and immunotherapy efficacy in breast cancer

Abstract

Background: Previous studies have shown that glycolysis-related genes (GRGs) are associated with the development of breast cancer (BC), and the prognostic significance of GRGs in BC has been reported. Considering the heterogeneity of BC patients, which makes prognosis difficult to predict, and the fact that glycolysis is regulated by multiple genes, it is important to establish and evaluate new glycolysis-related prediction models in BC.

Methods: In total, 170 GRGs were selected from the GeneCards database. We analyzed data from the Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA) database as a training set and data from the Gene Expression Omnibus (GEO) database as a validation cohort. Based on the overall survival data and the expression levels of GRGs, Cox regression analyses were applied to develop a glycolysis-related prognostic gene (GRPGs)-based prediction model. Kaplan (KM) survival and ROC analyses were performed to assess the performance of this model. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used to identify the potential biological functions of GRPGs. cBioPortal database was used to explore the tumor mutation burden (TMB). The tumor immune dysfunction and exclusion indicator (TIDE) was used to estimate the patient response to immune checkpoint blockade (ICB). The levels of tumor-infiltrating immune cells (TICs) and stromal cells were quantitatively analyzed based on gene expression profiles.

Results: We constructed a prediction model of 10 GRPGs (ADPGK, HNRNPA1, PGAM1, PIM2, YWHAZ, PTK2, VDAC1, CS, PGK1, and GAPDHS) to predict the survival outcomes of patients with BC. Patients were divided into low- and high-risk groups based on the gene signature. The AUC values of the ROC curves were 0.700 (1-year OS), 0.714 (3-year OS), 0.681 (5-year OS). TMB and TIDE analyses showed that patients in the high-risk group might respond better to ICB. Additionally, by combining the GRPGs signature and clinical characteristics of patients, a novel nomogram was constructed. The AUC values for this combined prediction model were 0.827 (1-year OS), 0.792 (3-year OS), and 0.783 (5-year OS), indicating an outstanding predictive performance.

Conclusion: A new GRPGs based prediction model was built to predict the OS and immunotherapeutic response of patients with BC.

Keywords: bioinformatics; breast cancer; glycolysis; prognostic signature; the cancer genome atlas.

Figure 1

Overall workflow of this study.

Figure 2

GRPGs selection using univariate and multivariate Cox regression analyses. (A) The Venn diagram displayed how 169 GRGs were selected. (B) Univariate Cox regression analysis selected 10 GRPGs correlated with OS. (C) Multivariate Cox regression analysis result was shown by the forest plot. (D) The risk factor diagram showed the risk score distribution, the survival status of BC patients, and the gene expression levels of 10 GRPGs.

Figure 3

Differential expression and GO/KEGG enrichment analyses of GRPGs. Differential expression analyses of 10 GRPGs between BC tumor and adjacent normal tissues were performed in TCGA-BRCA (A), GSE42568 (B), GSE29044 (C). Biological process (BP) and KEGG enrichment analyses of 10 GRPGs were shown in histogram (D) and network diagrams (E, F). **: P value<0.01, ***: P value<0.001.

Figure 4

Clinicopathological and survival information of the low-risk and high-risk group for TCGA-BRCA. (A) T stage, (B) N stage, (C) M stage, (D) pathologic stage, (E) age, (F) OS, (G) DSS, (H) PFI, (I)ER status, (J) PR status, (K) HER2 status, (L) TNBC or non-TNBC.

Figure 5

Kaplan-Meier survival analyses in BC patients based on risk stratification and the expression level of each GRPG. The OS difference between low-risk and high-risk group was shown in (A). Kaplan-Meier survival analyses were based on the expression levels of CS (B), PIM2 (C), PGK1 (D), GAPDHS (E), HNRNPA1 (F), ADPGK (G), YWHAZ (H), PTK2 (I), PGAM1 (J), VDAC1 (K) in TCGA-BRCA. (L) The OS difference between low-risk and high-risk group in GSE20685 was displayed by the KM curves. (M) AUC values were calculated in ROC analysis for risk scores predicting the OS from TCGA-BRCA. (N) Differential expression analyses of 10 GRPGs between low-risk group and high-risk group were performed in TCGA-BRCA. (O) The expression patterns of 10 GRPGs were shown in the heatmap.

Figure 6

Differentially expressed genes (DEGs) of low-risk group versus high-risk group. (A) 1148 DEGs were shown in the Volcano plot (|logFC| > 0.5 and adjusted P<0.05). (B) 20 enriched biological functions obtained by GSVA analysis were shown in the heatmap. (C) Mountain plot showed the four main biological features of DEGs achieved by GSEA enrichment analysis. DEGs were significantly enriched in oxidative stress induced senescence (D), cellular senescence (E), folate metabolism (F) and primary immunodeficiency (G).

Figure 7

Estimation of TMB, TIDE, stromal score, immune score, ESTIMATE score and tumor purity score in low-risk and high-risk group. (A) The histogram showed the differential TMB scores between the low-risk and high-risk group. (B) A positive correlation between risk scores and TMB scores was found by Spearman correlation test. (C) The histogram showed the differential TIDE scores between the low-risk and high-risk group. (D) A negative correlation between risk scores and TMB scores was found by Spearman correlation test. Violin plots showed the differential stromal score (E), immune score (F), ESTIMATE score (G) and tumor purity score (H) between the high-risk and low-risk groups. Risk scores were negatively correlated with stromal score (I), immune score (J) and ESTIMATE score (K) but were positively correlated with tumor purity score (L). ***: P value<0.001.

Figure 8

Estimation of TICs in low-risk and high-risk group. (A) The differences of 28 TICs between low-risk and high-risk group were evaluated by ssGSEA algorithm. The correlation analyses of TICs were conducted in low-risk group (B) and high-risk group (C). The dot plots showed the correlation between the abundance of TICs and the expression levels of GRPGs in low-risk group (D) and high-risk group (E). *: P value<0.05, ***: P value<0.001, ns: P values≥0.05.

Figure 9

Genetic alterations of GRPGs in low-risk and high-risk group. The combined graph displayed the variant classification, variant types, single nucleotide variations (SNV) class, the number of variants and top 10 mutated genes in low-risk group (A) and high-risk group (B). The waterfall plots showed the genetic alterations of GRPGs sorted by mutation rate in low-risk group (C) and high-risk group (D). The histograms showed the top 20 gene amplifications in low-risk group (E) and high-risk group (F). The histograms showed the top 20 gene deletions in low-risk group (G) and high-risk group (H).

Figure 10

Specific genetic alterations and corresponding number of samples depending on different molecular types of BC. (A) ER-, (B) ER+, (C) HER2-, (D) HER2+, (E) non-ER+HER2+, (F) ER+HER2+, (G) non-TNBC, (H) TNBC.

Figure 11

Differential expression of 10 GRPGs and risk scores in different BC subgroups. Differential analyses of the expression of 10 GRPGs was conducted in ER+/- BC (A), HER2+/- BC (B), ER+HER2+/non-ER+HER2+ BC (C) and TNBC/non-TNBC (D). The heatmaps showed the expression patterns of 10 GRPGs in ER+/- BC (E), HER2+/- BC (F), ER+HER2+/non-ER+HER2+ BC (G) and TNBC/non-TNBC (H). The histograms showed risk scores in ER+/- BC (I), HER2+/- BC (J), ER+HER2+/non-ER+HER2+ BC (K) and TNBC/non-TNBC (L). *: P value<0.05, **: P value<0.01, ***: P value<0.001.

Figure 12

Survival subgroup analyses based on the clinicopathological variables. (A) T1-T2 stage, (B) T3-4 stage, (C) N0 stage, (D) N (+) stage (N1-N3), (E) M0 stage, (F) M1 stage, (G) pathologic stage I-II, (H) pathologic stage III-IV, (I) age ≤ 60, (J) age>60, (K) ER-, (L) ER+, (M) PR-, (N) PR+, (O) HER2-, (P) HER2+, (Q) non-ER+HER2+, (R) ER+HER2+, (S) non-TNBC, (T) TNBC.

Figure 13

Construction a nomogram in TCGA-BRCA. Univariate (A) and multivariate (B) Cox analyses were performed to analyze several clinicopathological data. (C) The nomogram consisted of TNM stage, age, ER status, PR status, HER2 status and risk scores to predict the probability of 1-year, 3-year and 5-year OS.

Figure 14

Evaluation of this nomogram. Calibration curves of 1-year (A), 3-year (B) and 5-year (C) OS predicted by the nomogram showed the relationship between predicted survival probability and observed fraction survival probability. DCA curves of 1-year (D), 3-year (E) and 5-year (F) OS prediction showed the clinical predictive effects of this combined prediction model. The OS difference between combined low-risk and high-risk group was shown in (G). (H) AUC values were calculated in ROC analysis for combined risk scores predicting the OS from TCGA-BRCA.