Identification and validation of real hub genes in hepatocellular carcinoma based on weighted gene co-expression network analysis

Abstract

BACKGROUND:

Hepatocellular carcinoma (HCC) is one of the most common liver malignancies in the world. With highly invasive biological characteristics and a lack of obvious clinical manifestations, hepatocellular carcinoma usually has a poor prognosis and ranks fourth in cancer mortality. The etiology and exact molecular mechanism of primary hepatocellular carcinoma are still unclear.

OBJECTIVE:

This work aims to help identify biomarkers of early HCC diagnosis or prognosis based on weighted gene co-expression network analysis (WGCNA).

METHODS:

Expression data and clinical information of HTSeq-Counts were downloaded from The Cancer Genome Atlas (TCGA) database, and gene expression map GSE121248 was downloaded from Gene Expression Omnibus (GEO). By differentially expressed genes (DEGs) and weighted gene co-expression network analysis (WGCNA) searched for modules in the two databases that had the same effect on the biological characteristics of HCC, and extracted the module genes with the highest positive correlation with HCC from two databases, and finally obtained overlapping genes. Then, we performed functional enrichment analysis on the overlapping genes to understand their potential biological functions. The top ten hub genes were screened according to MCC through the string database and Cytoscape software and then subjected to survival analysis.

RESULTS:

High expression of CDK1, CCNA2, CDC20, KIF11, DLGAP5, KIF20A, ASPM, CEP55, and TPX2 was associated with poorer overall survival (OS) of HCC patients. The DFS curve was plotted using the online website GEPIA2. Finally, based on the enrichment of these genes in the KEGG pathway, real hub genes were screened out, which were CDK1, CCNA2, and CDC20 respectively.

CONCLUSIONS:

High expression of these three genes was negatively correlated with survival time in HCC, and the expression of CDK1, CCNA2, and CDC20 were significantly higher in tumor tissues of HCC patients than in normal liver tissues as verified again by the HPA database. All in all, this provides a new feasible target for early and accurate diagnosis of HCC, clinical diagnosis, treatment, and prognosis.

1. Introduction

Hepatocellular carcinoma (HCC) is currently the fifth most prevalent and fourth most deadly disease in the world [1]. It increases year by year, and the incidence of males is much higher than that of females, and the incidence of males and females is 2.8:1 [2]. Globally, East Asia and Africa have the highest incidence of HCC, but the incidence and mortality of HCC are also increasing in Europe and the United States. Hepatitis B virus (HBV) infection accounts for 60% in Asia, and cirrhosis is the strongest risk factor [3]. Among the new cancer cases in China in 2020, liver cancer added about 302,598 among men, ranking fourth and accounting for 47.86% of the global incidence. Liver cancer added 107,440 among women, ranking 7th and accounting for 39.30% of the global incidence rate. In China, the death toll of liver cancer ranked second among men and fifth among women in 2020, accounting for 49.89% and 40.78% of the global death toll of liver cancer, respectively [4]. According to the SEER database of The National Cancer Institute, the average five-year survival rate of HCC patients in the United States is 19.6%, but for advanced and metastatic diseases, the survival rate may be as low as 2.5% [1]. Currently, there are surgical and non-surgical treatments for HCC. Among them, surgical treatment includes partial liver resection and liver transplantation, which are suitable for early patients. Non-surgical treatment includes tumor radiotherapy, chemotherapy, biotherapy, immunotherapy, etc., both surgical and non-surgical treatment have a significant effect on early patients, but the treatment effect on advanced patients is not good. Therefore, how to improve the early diagnosis rate of HCC patients is the focus of our attention. The current outlook for HCC has mentioned that screening with novel serum biomarkers may help in early patient diagnosis and more effective selection of appropriate treatments for the disease, thus effectively improving patient survival cycles, survival rates, etc. [3].

In recent years, with the rapid development of high-throughput measurement technologies, bioinformatics analysis of expression profiles has been widely used to find potential targets for cancer and biomarkers for patient prognostic analysis [5]. Differential gene expression analysis between different clinical phenotypes is a common method to search for potential biomarkers of diseases. Weighted gene co-expression network analysis allows the construction of co-expressed gene modules by finding coordinately expressed gene modules and analyzing the relationship between different modules and specific features [6].

In our study, overlapping genes were obtained by using a combination of differential gene analysis and WGCNA on HCC samples from the TCGA and GEO databases. The real hub genes were obtained by further analysis of biological functions and clinical characteristics of overlapping genes. The final results are helpful for us to further understand the mechanism of occurrence and development of HCC and provide new molecular targets for early detection and treatment of HCC.

2. Methods

The flow chart of this research is shown in Fig. 1.

Figure 1.

The flow chart of present study.

2.1 Data collection and data processing

The Cancer Genome Atlas (TCGA) [7] database downloaded 424 HTSeq-Counts expression data on HCC and clinical information (https://portal.gdc.cancer.gov/). They included 374 cancerous tissue samples and 50 solid tissue normal samples. The genetic transformation was performed using Perl language, and, 19600 genes were analyzed. Downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/) on human liver cancer cell gene expression profile, select sample than 50 sets of data, eventually, take GSE121248 [8] as the data source. It was provided by Hui KM based on the GPL570 platform ([HG-U133_PLus_2] Affymetrix Human Genome U133 Plus 2.0 Array) and included 37 para cancer tissue samples and 70 tumor tissue samples. Probes were converted to gene symbols according to the annotation file, and duplicate probes for the same genes were removed by determining the median expression values of all corresponding probes. The results resulted in a total of 21,654 genes screened for subsequent analysis.

2.2 Differentially expressed genes identification

Significance analysis of differential gene expressions (DEGs) between tumor tissue and normal tissue samples was performed using the Limma package (https://www.bioconductor.org/) in R language [9], with set criteria of $\mid$ logFC $\mid$ $⩾$ 1 and adjusted $p$ -value $<$ 0.05. Two R language packages, “pheatmap” and “ggplot2”, were used to create the heatmap and volcano map.

2.3 Weighted gene co-expression network analysis and overlapping genes construction

WGCNA is an analytical method for analyzing the gene expression profiles of multiple samples. It allows clustering of genes with similar expression patterns and analysis of the relationship between modules and specific traits or phenotypes. Distances between different transcripts were calculated using Pearson correlation coefficients. In our study, the “WGCNA” [6] software package of R language was used to implement weighted gene co-expression network analysis (WGCNA) and construct a weighted co-expression network for gene expression data of HCC. The weighted adjacency matrix is transformed into a topological overlap matrix (TOM) to detect network interactions and a hierarchical clustering approach is used to build a clustering tree structure for the TOM matrix. Eventually, genes with similar expression patterns are combined into the same module. The “VennDiagram” R language package [10] was used to draw the Venn diagram of TCGA and modules in GSE121248 that are positively correlated with tumor tissues.

Figure 2.

The volcano plot of DEGs between HCC and normal tissue. (A) Volcano plot of DEGs in the TCGA-LIHC dataset. (B) Volcano plot of DEGs in the GSE121248 dataset.

2.4 GO and KEGG function enrichment analysis of gene modules

GO enrichment analysis can effectively predict the biological characteristics of overlapping genes, and KEGG can effectively analyze the pathways of gene function. Functional annotation of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) based on R language package clusterProfiler [11] pathway enrichment analysis. Adjusted $p$ -value $<$ 0.05 was considered to have statistical significance.

2.5 Identification of hub genes in functional modules

The intersection of TCGA-LIHC and tumor tissue positively related modules in GSE121248 is obtained to obtain overlapping genes, which are imported into STRING 11.5 (https://cn.string-db.org/) [12] to construct Protein-Protein Interaction (PPI) network diagram. Set confidence score $⩾$ 0.90 and import the results into Cytoscape 3.9.0 (http://www.cytoscape.org) to build network visualization images. The CytoHubba plugin in Cytoscape was used to calculate each of its nodes, where the maximum cluster centrality (MCC) algorithm was considered to be the most efficient way to find the central node, and the top 10 genes with the largest weight coefficients were obtained for subsequent analysis [13].

2.6 Hub genes clinical trait validation and survival analysis

To better verify the relationship between hub genes and HCC patients, the Transcripts Per Million (TPM) was used as the standardized value for sample difference analysis, and the expression point map of hub genes in normal tissues, as well as tumor tissues, were plotted. To obtain the relationship between hub genes and clinical characteristics, a survival R language package was used to analyze the clinical information of HCC patients downloaded from TCGA with hub genes. The R package “Survival” and “SurvMiner” used in survival analysis need to perform the log-rank test and set $p$ -value $<$ 0.05. Kaplan-Meier curve was plotted by expression spectrum, and overall survival (OS) curves were plotted for HCC patients by dividing them into two groups based on the median expression values of the hub genes. The online tool Gene Expression Profilings Interactive Analysis (http://gepia.cancer-pku.cn/) [14] was used to plot disease-free survival (DFS). The genes related to the survival of HCC patients were obtained, and then compared with the genes contained in the top 10 pathways of KEGG enrichment analysis, Finally, 3 real hub genes were identified for subsequent analysis.

2.7 Validation of hub genes protein expression by the human protein atlas database

The Human Protein Atlas database (HPA) (https://www.proteinatlas.org/), which is based on proteomics, transcriptome, and systems biology data on tissue, cells, and organs such as mapping [15]. Not only tumor tissue but also normal tissue protein expression were included. We used the HPA database to validate the protein expression of real hub genes in normal and HCC tissues based on the immunohistochemical (IHC) method.

Figure 3.

The heatmaps of DEGs between HCC and normal tissues. (A) The heatmap of DEGs obtained from TCGA-LIHC dataset. (B) The heatmap of DEGs obtained from GSE121248 dataset.

Figure 4.

Weighted gene co-expression modules construction and selection from TCGA-LIHC. (A) The left graph shows the scale free fit index for soft-threshold powers. The right graph displays the mean connectivity for soft-threshold powers. (B) Hierarchical clustering analyses in TCGA-LIHC. (C) Dendrogram of all co-expression genes modules were clustered based on a dissimilarity measure (1-TOM).

Figure 5.

Weighted gene co-expression modules construction and selection from GES121248. (A) The left graph shows the scale free fit index for soft-threshold powers. The right graph displays the mean connectivity for soft-threshold powers. (B) Hierarchical clustering analyses in GES121248. (C) Dendrogram of all co-expression genes modules were clustered based on a dissimilarity measure (1-TOM).

Figure 6.

The heatmaps of correlation between clinical traits and the modular genes. (A) The heatmap of traits and modular genes constructed from TCGA-LIHC. (B) The heatmap of traits and modular genes constructed from GSE121248.

Figure 7.

MM and GS correlation analysis scatterplots. (A) Scatterplot of MM and GS correlation analysis of yellow module gene in TCGA-LIHC. (B) Scatterplot of MM and GS correlation analysis of turquoise module gene in GSE121248.

Figure 8.

The Venn diagram of DEGs and co-expression genes. As a result, 73 overlapping genes were obtained from the intersection of DEGs and co-expression.

Figure 9.

Functional enrichment analysis of the 73 overlapping genes. (A) GO enrichment analysis. (B) KEGG pathway enrichment analysis. The size of the bubble represents the number of enriched genes and the color of the bubble represents the $p$ -value.

Figure 10.

Construction of PPI network and screening of hub genes based on CytoHubba. (A) Construct overlapping genes PPI network based on STRING. (B) PPI network of overlapping genes was constructed based on Cytoscape. (C) The top 10 genes with the highest MCC score were screened from PPI.

3. Results

3.1 Identification of differentially expressed genes

Data information of HCC patients from different platforms was analyzed using the Limma software package and set to $\mid$ logFC $\mid$ $⩾$ 1 and adjusted $p$ -value $<$ 0.05. A total of 2703 DEGs were screened from HCC and normal control samples from The TCGA-LIHC database, among which 1023 genes were up-regulated and 1680 genes were down-regulated in tumor tissues (Fig. 2A). A total of 580 DEGs were screened from the GSE121248 dataset, of which 167 were up-regulated genes and 413 were down-regulated genes (Fig. 2B). In addition, heat maps of DEGs of two data sets were constructed using the R language, as shown in Fig. 3A and B.

3.2 Weighted gene co-expression network construction

The data sets of TCGA-LIHC and GSE121248 were used to construct the sample clustering tree of HCC and normal samples through WGCNA data packets, and the soft power threshold $\beta$ was selected as 3 and 10, respectively. To ensure that the fitting index (R2) and average connectivity of the scale-free topological model reaches a stable state, as shown in Figs 4A and 5A. Then, gene modules were detected based on the TOM matrix, and genes with similar expression patterns were placed into different modules through average linkage clustering. Hierarchical clustering analysis is shown in Figs 4B and 5B. Finally, after assigning colors to each module, 9 modules in TCGA-LIHC (Fig. 4C) and 8 modules in GSE121248 (Fig. 5C) were identified respectively.

3.3 Identification of DEGs overlapping genes in TCGA-LIHC and GSE121248 dataset

Subsequently, hierarchical cluster analysis was further evaluated to construct the correlation heat maps between different modules and clinical features. The results showed that the yellow color in TCGA-LIHC and the turquoise color in GSE121248 had the highest correlation with tumor tissues (Fig. 6A and B) (yellow module: $R=$ 0.47, $p=$ 5e-25; turquoise module: $r=$ 0.72, $P=$ 3e-18). The yellow module and the turquoise module selected above were subjected to the next step of the analysis, and the scatter plots of individual genes within the modules were obtained by calculating the gene significance (GS) and module membership (MM) analysis within the modules (Fig. 7A and B). Finally, the “VennDiagram” software package was used to intersect the DEGs of TCGA-LIHC and GSE121248 with the selected yellow and turquoise module genes to obtain the Venn diagram. For the final 2,703 DEGs from TCGA and 580 DEGs from GEO, 73 co-expressed genes were extracted from the yellow module of TCGA-LIHC and the turquoise module of GSE121248, respectively, with 1235 co-expressed genes and 734 co-expressed genes (Fig. 8).

3.4 GO and KEGG functional enrichment analysis of the 73 overlapping genes

The “clusterProfiler” R language package was used to enrich the functions of 73 overlapping genes, map their GO and KEGG pathway bubbles diagram (Fig. 9A and B), to predict the potential biological functions and related pathways of these 73 genes. GO enrichment analysis showed that these genes were significantly enriched in the nuclear division, organelle fission, spindle, chromosomal region, microtubule-binding, and tubulin binding ( $p<$ 0.05). In the KEGG enrichment pathway, these genes are significantly enriched in the cell cycle, DNA replication, cellular senescence, p53 signaling pathway, viral carcinogenesis and mismatch repair, and nucleotide excision repair ( $P<$ 0.05).

3.5 Identification of hub genes from overlapping genes

To further select hub genes in the overlapping genes obtained by WGCNA analysis, 73 genes were imported into the STRING 11.5 database in this study to construct a PPI network (Fig. 10A). The results were then imported into the visualization software Cytoscape (Fig. 10B), and the MCC algorithm of the CytoHubba plug-in was used to screen out the top 10 nodes (Fig. 10C). These 10 genes are hub genes, namely CDK1, TOP2A, DLGAP5, CDC20, ASPM, CEP55, TPX2, KIF11, KIF20A, and CCNA2.

3.6 Correlation analysis of hub genes expression and survival

To verify the authenticity of the 10 screened hub genes, we standardized the target genes based on TCGA-LIHC data set through the TPM method, and finally drew the expression point map of hub genes in normal tissues and liver cancer tissues (Fig. 11). Overall survival (OS) of these 10 genes was analyzed and Kaplan-Meier curves were drawn by the “Survival” and “SurvMiner” R language package (Fig. 12). Disease-free survival (DFS) curves of 10 genes were plotted using the online website GEPIA2 (Fig. 13), and their log-rank $P$ values were obtained.

3.7 Validation of the real hub genes based on HPA

By analyzing the difference in expression of hub genes in normal and hepatocellular carcinoma tissues, combined with the OS curve and DFS curve of hub genes. We found that the expression levels of CDK1, DLGAP5, CDC20, ASPM, CEP55, TPX2, KIF11, KIF20A, and CCNA2 were significantly correlated with OS and DFS in HCC patients ( $P<$ 0.05). Combined with the distribution of these 9 genes in KEGG enrichment pathways (Table 1), it can be concluded that the pathways enriched by CDK1, CDC20, and CCNA2 account for the highest proportion and are the real hub genes we need to find. It is mainly involved in the cell cycle, DNA replication, oocyte meiosis, and other related pathways. The immunohistochemical images of real hub genes were obtained from the HPA database (Fig. 14), and the expression of real hub genes in different tissue samples could be directly reflected by immunohistochemistry.

Table 1

The KEGG pathway of the hub genes

KEGG ID	Description	$P$ value	$P$ adjust	Gene ID
hsa04110	Cell cycle	2.45E-12	1.13E-10	BUB1B/PTTG1/CDK1/CCNB1/MCM3/TTK/
				CDC20/CCNA2/CDC7/CDC6
hsa03030	DNA replication	4.55E-06	1.05E-04	RFC4/MCM3/PRIM1/POLE2
hsa04114	Oocyte meiosis	7.37E-04	1.13E-02	PTTG1/CDK1/CCNB1/CDC20
hsa04218	Cellular senescence	1.41E-03	1.44E-02	CDK1/CCNB1/CCNA2/FOXM1
hsa04115	p53 signaling pathway	1.57E-03	1.44E-02	RRM2/CDK1/CCNB1
hsa03430	Mismatch repair	2.40E-03	1.84E-02	MSH2/RFC4
hsa04914	Progesterone-mediated oocyte maturation	4.07E-03	2.67E-02	CDK1/CCNB1/CCNA2
hsa05166	Human T-cell leukemia virus 1 infection	5.08E-03	2.92E-02	BUB1B/PTTG1/CDC20/CCNA2
hsa03420	Nucleotide excision repair	9.77E-03	4.68E-02	RFC4/POLE2
hsa04120	Ubiquitin mediated proteolysis	1.02E-02	4.68E-02	CDC20/UBE2S/UBE2C
hsa01524	Platinum drug resistance	2.26E-02	9.44E-02	TOP2A/MSH2
hsa05203	Viral carcinogenesis	2.66E-02	1.02E-01	CDK1/CDC20/CCNA2
hsa05170	Human immunodeficiency virus 1 infection	2.94E-02	1.04E-01	CDK1/CCNB1/APOBEC3B

Figure 11.

The expression level of hub genes in HCC and normal tissues was verified based on TCGA-LICH. (A) ASPM, (B) CCNA2, (C) CDC20, (D) CDK1, (E) CEP55, (F) DLGAP5, (G) KIF11, (H) KIF20A, (I) TOP2A, (J) TPX2.

Figure 12.

Validation of overall survival (OS) curve of hub genes. (A) ASPM, (B) CCNA2, (C) CDC20, (D) CDK1, (E) CEP55, (F) DLGAP5, (G) KIF11, (H) KIF20A, (I) TOP2A, (J) TPX2.

Figure 13.

Disease-free survival analysis of hub genes. (A) ASPM, (B) CCNA2, (C) CDC20, (D) CDK1, (E) CEP55, (F) DLGAP5, (G) KIF11, (H) KIF20A, (I) TOP2A, (J) TPX2.

Figure 14.

Immunohistochemistry of the real hub genes in tumor (top) and normal tissues (base) from the human protein atlas (HPA) database. (A) CDK1, (B) CDC20, (C) CCNA2. I represents liver tumor tissues. II represents normal tissues.

4. Discussion

In the current clinical diagnosis and treatment, the continuous increase of serum alpha-fetoprotein (AFP) is often used as an important indicator of liver cancer, but the increase of AFP is also common in some other liver diseases, such as cirrhosis [16], indicating that the specificity and sensitivity of AFP for the diagnosis of liver cancer is not enough. Therefore, biomarkers with higher specificity and sensitivity are needed to diagnose early HCC and effectively judge the future survival of HCC patients. At present, many new and more effective biomarkers are being discovered and applied in clinics. For example, silencing DUSP21 significantly inhibited cell proliferation, colony formation, and tumorigenicity of HCC cells in vivo. DUSP21 plays a significant role in maintaining hepatocellular carcinoma cell proliferation, therefore it could be a potential therapeutic target for liver cancer [17]. High expression of cholesterol esterase (SOAT1) was associated with the worst prognostic risk. By inhibiting the cholesterol esterification enzyme SOAT1, cholesterol levels on the cytoplasmic membrane can be reduced, effectively inhibiting tumor cell proliferation and migration [18]. Therefore, in our study, DEGs and WGCNA were used for the biological analysis of HCC in TCGA and GEO databases. As a result, 73 genes were obtained; PPI and Cytoscape were used for further analysis of overlapping genes to obtain 10 hub genes. Namely, CDK1, TOP2A, DLGAP5, CDC20, ASPM, CEP55, TPX2, KIF11, KIF20A, CCNA2. Combining the enrichment of OS, DFS, and KEGG pathways of 10 genes, ultimately we found that high expression of CDK1, CDC20, and CCNA2 was significantly associated with OS in HCC patients.

Cyclin Dependent Kinase 1 (CDK1) is the protein-encoding gene for a protein that is a member of the Serine/threonine-protein kinase family. The protein is a catalytic subunit of a highly conserved protein kinase complex called the M-phase promoting factor (MPF), which is required for the G1/S and G2/M phase transitions of the eukaryotic cell cycle. It was found that miR-582-5p was down-regulated in hepatocellular carcinoma tissues and up-regulation of miR-582-5p inhibited cell proliferation and blocked the cell cycle in G0/G1 phase. miR-582-5p directly inhibited the expression of CDK1 and AKT3, and indirectly inhibited the expression of cyclinD1, thereby regulating HCC [19]. Experiments have also confirmed that dihydroartemisinin (DHA) can reduce the expression level of CDK1, thereby inhibiting the proliferation of liver cancer cells [20]. At present, CDK1 is highly expressed in various tumor tissues to varying degrees, and it has been proven that CDK1 can play an important role in lung cancer [21], breast cancer [22], pancreatic cancer [23], colorectal cancer [24] and so on by regulating the cell cycle. In the present study, high CDK1 expression was strongly associated with poor prognosis in HCC patients. KEGG pathway also showed that CDK1 was mainly involved in cell cycle, oocyte meiosis, cellular senescence, P53 signaling pathway, and other signaling pathways. Studies have found that CCNB1/CDK1 can mediate mitochondria, coordinate ATP output, and provide biological energy for cell G2/M transformation, thus interrupting the cell replication cycle [25, 26]. Therefore, CDK1 plays an essential role in the cell cycle and cell proliferation.

Anaphase Promoting Complex (APC) was activated by cell division cycle 20 (CDC20) and CDH1 to form two different E3 ubiquitin ligase complex sub-complexes, APC ${}^{\text{CDC20}}$ and APC ${}^{\text{CDH1}}$ . CDC20 is an oncogenic factor and CDH1 is a tumor suppressor factor [27]. It is currently believed that the APC ${}^{\text{CDC20}}$ complex can disrupt the ubiquitination of downstream cell cycle regulators securin and cyclin B [28] and regulate the activity of downstream pluripotent-related transcription factors SOX2 to promote the invasion and renewal of glioma stem cells [29]. The degree of CDC20 expression has now been found to correlate significantly with tumor prognosis, with higher CDC20 expression associated with lower 5-year survival rates for patients. Studies have shown that the expression level of CDC20 increases in non-small cell lung cancer (NSCLC), resulting in the more significant pleural invasion of cancer cells, while the growth rate of lung cancer cells is significantly slowed down after the expression level of CDC20 is down-regulated [30]. In our study, the high expression of CDC20 also predicted the poor prognosis of HCC patients, and CDC20 was mainly enriched in the cell cycle, oocyte meiosis, and other pathways. A related study found that CDC20 was overexpressed in 68% of hepatocellular carcinoma tissues [31], that overall survival was lower in HCC patients with high CDC20 expression, and that CDC20 was predictive of prognosis in HCC patients [32]. In addition, numerous studies have found that anti-mitotic agents such as paclitaxel can inhibit APC ${}^{\text{CDC20}}$ , and CDC20 gene knockout can lead to a checkpoint-independent mitotic block, which can effectively kill cancer cells [33, 34]. All of these suggest that CDC20 is a potential target for the treatment of HCC.

The CCNA2 gene encodes a protein that belongs to a family of highly conserved cell cycle proteins, whose members function to regulate the cell cycle. CCNA2 binds and activates cyclin-dependent kinase 2, thereby facilitating the G1/S and G2/M transitions. Studies have shown that p53-P21-DREAM-CDE/CHR pathway can down-regulate the CCNA2 gene, thereby blocking the G2/M cell cycle and controlling the growth of cancer cells [35]. This article also found that high CCNA2 expression was significantly associated with the prognosis of HCC patients and that the proliferation of hepatocellular carcinoma cells could be inhibited by drugs that inhibit the FXR-miR-22-CCNA2 signaling pathway [36]. The current studies have found that CCNA2 is highly expressed in all types of tumors. For example, E2F1 can promote the proliferation of triple-negative breast cancer (TNBC) through up-regulation of CCNA2 [37]. miR-545 can inhibit proliferation by inhibiting CCNA2 expression in osteosarcoma cells [38]. In addition, the high expression of CCNA2 has been reported in lung cancer [39], gastric cancer [40], ovarian cancer [41], and esophageal cancer [42].

Studies have begun to target CDK1, CDC20, and CCNA2 as potential targets for liver cancer treatment. The therapeutic effect of sorafenib in hepatocellular carcinoma can be improved by blocking the CDK1/PDK1/DK1/ $\beta$ -Catenin signaling pathway through the CDK1 inhibitor RO3306. It was found that the combination of different doses of RO3306 and sorafenib in hepatocellular carcinoma cells showed synergistic effects and significantly inhibited tumor self-renewal. Low-dose combinatorial inhibition decreased S-phase and increased cell arrest in the G2/M phase. Both the tumor incidence and the mean tumor volume and weight were significantly lower in the combination group than in the control group. This may be an important personalized treatment option for patients with high CDK1 expression [43]. Meanwhile, it was found that the growth of liver cancer cells could be inhibited by inhibition of CDC20 which resulted in G2/M phase arrest [44]. And SIRT2 regulated APC E3 ligase activity partly through deacetylating CDC20, and Sirt2-deficient male mice developed hepatocellular carcinoma [45]. There are already different drugs that act on CDC20, such as Apcin, which is a specific inhibitor that acts directly on the APC ${}^{\text{CDC20}}$ complex, as well as natural drugs such as withaferin, GDNT, and genistein, which inhibit CDC20 and may be a safer way to better treat human cancers [46]. Arylbenzofuran (Ary) is a novel natural product that induces apoptosis in cancer cells via the mitochondrial membrane potential and acts to inhibit the growth of cancer cells by activating ERK1/2 to localize the nucleolus and upregulating the expression of cyclin A2 and Cdk2 to block the cells in the G1/S phase [47]. From the current study, CDK1, CDC20, and CCNA2 have been found to be effective in cancer treatment. Therefore, we have more reasons to believe that these three hub genes are the future research direction for the treatment of HCC.

Multiple bioinformatics methods were used to jointly screen out potential targets for HCC. This study indicated that CDK1, CDC20, and CCNA2 show high expression in hepatocellular carcinoma tissues and that their high expression is negatively correlated with the survival rate of HCC patients, and these targets are mainly involved in regulating the cell cycle and thus controlling the proliferation of tumor cells. Therefore, in the future, we may be able to use them in the clinical diagnosis of early-stage liver cancer and also to predict the survival time of patients with liver cancer. The existing studies corroborate our view. Studies on different hepatocellular carcinoma samples have found that CDK1 is highly expressed in hepatocellular carcinoma cells and can influence tumorigenesis and proliferation by regulating the cell cycle [48, 49]. Other studies have also found that Tax protein can alter mitotic progression and mitotic defects by activating the CDC20 regulatory complex APC ${}^{\text{CDC20}}$ [49, 50, 51]. This further confirms the high reliability of our screened genes and lays the foundation for further studies. These newly discovered potential targets may provide new ideas and new research directions for the treatment of HCC and have important clinical significance.

However, our study also had many limitations. First, the sample size we selected was not large enough for a specific analysis of subtypes in HCC, which may lead to inaccurate hub gene screening and insufficient coverage. Secondly, the reliability of hub genes has not been verified experimentally in this paper, and the specific mechanism of hub genes regulating HCC is not fully understood, which needs further exploration in the future. Finally, the correlation between the high expression of these targets and the clinical manifestations and staging of hepatocellular carcinoma was not thoroughly clarified, which may be directly related to the insufficient size of our sample data. In the future direction of our study, we will further incorporate more diverse data to make our analysis more comprehensive.

In the future, we will include a more diverse and richer sample of liver cancer and will analyze it for different stages of liver cancer and those caused by different etiologies. These will provide further support for early liver cancer diagnosis and targeted treatment of liver cancer.

5. Conclusion

In summary, our combined analysis of the WGCNA and TCGA, and GEO databases has uncovered CDK1, CDC20, and CCNA2 as three potential targets that may be able to diagnose, analyze prognosis and even treat liver cancer. We also further validated the expression of real hub genes in HCC tissues and normal tissues through tissue expression differences and immunohistochemistry. We also analyzed that CDK1, CDC20, and CCNA2 control tumor growth mainly by regulating the cell cycle, but we lack further experimental validation. These findings may provide ideas for our subsequent studies and provide a reference for the more complex study of hepatocellular carcinoma to investigate this complex mechanism of hepatocellular carcinoma treatment.