For screening DEGs, three high-throughput RNA-seq-based raw datasets, namely GSE25599 (Huang et al., 2011), GSE94660 (Yoo et al., 2017), and GSE77509 (Yang et al., 2017), which comprised patients with HBV infection, were downloaded from the NCBI GEO database¹. GSE94660 (21 paired normal and HBV-related HCC tissue samples) and GSE77509 (16 paired normal and HBV-related HCC tissue samples) were established using GPL16791 Illumina HiSeq 2500 (Homo sapiens), while GSE25599 (10 paired normal and HBV-related HCC tissue samples) was established using the GPL9052 Illumina Genome Analyzer. RNA-seq raw data and clinical data of 50 normal samples and 374 HCC samples² were downloaded from TCGA³. For the validation of independent prognostic genes, a dataset named LIRI-JP⁴, including 202 normal and 243 HCC tissue samples, was downloaded from the International Cancer Genome Consortium (ICGC)⁵. For comparing the expression levels of independent prognostic genes between HCV- and HBV-related HCC, the GSE69715 dataset [66 normal and 37 HCV-related HCC tissue samples established using GLP570 (HG-U133_Plus_2)] was downloaded from the GEO database.

Gene expression profile matrix files of GSE25599, GSE77509, and GSE94660 were obtained from raw datasets using Perl (de Hoon et al., 2004). Nevertheless, the gene expression profile matrix data of TCGA was acquired using the “TCGAbiolinks” R package (Colaprico et al., 2016) and Perl. Genes that were differentially expressed between normal and HBV-related HCC tissue samples were screened by the limma R (Ritchie et al., 2015) and edgeR R packages (Robinson et al., 2010). | Log₂(FC)| ≥ 1.0, p-value ≤ 0.05, and FDR ≤ 0.05 were set as the cutoff criteria for DEGs screening after background correction and data normalization. Overlapped DEGs among GSE25599, GSE77509, GSE94660, and TCGA were identified using the VennDiagram R package (Chen and Boutros, 2011). The heatmaps of DEGs, which could be divided into up- and down-regulated groups, were drawn using the “pheatmap” R package (Galili et al., 2018).

The names of DEGs were translated into gene IDs using the R programming language. To investigate the biological pathways that might be involved in the occurrence and development of HBV infection and HCC, candidate DEGs were segregated into up- and down-regulated groups and subjected to pathway enrichment analysis. Gene Ontology (GO) analysis, which involved three categories, namely molecular functions (MF), cellular components (CC), and biological processes (BP), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed with the threshold of FDR-value < 0.05 using the clusterProfiler R package (Yu et al., 2012), which facilitated biological terminology classification and gene cluster enrichment.

A protein–protein interaction (PPI) network of DEGs was constructed using STRING⁶ (Szklarczyk et al., 2019) and visualized with Cytoscape v3.6.1 (Shannon et al., 2003). DEGs that consisted of several important nodes with many other interaction partners were analyzed using Molecular Complex Detection (Bader and Hogue, 2003) and CytoHubba (Chin et al., 2014; Tang et al., 2019). Subnets of the vast protein interaction network were extracted by calculating the degree of nodes, and highly connected nodes with a degree score of >45 and p-value < 0.05 were identified as hub genes.

Survival analysis were primarily performed using clinical data from TCGA to predict the prognostic value of hub genes. Kaplan–Meier survival curves of hub genes were plotted using the survival R package⁷ and differences in survival rate were evaluated with a log-rank test threshold of p-value < 0.05. To evaluate the accuracy of the survival curves, receiver operating characteristic (ROC) curves were then constructed using the “survival ROC” R package (Huang et al., 2017) with the threshold of AUC ≥ 0.6. Next, Cox proportional-hazards models were used to estimate the effects of prognostic factors on survival using the survival R and “survminer” R packages⁸ with the threshold of p-value < 0.05. Univariate Cox analysis was first performed to screen for genes significantly associated with overall survival rate, and multivariate Cox analysis was then performed to identify independent prognostic genes (Orimo et al., 2008; Uemura et al., 2009).

To validate independent prognostic genes for HCC screened by survival analyses, the aforementioned TCGA clinical data were used to analyze individual gene expression levels between normal and HCC tissue specimens at different stages of tumor progression using the “ggpubr” R package⁹. Data pertaining to normal and HCC tissue samples were compared using Wilcoxon test, and those pertinent to multiple samples from different stages of tumor progression were compared using the Kruskal–Wallis test, with the threshold of p-value < 0.05.

To further evaluate the clinical value of the independent prognostic genes, a dataset of patients with HCC was downloaded from the ICGC portal (see text footnote 4) for survival and ROC curve analyses; for this purpose, we used the survival R package, survival ROC R package, and Perl. Ultimately, a meta-analysis of the independent prognostic genes in Oncomine¹⁰ (Rhodes et al., 2004), a cancer-profiling database containing published data and listing differential gene expression analyses, were performed to verify their expression levels in patients with HCC using four published data (Chen et al., 2002; Wurmbach et al., 2007; Roessler et al., 2010).

To analyze the potential relationships among the four independent prognostic genes in HCC occurrence and development, TCGA dataset was subjected to correlation analyses using the corrplot¹¹ R software. The correlation coefficient (Cor), ranging from −1 (perfect negative correlation) to +1 (perfect positive correlation), indicated how closely data in a scatterplot were arranged along a straight line. p-value < 0.05 for the coefficients indicates a statistically significant relationship.

Since viral Hepatitis B and Hepatitis C are the most commonly implicated risk factors for HCC, to compare the expression levels of the independent prognostic genes between HCV-related HCC and HBV-related HCC, the GSE69715 dataset for HCV and the GSE94660 dataset for HBV were analyzed using gglpot2¹², cowplot¹³, and ggpubr¹⁴ package. The relative expression levels (i.e., fold change) of these four genes in the tumor tissues comparing with normal tissues were calculated. Wilcoxon test was carried out between the HCC and normal tissues. p-value < 0.05 indicate statistical significance.

Differentially expressed genes were identified from three raw datasets, namely GSE25599, GSE94660, and GSE77509, downloaded from the NCBI GEO database and one downloaded dataset from TCGA database using the limma R package and edgeR R package. The cutoff criteria were | log₂(FC)| ≥ 1.0, p-value ≤ 0.05 and FDR ≤ 0.05. In total, 230 DEGs were overlapping among the four datasets, of which 103 were up-regulated (Figure 1A) and 127 were down-regulated (Figure 1B). Cluster heatmaps showing the expression levels of the 230 DEGs in each of the four datasets were generated (Figures 1C–F). Details of the top 20 up- and 20 down-regulated DEGs in HBV-related HCC are shown in Supplementary Table S1.

To investigate the functional annotation of DEGs, GO, and KEGG pathway enrichment analyses were performed. The results were considered to be statistically significant if FDR value was <0.05. The top 15 GO terms of up-regulated genes are listed in Supplementary Table S2. As evident from Figure 2A and Supplementary Table S2, in the MF, CC, and BP categories, the up-regulated genes were significantly enriched in nuclear division, organelle fission, and mitotic nuclear division; spindle, chromosomal region, and spindle pole; and protein kinase binding, enzyme binding, and chromatin binding, respectively. Further, KEGG pathway analysis of the up-regulated genes indicated that they were primarily enriched in cell cycle, p53 signaling pathway, and oocyte meiosis (Figure 2C and Supplementary Table S3).

As evident from Figure 2B,D and Supplementary Tables S4, S5, in the MF, CC, and BP categories, the down-regulated genes were mainly involved in small molecule catabolic process, organic acid catabolic process, carboxylic acid catabolic process, extracellular matrix, collagen-containing extracellular matrix, collagen trimer, cofactor binding, iron ion binding, and monooxygenase activity, respectively. Moreover, KEGG pathway analysis of the down-regulated genes indicated that they were enriched in retinol metabolism, arachidonic acid metabolism, and drug metabolism–cytochrome P450.

A PPI network of DEGs (Figure 3A) containing 230 nodes and 1189 edges was constructed by STRING and visualized by Cytoscape, which provides critical assessment and integration of protein-protein interactions, including direct (physical) and indirect (functional) correlations. Pivotal modules of the network were obtained using Molecular Complex Detection, and the degree of nodes was calculated using CytoHubba. In the PPI network, the number of edges involved determines the degree score of nodes; the nodes with high degree scores were considered to be hub genes (Chin et al., 2014). 54 DEGs with a degree score of >10 and p-value < 0.05 are listed in Supplementary Figure S1. There were 39 genes with degree scores of >30, 24 genes with degree scores of >40, 12 genes with degree scores of >45, and only one gene with degree score of >50, and all of these genes meet the requirements of p-value < 0.05. The modules with 39 nodes and 698 edges (degree score >30 and p-value < 0.05) were extracted to construct a subnet (Figure 3B). The most significant modules of 12 genes (degree score >45 and p-value < 0.05; Figure 3C) were identified as hub genes. The names, abbreviations, and scores of hub genes are summarized in Supplementary Table S6. The top five hub genes with the highest interaction node degrees were CDK1, CCNB1, CCNA2, BUB1B, and CCNB2, implying their potential roles in the development of HBV-related HCC.

Gene Ontology and KEGG pathway enrichment analyses were utilized to investigate the functional enrichment of the 12 hub genes. In the MF, CC, and BP categories, these 12 hub genes were mainly enriched in nuclear division, organelle fission, mitotic nuclear division, spindle, spindle pole, condensed chromosome, protein kinase binding, protein C-terminus binding, and protein serine/threonine kinase activity, respectively (Figure 3D and Supplementary Table S7). Further, KEGG pathway analysis for the hub genes (Figure 3E and Supplementary Table S8) indicated that they were primarily enriched in cell cycle, progesterone-mediated oocyte maturation, and oocyte meiosis.

It is noteworthy that all the 12 hub genes were up-regulated in patients with HBV-related HCC. To explore their prognostic importance, all of them were evaluated using the Kaplan–Meier plot and ROC curve with clinical and expression data from TCGA. As shown in Figure 4, based on their expression levels, AUC values of the 12 hub genes (BUB1B, CCNA2, CCNB1, CCNB2, CDC20, CDK1, KIF11, KIF20A, MAD2L1, PLK1, TOP2A, and TPX2), ranged from 0.6 to 0.7, while log-rank test showed p-value < 0.05 in all of the survival curves. Therefore, we considered that all the 12 hub genes appeared to be capable of survival prediction with the thresholds of p-value < 0.05 and AUC ≥ 0.6. Patients with HCC and up-regulation of these genes showed worse survival rate.

Further, univariate and multivariate Cox regression analyses were performed to analyze their independent prognostic importance in patients with HCC. As indicated in Table 1, univariate Cox regression analysis showed that all the 12 hub genes were high-risk genes (hazard ratio >1, p-value < 0.05); however, multivariate Cox regression analysis suggested that only KIF11, TPX2, KIF20A, and CCNB2 were independent prognostic genes in case of patients with HCC (hazard ratio >1, p-value < 0.05).

The Cancer Genome Atlas dataset of normal and HCC tissue samples were subjected to Wilcoxon test; KIF11, TPX2, KIF20A, and CCNB2 were found to have higher mRNA expression levels in HCC than in normal tissue samples (Figures 5A–D). Furthermore, the expression levels of KIF11, TPX2, KIF20A, and CCNB2 in multiple samples from different stages (I–IV) of tumor progression were compared using the Kruskal–Wallis test, and the results revealed that in comparison with normal tissue samples, the expression levels of KIF11, TPX2, KIF20A, and CCNB2 were higher at each stage of HCC (Figures 5E–H). These findings indicated the potential roles of these genes for diagnostic and prognostic prediction.

Another dataset with RNA-Seq mRNA expression data and clinical pathological data were obtained from the ICGC portal as an independent validation cohort to verify the prognostic potential of KIF11, TPX2, KIF20A, and CCNB2 in HBV-related HCC. Overall survival rate analysis of these four genes was performed using the Kaplan–Meier plot and ROC curves. The results were consistent with those obtained from TCGA datasets, revealing that patients with up-regulated KIF11, TPX2, KIF20A, and CCNB2 genes showed worse survival rate (Figures 6A–D, p-value < 0.05 and AUC ≥ 0.6). Notably, as indicated in Figure 6, AUC values calculated using ICGC data were a little higher (from 0.7 to 0.8) compared with those using the TCGA data. In general, an AUC of 0.7 to 0.8 is considered to be acceptable (Mandrekar, 2010). In addition, meta-analysis in Oncomine showed that KIF11, TPX2, KIF20A, and CCNB2 were highly expressed in HCC comparing with normal tissues samples (Figure 6E–H, p-value < 0.05). Correlation analyses (Supplementary Figure S2) revealed the potential relationships among these four independent prognostic genes, implying that these four genes have combined effects in HCC occurrence and development.

To deternmine whether KIF11, TPX2, KIF20A, and CCNB2 are specific to HBV-induced HCC comparing with HCV-induced HCC, HCV-related HCC (GSE69715) and HBV-related HCC (GSE94660) datasets were used to analyze the relative expression levels (i.e., fold change) of these 4 genes in HCC and normal tissue samples. As indicated in Figure 7 and Table 2, log₂(FC) of KIF11, TPX2, KIF20A, and CCNB2 in the HBV-related HCC dataset ranges from 1.38 to 2.66, while log₂(FC) of these genes in the HCV-related HCC dataset ranges from 0.09 to 0.35. Although all of the p-values in the Figure 7 are less than 0.05 (p-value < 0.05), which are statistically significant, we do not believe that they are biologically significant, because the expression level of these four genes in HCV-related HCC showed only a minor increase as compared with that in HBV-related HCC.