Increased tumor glycolytic activity indicates inferior survival in various cancers
A total of 23 types of cancer in 8633 cases were identified in the survival analysis of glycolysis score and OS. A 72-gene expression signature was used to calculate the glycolysis scores. Univariate Cox analysis showed that increased tumor glycolytic activity was significantly associated with lower OS in six types of cancers: liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), head and neck squamous cell carcinoma (HNSC), sarcoma (SARC), brain lower grade glioma (LGG), and pancreatic adenocarcinoma (PAAD) (Table 1). A list of the TCGA cancer type abbreviations is available at https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/tcga-study-abbreviations. To evaluate the prognostic value of glycolysis in pan-cancer patients, we combined HRs for OS. Increased tumor glycolytic activity correlated with inferior OS (HR: 1.70, 95% confidence interval: 1.20–2.40, P = 0.003; Figure 1A). Furthermore, we conducted subgroup analysis to observe the prognostic value of glycolysis in different clinicopathological parameters. We also performed analyses of the effects of glycolysis on OS by subgroups defined by age (<60, ≥60 years), gender (female, male), and pathological stage (early stage, advanced stage). Subgroup analysis revealed that glycolysis is widespread and acts as an unfavorable factor of tumor patients’ prognosis (Table 2).
Table 1. Univariate Cox analysis of glycolysis score.
Study
|
N
|
Hazard ratio
|
95% confidence interval
|
P-value
|
BLCA
|
399
|
1.793
|
0.954-3.367
|
0.070
|
BRCA
|
1052
|
1.894
|
0.974-3.684
|
0.060
|
CESC
|
273
|
2.591
|
0.931-7.214
|
0.068
|
COAD
|
426
|
0.672
|
0.277-1.631
|
0.380
|
ESCA
|
178
|
0.985
|
0.404-2.399
|
0.973
|
GBM
|
146
|
1.034
|
0.400-2.676
|
0.945
|
HNSC
|
512
|
2.819
|
1.553-5.117
|
0.001
|
KIRC
|
518
|
0.972
|
0.460-2.053
|
0.940
|
KIRP
|
278
|
3.941
|
0.840-18.483
|
0.082
|
LGG
|
481
|
3.276
|
1.161-9.243
|
0.025
|
LIHC
|
343
|
8.862
|
3.471-22.628
|
5.00E-06
|
LUAD
|
492
|
3.460
|
1.940-6.173
|
2.60E-05
|
LUSC
|
474
|
0.714
|
0.401-1.271
|
0.252
|
OV
|
295
|
0.724
|
0.316-1.658
|
0.444
|
PAAD
|
172
|
3.034
|
1.113-8.271
|
0.030
|
PCPG
|
172
|
10.703
|
0.112-1020.485
|
0.308
|
PRAD
|
494
|
20.645
|
0.705-604.846
|
0.079
|
READ
|
153
|
0.243
|
0.042-1.403
|
0.114
|
SARC
|
255
|
3.723
|
1.571-8.825
|
0.003
|
STAD
|
375
|
0.609
|
0.341-1.087
|
0.093
|
TGCT
|
130
|
0.787
|
0.015-40.206
|
0.905
|
THCA
|
502
|
4.915
|
0.402-60.072
|
0.213
|
THYM
|
118
|
0.112
|
0.008-1.580
|
0.105
|
UCEC
|
513
|
1.389
|
0.622-3.100
|
0.423
|
Table2. Subgroup analysis for the prognostic value of glycolysis score.
Parameters
|
Types of cancer
|
No. of patients
|
HR (95%CI)
|
P-value
|
Model
|
Age
|
|
|
|
|
|
Age <60
|
23
|
4034
|
1.82 (1.08-3.08)
|
0.025
|
Random effects
|
Age ≥60
|
22
|
4575
|
1.50 (1.01-2.22)
|
0.044
|
Random effects
|
Gender
|
|
|
|
|
|
Female
|
21
|
4477
|
1.62 (1.25-2.11)
|
<0.001
|
Fixed effects
|
Male
|
20
|
4156
|
1.81 (1.11-2.93)
|
0.016
|
Random effects
|
Stage
|
|
|
|
|
|
I/II
|
15
|
3631
|
1.82 (1.07-3.08)
|
0.027
|
Random effects
|
III/IV
|
14
|
2078
|
1.56 (1.17-2.10)
|
0.003
|
Fixed effects
|
To identify the clinical parameters associated with glycolysis status in cancer, we compared the glycolysis scores between different parameters. In many types of cancers, the glycolysis scores displayed a different status for different clinical features (Figure 1B). For example, there were significant differences between early and advanced stage in seven types of cancer. We further explored the status of glycolysis in different molecular subtypes of tumors (Figure 1C–F). In the case of breast cancer, the glycolysis score was lowest in the luminal A subtype, which is the subtype with the best prognosis.
The prognostic value of glycolysis was most significant in LIHC. Hence, three independent cohorts were used to validate the results: GSE14520 (N = 242)(25) and GSE54236 (N = 78)(26), acquired from the Gene Expression Omnibus (GEO) repository, and the LIRI-JP dataset (N = 229)(27), downloaded from the International Cancer Genome Consortium (ICGC) database. The glycolysis score was calculated using GSVA in hepatocellular carcinoma patients with OS of no less than 30 days. Based on the median value of the glycolysis score, LIHC patients were divided into high- and low-glycolysis groups. The patients in the high glycolysis group exhibited lower OS (Figure 2).
Global patterns of glycolysis signatures across cancer types
Based on the tumor expression levels of 72 glycolysis genes, 2047 patients were assigned to a group of high expression of glycolysis-associated genes, and 6585 patients were assigned to a group of low expression of glycolysis-associated genes (Figure 3A). The PCA plot showed that the two clusters were markedly different (Figure 3B). The K–M survival plot showed that patients in the high glycolysis group had inferior OS than patients grouped into the low glycolysis group (Figure 3C). The global expression profiles also indicated differences between the two groups; different glycolysis levels significantly correlated with different clinical parameters (Figure 3D). Notably, different types of tumors accounted for different high-to-low glycolysis group ratios.The ratio of high to low glycolysis group was highest in kidney renal clear cell carcinoma (KIRC) patients and lowest in prostate adenocarcinoma (PRAD) and thyroid cancer (THCA) patients (Figure 3E).
Glycolysis effects on multidimensional molecular factors
Signaling pathways are characterized by frequent somatic alterations and complex interactions with each other. A systematic characterization and exploration of the relationships between oncogenic signaling pathway alterations and glycolysis is becoming meaningful. The glycolysis scores were significantly different between many signaling pathways with and without alterations in different types of cancers (Figure 4A). Glycolysis was upregulated in the PI3K signaling pathway alterations group compared with the no alterations group in 12 types of cancer, and in the cell cycle, tumor protein 53 (TP53), and Hippo signaling pathway alterations group compared with the no alterations group in 10 types of cancer (Figure 4B).
We also explored the glycolysis-associated mRNA expression profiles, as mRNA expression serves as a link between genetic alterations and protein actions. Genes positively correlated with glycolysis in more than 15 cancer types were subjected to functional enrichment analysis. We found that these genes were mainly involved in cell proliferation–associated biological processes and energy-related molecular functions (Figure 5A). The pathway annotations indicated that “Cell cycle” was the most significant relevant term related to glycolysis alterations (Figure 5B). The PPI network analysis also suggested that some genes important in the cell cycle, such as CDK1 and PLK1, were the cores in the network (Figure 5C).
At the protein level, some findings verified the results based on transcriptomics data and led to further discoveries (Figure 6A). Cyclin B1 was significantly positively associated with the glycolysis score in 16 types of tumors. FOXO3A PS318S321, and P27 were markedly negatively related to the glycolysis score in seven types of tumors (Figure 6B).
Glycolytic activity likely positively correlates with immune evasion
Tumor infiltration of lymphocytes is one of the key mechanisms in cancer progression and therapeutic responses. Spearman’s correlation analysis showed that the glycolysis scores significantly negatively correlated with the infiltration of lymphocytes and positively correlated with the infiltration of macrophages and neutrophils (Figure 6). Programmed death-ligand 1 (PD-L1) expression was higher in patients with high glycolysis scores than those with low scores. Increased glycolysis generally correlated with increased T helper (TH) 2 cell infiltration, and TH1 cell infiltration showed the opposite correlation in most types of cancers (Figure 7).
We also explored the relationship between glycolytic activity and genomic instability. Glycolysis was related to measures of DNA damage, including copy number variation (CNV) burden, aneuploidy, homology directed repair (HRD), and intratumor genetic heterogeneity (ITH). Neoantigen load was also associated with higher glycolysis. All these findings suggested that glycolysis was markedly associated with high genomic instability (Figure 8).
Compounds potentially capable of targeting glycolysis
The Cmap online tool was used to discover relationships between genes, compounds, and biological processes and search for valuable targeted molecular compounds for glycolysis. We identified 75 that closely related with glycolysis in not less than 10 cancer types. Sirolimus and LY-294002 as mammalian target of rapamycin (mTOR) inhibitors were found to target glycolysis in all 23 types of cancers (Figure 9A). Cmap mode-of-action (MoA) analysis revealed that 42 mechanisms of action were shared by 44 of the compounds (Figure 9B). Nine compounds shared the MoA of dopamine receptor antagonists.