Abstract

Hepatocellular carcinoma (HCC) is a malignant tumor with high mortality. The abnormal expression of genes is significantly related to the occurrence of HCC. The aim of this study was to explore the differentially expressed genes (DEGs) of HCC and to provide bioinformatics basis for the occurrence, prevention and treatment of HCC. The DEGs of HCC and normal tissues in GSE102079, GSE121248, GSE84402 and GSE60502 were obtained using R language. The GO function analysis and KEGG pathway enrichment analysis of DEGs were carried out using the DAVID database. Then, the protein–protein interaction (PPI) network was constructed using the STRING database. Hub genes were screened using Cytoscape software and verified using the GEPIA, UALCAN, and Oncomine database. We used HPA database to exhibit the differences in protein level of hub genes and used LinkedOmics to reveal the relationship between candidate genes and tumor clinical features. Finally, we obtained transcription factor (TF) of hub genes using NetworkAnalyst online tool. A total of 591 overlapping up-regulated genes were identified. These genes were related to cell cycle, DNA replication, pyrimidine metabolism, and p53 signaling pathway. Additionally, the GEPIA database showed that the CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 were associated with the poor survival of HCC patients. UALCAN, Oncomine, and HPA databases and qRT-PCR confirmed that these genes were highly expressed in HCC tissues. LinkedOmics database indicated these genes were correlated with overall survival, pathologic stage, pathology T stage, race, and the age of onset. TF analysis showed that MYBL2, KDM5B, MYC, SOX2, and E2F4 were regulators to these nine hub genes. Overexpression of CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 in tumor tissues predicted poor survival in HCC. They may be potential therapeutic targets for HCC.

1. Introduction

Currently, liver cancer is the fifth most common malignancy and the second most common cause of cancer death [1, 2]. Hepatocellular carcinoma (HCC) accounts for more than 90% of primary liver cancer, killing about 750,000 people worldwide each year [3]. In particular, the vast majority of cases (83%) have occurred in less developed regions of the world, causing a major health crisis in Asia [4]. Because 80–90% of liver cancer cannot be completely resected, and its' prognosis is very poor, it seriously threatens people's physical and mental health. According to the epidemiological investigation, it may be related to viral hepatitis [5], and drinking, moldy food, toxicants and genetic factors [6]. At present, more and more studies believe that the occurrence and poor prognosis of HCC are related to the abnormal expression of genes [7, 8]. However, since multiple genes are often involved in the process of cell carcinogenesis, and these genes can interact with each other and function through the regulatory network [9], the specific pathogenesis of HCC is still unclear. Although great efforts have been made to search for biomarkers of tumor prognosis or diagnosis, it is estimated that less than 1% of biomarkers are used in clinical practice [10]. Therefore, it is especially important to assess the relationship between biomarkers and diseases at the genetic level, protein levels and clinical factors.

In recent years, gene chip technology has played an important role in studying tumor gene expression profiles and searching for tumor key genes [11]. This study aims to provide bioinformatics basis for further research on the molecular mechanism of HCC, and provide a new way to carry out individualized treatment on genes level.

2. Materials and Methods

2.1. Microarray Data

In this study, the gene expression profile datasets (GSE102079, GSE121248, GSE84402, and GSE60502) were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) of NCBI. The microarray data from GSE102079, GSE121248, and GSE84402 were based on GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array. The GSE102079 profile was composed of 152 HCC tissues and 91 nontumorous tissues. GSE121248 included 70 chronic hepatitis B-induced HCC tissues and 39 adjacent normal tissues. GSE84402 included 14 pairs of HCC tissues and corresponding nontumorous tissues. Totally, 18 pairs of HCC tissues and adjacent nontumorous tissues were enrolled in GSE60502 which was based on GPL96[HG-U133A] Affymetrix Human Genome U133A Array.

2.2. Screening Overlapping Up-Regulated DEGs in HCC

We used R software to analyze GSE102079, GSE121248, GSE84402, and GSE60502 raw data of the CEL file for identifying DEGs. RMA package was used for data normalization processing. Affy package was used for quality assessment of samples in each GEO dataset. The Limma package was used to identify DEGs. The criterion for selection of DEGs was set as |log2FC| > 1 and value <0.05 for each GEO dataset. To identify shared up-regulated DEGs among GSE102079, GSE121248, GSE84402 and GSE60502, we used R software to generate a Venn diagram.

2.3. Function and Pathway Enrichment Analyses of Common Up-Regulated DEGs

To investigate the function of 591 common up-regulated DEGs, we used the Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) online tool to perform functional and pathway enrichment analysis. Gene ontology (GO) analysis, including the biological process (BP), cellular component (CC), and molecular function (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were conducted for the selected common up-regulated DEGs using DAVID. was regarded as statistical significance. Subsequently, the shared up-regulated DEGs enriched in top ten KEGG pathways/GO functions were determined for the Venn diagram using R software.

2.4. Protein–Protein Interaction (PPI) Network Construction and Hub Genes Selection

We used the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/) online database to construct PPI network. CytoHubba was used to get the top 10 hub genes with the highest degree in the PPI network in Cytoscape software.

2.5. Validation of the Hub Genes

We used Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn) to identify potential candidate biomarkers for overall survival (OS) and disease-free survival (DFS) in liver hepatocellular carcinoma (LIHC) patients. Genes that are significantly associated with OS and DFS were considered as potential biomarkers for LIHC prognosis. Furthermore, to evaluate mRNA expression of hub genes, we used UALCAN (http://ualcan.path.uab.edu/analysis.html) database and Oncomine (https://www.oncomine.org/) database to differentiate expression of hub genes in LIHC tissues and normal tissues.

2.6. Quantitative Real-Time PCR

HCC and adjacent tissues were taken from Cuiying center sample library of Lanzhou University Second Hospital. The total RNA was extracted by using TRNzol Reagent, and was reverse-transcribed with FastKing gDNA Dispelling RT SuperMix (TIANGEN, Beijing, China). All qRT-PCR reactions were conducted with Rotor-Gene 6000 PCR system (Qiagen) and performed with SsoFast EvaGreen Supermix (Bio-Rad) in 20 μl volume containing 10 μl of 2× SsoFast EvaGreen Supermix, 1 μl of each 10 μM forward and reverse primer, 1 μl of cDNA sample, and nuclease-free water up to 20 μl. Amplification was carried out according to the following conditions: initial denaturation 95°C 5 min, followed by 45 cycles of denaturation 95°C 10 s, annealing 57°C 15 s, extension 72°C 15 s. The relative expression of the gene was calculated by the 2−ΔΔCt method. The primers are listed in Table 1.

2.7. Evaluation of Immunohistochemical Staining

To verify the protein expression level of candidate genes in HCC tissues, we used Human Protein Atlas (HPA, https://www.proteinatlas.org/) database to obtain immunohistochemical staining.

2.8. Relationship between Candidate Genes and Clinical Features in HCC Patients

To further explore the relationship between candidate genes and tumor clinical features, we analyzed the TCGA clinical data using LinkedOmics (http://www.linkedomics.org/) database.

2.9. Transcription Factor (TF) and Expression Correlation Analyses

TF of hub genes was explored using NetworkAnalyst (http://www.networkanalyst.ca). Expression correlation analysis based on TCGA samples was conducted in GEPIA.

2.10. Statistical Analysis

Statistical analysis and graphs were performed with GraphPad Prism 7.00 sofware. Data were presented as the mean ± SD. The test was used for comparison between the two groups.

3. Results

3.1. Screening Overlapping Up-Regulated DEGs

We processed the data of four chips with R language, and set the cut-off criteria as |log2FC| > 1, the -value <0.05 to screen the DEGs. A total of 591 overlapping up-regulated genes (2766 in GSE102079, 2483 in GSE121248, 2448 in GSE84402, and 3284 in GSE60502) were identified using a Venn diagram (Figure 1).

3.2. Functions and Pathways of Up-Regulated Genes

We presented the top ten pathways/GO functions in this study (Table 2). DEGs were mainly involved in biological processes such as cell division, sister chromatid cohesion, DNA replication, mitotic nuclear division, and G1/S transition of the mitotic cell cycle. Cytological composition analysis showed that most of these genes were involved in the composition of nucleoplasm, cytosol, nucleus, cytoplasm and membrane. The molecular functions were mainly concentrated in protein binding, poly (A) RNA binding, and ATP binding. KEGG pathway showed that the DEGs were mainly involved in the cell cycle, DNA replication, pyrimidine metabolism, and p53 signaling pathway.

Additionally, 104 genes were enriched in KEGG pathways, 133 genes in biological processes, 500 genes in cellular component, and 485 genes in molecular function. Subsequently, we generated a Venn diagram for pathways/GO functions and obtained 51 overlapping genes. AURKA, BUB1, BUB1B, BUB3, CCNA2, CCNB1, CCNB2, CCNE1, CCNE2, CDC20, CDC25C, CDC45, CDC6, CDC7, CDK1, CDK4, CDK7, CDKN2A, CHEK1, DUT, FEN1, MAD2L1, MCM2, MCM3, MCM4, MCM5, MCM6, MCM7, NUP107, NUP133, NUP155, NUP205, NUP37, NUP43, NUP85, NUPL2, PCNA, POLA1,POLE2, PTTG1, RAD21, RAN, RFC1, RFC3, RFC4, RRM1, RRM2, SMC3, TACC3,TPR, and XPO1 were shared in the four datasets (Figure 2).

3.3. PPI Network Construction and Hub Genes Selection

We constructed the PPI network for 51 overlapping up-regulated genes among KEGG pathway and GO analysis. The PPI included 51 nodes and 836 edges (Figure 3(a)). We used cytoHubba to get the top 10 hub genes with the highest degree of connectivity in the PPI network, and the top 10 hub genes included CDK1, CCNB2, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 (Figure 3(b)). The connectivity degree of top 10 hub genes is shown in Table 3.

3.4. Survival Analysis of Top Ten Up-Regulated Genes

We used GEPIA database to get the survival curves of 182 pairs of HCC tissues with high expression and low expression of hub genes. LIHC patients with high CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 experienced poor OS; there was no statistical difference between CCNB2 expression and OS in LIHC patients (Figure 4). Similarly, compared with the low expression of these 10 genes, overexpression of CDK1, CCNB2, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 in tumors was significantly associated with DFS in LIHC patients (Figure 5). Therefore, CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 were considered as potential biomarkers.

3.5. Validation of Selected Up-Regulated Genes in HCC

Using the TCGA data in UALCAN online tool, we analyzed the expression of the nine selected up-regulated genes in LIHC tissues (371 cases) and normal tissues (50 cases). The results showed that the CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 were highly expressed in LIHC tissues, and the differences were statistically significant (Figure 6). Similarly, we also investigated the transcriptional levels of these genes in liver cancer and normal samples by using the Oncomine database. The mRNA expression levels of these genes were significantly up-regulated in liver cancer tissues compared with normal tissues in several datasets (Table 4). Furthermore, these 9 genes were also highly expressed in various grades of HCC compared with the normal group. In addition, overexpression of these 9 genes was also related to advanced tumor grade (Figure 7).

To further validate the data mining results, we performed qRT-PCR with paired tumor and adjacent tissues borrowed from the Cuiying center sample library of Lanzhou University Second Hospital. Although the sample was limited, except for CNB1, the other 8 candidate genes were highly expressed in tumor tissues (Figure 8).

3.6. Differences of Selected Up-Regulated Genes in Protein Level between HCC and Normal Tissues

We used the HPA database to exhibit the differences in protein level of CDK1, CCNB1, CDC20, MAD2L1, MCM3, MCM2, and RFC4. The results showed the immunohistochemical staining of CDK1, CCNB1, CDC20, MAD2L1, MCM3, MCM2, and RFC4 was negative staining in normal tissues and positive in HCC tissues, demonstrating that these genes were significantly expressed in HCC tissues than in normal liver tissues. The immunohistochemical staining is displayed in Figure 9.

3.7. Expression of Selected Up-Regulated Genes and Its Clinical Significance in LIHC Patients

Downloading the TCGA clinical data in LinkedOmics online tool, we analyzed the relationship between selected up-regulated genes and clinical features in LIHC patients. The results showed that all 9 candidate genes were significantly correlated with overall survival, pathologic stage, and pathology T stage, indicating that high expression of candidate genes predicted poor survival and tumor progression. The CDK1, CCNB1, CDC20, MCM3, BUB1B, MCM2, and RFC4 in LIHC patients were significantly correlated with race, which were significantly higher in Asian, Black, or African American and White than in American Indian or Alaska native, while BUB1, MAD2L1 were not significantly different in race. There was a significant correlation between CDK1, MCM3, BUB1B, RFC4 and the age of onset. The expression levels of the 9 genes in pathology N stage, pathology M stage, histological type, ethnicity, residual tumor, radiation therapy, and tumor purity were not statistically different (Table 5).

3.8. TF Analysis for Selected Up-Regulated Genes

We further investigated the molecular that can regulate CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4. We used NetworkAnalyst tool to predict the TFs that can regulate the expression of these nine genes. We found five TFs, MYBL2, KDM5B, MYC, SOX2, and E2F4, that can regulate these nine hub genes expression (Figure 10(a)). Correlation analysis showed MYBL2, KDM5B, SOX2, and E2F4 were positively correlated with these nine genes expressions. MYC was positively correlated with CCNB1, BUB1, MAD2L1, MCM3, and BUB1B expressions (Figure 10(b)).

4. Discussion

The occurrence of HCC is a complex biological process. In recent years, a large number of biomarkers have been used in the early diagnosis of HCC [15]. Many anti-HCC mechanisms have also been discovered [16]. However, there is still very little study at the multiple gene levels. The researches at the multi-gene levels can contribute to explore the pathogenesis of cancer. In this study, the data of four gene chips in HCC were analyzed by bioinformatics method. Finally, it was found that CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 were related to the poor survival of HCC patients and the advanced tumor grade. A study had similar results [17]. The difference was that this study was validated at the transcriptional level and protein level. In addition, the study also analyzed the relationship between clinical features and biomarkers: these genes were correlated with overall survival, pathologic stage, pathology T stage, race and the age of onset. To explore the molecular mechanism of HCC, TF-hub genes regulatory network was also constructed. We identified 5 TFs, MYBL2, KDM5B, MYC, SOX2, and E2F4, all of which can regulate the expression of these 9 hub genes and provide more evidences for the elucidation of the mechanism of HCC progression.

Cyclin-dependent kinase 1 (CDK1) belongs to serine/threonine protein kinase family. A recent study has found that metformin can significantly inhibit the proliferation of HCC cells by inducing G2/M arrest and can effectively reduce the expression of CDK1 [18]. This result suggested that CDK1 may be involved in the process of cell proliferation in the cell cycle of HCC. Another study showed the miR-582-5p regulated the progression of HCC through directly inhibiting the expression of CDK1 and AKT3, and indirectly inhibiting the expression of cyclinD1 [19], and supported this theory. In addition, CDK1 is also expressed in other tumors. Studies have shown that CDK1 is active in the cell cycle of several tumor-regulating cell adhesion [20] and can be used as clinical prognostic biomarkers for nonsmall cell lung cancer [21], colon cancer [22], breast cancer [23], and ovarian cancer [24]. In this study, high expression of CCNB1 was closely associated with poor prognosis in HCC patients. This conclusion is further confirmed by a study that knockdown of CCNB1 regulated by microRNA-144 significantly inhibited cell proliferation, migration, and invasion in HCC [25]. In this study, the KEGG pathway showed that the DEGs were mainly involved in the cell cycle and DNA replication. Previous study [26] has shown that CCNB1/CDK1-mediated phosphorylation provides cells with efficient bioenergy for G2/M transition and shortens the overall cell-cycle time. Therefore, CCNB1/CDK1 plays an important role in the cell cycle and cell proliferation.

Overexpression of cell division cycle 20 (CDC20) is associated with poor prognosis of prostate cancer [27], breast cancer [28], and colon cancer [29]. However, the expression of CDC20 in HCC still lacks sufficient experimental data. In cutaneous squamous cell carcinoma, CDC20 promotes cell proliferation and migration through the Wnt/-catenin signaling pathway [30]. CDC20 can contribute to cardiac hypertrophy by promoting LC3 degradation and inhibiting autophagy [31]. High expression of BUB1B can increase proliferation, migration, and invasion of prostate cancer cells [32]. MiR-200c-5p inhibits the proliferation, migration, and invasion of HCC cells by down-regulating MAD2L1 [33], which suggested that the expression of MAD2L1 is significantly higher in HCC and related to the poor prognosis of HCC. Moreover, it also indicated that MAD2L1 can be used as a prognostic and therapeutic target in HCC patients. The minichromosome maintenance (MCM) participates in DNA synthesis [34] and can be used as a biomarker of oral squamous cell carcinoma [35], melanoma [36], glioma [37] and colon cancer [38]. Replication factor C (RFC) plays an important role in DNA repair activities following DNA damage [39]. Targeted therapy of RFC3 can inhibit the proliferation and survival of HCC cells [40].

In this study, we identified DEGs in HCC by bioinformatics analysis and found that overexpression of CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 in tumor tissues predicted poor survival in HCC. We hypothesized that CDK1, CCNB1, CDC20, BUB1, MAD2L1, MCM3, BUB1B, MCM2, and RFC4 may be potential therapeutic targets for HCC. We analyzed these genes at the transcriptional and protein levels, verified with qRT-PCR, explored the relationship between candidate genes and clinical factors, and constructed TF-hub genes regulatory network. Therefore, there are some advantages in this study.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors’ Contributions

Wan-Xia Yang and Chong-Ge You conceived and designed the study. Wan-Xia Yang and Yun-Yan Pan made the diagrams and tables of the article. Wan-Xia Yang and Chong-Ge You wrote the paper. All the authors read and approved the manuscript.

Acknowledgments

This study was supported by Cuiying Scientific and Technological Program of Lanzhou University Second Hospital (CY2018-MS10) and Gansu Province Science and Technology Plan Key Research and Development Project (18YF1FA108).