Research Paper Volume 13, Issue 9 pp 13318—13332

Development of a prognostic model for hepatocellular carcinoma using genes involved in aerobic respiration

Jiawei Rao1,2,3, *,#, , Xukun Wu4, *,#, , Xiaozhuan Zhou5, , Ronghai Deng1,2,3, , Yi Ma1,2,3, ,

  • 1 Organ Transplant Center, The First Affiliated Hospital, Sun Yat-Sen University, Guangzhou 510080, China
  • 2 Guangdong Provincial Key Laboratory of Organ Donation and Transplant Immunology, The First Affiliated Hospital, Sun Yat-Sen University, Guangzhou 510080, China
  • 3 Guangdong Provincial International Cooperation Base of Science and Technology (Organ Transplantation), The First Affiliated Hospital, Sun Yat-Sen University, Guangzhou 510080, China
  • 4 Department of Hepatology Surgery, The First Affiliated Hospital of Sun Yat-Sen University, Guangzhou, China
  • 5 Department of Gastroenterology, The First Affiliated Hospital of Sun Yat-Sen University, Guangzhou, China
* Equal contribution
# Co-first authors

Received: August 17, 2020       Accepted: October 3, 2020       Published: April 26, 2021      

https://doi.org/10.18632/aging.203021
How to Cite

Copyright: © 2021 Rao et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Hepatocellular carcinoma (HCC) is the second leading cause of cancer-related death worldwide. Currently, recent risk stratification has only focused on liver function and tumor characteristics. Thus, the purpose of this study was to develop a prognostic model based on genes involved in aerobic respiration. Matched tumor and normal tissues from TCGA and ICGC cohorts were analyzed to identify 15 overlapping differential expressed genes. Cox univariate analysis of the 15 genes in the TCGA cohort revealed they were all associated with disease-specific survival (DSS) in HCC patients. Using LASSO estimation and the optimal value for penalization coefficient lambda 12 genes were selected for the prognostic model, and then HCC patients in the TCGA cohort were dichotomized into low-risk and high-risk groups. Univariate and multivariate Cox analysis demonstrated patients in low-risk group had better survival. Validation of the risk score model with the ICGC cohort produces results consistent with those of the TCGA cohort. In conclusion, this study developed and validated a prognostic model of HCC through a comprehensive analysis of genes involved in aerobic respiration. This model may help develop personalized treatments for patients with HCC.

Introduction

Hepatocellular carcinoma (HCC) is the second leading cause of cancer-related death worldwide, and accounts for around 800,000 deaths yearly [3]. HCC is a heterogeneous disease with a number of risk factors including hepatitis B virus (HBV) infection, hepatitis C virus (HCV) infection, obesity, and alcohol abuse. Treatments for HCC include hepatectomy, liver transplantation, radiotherapy, transcatheter arterial chemoembolization (TACE), and chemotherapy [4]; however, the survival of patients with advanced HCC is poor. The assessment of prognosis is important for patients with HCC, as an adequate assessment can prevent low-risk patients from receiving unnecessary treatments, and high-risk patients from relapse or metastasis due to inadequate treatment. Current well-known staging strategies such as tumor, node, metastasis (TNM) and the Barcelona Clinic Liver Cancer (BCLC) system only focus on tumor diameter, tumor metastasis, and liver function. However, the important role of metabolic subtypes is ignored. Thus, the addition of genetic factors to traditional prognostic models may improve risk assessment.

Oxidative phosphorylation (OXPHOS) and the tricarboxylic acid cycle (TCA) in mitochondria are important methods for the production of ATP, and OXPHOS accounts for 80% of the ATP production in mammalian cells. In the 1920s, Otto Warburg observed that cancer cells have high glucose consumption and high lactate production, even in the presence of abundant oxygen, indicating a switch from OXPHOS to glycolysis [1]. The switch may be due to the faster rate of ATP production by glycolysis, which supports the rapid proliferation of malignant cells [2]. To facilitate glycolytic efficiency, the activity of pyruvate dehydrogenase (PDH) is usually suppressed which prevents pyruvate from being transformed to acetyl coenzyme A (acetyl-CoA) for subsequent aerobic mitochondrial metabolism [5]. Thus, the TCA cycle in cancer cells is bypassed. Defects of OXPHOS might be associated with decrease levels of complex I or autophagic degradation, which is caused by oncogenic K-Ras activation [6, 7]. Lee et al. reported that exogenous lactate released by cancer cells can inhibit OXPHOS in normal cells by reducing MRPL13 expression, thus forming a vicious circle [5].

A number of studies have confirmed the relations between OXPHOS, the TCA cycle, and HCC. A study showed that inactivation of OXPHOS, rather than increased activation of glycolysis, was responsible for both inherent and acquired sorafenib resistance of HCC cells [8]. Another study showed that citramalic acid in tumor tissues of HCC patients who were HBV surface (HBsAg) positive was significantly lower than in patients who were HBsAg negative, which indicated an association of TCA cycle suppression and inflammatory status [9].

The TCA cycle and OXPHOS are tightly controlled by the expression levels of associated genes; however, studies have not examined the relations of expression profiles of these genes and HCC prognosis. Thus, the objective of this study was to develop a HCC prognostic model based on the expression of genes associated with the TCA cycle and OXPHOS.

Results

Identification of DEGs associated with aerobic respiration

We obtained RNA-sequencing data of 50 paired HCC tissues and adjacent normal tissues from TCGA. RNA-sequencing data of 203 HCC tissues and 201 matched adjacent normal tissues from ICGC were also included. Using a threshold |log2 fold-change| > 1.0 and false discovery rate (FDR) adjusted to P < 0.05, 3 up-regulated genes (ATP6V1C1, SLC16A3, ME1) and 12 down-regulated genes (GSTZ1, ETFDH, ACAA1, ACADSB, ACAT1, BDH2, PCK1, ACAA2, ALDH6A1, ECHS1, OGDHL, PCK2) were identified (Figure 1, Figure 2).

Volcano plot showing differentially expressed genes in aerobic respiration between hepatocellular carcinoma and normal tissues. Red dots represent significantly up-regulated genes, blue dots represent significantly down-regulated genes, and gray dots represent no differences between genes.

Figure 1. Volcano plot showing differentially expressed genes in aerobic respiration between hepatocellular carcinoma and normal tissues. Red dots represent significantly up-regulated genes, blue dots represent significantly down-regulated genes, and gray dots represent no differences between genes.

Heat map showing differentially expressed genes in aerobic respiration between 50 matched hepatocellular carcinoma and normal tissue pairs in the TCGA cohort.

Figure 2. Heat map showing differentially expressed genes in aerobic respiration between 50 matched hepatocellular carcinoma and normal tissue pairs in the TCGA cohort.

Functional enrichment analysis of the DEGs

To explore the biological function of the 15 genes, we conducted functional enrichment analysis via GO annotation and KEGG pathway analyses. Top 30 genes identified by GO analysis are shown in Figure 3, and the results demonstrated that the most frequent biological process category was organic acid catabolic processes. The top 30 enrichment pathways are summarized in Figure 4. The results indicated that the most enriched pathway was fatty acid degradation. Protein-protein network analysis showed that ACAA1, ACAA2, ECHS1, ACAT1 and ALDH6A1 were the core genes (Figure 5).

GO analysis showing the biological processes and molecular functions associated with the differentially expressed genes in aerobic respiration.

Figure 3. GO analysis showing the biological processes and molecular functions associated with the differentially expressed genes in aerobic respiration.

KEGG analysis showing the signaling pathway involved in the differentially expressed genes in aerobic respiration.

Figure 4. KEGG analysis showing the signaling pathway involved in the differentially expressed genes in aerobic respiration.

Protein-protein interaction (PPI) network of 15 differentially expressed genes in aerobic respiration.

Figure 5. Protein-protein interaction (PPI) network of 15 differentially expressed genes in aerobic respiration.

Development of a prognostic model and model validation

In the 357 patients from TCGA cohort, univariate Cox regression analysis indicated that 15 genes were significantly related to DSS. LASSO Cox regression analysis indicated that 12 of the genes had maximum prognostic value (Figure 6A, 6B). These genes were then used to develop a risk score model to predict the prognostic of HCC patients.

(A) LASSO coefficient profiles of the 12 genes of high prognostic value. (B) The optimal values of the penalty parameter were determined by 10-fold cross-validation.

Figure 6. (A) LASSO coefficient profiles of the 12 genes of high prognostic value. (B) The optimal values of the penalty parameter were determined by 10-fold cross-validation.

Riskscore = ACAA1 × (–0.8576835) + ACAT1 × 0.1040123 + ATP6V1C1 × (–0.3383703) + BDH2 × (–0.1431161) + ECHS1 × 0.1992799 + ETFDH × 0.2083975 + GSTZ1 × 0.1662454 + ME1 × (–0.8860555) + OGDHL × 0.4802452 + PCK1 × 0.2184216 + PCK2 × 0.2346868 + SLC16A3 × (–0.2621411).

Using the risk score, patients from TCGA were categorized as high-risk and low-risk. The area under the receiver operating characteristic (ROC) curve (AUC) of the prognostic model for DSS was 0.72 at 1 year, 0.762 at 3 years, and 0.745 at 5 years (Figure 7). Kaplan-Meier analysis showed a significant survival difference between the 2 groups (P < 0.001) (Figure 8). The riskscore was significantly higher (P < 0.05) in patients with higher grade tumors (Figure 9A), more advanced T stage (Figure 9B), and more advanced American Joint Committee on Cancer (AJCC) stage (Figure 9C), but was not associated with HBV status (Supplementary Figure 1).

Survival-dependent receiver operating characteristic (ROC) curves showing the prognostic value of the model.

Figure 7. Survival-dependent receiver operating characteristic (ROC) curves showing the prognostic value of the model.

Hepatocellular carcinoma patients from the TCGA cohort in the low-risk group had better disease-specific survival (DSS) (P

Figure 8. Hepatocellular carcinoma patients from the TCGA cohort in the low-risk group had better disease-specific survival (DSS) (P < 0.001).

Significant differences of risk score were found between: (A) different tumor grades; (B) different T stages; and (C) different AJCC stages.

Figure 9. Significant differences of risk score were found between: (A) different tumor grades; (B) different T stages; and (C) different AJCC stages.

After controlling for the impact of clinical factors such as pathological stage and histological grade via multivariate Cox regression analysis, the prognosis of the low-risk group remained better than that of the high-risk group (hazard ratio [HR] = 0.24; 95% confidence interval [CI]: 0.13–0.44; P < 0.001) (Figure 10).

Low-risk of HCC patients from the TCGA cohort is an independent prognostic factor for disease-specific survival (DSS).

Figure 10. Low-risk of HCC patients from the TCGA cohort is an independent prognostic factor for disease-specific survival (DSS).

We also constructed a nomogram that integrated risk score and clinical data including pathological stage and histological grade to predict the prognosis of HCC patients in the TCGA cohort more precisely, with a C-index = 0.768(95% CI: 0.712–0.824). Based on the score for each variable on the point scale of this nomogram, the probability of survival at 1, 3, and 5 years can be determined by calculating the total score (Figure 11).

Nomogram to predict the probability of disease-specific survival (DSS) at 1, 3, and 5 years for HCC patients in the TCGA cohort.

Figure 11. Nomogram to predict the probability of disease-specific survival (DSS) at 1, 3, and 5 years for HCC patients in the TCGA cohort.

To examine whether the riskscore model was accurate in different cohorts, the same cut-off values of the 12 genes were used to calculate the risk scores of the 235 HCC patients in the ICGC cohort. The patients were divided into high-risk and low-risk groups. Consistent with the findings of the TCGA cohort, patients in the low-risk group had a better overall survival (OS) than patients in the high-risk group (P = 0.0048) (Figure 12). The difference remained significant after controlling for clinical factors (HR = 0.44, 95% CI: 0.23–0.86, P = 0.016) (Figure 13).

Validation of the prognostic model in the ICGC cohort by univariate Cox analysis.

Figure 12. Validation of the prognostic model in the ICGC cohort by univariate Cox analysis.

Validation of the prognostic model in the ICGC cohort by multivariate Cox analysis.

Figure 13. Validation of the prognostic model in the ICGC cohort by multivariate Cox analysis.

Taken together, the results indicate that the risk score can be used as an independent prognostic factor for patients with HCC.

Discussion

In recent years, more research has been devoted to examining the relations between energy metabolism and human malignancy. Numerous studies have focused on the impact of glycolysis on HCC progression, while the clinical significance of aerobic respiration has not been deeply explored. In this study, we compared gene expression from matched tumor and normal tissues in TCGA and ICGC cohorts, and identified 15 DEGs. After setting optimal cut-off points, univariate Cox regression analysis showed that all 15 genes were significantly associated with DSS in HCC patients from the TCGA cohort. Finally, based on the LASSO estimation [10], we chose 12 genes to construct a prognostic model.

Many studies have attempted to develop prognostic models for HCC patients. Chen et al. developed a 4-gene signature model to predict the prognosis of HCC patients, with a model formula that is more concise than ours [11]. However, the AUC or C-index was not provided in that study, and thus it is difficult to evaluate the quality of the model. A prognostic model for patients with HCC undergoing radiofrequency ablation (RFA) was recently reported [12]. The model incorporated tumor size and number, alpha fetoprotein (AFP) level (ng/ml), protein induced by vitamin K absence (PIVKA) level (mAU/ml), lymphocyte count (/ul), serum albumin level (g/dl), and ascites, and a nomogram was developed to predict patient survival. However, the C-index of the model (0.759) was a little lower than that of our model (0.768). As one of main metabolic features of tumor cells, genes involved in glycolysis have been screened to establish a six gene signature for HCC patients [13]. The diagnostic efficacy of the prognostic model was good (AUC = 0.765), but the prognostic value was not evaluated. We also tried to develop a prognostic model based on both glycolysis and aerobic respiration, but the result seemed to be unsatisfactory (Supplementary Figure 2), which might be mainly explained by the high heterogeneity of aerobic respiration and glycolysis pathways.

Most of the genes included in our prognostic model have been reported to be associated with survival in many different malignancies. ACAA1 may be a marker for non-small cell lung cancer (NSCLC) diagnosis and prognosis, and may provide new insights for NSCLC treatment [14]. Up-regulation of ACAT1 expression is involved in the progression of colorectal cancer (CRC) [15, 16]. Breast cancer growth and bone metastasis are prevented by silencing of ATP6V1C1, suggesting a potential therapy target [17]. BDH2 is a tumor suppressor gene in HCC that regulates apoptosis and autophagy [18]. Up-regulation of ECHS1 may promote cell proliferation in HCC [19], which is different than the results or our study. The difference of results may be explained by the sample divergence: ETDFH expression level has been found to be associated with AFP level, and it may also be an independent prognostic factor in HCC patients [20]. GSTZ1 acts as a tumor suppressor gene in HCC, and GSTZ1 deficiency might accelerate HCC progression [21]. ME1 has been reported to be associated with poor survival in gastric cancer patients, and it may be an oncogene [22]. Down-regulation of OGDHL expression is correlated with poor prognosis in patients with pancreatic cancer [23]. PCK1 suppressed HCC through promoting TCA cataplerosis, oxidative stress, and apoptosis in HCC cells [24]. Low expression of PCK2 has been observed in osteosarcoma patients with metastasis and recurrence, and is associated with poorer survival [25]. SLC16A3 has been shown to be a possible biomarker for prognosis of pancreatic cancer [26].

In this study, HCC patients identified in TCGA were dichotomized into a high-risk group and low-risk group by the risk score. Univariate and multivariate Cox regression analysis both indicated that the low-risk group had better DSS. We also validated the prognostic model in an ICGC cohort and obtained similar results. Based on the development of the prognostic model, low-risk was defined as a risk score < –0.6943, and high-risk was defined as a risk score > –0.6943. Overall, patients in the low-risk group had a better prognosis than patients in high-risk group, and a relatively precise survival probability could be calculated for each patient using the nomogram developed in the study. As such, the prognostic model may be used for the risk stratification in HCC patients, and potentially guide individualized patient treatments.

Several factors may have potentially influenced the outcome of this research. First, batch effects are important source of error in research based on RNA-sequencing data. In this study, we convert the format of RNA-sequencing data from fragments per kilobase of exon model per million mapped fragments (FPKM) to transcripts per kilobase of exon model per million mapped reads (TPM), which removed the impact of different sequencing depths. Second, and importantly, we identified DEGs from paired tumor and normal tissues, because we thought tumor and normal tissues from the same patient were comparable. Third, we used DSS as the prognostic indicator, which removed the impact of death from other factors to some extent. In the ICGC cohort, only 1 survival indicator was provided, so we empirically excluded patients with survival < 30 days to remove patients who died from other causes.

Some limitations of this study must be pointed out. We did not explore the possible mechanisms by which the genes identified serve a prognostic role. Besides, there is the potential of different sequencing data from different sequencing platforms; however, the prognostic model exhibited ideal performance based on Illumina platform. While the results of this study are promising, they need to be validated in other studies with larger numbers of patients and different sequencing platform. Finally, our study was a retrospective research, so it was difficult to prevent some inherent bias.

In conclusion, this study constructed and validated a prognostic model of HCC through a comprehensive analysis of genes involved in aerobic respiration. Based on risk score, HCC patients could be stratified into high-risk and low-risk groups, which may assist in developing individualized treatments for patients with HCC.

Methods

Data acquisition and processing

RNA-sequencing and clinical data were downloaded from the LIHC project of The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/), and the LIRI-JP project of International Cancer Genome Consortium (ICGC) data portal (https://dcc.icgc.org/releases/current/Projects/LIRI-JP).

A total of 285 genes were obtained for the citric acid cycle (TCA) and for respiratory electron transport (R-HSA-1428517) in Reactome (https://reactome.org/PathwayBrowser/#/R-HSA-1428517&PATH=R-HSA-1430728&DTAB=MT). Oxidative phosphorylation hallmark gene sets and KEGG_CITRATE_CYCLE_TCA_CYCLE of KEGG gene sets were obtained from the molecular signatures database (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#C5).

For calculations, the format of RNA-sequencing data was converted from fragments per kilobase of exon model per million mapped fragments (FPKM) to transcripts per kilobase of exon model per million mapped reads (TPM).

Differential and functional analysis

The “edgR” package in R software was used to analyze the raw RNA-sequencing data counts from TCGA and the ICGC with the criteria of: |log2 fold-change| > 1.0, and false discovery rate (FDR) adjusted to P < 0.05. A total of 15 overlapping differentially expressed genes (DEGs) associated with aerobic respiration were identified. Volcano plots and heat maps of the DEGs were generated to visualize the results using the “ggpubr” package and the “pheatmap” package. To understand the biological processes of the DEGs, functional enrichment analysis of the DEGs was performed with the “clusterProfile” package [27]. Gene Ontology (GO) and Kyoto Gene and Genomic Encyclopedia (KEGG) analyses were conducted to assess relevant functional categories based on the threshold of P < 0.05.

Construction of the prognostic model and model validation

To eliminate patients who might die of non-cancer-related diseases, we included 357 patients from TCGA who were followed-up for more than 1 day with clear disease-specific survival (DSS) data, and 235 patients from the ICGC who were followed-up for more than 30 days with clear survival status data. The optimum cut-off point of DEG expression was determined using the “surv_cutpoint” function of the “survminer” package in R. Ideal DEGs with a significance level of P < 0.05 were selected through univariate Cox regression analysis using the “survival” package in R. The LASSO function of the “glmnet” package in R was used to select the most useful prognostic genes among the DEGs. The most useful prognosis-related DEGs were selected based on LASSO estimation, and the optimal value for penalization coefficient lambda was chosen after running cross-validation likelihood (cvl) 1,000 times [10]. Subsequently, the Cox coefficients and expression values of the DEGs were extracted to calculate the risk score. The expression value of a gene, which was considered highly expressed, was defined as 0, and the expression value of a gene which was considered to have low expression was defined as 1.

A formula for the prognostic model was established to predict patient survival based on the formula:

Risk score = ∑Cox coefficient of gene Xi × expression value of gene Xi. The sensitivity and specificity of survival prediction based on the risk score was calculated, and receiver operating characteristic (ROC) curve analysis was performed using the “survivalROC” package in R. We also used multivariate Cox regression analysis to control for the impact of potential confounding factors, and Forest plots and a nomogram were constructed using the “rms” package and “ggforest” functions in the “survminer” package of R. The C-index of the prognostic model was calculated using the “rcorrcens” function. The prognostic model was validated using the ICGC cohort and the same cut-off points.

Overall, we developed a flow chart to help readers understand the steps of our manuscript (Supplementary Figure 3).

Supplementary Materials

Supplementary Figures

Acknowledgments

Conceptualization of the study: Yi Ma; Investigations and statistical analysis: Jiawei Rao, Xukun Wu, Xiaozhuan Zhou; Writing the original manuscript: Jiawei Rao, Xukun Wu; Writing review and editing: Jiawei Rao, Ronghai Deng; Funding acquisition: Yi Ma.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Funding

Funding for this research was obtained from the following sources: National Natural Science Foundation of China, China (81873591, and 81670591); Guangdong Natural Science Foundation, China (2016A030311028); Science and Technology Planning Project of Guangdong Province, China (2018A050506030); Science and Technology Program of Guangzhou, China (201704020073); Guangdong Provincial Key Laboratory Construction Projection on Organ Donation and Transplant Immunology (2013A061401007 and 2017B030314018); Guangdong Provincial International Cooperation Base of Science and Technology (Organ Transplantation) (2015B050501002).

References

  • 1. Warburg O. On the origin of cancer cells. Science. 1956; 123:309–14. https://doi.org/10.1126/science.123.3191.309 [PubMed]
  • 2. Pfeiffer T, Schuster S, Bonhoeffer S. Cooperation and competition in the evolution of ATP-producing pathways. Science. 2001; 292:504–07. https://doi.org/10.1126/science.1058079 [PubMed]
  • 3. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, Parkin DM, Forman D, Bray F. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer. 2015; 136:E359–86. https://doi.org/10.1002/ijc.29210 [PubMed]
  • 4. De Lorenzo S, Tovoli F, Barbera MA, Garuti F, Palloni A, Frega G, Garajovà I, Rizzo A, Trevisani F, Brandi G. Metronomic capecitabine vs. best supportive care in Child-Pugh B hepatocellular carcinoma: a proof of concept. Sci Rep. 2018; 8:9997. https://doi.org/10.1038/s41598-018-28337-6 [PubMed]
  • 5. Lee YK, Lim JJ, Jeoun UW, Min S, Lee EB, Kwon SM, Lee C, Yoon G. Lactate-mediated mitoribosomal defects impair mitochondrial oxidative phosphorylation and promote hepatoma cell invasiveness. J Biol Chem. 2017; 292:20208–17. https://doi.org/10.1074/jbc.M117.809012 [PubMed]
  • 6. Kim JH, Kim HY, Lee YK, Yoon YS, Xu WG, Yoon JK, Choi SE, Ko YG, Kim MJ, Lee SJ, Wang HJ, Yoon G. Involvement of mitophagy in oncogenic K-Ras-induced transformation: overcoming a cellular energy deficit from glucose deficiency. Autophagy. 2011; 7:1187–98. https://doi.org/10.4161/auto.7.10.16643 [PubMed]
  • 7. Baracca A, Chiaradonna F, Sgarbi G, Solaini G, Alberghina L, Lenaz G. Mitochondrial Complex I decrease is responsible for bioenergetic dysfunction in K-ras transformed cells. Biochim Biophys Acta. 2010; 1797:314–23. https://doi.org/10.1016/j.bbabio.2009.11.006 [PubMed]
  • 8. Shen YC, Ou DL, Hsu C, Lin KL, Chang CY, Lin CY, Liu SH, Cheng AL. Activating oxidative phosphorylation by a pyruvate dehydrogenase kinase inhibitor overcomes sorafenib resistance of hepatocellular carcinoma. Br J Cancer. 2013; 108:72–81. https://doi.org/10.1038/bjc.2012.559 [PubMed]
  • 9. Huang Q, Tan Y, Yin P, Ye G, Gao P, Lu X, Wang H, Xu G. Metabolic characterization of hepatocellular carcinoma using nontargeted tissue metabolomics. Cancer Res. 2013; 73:4992–5002. https://doi.org/10.1158/0008-5472.CAN-13-0308 [PubMed]
  • 10. Tibshirani R. The lasso method for variable selection in the Cox model. Stat Med. 1997; 16:385–95. https://doi.org/10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3 [PubMed]
  • 11. Chen PF, Li QH, Zeng LR, Yang XY, Peng PL, He JH, Fan B. A 4-gene prognostic signature predicting survival in hepatocellular carcinoma. J Cell Biochem. 2019; 120:9117–24. https://doi.org/10.1002/jcb.28187 [PubMed]
  • 12. Kim CG, Lee HW, Choi HJ, Lee JI, Lee HW, Kim SU, Park JY, Kim DY, Ahn SH, Han KH, Kim HS, Kim KH, Choi SJ, et al. Development and validation of a prognostic model for patients with hepatocellular carcinoma undergoing radiofrequency ablation. Cancer Med. 2019; 8:5023–32. https://doi.org/10.1002/cam4.2417 [PubMed]
  • 13. Jiang L, Zhao L, Bi J, Guan Q, Qi A, Wei Q, He M, Wei M, Zhao L. Glycolysis gene expression profilings screen for prognostic risk signature of hepatocellular carcinoma. Aging (Albany NY). 2019; 11:10861–82. https://doi.org/10.18632/aging.102489 [PubMed]
  • 14. Zhang X, Yang H, Zhang J, Gao F, Dai L. HSD17B4, ACAA1, and PXMP4 in Peroxisome Pathway Are Down-Regulated and Have Clinical Significance in Non-small Cell Lung Cancer. Front Genet. 2020; 11:273. https://doi.org/10.3389/fgene.2020.00273 [PubMed]
  • 15. Chen X, Liang H, Song Q, Xu X, Cao D. Insulin promotes progression of colon cancer by upregulation of ACAT1. Lipids Health Dis. 2018; 17:122. https://doi.org/10.1186/s12944-018-0773-x [PubMed]
  • 16. Ye K, Wu Y, Sun Y, Lin J, Xu J. TLR4 siRNA inhibits proliferation and invasion in colorectal cancer cells by downregulating ACAT1 expression. Life Sci. 2016; 155:133–39. https://doi.org/10.1016/j.lfs.2016.05.012 [PubMed]
  • 17. Feng S, Zhu G, McConnell M, Deng L, Zhao Q, Wu M, Zhou Q, Wang J, Qi J, Li YP, Chen W. Silencing of atp6v1c1 prevents breast cancer growth and bone metastasis. Int J Biol Sci. 2013; 9:853–62. https://doi.org/10.7150/ijbs.6030 [PubMed]
  • 18. Liang H, Xiong Z, Li R, Hu K, Cao M, Yang J, Zhong Z, Jia C, Yao Z, Deng M. BDH2 is downregulated in hepatocellular carcinoma and acts as a tumor suppressor regulating cell apoptosis and autophagy. J Cancer. 2019; 10:3735–45. https://doi.org/10.7150/jca.32022 [PubMed]
  • 19. Lin BY, Xiao CX, Zhao WX, Xiao L, Chen X, Li P, Wang XM. Enoyl-coenzyme A hydratase short chain 1 silencing attenuates the proliferation of hepatocellular carcinoma by inhibiting epidermal growth factor signaling in vitro and in vivo. Mol Med Rep. 2015; 12:1421–28. https://doi.org/10.3892/mmr.2015.3453 [PubMed]
  • 20. Wu Y, Zhang X, Shen R, Huang J, Lu X, Zheng G, Chen X. Expression and significance of ETFDH in hepatocellular carcinoma. Pathol Res Pract. 2019; 215:152702. https://doi.org/10.1016/j.prp.2019.152702 [PubMed]
  • 21. Li J, Wang Q, Yang Y, Lei C, Yang F, Liang L, Chen C, Xia J, Wang K, Tang N. GSTZ1 deficiency promotes hepatocellular carcinoma proliferation via activation of the KEAP1/NRF2 pathway. J Exp Clin Cancer Res. 2019; 38:438. https://doi.org/10.1186/s13046-019-1459-6 [PubMed]
  • 22. Shi Y, Zhou S, Wang P, Guo Y, Xie B, Ding S. Malic enzyme 1 (ME1) is a potential oncogene in gastric cancer cells and is associated with poor survival of gastric cancer patients. Onco Targets Ther. 2019; 12:5589–99. https://doi.org/10.2147/OTT.S203228 [PubMed]
  • 23. Liu Y, Meng F, Wang J, Liu M, Yang G, Song R, Zheng T, Liang Y, Zhang S, Yin D, Wang J, Yang H, Pan S, et al. A Novel Oxoglutarate Dehydrogenase-Like Mediated miR-214/TWIST1 Negative Feedback Loop Inhibits Pancreatic Cancer Growth and Metastasis. Clin Cancer Res. 2019; 25:5407–21. https://doi.org/10.1158/1078-0432.CCR-18-4113 [PubMed]
  • 24. Liu MX, Jin L, Sun SJ, Liu P, Feng X, Cheng ZL, Liu WR, Guan KL, Shi YH, Yuan HX, Xiong Y. Metabolic reprogramming by PCK1 promotes TCA cataplerosis, oxidative stress and apoptosis in liver cancer cells and suppresses hepatocellular carcinoma. Oncogene. 2018; 37:1637–53. https://doi.org/10.1038/s41388-017-0070-6 [PubMed]
  • 25. Zhang Y, Zhao HJ, Xia XR, Diao FY, Ma X, Wang J, Gao L, Liu J, Gao C, Cui YG, Liu JY. Hypoxia-induced and HIF1α-VEGF-mediated tight junction dysfunction in choriocarcinoma cells: Implications for preeclampsia. Clin Chim Acta. 2019; 489:203–11. https://doi.org/10.1016/j.cca.2017.12.010 [PubMed]
  • 26. Yu S, Wu Y, Li C, Qu Z, Lou G, Guo X, Ji J, Li N, Guo M, Zhang M, Lei L, Tai S. Comprehensive analysis of the SLC16A gene family in pancreatic cancer via integrated bioinformatics. Sci Rep. 2020; 10:7315. https://doi.org/10.1038/s41598-020-64356-y [PubMed]
  • 27. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]