# Identification of Hypoxia-Related Differentially Expressed Genes and Construction of the Clinical Prognostic Predictor in Hepatocellular Carcinoma by Bioinformatic Analysis

Authors: Hao Guo; Jing Zhou; Yanjun Zhang; Zhi Wang; Likun Liu; Weiwei Zhang
Journal: BioMed Research International
 (2021)

---

## Abstract

Background. Hypoxia closely relates to malignant progression and appears to be prognostic for outcome in hepatocellular carcinoma (HCC). Our research is aimed at mining the hypoxic-related genes (HRGs) and constructing a prognostic predictor (PP) model on clinical prognosis in HCC patients. Methods. RNA-sequencing data about HRGs and clinical data of patients with HCC were obtained from The Cancer Genome Atlas (TCGA) database portal. Differentially expressed HRGs between HCC and para-carcinoma tissue samples were obtained by applying the Wilcox analysis in R statistical software. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes were used for gene functional enrichment analyses. Then, the patients who were asked to follow up for at least one month were enrolled in the following study. Cox proportional risk regression model was applied to obtain key HRGs which related to overall survival (OS) in HCC. PP was constructed and defined, and the accuracy of PP was validated by constructing the signature in a training set and validation set. Connectivity map (CMap) was used to find potential drugs, and gene set cancer analysis (GSCA) was also performed to explore the underlying molecular mechanisms. Results. Thirty-seven differentially expressed HRGs were obtained. It contained 28 upregulated and 9 downregulated genes. After the univariate Cox regression model analysis, we obtained 27 prognosis-related HRGs. Of these, 25 genes were risk factors for cancer, and 2 genes were protective factors. The PP was composed by 12 key genes (HDLBP, SAP30, PFKP, DPYSL4, SLC2A1, HMOX1, PGK1, ERO1A, LDHA, ENO2, SLC6A6, and TPI1). GSCA results showed the overall activity of these 12 key genes in 10 cancer-related pathways. Besides, CMap identified deferoxamine, crotamiton, talampicillin, and lycorine might have effects with HCC. Conclusions. This study firstly reported 12 prognostic HRGs and constructed the model of the PP. This comprehensive research of multiple databases helps us gain insight into the biological properties of HCC and provides deferoxamine, crotamiton, talampicillin, and lycorine as potential drugs to fight against HCC.


---

## Body

## 1. Introduction

Hypoxia is a hallmark of the tumor microenvironment and contributes to tumor malignant phenotypes. It presents in the majority of tumors, and the distribution of hypoxic regions in tumors is heterogeneous [1]. Hypoxia has been considered as a vital sign of poor prognosis in almost all solid tumors, including the hepatocellular carcinoma (HCC) [2]. Hypoxia can lead to changes in the proteome and/or genome. The mechanism may be related to cells’ continued access to nutrients, escape from unsatisfactory microenvironments, and promote unrestricted growth to accelerate tumor progression [3]. Tumor cells in hypoxic regions are closely associated with treatment resistance, local recurrence, and distant metastasis [4]. All these can be regarded as important factors for poor prognosis of patients.Liver cancer is the fourth leading cause of cancer death in 2018. It kills about 782,000 people annually [5, 6]. It has been predicted that liver cancer may become the sixth most commonly diagnosed cancer in the future. HCC comprises 75%–85% of primary cancer of the liver. Unfortunately, the exact mechanism and pathways leading to HCC development are still unclear explanation [7]. In particular, the absence of indicators is related to the prognosis of HCC. The 5-year survival rate is still very low [8, 9]. Therefore, finding the key molecular biomarkers focused on hypoxia which is related to the occurrence and development of HCC patients has certain reference value for reliable estimation of the deterioration of HCC and may be an effective measure to against HCC.In recent years, progress using high-throughput sequencing technologies has greatly expanded the research of cancer genome. Here, we compared the expression profile of hypoxic-related genes (HRGs) between HCC tissue samples and para-carcinoma tissue samples which were collected from The Cancer Genome Atlas (TCGA) database [10]. Furthermore, in order to clarify the correlation between the HRGs and clinical outcomes in HCC patients, the prognostic predictor (PP) was developed and validated as an independent indicator to assess the overall survival (OS) of patients. We also tried to find potential molecular drugs to fight against HCC. These findings may provide new research targets and therapeutic strategies for HCC patients and provide a reliable theoretical basis for judging treatment outcomes and assessing prognosis.


## 2. Materials and Methods

### 2.1. Data Acquisition

The Molecular Signatures Database (MSigDB,http://software.broadinstitute.org/gsea/msigdb/index.jsp) is a collection of annotated gene sets for GSEA software. A variety of HRGs were extracted from the Hallmark hypoxic gene sets in the MSigDB [11, 12], which contains a list of 200 HRGs involved in hypoxic regulatory pathways reported in the literature. RNA-sequencing data about HRGs and clinical data of patients with HCC were obtained from the TCGA data portal. The dataset consists of 374 HCC tissue samples and 50 para-carcinoma tissue samples. We performed on these samples to obtain the differentially expressed HRGs. Then, in order to build a model related to the prognosis of HCC, a total of 349 patients who were asked to follow up for at least one month were enrolled in our study. The follow-up periods of included patients were for a range of 30 days to 3675 days.


### 2.2. Differentially Expressed HRG Enrichment Analysis

We screened for differences in HRG expression between HCC tissue samples and para-carcinoma tissue samples by applying the Wilcox analysis in R statistical software. The significantly differentially expressed HRGs were defined as thep value is less than 0.05; the gene expression changes at least 1.5-fold. To explain the main biological properties of HRGs, the tools of Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used for gene functional enrichment analyses. GO and KEGG terms with a p value and q value both smaller than 0.05 were considered significant categories. To take full advantage of the well-known GO and KEGG themes and highlight the most relevant GO terms associated with a given gene list, these results of gene enrichment were visualized by R language clustering package [13].


### 2.3. Construction of the Individualized Prognostic Predictor Model Based on HRGs

We extracted the expression profile of HRGs downloaded from TCGA in the FPKM format. Univariate Cox regression analysis was performed to screen HRGs whose expression profile was significantly correlated with OS in HCC patients from 200 hypoxic-related genes (p<0.01). Multivariate Cox regression analysis was performed on these survival-related genes to remove genes that might not be independent indicators of prognosis. We got the optimal solution of the model by narrowing down the relevant genes step by step. Finally, the PP was constructed and significantly divided patients of HCC into high- and low-risk groups according to OS. Kaplan-Meier (K-M) method was used to plot the survival curve. The survival rates of the two groups were tested by log-rank test.To further verify the hypoxia-related gene model and whether the PP could be used as an independent predictor of prognosis for HCC patients, we divided the TCGA data set into a training set and a validation set, and performed survival analysis on the hypoxia-related signature in the two data sets to illustrate the robustness of the PP. Besides, we introduced several following indicators closely related to the clinic into the model as covariables for multivariate Cox regression analysis: PP and age were coded as continuous variables, gender (male or female), tumor grade (high or low), AJCC TNM stage (I=1, II=2, III=3, and IV=4), tumor stage [1–4], lymph node metastasis (positive or negative), and distant metastasis (positive or negative).


### 2.4. Identification of Potential Therapeutic Agents

Connectivity map (CMap) was used to identify potential therapeutic agents for HCC. It is a program that predicts potential drugs and produces biological states by encoding specific gene expression characteristics [14]. 28 upregulated and 9 downregulated HRGs were used to check the CMap database. Finally, the enrichment score representing similarity was calculated, ranging from -1 to 1. The positive connectivity score indicates that drugs can induce the biological phenomena queried in human cell lines. In contrast, a negative connectivity score indicates that the drug reverses the requested biological characteristics and has potential therapeutic value. With p<0.05 as the cut-off value, the connectivity scores of different instances were filtered out. We carried out studies on these potentially effective drugs in PubChem.


### 2.5. Gene Set Cancer Analysis (GSCA)

To further explore the underlying molecular mechanisms involved in these modeling genes, GSCA (http://bioinfo.life.hust.edu.cn/web/GSCALite/) [15] was used to obtain the activity of 10 cancer-related pathways. Genes were divided into high and low group by median expression. Student’s T test was used to define the difference of pathway activity score (PAS) between groups. p value was adjusted by FDR (FDR≤0.05 is considered as significant).


### 2.6. Statistical Analysis

All statistical analyses were conducted using SPSS 22.0 (Chicago, IL, USA) and R 3.6.0 (https://www.r-project.org/); Adobe Illustrator CS 5 (San Jose, California, USA) was performed to draw plots. Univariate Cox regression analysis was used for illustrating the correlation between genes expression profiles and OS. The multivariate Cox proportional hazard regression model was used for performing multivariate analysis and construct prognostic models. Receiver-operating characteristic (ROC) curve and the corresponding area under the ROC curve (AUC) were applied for verifying that PP was accurate. All tests were performed bilaterally, p<0.05 indicates a difference with statistical significance.


## 2.1. Data Acquisition

The Molecular Signatures Database (MSigDB,http://software.broadinstitute.org/gsea/msigdb/index.jsp) is a collection of annotated gene sets for GSEA software. A variety of HRGs were extracted from the Hallmark hypoxic gene sets in the MSigDB [11, 12], which contains a list of 200 HRGs involved in hypoxic regulatory pathways reported in the literature. RNA-sequencing data about HRGs and clinical data of patients with HCC were obtained from the TCGA data portal. The dataset consists of 374 HCC tissue samples and 50 para-carcinoma tissue samples. We performed on these samples to obtain the differentially expressed HRGs. Then, in order to build a model related to the prognosis of HCC, a total of 349 patients who were asked to follow up for at least one month were enrolled in our study. The follow-up periods of included patients were for a range of 30 days to 3675 days.


## 2.2. Differentially Expressed HRG Enrichment Analysis

We screened for differences in HRG expression between HCC tissue samples and para-carcinoma tissue samples by applying the Wilcox analysis in R statistical software. The significantly differentially expressed HRGs were defined as thep value is less than 0.05; the gene expression changes at least 1.5-fold. To explain the main biological properties of HRGs, the tools of Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used for gene functional enrichment analyses. GO and KEGG terms with a p value and q value both smaller than 0.05 were considered significant categories. To take full advantage of the well-known GO and KEGG themes and highlight the most relevant GO terms associated with a given gene list, these results of gene enrichment were visualized by R language clustering package [13].


## 2.3. Construction of the Individualized Prognostic Predictor Model Based on HRGs

We extracted the expression profile of HRGs downloaded from TCGA in the FPKM format. Univariate Cox regression analysis was performed to screen HRGs whose expression profile was significantly correlated with OS in HCC patients from 200 hypoxic-related genes (p<0.01). Multivariate Cox regression analysis was performed on these survival-related genes to remove genes that might not be independent indicators of prognosis. We got the optimal solution of the model by narrowing down the relevant genes step by step. Finally, the PP was constructed and significantly divided patients of HCC into high- and low-risk groups according to OS. Kaplan-Meier (K-M) method was used to plot the survival curve. The survival rates of the two groups were tested by log-rank test.To further verify the hypoxia-related gene model and whether the PP could be used as an independent predictor of prognosis for HCC patients, we divided the TCGA data set into a training set and a validation set, and performed survival analysis on the hypoxia-related signature in the two data sets to illustrate the robustness of the PP. Besides, we introduced several following indicators closely related to the clinic into the model as covariables for multivariate Cox regression analysis: PP and age were coded as continuous variables, gender (male or female), tumor grade (high or low), AJCC TNM stage (I=1, II=2, III=3, and IV=4), tumor stage [1–4], lymph node metastasis (positive or negative), and distant metastasis (positive or negative).


## 2.4. Identification of Potential Therapeutic Agents

Connectivity map (CMap) was used to identify potential therapeutic agents for HCC. It is a program that predicts potential drugs and produces biological states by encoding specific gene expression characteristics [14]. 28 upregulated and 9 downregulated HRGs were used to check the CMap database. Finally, the enrichment score representing similarity was calculated, ranging from -1 to 1. The positive connectivity score indicates that drugs can induce the biological phenomena queried in human cell lines. In contrast, a negative connectivity score indicates that the drug reverses the requested biological characteristics and has potential therapeutic value. With p<0.05 as the cut-off value, the connectivity scores of different instances were filtered out. We carried out studies on these potentially effective drugs in PubChem.


## 2.5. Gene Set Cancer Analysis (GSCA)

To further explore the underlying molecular mechanisms involved in these modeling genes, GSCA (http://bioinfo.life.hust.edu.cn/web/GSCALite/) [15] was used to obtain the activity of 10 cancer-related pathways. Genes were divided into high and low group by median expression. Student’s T test was used to define the difference of pathway activity score (PAS) between groups. p value was adjusted by FDR (FDR≤0.05 is considered as significant).


## 2.6. Statistical Analysis

All statistical analyses were conducted using SPSS 22.0 (Chicago, IL, USA) and R 3.6.0 (https://www.r-project.org/); Adobe Illustrator CS 5 (San Jose, California, USA) was performed to draw plots. Univariate Cox regression analysis was used for illustrating the correlation between genes expression profiles and OS. The multivariate Cox proportional hazard regression model was used for performing multivariate analysis and construct prognostic models. Receiver-operating characteristic (ROC) curve and the corresponding area under the ROC curve (AUC) were applied for verifying that PP was accurate. All tests were performed bilaterally, p<0.05 indicates a difference with statistical significance.


## 3. Results

### 3.1. Differentially Expressed HRGs

A total of 374 HCC tissue samples and 50 para-carcinoma tissue specimens were obtained for RNA-seq and clinical data from TCGA. A list of 200 HRGs involved in hypoxic regulatory pathways reported in the literature was extracted from the Hallmark hypoxic gene sets in the MSigDB database. According to the limitation of aFDR<0.05 and ∣log2FoldChange∣>1.5, our study screened 37 genes with significant differential expression, of which 28 were upregulated and 9 were downregulated (Figures 1(a) and 1(b)). Then, we visualized the expression patterns of these HRGs between HCC and normal tissues by making boxplot (Figure 1(c)). The boxplot showed expression patterns that 28 upregulated genes (ALDOA, ANKZF1, ANXA2, B4GALNT2, BCAN, CA12, COL5A1, DDIT3, DPYSL4, DTNA, EFNA3, GPC3, HOXB9, INHA, KDELR3, KIF5A, LOX, NDRG1, P4HA2, PDGFB, PFKP, PGF, PPFIA4, RRAGD, SLC2A1, SLC2A5, STC1, and STC2) and 9 downregulated genes (DCN, FBP1, FOS, MT1E, MT2A, PCK1, PLAC8, SERPINE1, and ZFP36).Figure 1
Differentially expressed HRGs between HCC and para-carcinoma tissues. (a) The volcano plot for the 200 HRGs from the TCGA data portal. 37 significantly different genes were labeled in the figure. Red indicates high expression and green low expression. (b) Hierarchical clustering of differentially expressed HRG expression levels. (c) The expression patterns of box plots displayed 37 differentially expressed HRGs in HCC types and paired nontumor samples. Each red dot represents a distinct tumor sample and green represents a nontumor sample.
(a)(b)(c)


### 3.2. Functional Annotation and Enrichment Analysis of Differentially Expressed HRGs

To further understand the biologically relevant information of these 37 differential HRGs, we performed functional enrichment and enrichment analysis. The GO term function (Figure2) and KEGG pathway enrichment (Figure 3) of these genes were reviewed, respectively. We found the results that the GO terms for biological processes were carbohydrate biosynthetic process, NADH regeneration, canonical glycolysis, glucose catabolic process to pyruvate, and so on. The detailed results of all the above gene enrichment are shown in Figure 2.Figure 2
The circle plot of the top 10 GO terms function enrichment of differentially expressed HRGs.Figure 3
The circle plot of the top 10 enriched KEGG pathways of differentially expressed HRGs.Through the enrichment analysis function in the KEGG pathway, we found many important pathways associated with these genes such as HIF-1 signaling pathway, glycolysis/gluconeogenesis, carbon metabolism, MAPK signaling pathway, and PI3K-Akt signaling pathway (Figure3). As shown in the circle plot, the change from the Z-score indicated mostly related signaling pathways were more inclined to be increased.


### 3.3. Identification of Prognostic HRGs

In order to build a model related to the prognosis of HCC, we selected patients with follow-up period of more than 1 month from this data set. Finally, we obtained 349 patients for the analysis and establishment of the prognostic predictor (PP). Prognostic HRGs were computed with univariate Cox regression model by SPSS 22.0. The relationship between the expression profiles of HRGs and OS was evaluated using TCGA data and the Hallmark hypoxic gene sets in the MSigDB database, resulting in 27 prognosis-related HRGs (HR>1 or HR<1, p<0.01) in Table 1.Table 1
The associated HRGs with OS for univariate Cox regression model in HCC patients.
IDHRHR.95% CIp valueHDLBP1.761.19-2.60<0.01SAP301.921.47-2.51<0.01NDRG11.341.17-1.54<0.01ADORA2B1.561.16-2.09<0.01PFKP1.261.13-1.42<0.01DPYSL41.551.12-2.14<0.01PFKFB31.211.07-1.37<0.01SLC2A11.621.38-1.90<0.01VEGFA1.431.14-1.80<0.01HMOX11.261.09-1.45<0.01SDC31.371.11-1.70<0.01ATP7A2.011.28-3.16<0.01PGK11.451.17-1.79<0.01ERO1A1.481.17-1.87<0.01XPNPEP12.031.32-3.12<0.01LDHA2.091.60-2.73<0.01ENO21.291.10-1.52<0.01GRHPR0.730.60-0.91<0.01UGP20.740.60-0.92<0.01SLC6A61.291.08-1.53<0.01MAP3K11.461.11-1.92<0.01BNIP3L1.431.12-1.82<0.01TPI11.611.24-2.08<0.01RRAGD1.321.07-1.62<0.01TES1.281.07-1.54<0.01SERPINE11.161.05-1.28<0.01JMJD61.771.32-2.37<0.01HR: hazard ratio; CI: confidence interval.By mapping the forest of HRGs and OS (Figure4), we can intuitively understand the role of these genes in cancer prognosis. Of these, 25 genes (HDLBP, SAP30, NDRG1, ADORA2B, PFKP, DPYSL4, PFKFB3, SLC2A1, VEGFA, HMOX1, SDC3, ATP7A, PGK1, ERO1A, XPNPEP1, LDHA, ENO2, SLC6A6, MAP3K1, BNIP3L, TPI1, RRAGD, TES, SERPINE1, and JMJD6) were risk factors for cancer, and 2 genes (GRHPR and UGP2) were protective factors.Figure 4
The relationships between the expression profiles of related genes (HRGs) and OS (overall survival). In these 27 HRGs, red squares display risk factors, and the other 2 green squares display protective factors.In order to better evaluate the prognosis and survival time of patients, an optimal equation for the prognostic predictor (PP) was further conducted by multivariate Cox proportional hazard regression model. Finally, 12 genes (HDLBP, SAP30, PFKP, DPYSL4, SLC2A1, HMOX1, PGK1, ERO1A, LDHA, ENO2, SLC6A6, and TPI1) were obtained and incorporated into the final model.Among these genes, 9 genes (HDLBP, SAP30, PFKP, DPYSL4, SLC2A1, HMOX1, ERO1A, LDHA, and TPI1) were risk factors for HCC patients, and 3 genes (PGK1, ENO2, and SLC6A6) were protective factors (Table 2).Table 2
The associated HRGs with OS for multivariate Cox proportional hazard regression model.
IDCoefficientHRHR.95% CIp valueHDLBP0.50031.64921.0947-2.48460.0167SAP300.41921.52081.0984-2.10560.0116PFKP0.23451.26431.0443-1.53050.0162DPYSL40.37181.45031.0289-2.04420.0338SLC2A10.41661.51681.1191-2.05570.0072HMOX10.13111.14010.9705-1.33940.1105PGK1-0.42440.65420.4637-0.92300.0157ERO1A0.30851.36140.9634-1.92390.0803LDHA0.53511.70771.2296-2.37160.0014ENO2-0.51500.59750.4269-0.83620.0027SLC6A6-0.22780.79630.6013-1.05470.1121TPI10.36581.44171.0196-2.03870.0385HR: hazard ratio; CI: confidence interval.


### 3.4. Construction and Definition of the PP

In the multivariate Cox proportional hazard regression model,Y (survival time and status) was obtained by multiplying the contribution weight of every gene and X (expression value of each gene) and then adding them. The formula of PP expressed as follows: PP=0.5003×expressionvalueofHDLBP+0.4192×expressionvalueofSAP30+0.2345×expressionvalueofPFKP+0.3718×expressionvalueofDPYSL4+0.4166×expressionvalueofSLC2A1+0.1311×expressionvalueofHMOX1+−0.4244×expressionvalueofPGK1+0.3085×expressionvalueofERO1A+0.5351×expressionvalueofLDHA+−0.5150×expressionvalueofENO2+−0.2278×expressionvalueofSLC6A6+0.3658×expressionvalueofTPI1. It was worth noting that the coefficients of the genes were positive; the expression of these genes may shorten the OS of patients with HCC. Conversely, the negative coefficient of the gene means that these genes were beneficial for increasing the OS of HCC patients.


### 3.5. Validation the Value of PP in Evaluating Patients’ Clinical Prognosis

To validate the performance of PP in evaluating patients’ clinical prognosis, we divided patients into the high- and low-risk group with the median risk score as the cut-off point according to the risk score formula (Figure5(a)). The K-M plots were plotted to analyze the different survival time between the two groups. The K-M analysis showed that patients in the high-risk group suffered significantly worse survival than those in the low-risk group (Figure 5(b)).Figure 5
Hypoxia-related prognostic predictor (PP) of hepatocellular carcinoma patients. (a) High- and low-risk group with the median risk score as the cut-off point according to the risk score formula. Red indicates high risk and green low risk. (b) Kaplan-Meier plot represents that patients in the high-risk group had significantly shorter overall survival time than those in the low-risk group. (c) The overall survival of patients in the high- and low-risk group. The dotted line was the dividing line between high- and low-risk group of HCC patients. The right side of the dotted line represented the high-risk group, the left side represented the low-risk group, the red dots represented patients who had died, and the green dots represented patients who had still survived (d) The heat map shows the distribution of these 12 key gene expression profiles in the TCGA dataset.
(a)(b)(c)(d)All these results indicated that the prediction ability of the equation was very well. Furthermore, the distribution of patients according to their prognostic risk was also showed (Figures5(c) and 5(d)). We could see intuitively that the high-risk group had a much higher number of deaths than the low-risk group.


### 3.6. The Accuracy of PP Was Verified by Combining with Clinical Parameters

To further validate the accuracy of PP in predicting OS in HCC patients, our research performed the 1-, 3-, and 5-year ROC curves by integrating age, gender, tumor grade, and AJCC TNM stage. We got the AUC value of the 1-, 3-, and 5-year risk score were all >0.7, which were much bigger than other clinical parameters (Figure6(a)). Thus, compared with other indicators, PP had important reference value in predicting patients’ prognosis.Figure 6
The accuracy of hypoxia-related prognostic predictor (PP) in hepatocellular carcinoma patients. (a) The 1-, 3-, and 5-year risk score ROC (receiver-operating characteristic) curve analysis by integrating age, gender, tumor grade, and AJCC TNM stage. Kaplan-Meier (K-M) method was used to plot the survival curve for the training set (b) and the validation set (c). Risk score with clinicopathological features in univariate analyses (d) and in multivariate analyses (e). Red squares displayed risk factors, and green squares displayed protective factors.
(a)(b)(c)(d)(e)Furthermore, we used four-fifths of the samples that met the inclusion criteria from the TCGA data set as the training set, and the remaining one-fifth samples as the validation set. The results of the two data sets with survival analysis showed that the survival probability of the high-risk group was significantly lower than that of the low-risk group, and the difference was significant (Figures6(b) and 6(c)). Risk score remained as an independent prognostic indicator for HCC patient in univariate and multivariate analyses, after adjusting for clinicopathological features such as age, gender, tumor grade, and AJCC TNM stage (HR=1.484, 95%CI=1.342–1.642, p<0.001, and HR=1.414, 95%CI=1.258–1.588, p<0.001, Figures 6(d) and 6(e)).


### 3.7. Relationship between PP and Clinical Parameters

To further clarify the important value of PP in clinical application, we explored the relationship between risk score and clinical indicators by using the “Beeswarm” software package. The relationship between risk score of the high- and low-risk group of HCC patients and the six clinical indicators (age, gender, tumor grade, and AJCC TNM stage) was verified by the “Beeswarm” software package;p<0.05 was considered to be statistically significant. The final results showed that the risk score calculated by the PP model was consistent with the change in tumor stage (Figure 7(a)) and T stage (Figure 7(b)). The higher the patient’s risk score, the higher the corresponding tumor stage and T stage, and the patient’s prognosis was poor.Figure 7
The relationship between risk score and clinical indicators. (a) The boxplot showed the relationship between risk score and tumor stage. (b) The boxplot showed the relationship between risk score and T stage.
(a)(b)


### 3.8. Potential Therapeutic Drug Screening

In order to screen out potential therapeutic drugs, we analyzed 28 upregulated and 9 downregulated HRGs with CMap. We got the top 4 of small molecules: deferoxamine, crotamiton, and talampicillin showed a highly positive correlation with HCC; lycorine showed a highly negative correlation with HCC. They all might become effective drugs for treating HCC. The 2D chemical structural formulas of the top 4 most promising molecular drugs from PubChem were shown in Figures8(a)–8(d).Figure 8
The structural formulas of the top 4 molecular drugs: (a) deferoxamine, (b) crotamiton, (c) talampicillin, and (d) lycorine.
(a)(b)(c)(d)


### 3.9. Gene Set Cancer Analysis (GSCA)

In order to further explore the underlying molecular mechanisms of the 12 core genes involved in the PP, GSCA was used to obtain the activity of 10 cancer-related pathways. Figure9 showed the overall activity of these 12 key genes in 10 cancer-related pathways.Figure 9
The activity of 10 cancer-related pathways with GSCA.


## 3.1. Differentially Expressed HRGs

A total of 374 HCC tissue samples and 50 para-carcinoma tissue specimens were obtained for RNA-seq and clinical data from TCGA. A list of 200 HRGs involved in hypoxic regulatory pathways reported in the literature was extracted from the Hallmark hypoxic gene sets in the MSigDB database. According to the limitation of aFDR<0.05 and ∣log2FoldChange∣>1.5, our study screened 37 genes with significant differential expression, of which 28 were upregulated and 9 were downregulated (Figures 1(a) and 1(b)). Then, we visualized the expression patterns of these HRGs between HCC and normal tissues by making boxplot (Figure 1(c)). The boxplot showed expression patterns that 28 upregulated genes (ALDOA, ANKZF1, ANXA2, B4GALNT2, BCAN, CA12, COL5A1, DDIT3, DPYSL4, DTNA, EFNA3, GPC3, HOXB9, INHA, KDELR3, KIF5A, LOX, NDRG1, P4HA2, PDGFB, PFKP, PGF, PPFIA4, RRAGD, SLC2A1, SLC2A5, STC1, and STC2) and 9 downregulated genes (DCN, FBP1, FOS, MT1E, MT2A, PCK1, PLAC8, SERPINE1, and ZFP36).Figure 1
Differentially expressed HRGs between HCC and para-carcinoma tissues. (a) The volcano plot for the 200 HRGs from the TCGA data portal. 37 significantly different genes were labeled in the figure. Red indicates high expression and green low expression. (b) Hierarchical clustering of differentially expressed HRG expression levels. (c) The expression patterns of box plots displayed 37 differentially expressed HRGs in HCC types and paired nontumor samples. Each red dot represents a distinct tumor sample and green represents a nontumor sample.
(a)(b)(c)


## 3.2. Functional Annotation and Enrichment Analysis of Differentially Expressed HRGs

To further understand the biologically relevant information of these 37 differential HRGs, we performed functional enrichment and enrichment analysis. The GO term function (Figure2) and KEGG pathway enrichment (Figure 3) of these genes were reviewed, respectively. We found the results that the GO terms for biological processes were carbohydrate biosynthetic process, NADH regeneration, canonical glycolysis, glucose catabolic process to pyruvate, and so on. The detailed results of all the above gene enrichment are shown in Figure 2.Figure 2
The circle plot of the top 10 GO terms function enrichment of differentially expressed HRGs.Figure 3
The circle plot of the top 10 enriched KEGG pathways of differentially expressed HRGs.Through the enrichment analysis function in the KEGG pathway, we found many important pathways associated with these genes such as HIF-1 signaling pathway, glycolysis/gluconeogenesis, carbon metabolism, MAPK signaling pathway, and PI3K-Akt signaling pathway (Figure3). As shown in the circle plot, the change from the Z-score indicated mostly related signaling pathways were more inclined to be increased.


## 3.3. Identification of Prognostic HRGs

In order to build a model related to the prognosis of HCC, we selected patients with follow-up period of more than 1 month from this data set. Finally, we obtained 349 patients for the analysis and establishment of the prognostic predictor (PP). Prognostic HRGs were computed with univariate Cox regression model by SPSS 22.0. The relationship between the expression profiles of HRGs and OS was evaluated using TCGA data and the Hallmark hypoxic gene sets in the MSigDB database, resulting in 27 prognosis-related HRGs (HR>1 or HR<1, p<0.01) in Table 1.Table 1
The associated HRGs with OS for univariate Cox regression model in HCC patients.
IDHRHR.95% CIp valueHDLBP1.761.19-2.60<0.01SAP301.921.47-2.51<0.01NDRG11.341.17-1.54<0.01ADORA2B1.561.16-2.09<0.01PFKP1.261.13-1.42<0.01DPYSL41.551.12-2.14<0.01PFKFB31.211.07-1.37<0.01SLC2A11.621.38-1.90<0.01VEGFA1.431.14-1.80<0.01HMOX11.261.09-1.45<0.01SDC31.371.11-1.70<0.01ATP7A2.011.28-3.16<0.01PGK11.451.17-1.79<0.01ERO1A1.481.17-1.87<0.01XPNPEP12.031.32-3.12<0.01LDHA2.091.60-2.73<0.01ENO21.291.10-1.52<0.01GRHPR0.730.60-0.91<0.01UGP20.740.60-0.92<0.01SLC6A61.291.08-1.53<0.01MAP3K11.461.11-1.92<0.01BNIP3L1.431.12-1.82<0.01TPI11.611.24-2.08<0.01RRAGD1.321.07-1.62<0.01TES1.281.07-1.54<0.01SERPINE11.161.05-1.28<0.01JMJD61.771.32-2.37<0.01HR: hazard ratio; CI: confidence interval.By mapping the forest of HRGs and OS (Figure4), we can intuitively understand the role of these genes in cancer prognosis. Of these, 25 genes (HDLBP, SAP30, NDRG1, ADORA2B, PFKP, DPYSL4, PFKFB3, SLC2A1, VEGFA, HMOX1, SDC3, ATP7A, PGK1, ERO1A, XPNPEP1, LDHA, ENO2, SLC6A6, MAP3K1, BNIP3L, TPI1, RRAGD, TES, SERPINE1, and JMJD6) were risk factors for cancer, and 2 genes (GRHPR and UGP2) were protective factors.Figure 4
The relationships between the expression profiles of related genes (HRGs) and OS (overall survival). In these 27 HRGs, red squares display risk factors, and the other 2 green squares display protective factors.In order to better evaluate the prognosis and survival time of patients, an optimal equation for the prognostic predictor (PP) was further conducted by multivariate Cox proportional hazard regression model. Finally, 12 genes (HDLBP, SAP30, PFKP, DPYSL4, SLC2A1, HMOX1, PGK1, ERO1A, LDHA, ENO2, SLC6A6, and TPI1) were obtained and incorporated into the final model.Among these genes, 9 genes (HDLBP, SAP30, PFKP, DPYSL4, SLC2A1, HMOX1, ERO1A, LDHA, and TPI1) were risk factors for HCC patients, and 3 genes (PGK1, ENO2, and SLC6A6) were protective factors (Table 2).Table 2
The associated HRGs with OS for multivariate Cox proportional hazard regression model.
IDCoefficientHRHR.95% CIp valueHDLBP0.50031.64921.0947-2.48460.0167SAP300.41921.52081.0984-2.10560.0116PFKP0.23451.26431.0443-1.53050.0162DPYSL40.37181.45031.0289-2.04420.0338SLC2A10.41661.51681.1191-2.05570.0072HMOX10.13111.14010.9705-1.33940.1105PGK1-0.42440.65420.4637-0.92300.0157ERO1A0.30851.36140.9634-1.92390.0803LDHA0.53511.70771.2296-2.37160.0014ENO2-0.51500.59750.4269-0.83620.0027SLC6A6-0.22780.79630.6013-1.05470.1121TPI10.36581.44171.0196-2.03870.0385HR: hazard ratio; CI: confidence interval.


## 3.4. Construction and Definition of the PP

In the multivariate Cox proportional hazard regression model,Y (survival time and status) was obtained by multiplying the contribution weight of every gene and X (expression value of each gene) and then adding them. The formula of PP expressed as follows: PP=0.5003×expressionvalueofHDLBP+0.4192×expressionvalueofSAP30+0.2345×expressionvalueofPFKP+0.3718×expressionvalueofDPYSL4+0.4166×expressionvalueofSLC2A1+0.1311×expressionvalueofHMOX1+−0.4244×expressionvalueofPGK1+0.3085×expressionvalueofERO1A+0.5351×expressionvalueofLDHA+−0.5150×expressionvalueofENO2+−0.2278×expressionvalueofSLC6A6+0.3658×expressionvalueofTPI1. It was worth noting that the coefficients of the genes were positive; the expression of these genes may shorten the OS of patients with HCC. Conversely, the negative coefficient of the gene means that these genes were beneficial for increasing the OS of HCC patients.


## 3.5. Validation the Value of PP in Evaluating Patients’ Clinical Prognosis

To validate the performance of PP in evaluating patients’ clinical prognosis, we divided patients into the high- and low-risk group with the median risk score as the cut-off point according to the risk score formula (Figure5(a)). The K-M plots were plotted to analyze the different survival time between the two groups. The K-M analysis showed that patients in the high-risk group suffered significantly worse survival than those in the low-risk group (Figure 5(b)).Figure 5
Hypoxia-related prognostic predictor (PP) of hepatocellular carcinoma patients. (a) High- and low-risk group with the median risk score as the cut-off point according to the risk score formula. Red indicates high risk and green low risk. (b) Kaplan-Meier plot represents that patients in the high-risk group had significantly shorter overall survival time than those in the low-risk group. (c) The overall survival of patients in the high- and low-risk group. The dotted line was the dividing line between high- and low-risk group of HCC patients. The right side of the dotted line represented the high-risk group, the left side represented the low-risk group, the red dots represented patients who had died, and the green dots represented patients who had still survived (d) The heat map shows the distribution of these 12 key gene expression profiles in the TCGA dataset.
(a)(b)(c)(d)All these results indicated that the prediction ability of the equation was very well. Furthermore, the distribution of patients according to their prognostic risk was also showed (Figures5(c) and 5(d)). We could see intuitively that the high-risk group had a much higher number of deaths than the low-risk group.


## 3.6. The Accuracy of PP Was Verified by Combining with Clinical Parameters

To further validate the accuracy of PP in predicting OS in HCC patients, our research performed the 1-, 3-, and 5-year ROC curves by integrating age, gender, tumor grade, and AJCC TNM stage. We got the AUC value of the 1-, 3-, and 5-year risk score were all >0.7, which were much bigger than other clinical parameters (Figure6(a)). Thus, compared with other indicators, PP had important reference value in predicting patients’ prognosis.Figure 6
The accuracy of hypoxia-related prognostic predictor (PP) in hepatocellular carcinoma patients. (a) The 1-, 3-, and 5-year risk score ROC (receiver-operating characteristic) curve analysis by integrating age, gender, tumor grade, and AJCC TNM stage. Kaplan-Meier (K-M) method was used to plot the survival curve for the training set (b) and the validation set (c). Risk score with clinicopathological features in univariate analyses (d) and in multivariate analyses (e). Red squares displayed risk factors, and green squares displayed protective factors.
(a)(b)(c)(d)(e)Furthermore, we used four-fifths of the samples that met the inclusion criteria from the TCGA data set as the training set, and the remaining one-fifth samples as the validation set. The results of the two data sets with survival analysis showed that the survival probability of the high-risk group was significantly lower than that of the low-risk group, and the difference was significant (Figures6(b) and 6(c)). Risk score remained as an independent prognostic indicator for HCC patient in univariate and multivariate analyses, after adjusting for clinicopathological features such as age, gender, tumor grade, and AJCC TNM stage (HR=1.484, 95%CI=1.342–1.642, p<0.001, and HR=1.414, 95%CI=1.258–1.588, p<0.001, Figures 6(d) and 6(e)).


## 3.7. Relationship between PP and Clinical Parameters

To further clarify the important value of PP in clinical application, we explored the relationship between risk score and clinical indicators by using the “Beeswarm” software package. The relationship between risk score of the high- and low-risk group of HCC patients and the six clinical indicators (age, gender, tumor grade, and AJCC TNM stage) was verified by the “Beeswarm” software package;p<0.05 was considered to be statistically significant. The final results showed that the risk score calculated by the PP model was consistent with the change in tumor stage (Figure 7(a)) and T stage (Figure 7(b)). The higher the patient’s risk score, the higher the corresponding tumor stage and T stage, and the patient’s prognosis was poor.Figure 7
The relationship between risk score and clinical indicators. (a) The boxplot showed the relationship between risk score and tumor stage. (b) The boxplot showed the relationship between risk score and T stage.
(a)(b)


## 3.8. Potential Therapeutic Drug Screening

In order to screen out potential therapeutic drugs, we analyzed 28 upregulated and 9 downregulated HRGs with CMap. We got the top 4 of small molecules: deferoxamine, crotamiton, and talampicillin showed a highly positive correlation with HCC; lycorine showed a highly negative correlation with HCC. They all might become effective drugs for treating HCC. The 2D chemical structural formulas of the top 4 most promising molecular drugs from PubChem were shown in Figures8(a)–8(d).Figure 8
The structural formulas of the top 4 molecular drugs: (a) deferoxamine, (b) crotamiton, (c) talampicillin, and (d) lycorine.
(a)(b)(c)(d)


## 3.9. Gene Set Cancer Analysis (GSCA)

In order to further explore the underlying molecular mechanisms of the 12 core genes involved in the PP, GSCA was used to obtain the activity of 10 cancer-related pathways. Figure9 showed the overall activity of these 12 key genes in 10 cancer-related pathways.Figure 9
The activity of 10 cancer-related pathways with GSCA.


## 4. Discussion

HCC is a major contributor of death caused by cancer. Previously, many progresses have been made in understanding the epidemiology, risk factors, and molecular profiles of HCC [16]. However, incidence and HCC-specific mortality still continue to increase [17–19]. Although some progress has been made in molecular targeted therapy for HCC, the results are still unsatisfactory [20]. It is urgently needed to better understand the molecular mechanism which leads to the development of HCC and to search for targeted genes that play a key role in the prognosis of HCC for early intervention.Exploration of hypoxia mechanism opens new perspectives for HCC [21]. In particular, the announcement of the 2019 Nobel Prize in physiology or medicine has increased researchers’ enthusiasm for exploring the mechanism of hypoxia in the field of cancer. Increasing evidence indicates that hypoxia is closely related to the migration and malignant progression of HCC [8, 22]. However, most studies have only focused on hypoxia by studying signaling genes [23, 24]. With the development of bioinformatics and the continuous improvement of high-throughput sequencing technology, some large-scale databases, including TCGA and GEO, have emerged. These databases provide effective means for sorting out gene signatures. In the current study, we deeply mined all the differential genes which are significantly expressed in tumors and normal tissues. By integrating them with a large sample of clinical data, we aimed to select molecular biomarkers for detecting the prognosis of HCC patients.We first screened 37 differentially expressed HRGs between HCC and para-carcinoma tissues. To discover the role of genes in HCC progression, GO and KEGG analysis of the differential expression genes was performed. Gene functional enrichment analysis suggested that these genes were mainly involved in carbohydrate biosynthetic process and HIF-1 signaling pathway. The process of carbohydrate biosynthesis was mainly related to the energy metabolism of cancer cells. HIF-1 and its related genes actively promoted HCC growth, HCC cell proliferation, and aggressive behavior, which was positively associated with the presence of intrahepatic metastasis and the histological grade of HCC [8, 25]. Our analysis confirmed that HIF-1 signaling pathway was the most critical process for hypoxia leading to poor prognosis of HCC, which is consistent with the findings in most literatures [26, 27]. As such, these genes could also act as clinical biomarkers for monitoring metastasis, assessing survival, and potential drug targets. Bioinformatics analysis provided several clues intervening the occurrence and development of HCC via several hypoxic pathways. The univariate Cox regression analysis revealed that 27 HRGs were associated with OS in the TCGA database. Further multivariate Cox regression analysis helped us determine 12 vital prognostic HRGs (HDLBP, SAP30, PFKP, DPYSL4, SLC2A1, HMOX1, PGK1, ERO1A, LDHA, ENO2, SLC6A6, and TPI1) to construct of the PP. There are already several existed reports about partial of these genes.HDLBP (high-density lipoprotein-binding protein), a positively regulated gene, is closely related to the cancer process [28]. HDLBP can promote HCC cell proliferation and tumorigenesis. Consistent with many other studies, HDLBP was overexpressed in other types of cancers [29]. Intriguingly, in breast cancer, HDLBP may be a tumor suppressor to accelerate the degradation and inhibit the translation of the c-fms protooncogene mRNA [30]. These observations suggested that the mechanism of HDLBP (either inhibit or promote carcinogenesis) was still unclear. SAP30 (serum amyloid P30) was a very sensitive indicator of liver disease [31]. It has been demonstrated that SAP upregulated in serum samples from HCC patients. HMOX1 (heme oxygenase 1) is an isoenzyme of heme oxygenase (HO). HO is a microsomal enzyme system that degrades heme. Their main biological role is to degrade heme to produce carbon monoxide (CO) and bilivertin and release iron ions. HMOX1 has been proved to be an oncogene and plays an important role in the development of various cancers. Inhibition of HMOX1 degradation can definitely promote HCC proliferation. LDHA (lactate dehydrogenase A) is a vital enzyme responsible for cancer growth and energy metabolism via the aerobic glycolytic pathway. Inhibition of LDHA could inhibit tumor-initiating cell survival and proliferation, which indicated that LDHA may be a potential therapeutics target [32]. PGK1 (phosphoglycerate kinase 1) is an important enzyme in the metabolic glycolysis pathway. Many articles have observed a significant overexpression of PGK1 in HCC tissues and a negative correlation between PGK1 expression and HCC patient survival. Depletion of PGK1 dramatically reduced cancer cell proliferation and tumorigenesis [33, 34]. These findings demonstrated a novel role of PGK1 in HCC progression. ENO2 (Enolase 2) encodes an enolase isoenzyme which is regarded as a sensitive biomarker for neuroendocrine tumors [35]. Additionally, ENO2 was believed to be an indicator in the diagnosis and prognosis rather than a potential target for therapy in renal tumor patients [36]. In our study, ENO2 was defined as a negative regulatory gene. It provided new therapeutic targets and biomarkers for HCC prognosis through inhibiting its expression to improve the prognosis of HCC. Other genes were rarely reported; their roles and underlying mechanism in cancer still need to be further explored.According to the current study, some features of HCC prognosis based on expression profile have been found by large public databases. However, the purpose of these studies was often limited to mining new molecular markers but neglecting routine clinical parameters [37, 38]. Thus, we integrated patient clinical information to establish a PP model for assessing prognosis. Firstly, we verified the accuracy of the model through the K-M survival analysis. The risk score was calculated by PP formula, and the median was taken as the cut-off point dividing the high- and low-risk group, and the survival difference was significant (Figure 5(a)). Secondly, we performed the 1-, 3-, and 5-year risk score ROC curve analysis by integrating age, gender, tumor grade, and AJCC TNM stage and found that the AUC value of 1-, 3-, and 5-year risk score were all much bigger than other clinical parameters (Figure 6(a)). Besides, we divided the TCGA data set into a training set and a validation set and performed survival analysis on the hypoxia-related signature in the two data sets to illustrate the robustness of the PP. The results of the two data sets with survival analysis showed that the survival probability of the high-risk group was significantly lower than that of the low-risk group, and the difference was significant (Figures 6(b) and 6(c)). Similarly, we found that the risk score was still an independent indicator for evaluating the prognosis of patients with HCC through the process of univariate and multivariate analysis (Figures 6(d) and 6(e)). At last, by comprehensively analyzing the risk score and the patient’s tumor TNM stage, we found that the greater the risk score calculated by PP, the higher the patient’s tumor grade and the worse the prognosis. All of the above fully demonstrated the accuracy of our constructed model of PP and its important value in evaluating patients’ prognosis.The results of CMap demonstrated the therapeutic agents closely associated with HCC. Substantial evidence suggests that iron overload enhances cancer growth and metastasis. As a commonly used iron chelator, high-dose deferoxamine therapy could inhibit tumor migration by destroying iron homeostasis, arresting the cell cycle and enhancing apoptosis of hepatocellular carcinoma and hepatoblastoma cell lines [39]. The iron chelator deferoxamine has been proved that it could prevent liver injury as well as the development of preneoplastic lesions in rats. In particular, as a potential anticancer drug, it has been used in clinical studies to confirm its efficacy in patients with advanced hepatocellular carcinoma [40]. Lycorine is a multifunctional bioactive compound and powerful anticancer agent against various cancer cell lines. Focusing on hepatocellular carcinoma, lycorine can significantly induce apoptosis and autophagy of hepatocellular carcinoma cells both in vivo and in vitro [41]. The possible molecular mechanism involved the TCRP1/Akt/mTOR pathway. Lycorine inhibited TCRP1 to reduce the phosphorylation level of Akt and inhibited Akt/mTOR signal transduction and finally mediated the apoptosis and autophagy of xenograft hepatocellular tumors without remarkable toxicity. Similarly, research by Hu et al. suggested that lycorine could be used to treat colorectal cancer by inducing autophagy-associated apoptosis [42]. Additionally, lycorine could inhibit the proliferation of SGC-7901 gastric cancer cells and promote cell apoptosis [43]. The specific mechanism of treatment was exactly consistent with the results of GSCA for modeling gene pathway analysis. However, there were few reports about crotamiton and talampicillin in the treatment of cancer. Their anticancer mechanisms need to be further explored and studied.In conclusion, through a comprehensive bioinformatics analysis with HRG expression profiles, our study firstly reported 12 prognostic HRGs and constructed the model of the PP. This comprehensive research of multiple databases helps us gain insight into the biological properties of HCC and provides deferoxamine, crotamiton, talampicillin, and lycorine as potential molecular drugs to fight against HCC.



---

*Source: 7928051-2021-02-08.xml*