First, RNA-Seq and corresponding clinical data for liver hepatocellular carcinoma (LIHC) were obtained from the TCGA database (https://portal.gdc.cancer.gov/). A list of 424 samples was obtained, including 374 LIHC and 50 normal liver tissues, and RNA-seq count data on 19645 genes were obtained. The Illumina HiSeq 2000 platform was used to generate and annotate all data to a reference transcript set of the human hg38 gene standard track. The edgeR package tutorial suggested that genes with low read counts do not merit further analysis (11). Hence, genes with a count per million (CPM) <1 were omitted from this analysis. Next, the function rpkm in the edgeR package was adapted for further filtering. Consequently, 13,924 genes were acquired for subsequent analysis. Second, the RNA-Seq data and clinical data of HCC patients were acquired from the ICGC database (https://dcc.icgc.org/). A total of 260 HCC samples, which mainly originated from the Japanese population with HBV or HCV infection, were acquired (12). We chose the normalized read count values of the ICGC-LIRI-JP cohort. As a result, 22,913 genes were obtained for the next analysis.

The gene coexpression network of the TCGA-LIHC dataset was built through the WGCNA package (8). To build a scale-free network, a soft-power β = 7 ( Figures 2A, B ) was used in the TCGA-LIHC dataset. Next, the adjacency matrix was created according to the formula aij = |Sij|^β (aij: adjacency matrix between gene i and gene j, Sij: similarity matrix made by Pearson’s correlation coefficient of all gene pairs, as well as β: soft-power value). Subsequently, we converted this matrix into a topological overlap matrix (TOM) and the corresponding dissimilarity (1-TOM). The hierarchical clustering dendrogram of the 1-TOM matrix was established to aggregate the genes with similar expression patterns into the same coexpression module. Afterward, the module-trait relations between modules and external traits were analyzed to identify functional modules from the coexpression network. Hence, the modules with the largest correlation coefficients were regarded as modules that highly correlated with clinical traits. We chose the module that was positively associated with LIHC for our subsequent analysis.

The limma package is often used to perform differential gene expression analysis of gene expression profiles and RNA-Seq datasets (13). Here, we applied the limma package in the differential expression analysis of the TCGA-LIHC dataset to identify differentially expressed genes (DEGs) between LIHC and nontumorous tissues. To minimize the false discovery rate (FDR) to the greatest extent possible, we adjusted the P-value with the Benjamini–Hochberg (BH) method. The filtering criteria for DEGs were |logFC|>1 and adj. P <0.05. Afterward, we took the intersection of genes between DEGs and coexpressed genes to improve the reliability of screening closely related genes, and these differentially coexpressed genes were used for the next analysis.

The PPI network of differentially coexpressed genes was built through the Search Tool for the Retrieval of Interacting Genes (STRING) database (14). Then, we established a visual network of molecular interactions with combined scores ≥0.7 using Cytoscape (15). In addition, the degree values of all nodes in the PPI network were calculated using the CytoHubba plugin (16). The top 40 nodes with the highest degree scores were selected and regarded as hub genes associated with LIHC. The forty hub genes related to LIHC were displayed using the CytoHubba plug-in. In addition, we conducted gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the 40 hub genes to explore their biological functions. Adj. P values <0.05 were considered significant.

To analyze the prognostic roles of the differentially coexpressed hub genes in LIHC, we performed univariate Cox regression analysis of OS using the survival package based on the TCGA-LIHC dataset. LIHC patients without follow-up information or a survival time=0 days were excluded from our analysis, and the other patients in the TCGA-LIHC dataset were classified into two groups considering the median expression levels of the differentially coexpressed hub genes. Log-rank P<0.01 was considered significant. Additionally, the correlation network of these differentially coexpressed hub genes was established through the igraph package. The filtering criterion was a cutoff >0.75.

To decrease the risk of overfitting to the greatest extent possible, we used least absolute shrinkage and selection operator (LASSO) Cox regression analysis to build the prognostic module of LIHC (17, 18). The LASSO algorithm is widely utilized to select and shrink variables using the glmnet package. We used the expression matrix of the differentially coexpressed hub genes with prognostic value as the independent variable, while the OS and status of patients in the TCGA-LIHC dataset were used as the response variables. Then, we determined the penalty parameter (λ) of this module using tenfold cross-validation following the minimum criteria, namely, the λ value corresponding to the minimum partial likelihood deviance.

We calculated the risk scores of all LIHC patients using the expression level of every gene and the corresponding regression coefficient. The following formula was used: score= e^{sum (every gene’s expression × corresponding coefficient)}. Then, LIHC patients were divided into high- and low-risk cohorts based on the median value of the risk score. Subsequently, we constructed a nomogram of the prognostic signature to predict the survival of LIHC patients. Furthermore, we built calibration curves and time-dependent receiver operating characteristic (ROC) curves to evaluate the discrimination and accuracy of the prognostic multigene signature. The GSE112790 dataset was used to validate the expression patterns of the genes in the signature between LIHC and nontumorous tissues.

To analyze the prognostic value of the gene signature, we performed Kaplan-Meier survival analysis between the low- and high-risk groups using the survminer package based on the TCGA-LIHC and ICGC-LIRI-JP datasets. Additionally, to explore distribution in the low- and high-risk cohorts, we performed principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) on the TCGA-LIHC and ICGC-LIRI-JP datasets using the stats and Rtsne packages, respectively. To determine whether the risk score acts as an independent indicator of the prognosis of LIHC patients, we performed univariate and multivariate Cox regression analyses among all available variables using the TCGA-LIHC and ICGC-LIRI-JP datasets.

To acquire the DEGs between the low- and high-risk groups, we performed differential gene expression analysis using the limma package in the TCGA-LIHC and ICGC-LIRI-JP datasets. The P-value was adjusted using the BH method. The filtering criteria for DEGs were |logFC|>2 and adj. P <0.05. Afterward, we conducted GO and KEGG pathway analyses of the DEGs between the low- and high-risk groups in the TCGA-LIHC and ICGC-LIRI-JP datasets. To further analyze the relationship between the risk score and immune status, we calculated the infiltrating scores of 16 immune cells and 13 immune-related functions or pathways using single-sample gene set enrichment analysis (ssGSEA) (19).

To find the pivotal module in LIHC, the gene coexpression network was established in the TCGA-LIHC dataset. A list of 11 modules was generated ( Figure 2C ). Next, the heatmap revealed the correlations between the modules and clinical traits (normal and LIHC) in the TCGA-LIHC dataset ( Figure 2D ). Furthermore, the yellow module of the TCGA-LIHC dataset positively correlated with LIHC tissues (r=0.57, P=1e-37) and was used for our next analysis.

The heatmap displayed the expression patterns of fifty upregulated and fifty downregulated genes in the TCGA-LIHC dataset ( Figure 3A ). The volcano plot indicated that 2708 DEGs had a conspicuous dysregulation between LIHC and nontumorous tissues in the TCGA-LIHC dataset ( Figure 3B ). The Venn diagram showed the intersection of coexpressed genes ( Table S1 ) and DEGs ( Table S2 ); namely, 393 differentially coexpressed genes were identified ( Figure 3C ).

Figure 4A displays the PPI network of the differentially coexpressed genes with 241 nodes and 4792 edges. Subsequently, we quantified the degree scores of all nodes in this PPI network through the CytoHubba plugin ( Table S3 ) and chose the top 40 nodes as hub genes that are closely correlated with LIHC ( Figure 4B ). In addition, GO analysis showed significant enrichment in the mitotic nuclear division, organelle fission and spindles terms ( Figure S1A ). KEGG pathway analysis showed enrichment in the cell cycle and oocyte meiosis pathways ( Figure S1B ).

Univariate Cox regression analysis of the differentially coexpressed hub genes demonstrated that 38 hub genes were closely associated with the survival of LIHC patients ( Figure 5A ). The heatmap revealed that the 38 hub genes with prognostic value were significantly overexpressed in LIHC tissues ( Figure 5B ). Additionally, the correlation network suggested that the differentially coexpressed hub genes closely interact with each other ( Figure 5C ).

We used the LASSO Cox regression module to build a prognostic signature based on the expression matrix of the 38 differentially coexpressed hub genes. Consequently, we identified a four-gene signature module according to the optimal λ value ( Figures 6A, B ). In addition, we calculated the risk scores of LIHC patients using the following formula: score= e ^{(0.225*expression value of CDCA8+0.124*expression value of KIF20A+0.012*expression value of KIF2C+0.144*expression value of CEP55)}. Then, we established a nomogram to predict the 1-, 2-, and 3-year OS probability of LIHC patients ( Figure 6C ). The calibration curves of 1-, 2-, and 3-year OS probability showed satisfactory calibration of this nomogram ( Figures 6D–F ). Moreover, based on GSE112790, we confirmed that CDCA8, KIF20A, KIF2C and CEP55 were significantly overexpressed in LIHC tissues compared with nontumorous tissues ( Figure S2 ). Furthermore, the ROC curves suggested acceptable accuracy of this nomogram (area under the curve [AUC] of 1-year survival: 0.739; AUC of 2-year survival: 0.714; and AUC of 3-year survival: 0.673) ( Figure 7A ). Afterward, all LIHC patients were divided into a low-risk cohort (n=183) and a high-risk cohort (n=182) based on the median risk score ( Figure 7B ). The high-risk cohort in the TCGA-LIHC dataset had more deaths ( Figure 7C ), a poorer tumor grade, a higher clinical stage and a higher T stage ( Table 1 ). Consistently, Kaplan-Meier survival analysis showed that LIHC patients in the high-risk cohort experienced shorter survival than those in the low-risk cohort ( Figure 7D , P=1.14e-4). In the PCA of the TCGA-LIHC dataset, the first principal component (PC1) could explain 88.6% of total variance, and the PC1 scores were negatively correlated with the risk scores of patients ( Figure 7E ), while the second principal component (PC2) could explain 5.4% total variance ( Figure S3A ). Moreover, PCA and t-SNE analysis revealed that most LIHC patients in the high- and low-risk cohorts were distributed in two different directions ( Figure 7F ).

To validate the robustness of the four-gene signature module from the TCGA-LIHC dataset, we chose the ICGC-LIRI-JP dataset for further verification. First, we stratified LIHC patients from the ICGC-LIRI-JP dataset into high-risk and low-risk cohorts according to the median value of the risk score, which was calculated using the formula mentioned above. Consistent with the outcomes from the TCGA-LIHC dataset, the four-gene signature had an excellent AUC ( Figure 8A , 1-year survival: 0.752; 2-year survival: 0.751; and 3-year survival: 0.782). Moreover, the high-risk group correlated with a higher rate of mortality ( Figures 8B, C ). Additionally, patients from the high-risk cohort experienced significantly shorter survival than those in the low-risk cohort ( Figure 8D , P=1.24e-3). In the PCA of the ICGC-LIRI-JP dataset, the PC1 could explain 79% of total variance, and the PC1 scores were positively correlated with the risk scores of patients ( Figure 8E ), whereas the PC2 could explain 13% total variance ( Figure S3B ). In addition, t-SNE analysis validated that most patients in the high- and low-risk cohorts were distributed in two different directions ( Figure 8F ). In general, these outcomes in the ICGC-LIRI-JP dataset were similar to those in the TCGA-LIHC dataset.

To determine whether the risk score plays an independent prognostic role, we performed univariate and multivariate Cox regression analyses of the survival of LIHC patients. The univariate Cox regression analysis indicated that a higher risk score was closely correlated with worse survival in LIHC patients using the TCGA-LIHC ( Figure 9A , hazard ratio [HR]=3.324, 95% confidence interval [CI]: 2.181–5.066, P<0.001) and ICGC-LIRI-JP ( Figure 9B , HR=1.413, 95% CI: 1.243–1.607, P<0.001) datasets. Similar to the results of the univariate Cox regression analysis, the multivariate Cox regression analysis still suggested the risk score as an independent indicator for the survival of LIHC patients using the TCGA-LIHC ( Figure 9C , HR=3.041, 95% CI: 1.930–4.790, P<0.001) and ICGC-LIRI-JP ( Figure 9D , HR=1.378, 95% CI: 1.210–1.569, P<0.001) datasets.

Differential gene expression analyses were conducted in the TCGA-LIHC and ICGC-LIRI-JP datasets, and 499 and 185 DEGs between the high- and low-risk groups were obtained ( Tables S4 , S5 ), respectively. To explore the biological functions of the DEGs in the high- and low-risk groups, we again performed GO enrichment and KEGG pathway analyses. In the TCGA-LIHC dataset, GO enrichment analysis indicated significant enrichment in the organelle fission, nuclear division, chromosomal region and ATPase activity terms ( Figure 10A ). GO enrichment analysis of the ICGC-LIRI-JP dataset showed similar outcomes to the TCGA-LIHC dataset ( Figure 10B ). Additionally, KEGG pathway analysis of the TCGA-LIHC dataset showed significant enrichment in the cell cycle, oocyte meiosis and progesterone-medicated oocyte maturation pathways ( Figure 10C ). In the ICGC-LIRI-JP dataset, KEGG pathway analysis also demonstrated the analogical outcomes of the TCGA-LIHC dataset ( Figure 10D ).

To explore the correlation between the risk score and immune status, we calculated the infiltrating scores of 16 immune cells and 13 immune-related functions or pathways using ssGSEA. The scores of activated dendritic cells (aDCs), mast cells and follicular helper cells (Tfhs) were notably different between the high- and low-risk groups in the TCGA-LIHC dataset (all adj. P<0.001, Figure 11A ). In the TCGA-LIHC dataset, the scores of cytolytic activity, type I interferon (IFN) response and type II IFN response were obviously higher in the low-risk group, while the score of MHC class I was lower in the low-risk group (all adj. P<0.01, Figure 11B ). Moreover, the ICGC-LIRI-JP dataset showed that aDCs, mast cells, MHC class I and type II IFN responses were significantly different between the two risk cohorts ( Figures 11C, D ), which is consistent with the results of the TCGA-LIHC dataset.