Pan-Cancer analysis of clinical signi cance and associated molecular features of glycolysis

ABSTRACT Tumor glycolysis is a major promoter of carcinogenesis and cancer progression. Given its complex mechanisms and interactions, comprehensive analysis is needed to reveal its clinical significance and molecular features. On the basis of a well-established glycolysis gene expression signature, we quantified 8633 patients with different cancer types from the Cancer Genome Atlas (TCGA) and evaluated their prognostic associations. High tumor glycolytic activity correlated with inferior overall survival in the pan-cancer patients (hazard ratio: 1.70, 95% confidence interval: 1.20–2.40, P = 0.003). The prognostic value of glycolysis correlated with the molecular subtypes and was stable regardless of clinical parameters. The prognostic significance of glycolysis was validated using three independent datasets. In addition, genome, transcriptome, and proteome profiles were utilized to characterize the distinctive molecular features associated with glycolysis. Mechanistically, glycolysis fulfilled the fundamental needs of tumor proliferation in multiple ways. Exploration of the relationships between glycolysis and tumor-infiltrating immune cells showed that glycolysis enabled the immune evasion of tumor cells. Mammalian target of rapamycin (mTOR) inhibitors and dopamine receptor antagonists can effectively reverse the glycolytic status of cancers. Overall, our study provides an in-depth molecular understanding of tumor glycolysis and may have practical implications for clinical cancer therapy.


Background
Teleologically, carcinogenesis and tumor cell development are linked with and dependent on metabolic reprogramming to meet their needs in energy and macronutrient requirements (1,2). In the alterations of metabolic reprogramming, the Warburg effect refers to the tendency of cancer cells to utilize glucose via glycolysis irrespective of oxygen availability (3,4). Tumor glycolysis is a biological phenotype of most tumors and has served as the basis for cancer detection with positron emission tomography. Glycolysis has been associated with advanced tumor progression, treatment resistance, and poor clinical outcomes (5,6). More importantly, it has also signi cantly been linked to many cancer molecular characteristics, including proliferation, angiogenesis, and immune evasion (7,8). Hence, targeting tumor glycolysis is viewed as of major importance in cancer therapy (9). Nevertheless, a systematic characterization of the clinical and molecular characteristics of glycolysis is still needed.
In recent years, The Cancer Genome Atlas (TCGA) multi-omics data has driven the understanding of the molecular landscape of primary beyond individual molecular platforms by integrating genomic, transcriptomic, and proteomic characteristics and clinicopathological parameters (10). The tumor glycolysis status has also been linked to multiple layers of molecular alterations (11)(12)(13)(14). Many excellent studies on the mechanisms of glycolysis have presented a clearer picture of tumor characteristics. For example, upregulated tumor glycolysis was found to modulate T cell-mediated antitumor activity, thereby inhibiting melanoma patients' response to adoptive T cell therapy (15). Glycolysis could also endow maintaining strong tumorigenic activity on transcription factors such as YAP/TAZ (16).
Furthermore, metabolism of glycolytic modulates the translation of hypoxia-inducible factor 1-alpha (HIF1A) to control T cell responses to hypoxia (17). Taken together, these studies have established that glycolysis can lead to multiple layers of molecular changes and thus plays a pivotal role in cancer development. Yet, the underlying mechanistic details pertinent to the cancer glycolysis phenotype remain unclear.
With the availability of multi-omics data, we aimed to quantify glycolysis and reveal its prognostic value from pan-cancer perspectives. By multi-omics data integration of their genome, transcriptome, and proteome, we also aimed to identify molecular characteristics related to glycolysis. Finally, we aimed to identify speci c molecular compounds that may eventually result in novel therapy strategies that can reverse the glycolytic status of cancers. These approaches provide an unprecedented opportunity to explore the biological characteristics of tumors in great depth.

Materials And Methods
Classi cation of tumor glycolysis scores across various cancer types Molecular pro les including oncogenic signaling pathway, messenger RNA (mRNA) expression, and protein expression pro les were acquired from the TCGA pan-cancer project (https://gdc.cancer.gov/about-data/publications/pancanatlas). Genome-wide RNA sequencing (RNAseq) data pro les were obtained from the TCGA dataset using an Illumina HiSeq 2000 RNA Sequencing System. Cancer-relevant proteins and phosphoproteins were detected by reverse-phase protein arrays.
The gene signature for glycolysis (REACTOME_GLYCOLYSIS) was de ned based on 72 genes involved in glycolysis from the Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb). To obtain more reliable results for survival analysis, only primary tumor samples with overall survival (OS) of no less than 30 days were included in further analysis. Finally, we analyzed 23 TCGA non-hematologic cancer types with a sample size ≥100. The glycolysis score for each tumor sample across cancer types was calculated using Gene Set Variation Analysis (GSVA) (18).
Estimation of the prognostic value of glycolysis in pan-cancer Univariate Cox analysis was performed to explore the relationship between glycolysis scores and tumor patients' OS. Then, meta-analysis was performed to integrate the univariate Cox survival analysis results into each type of cancer. Heterogeneity between individual studies was evaluated by the Q test: I2 > 50% and/or P ≤ 0.05 indicated signi cant heterogeneity. Random effects or xed effects model was identi ed based on the heterogeneity analysis results. Meta-analysis was performed using Stata 14.0 (Stata Corporation, College Station, TX, USA). To further explore variations in the effect of glycolysis on OS, subgroup analyses strati ed by age, gender, and tumor stage were performed.
We further performed consensus cluster analysis to identify a global pattern of glycolysis across the 23 types of cancers using the ConsensusClusterPlus package in R software (19). Patients were grouped into two clusters according to glycolysis-related genes using two-class K-means clustering with the Euclidean distance. Principal component analysis (PCA) was used to observe the patients in the two clusters. The OS difference between the clusters was calculated by the Kaplan-Meier (K-M) method and the log-rank test.
We also computed the correlation between glycolysis scores and mRNA and protein expression pro les.
Spearman's rank correlation was used to assess the correlation between glycolysis scores and molecular factors. Glycolysis-associated genes and proteins in each cancer type were determined as follows: Spearman's correlation |coe cient| > 0.3 and P < 0.05. We found that genes signi cantly positively correlated with glycolysis scores in at least 15 cancer types were glycolysis-associated genes. Glycolysisassociated genes were further subjected to a process of gene ontology (GO) classi cation analysis using the R package clusterPro ler (21). To ascertain functional interactions between these glycolysis-related genes and found hub genes, protein-protein interaction (PPI) networks were constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database version 11.0 (required interaction score: >0.9)(22) .

Exploration of the relation between increased glycolysis and immune evasion
Growing evidence suggests that increased glycolysis actively interferes with immune cell functions. Thus, we aimed to unveil the association between glycolysis and immune status. The tumor immune microenvironment is an executor of immunotherapy. Using CIBERSORT, we obtained 22 types of immune cells in 6 classes from the TCGA pan-cancer project: (1) lymphocytes (B cells naïve, B cells memory, T cells CD4 naïve, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory Tregs, T cells gamma delta, T cells CD8, NK cells resting, NK cells activated, plasma cells); (2) macrophages (monocytes, macrophages M0, macrophages M1, macrophages M2); (3)
Antigen-speci c T cell receptor and B cell receptor repertoires are critical for the the process of malignant cells recognition and may re ect a robust antitumor response comprising a large number of antigenspeci c adaptive immune cells.
Owing to the complexity of immune evasion, we also explored the associations between glycolysis and DNA damage, including copy number variation (CNV) burden (in terms both of number of segments and of fraction of genome alterations), aneuploidy, loss of heterozygosity (LOH), homologous recombination de ciency (HRD), predicted neoantigen number (including single-nucleotide variant [SNV] and indel neoantigens), and intratumor heterogeneity (ITH). The detailed process of quantitative analysis of these indicators was described in a previous study (29628290).

Compounds targeting cancer glycolysis
To explore which molecules could be effective in cancer-inhibitory glycolytic activity, we performed Connectivity Map (Cmap) analysis (24). Cmap is a public online database similar to Gene Set Enrichment Analysis (GSEA), which predicts target drugs based on a query signature. Gene symbols were mapped to the HG-U133A probe set GPL96 platform ID. We identi ed 500 genes most positively associated and 500 genes most negatively associated with glycolysis scores to query drugs matching the "reference" signature. Compounds with an enrichment score <0 and a P-value <0.05 were identi ed as glycolysis antagonist drugs.

Increased tumor glycolytic activity indicates inferior survival in various cancers
A total of 23 types of cancer in 8633 cases were identi ed 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 signi cantly 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% con dence 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 de ned 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). 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 signi cant 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 signi cant 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 pro les also indicated differences between the two groups; different glycolysis levels signi cantly 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 signi cantly 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 pro les, 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 signi cant 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 ndings veri ed the results based on transcriptomics data and led to further discoveries ( Figure 6A). Cyclin B1 was signi cantly 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 in ltration of lymphocytes is one of the key mechanisms in cancer progression and therapeutic responses. Spearman's correlation analysis showed that the glycolysis scores signi cantly negatively correlated with the in ltration of lymphocytes and positively correlated with the in ltration 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 in ltration, and TH1 cell in ltration 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 ndings 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 identi ed 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.

Discussion
In this study, we systematically analyzed the prognostic value and effects of glycolysis on molecular signatures in 8633 primary human tumors across 23 cancer types. We gained insight into the relationship between glycolysis and low pan-cancer tumor survival. Our study provides a comprehensive view of glycolysis-associated molecular signatures, including genome, transcriptome, and proteome data. We provided a description of the relation between glycolysis and immune evasion from immune cell in ltrations and immune-related indicators. And our exploration led to potentially actionable compounds as candidates for alternative metabolism therapy for solid tumors. Our integrative analysis further suggests that glycolysis can impact tumors in many ways.
Metabolic reprogramming ful lls tumors' energy/nutrient requirements and is considered a hallmark of cancer (2). Previous studies found that a higher glycolysis level signi cantly correlates with a higher risk of adverse events or death in many types of cancers, including head and neck (28), lung (29), and esophageal cancer (30). With advancements in next-generation technologies and public cancer databases, it is now possible to comprehensively analyze the prognostic value of glycolysis in pancancer. We found that increased tumor glycolytic activity is associated with poor clinical outcomes in pan-cancer samples. This nding is consistent with recent work demonstrating that high glycolytic activity positively correlates with higher death risk. The prognostic signi cance was stable in patients with different clinicopathological. These results indicate fundamental and universal traits of glycolysis in cancer prognosis.
By integrating multi-omics data, we observed that oncogenic glycolysis to be related with some molecular characteristics. Multi-omics analysis can link and integrate multilayered information and provide more reliable and accurate results than single platform analysis. The rst one is perturbations in oncogenic signaling pathways. Several important signaling pathways have been identi ed as frequently genetically altered in cancer (31). Patients with PI3K, cell cycle, and TP53 pathway alterations were found to have high glycolysis scores in multiple tumor types. At the transcriptome level, we identi ed glycolysisassociated genes, and subsequent gene functional enrichment analysis also indicated essential cell cycle and proliferation processes. Proteomics analysis revealed that cell cycle regulatory protein cyclin B1 is an important participant in glycolysis process. Different multi-omics data suggested that glycolysis is inextricably related with tumor proliferation. The most basic characteristic of cancer cells is continuous proliferation (2). It can be said that glycolysis rst meets the most basic requirements of tumors.
Carcinogenesis and cancer progression modulated by tumor glycolysis have aroused considerable interest. However, the immune regulatory role of glycolysis is still unclear. The high concentration of lactate produced by the glycolytic process in TME suppresses the anticancer immune cells by disturbing their intracellular pH (32,33). Furthermore, cancer glycolytic activity competes with immune cells for glucose uptake (34,35), whereby cancer cells effectively achieve immune evasion. We found that increased tumor glycolytic activity is related to poor tumor lymphocyte in ltration in the majority of cancer types. Previous studies also found that it is inversely associated with tumor in ltration of T cells in melanoma, non-small-cell lung carcinoma (NSCLC), and HNSC samples (15,36). As for T helper cells, increased glycolysis scores strongly correlated with high TH2 and low TH1 cell in ltration, which is an interesting phenomenon. Preferential accumulation of immunosuppressive and TH2 cells, rather than antitumor TH1 cells, have been found to be indispensable in tumor immune evasion (37,38). Similarly, the relationships between glycolytic activity and the PD-L immune checkpoint and TCR richness also suggest that increased glycolysis promotes tumor proliferation by immune evasion.
We queried Cmap using the gene expression signatures from glycolysis-associated gene analysis.
Perhaps surprisingly, the Cmap analysis, which is based on a limited number of treated cell lines, very precisely selected drugs that have been shown to affect cancer metabolism with speci city. MTOR inhibitors were signi cantly enriched in all cancer types and have been reported to inhibit glycolysisrelated tumorigenicity (39). Interestingly, mTOR actively participates in T cell function and activity (40). Dopamine receptor antagonists, the most frequent compounds in the Cmap analysis, also showed a capacity to inhibit tumor progression (41). These translational analyses may provide alternative approaches for metabolic cancer treatment.

Conclusion
In summary, glycolysis could be a reliable predictor of clinical outcomes of tumors. This study also provides a comprehensive catalogue of molecular alterations associated with glycolysis and contributes to our understanding of glycolysis, perhaps leading to effective survival prediction, treatment decisionmaking, and target therapy identi cation.

Declarations
Ethics approval and consent to participate This study was approved by the First A liated Hospital of Guangxi Medical University as a retrospective study, and the requirement for informed consent was waived. All procedures performed in studies involving human participants were in accordance with the ethical standards of the1964 Helsinki declaration and its later amendments. and interpretation: Yi-chen Liu, Peng Lin and Lin-yong Wu and Yu-jia Zhao; (V) Manuscript writing: All authors; (VI) Final approval of manuscript: All authors.       Associations between glycolysis and immune cell in ltrations.

Figure 8
Associations between glycolysis and immune-related index and genomic instability indicators.