Identification of prognostic biomarkers for patients with hepatocellular carcinoma after hepatectomy.

Hepatocellular carcinoma (HCC) is a lethal malignancy with high morbidity and mortality rates worldwide. The identification of prognosis‑associated biomarkers is crucial to improve HCC patient survival. The present study aimed to explore potential predictive biomarkers for HCC. Differentially expressed genes (DEGs) were analyzed in the GSE36376 dataset using GEO2R. Hub genes were identified and further investigated for prognostic value in HCC patients. A risk score model and nomogram were constructed to predict HCC prognosis using the prognosis‑associated genes and clinical factors. Pearson's correlation was employed to show interactions among hub genes. Gene enrichment analysis was performed to identify detailed biological processes and pathways. A total of 71 DEGs were obtained and seven (ADH4, CYP2C8, CYP2C9, CYP8B1, SLC22A1, TAT and HSD17B13, all adjusted P≤0.05) of the 10 hub genes were identified as prognosis‑related genes for survival analysis in HCC patients, including alcohol dehydrogenase 4 (class II), pi polypeptide (ADH4), cytochrome p450 family 2 subfamily C member 8 (CYP2C8), cytochrome P450 family 2 subfamily C member 9 (CYP2C9), cytochrome P450 family 8 subfamily B member 1 (CYP8B1), solute carrier family 22 member 1 (SLC22A1), tyrosine aminotransferase (TAT) and hydroxysteroid 17‑β dehydrogenase 13 (HSD17B13). The risk score model could predict HCC prognosis and the nomogram visualized gene expression and clinical factors of probability for HCC prognosis. The majority of genes showed significant Pearson's correlations with others (41 Pearson correlations P≤0.01, four Pearson correlations P>0.05). GO analysis revealed that terms such as 'chemical carcinogenesis' and 'drug metabolism‑cytochrome P450' were enriched and may prove helpful to elucidate the mechanisms of hepatocarcinogenesis. Hub genes ADH4, CYP2C8, CYP2C9, CYP8B1, SLC22A1, TAT and HSD17B13 may be useful as predictive biomarkers for HCC prognosis.

Introduction

Liver cancer is more common in men than in women, and is the second leading cause of cancer-associated mortality worldwide in men living in less developed countries (1). It is estimated that ~50% of newly diagnosed cancer cases and mortalities globally occurred in China alone in 2012, which equates to approximately 391,250 new liver cancer cases and 372,750 mortalities (1). Hepatocellular carcinoma (HCC) accounts for the majority (70–90%) of primary liver cancer cases worldwide (2). Previous epidemiological surveys have revealed that chronic hepatitis B or C viral infection, cirrhosis, exposure to aflatoxin B1, obesity, chronic alcohol consumption and diabetes mellitus, as well as metabolic abnormalities including haemochromatosis and α1-antitrypsin deficiency, are common and significant risk factors for HCC development (3–5). Furthermore, recent studies indicate that approximately 1,000,000 new HCC cases are diagnosed each year worldwide, with the same incidence and morbidity rate, indicating that HCC diagnosis is typically at the advanced stage and prognosis remains poor (6,7). Despite advances in diagnosis, prevention and treatment, including ultrasonography, multiphase computerized tomography, magnetic resonance imaging, surgical resection, liver transplantation, transarterial chemoembolization, radiofrequency ablation and transarterial radiation, percutaneous ablation and systemic therapy, the prognosis for HCC patients remains unsatisfactory (8). The 5-year relative survival rate is approximately 7% (9). Therefore, identification of novel biomarkers for the early diagnosis of HCC is vital and may improve HCC prognosis.

At present, many reports have focused on the identification of prognosis-associated biomarkers The role of astrocytic phosphoprotein PEA-15 in HCC has been evaluated and may be a novel target for HCC treatment, as well as a predictive biomarker for HCC patient prognosis (10). Furthermore, it has been suggested that CKLF-like MARVEL transmembrane domain-containing protein 5 may function as a tumor suppressor in human HCC and represent a valuable therapeutic target (11). microRNA (miR)-182-5p is recognized as a potential predictive biomarker for the early recurrence of HCC (12). Currently, DNA microarray approach has been used to investigate the genetic features of HCC in molecular biology (13). Accumulating evidence regarding gene expression in HCC has demonstrated that a variety of differentially expressed genes (DEGs) may be involved in the process of hepatocarcinogenesis (14). The combination of microarray techniques and bioinformatic analysis makes it conceivable to use DEGs in one or more chips to detect potential predictive biomarkers for several types of malignancies (15). Therefore, the present study focused on DEGs in microarrays and aimed to identify potential predictive biomarkers for patients with HCC.

Materials and methods

Data collection and processing

GEO2R (www.ncbi.nlm.nih.gov/geo/geo2r/) was initially used to identify DEGs in the GSE36376 gene expression profile in the Gene Expression Omnibus (GEO) database (16). The GSE36376 dataset (www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36376) is embedded in the platform GPL10558 (Illumina HumanHT-12 v4 Expression Beadchip) and includes 240 tumor tissues and 193 adjacent non-tumor tissues of two gene expression by array sets: A training set and validation set (17).

Next, DEGs from the GSE36376 profile were acquired using the GEO2R tool with the criteria of |log fold change|≥2 and P≤0.05. Furthermore, the Cytoscape software (version 3.6.0) plugin, CentiScaPe (version 2.2), was employed to identify the top 10 hub genes in these DEGs on the basis of the following centralities: degree, betweenness and closeness (18).

Hub gene expression, survival and stratified analysis

The expression levels of 10 hub genes in multiple organs, as well as in tumor and non-tumor tissues were obtained from the GTEx portal (www.gtexportal.org/home) and the Gene Expression Profiling Interactive Analysis (GEPIA; gepia.cancer-pku.cn/index.html) server (19). Protein expression levels of the hub genes in liver tissue were obtained from The Human Protein Atlas (www.proteinatlas.org) (20). A co-expression matrix was constructed using R (version 3.5.0; www.r-project.org) to show Pearson correlations between two of these genes.

In addition, hub genes were further analyzed for their prognostic values using The Cancer Genome Atlas database (cancergenome.nih.gov/). The clinical characteristics of 360 HCC patients, including age, sex, race and tumor stage, were obtained and used in the analysis. Gene expression data were downloaded from the OncoLnc website (www.oncolnc.org) at median cut off. In addition, clinical data with statistically significant P-values (≤0.05) were adjusted for further analysis to identify prognosis-associated genes. For the above identified genes, clinical data were stratified for further analysis.

Hub gene mutational and transcriptional analysis

The mutation status of hub genes using the cBioPortal website (www.cbioportal.org/) (21,22). Furthermore, transcripts per million (TPM) of the hub genes in liver tissues were analyzed to identify transcription levels at the log scale [Log2 (TPM+1)] using the GEPIA website (gepia.cancer-pku.cn/index.html/) (19).

Construction of risk score model and nomogram

A risk score model was constructed based on the expression levels of prognosis-associated genes and the contribution coefficient (β) of the multivariate Cox proportional hazards regression model. The formula of the model was as follows: risk score=expression of gene1 × β1 of gene1 + expression of gene2 × β2 of gene2 + … + expression of genen × βn of genen. The risk score was then divided into high and low risk groups with the same cut-off criteria of gene expression. A Kaplan-Meier survival curve was drawn to predict patient survival. Prognostic receiver operating characteristic (ROC) curves were then created at 1–5 years.

A nomogram was also constructed based on seven prognosis-associated hub genes, clinical factors and tumor stage, to obtain patient survival at 1, 3 and 5 years. The contribution of each factor was limited to a maximum of 100 points.

Gene-gene and protein-protein interactions (PPI) network construction and functional enrichment analysis

To identify the biological processes and metabolic pathways of these hub genes, enrichment analysis of Gene Ontology (GO) was performed, including biological process (BP), cellular component (CC), and molecular function (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using The Database for Annotation, Visualization and Integrated Discovery (DAVID; david.ncifcrf.gov/) (15,23). In addition, visualized GO results were presented by the Cytoscape BinGO plugin version 3.6.0 (18,24). Gene-gene interaction networks were constructed using the GeneMANIA plugin (version 3.4.1) in Cytoscape software version 3.6.0 (25). The PPI network depicted interactions among these proteins and was constructed using the Search Tool for the Retrieval of Interacting Genes website (string-db.org/cgi/input.pl) (26).

Statistical analysis

Kaplan-Meier survival analysis by log rank test was used to calculate the median survival time (MST). Univariate and multivariate Cox proportional hazards models were used to calculate the hazard ratio (HR) and 95% confidence interval (CI). P≤0.05 was considered to indicate a statistically significant difference. Box plots by unpaired t-test, survival curves and ROC curves including area under the curve (AUC) were produced using GraphPad version 7.0 (GraphPad Software, Inc., La Jolla, CA, USA). Statistical analysis was performed using SPSS software version 16.0 (SPPS, Inc., Chicago, IL, USA).

Results

Processing of DEGs and hub genes

A total of 71 DEGs were obtained from the GSE36376 dataset and the top 10 hub genes were identified. These 10 hub genes included cytochrome P450 family 2 subfamily E member 1 (CYP2E1), tyrosine aminotransferase (TAT), cytochrome P450 family 2 subfamily A member 6 (CYP2A6), cytochrome P450 family 8 subfamily B member 1 (CYP8B1), cytochrome P450 family 2 subfamily C member 9 (CYP2C9), hydroxysteroid 11-β dehydrogenase 1 (HSD11B1), hydroxysteroid 17-β dehydrogenase 13 (HSD17B13), solute carrier family 22 member 1 (SLC22A1), cytochrome p450 family 2 subfamily C member 8 (CYP2C8) and alcohol dehydrogenase 4 (class II), pi polypeptide (ADH4). Detailed results are shown in Table I.

Table I.

Identified hub genes in the GSE36376 dataset and their screening centralities.

Genes	Degree centralities	Betweenness centralities	Closeness centralities
CYP2E1	37	362.1080892	0.009345794
TAT	32	266.8918244	0.008928571
CYP2A6	32	97.3346815	0.008064516
CYP8B1	32	109.0359712	0.008064516
CYP2C9	31	179.8498746	0.007874016
HSD11B1	31	333.8964087	0.008849558
HSD17B13	31	147.7138107	0.008403361
SLC22A1	30	84.7381197	0.008403361
CYP2C8	29	36.7257238	0.007751938
ADH4	28	52.5113332	0.007575758

CYP2E1, cytochrome P450 family 2 subfamily E member 1; TAT, tyrosine aminotransferase; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Analysis of hub gene expression, mutation and transcription

Hub gene expression in tumor tissues was downregulated compared with normal tissues in all cases (Fig. 1). However, only CYP2A6, CYP2C8, CYP2E1, HSD17B13 and SLC22A1 had a significant alteration in expression (P<0.05). Tissue expression data indicated that all of the hub genes were highly expressed in liver tissue (Fig. 2). Mutation analysis of hub genes revealed that mutations were present in different ratios (Fig. 3A). At 11%, HSD11B1 was the most significant in terms of mutation ratio, of which the majority were amplification mutations. Transcriptional analysis (Fig. 3B) demonstrated that all hub genes were differentially transcribed in tumor and normal tissues. Furthermore, normal tissues had consistently high TPM in the 10 hub genes. Detailed results are shown in Fig. 3.

Figure 1.

Expression levels of hub genes ADH4, CYP2A6, CYP2C8, CYP2C9, CYP2E1, CYP8B1, HSD11B1, HSD17B13, SLC22A1 and TAT in tumor and normal tissues. *P<0.05 vs. normal tissue. LIHC, liver hepatocellular carcinoma; CYP2E1, cytochrome P450 family 2 subfamily E member 1; TAT, tyrosine aminotransferase; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Figure 2.

Expression levels of hub genes (A) ADH4, (B) CYP2A6, (C) CYP2C8, (D) CYP2C9, (E) CYP2E1, (F) CYP8B1, (G) HSD11B1, (H) HSD17B13, (I) SLC22A1 and (J) TAT in different organs. High liver expression was noted in all hub genes. TPM, transcripts per million; CYP2E1, cytochrome P450 family 2 subfamily E member 1; TAT, tyrosine aminotransferase; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Figure 3.

Mutational and transcriptional analysis of hub genes. (A) Mutational analysis including missense mutation, deep deletion and amplification of the 10 hub genes. (B) Transcriptional analyses (TPM) of the 10 hub genes in tumor and normal tissues. TPM, transcripts per million; CYP2E1, cytochrome P450 family 2 subfamily E member 1; TAT, tyrosine aminotransferase; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Protein expression data indicated that eight proteins were highly expressed in liver tissue, excluding TAT and CYP8B1 (data not shown; Fig. 4), which was similar to the tissue expression results.

Figure 4.

Immunohistochemistry results of protein expression levels of hub genes ADH4, CYP2A6, CYP2C8, CYP2C9, CYP2E1, HSD11B1, HSD17B13 and SLC22A1 in normal liver tissues. ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP2C8, cytochrome p450 family 2 subfamily C member 8; CYP2C9, cytochrome P450 family 2 subfamily C member 9; CYP2E1, cytochrome P450 family 2 subfamily E member 1; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1.

Demographic and clinicopathological characteristics

A total of 360 HCC patients were included in the dataset. Survival analysis indicated that tumor stage showed statistical significance (log-rank P<0.0001), but all other factors were not significant (log-rank P>0.05; Table II).

Table II.

Demographic characteristics of HCC patients.

Variables	Patients (n=360)	Number (%)	MST (days)	HR (95% CI)	Log-rank P-value
Race
Asian	155	44 (28.4)	NA	Ref.
Non-Asian	196	78 (39.8)	1,397	1.29 (0.89–1.87)	0.184
Sex
Male	244	78 (32.0)	2,486	Ref.
Female	116	48 (41.4)	1,560	1.21 (0.84–1.73)	0.308
Age (years)
≤61	186	59 (31.7)	2,116	Ref.
>61	171	65 (38.0)	1,622	1.18 (0.83–1.69)	0.349
Tumor stage
I and II	252	66 (26.2)	2,532	Ref.
III and IV	87	48 (55.2)	770	2.50 (1.72–3.63)	<0.0001

Data were unavailable for race in 9 patients, age in 3 patients and tumor stage in 21 patients. MST, median survival time; HR, hazard ratio; CI, confidence interval; NA, not available; Ref., reference; Tumor stage was determined with the Tumor, Node and Metastasis staging system version 7.

Expression levels and survival analysis of hub genes

In the univariate survival analysis, ADH4, CYP2C8, CYP2C9, CYP2E1, CYP8B1, HSD17B13, SLC22A1 and TAT were significant (P≤0.05; Table III and Fig. 5). Following adjustment for tumor stage, ADH4, CYP2C8, CYP2C9, CYP8B1, SLC22A1, TAT and HSD17B13 were significant (P≤0.05; Table III) while CYP2A6, CYP2E1, HSD11B1 were not significant (P>0.05). All hub genes were significantly different when comparing high and low expression levels (P<0.0001; Fig. 6A).

Table III.

Survival analysis of HCC patient prognosis.

Gene expression	Patients/events	MST (days)	Crude HR (95% CI)	Crude P-value	Adjusted HR (95% CI)	Adjusted P-value
ADH4
Low	180/80	1,372	Ref.		Ref.
High	180/46	2,456	0.52 (0.36–0.75)	<0.001[a]	0.55 (0.38–0.81)	0.002[a]
CYP2A6
Low	180/66	1,694	Ref.		Ref.
High	180/60	1,791	0.83 (0.58–1.18)	0.291	0.83 (0.57–1.20)	0.321
CYP2C8
Low	180/75	1,229	Ref.		Ref.
High	180/51	2,456	0.56 0.39–0.79)	0.001[a]	0.56 (0.38–0.83)	0.003[a]
CYP2C9
Low	180/74	1,271	Ref.		Ref.
High	180/52	2,456	0.56 (0.39–0.80)	0.001[a]	0.64 (0.43–0.93)	0.020[a]
CYP2E1
Low	180/75	1,490	Ref.		Ref.
High	180/51	2,456	0.67 (0.47–0.96)	0.029[a]	0.74 (0.51–1.08)	0.119
CYP8B1
Low	180/73	1,372	Ref.		Ref.
High	180/53	2,131	0.63 (0.44–0.90)	0.010[a]	0.61 (0.42–0.89)	0.011[a]
HSD11B1
Low	180/71	1,397	Ref.		Ref.
High	180/55	2,131	0.77 (0.54–1.09)	0.142	0.81 (0.56–1.18)	0.279
SLC22A1
Low	180/79	1,149	Ref.		Ref.
High	180/47	2,456	0.49 (0.34–0.70)	<0.0001[a]	0.51 (0.35–0.75)	0.001[a]
TAT
Low	180/74	1,372	Ref.		Ref.
High	180/52	2,131	0.56 (0.39–0.80)	0.001[a]	0.53 (0.36–0.78)	0.001[a]
HSD17B13
Low	180/74	1,372	Ref.		Ref.
High	180/52	2,456	0.57 (0.40–0.81)	0.001[a]	0.56 (0.39–0.82)	0.003[a]

P<0.05. MST, median survival time; HR, hazard ratio; CI, confidence interval; CYP2E1, cytochrome P450 family 2 subfamily E member 1; TAT, tyrosine aminotransferase; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Figure 5.

Kaplan-Meier plots of hub genes ADH4, CYP2A6, CYP2C8, CYP2C9, CYP2E1, CYP8B1, HSD11B1, HSD17B13, SLC22A1 and TAT. CYP2E1, cytochrome P450 family 2 subfamily E member 1; TAT, tyrosine aminotransferase; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Figure 6.

Expression levels, Pearson's correlations and expression network analysis of hub genes. (A) Relative mRNA expression levels in low expression and high expression groups. (B) Pearson correlation analysis of hub genes. Blue indicates positive correlation and red indicates negative correlation. (C) Protein-protein interaction network and (D) Gene-gene co-expression network of hub genes. **P<0.01, ***P<0.001, ****P<0.0001. CYP2E1, cytochrome P450 family 2 subfamily E member 1; TAT, tyrosine aminotransferase; CYP2A6, cytochrome P450 family 2 subfamily A member 6; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD11B1, hydroxysteroid 11-β dehydrogenase 1; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Stratified analysis of prognosis-associated genes

In the stratification of tumor stage, high expression levels of CYP2C8, CYP2C9, SLC22A1, TAT and HSD17B13 had tumor suppressor roles in stages I and II (P≤0.05) while high expression levels of ADH4, CYP2C8, CYP8B1, SLC22A1, TAT and HSD17B13 had tumor suppressor roles in stages III and IV (P≤0.05). Detailed results are presented in Tables IV and V.

Table IV.

Stratified analysis of HSD17B13, SLC22A1 and TAT for HCC patients (n=360) in terms of prognosis.

	HSD17B13				SLC22A1				TAT
Characteristics	Low	High	Adjusted HR (95% CI)	Adjusted P-value	Low	High	Adjusted HR (95% CI)	Adjusted P-value	Low	High	Adjusted HR (95% CI)	Adjusted P-value
Race
Asian	87	68	0.48 (0.24–0.93)	0.03[a]	88	67	0.32 (0.15–0.68)	0.003[a]	92	63	0.38 (0.19–0.78)	0.009[a]
Non-Asian	89	107	0.55 (0.34–0.89)	0.014[a]	89	107	0.64 (0.40–1.05)	0.076	82	114	0.66 (0.41–1.07)	0.092
Sex
Male	128	116	0.47 (0.28–0.78)	0.003[a]	111	133	0.48 (0.29–0.79)	0.004[a]	115	129	0.44 (0.27–0.72)	0.001[a]
Female	52	64	0.71 (0.39–1.28)	0.258	69	47	0.60 (0.31–1.13)	0.113	65	51	0.79 (0.43–1.45)	0.45
Age (years)
≤61	100	86	0.42 (0.24–0.74)	0.003[a]	104	82	0.29 (0.16–0.55)	<0.001[a]	103	83	0.37 (0.21–0.67)	0.001[a]
>61	80	91	0.66 (0.39–1.12)	0.121	74	97	0.76 (0.45–1.31)	0.325	74	97	0.73 (0.43–1.24)	0.244
Tumor stage
I and II	121	131	0.61 (0.37–0.99)	0.047[a]	116	136	0.55 (0.34–0.91)	0.019[a]	121	131	0.61 (0.37–1.00)	0.048[a]
III and IV	50	37	0.47 (0.26–0.86)	0.015[a]	58	29	0.40 (0.20–0.79)	0.009[a]	53	34	0.41 (0.21–0.77)	0.006[a]

Data were unavailable for race in 9 patients, age in 3 patients and tumor stage in 21 patients. Tumor stage was determined with the Tumor, Node and Metastasis staging system version 7.

P<0.05. HR, hazard ratio; 95% CI, 95% confidence interval.

Table V.

Stratified analysis of ADH4, CYP2C8, CYP2C9 and CYP8B1 for HCC patients (n=360) in terms of prognosis.

	ADH4				CYP2C8				CYP2C9				CYP8B1
Characteristics	Low	High	Adjusted HR (95% CI)	Adjusted P-value	Low	High	Adjusted HR (95% CI)	Adjusted P-value	Low	High	Adjusted HR (95% CI)	Adjusted P-value	Low	High	Adjusted HR (95% CI)	Adjusted P-value
Race
Asian	91	64	0.35 (0.17–0.74)	0.005[a]	85	70	0.40 (0.20–0.80)	0.009[a]	79	76	0.40 (0.21–0.75)	0.004[a]	93	62	0.28 (0.13–0.59)	0.001[a]
Non-Asian	85	111	0.62 (0.38–1.01)	0.057	92	104	0.70 (0.43–1.13)	0.143	97	99	0.75 (0.45–1.23)	0.251	84	112	0.90 (0.55–1.45)	0.658
Sex
Male	107	137	0.51 (0.31–0.83)	0.006[a]	120	124	0.48 (0.29–0.80)	0.005[a]	108	136	0.54 (0.33–0.89)	0.015[a]	109	135	0.48 (0.30–0.78)	0.003[a]
Female	73	43	0.65 (0.34–1.25)	0.196	60	56	0.75 (0.42–1.37)	0.355	72	44	0.85 (0.46–1.58)	0.613	71	45	1.00 (0.54–1.85)	0.998
Age (years)
≤61	103	83	0.43 (0.24–0.77)	0.005[a]	106	80	0.54 (0.30–0.98)	0.044[a]	98	88	0.73 (0.42–1.27)	0.269	115	71	0.27 (0.14–0.53)	<0.001[a]
>61	75	96	0.63 (0.36–1.09)	0.096	71	100	0.63 (0.37–1.06)	0.08	79	92	0.62 (0.37–1.06)	0.082	63	108	0.99 (0.58–1.69)	0.973
Tumor stage
I and II	117	135	0.62 (0.38–1.01)	0.055	113	139	0.58 (0.36–0.95)	0.029[a]	113	139	0.61 (0.38–1.00)	0.048[a]	121	131	0.79 (0.48–1.28)	0.334
III and IV	53	34	0.45 (0.24–0.84)	0.013[a]	59	28	0.50 (0.27–0.95)	0.033[a]	56	31	0.65 (0.35–1.20)	0.171	51	36	0.38 (0.20–0.71)	0.003[a]

Data were unavailable for race in 9 patients, age in 3 patients and tumor stage in 21 patients. Tumor stage was determined with the Tumor, Node and Metastasis staging system version 7.

P<0.05. HR, hazard ratio; 95% CI, 95% confidence interval.

Hub gene co-expression and Pearson correlation analysis

The majority of hub genes showed significant Pearson correlations with other genes. For example, CYP2E1 was positively correlated with CYP2C8, CYP2C9, CYP8B1, ADH4 and HSD11B1, and negatively correlated with CYP2A6 and HSD17B13. Detailed results are presented in Fig. 6B.

In the PPI network, most of the hub genes exhibited complicated interactive co-expression relationships, with the exception of HSD11B1, which was co-expressed with HSD17B13 alone (Fig. 6C). In the gene-gene interaction network, each gene was co-expressed with at least two other genes (Fig. 6D).

Risk score model and nomogram construction

A risk score model was constructed using the aforementioned formula, which contained risk score ranking, survival status and heat maps of gene expressions (Fig. 7A). The Kaplan-Meier plot revealed that the percent survival difference between the low and high risk groups was significant (Fig. 7B). ROC curves were then constructed to evaluate the prognostic values of the model. In the 1–5 year ROC curves, all AUCs were above 0.6 (Fig. 7C), which shows that this model was useful for prognosis prediction. In the comparison of high and low risk score groups, all of the comparisons showed significant P-values (all P<0.001; Fig. 7D).

Figure 7.

Risk score model, survival plot, ROC curves and boxplot. (A) Risk score model plot including risk score ranking, survival status and heatmap. (B) Kaplan-Meier plot of the risk score model. (C) ROC curves for 1-, 2-, 3-, 4- and 5-year survival rates from the risk score model. (D) Boxplot of risk score groups. ****P<0.0001. ROC, receiver operating characteristic; TAT, tyrosine aminotransferase; CYP8B1, cytochrome P450 family 8 subfamily B member 1; CYP2C9, cytochrome P450 family 2 subfamily C member 9; HSD17B13, hydroxysteroid 17-β dehydrogenase 13; SLC22A1, solute carrier family 22 member 1; CYP2C8, cytochrome p450 family 2 subfamily C member 8; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

The contributions of each factor were present in the nomogram (Fig. 8). In detail, tumor stages III and IV showed the maximum 100 points. Unlike the other genes, high expression levels of CYP8B1 had a high number of points. High points typically indicated low survival probability at 1, 3 and 5 years. As expected, there was a high probability of survival prediction at 1 year compared with 5 years.

Figure 8.

Nomogram for predicting HCC prognosis. TAT, tyrosine aminotransferase; CYP2C8, cytochrome p450 family 2 subfamily C member 8; CYP2C9, cytochrome P450 family 2 subfamily C member 9; CYP2C8, cytochrome p450 family 2 subfamily C member 8; CYP8B1, cytochrome P450 family 8 subfamily B member 1; SLC22A1, solute carrier family 22 member 1; ADH4, alcohol dehydrogenase 4 (class II), pi polypeptide.

Hub genes enrichment analysis

GO analysis revealed that genes were significantly enriched in the terms such as ‘oxidation-reduction process’, the ‘epoxygenase P450 pathway’, ‘metbolism of xenobiotics by cytochrome P450’, ‘chemical carcinogenesis’, ‘drug metabolic processes’ and ‘organelle membranes’ (Fig. 9).

Figure 9.

Enriched Gene Ontology and KEGG pathways of hub genes. Genes were analyzed in terms of biological process, cellular component and molecular function. Significantly enriched pathways were identified in the KEGG pathway analysis. KEGG, Kyoto Encyclopedia of Genes and Genomes.

BP terms, including ‘exogenous drug catabolic process’, ‘drug catabolic process’, ‘secondary metabolic process’ and ‘oxidation reduction’, were significantly enriched (Fig. 10). In CC terms, ‘subsynaptic reticulum’, ‘endoplasmic reticulum’, ‘microsome’ and ‘vesicular fraction’, were enriched (Fig. 11). ‘Oxygen binding’, ‘iron ion binding’, ‘heme binding’, ‘electron carrier activity’ and ‘oxidoreductase activity’, were enriched in terms of MF (Fig. 12).

Figure 10.

Enriched BP interaction network of hub genes. BP terms, including ‘exogenous drug catabolic process’, ‘drug catabolic process’, ‘secondary metabolic process’ and ‘oxidation reduction’, were significantly enriched. BP, biological process.

Figure 11.

Enriched CC interaction network of hub genes. In CC terms, ‘subsynaptic reticulum’, ‘endoplasmic reticulum’, ‘microsome’ and ‘vesicular fraction’, were enriched. CC, cellular component.

Figure 12.

Enriched MF interaction network of hub genes.‘Oxygen binding’, ‘iron ion binding’, ‘heme binding’, ‘electron carrier activity’ and ‘oxidoreductase activity’, were enriched in terms of MF. MF, molecular function.

Discussion

In the present study, DEGs were analyzed in microarrays to identify hub genes. The top 10 hub genes in the GSE36376 dataset were examined for their prognostic prediction value in HCC. A total of seven genes, including ADH4, CYP2C8, CYP2C9, CYP8B1, SLC22A1, TAT and HSD17B13, were identified as potential prognostic biomarkers. In addition, high expression levels of these hub genes were associated with tumor suppressing roles in HCC. Stratified analysis of clinical factors further revealed their prognostic values in subgroups. A risk score model for patient survival was constructed and evaluated, which confirmed its value for assessing prognosis. A nomogram was created to identify the degree of the contribution made by each factor. Enrichment analysis of the hub genes highlighted the metabolic pathways and biological processes that the hub genes were involved in, which may provide clues into the exact mechanisms of HCC development.

Valuable data in gene microarrays may be lost due to potentially unpredictable problems with the samples when the results of a single piece of research are analyzed (27). Furthermore, using a Student's t-test to analyze microarray data has several limitations (27). Small sample sizes may lead to unreliable variance estimation, leading to a high false-positive rate, while some significant and reliable differences in expression may be missed (28). However, in the present study, analyzing microarray data from a study with a large sample size, allowed the acquisition of potentially useful information for further analysis. In total, 433 samples from the GSE36376 dataset were analyzed to obtain DEGs, in order to determine potential serum biomarkers for HCC diagnosis and prognosis. Focused on an Asian population, the present study also searched for DEGs using the GEO2R online resource. The top 10 hub genes in these DEGs were selected for further analysis, and seven of these hub genes, ADH4, CYP2C8, CYP2C9, CYP8B1, SLC22A1, TAT and HSD17B13, were confirmed to have prognostic value in HCC.

The human ADH4 enzyme is encoded by the ADH4 gene, which maps to 4q22 within the ADH gene cluster (29). Previous studies have revealed that a ADH4 gene variant confers risk for alcohol dependence (AD) and related traits in European Americans and African Americans (29). Edenberg et al (30) reported that 16 single nucleotide polymorphisms (SNP) of ADH4, including rs2226896, are associated with AD, in an independent collaborative study on the genetics of alcoholism. Edenberg et al (31) showed that ADH4 promoter variant −75A/C (rs800759), which could alter ADH4 enzyme expression levels significantly, as well as the 159A/G variant were significantly linked to AD in European Americans and African Americans in a Brazilian population (31). Previous studies have reported that ADH4 may be associated with cluster headaches and personality traits such as agreeableness and extraversion (32,33). In addition, polymorphisms in the ADH4 gene are associated with a decreased risk of ovarian cancer (34) and an increased risk of upper aerodigestive tract cancer (35), which suggests that the ADH4 may be involved in tumorigenesis. Wei et al (36) found that ADH4 mRNA expression in HCC is significantly lower than that in non-cancerous tissue, and ADH4 protein expression is also reduced in HCC, which indicates that ADH4 may serve as a tumor suppressor. These findings are consistent with the results of the present study, where ADH4 was identified as a potential prognostic biomarker for HCC.

The CYP2 family contains many subfamilies, including CYP2A, CYP2B, CYP2C, CYP2D, CYP2E and CYP2F (37). CYP2C8 and CYP2C9 are members of the CYP2C subfamily that are localized in a single gene locus on chromosome 10 (38,39). CYP2C8 shares sequence homology with CYP2C9, that metabolizes several drugs including analgesics (40), antidiabetic and cholesterol-lowering drugs (41). CYP2C9 metabolizes the majority of angiotensin II type 1 receptor blockers (42) and neurological drugs (43). In addition, CYP2C8 has been associated with an increased risk of essential hypertension and coronary artery disease in Bulgarian patients (44), anemia (45), vascular inflammatory disease (46) and breast cancer (47). It has been reported that CYP2C9 downregulation by miR-128-3p is associated with HCC (48). In addition, we have previously reported that analysis of CYP2C8 and CYP2C9 expression in combination is better than analyzing them in isolation (37).

CYP8B1 is predominantly expressed in hepatocytes in a homogenous pattern (49) and is involved in bile acid synthesis of bile acids (50). Overexpression of CYP8B1 alone or in combination with CYP7A1, but not of CYP7A1 alone, reverses obeticholic acid-induced alterations in bile acid levels (taurocholic acid), bile acid composition (taurocholic acid and α/β-muricholic acids) and cholesterol absorption (51). The SNP rs3732860 in the 3′-untranslated region of the CYP8B1 gene is linked to gallstone disease risk in the Chinese Han population (52,53). However, the function of CYP8B1 in cancer remains elusive. The present study indicated that CYP8B1 may serve as a potential biomarker for HCC, and may be involved in tumor initiation and development.

SLC22A1 has the ability to encode solute carrier family 22 member 1, which is not only an uptake transporter but also has a predictive value for the molecular response to imatinib mesylate therapy (54). Patients who experience a major molecular response have higher SLC22A1 expression compared with those without a major molecular response (55). SLC22A1-ABCB1 haplotype profiles can predict imatinib pharmacokinetics in Asian patients with chronic myeloid leukemia (56). Indirect SLC22A1 gene upregulation by dexamethasone may be caused by glucocorticoid receptor-induced hepatocyte nuclear factor-4α expression in primary human hepatocytes, but not in hepatocyte-derived tumor cell lines (57). HCC and cholangiocarcinoma development is typically accompanied by decreased SLC22A1 expression, which may significantly alter the ability of sorafenib to reach active intracellular concentrations in these tumors (58). Downregulated SLC22A1 expression is associated with tumor progression and decreased survival in patients with cholangiocellular carcinoma (59). However, to the best of our knowledge, the relationship between SLC22A1 and HCC patient prognosis has not been reported. The present study demonstrated that SLC22A1 may be a predictive biomarker for HCC.

TAT is associated with the catalyzing the transamination of tyrosine and other aromatic amino acids and plays a role in recovery from tyrosinemia type II, hepatitis and hepatic carcinoma (60). Deficiency of TAT causes marked hypertyrosinemia, which leads to painful palmoplantar hyperkeratosis, pseudodendritic keratitis, and variable mental retardation (61). Recurrent mutation of the TAT gene has been reported in those affected by Richner-Hanhart syndrome (62). Fu et al (63) reported that downregulation of TAT at a frequently deleted region, 16q22, contributes to the pathogenesis of HCC, and it was demonstrated that TAT is a novel tumor suppressor gene. This finding is consistent with the results of the present study, and may be used as an effective serum biomarker for HCC.

Su et al (64) reported that HSD17B13 is upregulated in the livers of patients with non-alcoholic fatty liver disease. HSD17B13 expression is localized to liquid droplets (65). In addition, Chen et al (66) reported that HSD17B13 is downregulated in HCC, and has a tumor suppressor role via inhibition of HCC progression and recurrence (66). This finding is in accordance with the results of the present study, where it was concluded that HSD17B13 may serve as a potential predictive biomarker for HCC.

In regards to metabolic pathways, our previous study demonstrated that the CYP2C subfamily members are involved in chemical carcinogenesis (37). The formation of DNA adducts, dG-C8-IQ, dG-N-IQ, dG-C8-MeIQx and dG-N-MeIQx, may induce liver, colon lung and breast cancer tumorigenesis (37). The present study found that the identified hub genes were also enriched in chemical carcinogenesis. Additionally, gene were enriched in ‘drug metabolism-cytochrome P450’, ‘metabolism of xenobiotics by cytochrome P450’, ‘retinol metabolism’, ‘linoleic acid metabolism’, ‘arachidonic acid metabolism’, ‘tyrosine metabolism’, and ‘steroid hormone biosynthesis’. These metabolic processes and pathways provide evidence that the hub genes may be involved in hepatocarcinogenesis.

There are some limitations to the present study. First, larger samples that include other populations are required to validate the findings of the present study. Second, an increased number of valid clinical factors, such as race, drinking status, smoking status, cirrhosis, Barcelona-Clinic liver cancer staging, hepatitis infection status, antiviral therapy, α-fetoprotein levels and microvascular invasion, should be included in the analysis. Third, functional validation in a well-designed clinical trial is required to examine the biological behavior of prognosis-associated genes on HCC initiation and progression. As the present study explored the potential prognostic biomarkers for HCC, the identification and clinical significance of targeted drugs was not investigated in the present study. Thus, it is crucial to focus future studies on these topics.

The present study indicated that low gene expression of ADH4, CYP2C8, CYP2C9, CYP8B1, SLC22A1, TAT and HSD17B13 are predictors of poor prognosis in HCC. Further functional trials and identification of targeted drugs for these hub genes is warranted to determine clinical application. In detail, trials in vivo and in vitro should be conducted to explore biological behavior, such as invasion, metastasis and proliferation ability. Then, the influence potential drugs on their target genes should be validated, to determine whether targeted overexpression of these genes improve HCC prognosis.