OL Oncology Letters 1792-1074 1792-1082 D.A. Spandidos 10.3892/ol.2020.11752 OL-0-0-11752 Articles Identification of hub genes involved in the occurrence and development of hepatocellular carcinoma via bioinformatics analysis MiNingning 1 2 * CaoJie 1 2 3 4 5 6 * ZhangJinduo 2 3 4 5 FuWenkang 1 2 3 4 5 HuangChongfei 1 2 3 4 5 GaoLong 1 2 3 4 5 YuePing 2 3 4 5 BaiBing 2 3 4 5 LinYanyan 1 2 3 4 5 MengWenbo 1 2 3 4 5 LiXun 4 5 7 The First Clinical Medical School, Lanzhou University, Lanzhou, Gansu 730000, P.R. China Department of Special Minimally Invasive Surgery, The First Hospital of Lanzhou University, Lanzhou, Gansu 730000, P.R. China Institute of Genetics, School of Basic Medical Sciences, Lanzhou University, Lanzhou, Gansu 730000, P.R. China Gansu Province Institute of Hepatopancreatobiliary, Lanzhou, Gansu 730000, P.R. China Gansu Province Key Laboratory Biotherapy and Regenerative Medicine, Lanzhou, Gansu 730000, P.R. China Laboratory Department, The First Hospital of Lanzhou University, Lanzhou, Gansu 730000, P.R. China The Fifth Department of General Surgery, The First Hospital of Lanzhou University, Lanzhou, Gansu 730000, P.R. China Correspondence to: Dr Yanyan Lin or Professor Wenbo Meng, Department of Special Minimally Invasive Surgery, The First Hospital of Lanzhou University, 1 Donggang West Road, Lanzhou, Gansu 730000, P.R. China, E-mail: ldyy_linyy@lzu.edu.cn, E-mail: mengwb@lzu.edu.cn

Contributed equally

 08 2020 17 06 2020 20 2 1695 1708 28092019 07052020 Copyright: © Mi et al. 2020 This is an open access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

Hepatocellular carcinoma (HCC) is a heterogeneous malignancy, which is a major cause of cancer morbidity and mortality worldwide. Thus, the aim of the present study was to identify the hub genes and underlying pathways of HCC via bioinformatics analyses. The present study screened three datasets, including GSE112790, GSE84402 and GSE74656 from the Gene Expression Omnibus (GEO) database, and downloaded the RNA-sequencing of HCC from The Cancer Genome Atlas (TCGA) database. The differentially expressed genes (DEGs) in both the GEO and TCGA datasets were filtered, and the screened DEGs were subsequently analyzed for functional enrichment pathways. A protein-protein interaction (PPI) network was constructed, and hub genes were further screened to create the Kaplan-Meier curve using cBioPortal. The expression levels of hub genes were then validated in different datasets using the Oncomine database. In addition, associations between expression and tumor grade, hepatitis virus infection status, satellites and vascular invasion were assessed. A total of 126 DEGs were identified, containing 70 upregulated genes and 56 downregulated genes from the GEO and TCGA databases. By constructing the PPI network, the present study identified hub genes, including cyclin B1 (CCNB1), cell-division cycle protein 20 (CDC20), cyclin-dependent kinase 1, BUB1 mitotic checkpoint serine/threonine kinase β (BUB1B), cyclin A2, nucleolar and spindle associated protein 1, ubiquitin-conjugating enzyme E2 C (UBE2C) and ZW10 interactor. Furthermore, upregulated CCNB1, CDC20, BUB1B and UBE2C expression levels indicated worse disease-free and overall survival. Moreover, a meta-analysis of tumor and healthy tissues in the Oncomine database demonstrated that BUB1B and UBE2C were highly expressed in HCC. The present study also analyzed the data of HCC in TCGA database using univariate and multivariate Cox analyses, and demonstrated that BUB1B and UBE2C may be used as independent prognostic factors. In conclusion, the present study identified several genes and the signaling pathways that were associated with tumorigenesis using bioinformatics analyses, which could be potential targets for the diagnosis and treatment of HCC.

 hepatocellular carcinoma hub genes prognosis bioinformatics microarray Introduction

Hepatocellular carcinoma (HCC) is one of the few cancer types associated with a rising incidence in recent years, with the incidence in Mongolia being as high as 93.7 per 100,000 individuals (1). Certain aspects of the etiology and development of HCC have been reported, and it is speculated that the complex pathogenesis of HCC is associated with cirrhosis, viral hepatitis, aflatoxins, specific chemical carcinogens and abnormal regulation of hormones (2,3). Among these, chronic liver disease is the most important factor leading to the occurrence of liver cancer (4). The infection of hepatitis B virus (HBV) accounts for 60% of HCC cases in developing countries and for 23% of cases in developed countries, while infection of hepatitis C virus (HCV) accounts for 23% of HCC cases in developing countries and 20% of cases in developed countries (5). Moreover, sub-Saharan Africa, East Asia and Southeast Asia are the regions with the highest incidence of HBV, while USA, Europe and Japan have high incidences of HCV (4). The incidence of non-alcoholic fatty liver disease (NAFLD) is increasing with the rise in obesity and other metabolic syndromes, which seriously affects the health of the affected individual (6). According to data generated from a cross-sectional study, 25–30% of people with a Western lifestyle have a fat content in their liver that exceeds the normal level; of these, 2–5% suffer from NAFLD, and 1–2% of the affected adults may be at risk of developing non-alcoholic steatohepatitis cirrhosis (7).

The incidence of HCC is increasing worldwide; however, early diagnoses and treatment of HCC remain an issue (4). Due to the low-level medical and healthcare in developing countries, the incidence of HCC is rising, with an estimated global incidence of liver cancer per 100,000 individuals being 9.3 in 2018 (1), and the prognosis is poor (8). For patients without liver cirrhosis and distant metastasis of the tumor, hepatic resection is currently the first treatment choice and the most effective way to treat HCC (4). However, for patients with liver cirrhosis, surgical resection is a contraindication, and liver transplantation should be performed to achieve a favorable prognosis (4). The size of the tumor is not the definite limiting factor for surgery; the presence or absence of distant metastasis and vascular invasion are the main factors considered for surgical resection (9). In recent years, laparoscopic hepatectomy has become a major surgical method of radical resection of tumors (10). Furthermore, compared with traditional open surgery, laparoscopic surgery is minimally invasive, which has the advantages of small trauma, fewer perioperative complications and quick postoperative recovery, while having a similar prognosis to open surgery (10). For patients without indications for surgical resection, minimally invasive treatment can be applied according to a variety of factors such as tumor characteristics, complications, tumor site and performance status (4). Local treatment methods include radiofrequency ablation and microwave ablation, while other therapeutic methods are locoregional therapies such as chemoembolization and radioembolization (11). Moreover, systemic chemotherapy is a feasible method to improve the survival rate of patients with advanced liver cancer (3,4). Sorafenib is a first-line drug for the treatment of advanced HCC (12). With the development of sorafenib derivatives, additional satisfactory antitumor drugs are expected to be developed (12). In addition, gene chip technology for the study of liver cancer can accurately depict the molecular expression profile of HCC and identify the specific genes associated with HCC (13–15).

Microarray technology, which is a highly efficient and accurate transcriptional expression technology, has been successfully applied in the screening of molecular markers of almost all human malignant tumors, particularly liver cancer. The present study obtained GSE112790, GSE84402 and GSE74656 microarray datasets from the Gene Expression Omnibus (GEO) database. In addition, RNA-sequence (seq) data were downloaded from The Cancer Genome Atlas (TCGA) database and further analyzed via bioinformatics methods. These analyses resulted in the identification of a series of hub genes that may be associated with HCC, which might be used as molecular markers in the early diagnosis and treatment targets of HCC.

 Materials and methods <sec> <title>Data source

The present study downloaded three datasets, namely GSE112790 (16), GSE84402 (17) and GSE74656 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74656), from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), and the following conditions were met for each set of samples: i) Samples were obtained from human liver tissue; ii) chips contain both healthy samples and tumor tissue of HCC; iii) the total number of samples ≥10.

GSE112790 and GSE84402 datasets were based on GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array; Affymetrix; Thermo Fisher Scientific, Inc.; http://www.affymetrix.com/index.affx), while GSE74656 was based on GPL16043 (GeneChip® PrimeView™ Human Gene Expression Array with external spike-in RNAs; Affymetrix; Thermo Fisher Scientific, Inc.; http://www.affymetrix.com/index.affx). The GSE112790 dataset contained 15 healthy liver tissues and 183 tumor samples, while the GSE84402 dataset contained 14 healthy liver samples and 14 malignant tissues. Moreover, the GSE74656 dataset contained five healthy liver and five malignant liver tissues. The RNA-seq of HCC consists of 374 tumor liver samples and 50 healthy tissues that were downloaded from TCGA database (https://cancergenome.nih.gov/).

 Screening for differentially expressed genes (DEGs)

The raw data of the gene expression profiles were analyzed using R software v3.5.2 (https://www.r-project.org) and the Bioconductor package ‘Limma’ v3.36.5 (https://bioconductor.org/packages/limma/) (18). The robust multi-array average method (19) was used to complete background correction and normalization of all original microarray data. Subsequently, the present study used the combat function of the ‘sva’ package v3.30.1 (https://bioconductor.org/packages/sva/) (20) of R to remove any batch effect. DEGs were screened using the ‘Limma’ package, and the Bioconductor package ‘edgeR’ v3.24.3 (21) was used to analyze and process the data downloaded from TCGA database. The screening criteria for DEGs were both adjusted to P<0.01 and |log fold-change| >2.0.

 Functional enrichment analysis of DEGs

In order to perform Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp) pathway and Gene Ontology (GO; http://geneontology.org) function analyses, DEGs were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov/). Biological annotations were considered to be significantly enriched when P<0.05. DAVID bioinformatics resources include an integrated biological knowledge base and analytical tools designed to systematically extract biological significance from large gene lists (22). The results were visualized using the GOChord and GOCircle functions of ‘GOplot’ package v1.0.2 (https://CRAN.R-project.org/package=GOplot) (23).

 Construction of protein-protein interaction (PPI) network and analysis of modules

The Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) online software (24) and Cytoscape tools were utilized to construct a PPI regulatory network of DEGs. Cytoscape 3.6.0 (https://cytoscape.org/) is a free software that graphically displays PPI networks, and performs analysis and editing. Moreover, Cytoscape is an application that runs on a personal computer rather than on a web browser (25). Medium confidence (0.400) was implied in the STRING tool. The present study utilized the plug-in Molecular Complex Detection (MCODE 1.5.1) (26), an application of Cytoscape, to mine the core modules of PPI networks that were significantly associated with liver cancer, and further investigate the biological functions of the module. The criteria for screening are as follows: Degree cut-off=2; node score cut-off=0.2; maximum depth=100; and k-score=2. Hub proteins were selected based on their associations with other proteins, which were sorted by the degree value in the network.

 Selection and analysis of hub genes

The numbers of nodes and edges in the PPI network were counted in order to screen for hub genes. The present study used the Kaplan-Meier curve to determinate the overall and disease-free survival of the hub genes in the Liver Hepatocellular Carcinoma (TCGA, Firehose Legacy; http://www.cbioportal.org/study/summary?id=lihc_tcga) database of the cBioPortal for cancer genomics (v3.0.2; http://www.cbioportal.org/) (27,28). The study downloaded the original data of hub genes and used the Renyi test to determine if there is a significant difference in survival rate between two groups when their survival curves cross each other using SAS software v9.4 (SAS Institute, Inc.) (29).

For the purpose of assessing the expression levels of hub genes in different datasets, these were submitted to the Oncomine (https://www.oncomine.org/) database (30), which was also used to analyze the association between gene expression and tumor grade, hepatitis virus infection status, satellites and vascular invasion by analyzing the Wurmbach Liver dataset (31). The clinical data of HCC from TCGA database was then obtained, and the ‘survival’ package (v3.1–8; http://CRAN.R-project.org/package=survival) of R software was used to perform a univariate analysis and construct the multivariable Cox regression models on BUB1B and UBE2C, and other clinicopathological parameters, and calculated hazard ratios and corresponding 95% confidence intervals.

 Gene set enrichment analysis (GSEA)

In order to investigate the biological functions and pathways of different expression levels of core genes, GSEA was performed using GSEA software (v4.0.2; http://software.broadinstitute.org/gsea/index.jsp) (32) with gene set c2 (cp.kegg.v7.0.symbols.gmt) from the Molecular Signatures Database (MSigDB; http://www.gsea-msigdb.org/gsea/msigdb). The MSigDB of c2 is a pathway gene set, which was curated from published canonical pathways and experimental signatures (33). P<0.05 and false discovery rate (FDR) <0.25 were considered to indicate a statistically significant difference.

 Results <sec> <title>Screening for DEGs

In total, 126 genes that met the differential threshold were screened by bioinformatics methods. Overall, 83 genes were upregulated and 93 genes were downregulated in the GEO datasets, which contained a total of 34 healthy samples and 202 cancer samples. In total, 374 tumors and 50 healthy liver tissues were analyzed in TCGA database using the same bioinformatics methods, and 3,297 upregulated and 300 downregulated genes were identified. Hub genes were identified according to the adjusted P<0.01 and a |log fold-change| of >2.0.

The heatmaps and volcano plots of the DEGs obtained from the GEO and TCGA databases are shown in Figs. 1 and 2. Heatmaps were created to visualize the expression values of DEGs in the different samples, while volcano plots were plotted to present the distribution of each gene based on the criteria of an adjusted P<0.01 and a |log fold-change| of >2.0. Via comprehensive bioinformatics analysis of these DEGs, 126 DEGs that existed in both the GEO and TCGA databases were identified, including 70 upregulated and 56 downregulated genes (Table I).

 Functional enrichment analyses of DEGS

The present study used DAVID and R software to perform functional and pathway enrichment analyses for identifying the biological classification of DEGs. The following conclusions were drawn: i) The DEGs were enriched in 83 GO terms; ii) the biological processes (BP) of DEGs were primarily enriched in ‘negative regulation of growth’, ‘cell division’, ‘cell proliferation’, ‘cell migration’ and ‘cell cycle’; iii) molecular function (MF) was mainly enriched in ‘oxygen binding’, ‘monooxygenase activity’ and ‘protein binding’; and iv) the ‘midbody’, ‘cytosol’ and ‘nucleus’ were primarily enriched in the cell components (CC) function (Table SI). The top ten GO terms are presented in (Fig. 3A). The results of the KEGG pathway enrichment analysis indicated that DEGs were primarily enriched in 18 pathways, with the top six being ‘p53 signaling pathway’ (hsa04115), ‘cell cycle’ (hsa04110), ‘cellular senescence’ (hsa04218), ‘5’ AMP-activated protein kinase (AMPK) signaling pathway’ (hsa04152), ‘peroxisome proliferator-activated receptors (PPARs) signaling pathway’ (hsa03320) and ‘chemical carcinogenesis’ (hsa05204) (Fig. 3B).

 Construction of PPI network and analysis of modules

The PPI network of hub genes, which was composed of 102 nodes and 1,023 edges, was constructed using the STRING website and visualized by Cytoscape (Fig. 4A). The most essential module was obtained using MCODE, a Cytoscape plugin, and is presented in (Fig. 4B). DAVID was used to identify GO items of modules that were significant (P<0.05). Moreover, it was demonstrated that ‘cell division’, ‘cell cycle’, ‘cell division’, ‘DNA repair’ and ‘regulation of cell cycle’ were the five most significantly enriched BP of the module (data not shown).

 Selection and analysis of hub genes

The top eight genes that had >44 nodes in the PPI network were selected as hub genes. The names, coded proteins and functions of these hub genes are presented in Table II. A Kaplan-Meier curve for these eight hub genes was constructed in cBioPortal, and it was demonstrated that patients with alterations in cyclin B1 (CCNB1), cell-division cycle protein 20 (CDC20), BUB1 mitotic checkpoint serine/threonine kinase β (BUB1B) and ubiquitin-conjugating enzyme E2C (UBE2C) exhibited worse disease-free survival (Fig. 5A) and overall survival rate (Fig. 5B). Subsequently, the results of the expression of these four genes were validated in the Oncomine database, and it was indicated that BUB1B and UBE2C were significantly upregulated in HCC in different datasets (Fig. 6). By analyzing the Wurmbach liver dataset, it was speculated that high expression of these two genes is associated with tumor grade, hepatitis virus infection status, satellites and vascular invasion (Fig. 7).

After deleting the incomplete sample information of the clinical data, a univariate analysis was performed. The results suggested that high expression levels of BUB1B and UBE2C, tumor stage, T classification and metastasis predicted poorer survival, and multivariate Cox analyses demonstrated that BUB1B and UBE2C were independent prognostic factors for HCC (Tables III and IV). Collectively, the results suggested that the alteration of BUB1B and UBE2C may be a potential target for the diagnosis and treatment of HCC.

 GSEA

The GSEA results indicated that the high-expression phenotypes of BUB1B and UBE2C were significantly associated with tumor-associated pathways (P<0.05; FDR<0.25; Fig. 8), including ‘p53 signaling pathway’, ‘cell cycle’, ‘WNT signaling pathway’, ‘T-cell receptor signaling pathway’, ‘vascular endothelial growth factor (VEGF) signaling pathway’ and other pathways associated with malignant tumors.

 Discussion

Recently, with the improvement of genome sequencing, biochips and high-throughput sequencing technology, an increasing number of studies have used bioinformatics methods for chip datasets analysis, which provides a novel and effective way to identify promising markers for HCC diagnosis and treatment (13–15,34,35). In contrast to previous studies, in the present study three datasets (GSE112790, GSE84402 and GSE74656) were downloaded from GEO and RNA-seq data from TCGA, and were combined to determine the DEGs of tumor and healthy tissues via comprehensive bioinformatics analysis methods. A total of 126 hub genes were selected, including 56 downregulated and 70 upregulated genes. Subsequently, the present study performed GO and KEGG pathway analyses for these DEGs. It was found that BPs were mainly involved in ‘negative regulation of growth’, ‘cell division’, ‘cell proliferation’, ‘cell migration’ and ‘cell cycle,’ while CC primarily contained the ‘midbody’, ‘cytosol’ and ‘nucleus’. Moreover, MF of DEGs were primarily enriched in ‘oxygen binding’, ‘monooxygenase activity’ and ‘protein binding’. In addition, KEGG were mainly enriched in ‘p53 signaling pathway’, ‘cell cycle’, ‘cellular senescence’ and ‘AMPK signaling pathway’.

Among the 126 DEGs that were identified, a total of eight hub genes were screened, including CCNB1, CDC20, cyclin-dependent kinase 1 (CDK1), BUB1B, cyclin A2 (CCNA2), nucleolar and spindle-associated protein 1 (NUSAP1), UBE2C and ZW10 interactor (ZWINT). A survival analysis was performed on these hub genes, and it was demonstrated that alterations in CCNB1, CDC20, BUB1B and UBE2C indicated poor outcomes, while alterations in CDK1, CCNA2, NUSAP1 and ZWINT did not result in meaningful differences (data not shown). In addition, the present study compared four different groups of experimental results in Oncomine, and analysis of tumoral and healthy tissue identified that BUB1B and UBE2C were significantly upregulated and highly associated with tumor grade, hepatitis virus infection status, satellites and vascular invasion, thus suggesting that these selected genes may play a critical role in the initiation of HCC. Furthermore, univariate and multivariate Cox analysis results suggested that BUB1B and UBE2C were independent prognostic factors for HCC. In addition, GSEA analysis identified the pathways via which these two genes may function, including t ‘p53 signaling pathway’, ‘cell cycle’, ‘WNT signaling pathway’, ‘T-cell receptor’, ‘VEGF signaling pathway’ and other pathways associated with malignant tumors.

In the present study, CDC20 was primarily involved in ‘mitotic nuclear division’ and ‘cell division’ of BPs and the ‘spindle’ and ‘cytoplasm’ (data not shown) of the cell. Zhang et al (36), reported that high expression of CDC20 in prostate cancer often indicates a poor prognosis. Moreover, a previous study revealed that CDC20 is associated with tumor metastasis (37); thus, future research should focus on CDC20 target inhibitors to develop new drugs to inhibit tumor metastasis.

As a significant checkpoint of cell mitosis, BUB1B is involved in the normal process of cell mitosis (38), and upregulation of BUB1B is observed in various human malignancies, which are associated with the development and progression of the tumor, as well as aggressive biological features (39–42). However, low expression of BUB1B in colon adenocarcinoma (43) and lung tumor (44) causes metastasis and poorer survival.

UBE2C is a protein in the UBE2 family, which is involved in ubiquitination during protein modification, and degrades anaphase promoting complex/cyclosome substrates via proteasomes and promotes mitotic exit (45). Zhang et al (46), revealed that UBE2C is highly expressed in gastric cancer, which could lead to the occurrence and development of malignant tumors, and this upregulation is significantly associated with lymph node metastasis and tumor stage. In addition, high expression of UBE2C is also closely associated with breast cancer, colon cancer and other malignant tumors (47–49).

ZWINT is a member of the MIS12 kinetochore complex component, which is necessary for spindle checkpoint activity. It also participates in biological processes such as cell cycle and cell division (50), and may be a promising target for tumor therapy. Previous studies have reported that CCNA2 expression is elevated in human breast cancer and pancreatic ductal adenocarcinoma (35,51). Moreover, NUSAP1 may be involved in the regulation of chromosome alignment (52), and several studies have revealed that its expression in high levels is associated with prostate, colon and liver cancer (34,53,54). In addition, CDK1 may be an underlying therapeutic target for the treatment of cancer, as it is a principal regulatory agent for the cell cycle, and functions by controlling the centrosome cycle and the initiation of mitosis (55). The protein CCNB1 plays a pivotal role in tumor aggressiveness by modulating the cell cycle at the G2/M (mitosis) transition (56). Zhuang et al (57) also indicated that upregulation of CCNB1 may predict poor prognosis in patients with HCC.

In conclusion, the present study performed preliminary examinations to identify the mechanism underlying the initiation of cancer formation, as well as the development and progression of HCC. Comprehensive bioinformatics analyses were performed in order to identify 126 DEGs and eight hub genes. Furthermore, the results demonstrated that BUB1B and UBE2C may be potential targets for the diagnosis and treatment of HCC. However, the main limitation of the present study is the lack of quantitative PCR analysis to verify the expression levels of BUB1B and UBE2C. Collectively, the present study provides a new direction for further studies on HCC. The DEGs identified in HCC tissues may be involved in carcinogenesis and progression. Particularly, BUB1B and UBE2C, may serve as potential candidate biomarkers and therapeutic targets for HCC. Future studies should focus on the verification of the present findings in clinical practice and on identifying the molecular mechanisms involved in HCC.

 Supplementary Material Supporting Data Acknowledgements

Not applicable.

 Funding

This study was funded by the National Natural Science Foundation of China (grant no. 81872036) and The First Hospital of Lanzhou University Fund 2017 (grant no. ldyyyn2017-04).

 Availability of data and materials

The datasets analyzed during the present study are available in the GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA repositories (https://cancergenome.nih.gov/).

 Authors' contributions

NM and JC designed the study and wrote the initial draft of the manuscript. JZ, WF, CH and LG contributed to the analysis and interpretation of data. PY and BB helped to design the study and acquire data, YL, WM and XL participated in interpreting the results and revised the writing of manuscript. All authors have reviewed and approved the final manuscript.

 Ethics approval and consent to participate

Not applicable.

 Patient consent for publication

Not applicable.

 Competing interests

The authors declare that they have no competing interests.

 References McGlynnKAPetrickJLEl-SeragHBEpidemiology of hepatocellular carcinomaHepatologyApr222020(Epub ahead of print)10.1002/hep.3128832319693 BrayFFerlayJSoerjomataramISiegelRLTorreLAJemalAGlobal cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countriesCA Cancer J Clin68394424201810.3322/caac.2149230207593 MaluccioMCoveyARecent progress in understanding, diagnosing, and treating hepatocellular carcinomaCA Cancer J Clin62394399201210.3322/caac.2116123070690 FornerAReigMBruixJHepatocellular carcinomaLancet39113011314201810.1016/S0140-6736(18)30010-229307467 ParkinDMThe global health burden of infection-associated cancers in the year 2002Int J Cancer11830303044200610.1002/ijc.2173116404738 RinellaMENonalcoholic fatty liver disease: A systematic reviewJAMA31322632273201510.1001/jama.2015.537026057287 MirizziAFrancoILeoneCMBonfiglioCCozzolongoRNotarnicolaMGiannuzziVTutinoVDe NunzioVBrunoIEffects of some food components on non-alcoholic fatty liver disease severity: Results from a cross-sectional studyNutrients112744201910.3390/nu11112744 DimitroulisDDamaskosCValsamiSDavakisSGarmpisNSpartalisEAthanasiouAMorisDSakellariouSKykalosSFrom diagnosis to treatment of hepatocellular carcinoma: An epidemic problem for both developed and developing worldWorld J Gastroentero2352825294201710.3748/wjg.v23.i29.5282 FuksDDokmakSParadisVDioufMDurandFBelghitiJBenefit of initial resection of hepatocellular carcinoma followed by transplantation in case of recurrence: An intention-to-treat analysisHepatology55132140201210.1002/hep.2468021932387 HanHSShehtaAAhnSYoonYSChoJYChoiYLaparoscopic versus open liver resection for hepatocellular carcinoma: Case-matched study with propensity score matchingJ Hepatol63643650201510.1016/j.jhep.2015.04.00525872167 HartkeJJohnsonMGhabrilMThe diagnosis and treatment of hepatocellular carcinomaSemin Diagn Pathol34153159201710.1053/j.semdp.2016.12.01128108047 ZhangPYangYWenFHeXTangRDuZZhouJZhangJLiQCost-effectiveness of sorafenib as a first-line treatment for advanced hepatocellular carcinomaEur J Gastroenterol Hepatol27853859201510.1097/MEG.000000000000037325919775 LiLLeiQZhangSKongLQinBScreening and identification of key biomarkers in hepatocellular carcinoma: Evidence from bioinformatic analysisOncol Rep3826072618201710.3892/or.2017.594628901457 ZhangRYeJHuangHDuXMining featured biomarkers associated with vascular invasion in HCC by bioinformatics analysis with TCGA RNA sequencing dataBiomed Pharmacother118109274201910.1016/j.biopha.2019.10927431545220 ZhangCPengLZhangYLiuZLiWChenSLiGThe identification of key genes and pathways in hepatocellular carcinoma by bioinformatics analysis of high-throughput dataMed Oncol34101201710.1007/s12032-017-0963-928432618 ShimadaSMogushiKAkiyamaYFuruyamaTWatanabeSOguraTOgawaKOnoHMitsunoriYBanDComprehensive molecular and immunological characterization of hepatocellular carcinomaEBioMedicine40457470201910.1016/j.ebiom.2018.12.05830598371 WangHHuoXYangXRHeJChengLWangNDengXJinHWangNWangCSTAT3-mediated upregulation of lncRNA HOXD-AS1 as a ceRNA facilitates liver cancer metastasis by regulating SOX4Mol Cancer16136201710.1186/s12943-017-0680-128810927 RitchieMEPhipsonBWuDHuYLawCWShiWSmythGKLimma powers differential expression analyses for RNA-sequencing and microarray studiesNucleic Acids Res43e47201510.1093/nar/gkv00725605792 BolstadBMIrizarryRAAstrandMSpeedTPA comparison of normalization methods for high density oligonucleotide array data based on variance and biasBioinformatics19185193200310.1093/bioinformatics/19.2.18512538238 LeekJTJohnsonWEParkerHSJaffeAEStoreyJDThe sva package for removing batch effects and other unwanted variation in high-throughput experimentsBioinformatics28882883201210.1093/bioinformatics/bts03422257669 RobinsonMDMcCarthyDJSmythGKEdgeR: A Bioconductor package for differential expression analysis of digital gene expression dataBioinformatics26139140201010.1093/bioinformatics/btp61619910308 Huang daWShermanBTLempickiRASystematic and integrative analysis of large gene lists using DAVID bioinformatics resourcesNat Protoc44457200910.1038/nprot.2008.21119131956 WalterWSánchez-CaboFRicoteMGOplot: An R package for visually combining expression data with functional analysisBioinformatics3129122914201510.1093/bioinformatics/btv30025964631 SzklarczykDGableALLyonDJungeAWyderSHuerta-CepasJSimonovicMDonchevaNTMorrisJHBorkPSTRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasetsNucleic Acids Res47D607D613201910.1093/nar/gky113130476243 DemchakBHullTReichMLiefeldTSmootMIdekerTMesirovJPCytoscape: The network visualization tool for GenomeSpace workflowsF1000Res3151201410.12688/f1000research.4492.225165537 BaderGDHogueCWAn automated method for finding molecular complexes in large protein interaction networksBMC Bioinformatics42200310.1186/1471-2105-4-212525261 GaoJAksoyBADogrusozUDresdnerGGrossBSumerSOSunYJacobsenASinhaRLarssonEIntegrative analysis of complex cancer genomics and clinical profiles using the cBioPortalSci Signal6pl1201310.1126/scisignal.200408823550210 CeramiEGaoJDogrusozUGrossBESumerSOAksoyBAJacobsenAByrneCJHeuerMLLarssonEThe cBio cancer genomics portal: An open platform for exploring multidimensional cancer genomics dataCancer Discov2401404201210.1158/2159-8290.CD-12-009522588877 LiHHanDHouYChenHChenZStatistical inference methods for two crossing survival curves: A comparison of methodsPLoS One10e0116774201510.1371/journal.pone.011677425615624 RhodesDRKalyana-SundaramSMahavisnoVVaramballyRYuJBriggsBBBarretteTRAnstetMJKincead-BealCKulkarniPOncomine 3.0: Genes, pathways, and networks in a collection of 18,000 cancer gene expression profilesNeoplasia9166180200710.1593/neo.0711217356713 WurmbachEChenYBKhitrovGZhangWRoayaieSSchwartzMFielIThungSMazzaferroVBruixJGenome-wide molecular profiles of HCV-induced dysplasia and hepatocellular carcinomaHepatology45938947200710.1002/hep.2162217393520 SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMAPaulovichAPomeroySLGolubTRLanderESMesirovJPGene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profilesProc Natl Acad Sci USA1021554515550200510.1073/pnas.050658010216199517 LiberzonABirgerCThorvaldsdóttirHGhandiMMesirovJPTamayoPThe Molecular Signatures Database (MSigDB) hallmark gene set collectionCell Syst1417425201510.1016/j.cels.2015.12.00426771021 ZhouLDuYKongLZhangXChenQIdentification of molecular target genes and key pathways in hepatocellular carcinoma by bioinformatics analysisOnco Targets Ther1118611869201810.2147/OTT.S15673729670361 DongSHuangFZhangHChenQOverexpression of BUB1B, CCNA2, CDC20 and CDK1 in tumor tissues predicts poor survival in pancreatic ductal adenocarcinomaBiosci Rep39BSR20182306201910.1042/BSR2018230630765611 ZhangQHuangHLiuALiJLiuCSunBChenLGaoYXuDSuCCell division Cycle 20 (CDC20) drives prostate cancer progression via stabilization of β-catenin in cancer stem-like cellsEBioMedicine42397407201910.1016/j.ebiom.2019.03.03230904606 ChengSCastilloVSlivaDCDC20 associated with cancer metastasis and novel mushroom-derived CDC20 inhibitors with antimetastatic activityInt J Oncol5422502256201931081056 ShinHJBaekKHJeonAHParkMTLeeSJKangCMLeeHSYooSHChungDHSungYCDual roles of human BubR1, a mitotic checkpoint kinase, in the monitoring of chromosomal instabilityCancer Cell4483497200310.1016/S1535-6108(03)00302-714706340 AndoKKakejiYKitaoHIimoriMZhaoYYoshidaROkiEYoshinagaKMatumotoTMoritaMHigh expression of BUBR1 is one of the factors for inducing DNA aneuploidy and progression in gastric cancerCancer Sci101639645201010.1111/j.1349-7006.2009.01457.x20132214 YamamotoYMatsuyamaHChochiYOkudaMKawauchiSInoueRFuruyaTOgaANaitoKSasakiKOverexpression of BUBR1 is associated with chromosomal instability in bladder cancerCancer Genet Cytogen1744247200710.1016/j.cancergencyto.2006.11.012 LiuAWCaiJZhaoXLXuAMFuHQNianHZhangSHThe clinicopathological significance of BUBR1 overexpression in hepatocellular carcinomaJ Clin Pathol6210031008200910.1136/jcp.2009.06694419861558 FuXChenGCaiZDWangCLiuZZLinZYWuYDLiangYXHanZDLiuJCZhongWDOverexpression of BUB1B contributes to progression of prostate cancer and predicts poor outcome in patients with prostate cancerOnco Targets Ther922112220201627143916 ShichiriMYoshinagaKHisatomiHSugiharaKHirataYGenetic and epigenetic inactivation of mitotic checkpoint genes hBUB1 and hBUBR1 and their relationship to survivalCancer Res621317200211782350 ParkHYJeonYKShinHJKimIJKangHCJeongSJChungDHLeeCWDifferential promoter methylation may be a key molecular mechanism in regulating BubR1 expression in cancer cellsExp Mol Med39195204200710.1038/emm.2007.2217464181 WilliamsonAWickliffeKEMelloneBGSongLKarpenGHRapeMIdentification of a physiological E2 module for the human anaphase-promoting complexProc Natl Acad Sci USA1061821318218200910.1073/pnas.090788710619822757 ZhangHQZhaoGKeBMaGLiuGLLiangHLiuLRHaoXSOverexpression of UBE2C correlates with poor prognosis in gastric cancer patientsEur Rev Med Pharmacol Sci2216651671201829630110 MoCHGaoLZhuXFWeiKLZengJJChenGFengZBThe clinicopathological significance of UBE2C in breast cancer: A study based on immunohistochemistry, microarray and RNA-sequencing dataCancer Cell Int1783201710.1186/s12935-017-0455-129021715 BaviPUddinSAhmedMJehanZBuRAbubakerJSultanaMAl-SaneaNAbduljabbarAAshariLHBortezomib stabilizes mitotic cyclins and prevents cell cycle progression via inhibition of UBE2C in colorectal carcinomaAm J Pathol17821092120201110.1016/j.ajpath.2011.01.03421514426 DastsoozHCeredaMDonnaDOlivieroSA comprehensive bioinformatics analysis of UBE2C in cancersInt J Mol Sci202228201910.3390/ijms20092228 LinYTChenYWuGLeeWHHec1 sequentially recruits Zwint-1 and ZW10 to kinetochores for faithful chromosome segregation and spindle checkpoint controlOncogene2569016914200610.1038/sj.onc.120968716732327 GaoTHanYYuLAoSLiZJiJCCNA2 is a prognostic biomarker for ER+ breast cancer and tamoxifen resistancePLoS One9e91771201410.1371/journal.pone.009177124622579 LiCYXueCYangQLowBCLiouYCNuSAP governs chromosome oscillation by facilitating the Kid-generated polar ejection forceNat Commun710597201610.1038/ncomms1059726839278 GordonCAGongXGaneshDBrooksJDNUSAP1 promotes invasion and metastasis of prostate cancerOncotarget82993529950201710.18632/oncotarget.1560428404898 LiuZGuanCLuCLiuYNiRXiaoMBianZHigh NUSAP1 expression predicts poor prognosis in colon cancerPathol Res Pract214968973201810.1016/j.prp.2018.05.01729853313 WangQSuLLiuNZhangLXuWFangHCyclin Dependent Kinase 1 Inhibitors: A review of recent progressCurr Med Chem1820252043201110.2174/09298671179559011021517772 PetriETErricoAEscobedoLHuntTBasavappaRThe crystal structure of human cyclin BCell Cycle613421349200710.4161/cc.6.11.429717495533 ZhuangLYangZMengZUpregulation of BUB1B, CCNB1, CDC7, CDC20, and MCM3 in tumor tissues predicted worse overall survival and disease-free survival in hepatocellular carcinoma patientsBiomed Res Int20187897346201810.1155/2018/789734630363964 BrownNRLoweEDPetriESkamnakiVAntrobusRJohnsonLNCyclin B and Cyclin A confer different substrate recognition properties on CDK2Cell Cycle613501359200710.4161/cc.6.11.427817495531 QiaoMShapiroPFosbrinkMRusHKumarRPassanitiACell cycle-dependent phosphorylation of the RUNX2 transcription factor by cdc2 regulates endothelial cell proliferationJ Biol Chem28171187128200610.1074/jbc.M50816220016407259 PaganoMPepperkokRVerdeFAnsorgeWDraettaGCyclin A is required at two points in the human cell cycleEMBO J11961971199210.1002/j.1460-2075.1992.tb05135.x1312467 DavidYZivTAdmonANavonAThe E2 ubiquitin-conjugating enzymes direct polyubiquitination to preferred lysinesJ Biol Chem28585958604201010.1074/jbc.M109.08900320061386

Heat maps of DEGs between the two databases. (A) Heat map of GEO datasets. (B) Heat map of TCGA datasets. Red lines represent upregulation and blue lines represent downregulation. DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; TCGA, The Cancer Genome Atlas; N, normal tissues; T, tumor tissues.

Volcano plots of differentially expressed genes used to distinguish HCC from healthy tissue. Red dots represent genes that were significantly upregulated in HCC; blue dots represent genes that were significantly downregulated in HCC; grey dots represent genes that were not differentially expressed. (A) Gene Expression Omnibus datasets. (B) The Cancer Genome Atlas datasets. HCC, hepatocellular carcinoma.

Functional enrichment analyses of DEGs. (A) GOChord plot of GO analyses, with the top ten GO terms. The left side of the GOChord diagram represents logFC, and the right side represents different GO terms. The connecting bands indicate the corresponding pathways for each gene. (B) GOCircle plot of KEGG pathway analyses, including the top six KEGG pathways. The height of the inner ring indicates the significance of the term (−log10 P-value), and the color represents the z-score. The outer ring indicates the logFC of the genes in each GO term. DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and genomes.

PPI network. (A) PPI network of DEGs with red marking upregulated and blue marking downregulated genes. (B) The most significant module in the PPI network consisted of 45 nodes and 922 edges. PPI, protein-protein interaction; DEGs, differentially expressed genes.

Selection and analysis of hub genes. Analyses of (A) disease-free and (B) overall survival rate of hub genes. CDC20, cell division cycle 20; UBE2C, ubiquitin-conjugating enzyme E2 C; BUB1B, BUB1 mitotic checkpoint serine/threonine kinase β; CCNB1, cyclin-B1.

Comparison of the expression difference of two core genes between tumors and healthy tissues. UBE2C, ubiquitin-conjugating enzyme E2C; BUB1B, BUB1 mitotic checkpoint serine/threonine kinase β.

Association between BUB1B mRNA expression and (A) tumor grade, (B) hepatitis virus infection status, (C) satellites and (D) vascular invasion. Association between UBE2C mRNA expression and (E) tumor grade, (F) hepatitis virus infection status, (G) satellites and (H) vascular invasion. UBE2C, ubiquitin-conjugating enzyme E2C; BUB1B, BUB1 mitotic checkpoint serine/threonine kinase β.

Enrichment plots from gene set enrichment analysis. (A) Signaling pathways of BUB1B. (B) Signaling pathways of UBE2C. UBE2C, ubiquitin-conjugating enzyme E2 C; BUB1B, BUB1 mitotic checkpoint serine/threonine kinase β; KEGG, Kyoto Encyclopedia of Genes and Genomes; WNT, wingless-related integration; VEGF, vascular endothelial growth factor.

DEGs (n=151) in The Cancer Genome Atlas and Gene Expression Omnibus databases.

DEGs Genes
Downregulated MT1M, HAMP, CNDP1, SLC22A1, APOF, CLEC4M, FCN3, CLEC1B, CRHBP, KCNN2, OIT3, CYP1A2, CLEC4G, FOS, C9, MT1F, GYS2, CYP2C19, FOSB, TTC36, MT1H, HAO2, GBA3, SLCO1B3, CYP26A1, MT1E, NAT2, CYP2B6, PCK1, BCO2, MT1X, LPA, MARCO, CYP4A11, NPY1R, MT1G, SLC25A47, GHR LYVE1, ASPG PLAC8, CXCL14, IGFALS, CFP, TSLP, FCN2, IL1RAP, FBP1, DIRAS3, RDH16, PZP, MT2A, RND3, AADAT, CYP39A1, CETP
Upregulated KIF15, SHCBP1, MMP12, ESM1, STK39, UHRF1, CENPE, CEP55, COX7B2, TPX2, S100P, GMNN, ZWINT, CENPL, TRIM71, EZH2, BEX2, MDK, E2F8, TRIP13, CCNA2, CRNDE, NEK2, CD24, CDK1, KIF4A, DEPDC1B, PEG10, ACSL4, FOXM1, MELK, CDC20, FLVCR1, ZIC2, OIP5, FAM83D, NQO1, NUSAP1, RBM24, ANLN, SULT1C2, ASPM, NCAPG, NDC80, UBE2T, UBE2C, RAD51AP1, TOP2A, PRC1, DLGAP5, KIF20A, DTL, HMMR, TTK, NUF2, ROBO1, CENPW, BUB1B, GPC3, GINS1, PBK, CCNB1, CTHRC1, COL15A1, CDKN3, RACGAP1, CAP2, RRM2, AKR1B10, SPINK1

DEGs, differentially expressed genes.

Function of eight hub genes.

No. Gene symbol Protein Function (Refs)
1 CCNB1 G2/mitotic-specific cyclin-B1 Control of the cell cycle at the G2/M (mitosis) transition Brown et al, 2007 (58); Petri et al, 2007 (56).
2 CDK1 Cyclin-dependent kinase 1 Regulation of cell cycle progression, apoptosis and carcinogenesis of tumor cells Qiao et al, 2006 (59); Wang et al, 2011 (55).
3 BUB1B BUB1 mitotic checkpoint serine/ threonine kinase β Essential component of the mitotic checkpoint; mutated in colorectal cancer and other tumors Shin et al, 2003 (38); Yamamoto et al, 2007 (40); Liu et al, 2009 (41); Fu et al, 2016 (42).
4 CDC20 Cell division cycle 20 Regulatory protein that participates in cell cycle processes Zhang et al, 2019 (36); Cheng et al, 2019 (37).
5 CCNA2 Cyclin-A2 Controls the G1/S and the G2/M transition phases of the cell cycle Pagano et al, 1992 (60).
6 UBE2C Ubiquitin-conjugating enzyme E2 C Involved in ubiquitination during protein modification Williamson et al, 2009 (45); David et al, 2010 (61).
7 ZWINT ZW10 interactor Member of the MIS12 complex; participates in biological processes such as cell cycle and cell division Lin et al, 2006 (50).
8 NUSAP1 Nucleolar and spindle- associated protein 1 High expression may be associated with prostate, colon and liver cancer Gordon et al, 2017 (53); Liu et al, 2018 (54); Zhou et al, 2018 (34)

Univariate and multivariate analyses of the correlation of BUB1B expression with overall survival.

Univariate analysis Multivariate analysis
Parameter HR 95% CI P-value HR 95% CI P-value
Age 1.005 0.987–1.023 0.5912 1.013 0.993–1.033 0.2144
Sex 0.780 0.487–1.249 0.3012 1.105 0.654–1.869 0.7088
Histological grade 1.017 0.746–1.387 0.9143 1.003 0.715–1.406 0.9871
TNM stage 1.865 1.456–2.388 <0.0001 0.961 0.361–2.556 0.9366
T classification 1.804 1.434–2.270 <0.0001 1.690 0.702–4.066 0.2414
M classification 3.850 1.207–12.281 0.0228 2.313 0.580–9.220 0.2346
N classification 2.022 0.494–8.276 0.3276 1.690 0.266–10.728 0.5779
BUB1B 1.329 1.157–1.527 <0.0001 1.295 1.111–1.508 <0.0001

HR, hazard ratio; CI, confidence interval; TNM, tumor-node-metastasis; BUB1B, BUB1 mitotic checkpoint serine/threonine kinase β.

Univariate and multivariate analyses of the correlation of UBE2C expression with overall survival.

Univariate analysis Multivariate analysis
Parameter HR 95% CI P-value HR 95% CI P-value
Age 1.005 0.987–1.023 0.5912 1.014 0.994–1.035 0.1755
Sex 0.780 0.487–1.249 0.3012 1.123 0.666–1.894 0.6628
Histological grade 1.017 0.746–1.387 0.9143 0.983 0.703–1.375 0.9189
TNM stage 1.865 1.456–2.388 <0.0001 0.941 0.360–2.462 0.9011
T classification 1.804 1.434–2.270 <0.0001 1.728 0.729–4.093 0.2138
M classification 3.850 1.207–12.281 0.0227 1.970 0.507–7.656 0.3277
N classification 2.022 0.494–8.276 0.3275 1.904 0.314–11.554 0.4840
UBE2C 1.374 1.194–1.582 <0.0001 1.332 1.145–1.550 <0.0001

HR, hazard ratio; CI, confidence interval; TNM, tumor-node-metastasis; UBE2C, ubiquitin-conjugating enzyme E2C.