Screening of significant biomarkers with poor prognosis in hepatocellular carcinoma via bioinformatics analysis.

Hepatocellular carcinoma (HCC) is a malignant tumor with unsatisfactory prognosis. The abnormal genes expression is significantly associated with initiation and poor prognosis of HCC. The aim of the present study was to identify molecular biomarkers related to the initiation and development of HCC via bioinformatics analysis, so as to provide a certain molecular mechanism for individualized treatment of hepatocellular carcinoma.Three datasets (GSE101685, GSE112790, and GSE121248) from the GEO database were used for the bioinformatics analysis. Differentially expressed genes (DEGs) of HCC and normal liver samples were obtained using GEO2R online tools. Gene ontology term and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway analysis were conducted via the Database for Annotation, Visualization, and Integrated Discovery online bioinformatics tool. The protein-protein interaction (PPI) network was constructed by the Search Tool for the Retrieval of Interacting Genes database and hub genes were visualized by Cytoscape. Survival analysis and RNA sequencing expression were conducted by UALCAN and Gene Expression Profiling Interactive Analysis.A total of 115 shared DEGs were identified, including 30 upregulated genes and 85 downregulated genes in HCC samples. P53 signaling pathway and cell cycle were the major enriched pathways for the upregulated DEGs whereas metabolism-related pathways were the major enriched pathways for the downregulated DEGs. The PPI network was established with 105 nodes and 249 edges and 3 significant modules were identified via molecular complex detection. Additionally, 17 candidate genes from these 3 modules were significantly correlated with HCC patient survival and 15 of 17 genes exhibited high expression level in HCC samples. Moreover, 4 hub genes (CCNB1, CDK1, RRM2, BUB1B) were identified for further reanalysis of KEGG pathway, and enriched in 2 pathways, the P53 signaling pathway and cell cycle pathway.Overexpression of CCNB1, CDK1, RRM2, and BUB1B in HCC samples was correlated with poor survival in HCC patients, which could be potential therapeutic targets for HCC.

Introduction

Currently, the incidence and mortality rates of cancer overall are declining year by year, whereas for liver cancer, both incidence and mortality rates are increasing.[ Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer.[ Although great progress has been made in the treatment of HCC in recent years, the prognosis of HCC remains unsatisfactory. It is estimated that 5 million deaths worldwide between 2015 and 2030 will be attributable to HCC due to its difficulty in early diagnosis, relapse, and metastasis.[ Accumulating evidence indicated that the abnormal genes expression is associated with initiation and poor prognosis of HCC,[ and numerous studies have been performed to search for suitable biomarkers of early detection and prognosis.[ However, limit number of biomarkers are used in clinical practice.[ Therefore, it is especially important to explore more reliable biomarkers for improving the early detection, prognosis prediction, and better understanding the potential mechanism of HCC.

In the last decades, genome-wide molecular profiling has played crucial roles in quickly detecting differentially expressed genes that were involved in the tumorigenesis and progression and has proved to be a reliable technique to identify core genes.[ In this study, we aim to identify molecular biomarkers related to the initiation and development of hepatocellular carcinoma via bioinformatics analysis, so as to provide a certain molecular mechanism for individualized treatment of hepatocellular carcinoma in the future.

Methods

Microarray data

In this study, we obtained the 3 independent gene expression profiles of HCC (GSE101685, GSE112790, and GSE121248) from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). All of these 3 microarray data were based on GPL570 Platforms ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array), which included 24 HCC samples and 8 normal liver samples, 183 HCC samples, and 15 normal liver samples and 70 HCC samples and 37 normal liver samples, respectively. Ethical approval was not necessary as this study is bioinformatics analysis.

Identification of differentially expressed genes (DEGs)

GEO2R online tools[ were used to identify the DEGs between HCC samples and normal liver samples. The genes with |logFC (fold change)| > 2 and adjust P value < .05 were considered as the DEGs. Venny online software (https://bioinfogp.cnb.csic.es/tools/venny/) was used to identify the shared upregulated and downregulated DEGs among the 3 datasets.

Gene ontology term and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway analysis

To investigate the function of commonly DEGs, the Database for Annotation, Visualization, and Integrated Discovery online bioinformatics tool (DAVID, https://david.ncifcrf.gov/) were used to conduct functional and pathway enrichment analysis. GO analysis, including the biological process, cellular component, and molecular function, and KEGG pathway analysis, were performed for the shared DEGs via DAVID. P value less than .05 was considered statistical significance.

Protein–protein interaction (PPI) network and module analysis

To construct PPI network for the commonly DEGs, an online biological database, the Search Tool for the Retrieval of Interacting Genes (STRING; https://string-db.org/), was used. Then, the potential correlation between these DEGs was evaluated by STRING app in Cytoscape.[ Furthermore, the modules of the PPI network were explored by Molecular Complex Detection (MCODE) app in Cytoscape (degree cutoff = 2, K-Core = 4, Depth = 100, and node score cutoff = 0.2).

Survival analysis and RNA sequencing expression

The analysis of patient survival of the hub genes was performed via UALCAN (http://ualcan.path.uab.edu/), a comprehensive, widely used online bioinformatics tool, which provide graphs and plots depicting gene expression and patient survival information based on The Cancer Genome Atlas (TCGA) database and MET500 transcriptome sequencing.[ To validate these hub genes, Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) was applied to assess the RNA sequencing expression of candidate genes from the TCGA database and Genotype-Tissue Expression Projects.[

Results

Identification of DEGs in HCC

In this study, there were 277 HCC samples and 60 normal liver samples. By using GEO2R online tools, 459, 325, and 145 DEGs were extracted from GSE101685, GSE112790, and GSE121248, respectively. Then, Venny online software was used to identify the shared upregulated and downregulated DEGs among the 3 datasets. As shown in Fig. 1 and Table 1, a total of 115 shared DEGs were identified, including 30 upregulated genes and 85 downregulated genes in HCC samples.

Figure 1

Identification of common DEGs in 3 gene expression datasets (GSE101685, GSE112790, and GSE121248). A, Thirty upregulated DEGs in the 3 datasets. B, Eighty-five downregulated DEGs in the 3 datasets. DEGs = differentially expressed genes.

Table 1

All 115 differentially expressed genes (DGEs) were identified from 3 profile datasets, including 30 upregulated genes and 85 downregulated genes in the hepatocellular carcinoma tissues, compared with normal liver tissues.

Functional enrichment analysis of DEGs genes

To further demonstrate the biological functions of the DEGs, Gene ontology in DAVID software was performed. As shown in Table 2, for biological process, the upregulated DEGs were mainly involved in mitotic nuclear division, regulation of attachment of spindle microtubules to kinetochore, G2/M transition of mitotic cell cycle, negative regulation of ubiquitin–protein ligase activity involved in mitotic cell cycle, and positive regulation of ubiquitin–protein ligase activity involved in regulation of mitotic cell cycle. Whereas the downregulated DEGs were significantly involved in epoxygenase P450 pathway, oxidation-reduction process, exogenous drug catabolic process, xenobiotic metabolic process, and monocarboxylic acid metabolic process. For cellular component, the upregulated DEGs were mainly involved in midbody, cytoplasm, mitotic spindle, centralspindlin complex, and cytosol. Whereas the downregulated DEGs were significantly enriched in extracellular region, organelle membrane, insulin-like growth factor ternary complex, extracellular space, and insulin-like growth factor binding protein complex. About molecular function, the upregulated DEGs were mainly involved in protein binding, histone kinase activity, protein kinase activity, and protein serine/threonine kinase. Whereas the downregulated DEGs were significantly enriched in heme binding, iron ion binding, oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, arachidonic acid epoxygenase activity and oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as 1 donor, and incorporation of 1 atom of oxygen.

Table 2

The top 5 pathways in gene ontology analysis of DEGs in hepatocellular carcinoma.

KEGG analysis indicated that the upregulated DEGs were significantly enriched in P53 signaling pathway and cell cycle. While downregulated DEGs were mainly involved in retinol metabolism, caffeine metabolism, and drug metabolism—cytochrome P450, etc (Table 3).

Table 3

KEGG pathway analysis of DEGs in hepatocellular carcinoma.

PPI network and module analysis

To investigate the relationship between the identified DEGs in the HCC, we constructed the PPI network via online database STRING and displayed it via Cytoscape software.[ As shown in Fig. 2A, there were 105 nodes and 249 edges. In addition, 3 significant modules were identified via MCODE. All of the 15 nodes were upregulated DEGs in the module 1 (Fig. 2B). Whereas in the module 2, the other 8 nodes were downregulated DEGs expect SPP1 (Fig. 2C). In the module 3, all of the 6 nodes were downregulated DEGs (Fig. 2D). Interestingly, most of the upregulated genes have the highest connectivity with each other (Fig. 2B).

Figure 2

Protein–protein interaction (PPI) network construction and 3 significant modules. A, PPI relationships among common DEGs from the 3 datasets of GEO. B, Module 1. C, Module 2. D, Module 3. Red nodes represent upregulated DEGs and green nodes represent downregulated DEGs.

Survival analysis of hub genes

UALCAN was utilized to identify 30 hub genes survival data from 3 modules mentioned above. As shown in Fig. 3, 17 candidate genes were significantly correlated with HCC patient survival. Briefly, 14 of the 15 genes in module 1 experienced poor survival; no significant difference between PRC1 expression and survival in HCC patients. In the module 2, HCC patients with higher SPP1 or lower ESR1 expression exhibited worse survival. However, the upregulation of SERPINE1 gene is correlated with poor survival, which is contrary to the results of the microarray data. The other 6 genes had no significant. In the module 3, the upregulation of CYP2C9 gene is correlated with better survival, whereas the other 5 genes had no significant. Next, GEPIA was applied to assess the RNA sequencing expression of 17 candidate genes. Results demonstrated that 15 of 17 genes exhibited high expression level in HCC samples and 2 of 17 genes reflected low expression level in HCC samples compared with normal liver samples (Fig. 4), which is consistent with the microarray array data.

Figure 3

Overall survival validation of hub genes in HCC patients. HCC = hepatocellular carcinoma.

Figure 4

Relative expression of hub genes in HCC samples and normal liver samples. HCC = hepatocellular carcinoma.

Reanalysis of 17 selected genes for pathway enrichment

To investigate the possible pathway of these 17 hub genes, DAVID software was used to reanalyze the pathway enrichment. As shown in Table 4 and Fig. 5, 4 genes were enriched in 2 pathways, the P53 signaling pathway and cell cycle (all P value less than.05).

Table 4

Reanalysis of 17 selected genes for pathway enrichment.

Figure 5

Reanalysis of 17 hub genes via KEGG pathway enrichment. A, CCNB1, CDK1, and BUB1B are enriched in the Cell cycle. B, CCNB1, CDK1, and RRM2 are enriched in the P53 signaling pathway.

Discussion

In the present study, a total of 115 shared DEGs were identified from 3 independent gene expression profiles via bioinformatics analysis, including 30 upregulated genes and 85 downregulated genes in HCC samples. Upregulated DEGs were mainly enriched in mitotic nuclear division and cell cycle related functional terms, whereas downregulated DEGs were mainly enriched in multiple metabolic-related functional terms. Among these DEGs, 3 modules and 30 vital genes were screened from the PPI network complex. In addition, through survival analysis of these 30 vital genes, we found that 15 upregulated genes and 2 downregulated genes were significantly associated with worse survival of HCC patients. Furthermore, reanalysis of 17 selected genes for KEGG pathway enrichment clarified that 4 hub genes (CCNB1, cyclin-dependent kinases [CDK1], RRM2, BUB1B) enriched in P53 signaling pathway and cell cycle, and all of these 4 hub genes were upregulated in the HCC samples compared with normal liver samples, indicating the “driver” function of these genes in HCC development.

Previous studies indicated that cell cycle dysregulation is a hallmark of cancer,[ and target cell cycle could be an effective therapeutic strategy in the treatment of cancer. Increasing evidence indicated that cell cycle-related genes such as CCNB1, CDK1, and BUB1B, which were also screened in this study, are involved in the initiation and progression of cancers.[ CCNB1 is overexpressed in many cancers and correlated with worse survival such as pancreatic cancer,[ breast cancer,[ and prostate cancer[ and could promote the G2/M transition of the cell cycle by regulating and forming a complex with CDK1.[ The transcription factor p53, located on chromosome 17p, plays an essential role in regulating cell cycle, DNA repair, and apoptosis and acts as a powerful tumor suppressor.[ However, in most tumors, p53 gene is inactivated by mutation or degradation of its associated pathway.[ Zhang et al[ have demonstrated that CCNB1 silencing could suppress pancreatic cancer cell proliferation and promote apoptosis through activating the p53 signaling pathway. Chai et al[ have reported that FOXM1 could promote proliferation of HCC cells via transcriptional activation of CCNB1. These reports are in line with our present demonstration that CCNB1 was overexpressed in HCC samples, which was correlated with worse survival of HCC patients.

Accumulating studies demonstrated that CDKs play essential roles in the regulation of cell cycle.[ A recent study indicated that karyopherin subunit-α 2 could promote cell proliferation and induce cell cycle arrest via upregulating CDK1 and CCNB2 in HCC.[ Another report revealed that miR-582–5p inhibited the HCC progression in vitro via targeting CDK1 and AKT3.[ In other tumors, suppressing the expression of CDK1 could significantly inhibit cell proliferation and induce apoptosis in breast cancer.[ In addition, overexpression of miR-31–5p significantly inhibited renal cell carcinoma progression in vitro through targeting CDK1.[ Similar association was identified in the present study that CDK1 was upregulated in HCC and higher expression of CDK1 correlated with worse survival of HCC patients. Taken together, these data indicated that CDK1 may act as a potential biomarker for predicting the survival in HCC.

Mitotic checkpoint serine/threonine kinase B (BUB1B), known also as BubR1, encodes a spindle checkpoint-related kinase and plays an important role in the proper chromosome segregation.[ Accumulating data reveal that BUB1B acts as an essential role in the development and progression of many forms of cancer including breast,[ prostate,[ lung,[ and liver cancer.[ A recent investigation revealed that high level of BUB1B was associated with worse survival in HCC patients,[ which is in agreement with our present study.

Human ribonucleotide reductase (RR) is the rate-limiting enzyme catalyzing the conversion of ribonucleoside diphosphates to deoxyribonucleoside diphosphates in cells. RR is composed of 2 identical subunits, RRM1 and RRM2. The expression level of RRM1 is relatively stable during the course of cell cycle, while the expression of RRM2 only occurred in the late G1 and early S phase and functioned as an essential role in DNA synthesis and repair.[ Therefore, RRM2 is identified as an important anticancer target.[ Several investigations demonstrated that the hepatitis B virus (HBV) could induce RRM2 expression through activating DNA damage response and targeting RRM2 could effectively inhibit HBV replication.[ These studies indicated that RRM2 may be identified as a therapeutic target in HBV-related HCC. Indeed, upregulation of RRM2 was correlated with worse survival in HCC,[ which is in line with our study. Taken together, these data indicated that RRM2 may act as a potential therapeutic target for HCC.

In brief, our bioinformatics analysis identified 4 hub genes (CCNB1, CDK1, RRM2, BUB1B) which were correlated with poor survival in HCC patients. However, the specific mechanism of these hub genes involved in regulating the occurrence, development, and prognosis of HCC still needs to be further explored.

Author contributions

Conceptualization: Tongxin Liu.

Formal analysis: Peng Liu.

Funding acquisition: Tongxin Liu, Peng Liu.

Investigation: Bin Long.

Methodology: Quanquan Sun.

Writing-original draft: Quanquan Sun.

Writing-review & editing: Yuan Zhu.