OL Oncology Letters 1792-1074 1792-1082 D.A. Spandidos 10.3892/ol.2018.8221 OL-0-0-8221 Articles Identification and functional analysis of a core gene module associated with hepatitis C virus-induced human hepatocellular carcinoma progression BaiGaobo ZhengWenling MaWenli Institute of Genetic Engineering, Southern Medical University, Guangzhou, Guangdong 510515, P.R. China Correspondence to: Professor Wenli Ma, Institute of Genetic Engineering, Southern Medical University, 1838 North Guangzhou Avenue, Guangzhou, Guangdong 510515, P.R. China, E-mail: wenlima668@sohu.com 05 2018 09 03 2018 15 5 6815 6824 31012016 27022017 Copyright: © Bai et al. 2018 This is an open access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

Hepatitis C virus (HCV)-induced human hepatocellular carcinoma (HCC) progression may be due to a complex multi-step processes. The developmental mechanism of these processes is worth investigating for the prevention, diagnosis and therapy of HCC. The aim of the present study was to investigate the molecular mechanism underlying the progression of HCV-induced hepatocarcinogenesis. First, the dynamic gene module, consisting of key genes associated with progression between the normal stage and HCC, was identified using the Weighted Gene Co-expression Network Analysis tool from R language. By defining those genes in the module as seeds, the change of co-expression in differentially expressed gene sets in two consecutive stages of pathological progression was examined. Finally, interaction pairs of HCV viral proteins and their directly targeted proteins in the identified module were extracted from the literature and a comprehensive interaction dataset from yeast two-hybrid experiments. By combining the interactions between HCV and their targets, and protein-protein interactions in the Search Tool for the Retrieval of Interacting Genes database (STRING), the HCV-key genes interaction network was constructed and visualized using Cytoscape software 3.2. As a result, a module containing 44 key genes was identified to be associated with HCC progression, due to the dynamic features and functions of those genes in the module. Several important differentially co-expressed gene pairs were identified between non-HCC and HCC stages. In the key genes, cyclin dependent kinase 1 (CDK1), NDC80, cyclin A2 (CCNA2) and rac GTPase activating protein 1 (RACGAP1) were shown to be targeted by the HCV nonstructural proteins NS5A, NS3 and NS5B, respectively. The four genes perform an intermediary role between the HCV viral proteins and the dysfunctional module in the HCV key genes interaction network. These findings provided valuable information for understanding the mechanism of HCV-induced HCC progression and for seeking drug targets for the therapy and prevention of HCC.

 hepatitis C virus hepatocellular carcinoma co-expression differential co-expression interaction network Introduction

Hepatocellular carcinoma (HCC) is one of the most common types of cancer in the world (1), and it was the third leading cause of cancer-associated mortality worldwide in 2013 (2). A number of effective measures, including percutaneous ablation, surgical resection and liver transplantation, have been developed for HCC therapy (3). Although the short-term prognosis of patients with HCC has improved due to advances in early diagnosis and treatments, the long-term prognosis remains poor due to frequent recurrence or metastasis (4).

HCC is a heterogeneous cancer that usually develops in patients with chronic liver disease, particularly those with cirrhosis. The progression between chronic liver disease and HCC is a complex and multistep process. Chronic viral hepatitis caused by hepatitis C virus (HCV) infection has become one of the main causes of HCC. A typical HCV-induced HCC progression may undergo the following five successive stages: Normal; cirrhosis; dysplasia; early HCC; and advanced HCC (5).

In the development from hepatic lesions to HCC, not only does the expression level of numerous genes change, but dynamic changes may also appear in the regulation of different genes and protein-protein interaction (6). Classical approaches, including identification of differentially expressed genes, may not reveal the complex interactions and functional association among genes in biological processes. Alternatively, gene co-expressed network analysis may provide a powerful approach for elucidating the co-regulation and interaction between proteins in biological processes (7). Comparing with identification of differentially expressed genes network-based analysis may provide valuable information for understanding complex interaction among genes in biological processes (7,8). At present, a number of studies have identified the altered gene co-expression associated with tumors (9–11).

HCV is a type of single-stranded RNA virus that replicates in the cytoplasm of host hepatocytes. HCV viral proteins are mainly divided into two classes: Structural proteins (CORE, envelope proteins 1 and 2 and p7) and non-structural (NS) proteins (NS2, NS3, NS4B, NS5A and NS5B) (12). Multiple important HCV viral proteins, including CORE, NS3, NS4B and NS5A, were identified to target important cancer-associated proteins. These proteins can promote cell growth and cell cycle deregulation and destroy the stable structure of the host cell genome. These viral proteins were revealed to potentiate oncogenic transformation (12) and perform crucial roles in HCC (13).

However, the exact role of HCV in HCC progression remains to be determined. Therefore, studying the association between the deregulated biological networks and HCV may help to understand the molecular mechanisms of HCC development.

 Materials and methods <sec> <title>Microarray data and differentially expressed genes

The HCC gene expression dataset GSE6764 (CEL file) was downloaded (14) from the Gene Expression Omnibus database (15). GSE6764 was based on the GPL570 platform: Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA). The dataset contained 75 tissue samples representing the stepwise carcinogenic process from pre-neoplastic lesions (cirrhosis and dysplasia) to HCC. The raw sample groups were categorized into five groups: Normal; cirrhosis; dysplasia, including low-grade dysplastic nodules and high-grade dysplastic nodules; early HCC group, including early HCC and extremely early HCC; and advanced HCC group, including advanced HCC and extremely advanced HCC. The numbers of samples contained in these groups were 10, 13, 17, 18 and 17, respectively.

The Robust Multichip Analysis (RMA) algorithm in the ‘affy’ package (16) (http://www.bioconductor.org/packages/3.6/bioc/html/affy.html) in Bioconductor (17) (http://www.bioconductor.org/) was used to perform background correction, quantile normalization and summarization. For multiple probe sets matching to the same gene, the probes with the largest variation were selected. The linear models for microarray and RNA-sequencing data (LIMMA) software package (https://bioconductor.org/packages/release/bioc/html/limma.html, The Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia, Version 3.24.15) (18) in Bioconductor was utilized for linear models and testing for differential expression, and to adjust for multiple testing. t-tests and F-tests were performed on the matrix. The genes with 2-fold expression change (increase or decrease) (P<0.001) and a false discovery rate (FDR) <0.05, compared with the normal group, were selected as differential expression genes (DEGs) in each disease stage. To reveal the expression dynamic of DEGs across all stages, the DEGs in each stage were merged into a DEG set for subsequent analysis.

 Analysis of the modules significantly associated with disease progression

In order to identify the modules of highly correlated genes involved in disease progression, Weighted Gene Co-expression Network Analysis (WGCNA) software (https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/) (19) was implemented in R language (https://cran.r-project.org/, version 3.2.2). First, the WGCNA package computed co-expressed association among genes in the DEG set and inferred the co-expressed networks in 75 samples. Based on a scale-free topology criterion (20), β=12 was selected as the soft threshold power in the present study. Deep split was set as 2, cut height was set as 0.975 and ‘minModuleSize’ was set as 10, and other parameters were set at default levels.

Module eigengene (ME) is the first principle component and explains the majority of information of the module genes. To identify modules that were significantly associated with disease progression, Spearman's rank correlation coefficient was used to analyze the association between the ME of each module and disease stages. Modules that were significantly correlated with disease progression were labeled as candidate modules (absolute correlation coefficient |r|>0.8; P<0.01).

 Screening for the key genes significantly associated with disease progression

The present study aimed to identify key genes with functional similarity and high co-expression over all disease developmental stages that were significantly associated with disease progression. The expression pattern of these genes may reflect the trajectory of disease progression. In candidate modules, the genes with absolute gene significance value (GS) >0.75 and modular membership (MM) >0.8 were selected as genes that were highly associated with disease stages. To exclude the gene pairs with lower correlation in all disease stages, the correlation coefficient between all reserved gene pairs were examined. A threshold of |r|=0.85 (P<0.01) was set as the screening criterion. The gene pairs with an absolute correlation coefficient higher than this threshold were retained. The Database for Annotation, Visualization and Integrated Discovery 6.7 (DAVID, https://david-d.ncifcrf.gov/) (21) was utilized to assess the function of these genes. By Gene Ontology (GO) enrichment analysis, the genes with functions associated with disease progression were retained. The heatmap and self-organizing tree algorithm (SOTA) (22) in cogena package (http://www.bioconductor.org/packages/3.2/bioc/html/cogena.html) was implemented for clustering and visualization of expression patterns of the key genes.

 Analysis of the interactions among the key genes

GeneMANIA (http://genemania.org/) (23,24) is a web-based tool for the construction of the composite gene-gene functional interaction network and the prediction of protein function, based on multiple networks derived from different databases. The interactions among key genes, including physical and genetic interactions, pathways, co-expression, co-localization, shared protein domain and predicted interactions were analyzed using GeneMANIA.

 Differential co-expressed analysis of the key genes

The genes retained as seeds were used to search for highly co-expressed genes (absolute correlation coefficient |r|>0.85, P<0.01) in the DEG set in which seed genes were ruled out in each stage. The software package DiffCorr (https://sourceforge.net/projects/diffcorr/, Version 0.4.1, Slashdot media and dice, Inc., San Jose, CA, USA) (25) in R language was used to search for the differential correlation gene pairs between adjacent stages. Fisher's z-test was used to identify significant differences between adjacent stages. First, the correlation coefficients for each of the two adjacent stages were transformed into z by Fisher's transformation, separately as follows:

z=12log1+r1–r

Then, differences between the two correlation coefficients were tested using the equation:

Z=za–zb1na–3+1nb–3

na and nb represent the sample size for each of the disease stages for each gene pair. The gene pairs with FDR<0.05 and P<0.001 of differential correlations were considered as differential coexpression between adjacent phases.

Cytoscape software (http://www.cytoscape.org/, version 3.2; Cytoscape Consortium, La Jolla, CA, USA) (26) was used for the construction of a differentially co-expressed network of the key genes.

 Construction of interaction networks between HCV viral proteins and products of key genes

Search Tool for the Retrieval of Interacting Genes (STRING) (http://version10.string-db.org/, version 10.0) is an online protein-protein interaction association database curated from experimental interactions of proteins and comprehensive information (27). First, the key genes were directly mapped to the STRING database to get pairwise interactions (combined scores >0.4). Next, the HCV viral proteins and their directly targeted proteins in the key genes set were extracted from the literature and a comprehensive interaction dataset containing 481 interaction pairs between HCV viral proteins and human proteins from yeast two-hybrid experiments (28). Furthermore, the interaction network of the key genes and the interaction network between HCV and their targeted proteins were combined to make an interaction network between HCV viral proteins and key genes, by anchoring the overlapped proteins of these two networks. The interaction network between HCV viral proteins and key genes was visualized using Cytoscape software (26).

 Results <sec> <title>DEGs between disease groups and normal group

Subsequent to performing the RMA algorithm for data preprocessing, the LIMMA software package (18) was run to obtain 509, 156, 802 and 1,501 DEGs from the cirrhosis, dysplasia, early HCC and advanced HCC groups, respectively. The DEGs from each disease group were merged into a gene set, termed the DEG set, which contained a total of 1,973 genes.

Identified modules associated with the disease progression. It was hypothesized that there are core gene modules that can exhibit the disease progression dynamic. It is thought that genes with similar expression patterns and similar functions may be regulated by the same mechanisms (29). Therefore, the modules consisting of genes with similar expression patterns and similar functions were expected to be identified via WGCNA (19). The expression profile of the DEG set was used as the input of WGCNA.

Based on the dissimilarity corresponding to the Topological Overlap Matrix, hierarchical average linkage clustering produced a hierarchical clustering tree (dendrogram) of genes. A total of 10 modules corresponding to the 10 branches of the dendrogram were obtained (Fig. 1A), based on provided parameters. The modules were labeled by 10 colors. The sizes of the 10 modules were ranked between 11 and 1,018 genes. The grey module contained the genes with dissimilar expression patterns.

In order to identify the association between the modules and disease progression, MEs were utilized, instead of single genes, as the overall expression levels of the module. The association between MEs and disease progression obtained was presented as a heatmap (Fig. 1B). In the 10 modules, the most strongly associated with the disease stage modules (|r|>0.8; P<0.01), the turquoise and blue modules, were selected as the candidate modules. The turquoise module was positively associated with disease stages, while the blue module was negatively associated with disease stages.

 Key genes associated with disease development and progression

Since the genes in modules may not have equivalent contribution to the association of the module with disease stages, GS and MM were used to identify the genes highly associated with disease stages in the two candidate modules. The numbers of retained genes (|GS|>0.75 and |MM|>0.8) were 53 in the turquoise module and 48 in the blue module (Fig. 1C and D). To yield the consistent co-expression of genes in modules across all stages, the pairwise correlation coefficient among genes in the two candidate modules was examined. The genes that weakly correlated (|r|<0.85 or P>0.01) with any gene in the same module were removed. Finally, 44 genes in the turquoise module and 31 genes in the blue module were obtained. The significantly enriched GO terms are shown in Tables I and II.

For those retained genes in the blue module, the GO terms that were significantly enriched were associated with carbohydrate binding (GO:0030246), sugar binding (GO:0005529), protein processing (GO:0016485), protein maturation (GO:0051604), complement activation (GO:0006956) and activation of plasma proteins involved in the acute inflammatory response (GO:0002541). These genes may be involved in energy metabolism and protein synthesis.

The GO terms that the genes retained in the turquoise module significantly enriched included mitotic cell cycle (GO:0000278), M phase (GO:0000279), cell cycle phase (GO:0022403), cell cycle process (GO:0022402), mitosis (GO:0007067), nuclear division (GO:0000280) and cell division (GO:0051301). In previous studies the process of tumorigenesis, cell cycle regulation was significantly altered and mitosis was significantly faster (30). The protein associated with cell cycle and mitosis notably changed along with tumorigenesis (31). In the turquoise module, the functions of these genes involved in the process of mitotic cell cycle and cell division were associated with tumorigenesis.

The heatmap showed that the expression of certain genes changed gradually between the normal stage and HCC; specifically, between the dysplasia stage and advanced HCC (Fig. 2A). Furthermore, the expression intensity of those genes in early HCC was in a transitional state between the dysplasia stage and advanced HCC. The expression patterns and clustering of these genes in the turquoise module are shown in Fig. 2B. A total of 44 genes were grouped into five clusters, with an evident trend along with the disease progression between the normal stage and advanced HCC. As the majority of the 44 genes showed gradually increased expression, they were selected as key genes of the core module. The interactions between the 44 key genes were analyzed using GeneMANIA (23,24). The co-localization, co-expression, pathway, shared protein domains, genetic, physical and predictive interactions of the 44 key genes are shown in Fig. 2C.

 Altered association of co-expression with key genes in disease progression

As aforementioned, the changes in expression of the 44 key genes reflected the progression dynamic between the normal stage and advanced HCC. Therefore, it was inferred that co-expressed associations among these genes and others may be altered correspondingly in this process. The differential co-expression genes of 44 key genes in disease progression were investigated. The results showed that the association of co-expression of 9 gene pairs was enhanced (P<0.01, FDR<0.05), while for 76 gene pairs it was reduced (P<0.01, FDR<0.05), between the normal and cirrhosis stages (Fig. 3A). In the cirrhosis-dysplasia stage, there were only 4 pairs of differentially co-expressed genes. In total, 103 pairs of genes were differentially co-expressed in dysplasia-early HCC stages. Among those, the correlations of the 32 gene pairs demonstrated either significant decrease, including the changes from strong correlation to weak correlation or to no correlation, or a transformation from positive correlation to negative correlation (the difference of the correlation coefficients in the adjacent phases: r2-r1<0; P<0.01; FDR<0.05). By contrast, the 71 gene pairs exhibited either significant increase in correlation, including from weak correlation or no correlation to strong correlation, or the transformation from negative correlation to positive correlation (r2-r1>0; P<0.01; FDR<0.05; Fig. 3B). However, between the early HCC stage and the advanced HCC stage, no differentially co-expressed gene pairs were identified. This analysis showed that the association of co-expression among key genes and other genes had maximum ratio of change in normal-cirrhosis stages and dysplasia-early HCC stages.

In normal-cirrhosis stages, genes that were differentially co-expressed with key genes significantly enriched the GO terms: Organic acid biosynthetic process (GO:0016053) and carboxylic acid biosynthetic process (GO:0046394) (P<0.01). According to the functional annotation of the genes, these genes are involved in the biological synthesis of the liver. The differential co-expression indicated the disorder in some of biosynthesis functions of the liver during the transition from normal to cirrhosis. In dysplasia-early HCC stages, in cases of weakened correlation, no GO terms were significantly enriched in those genes (P<0.01). The genes that enhanced correlation with key genes, including aurora kinase A (AURKA), marker of proliferation Ki-67 (MKI67), baculoviral IAP repeat containing 5 (BIRC5), ZW10 interacting kinetochore protein (ZWINT), cell division cycle associated 3 (CDCA3), cyclin (CCN) E2, kinesin family member 18A (KIF18A), NUF2 and nucleolar and spindle-associated protein 1 (NUSAP1), were involved in cell cycle, nuclear division, mitosis, spindle, microtubule-based process and spindle (Table III), which were closely associated with cancer. The differential co-expression reflected the dysfunction of cell cycle, mitosis and other associated processes during the development from dysplasia to HCC.

 Interaction networks between HCV viral proteins and their targeted proteins in the key genes set

In the HCV-human interactome from Chassey's yeast two-hybrid experiments (28), HCV viral proteins NS3 and NS5A were revealed to target the proteins of NDC80 and cyclin-dependent kinase 1 (CDK1) separately, which were present in the key gene set. In addition, HCV NS5B was reported to specifically interact with CCNA2 (32) and Rac GTPase-activating protein 1 (RACGAP1) (33). The interaction networks between HCV viral proteins (NS3, NS5A and NS5B) and their targeted proteins in the key genes set are shown in Fig. 4.

 Discussion

HCV infection has become one of the predominant causes of HCC in patients with chronic viral hepatitis (34). The exploration of HCV-induced HCC progression may help to understand the common developmental mechanism of HCC induced by other risk factors, including HBV and alcoholic liver diseases. In the present study, the WGCNA tool was performed for screening of the module consisting of the key genes that exhibited HCC progression. The 44 genes identified as key genes were associated with disease progression. The majority of these genes had similar expression patterns and were associated and interconnected with each other. The enriched GO terms were mainly involved in the cell cycle process, mitotic cell cycle, M phase, cell cycle phase, mitosis, nuclear division, cell division and phosphorylation of proteins. This indicated that the dynamic of the module consisting of these genes may contribute to or be driven by disease progression. As cell cycle deregulation is one of the common hallmark traits of cancer (35), dynamic features of this module may also be associated with numerous other cancer progressions.

In the key genes set, genes including CDK1, cell-division cycle protein 20 (CDC20), CCNB2, NIMA related kinase 2 (NEK2) and CCNB1 were known to be involved in common regulatory processes of the cell cycle, mitosis and cell division. Additionally, a number of the 44 genes were closely associated with HCC. DLG associated protein 5 (DLGAP5), also termed HURP (hepatoma upregulated protein), is overexpressed in HCC (36,37). DLGAP5 activates p38/nuclear factor κ-light chain-enhancer of activated B cells (NF-κB), and combines NF-κB into the HURP/NF-κB complex to regulate CCNE1 expression (38). Previous studies have revealed that silencing or knockdown of DLGAP5 significantly inhibited the proliferation and invasion of HCC cells (36,37). The TTK gene encodes a dual specificity protein kinase with the ability to phosphorylate tyrosine, serine and threonine. TKK is essential for chromosome alignment at the centromere during mitosis and is required for centrosome duplication. TTK is upregulated in the majority of HCC specimens and its overexpression can promote cell proliferation, anchor-dependent colony formation and resistance to sorafenib of HCC cells (39). CDKN3 is frequently upregulated in HCC and is associated with poor prognosis of HCC. Overexpression of CDKN3 can stimulate the proliferation of HCC cells by promoting G1/S phase transition (40). In addition, certain genes are associated with multiple cancers. For example, enhancer of zeste homolog 2 (EZH2) promotes lung cancer progression via the vascular endothelial growth factor-A/AKT signaling pathway in non-small cell lung cancer (41), and it is also involved in aggressive breast cancer (42,43).

The analysis indicated that the changes in expression levels of key genes were associated with transition between chronic liver disease and HCC. These changes may contribute to or be driven by the dysregulation of the interaction or regulated associations between certain genes and these key genes in the disease progression.

Differentially co-expressed analysis identified several genes that showed enhanced association of co-expression with key genes in dysplastic-early HCC stages, including AURKA, MKI67, BIRC5, cell division cycle associated 3 (CDCA3), ZWINT, NUSAP1, ubiquitin conjugating enzyme E2 C (UBE2C), CCNE2, KIF18A and NUF2, which were highly associated with cancer. In other enhanced co-expression genes, the α-fetoprotein (AFP) gene, is a well-known marker for primary HCC.KIAA0101 (alternative symbol NS5ATP9) is highly co-expressed with several key genes. A previous study revealed that KIAA0101 may be involved in the pathogenesis of HCV-associated HCC and is upregulated by HCV NS5A protein (44). Furthermore, KIAA0101 potentially has an important role in NS5A-induced hepatocyte autophagy (45). Notably, the correlation coefficient between ZWINT and NEK2 in the early HCC stage (r=0.904; P=2.69×10−7) is significantly increased (P=1.19×10−5, FDR=0.022) compared with the dysplastic stage (r=−0.134). A previous study reported that the functions of the two genes were associated with chromosomal instability, which is a major characteristic of numerous cancers (46). Similarly, pituitary tumor transforming gene 1 (PTTG1), differentially co-expressed with ZWINT, is also associated with chromosomal instability (46). Therefore, these findings indicated that chromosomal instability of HCC may be associated with not only the increasing expression intensity of these genes, but also the altered interaction among them. Centrosomal protein of 55 kDa (CEP55) and a DNA helicase/putative stem cell marker HELLS (gene encoding lymphoid-specific helicase) are downstream targets of the key oncogene forkhead box M1 (FOXM1). Waseem et al suggested the proliferation-associated gene set consisting of FOXM1, CEP55 and HELLS shared a progressive expression pattern during the progression of head and neck squamous cell carcinoma (47). In this analysis, CEP55 and HELLS showed enhanced correlation in the early HCC stage, compared with the dysplastic stage. Whether these genes possess a similar function in HCC, as demonstrated in a previous study (47), is worth additional investigation. These altered interactive associations among important genes may be associated with tumorigenesis.

To study the interaction association between the key genes identified and HCV, an interaction network between HCV viral proteins and their targeted proteins, which were contained in the key genes set, was constructed by mining the HCV-human database and published literature. In the present study, CDK1, NDC80, CCNA2 and RACGAP1 were revealed to be targets of HCV NS5A, NS3 and NS5B, respectively. CDK1 is a key regulatory kinase of the cell cycle in the CDK family. A previous study demonstrated that CDK1 is upregulated in HCC (48), and serves a crucial role of the G2/M modulators in the cell cycle and cell proliferation of HCC (49). As CDK1 activation is required for mitosis, the increase of CDK1 activation in HCC may indicate the increase of cancer cell proliferation. NDC80 is a coiled-coil protein critical for cell mitosis, and performs essential roles in chromosome segregation by interacting with several proteins that modulate the G2/M phase through its coiled-coil domains (50). NDC80 is overexpressed in a variety of human cancers, such as gastric and breast cancer (51,52), and is associated with multiple cancers (53,54).

In the key genes set, CCNA2 is a cell cycle regulatory protein that can bind to and activate CDK1 or CDK2 kinases, and thus promote cell cycle G1/S and G2/M transitions. NS5B was reported to specifically interact with CCNA2 in vitro and in vivo (32). However, unlike NDC80 and CDK1, CCNA2 has dual regulatory roles in cell cycle and viral propagation. CCNA2 is required for HCV replication. Small interfering RNA-mediated depletion of CCNA2 may significantly inhibit HCV replication in HCV subgenomic replicon cells and HCVcc-infected cells (32). In addition, RACGAP1 belongs to the GTPase-activating protein family (55). A previous study suggested that the level of RACGAP1 is upregulated by HCV infection in human hepatoma cells to enhance HCV replication by binding to and affecting NS5B polymerase activity (33). The interaction network between the key genes and HCV showed that a number of important proteins directly interacted with CDK1, NDC80, CCNA2 and RACGAP1. Therefore, these proteins that are targeted by HCV may perform an intermediary role between the HCV viral proteins and the dysfunctional module. The results suggested that CDK1, NDC80, RACGAP1 and CCNA2 inhibition or interference by drugs may have potential as an effective method of therapy for HCC and prevention of the HCC progression.

Overall, in the present study, using the WGCNA tool, a core gene module was identified, in which the dynamic characteristics of the genes were associated with HCV-induced HCC progression. The marked changes in expression levels and association of co-expression of the key genes may contribute to or be driven by progression of HCV-induced HCC. The changes may be markers of transition from precancerous to HCC states. Particularly, in the identified core module, CDK1, NDC80, CCNA2 and RACGAP1 were revealed to be targeted by HCV viral proteins. These findings confirmed the close link between the identified module and HCV-induced HCC. The present study may be helpful for understanding the molecular mechanism that underlies the progression of HCV-induced HCC and provides valuable information for seeking the drug targets for the therapy of HCC and prevention of HCC development.

 Acknowledgements

Not applicable.

 Funding

No funding was received.

 Availability of data and materials

The datasets analyzed during the current study are available in the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6764).

 Author's contributions

GB, WZ and WM conceived and designed the experiments, GB performed the experiments and analyzed the data. GB, WZ and WM drafted the manuscript.

 Ethics approval and consent to participate

Not applicable.

 Consent for publication

Not applicable.

 Competing interests

The authors declare that they have no competing interests.

 References TorreLABrayFSiegelRLFerlayJLortet-TieulentJJemalAGlobal cancer statistics, 2012CA Cancer J Clin6587108201510.3322/caac.2126225651787 Global Burden of Disease Cancer CollaborationFitzmauriceCDickerDPainAHamavidHMoradi-LakehMMacIntyreMFAllenCHansenGWoodbrookRThe global burden of cancer 2013JAMA Oncol1505527201510.1001/jamaoncol.2015.073526181261 YamaneBWeberSLiver-directed treatment modalities for primary and secondary hepatic tumorsSurg Clin North Am8997113200910.1016/j.suc.2008.10.00419186233 IshikawaTStrategy for improving survival and reducing recurrence of HCV-related hepatocellular carcinomaWorld J Gastroenterol1961276130201310.3748/wjg.v19.i47.886124115808 YuHLinCCLiYYZhaoZDynamic protein interaction modules in human hepatocellular carcinoma progressionBMC Syst Biol57 SupplS2201310.1186/1752-0509-7-S5-S2 ChoiJKYuUYooOJKimSDifferential coexpression analysis using microarray data and its application to human cancerBioinformatics2143484355200510.1093/bioinformatics/bti72216234317 LeeHKHsuAKSajdakJQinJPavlidisPCoexpression analysis of human genes across many microarray data setsGenome Res1410851094200410.1101/gr.191090415173114 HeBZhangHShiTA comprehensive analysis of the dynamic biological networks in HCV induced hepatocarcinogenesisPLoS One6e18516201110.1371/journal.pone.001851621526182 VarelasXBouchieMPKukuruzinskaMAProtein N-glycosylation in oral cancer: Dysregulated cellular networks among DPAGT1, E-cadherin adhesion and canonical Wnt signalingGlycobiology24579591201410.1093/glycob/cwu03124742667 MentzenWIFlorisMde la FuenteADissecting the dynamics of dysregulation of cellular processes in mouse mammary gland tumorBMC Genomics10601200910.1186/1471-2164-10-60120003387 SouthworthLKOwenABKimSKAging mice show a decreasing correlation of gene expression within genetic modulesPLoS Genet5e1000776200910.1371/journal.pgen.100077620019809 BanerjeeARayRBRayROncogenic potential of hepatitis C virus proteinsViruses221082133201010.3390/v209210821994721 McGivernDRLemonSMVirus-specific mechanisms of carcinogenesis in hepatitis C virus associated liver cancerOncogene3019691983201110.1038/onc.2010.59421258404 WurmbachEChenYBKhitrovGZhangWRoayaieSSchwartzMFielIThungSMazzaferroVBruixJGenome-wide molecular profiles of HCV-induced dysplasia and hepatocellular carcinomaHepatology45938947200710.1002/hep.2162217393520 BarrettTTroupDBWilhiteSELedouxPRudnevDEvangelistaCKimIFSobolevaATomashevskyMMarshallKANCBI GEO: Archive for high-throughput functional genomic dataNucleic Acids Res37(Database Issue)D885D890200910.1093/nar/gkn76418940857 GautierLCopeLBolstadBMIrizarryRAaffy-analysis of Affymetrix GeneChip data at the probe levelBioinformatics20307315200410.1093/bioinformatics/btg40514960456 GentlemanRCCareyVJBatesDMBolstadBDettlingMDudoitSEllisBGautierLGeYGentryJBioconductor: Open software development for computational biology and bioinformaticsGenome Biol5R80200410.1186/gb-2004-5-10-r8015461798 SmythGKLinear models and empirical bayes methods for assessing differential expression in microarray experimentsStat Appl Genet Mol Biol3Article3200410.2202/1544-6115.102716646809 LangfelderPHorvathSWGCNA: An R package for weighted correlation network analysisBMC Bioinformatics9559200810.1186/1471-2105-9-55919114008 ZhangBHorvathSA general framework for weighted gene co-expression network analysisStat Appl Genet Mol Biol4Article17200510.2202/1544-6115.112816646834 Huang daWShermanBTLempickiRASystematic and integrative analysis of large gene lists using DAVID bioinformatics resourcesNat Protoc44457200910.1038/nprot.2008.21119131956 HerreroJValenciaADopazoJA hierarchical unsupervised growing neural network for clustering gene expression patternsBioinformatics17126136200110.1093/bioinformatics/17.2.12611238068 Warde-FarleyDDonaldsonSLComesOZuberiKBadrawiRChaoPFranzMGrouiosCKaziFLopesCTThe GeneMANIA prediction server: Biological network integration for gene prioritization and predicting gene functionNucleic Acids Res38W214W220201010.1093/nar/gkq53720576703 ZuberiKFranzMRodriguezHMontojoJLopesCTBaderGDMorrisQGeneMANIA prediction server 2013 updateNucleic Acids Res41W115W122201310.1093/nar/gkt53323794635 FukushimaADiffCorr: An R package to analyze and visualize differential correlations in biological networksGene518209214201310.1016/j.gene.2012.11.02823246976 ShannonPMarkielAOzierOBaligaNSWangJTRamageDAminNSchwikowskiBIdekerTCytoscape: A software environment for integrated models of biomolecular interaction networksGenome Res1324982504200310.1101/gr.123930314597658 SzklarczykDFranceschiniAKuhnMSimonovicMRothAMinguezPDoerksTStarkMMullerJBorkPThe STRING database in 2011: Functional interaction networks of proteins, globally integrated and scoredNucleic Acids Res39(Database Issue)D561D568201110.1093/nar/gkq97321045058 de ChasseyBNavratilVTafforeauLHietMSAublin-GexAAgauguéSMeiffrenGPradezynskiFFariaBFChantierTHepatitis C virus infection protein networkMol Syst Biol4230200810.1038/msb.2008.6618985028 AlloccoDJKohaneISButteAJQuantifying the relationship between co-expression, co-regulation and gene functionBMC Bioinformatics518200410.1186/1471-2105-5-1815053845 MasVRMalufDGArcherKJYanekKKongXKulikLFreiseCEOlthoffKMGhobrialRMMcIverPFisherRGenes involved in viral carcinogenesis and tumor initiation in hepatitis C virus-induced hepatocellular carcinomaMol Med158594200910.2119/molmed.2008.0011019098997 De GiorgiVBuonaguroLWorschechATorneselloMLIzzoFMarincolaFMWangEBuonaguroFMMolecular signatures associated with HCV-induced hepatocellular carcinoma and liver metastasisPLoS One8e56153201310.1371/journal.pone.005615323441164 PhamLVNgoHTLimYSHwangSBHepatitis C virus non-structural 5B protein interacts with cyclin A2 and regulates viral propagationJ Hepatol57960966201210.1016/j.jhep.2012.07.00622796893 WuMJKePYHorngJTRacGTPase-activating protein 1 interacts with hepatitis C virus polymerase NS5B to regulate viral replicationBiochem Biophys Res Commun4541924201410.1016/j.bbrc.2014.10.00825305482 KoikeKMoriyaKKimuraSRole of hepatitis C virus in the development of hepatocellular carcinoma: Transgenic approach to viral hepatocarcinogenesisJ Gastroenterol Hepatol17394400200210.1046/j.1440-1746.2002.02763.x11982718 ChanKSKohCGLiHYMitosis-targeted anti-cancer therapies: Where they standCell Death Dis3e411201210.1038/cddis.2012.14823076219 KuoTCChangPYHuangSFChouCKChaoCCKnockdown of HURP inhibits the proliferation of hepacellular carcinoma cells via downregulation of gankyrin and accumulation of p53Biochem Pharmacol83758768201210.1016/j.bcp.2011.12.03422230478 LiaoWLiuWYuanQLiuXOuYHeSYuanSQinLChenQNongKSilencing of DLGAP5 by siRNA significantly inhibits the proliferation and invasion of hepatocellular carcinoma cellsPLoS One8e80789201310.1371/journal.pone.008078924324629 ChenJMChiuSCWeiTYLinSYChongCMWuCCHuangJYYangSTKuCFHsiaJYYuCTThe involvement of nuclear factor-kappaB in the nuclear targeting and cyclin E1 upregulating activities of hepatoma upregulated proteinCell Signal272636201510.1016/j.cellsig.2014.09.02025289861 LiangXDDaiYCLiZYGanMFZhangSRYin-PanLuHSCaoXQZhengBJBaoLFExpression and function analysis of mitotic checkpoint genes identifies TTK as a potential therapeutic target for human hepatocellular carcinomaPLoS One9e97739201410.1371/journal.pone.009773924905462 XingCXieHZhouLZhouWZhangWDingSWeiBYuXSuRZhengSCyclin-dependent kinase inhibitor 3 is overexpressed in hepatocellular carcinoma and promotes tumor cell proliferationBiochem Biophys Res Commun4202935201210.1016/j.bbrc.2012.02.10722390936 GengJLiXZhouZWuCLDaiMBaiXEZH2 promotes tumor progression via regulating VEGF-A/AKT signaling in non-small cell lung cancerCancer Lett359275287201510.1016/j.canlet.2015.01.03125633838 CollettKEideGEArnesJStefanssonIMEideJBraatenAAasTOtteAPAkslenLAExpression of enhancer of zeste homologue 2 is significantly associated with increased tumor cell proliferation and is a marker of aggressive breast cancerClin Cancer Res1211681174200610.1158/1078-0432.CCR-05-153316489070 KleerCGCaoQVaramballySShenROtaITomlinsSAGhoshDSewaltRGOtteAPHayesDFEZH2 is a marker of aggressive breast cancer and promotes neoplastic transformation of breast epithelial cellsProc Natl Acad Sci USA1001160611611200310.1073/pnas.193374410014500907 ShiLZhangSLLiKHongYWangQLiYGuoJFanWHZhangLChengJNS5ATP9, a gene up-regulated by HCV NS5A proteinCancer Lett259192197200810.1016/j.canlet.2007.10.02918068894 QuanMLiuSLiGWangQZhangJZhangMLiMGaoPFengSChengJA functional role for NS5ATP9 in the induction of HCV NS5A-mediated autophagyJ Viral Hepat21405415201410.1111/jvh.1215524750205 BrendleABrandtAJohanssonREnquistKHallmansGHemminkiKLennerPFörstiASingle nucleotide polymorphisms in chromosomal instability genes and risk and clinical outcome of breast cancer: A Swedish prospective case-control studyEur J Cancer45435442200910.1016/j.ejca.2008.10.00119008095 WaseemAAliMOdellEWFortuneFTehMTDownstream targets of FOXM1: CEP55 and HELLS are cancer progression markers of head and neck squamous cell carcinomaOral Oncol46536542201010.1016/j.oraloncology.2010.03.02220400365 LiKKNgIOFanSTAlbrechtJHYamashitaKPoonRYActivation of cyclin-dependent kinases CDC2 and CDK2 in hepatocellular carcinomaLiver22259268200210.1046/j.0106-9543.2002.01629.x12100577 ItoYTakedaTSakonMMondenMTsujimotoMMatsuuraNExpression and prognostic role of cyclin-dependent kinase 1 (cdc2) in hepatocellular carcinomaOncology596874200010.1159/00001214010895070 ChenYRileyDJZhengLChenPLLeeWHPhosphorylation of the mitotic regulator protein Hec1 by Nek2 kinase is essential for faithful chromosome segregationJ Biol Chem2774940849416200210.1074/jbc.M20706920012386167 QuYLiJCaiQLiuBHec1/Ndc80 is overexpressed in human gastric cancer and regulates cell growthJ Gastroenterol49408418201410.1007/s00535-013-0809-y23591767 BiècheIVacherSLallemandFTozlu-KaraSBennaniHBeuzelinMDriouchKRouleauELereboursFRipocheHExpression analysis of mitotic spindle checkpoint genes in breast carcinoma: Role of NDC80/HEC1 in early breast tumorigenicity, and a two-gene signature for aneuploidyMol Cancer1023201110.1186/1476-4598-10-2321352579 ChenYRileyDJChenPLLeeWHHEC, a novel nuclear protein rich in leucine heptad repeats specifically involved in mitosisMol Cell Biol1760496056199710.1128/MCB.17.10.60499315664 GlinskyGVBerezovskaOGlinskiiABMicroarray analysis identifies a death-from-cancer signature predicting therapy failure in patients with multiple types of cancerJ Clin Invest11515031521200510.1172/JCI2341215931389 LiangYLiuMWangPDingXCaoYAnalysis of 20 genes at chromosome band 12q13: RACGAP1 and MCRS1 overexpression in nonsmall-cell lung cancerGenes Chromosomes Cancer52305315201310.1002/gcc.2203023225332

Weighted gene co-expression network analysis. (A) A hierarchical clustering tree of the differential expression genes set was constructed by average linkage of hierarchical clusters. Gene set is grouped into 10 co-expression modules. (B) Heatmap of correlation between MEs and disease stages. The first value of each rectangle represents Spearman's correlation coefficient between the MEs and stages. The second value in each rectangle is the P-value. (C) Blue spots are genes with |GS|>0.75 and |MM|>0.8 in the blue module. (D) Turquoise spots are genes with |GS|>0.75 and |MM|>0.8 in the turquoise module. In (C) and (D), red spots represent the genes with |GS|<0.75 or |MM|<0.8 in the modules. ME, module eigengene; GS, gene significance; MM, module membership.

Expression patterns of the 44 key genes. (A) Heatmap shows the gradual change in the expression level of key genes between the normal stage and hepatocellular carcinoma. (B) The 44 key genes are grouped into five clusters by self-organizing tree algorithm, and the genes show gradual expression trends along with the disease stages (|r|>0.8, P<0.01). (C) Interaction network of the 44 genes, based on the GeneMANIA online tool.

Differentially co-expressed network in normal-cirrhosis and dysplasia-early HCC stages. (A) Differentially co-expressed network in normal-cirrhosis stages. (B) Differentially co-expressed network in dysplasia-early HCC stages. Turquoise nodes represent key genes; pink nodes represent differential co-expression genes with key genes. Edges represent the association between genes. Solid lines and dashed lines represent enhanced and weakened association between consecutive stages, respectively. HCC, hepatocellular carcinoma.

The interaction networks between HCV viral proteins and their targets in the key genes set. The red nodes in the network are the HCV viral proteins. The green nodes are the proteins targeted by HCV viral proteins. The pink nodes are the other members in the key genes set. The grey edges represent interactions between proteins. HCV, hepatitis C virus.

Significant functions of genes retained in the blue module.

Term Name P-value
GO:0030246 Carbohydrate binding 1.78×10−6
GO:0005529 Sugar binding 1.96×10−5
GO:0016485 Protein processing 7.59×10−4
GO:0051604 Protein maturation 9.72×10−4
GO:0006956 Complement activation 2.09×10−3
GO:0002541 Activation of plasma proteins involved in acute inflammatory response 2.19×10−3
GO:0005509 Calcium ion binding 4.40×10−3
GO:0006959 Humoral immune response 7.21×10−3
GO:0044421 Extracellular region part 7.35×10−3
GO:0005615 Extracellular space 7.91×10−3
GO:0051605 Protein maturation by peptide bond cleavage 8.50×10−3
GO:0009611 Response to wounding 9.71×10−3
GO:0001867 Complement activation, lectin pathway 9.71×10−3

GO, Gene Ontology.

Top 30 significant functions of genes retained in the turquoise module.

Term Name P-value
GO:0000278 Mitotic cell cycle 1.58×10−23
GO:0000279 M phase 6.52×10−23
GO:0022403 Cell cycle phase 1.68×10−22
GO:0022402 Cell cycle process 3.53×10−21
GO:0007067 Mitosis 4.37×10−21
GO:0000280 Nuclear division 4.37×10−21
GO:0000087 M phase of mitotic cell cycle 5.96×10−21
GO:0048285 Organelle fission 8.73×10−21
GO:0007049 Cell cycle 1.35×10−19
GO:0051301 Cell division 6.69×10−19
GO:0015630 Microtubule cytoskeleton 1.60×10−15
GO:0005819 Spindle 4.93×10−15
GO:0044430 Cytoskeletal part 1.94×10−11
GO:0005856 Cytoskeleton 1.07×10−10
GO:0007346 Regulation of mitotic cell cycle 2.78×10−10
GO:0043228 Non-membrane-bounded organelle 6.31×10−10
GO:0043232 Intracellular non-membrane-bounded organelle 6.31×10−10
GO:0051726 Regulation of cell cycle 1.06×10−9
GO:0007093 Mitotic cell cycle checkpoint 1.84×10−9
GO:0007059 Chromosome segregation 2.30×10−9
GO:0000075 Cell cycle checkpoint 1.82×10−7
GO:0007017 Microtubule-based process 4.03×10−7
GO:0000922 Spindle pole 2.59×10−6
GO:0000070 Mitotic sister chromatid segregation 2.92×10−6
GO:0000819 Sister chromatid segregation 3.27×10−6
GO:0005694 Chromosome 4.74×10−6
GO:0051276 Chromosome organization 5.86×10−6
GO:0007051 Spindle organization 7.26×10−6
GO:0007052 Mitotic spindle organization 9.09×10−6

GO, Gene Ontology.

Top 30 significant functions of genes showed enhanced association with key genes in dysplasia-early hepatocellular carcinoma stages.

Term Name P-value
GO:0007049 Cell cycle 4.43×10−21
GO:0000279 M phase 4.62×10−20
GO:0022403 Cell cycle phase 2.29×10−18
GO:0022402 Cell cycle process 4.21×10−16
GO:0000280 Nuclear division 6.09×10−16
GO:0007067 Mitosis 6.09×10−16
GO:0000087 M phase of mitotic cell cycle 7.71×10−16
GO:0048285 Organelle fission 1.03×10−15
GO:0000278 Mitotic cell cycle 1.79×10−14
GO:0051301 Cell division 2.53×10−11
GO:0007059 Chromosome segregation 4.96×10−10
GO:0000775 Chromosome, centromeric region 9.18×10−10
GO:0007017 Microtubule-based process 3.42×10−9
GO:0005819 Spindle 3.56×10−9
GO:0000226 Microtubule cytoskeleton organization 3.28×10−8
GO:0000793 Condensed chromosome 3.88×10−8
GO:0007051 Spindle organization 5.16×10−8
GO:0005694 Chromosome 1.99×10−7
GO:0000779 Condensed chromosome, centromeric region 7.80×10−7
GO:0043228 Non-membrane-bound organelle 9.53×10−7
GO:0043232 Intracellular non-membrane-bound organelle 9.53×10−7
GO:0015630 Microtubule cytoskeleton 1.02×10−6
GO:0044427 Chromosomal part 5.86×10−6
GO:0000777 Condensed chromosome kinetochore 1.60×10−5
GO:0044430 Cytoskeletal part 2.14×10−5
GO:0007010 Cytoskeleton organization 4.63×10−5
GO:0000776 Kinetochore 4.92×10−5
GO:0005876 Spindle microtubule 5.99×10−5
GO:0005524 ATP binding 6.05×10−5
GO:0032559 Adenyl ribonucleotide binding 6.85×10−5

GO, Gene Ontology.