Alterations in glycolytic/cholesterogenic gene expression in hepatocellular carcinoma

Four metabolic subgroups of HCC were identified based on the dual analysis of glycolytic and cholesterogenic gene expression

In total, 610 HCC tumor samples were included in this study (The Cancer Genome Atlas [TCGA], n = 373, and International Cancer Genome Consortium [ICGC], n = 237). Samples with < 30% were excluded. Reactome gene set enrichment analysis was applied to obtain the “glycolysis” (n = 29) and “cholesterol biosynthesis” (n = 72) gene sets. Unsupervised consensus clustering analysis was implemented to classify the two subgroups of significantly expressed glycolytic (n = 9) and cholesterogenic (n = 11) genes for further metabolic subtyping (Figure 1A).

Figure 1. The metabolic gene landscape of HCC based on glycolytic and cholesterogenic clusters. (A) Heat map of consensus clustering (k=5) for glycolytic and cholesterogenic genes in resected and metastatic LIHC samples (n=610). (B) Scatter plot of the median levels of co-expressed glycolytic (x-axis) and cholesterogenic (y-axis) genes in each LIHC sample. Metabolic subgroups were assigned based on the relative levels of glycolytic and cholesterogenic genes. (C) Heat map of differential gene expression patterns in glycolytic and cholesterogenic gene clusters across subgroups. (D) Kaplan-Meier survival analyses of patients with all subtypes of LIHC; the log-rank test p value is shown. (E) Overall survival analyses in the metastatic subgroup of LIHC patients; the log-rank test p value is shown. (F) Overall survival analyses in the non-metastatic LIHC cohort; the log-rank test p value is shown.

We calculated the median levels of glycolytic and cholesterol-producing genes in each sample. Based on the co-expression of the two gene sets, we separated the gene profiles into four metabolic subtypes of HCC: the glycolytic subgroup (high expression of glycolytic genes and low expression of cholesterol synthesis genes), cholesterogenic (high expression of cholesterol synthesis genes and low expression of glycolytic genes), mixed (high expression of both cholesterol synthesis genes and glycolytic genes) and quiescent (low expression of both cholesterol synthesis genes and glycolytic genes) subgroups (Figure 1B). The levels of genes involved in glycolysis and cholesterol synthesis are illustrated in Figure 1C, along with the proportion of patients in each metabolic subgroup. The quiescent phenotype group contained the largest number of patients (332/610, 54.4%), followed by the cholesterogenic (164/610, 26.9%), glycolytic (65/610, 11%) and mixed subtypes (49/610, 8%).

To determine the relationship between the HCC metabolic subgroup and patient prognosis, we performed a statistical cluster analysis on HCC metabolic oncogene characteristics in metastatic and non-metastatic patients. There was no significant difference in the distribution of the metabolic subgroups between the metastatic and non-metastatic groups. However, there were significant differences in overall survival based on cholesterol-generating and glycolytic gene expression. The clinical outcomes were significantly worse in the cholesterol-generating subgroup than in the glycolytic group (p = 0.017) (Figure 1D). In both the non-metastatic and metastatic groups of HCC patients (Figure 1E and 1F, p = 0.049 and p = 0.027, respectively), there was a significant difference in survival between the cholesterogenic and glycolytic subgroups. Surprisingly, survival benefits were observed in patients with increased expression of glycolytic genes. These findings indicated that there are metabolic phenotypes associated with glycolysis/cholesterogenesis and prognosis in HCC.

Relationship between the tumor genome metabolic subtype and the HCC subtype

Metabolic reprogramming is now accepted as a hallmark of cancer [20, 21], and we found significant abnormalities in metabolic gene expression in HCC patients. We next investigated the frequency of single nucleotide variations (SNVs), insertion-deletion mutations (INDELs) and copy number variations (CNVs) [22] in genes associated with the different metabolic subtypes of liver hepatocellular carcinoma (LIHC) cohorts (Figure 2A). The mutational frequency of each gene did not differ significantly among the subtypes (Fisher’s exact test and Benjamini-Hochberg (BH) test-corrected, p > 0.05 after correction). However, the median levels of cholesterol synthesis genes were significantly greater in MYC-amplified and TP53-deleted samples than in samples without these alterations (Figure 2B). Cholesterogenic gene expression correlated positively with MYC mutations (p = 0.039, R = 0.1) (Figure 2C) and negatively with TP53 mutations (p = 0.037, R = -0.14) (Figure 2D). These findings are consistent with the high mutational frequency of TP53 and MYC, and support the notion that mutations in TP53 and MYC promote tumor progression by inducing abnormal glucose utilization. Thus, tumors may be vulnerable to changes in glycolysis.

Figure 2. Gene mutational landscape across metabolic subgroups of HCC. (A) Oncoprint analysis indicating the distribution of SNVs, INDELs and CNVs of frequently mutated genes in LIHC across the metabolic subtypes. (B) Box plot of the median expression of cholesterogenic genes in samples with CNVs in TP53 and/or MYC. (C) Scatter plot of the correlation between the median cholesterogenic gene expression and MYC expression. (D) Scatter plot of the relationship between the median cholesterogenic gene expression and TP53 expression.

The relationship between survival-related features and gene expression in HCC has been illustrated in previous studies [23, 24]. Hoshida’s classification, Budhu’s metastasis-averse/inclined microenvironment (MAM/MIM) and Chew’s classification have been used to classify liver cancer according to the clinical prognosis (see Methods for details) (Supplementary Table 1). To explore the association between the glycolytic/cholesterogenic metabolic phenotypes and the tumor subtypes from the above-mentioned prognostic classification methods, we identified the tumor subtype and metabolic phenotype of each sample (Figure 3A). Most of the patients in the quiescent group had favorable prognoses according to the Hoshida classification (75.4%), while fewer patients in the mixed group (29.9%, adjusted p = 0.0025) and the cholesterogenic group (34.8%, adjusted p = 0.022) had favorable prognoses (Figure 3B). The number of MAM/MIM samples according to the Budhu classification differed significantly between the quiescent group and the mixed group. The number of samples with a poor prognosis according to the Chew classification differed significantly between the quiescent group and the glycolytic group.

Figure 3. The alignment of LIHC metabolic subgroups with known gene expression subtypes. (A) Overlay of the metabolic gene profiles with LIHC expression subtypes based on the known classifications of Hoshida et al., Budhu et al. and Chew et al. (B) Bar plots of the proportion of LIHC expression subtypes in each metabolic subgroup. (C) Scatter plots depicting the correlations of glycolytic and cholesterogenic gene levels with prognostic gene levels in the Hoshida classification.

We also analyzed the correlation of cholesterogenic/glycolytic gene expression with prognostic gene expression according to the Hoshida classification. The expression of genes associated with a poor prognosis correlated positively with the expression of genes in the cholesterol synthesis pathway. Likewise, the expression of genes associated with a good prognosis correlated negatively with the expression of genes in the cholesterol synthesis pathway. Glycolytic gene expression correlated positively with prognostic gene expression in both the better and poorer prognosis groups, but the correlation was stronger in the poorer prognosis group (Figure 3C). This was consistent with the significant association between cholesterogenic gene expression and a poorer survival prognosis. These data indicated that glycolysis and cholesterol synthesis are potential metabolic targets in patients with different HCC subtypes.

Mitochondrial pyruvate carrier complex expression differed among the metabolic subtypes

The mitochondrial pyruvate carrier (MPC) complex in the inner membrane of the mitochondria transfers free pyruvate from the cytoplasm to the mitochondrial matrix [25]. Previous studies observed that MPC1 and MPC2 were different in metabolic pathways and promoted tumor glycolysis activity. This difference is important in lactic acid production [26, 27].

To explore the association of MPC1 and MPC2 with the different metabolic phenotypes, we assessed the metabolic subpopulations for SNVs, INDELs, CNVs and mRNA expression changes in these genes. The CNVs in MPC1 and MPC2 differed significantly among the metabolic subgroups; however, the CNVs in MPC1 were deletions in almost all of the HCC samples, while the CNVs in MPC2 were mostly amplifications (Figure 4A). MPC1 and MPC2 mRNA levels differed significantly among the metabolic subgroups. MPC2 expression was significantly lower in the glycolytic group than in the other groups, whereas MPC1 expression was significantly greater in the glycolytic group. MPC2 expression was significantly higher in the cholesterogenic group than in the quiescent group (Figure 4B). Thus, the dysregulation of mitochondrial pyruvate transport at the mRNA level may be associated with the metabolic tumor subtypes.

Figure 4. Association of MPC1 and MPC2 expression with LIHC metabolic subgroups and cell signaling pathways. (A) Oncoprint indicating the distribution of MPC1 and MPC2 SNVs and CNVs across the metabolic groups. Only one case was found with an SNV in MPC2. (B) Box plots of significant (p < 0.001) differences in MPC1 and MPC2 levels across the LIHC metabolic subgroups. (C) Scatter plot of the correlations of 25,483 genes with MPC1 (x-axis) and MPC2 (y-axis). In total, 168 genes correlated positively (Spearman correlation BH-adjusted p < 0.01) with both MPC1 and MPC2, while 14 genes correlated negatively with both MPC1 and MPC2 (adjusted p < 0.01). (D) The most significantly enriched (hypergeometric test BH-adjusted p < 0.05) gene sets among the genes positively (left) and negatively (right) associated with MPC1/2 expression.

To find cellular pathways associated with MPC1/2 expression, we performed a comprehensive correlation analysis between MPC1/2 and all the other tested genes (n = 25,483). In total, 168 genes correlated positively with MPC1/2, while 14 genes correlated negatively with MPC1/2 (Spearman correlation BH-corrected p < 0.01) (Figure 4C). The positively correlated genes were associated with extracellular matrix function (hypergeometric test, BH regulation p < 0.05). The pathways enriched in the negatively correlated genes were involved in carbon chain binding and phosphorylation (Figure 4D). These data suggested that MPC1 and 2 participate in cellular networks associated with tumor progression in HCC.

Correlation between glycolytic and cholesterogenic gene expression in various cancers

It is very important to identify cancers that have unique metabolic characteristics driven by their mutational environment and organ-specific enzyme expression [28]. To identify the dysregulated routes of metabolism in different cancers, we performed cluster analyses of glycolytic and cholesterogenic gene expression in 26 cancer types (tumor content ≥ 30%) from TCGA (Supplementary Table 2). Gene co-expression pathway-specific profiles and hub genes were identified via network topology analysis in 13 cancer types. Metabolism-related genes were co-expressed in most tumor types. However, due to the differential co-expression of glycolytic and cholesterogenic genes, some genes were only co-expressed in a few cancer types, indicating that certain genes contribute to the tumor metabolic program in a cell-type-specific manner (Figure 5A). Cholesterogenic gene expression correlated positively with poor prognostic gene expression in the Hoshida classification system, not only in HCC (Spearman correlation BH-corrected p < 0.05), but also in a range of other tumor types (Figure 5B). For some cancers, the median cholesterogenic gene expression correlated positively with KRAS expression (cervical squamous cell carcinoma [CESC], glioblastoma multiforme [GBM], kidney renal clear cell carcinoma [KIRC], brain lower grade glioma [LGG], lung squamous cell carcinoma [LUSC], ovarian serous cystadenocarcinoma [OV], pancreatic adenocarcinoma [PAAD], pheochromocytoma and paraganglioma [PCPG], prostate adenocarcinoma [PRAD] and stomach adenocarcinoma [STAD]) or MYC expression (PCPG, STAD and LUSC) (BH-corrected p < 0.05). The expression of MPC1 was significantly upregulated in the cholesterogenic group in KIRC, LGG, lung adenocarcinoma (LUAD) and LUSC, and the expression of MPC2 was significantly upregulated in the cholesterogenic group in PAAD, sarcoma (SARC) and CESC (Figure 5B), similar to our findings in HCC. Further research is needed to explore the mRNA levels of these genes.

Figure 5. The glycolytic and cholesterogenic gene profiles of other cancer types. (A) Heat map depicting that glycolytic and cholesterogenic genes were robustly co-expressed when consensus clusters were applied to each individual cancer type. (B) Top: Bar plots showing the proportions of the metabolic subgroups across the various cancer types. Bottom: The correlation between cholesterogenic gene expression and the expression of Hoshida poor prognostic genes, KRAS, MYC and MPC1/2 in each cancer type. Median glycolytic gene expression correlated positively (BH-adjusted p < 0.05) with basal-like gene expression in all cancer types. The correlation between MPC1/2 expression and glycolytic gene expression was validated using Wilcoxon rank sum tests and BH correction. (C) Kaplan-Meier survival analysis curves depicting the differences in median overall survival across the metabolic subgroups in CESC. (D) Kaplan-Meier survival analysis curves demonstrating the differences in median overall survival in KIRC.

The potential associations between specific metabolic genes and the clinical features of various cancers were analyzed. The survival rates differed significantly among patients in the four metabolic subtypes in CESC (p = 0.001) (Figure 5C) and KIRC (p < 0.001) (Figure 5D). In CESC, the overall survival rate was significantly lower in the mixed subgroup than in the cholesterogenic group. In KIRC, the prognosis was worse in the glycolytic and quiescent groups than in the mixed and cholesterogenic groups. Thus, our model has provided novel insights into the molecular features of multiple cancers and revealed the tumor specificity of metabolic subtype gene expression.

HCC dataset acquisition and processing

The HCC datasets and all corresponding clinical data were downloaded from the data portal of TCGA (Illumina HiSeq Systems; https://cancergenome.nih.gov/) and the ICGC (https://www.icgc.org) [45]. We also downloaded standardized RNA sequence data for all 423 available cases from TCGA and for 237 samples from the ICGC data portal. We used human genome reference sequence GRCh37 from the Genome Reference Consortium [46]. Samples from the ICGC were filtered to exclude those labeled as cell lines, xenografts, metastatic, normal or non-laser microscopy enrichment. After excluding these samples, we downloaded somatic mutational data (CNVs, SNVs and indexes) for all the screened samples.

Analysis of RNA sequence data

We normalized the RNA levels in each sample using the Transcripts Per Million algorithm and log-transformed log₁₀ ((normalized count*1e⁶¹⁴⁵) +1). To identify significantly differentially expressed RNAs, we used a standard screen of a log₂ fold change ≥ 1. All samples were screened to exclude those with tumor contents < 30% [47].

Batch profile estimation and correction

We analyzed data from TCGA and ICGC cohorts, applied a normalization strategy and removed unwanted variation as a correction. Batch corrections were performed in each cohort using a gene localization scale (Figure 6).

Figure 6. Batch correction of queued datasets from TCGA and ICGC. (A) Gene datasets were validated before normalization. (B) Gene datasets were illustrated after normalization.

Classification of metabolic gene subgroups

The Molecular Signatures Database is critical for insightful interpretation of large-scale genomic data [48]. We used the "REACTOME GLYCOLYSIS" (n = 29) to identify the glycolytic subgroup and the "REACTOME CHOLESTEROL BIOSYNTHESIS" (n = 72) to identify the cholesterogenic subgroup. Subtypes were separated based on the consistent clustering of glycolytic and cholesterol-generating genes in ConsensusClusterPlus v1.38, and their statistical significance was assessed with SigClust (parameters: 160 reps=100, p Item=0.8, p Feature=1) [49]. Ward.D2 and the Euclidean distance were respectively used as clustering algorithms and distance measurements (k = 5 for glycolytic and cholesterogenic genes in resected and metastatic LIHC samples [n = 610), the cluster1]), and clusters 1-5 were validated. All patients were divided into four groups based on the median levels of glycolytic and cholesterol-generating genes in each sample: the quiescent group (GLYCOLYSIS ≤ 0, CHOLESTEROL ≤ 0), the glycolytic group (GLYCOLYSIS > 0, CHOLESTEROL ≤ 0), the cholesterogenic group (GLYCOLYSIS ≤ 0, CHOLESTEROL > 0) and the mixed group (GLYCOLYSIS > 0, CHOLESTEROL > 0). For each gene cluster, the ratio of glycolytic genes and cholesterol-generating genes was calculated. Gene clusters containing > 90% CHOLESTEROL genes or > 30% GLYCOLYSIS genes were considered as “core” gene clusters.

Classification of pre-existing HCC subgroups

Consistent clustering was applied to classify samples based on the common tumor subtypes studied by Hoshida et al. [50], Budhu et al. [51] and Chew et al. [52]. The Hoshida subgrouping was based on the 186 genetic characteristics in the original publication, the Budhu subtyping was based on the 17 genetic characteristics in the original paper, and the Chew subtyping was based on the 14 genes in the original study. During each subtype classification, the samples were consistently clustered based on the genes of each classifier, and then semi-automatic subtype assignment was applied.

RNA expression analysis of MPC1/2

We used RNA-seq data to identify gene sets that correlated positively or negatively with MPC1/2 (Spearman correlation analysis), and performed BH correction for multiple test corrections. A significant correlation was established between two gene sets based on an adjusted p value < 0.01. A correlation coefficient of r > 0 was used to establish a positive correlation of a gene with MPC1/2, while a correlation coefficient of r < 0 was used to establish a negative correlation. To determine the pathway enrichment of genes that correlated positively and negatively with MPC1/2, we performed a comprehensive gene set enrichment analysis on the two sets of genes.

Calculation of tumor content

The purity of each tumor sample was limited due to the presence of various immune cells surrounding the tumor cells, as well as other cells in the tumor microenvironment. We estimated the tumor purity of the HCC samples using an R package.

Mutation analysis of HCC genes

We identified and analyzed gene sequences from the human genome assembly GRCh37/hg19 [53]. In order to identify oncogenic molecular events in the HCC metabolic subtypes, we investigated the frequency of SNVs, INDELs and CNVs in commonly mutated HCC genes [54] and explored their relationship with the HCC metabolic subtypes. With respect to tumor ploidy, DNA fragments with copy statuses ≥ 3 and ≤ 1 were considered amplified and deleted, respectively. Based on a previous study [55], we screened HCC copy number events with at least 10 support probes and fragment averages > 0.2 (amplified) or < -0.2 (deleted). The coordinates of the copy number events were mapped to the gene coding region with Bedtools v2.26, and the SNVs and CNVs of each gene were tested with a contingency analysis. The selected genes in each subgroup were tested, and Fisher’s exact test was used to determine whether there was a loss-of-function mutation or copy number amplification/deletion in each subgroup. BH correction was applied to the resulting p values.

Pan-cancer RNA-seq analysis

We downloaded the Transcripts Per Million data of all pan-cancer samples and screened the samples according to cancer type. There were at least 100 samples after the preliminary screening, and there were 26 matched cancer types. Logarithmic conversion, batch correction, grouping and cluster analysis of RNA expression values were performed as described above. For each gene cluster, we calculated the ratio of glycolytic and cholesterol-generating genes. Gene clusters containing > 90% cholesterogenic genes or > 30% glycolytic genes were considered as “core” gene clusters. For cancer types with multiple core clusters in the same gene set, the most homogeneous cluster was considered to be the core cluster. Cancer types that did not have at least 75% homogeneous core glycolytic and cholesterogenic clusters were omitted from further analyses (specifically: bladder urothelial carcinoma, breast invasive carcinoma, colon adenocarcinoma, head and neck squamous cell carcinoma, kidney renal papillary cell carcinoma, acute myeloid leukemia, rectum adenocarcinoma, skin cutaneous melanoma, testicular germ cell tumors, thyroid carcinoma, thymoma and uterine corpus endometrial carcinoma (Supplementary Figure 1)). The metabolic subtypes of patients with each cancer type were determined based on the median values of the respective core glycolytic and cholesterol-generating genes.

Survival analysis of HCC patients

Kaplan-Meier diagrams were generated using the R packages "survival" v.2.4.2 and "survminer" v.0.4.2 [56]. Patients with a total survival of less than one month were removed from the survival analysis.

Statistical analysis

We used SPSS 23.0 software (IBM Corp., Armonk, NY, USA) and GraphPad Prism 7 (San Diego, CA, USA) to analyze the data [57]. Pearson correlation coefficients were used to express the linear correlations between pairs of variables. Kaplan-Meier and log-rank tests were used to analyze patient survival. R Studio was used to obtain the best cut-off value for each gene and survival curve. Two-sided p values < 0.05 were considered significant.