Gene mutation and mRNA and miRNA expression data as well as clinical information were downloaded from TCGA¹ (Deng et al., 2016). In TCGA, liver hepatocellular carcinoma (LIHC) samples are divided into two groups: the first group contains HCV RNA or genotype or hepatitis C antibody in the patient’s clinical information, and the other group does not; here, we named the two groups HCC-HCV and HCC. The gene/miRNA microarray as verifying cohorts GSE154211 (Wang et al., 2021) and GSE119159 (Umezu et al., 2020) were downloaded from GEO database. The data were normalized, and R and its packages were employed in all analysis steps.

MuTect2 Somatic Mutation data, analyzed using MuTect2, were download from TCGA. TCGA provides somatic mutation data in the MAF format. Therefore, we visualized somatic mutations using the R package “maftools” (Mayakonda et al., 2018). In total, 96 HCC-HCV samples and 269 HCC samples were present in the dataset. We calculated the mutational status of genes using the algorithm in maftools, and the genes with p < 0.05, OR > 2, and number of mutations > 5 were selected as the significantly and differentially mutated genes.

According to the groupings, we performed normalization and differential gene expression analysis using the R package “edgeR.” False discovery rate (FDR) < 0.01 and | log₂ fold change (FC)| > 1 were used as cutoffs for identify differentially expressed genes (DEGs) for further analysis. Two R packages “pheatmap” and “ggplot2” were used for visualizing the heatmaps and volcano maps, respectively.

In total, 139 HCV patients were enrolled in the GSE119159. They included 99 patients who had not developed HCC and 40 who had developed HCC. A total of 10 samples (tumor and non-tumor regions) from two HCV-related HCC patients and three HCC patients were used to find the gene candidates in HCV-related HCC in the GSE154211. For differential expression analysis, we used the R package “limma.” | log2FC| > 1 and logFDR < 0.01 were used as cutoffs to identify DEGs for further analysis.

The human transcription factor (TF) and miRNA regulatory networks were constructed by integrating miRTarBase, TRANSFAC, and TransmiR (Vlachos et al., 2015; Chou et al., 2018; Tong et al., 2019). The three databases include curated interactions among human TFs, miRNAs, and target genes. We uniformly named the genes and miRNAs within the regulatory networks according to the National Center for Biotechnology Information (NCBI) and miRbase databases. Moreover, all regulatory relationships within the regulatory network were supported experimentally. In total, 888 TFs, 1,072 miRNAs, 3,150 target genes, and 18,056 edges were discovered in the regulatory network.

The key regulatory gene symbols were converted to Entrez ID using the R package “org.Hs.eg.db.” To identify the biological pathways involved in HCV-HCC occurrence and development, we employed Gene Ontology–biological process (GO-BP) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and visualized the results using the R packages “clusterProfiler” and “ggplot2.”

We constructed an HCV-HCC-related subnetwork and identified key regulators from the subnetworks. Next, we investigated whether the key regulators could distinguish HCC patients with good or poor outcomes. From these data, we obtained TCGA HCC dataset with mRNA/miRNA expression and clinical information. Then, we used the key regulator expression values and mutation information to cluster all patients into two groups. The differential survival of the two study groups was finally assessed using the log-rank test.

We downloaded and analyzed the somatic mutation data of 392 TCGA-LIHC samples. The mutation information of all genes in the samples is displayed as a waterfall diagram, with different colors representing different mutation types (Figure 1A). Further analysis showed that missense mutation, single-nucleotide polymorphisms (SNPs), and C > T accounted for the highest proportion of the variations (Figures 1B–D). The median number of variations in all HCC samples was 74.5, and the maximum number of variations in a single sample was 1,250 (Figure 1E). The number of variations in different classifications in all samples is shown in a box diagram (Figure 1F). The top 10 mutated genes in the 392 samples were TTN (25%), TP53 (28%), CTNNB1 (24%), MUC16 (16%), PCLO (11%), ALB (11%), RYR2 (10%), ABCA13 (9%), MUC4 (10%), and APOB (9%; Figure 1G). In total, 96 HCC-HCV samples and 269 HCC samples were present in TCGA dataset; the survival analysis indicated that HCC patients without HCV lived significantly longer than HCC-HCV patients (Figure 2B). With the use of the maftools algorithm, 41 differentially mutated genes were identified (Supplementary Table 1). The top 10 differential mutated genes between the two groups of patients were UNC5D (6–0), MYRF (5–0), PGLYRP4 (5–0), PREX2 (11–7), EPHA4 (9–5), HECTD4 (8–4), REV3L (8–4), HIPK2 (6–2), CHST3 (5–1), and TRO (5–1; Figure 2A). Moreover, HCC-HCV patients with mutations in some genes had a poor prognosis (Supplementary Figure 1).

Differentially expressed mRNAs and miRNAs were identified from two raw datasets: one containing GSE154211 and GSE119159 downloaded from the GEO database and another dataset from TCGA database. In total, 530 mRNA and 30 miRNA transcripts were observed to be expressed differentially in the HCC-HCV samples compared with HCC samples in TCGA dataset—including, respectively, 412 and 25 upregulated and 118 and five downregulated transcripts. Hierarchical clustering showed systematic variations in mRNA and miRNA expression in the HCC-HCV and HCC samples (Supplementary Figure 2). To identify the genes related to HCC-HCV in GSE154211, we first divided the expression data into four groups to identify DEGs between (A) HCC vs. HCC-HCV-adjacent, suggesting related to HCV-related carcinogenesis; (B) HCC-HCV vs. HCC, suggesting related HCV-related hepatocarcinogenesis; (C) HCC-HCV-adjacent vs. HCC-adjacent, suggesting related to HCV-related non-oncogenic effects; and (D) HCC vs. HCC-adjacent, suggesting related to non-HCV-related carcinogenesis. Four groups of data were then analyzed. Consequently, we identified 1,494 DEGs belonging to group A or B, but not group C or D, as genes with strong potential to be relevant to HCC-HCV (Supplementary Figure 3). In addition, 21 miRNA transcripts were observed to be differentially expressed in the developed HCC samples compared with the non-developed HCC samples in GSE119159, including nine upregulated and 12 downregulated transcripts.

To mine HCV-HCC-related regulatory relationships, we first constructed a TF–miRNA–mRNA regulatory network as a background network. Then, the HCV-induced HCC-related subnetwork was constructed by mapping DEGs into the background network. The nodes in the subnetwork contained DEGs and genes directly connected to the DEGs. In total, 359 TFs, 395 miRNAs, 739 target genes, and 2626 edges were present in the subnetwork.

We next mined the core regulatory module from the subnetwork by extracting the top 20 nodes ranked by closeness centrality and the edges among them. Notably, the regulatory relationships between these 20 nodes and differential mutated genes were added into the core regulatory module (Figure 3 and Supplementary Table 2). Finally, the core regulatory module contained 24 nodes and 36 edges.

To analyze the function of genes in the module, we conducted enrichment analysis of GO and KEGG, with FDR < 0.05 used as the cutoff to identify statistically significant GO terms and KEGG pathways. We found that many GO terms and KEGG pathways were implicated in the HCV-HCC processes in previous studies. As shown in Figure 4A and Supplementary Table 3, in the biological process and molecular function categories, the significantly enriched genes were for vasculature development regulation (Vescovo et al., 2016), ameboidal-type and epithelial cell migration (Khera et al., 2017), cell aging (Naggie, 2017), cell-matrix adhesion (Ninio et al., 2019), and melanocyte differentiation and angiogenesis involved in wound healing (Mohsen et al., 2014). Furthermore, KEGG pathway analysis showed that the significantly enriched genes were for small cell lung cancer, miRNAs in cancer, PI3K–Akt signaling pathway (Cheng et al., 2015), Ras and p53 signaling pathways (Vescovo et al., 2016), cellular senescence (Shiu et al., 2017), endocrine resistance, and advanced glycation end products (AGE)–receptor for AGE (RAGE) signaling pathway in diabetic complications (Hyogo and Yamagishi, 2008; Figure 4B and Supplementary Table 4).

We further analyzed the genes in the core regulatory module and found that expression of EXO1, VCAN, has-miR-200c-5p, BMI1, has-miR-204-5p, and KIT was significantly correlated with HCC prognosis in all patients; and thus, these genes were considered key regulators (Figure 5A). In particular, we found that the patients with low EXO1, VCAN, or KIT expression had an adverse outcome (HR < 1; Figure 5A). The HCC-HCV patients with mutations in these three genes have possibly also poor prognoses (Figures 5B–D). They may be potential biomarkers to predict the prognosis of patients at the genomic level. Moreover, we found that has-miR-200c-5P was significantly overexpressed in HCC-HCV samples (Figure 6A). The survival time of patients with high has-miR-200c-5P expression was significantly lower than that of patients with low expression (Figure 6B), suggesting that has-miR-200c-5P may be a potential biomarker to predict the prognosis of patients at the expression level.

Mutations in specific locations in EXO1 have been reported to inactivate proteins that increase cancer susceptibility (Welchew et al., 2002). KDR was also a significantly differential mutated gene in the module. KDR is the principal receptor that promotes the proangiogenic action of vascular endothelial growth factor and is involved in the tumorigenesis and progression of many malignancies, including HCC (Zheng et al., 2014). Moreover, BCL2 was the downstream gene in the core regulatory module. BCL2 can be functionally divided into antiapoptotic and proapoptotic groups. The balance between these two groups may determine the fate of a tumor cell. In HCC, this balance is often tilted toward the antiapoptotic members, leading to resistance to death and rapid proliferation in cancer cells (Alenzi et al., 2010). BCL2 expression in the HCC-HCV samples was lower than that in the HCC samples, but the difference was non-significant in TCGA data—which may be a reason for the worse prognosis of HCC-HCV.