Identification of Energy Metabolism-Related Gene Signatures From scRNA-Seq Data to Predict the Prognosis of Liver Cancer Patients

The increasingly common usage of single-cell sequencing in cancer research enables analysis of tumor development mechanisms from a wider range of perspectives. Metabolic disorders are closely associated with liver cancer development. In recent years, liver cancer has been evaluated from different perspectives and classified into different subtypes to improve targeted treatment strategies. Here, we performed an analysis of liver cancer from the perspective of energy metabolism based on single-cell sequencing data. Single-cell and bulk sequencing data of liver cancer patients were obtained from GEO and TCGA/ICGC databases, respectively. Using the Seurat R package and protocols such as consensus clustering analysis, genes associated with energy metabolism in liver cancer were identified and validated. An energy metabolism-related score (EM score) was established based on five identified genes. Finally, the sensitivity of patients in different scoring groups to different chemotherapeutic agents and immune checkpoint inhibitors was analyzed. Tumor cells from liver cancer patients were found to divide into nine clusters, with cluster 4 having the highest energy metabolism score. Based on the marker genes of this cluster and TCGA database data, the five most stable key genes (ADH4, AKR1B10, CEBPZOS, ENO1, and FOXN2) were identified as energy metabolism-related genes in liver cancer. In addition, drug sensitivity analysis showed that patients in the low EM score group were more sensitive to immune checkpoint inhibitors and chemotherapeutic agents AICAR, metformin, and methotrexate.

Introduction

The new incidence of liver cancer accounted for 4.7% of all malignant tumors, sixth among all malignancies, and deaths due to liver cancer were nearly 830,000, accounting for 8.3% of all cancer deaths (Siegel et al., 2020) according to the global cancer data released by WHO International Agency for Research on Cancer (IARC) in 2020. Since liver tumors are highly heterogeneous, a variety of molecular typing schemes for liver cancer based on gene mutations and the clinical characteristics of liver cancer have been developed previously. For example, Boyault constructed a G1-G6 typing scheme related to the clinical and genetic characteristics of liver cancer based on P53, PIK3CA, hepatitis B virus (HBV) copy numbers and other indicators (Boyault et al., 2007). Hoshida constructed an S1-S3 typing scheme according to the expression levels of MyC, EpCAM, and other genes (Hoshida et al., 2009). There are also differences in prognosis and sensitivity to drug treatment between liver cancer subtypes, which may assist the development of patient-specific treatment strategies.

Tumor metabolism has received increasing attention from researchers in recent years. Reprogramming of energy metabolism is a characteristic feature of tumor cells that promotes rapid cell growth and proliferation. Tumor cells actively ingest glucose through an unusual “anaerobic glycolysis” process (known as the Warburg effect), which not only provides energy for the tumor cell but also allows intermediates to enter side pathways of anabolism to maintain the de novo synthesis of nucleotides, lipids, and amino acids required for cell proliferation (Peter et al., 2019). In addition, organic acids produced due to the Warburg effect cause acidification of the extracellular environment and thus promote ECM degradation and facilitate invasion of tumor cells, resulting in poor prognosis of patients (Jianrong, 2019). However, studies have shown that tumor cells can inhibit the function of tumor-infiltrating immune cells through competitive uptake of nutrients, and metabolites such as lactic acid and cholesterol can further inhibit the function of immune cells, resulting in immune escape (Siska et al., 2020). The liver is the largest metabolic organ in the human body, and metabolic abnormalities are thus closely associated with the occurrence and development of liver cancer. Research on metabolism has also provided new directions for liver cancer diagnosis and treatment (Satriano et al., 2019).

scRNA-seq has been widely used to identify previously unknown tumor subtypes in recent years (Timour and James, 2017). Here, single-cell sequencing data were used to classify liver cancer to obtain genes related to energy metabolism. Finally, a liver cancer cell energy metabolism score (EM score) was constructed to accurately predict the prognosis of patients.

Materials and Methods

Data Collection

Single-cell RNA sequencing data (scRNA-seq) of liver cancer patients were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/, accession number GSE146115). Avereps function in R package Limma was used to average the expression data of the same gene names. Data for liver cancer patients were obtained from the TCGA database (https://portal.gdc.cancer.gov/), including the transcriptome data [RNA-seq; in fragments per kilobase million (FPKM)] and other relevant clinical information. RNA-seq and related clinical data of liver cancer patients were obtained from the ICGC database (https://dcc.icgc.org/), and used as the validation set.

Identification of Liver Cancer Cell Subtypes From scRNA-Seq Data

The scRNA-seq data were analyzed using the Seurat R package. Samples with mitochondrial gene percentages greater than five were excluded. The NormalizeData function was used to standardize the data and extract 1,500 genes with a high coefficient of variation between cells. PCAP was then performed, and p values of each principal component were computed. Principal components with p values less than 0.05 were selected for the subsequent t-SNE analysis, where cells were divided into different clusters. The SingleR R package was used to annotate the cell types of the obtained cell clusters, and a liver cell expression matrix was extracted for further analysis (Aran et al., 2019). The Seurat package was used to perform the same analysis on the liver cell expression matrix, and the liver cancer cells were divided into different clusters. To analyze the differences between different cluster genes, the genes with p values less than 0.05, and log2 |FC| greater than 1 were considered markers of the cluster. Time-series (trajectory) analysis of liver cancer cells was performed using monocle R package to determine the differentiation direction of liver cancer cells.

Determination of Energy Metabolism Score in Different Subtypes of Liver Cancer Cells

Annotated gene information on human energy metabolism pathways was downloaded from the Reactome database (https://reactome.org/). ssGSEA analysis was then performed on the scRNA-seq of liver cancer cells using the annotated gene information, and the scores of each cell in different energy metabolism pathways were obtained. Heat maps were constructed using the pheatmap R package, and cluster assignments were added to reflect the metabolic status of different clusters on the heatmaps. Marker genes of the cluster with the highest metabolic score were then selected for further analysis.

Identification of Key Genes From the TCGA-LIHC RNA-Seq Data

Differential expression analysis of marker genes (Normal tissue vs. Tumor tissue) in the TCGA-LIHC RNA-seq data was performed using the following criteria to determine differentially expressed genes: p < 0.05, and log2 |FC|>1. A univariate Cox analysis was also performed to analyze prognosis-associated marker genes. The intersection of the two gene lists obtained from these procedures was determined, and genes with expression levels inconsistent with the prognosis were removed. TCGA FPKM data were converted to TPM format, and the ConsensusclusterPlus package was used to perform consistent cluster analysis for genes within the intersection. TCGA-LIHC patients were thus divided into different clusters, and K-M survival analysis was used to determine the difference in survival times between the different clusters. Subsequently, 1,000 times lasso analyses were conducted to screen out the most stable genes as key genes, and consistency cluster analysis was conducted once again for the selected genes to determine the survival differences among different clusters.

Validation of Selected Key Genes

PCA was applied to verify the results of the consistency analysis and demonstrate that the identified key genes corresponded to different groups of liver cancer patients. Gene annotations of 25 metabolism-related pathways in the KEGG database (https://www.genome.jp/kegg/) were downloaded, and ssGSEA was used to calculate TCGA-LIHC patient-related scores to analyze the different metabolic pathway scores between clusters. Then, liver cancer patient data was selected from the ICGC database to verify the key gene-based grouping and prognostic correlation.

Determination of High and Low Score Groups Based on Energy Metabolism Score Calculation

According to the previous typing results, if the typing value was found to increase with decreasing gene expression, the gene was classified as “signal A” gene. Vice versa, if gene expression was found to increase with increasing typing values, the gene was classified as “signal B” gene. PCA was used to calculate the energy metabolism score (EM score) of each sample using the following equation:

$\text{EM}\text{Score}={\displaystyle \sum \text{PC}1\text{A}}-{\displaystyle \sum \text{PC}1\text{B}}$

Used the SURv_cutpoint function to filter the best truncation value, and the patients were divided into high-and low-score groups. The difference in survival between the high-and low-score groups was subsequently analyzed.

Analysis of the Association Between Energy Metabolism Score and Immune Cells/Functions

The annotated gene information for 16 types of immune cells and 13 types of immune functions was downloaded from the GSEA database (http://www.gsea-msigdb.org/). The scores of immune cells and immune function of patients in the different evaluation groups were calculated using ssGSEA, and the differences in ssGSEA scores between the high and low EM score groups were analyzed. Additionally, immune cell and stromal cell scores in TCGA-LIHC patients were calculated using the Estimate package, and the correlation between the EM score and the estimated immune score was calculated using Spearman’s correlation coefficient. Expression data of immune checkpoint genes were extracted from patients, and the differences in immune checkpoint levels between the high-and low-score groups were calculated. Statistical significance was set at p < 0.05.

Determination of the Sensitivity of Patients in High and Low Groups to Tumor Treatment Drugs

pRRophetic is an R package created from gene expression and drug sensitivity data of cell lines from the Cancer Genome Project to predict clinical chemotherapy response from the perspective of tumor gene expression levels (Geeleher et al., 2017). Sensitivity to metabolism-related chemotherapeutic agents in patients with liver cancer was assessed using this package, and differences in sensitivity between the groups were analyzed (p < 0.05 was considered statistically significant). The sensitivity differences of TCGA-LIHC patients to PD-1 and CTLA4 inhibitors were analyzed using immunotherapy sensitivity data from the Cancer Immunochromatographic Database (TCIA) (https://tcia.at/).

Results

Data Collection

The scRNA-seq data used in this study was obtained from 3,200 cells of four liver cancer patients. Bulk RNA-Seq was obtained from TCGA and ICGC databases. The TCGA database contains 374 tumor samples and 50 paracancer normal tissue samples. The ICGC database contains gene expression data and clinical information for 231 patients with liver cancer. Patient information related to TCGA and ICGC is shown in Table 1.

TABLE 1

TABLE 1. Clinical characteristics of the liver cancer patients used in this study.

Identification of Liver Cancer Cell Subtypes

The Seurat R package was used to analyze scRNA-seq liver cancer cells. The FindVariableFeatures function of Seurat was used to extract the top 1,500 genes with the largest coefficients of variation for PCA. The p values of the top 20 principal components are shown in Figure 1A. The identified top 20 principal components were incorporated into the t-SNE analysis, where 3,200 cells were divided into 14 clusters (cluster 0–13, Figure 1B). Subsequently, the differentially expressed genes within each cluster were identified (p < 0.05). These differentially expressed genes were defined as marker genes for the corresponding clusters (Supplementary Table S1). To investigate the characteristics of the cells within each cluster, cell type annotations were obtained using the SingleR R package. The results revealed 14 clusters which included four cell types, with cells in cluster 5 containing monocytes, cells in clusters 6, 7, and 9 T cells, cells in cluster 13 NK cells, and cells in the remaining clusters hepatocytes (Figure 1C). GSEA analysis of differentially expressed genes yielded higher and lower metabolism- and immune system-related scores in hepatocyte cells, respectively. On the other hand, the opposite was observed for immune cells (Supplementary Figure S1), which highlighted changes in metabolism of liver cancer cells. To further explore liver cancer cell subtypes, we extracted the gene expression matrix of hepatocyte cells (in total 2,302 cells), and performed the same analysis as before using the Seurat package. PCA revealed a total of 15 principal components with p-values less than 0.05 (Figure 1D). These principal components were incorporated into the t-SNE analysis (Figure 1E). Accordingly, the hepatocytes were divided into nine clusters (0–8), and the marker genes of each cluster are shown in Supplementary Table S2. The monocle R package was used to perform a time-series analysis of the these hepatocytes to explore their differentiation directions. The results showed that hepatocytes gradually followed three differentiation directions, which are termed Branch 1, Branch 2 and Branch 3 (Figure 1F). The differentiation directions of hepatocytes were combined with the cluster information, and calculated for the major components in different Branch clusters. Figure 1G shows that the cells that account for the major component in Branch 1 stemmed from cluster 0 (688/809), those in Branch2 from cluster 4 (352/537), and those in Branch 3 mainly from cluster 5 (138/214).

FIGURE 1

FIGURE 1. Screening of high-energy metabolism hepatocellular carcinoma subtypes using scRNA-seq data. (A) Principal component analysis of all cells and p-value of each principal component. (B) tSNE algorithm divides cells into 14 clusters. (C) Cell types of different clusters. (D) Principal component analysis of hepatocellular carcinoma cells and the p-value of each principal component. (E) tSNE algorithm divides hepatocellular carcinoma cells into 9 clusters. (F) Proposed time series analysis of hepatocellular carcinoma cells. (G) Proposed time series analysis of different clusters of hepatocellular carcinoma cells.

Identification of Hepatocyte Cancer Cell Subtypes With High Energy Metabolism

Based on the annotated genes of the 11 energy metabolism-related pathways downloaded from the Reactome database (see Table 2 for details), ssGSEA analysis was performed to obtain the score of each cell in these pathways. Combining this scoring and cluster information of each cell, the pheatmap R package was used to plot a heat map that shows the scoring status of different clusters. As shown in Figure 2A, clusters 4 and 5 yielded the highest and lowest energy metabolism related scores, respectively. Cluster 0, the largest cluster, yielded a moderate metabolism score. Taken together with the results of the time-series analysis in the previous section, liver cancer cells were concluded to eventually differentiate into three different states of subtypes in terms of energy metabolism: high, medium and low, with Branch 2 showing the highest level of energy metabolism.

TABLE 2

TABLE 2. Metabolic pathways.

FIGURE 2

FIGURE 2. Screening of genes related to high energy metabolism. (A) Energy metabolism scores of differentcluster hepatocellular carcinoma cells. (B) Marker genes associated with survival in high energy metabolism cluster hepatocellular carcinoma cells. (C) Intersection of prognosis-related marker genes and differentially expressed marker genes.

Identification of a TCGA-LIHC High-Energy Metabolic Subtype

Marker gene expression data of cluster 4 were extracted from TCGA-LIHC data, and analyzed to determine differentially expressed genes between cancerous and normal tissues. A total of 101 differentially expressed genes were obtained. Univariate Cox regression analysis was performed in combination with the overall survival time of patients, and a total of 48 genes associated with OS were obtained (Figure 2B). The intersection of differentially expressed and prognosis-related genes was determined. A gene was excluded if its expression level in the tumor was found to be inconsistent with prognosis as described in Methods section. As a result, a total of 45 genes were identified from the initial intersection list (Figure 2C).

Based on the expression data of these 45 genes, a clustering analyses was performed. Based on the cumulative distribution function (CDF), we chose k = 2 as the optimal clustering parameter that yielded consistent results, and thus divided TCGA-LIHC patients into two clusters named cluster A and cluster B (Figure 3A). Survival analysis revealed that the OS of patients in clusters A and B were significantly different (p = 0.002, Figure 3B). Thus, we inferred that these 45 genes could influence the OS of liver cancer patients by affecting their metabolic status. Furthermore, 1000 lasso analysis also identified five genes (Figure 3C): ADH4, AKR1B10, CEBPZOS, ENO1, and FOXN2, as the most stable prognosis-related genes associated with energy metabolism in liver cancer. The correlation results among the five genes were shown in Supplementary Table S5 and Supplementary Figure S2 (Correlations between genes greater than 0.2 were shown in Supplementary Figure S2). The immunohistochemical results of the expression of the above five genes in liver cancer and normal tissues in HPA database were shown in Supplementary Figure S3.

FIGURE 3

FIGURE 3. Identification of TCGA-LIHC energy metabolism subtypes. (A)Typing of TCGA-LIHC based on 45 intersecting genes. (B) Survival of TCGA-LIHC patients within different subtypes based on 45 intersecting genes. (C) Lasso analysis to screen for the most stable energy metabolism-related genes. (D) TCGA-LIHC typing based on the five most stable energy metabolism-related genes. (E) Survival of TCGA-LIHC patients within different subtypes based on five key genes.

To explore the significance of these five genes, another clustering analysis was performed to divide TCGA-LIHC patients into two clusters, cluster A and cluster B (Figure 3D). Survival analysis revealed that the OS of cluster B was significantly lower than that of cluster A, and the difference was statistically significant (p = 0.013, Figure 3E).

Validation of the Identified Key Genes

ssGSEA results revealed that multiple higher metabolic pathway scores were higher in cluster B compared to cluster A (Figure 4A). The same analysis was performed on the liver cancer data obtained from the ICGC database to test whether choice of the dataset affected our conclusions. Based on ICGC data, the same five genes were able to classify liver cancer patients into two clusters (Figure 4B), and the survival analysis between the two groups showed a statistically significant difference with respect to OS (Figure 4C).

FIGURE 4

FIGURE 4. Validation of TCGA-LIHC typing results based on 5 key genes. (A) Differences in metabolism-related pathway scores among different subtypes. (B) Staging results of patients with hepatocellular carcinoma in ICGC database. (C) Survival of patients within different subtypes in the ICGC database.

Determination of Energy Metabolism Scores

To determine an energy metabolism score for TCGA-LIHC patients, we used principal component analysis (PCA) to calculate PC1 for genes within gene tags A and B, and calculated the sum of PC1 for both tags A and B (sPC1A and sPC1B), respectively. Subsequently, the difference between sPC1A and sPC1B was used as the energy metabolism score (EM score). Using the best cut-off value, TCGA cohort patients were divided into groups of high or low EM scores. Survival analysis showed that patients in the high score group had significantly worse survival than those in the low score group (p < 0.001) (Figure 5A). Thus, the EM score was successfully used to accurately classify patients based on overall survival.

FIGURE 5

FIGURE 5. Construction of EM score and its relationship with immunity. (A) Survival of patients within different EM score groups. (B) Correlation between EM score and Estimate ImmuneScore. (C) Differences in immune function scores of patients within different EM score groups. (D) Differences in immune cell scores of patients within different EM score groups.

The Relationship Between Energy Metabolism Score and Immunity

There is substantial evidence in literature on competition between tumor and immune system cells for nutrients, and inhibition of immune cell function by tumor cell metabolites (Xing et al., 2015). To explore the relationship between tumor energy metabolism and immunity in liver cancer, we analyzed 16 immune cell and 13 immune function scores in TCGA-LIHC patients by ssGSEA, and used the Wilcoxon test to determine the significance of the differences in scores of high and low EM Score groups. A total of 11 immune cell and 12 immune function scores were found to differ between the high and low EM Score groups (Figures 5B,C), with the majority being higher in the low EM Score group. Estimate analysis showed a decreasing trend of ImmuneScore (R = -0.17, p = 0.001) in TCGA-LIHC patients with increasing EM Score, and the results were statistically significant (Figure 5D). In addition, 25 of the 38 immune checkpoints differed between the high and low EM score groups, with most being highly expressed in the high EM score group (Figure 6). These results suggest that EM score can be used as an indicator of an influence on immunity in liver cancer patients.

FIGURE 6

FIGURE 6. Differences in immune checkpoint levels between different EM score groups.

The Role of Energy Metabolism Score in Predicting Drug Sensitivity in Cancer Therapy

Using the pRRophetic R package, the sensitivity (IC50) of patients to 138 chemotherapeutic agents, including AICAR, metformin, and other metabolism-related anticancer drugs, were determined. Accordingly, liver cancer patients in the low EM score group had higher sensitivity (lower IC50) to three metabolism-related anticancer drugs: AICAR (p = 2.22e-16), metformin (p = 4.3e-15), and methotrexate (p = 0.001) (Figures 7A–C). The TICA database data showed that liver cancer patients within the low EM score group had higher sensitivity to CTLA4 inhibitors (p = 0.00015), PD-1 inhibitors (p = 0.0021), and PD-1 inhibitors combined with CTLA4 inhibitors (p = 0.0066) (Figures 7D–F). Taken together, these data suggest that the EM score may be associated with response to chemotherapeutic agents in liver cancer patients.

FIGURE 7

FIGURE 7. Relationship between EM score and sensitivity of chemotherapy drugs and immune checkpoint inhibitors. (A–C) Differences in sensitivity to chemotherapeutic agents AICAR, Metformin and Methotrexate among patients in different EM score groups. (D–F) Differences in sensitivity to PD1 inhibitors, CTLA-4 inhibitors and the combination of these two inhibitors in different EM score groups.

Discussion

Liver cancer is a disease that seriously threatens human health. The disease ranked sixth in new cases and third in deaths among all malignant tumors in 2020. Therefore, in-depth investigations of the mechanism of liver cancer development has been prioritized by cancer researchers. The liver is the largest metabolic organ in the human body, and alterations in its metabolic profile are often closely related to the development of a primary tumor (Pope et al., 2019). Liver cancer cells in a hypermetabolic state are not only able to compete with normal cells for nutrients, but their metabolites can also affect the function of immune cells and promote tumor metastasis (Hanahan and Weinberg, 2011; Wu, 2017). Here, we focused on the metabolic status of different subtypes of liver cancer, and constructed and validated an energy metabolism (EM) score to quantify the metabolic status of liver cancer patients based on scRNA sequencing, TCGA-LIHC, and ICGC database data. Our findings showed that the EM score could be used to accurately predict the prognosis of liver cancer patients, and has potential to assist the selection of drugs for patient-specific treatment options.

Many studies have previously shown that the development of liver cancer is accompanied by disturbances in energy metabolism, which not only suppresses immune response to tumors, but also promotes tissue invasion and metastasis. Therefore, it is highly relevant to study the genes related to energy metabolism in liver cancer, and developed a score to characterize the status of energy metabolism in liver cancer patients. Here, we obtained liver cancer single-cell sequencing data from the GEO database, clustered and annotated the cells, and used ssGSEA to identify a subgroup of liver cancer cells with a high energy metabolism level. Based on the marker genes in this subgroup, we used univariate and Lasso Cox regression analysis to screen the TCGA database data for the most accurate prognostic biomarkers to establish an energy metabolism-related gene signature. As a result, we obtained a signature containing five genes: ADH4, AKR1B10, CEBPZOS, ENO1, and FOXN2. ADH4 encodes the class II ethanol dehydrogenase 4 subunit, which is a dimer consisting of two subunits with high oxidative activity towards long-chain fatty alcohols and aromatic alcohols. ADH4 has been associated with various disorders such as alcohol dependence (Chowdhury et al., 2020). Differential expression of ADH4 between tumor and normal tissues in patients with non-small cell lung cancer has been observed, and can be used as a prognostic marker in patients with esophageal cancer (Peng et al., 2018). Wei et al. showed that ADH4 was expressed at lower levels in the tumor tissues of liver cancer patients, and patients with low ADH4 expression levels had significantly shorter survival times (Wei et al., 2012).

AKR1B10 encodes a member of the aldehyde/ketone reductase superfamily, which is overexpressed in various solid tumors, and a potential diagnostic marker for tumors (Endo et al., 2021). In breast cancer, this gene enhances fatty acid utilization by tumor cells and activates ERK, Wnt, and other pathways to promote metastasis (Qu et al., 2020). In liver cancer, AKR1B10 mediates the proliferation of liver cancer cells through sphingosine-1-phosphate (Jin et al., 2016). In addition, AKR1B10 induces cellular resistance to erythromycin and nortriptyline by reducing the C13 ketone moiety, which leads to poor patient prognosis (Zhong et al., 2011).

CEBPZOS is the antisense chain of CEBPZ, and is expressed in various tissues. Only several studies have been conducted on CEBPZOS, yet it has been shown to be a novel mitochondrial protein (Hung et al., 2014). Due to their key role in intracellular energy generation, mitochondria have a great impact on cell growth, apoptosis, and maintenance of redox homeostasis. In tumors, oxidative phosphorylation of the mitochondria is very active, and mitochondria are thus primary sites where the Warburg effect is strong (Wallace, 2012). The role of CEBPZOS as a mitochondrial protein may also affect the energy metabolism of tumor cells in liver cancer, and thereby affect tumor development as well.

ENO1 encodes α-enolase, which is one of the three enolase isozymes found in mammals (Zhang et al., 2020). ENO1 is overexpressed in a variety of GI tumors, including colon, pancreatic, and gastric cancers (Hang et al., 2018), and promotes tumor development (Cheng et al., 2019; Xu et al., 2019). In addition, ENO1 is highly expressed in liver cancer, promotes glucose uptake and lactate production by tumor cells, and is involved in tumor cell division, proliferation, apoptosis, metastasis, immunomodulation, and chemoresistance (Qiao et al., 2018; Jiang et al., 2020).

FOXN2 plays an important role in the development of several tumors (Ye and Duan, 2019), including breast and cervical cancers (Cui et al., 2019). FOXN2 has been shown to be associated with tumor cell proliferation and invasion in liver cancer (Liu et al., 2021).

Based on these five genes, patients were divided into two groups, cluster A and cluster B. The prognosis of the two groups was found to be significantly different from each other. The energy metabolism of the patients was also analyzed using ssGSEA, and the results yielded significant differences in the energy metabolism-related pathway scores between the two groups. We conclude that these five genes are associated with the energy metabolism status of liver cancer patients, and can affect their overall survival. In addition, an external validation set from the ICGC database validated the accuracy of gene signatures consisting of these five energy metabolism-related genes for patient prognosis prediction.

Considering the variability in metabolic status among patients, it is important to quantify the status of energy metabolism of each patient. Identification of tumor subtypes based on different biomarkers has been used to improve the accuracy of patient prognosis for various tumors, such as breast cancer (Callari et al., 2016). Here, we used the above mentioned key genes as potential “subtype biomarkers,” and established an EM score to quantify the energy metabolism status of each sample. Survival analysis showed that higher EM scores were associated with lower survival rates.

Tumors are known to compete for glucose with immune system cells, and thus inhibit immune cell function. The glycolytic activities of tumor cells may also limit glucose consumption by immune cells, and thereby lead to immune escape (Fernández et al., 2019). In addition, metabolites produced by tumor cells can have profound effects on immune cells in the tumor microenvironment (TME). For example, aberrant glycolysis in tumor cells leads to the production of large amounts of lactic acid, resulting in an acidic TME. The low pH of the TME was shown to favor more aggressive tumor cells, and suppress immune response to tumors (Terrén et al., 2019). Arachidonic acid can also suppress immune responses by inhibiting Th1 differentiation, NK cell function, and T-cell activation (Pidgeon et al., 2007; Juman et al., 2013). Thus, we analyzed the differences in immune cell and immune function scores between the high and low EM score groups using ssGSEA to further explore the relationship between EM scores and immunity of liver cancer patients. In total, 11 immune cell and 12 immune function scores were found to differ between high and low EM score groups. Among them, B, T, and NK cells, which have anti-tumor effects, yielded higher scores in the low EM score group. However, the scores of immune checkpoints and immune functions with pro-cancer effects, such as CCR, were higher in the high EM score group. In addition, the results of the Estimate analysis showed that the immune and stromal scores of TCGA-LIHC patients decreased with increasing EM scores. In addition, the expression levels of most immune checkpoints were higher in liver cancer patients with high EM scores. These results suggest that the EM score reflects the status of the immune systems of liver cancer patients, with high EM scores representing lower immunity levels.

We also found that liver cancer patients in the low EM score group were more sensitive to three metabolism-related anticancer drugs, AICAR, metformin, and methotrexate, based on the results of the pRRophetic analysis. AICAR is an activator of AMP-activated protein kinase (AMPK) which can permeabilize cell membranes. Animal experiments have shown that AICAR significantly inhibits fatty acid and sterol synthesis in mouse hepatocytes (Gao et al., 2018), whereas cellular experiments have shown that AICAR inhibits the proliferation of liver cancer cells, and induces cell cycle arrest at the G1-S checkpoint (Cheng et al., 2014). Metformin is a drug related to glucose metabolism, and many studies have shown that metformin is not only effective against diabetes but also has therapeutic effects on a variety of tumors, including lung, pancreatic and breast cancers (Rizos and Elisaf, 2013). Wang et al. showed that in esophageal cancer, low doses of metformin reprogrammed the tumor immune microenvironment (TIME) in an anti-cancer direction by increasing the proportion of CD8⁺ cells (Wang et al., 2020). The combination of aloin and metformin was found to inhibit the growth and invasion of liver cancer through the PI3K/AKT/mTOR pathway and induce apoptosis and autophagy, thereby exerting anti-tumor effects (Sun et al., 2020). In addition, metformin and dichloroacetate inhibited the proliferation of liver cancer cells by inhibiting mTOR complex 1. Methotrexate is an anti-folate antitumor agent that inhibits the growth and multiplication of tumor cells by hindering their synthesis, mainly through the inhibition of dihydrofolate reductase (Rosh, 2022). Folic acid deficiency was found to cause defects in oxidative phosphorylation in human cells due to impaired mitochondrial translation, and thus affect energy metabolism (Morscher et al., 2018). Methotrexate has also been used as a therapeutic agent against lymphoma, osteosarcoma, and many other tumor types (Jaffe, 2009; Gaynon et al., 2010). In this study, patients with low EM scores showed higher sensitivity to the three chemotherapeutic agents mentioned above, suggesting that treatment with these three chemotherapeutic agents may be effective for patients with this type of liver cancer.

Kimiteru et al. showed that a CLTA-4 inhibitor (ipilimumab) was less effective in melanoma patients in a hypermetabolic state, and lead to shorter survival after treatment (Ito et al., 2018). In non-small cell lung cancer, patients with a larger total metabolic tumor volume (TMTV) have a shorter survival time after receiving immune checkpoint inhibitors (Seban et al., 2020). This suggests that EM scores may be related to the efficacy of immune checkpoint inhibitors in liver cancer treatment. Therefore, we analyzed the differences in sensitivity to PD-1 and CTLA-4 inhibitors between patients in the high and low EM score groups based on relevant data within the TCIA database. We found that patients in the low EM score group had higher sensitivity to PD-1 inhibitors, CTLA-4 inhibitors, and PD-1 inhibitors combined with CTLA4 inhibitors. Hence, immune checkpoint inhibition therapy may provide with a treatment option for patients with high EM scores.

Conculsion

In this study, we identified a class of liver cancer cells with high energy metabolism from single-cell sequencing data, and determined five key genes related to energy metabolism on marker genes of this class of liver cancer cells combined with bulk sequencing data. We then constructed an EM score based on expression levels of these five genes to predict patient survival. In addition, we found that these five genes may be related to the sensitivity of liver cancer patients to chemotherapeutic drugs and immune checkpoint inhibitors. On the other hand, our findings are based on the exploration of public databases, and experimental validation is necessary. Moreover, specific mechanisms of regulation of energy metabolism in hepatocellular carcinoma cells by these signature genes need to be explored as well.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

MS and LG conceived the study. BX conducted all statistical analyses. ZP and YA reviewed relevant literature and drafted the manuscript. GY and XY provided guidance to the study. All authors read and approved the final manuscript.

Funding

This work was supported by grants from the National Key R&D Program of China (No. 2016YFC1303601).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.858336/full#supplementary-material

Abbreviations

CDF, cumulative distribution function; CTLA-4, cytotoxic T-lymphocyte–associated antigen 4; EM Score, energy metabolism score; GEO, gene expression omnibus; HBV, hepatitis B virus; PD-1, programmed cell death protein 1; PD-L1, programmed death-ligand 1; ssGSEA, single-sample gene set enrichment analysis; SPC, sum of PCI; TCGA, the cancer genome atlas; TCIA, the cancer immunome atlas; TML, tumor mutation load; IARC, international agency for research on cancer; scRNA-seq: single-cell RNA sequencing data; TIME, tumor immune microenvironment; TMTV, total metabolic tumor volume.

References

Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V., Hsu, A., et al. (2019). Reference-based Analysis of Lung Single-Cell Sequencing Reveals a Transitional Profibrotic Macrophage. Nat. Immunol. 20, 163–172. doi:10.1038/s41590-018-0276-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyault, S., Rickman, D. S., de Reyniès, A., Balabaud, C., Rebouissou, S., Jeannot, E., et al. (2007). Transcriptome Classification of HCC Is Related to Gene Alterations and to New Therapeutic Targets. Hepatology (Baltimore, Md.) 45 (1), 42–52. doi:10.1002/hep.21467

PubMed Abstract | CrossRef Full Text | Google Scholar

Callari, M., Cappelletti, V., D'Aiuto, F., Musella, V., Lembo, A., Petel, F., et al. (2016). Subtype-Specific Metagene-Based Prediction of Outcome after Neoadjuvant and Adjuvant Treatment in Breast Cancer. Clin. Cancer Res. 22, 337–345. doi:10.1158/1078-0432.ccr-15-0757

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, J., Huang, T., Li, Y., Guo, Y., Zhu, Y., Wang, Q., et al. (2014). AMP-activated Protein Kinase Suppresses the In Vitro and In Vivo Proliferation of Hepatocellular Carcinoma. Plos One 9, e93256. doi:10.1371/journal.pone.0093256

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, Z., Shao, X., Xu, M., Zhou, C., and Wang, J. (2019). ENO1 Acts as a Prognostic Biomarker Candidate and Promotes Tumor Growth and Migration Ability through the Regulation of Rab1A in Colorectal Cancer. Cmar 11, 9969–9978. doi:10.2147/cmar.s226429

PubMed Abstract | CrossRef Full Text | Google Scholar

Chowdhury, N. P., Moon, J., and Müller, V. (2020). Adh4, an Alcohol Dehydrogenase Controls Alcohol Formation within Bacterial Microcompartments in the Acetogenic Bacterium Acetobacterium Woodii. Environ. Microbiol. 23 1, 499–511. doi:10.1111/1462-2920.15340

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, L., Nai, M., Zhang, K., Li, L., and Li, R. (2019). lncRNA WT1-AS Inhibits the Aggressiveness of Cervical Cancer Cell via Regulating P53 Expression via Sponging miR-330-5p. Cmar 11, 651–667. doi:10.2147/cmar.s176525

PubMed Abstract | CrossRef Full Text | Google Scholar

Endo, S., Matsunaga, T., and Nishinaka, T. (2021). The Role of AKR1B10 in Physiology and Pathophysiology. Metabolites 11. doi:10.3390/metabo11060332

CrossRef Full Text | Google Scholar

Gao, J., Xiong, R., Xiong, D., Zhao, W., Zhang, S., Yin, T., et al. (2018). The Adenosine Monophosphate (AMP) Analog, 5-Aminoimidazole-4-Carboxamide Ribonucleotide (AICAR) Inhibits Hepatosteatosis and Liver Tumorigenesis in a High-Fat Diet Murine Model Treated with Diethylnitrosamine (DEN). Med. Sci. Monitor Int. Med. J. Exp. Clin. Res. 24, 8533–8543. doi:10.12659/msm.910544

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaynon, P. S., Angiolillo, A. L., Carroll, W. L., Nachman, J. B., Trigg, M. E., Sather, H. N., et al. (2010). Long-term Results of the Children's Cancer Group Studies for Childhood Acute Lymphoblastic Leukemia 1983-2002: A Children's Oncology Group Report. Leukemia 24, 285–297. doi:10.1038/leu.2009.262

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2017). pRRophetic: an R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels. PLoS ONE 9, e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of Cancer: the Next Generation. Cell 144, 646–674. doi:10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoshida, Y., Nijman, S. M. B., Kobayashi, M., Chan, J. A., Brunet, J-P., Chiang, D. Y., et al. (2009). Integrative Transcriptome Analysis Reveals Common Molecular Subclasses of Human Hepatocellular Carcinoma. Cancer Res. 69 (18), 7385–7392. doi:10.1158/0008-5472.CAN-09-1089

PubMed Abstract | CrossRef Full Text | Google Scholar

Hung, V., Zou, P., Rhee, H.-W., Udeshi, N. D., Cracan, V., Svinkina, T., et al. (2014). Proteomic Mapping of the Human Mitochondrial Intermembrane Space in Live Cells via Ratiometric APEX Tagging. Mol. Cel 55, 332–341. doi:10.1016/j.molcel.2014.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ito, K., Schöder, H., Teng, R., Humm, J. L., Wolchok, J. D., and Weber, W. A. (2018). Prognostic Value of Baseline Metabolic Tumor Volume Measured on 18F-Fluorodeoxyglucose Positron Emission Tomography/computed Tomography in Melanoma Patients Treated with Ipilimumab Therapy. Eur. J. Nucl. Med. Mol. Imaging. 46 4, 930–939. doi:10.1007/s00259-018-4211-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaffe, N. (2009). Osteosarcoma: Review of the Past, Impact on the Future. The American Experience. Cancer Treat. Res. 152, 239–262. doi:10.1007/978-1-4419-0284-9_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, K., Dong, C., Yin, Z., Li, R., Mao, J., Wang, C., et al. (2020). Exosome-derived ENO1 Regulates Integrin α6β4 Expression and Promotes Hepatocellular Carcinoma Growth and Metastasis. Cell Death Dis 11, 972. doi:10.1038/s41419-020-03179-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jianrong, L. (2019). The Warburg Metabolism Fuels Tumor Metastasis. Cancer metastasis Rev. 38. doi:10.1007/s10555-019-09794-5

CrossRef Full Text | Google Scholar

Jin, J., Liao, W., Yao, W., Zhu, R., Li, Y., and He, S. (2016). Aldo-keto Reductase Family 1 Member B 10 Mediates Liver Cancer Cell Proliferation through Sphingosine-1-Phosphate. Sci. Rep. 6, 22746. doi:10.1038/srep22746

PubMed Abstract | CrossRef Full Text | Google Scholar

Juman, S., Hashimoto, M., Katakura, M., Inoue, T., Tanabe, Y., Arita, M., et al. (2013). Effects of Long-Term Oral Administration of Arachidonic Acid and Docosahexaenoic Acid on the Immune Functions of Young Rats. Nutrients 5, 1949–1961. doi:10.3390/nu5061949

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X. H., Liu, L. P., Xu, X. M., Hua, M., Kang, Q., Li, A., et al. (2021). FOXN2 Suppresses the Proliferation and Invasion of Human Hepatocellular Carcinoma Cells. Eur. Rev. Med. Pharmacol. Sci. 25, 731–737. doi:10.26355/eurrev_202101_24634

PubMed Abstract | CrossRef Full Text | Google Scholar

Morscher, R. J., Ducker, G. S., Hsin-Jung Li, S., and Mayer, J. A. (2018). Mitochondrial Translation Requires Folate-dependent tRNA Methylation. Nature 554 7690, 128–132. doi:10.1038/nature25460

PubMed Abstract | CrossRef Full Text | Google Scholar

Peter, V., Heinz, S., and Arnulf, M. (2019). The Warburg Effect: Essential Part of Metabolic Reprogramming and central Contributor to Cancer Progression. Int. J. Radiat. Biol. 95.

PubMed Abstract | Google Scholar

Pidgeon, G. P., Lysaght, J., Krishnamoorthy, S., Reynolds, J. V., O’Byrne, K., Nie, D., et al. (2007). Lipoxygenase Metabolism: Roles in Tumor Progression and Survival. Cancer Metastasis Rev. 26, 503–524. doi:10.1007/s10555-007-9098-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Piñeiro Fernández, J., Luddy, K. A., Harmon, C., and O'Farrelly, C. (2019). Hepatic Tumor Microenvironments and Effects on NK Cell Phenotype and Function. Int. J. Mol. Sci. 20, 4131. doi:10.3390/ijms20174131

CrossRef Full Text | Google Scholar

Pope, E. D., Kimbrough, E. O., Vemireddy, L. P., Surapaneni, P. K., Copland, J. A., and Mody, K. (2019). Aberrant Lipid Metabolism as a Therapeutic Target in Liver Cancer. Expert Opin. Ther. Targets 23, 473–483. doi:10.1080/14728222.2019.1615883

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiao, G., Xu, H., Li, C., Li, X., Farooqi, A. A., Zhao, Y., et al. (2018). Granulin A Synergizes with Cisplatin to Inhibit the Growth of Human Hepatocellular Carcinoma. Int. J. Mol. Sci. 19. doi:10.3390/ijms19103060

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, J., Li, J., He, R., Liu, X., and Luo, D. 2020. AKR1B10 Promotes Breast Cancer Cell Proliferation and Migration via the PI3K/AKT/NF-κB Signaling Pathway. Cel Biosci. 11, 163.

Google Scholar

Rizos, C. V., and Elisaf, M. S. (2013). Metformin and Cancer. Eur. J. Pharmacol. 705, 96–108. doi:10.1016/j.ejphar.2013.02.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosh, J. (2022). Methotrexatein StatPearls [Internet] (Treasure Island (FL): StatPearls Publishing).

Google Scholar

Satriano, L., Lewinska, M., Rodrigues, P. M., Banales, J. M., and Andersen, J. B. (2019). Metabolic Rearrangements in Primary Liver Cancers: Cause and Consequences. Nat. Rev. Gastroenterol. Hepatol. 16, 748–766. doi:10.1038/s41575-019-0217-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Seban, R. D., Mezquita, L., Berenbaum, A., Dercle, L., Botticella, A., Le Pechoux, C., et al. (2020). Baseline Metabolic Tumor burden on FDG PET/CT Scans Predicts Outcome in Advanced NSCLC Patients Treated with Immune Checkpoint Inhibitors. Eur. J. Nucl. Med. Mol. Imaging 47, 1147–1157. doi:10.1007/s00259-019-04615-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer Statistics, 2020. CA: A Cancer J. Clinicians 70. doi:10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

Siska, P. J., Singer, K., Evert, K., Renner, K., and Kreutz, M. (2020). The Immunological Warburg Effect: Can a Metabolic-Tumor-Stroma Score (MeTS) Guide Cancer Immunotherapy? Immunol. Rev. 295, 187–202. doi:10.1111/imr.12846

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, R., Zhai, R., Ma, C., and Miao, W. (2020). Combination of Aloin and Metformin Enhances the Antitumor Effect by Inhibiting the Growth and Invasion and Inducing Apoptosis and Autophagy in Hepatocellular Carcinoma through PI3K/AKT/mTOR Pathway. Cancer Med. 9, 1141–1151. doi:10.1002/cam4.2723

PubMed Abstract | CrossRef Full Text | Google Scholar

Terrén, I., Orrantia, A., VITALLé, J., Zenarruzabeitia, O., and Borrego, F. (2019). NK Cell Metabolism and Tumor Microenvironment. Front. Immunol. 10, 2278. doi:10.3389/fimmu.2019.02278

PubMed Abstract | CrossRef Full Text | Google Scholar

Timour, B., and James, H. (2017). Unravelling Biology and Shifting Paradigms in Cancer with Single-Cell Sequencing. Nat. Rev. Cancer 17. doi:10.1038/nrc.2017.58

CrossRef Full Text | Google Scholar

Wallace, D. C. (2012). Mitochondria and Cancer. Nat. Rev. Cancer 12, 685–698. doi:10.1038/nrc3365

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P., Zhang, L., Huang, C., Huang, P., and Zhang, J. (2018). Distinct Prognostic Values of Alcohol Dehydrogenase Family Members for Non-small Cell Lung Cancer. Med. Sci. Monit. 24, 3578–3590. doi:10.12659/MSM.910026

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Lin, Y., Xiong, X., Wang, L., and Zhang, H. (2020). Low-Dose Metformin Reprograms the Tumor Immune Microenvironment in Human Esophageal Cancer: Results of a Phase II Clinical Trial. Clin. Cancer Res. 26, 4921–4932. doi:10.1158/1078-0432.CCR-20-0113

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, R.-R., Zhang, M.-Y., Rao, H.-L., Pu, H.-Y., Zhang, H.-Z., and Wang, H.-Y. (2012). Identification of ADH4 as a Novel and Potential Prognostic Marker in Hepatocellular Carcinoma. Med. Oncol. 29, 2737–2743. doi:10.1007/s12032-011-0126-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, D. (2017). Innate and Adaptive Immune Cell Metabolism in Tumor Microenvironment. Adv. Exp. Med. Biol. 1011, 211–223. doi:10.1007/978-94-024-1170-6_7

PubMed Abstract | CrossRef Full Text | Google Scholar

Xing, Y., Zhao, S., Zhou, B. P., and Mi, J. (2015). Metabolic Reprogramming of the Tumour Microenvironment. FEBS J. 282, 3892–3898. doi:10.1111/febs.13402

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, X., Chen, B., Zhu, S., Zhang, J., He, X., Cao, G., et al. (2019). Hyperglycemia Promotes Snail-Induced Epithelial-Mesenchymal Transition of Gastric Cancer via Activating ENO1 Expression. Cancer Cel Int 19, 344. doi:10.1186/s12935-019-1075-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, H., and Duan, M. (2019). FOXN2 Is Downregulated in Breast Cancer and Regulates Migration, Invasion, and Epithelial- Mesenchymal Transition through Regulation of SLUG. Cancer Manag. Res. 11, 525–535. doi:10.2147/CMAR.S176938

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W.-Y., Jin, Z.-D., Liu, F., Yuan, H.-H., and Jiang, B. (2018). Antitumor Activity of Intratumoral Ethanol Injection in an Orthotopic Pancreatic Cancer Cell Mouse Xenograft Model. Gastroenterol. Res. Pract. 2018 (-2-1), 1–8. doi:10.1155/2018/7149565

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Gao, J., Shen, F., Ma, X., Wang, Z., Hou, X., et al. (2020). Cinnamaldehyde Changes the Dynamic Balance of Glucose Metabolism by Targeting ENO1. Life Sci. 258, 118151. doi:10.1016/j.lfs.2020.118151

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, L., Shen, H., Huang, C., Jing, H., and Cao, D. (2011). AKR1B10 Induces Cell Resistance to Daunorubicin and Idarubicin by Reducing C13 Ketonic Group. Toxicol. Appl. Pharmacol. 255, 40–47. doi:10.1016/j.taap.2011.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar