Identification and Validation of Hub Genes with Poor Prognosis in Hepatocellular Carcinoma by Integrated Bioinformatical Analysis

Introduction

Hepatocellular carcinoma (HCC) is the sixth most familiar malignancy and cause of the world’s second largest cancer-related death.¹ It is estimated that there are about 841,000 new cases of liver cancer worldwide in 2018, with 782,000 deaths, causing serious economic burden.² After a great progress has been made in the prevention, diagnosis and treatment of HCC, a significant decrease occurred in the number of viral cases and this was compensated by the progressive expansion of non-viral cases. The improved overall survival was observed due to the wider use of semi-annual surveillance, expanding the proportion of tumors that qualified for curative treatments, and the improved outcome of loco-regional treatments.³ In recent years, more evidence shows that the occurrence and development of HCC is relevant to the abnormal expression of various oncogenes and the inactivation of tumor suppressor genes; these genes have also been increasingly used in the early diagnosis and treatment.^4,5 Thus, it’s crucial to study the target molecules and molecular mechanisms of HCC occurrence and development for formulating more effective diagnosis and treatment strategies.

Gene mutation, cell environment and other factors are bound up with the occurrence, development and metastasis of tumors.⁶ A unique immune response profile of the liver microenvironment where CD4+CD25+ Foxp3 regulatory T-cells (Tregs) also play a crucial role through their immunosuppressive role in HCC development and progression.⁷ During the slow development of HCC, a large number of genomes changed, and accumulated into the phenotypic changes of hepatocytes, resulting in cell intermediates and multiple monoclones, which eventually evolved into HCC.^8,9 High-throughput technology and gene chip have become quicker methods to identify differential expressed genes (DEGs) and functional pathways in the occurrence and development of different diseases over recent years. Hence, bioinformatics has become an essential way to assess gene expression profile. Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) are the two largest public resource gene expression databases in the world. Molecular variations in HCC can be found through descriptive data generated by different microarray platforms.^10–12 Pathway analysis shows that the changes of some important cellular signaling pathways are related to the main pathogenic mechanism.¹³ A few hub genes were classified as pivotal regulators of HCC metastasis through the protein-protein interaction (PPI) network.¹⁴ But owing to the heterogeneity, small sample size and different statistical methods of independent studies, the number of identified functional genes is not nearly enough for verifying HCC pathogenesis, and the most remarkable maladjusted genes in previous studies are not consistent.

Here, in this study, five microarray data sets GSE45267, GSE101685, GSE84402, GSE62232 and GSE45267 downloaded from NCBI-GEO included 165 HCC tissues and 73 normal tissues, and DEGs between them were identified by GEO2R. Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and PPI network analysis were used to identify DEGs, and to evaluate the clinical significance of DEGs in the prognosis of HCC. Eight hub genes were identified as latent molecular targets for exploring new HCC intervention strategies. This work will further explore the development of HCC at molecular level.

Materials and Methods

Microarray Data

The gene expression profiles GSE45267, GSE101685, GSE84402, GSE62232 and GSE45267 were downloaded from GPL570 platform originate from the comprehensive database of gene expression of National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/geo/). According to the labeling information on the platform, the probes were transformed to corresponding gene symbols. After normalizing the gene expression data, log₂ transformation was performed in limma R package (http://www.bioconductor.org/).¹⁵

Recognition of DEGs

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) was employed to define DEGs between HCC and normal tissues. GEO2R is a software for differential analysis of expression microarray based on GEO database, which allows users to compare two or more data sets in GEO series to select different genes expressed under different experimental conditions. FDR and Benjamini-Hochberg control the possibility of large “false positive”. Probe sets without corresponding gene symbols were filtered out, and the average value of log₂FC was used when one gene corresponds to multiple probe sets. LogFC (fold change) >2 and P-value <0.01 were statistically remarkable.

Enrichment of DEGs by GO and KEGG Pathways

To better assess the selected DEGs, GO and KEGG pathway enrichment analysis were performed by David online tool (http://david.ncifcrf.gov) (Version 6.7); a P-value of <0.01 was statistically marked. David is a website for gene annotation, visualization and integrated discovery, which offers a comprehensive set of gene and protein annotation information to extract biological information. GO is the main bioinformatics tool for annotating genes and analyzing their biological functions.¹⁶ The GO mainly describes genes and their products from three aspects: molecular function (MF), biological process (BP) and cellular component (CC). KEGG is a database resource used for the advanced functions and utilities of biological systems from genomic and molecular information.¹⁷

Construction of PPI Network and Analysis of Selected Modules

To understand the relationship between DEGs, PPI network is established by using STRING online database (http://string-db.org) and Module analysis, and a confidence score of ≥0.7 is considered to be remarkable. STRING is a database for searching the interplay between known protein and predicting protein.¹⁸ Furthermore, PPI network is visualized by using Cytoscape (Version 3.4.0) (http://www.cytoscape.org/) (degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max depth = 100). By using the potential information of DEGs in GO and KEGG analysis module, genes with intersection of GO and KEGG ≥10 were selected as hub genes.

Verification of Hub Genes

The expression of hub genes in HCC and their relationship with prognosis were obtained from TCGA database. The expression of hub genes in HCC was analyzed by online database Oncomine, and Kaplan-Meier curve was employed to assess the hub genes and the overall survival and disease free survival rate of HCC patients on the online platform of cBioPortal (http://www.cbioportal.org) and GEPIA2 (http://gepia2.cancer-pku.cn/#index).

Results

Recognition of DEGs in HCC

The gene chips selected this time are GSE45267, GSE101685, GSE84402 and GSE62232 of GPL570 platform. GSE45267 contains 46 HCC and 41 normal samples. GSE101685 consists of 24 HCC and 8 normal samples. GSE84402 includes 14 HCC and 14 normal samples. GSE62232 involves 81 HCC and 10 normal samples. With P<0.01 and [logFC]>2 as the screening criteria, we extracted 308, 459, 463 and 292 DEGs from GSE45267, GSE101685, GSE84402 and GSE62232 respectively, and found 152 overlapping DEGs (Figure 1). Finally, it was confirmed that 39 genes were up-regulated and 113 genes were down-regulated in HCC (Table 1).

Table 1 All 152 Differentially Expressed Genes (DEGs)

Figure 1 Venn diagrams of DEGs in the four datasets. (A) Up-regulated DEGs; (B) down-regulated DEGs.

GO Functional Annotation of DEGs

To explore the biological functions of 152 DEGs, we used online database DAVID 6.7 to do GO analysis, and divided them into three functional groups: molecular function (MF), cellular component (CC) and biological process (BP). According to the survey, the up-regulated genes in HCC were concentrated in BP and CC. BP included mitosis, cell division, G2/M transition of mitotic cell cycle, protein ubiquitination involved in ubiquitin-dependent protein catabolism, complex-dependent catabolism, mitotic spindle tissue and the regulation of spindle microtubules attachment to centromere, etc. CC mainly includes mitotic spindle, cytoplasm, motile, condensed chromosome centromere, nucleus, spindle, spindle microtubules, spindle central region, centrosome, spindle pole and so on. In MF, the genes up-regulated in HCC mainly showed calmodulin-binding protein kinase, binding protein, microtubule-binding region, protein serine/threonine kinase activity, etc. (Table S1). The down-regulated genes in HCC are mainly concentrated in BP, including redox process, cyclooxygenase P450 pathway, exogenous drug catabolism process, cell response to cadmium ion, cell response to zinc ion, negative growth regulation, metabolic process of heterotypic biomass, etc. In CC, there are organelle membrane, extracellular region, blood particles, high density lipoprotein and endoplasmic reticulum membrane. The MF mainly includes oxygen binding, heme binding site, oxidoreductase activity, iron ion binding, monooxygenase activity and so on (Table S2).

KEGG Pathway Enrichment Analysis of DEGs

KEGG is a database resource used for the advanced functions and biological systems of large-scale molecular data sets generated from high-throughput experimental techniques. The up-regulated DEGs are accumulated on cell cycle, oocyte meiosis and p53 signaling pathway, while those down-regulated are rich in retinol metabolism, drug metabolism-cytochrome P450, chemical carcinogenesis, metabolism of cytochrome P450 to exogenous drugs, metabolic pathway, etc. (Table S3).

Construction of PPI Network and Analysis of Selected Modules

A PPI network of 152 DEGs (39 up-regulated and 113 down-regulated genes) was established by STRING. There are 7 up-regulated genes: RRM2, CCNB1, CDK1, AURKA, BUB1B, MAD2L1, TTK, and 25 down-regulated genes: CYP1A2, MT1G, MT1E, CNDP1, CYP26A1, CYP2A6, NAT2, AADAT, CYP2C8, ADH4, RDH16, AKR1D1, CYP8B1, MT1M, CYP2B6, CYP3A4, MT1F, KMO, ADH1C, CYP2C9, ADH1B, IDO2, MT1X, MT1H, TDO2 (Figure 2A). With Module analysis, the genes obtained above are intersected with the genes obtained from GO analysis and KEGG analysis (degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max depth = 100); 4 clusters of genes are obtained (Figure 2B–E, Table S4). At last, 19 hub genes are selected, that is, CCNB1, CDK1, AURKA, BUB1B, MAD2L1, TTK, CYP1A2, MT1E, CYP26A1, CYP2A6, NAT2, CYP2C8, MT1M, CYP2B6, CYP3A4, MT1F, CYP2C9, MT1X and MT1H.

Figure 2 PPI network of DEGs using STRING online database and Module analysis. (A) PPI of all genes. Nodes represent protein, and edges indicate interaction of protein. (B–E) The hub genes were screened from the PPI network in Module analysis.

Expression of Hub Genes in HCC and Its Relationship with Prognosis of Patients

We obtained the expression of hub genes in HCC and their relationship with prognosis from TCGA database. It was found that only 8 genes, CDK1, CYP2C8, CCNB1, AURKA, CYP2C9, BUB1B, MAD2L1 and TTK, were significantly different in HCC expression (P<0.05) (Figure 3), among which the up-regulated genes CDK1, CCNB1, AURKA, BUB1B, MAD2L1, TTK and down-regulated gene CYP2C8, CYP2C9 were interrelated to the overall survival (Figure 4) and disease-free survival (Figure S1) of patients.

Figure 3 Eight significantly expressed genes in TCGA data of HCC. (A) Expression of CDK1 in LIHC based on Sample types. (B) Expression of CCNB1 in LIHC based on Sample types. (C) Expression of AURKA in LIHC based on Sample types. (D) Expression of BUB1B in LIHC based on Sample types. (E) Expression of MAD2L1 in LIHC based on Sample types. (F) Expression of TTK in LIHC based on Sample types. (G) Expression of CYP2C8 in LIHC based on Sample types. (H) Expression of CYP2C9 in LIHC based on Sample types.

Figure 4 The Kaplan-Meier curve of the overall survival between the high-risk and low-risk groups of the eight hub genes. (A) Effect of CDK1 expression level on the survival time of patients with HCC. (B) Effect of CYP2C8 expression level on the survival time of patients with HCC. (C) Effect of CCNB1 expression level on the survival time of patients with HCC. (D) Effect of AURKA expression level on the survival time of patients with HCC. (E) Effect of CYP2C9 expression level on the survival time of patients with HCC. (F) Effect of BUB1B expression level on the survival time of patients with HCC. (G) Effect of MAD2L1 expression level on the survival time of patients with HCC. (H) Effect of TTK expression level on the survival time of patients with HCC.

Discussion

HCC, accounting for 90% of liver cancers, is one of the most familiar and deadly malignancies in those tumors.¹⁹ At present, the molecular mechanism of HCC is not completely clear. Precise diagnosis and prognosis evaluation are still enormous challenges in treatment. Hence, discovering new functional genes is helpful to assess the pathogenesis of HCC and enhance its diagnosis and prognosis. In order to identify new functional genes involved in HCC, we screened the differential genes tied to HCC prognosis through bioinformatics analysis.

Totally, 152 genes differentially expressed in HCC and normal tissues were determined. GO and KEGG enrichment analysis found that 39 up-regulated genes were concentrated on mitosis, cell cycle and oocyte meiosis, while those down-regulated (113) were rich in exogenous drug catabolism, drug metabolism-cytochrome P450, and the metabolism of cytochrome P450 on exogenous drugs. Cell cycle and mitosis crucial in tumor occurrence and development, and are also one of the main targets of tumor drug therapy.^20–22 In the treatment of anti-tumor drugs, the gradual emergence of drug resistance is a crucial reason leading to poor prognosis of patients. Thus, at first, we selected 19 hub genes by PPI network and Module analysis; then, we verified them by TCGA database, and finally discovered 8 hub genes. Those genes, CDK1, CYP2C8, CCNB1, AURKA, CYP2C9, BUB1B, MAD2L1 and TTK, are relevant to the overall survival rate of HCC patients. These genes may be promising biomarkers for HCC prognosis.

Most of these genes participate in the process of cell cycle and cell mitosis, such as CDK1, CCNB1, BUB1B, AURKA, MAD2L1 and TTK. CDK1 can form a complex with cyclin B1 (CCNB1) and cyclin B2 (CCNB2), regulate G2/M phase of mammalian cell cycle, and exert important effects in mitosis.²³ Recent studies have suggested that CDK1 suppresses the proliferation, migration and invasion of HCC cells.^24,25 AURKA, also described as Aurora kinase A, is a vital serine/threonine kinase responsible for regulating cell mitosis, which has a major part to play in centrosome replication and separation, spindle assembly, maturation, chromosome arrangement, spindle assembly checkpoint and cytoplasm division.²⁶ The increased expression of AURKA may bring about chromosome instability, transformation and centrosome amplification in mammalian cells, and promote the carcinogenesis of c-myc.^27,28 AURKA is over-expressed in many tumors, especially in liver cancer.^29,30 BUB1B is an integral part to spindle assembly checkpoint, including MAD1, MAD2, MAD3/Bub1b, Mps1, Bub1 and Bub3.³¹ In one study, Qiu et al³² verified the carcinogenicity of BUB1B in HCC, including promoting the progression of G0/G1 cell cycle, inhibiting the apoptosis of HCC cells and causing poor prognosis. MAD2L1 is effective in spindle mitosis, and the imbalance of MAD2L1 brings about chromosome instability and chromosome aneuploidy.³³ At the moment, many evidences demonstrated that MAD2L1 can be employed as a biomarker of poor prognosis of lung cancer,^34,35 and MAD2L1 gene variation can result in the decrease of shuttle checkpoint function and increase the risk of lung cancer.³⁶ The function of MAD2L1 in HCC needs to be verified by more clinical samples. In a basic experiment, miR-200c-5p suppressed HCC proliferation and metastasis by inhibiting MAD2L1.³⁷ The human monopolar spindle 1 (hMps1/TTK) gene (NM_003318) locates on chromosome 6q13-q21 and encodes a bis-serine/threonine and tyrosine protein kinase.³⁸ TTK is necessary for mitotic checkpoint and wrong chromosome attachment.³⁹ Elevated TTK level will result in centrosome amplification, over-activation of spindle assembly checkpoint and chromosome instability, thus promoting tumor occurrence.⁴⁰ According to Liang et al, TTK gene is a potential therapeutic target for anti-sorafenib resistance in HCC patients.⁴¹ CYP2C8 and CYP2C9 are members of CYP450 gene family, both of which are down-regulated in HCC tissues; the low level of CYP2C8 has been verified to be related to poor overall survival rate and disease-free survival rate.⁴² Overexpression of CYP2C8 suppressed the proliferation, clonality, migration, invasion and cell cycle of HCC cells via PI3K/Akt/p27 ^kip1 Axis. Members of CYP2C subfamily participate in the metabolism of many endogenous and exogenous substances, including drug metabolism. It may be the reason that low expression of CYP2C8 and CYP2C9 were associated with poor survival rate of HCC in part.⁴³

All these theories are in conformity with our results. These genes may play a central role in the occurrence, development and drug resistance of HCC. There are some limitations. One of the main problems is that we did not collect clinical samples for verification. What’s more, we only analyzed the overall survival rate of patients, but did not evaluate the disease-free survival rate, clinical characteristics and other indicators. Others, after the functional enrichment analysis of DEGs, we ignored that in signaling pathway, so the results could be further explored.

Conclusion

To summarize, we identified 152 DEGs with differential expression in HCC through bioinformatics analysis, and verified them with TCGA database, and identified 8 hub genes associated with the prognosis of patients. These hub genes are accumulated on cell cycle and drug catabolism, which may be latent targets for HCC, and supply further strategies for clinical diagnosis, accurate treatment and prognosis analysis.