Introduction
Hepatocellular carcinoma (HCC) represents a leading cause of cancer-related death worldwide.1 Currently, though several staging systems can stratify the HCC patients into appropriate risk categories, a great deal of divergence remains within each risk group, due to the molecular heterogeneity of tumor cells and microenvironment. One way to describe the progression of HCC is tumor growth, since it might hold for the long-held view that rapidly growing tumors are more likely to metastasize and become lethal than slow-growing tumors;2 although, the fundamental question regarding the ability of tumor cells to rapidly grow remains to be answered.
In fact, tumor cells adapt to changing environmental conditions and profoundly shape the dependencies of individual cells. For instance, through aerobic glycolysis, cancer cells produce energy by taking up glucose at much higher rates than other cells, while, at the same time, using a smaller fraction of the glucose for energy production. This allows cancer cells to function more like fetal cells, promoting extremely rapid growth.3 However, the underlying molecular basis of the intertwined interactions among tumor immunology, oncogenic signaling, and tissue/biochemical context, upon tumor growth remains largely unknown.4
Currently, a robust gene coexpression network for the mining of hub genes that drive pivotal signaling pathways in terms of large-scale gene expression profiles can be built through weighted gene coexpression network analysis (WGCNA).3 Previous studies have applied WGCNA to provide functional explanations of systems biology, proposing candidate therapeutic targets or diagnostic biomarkers for cancer in recent years.5 For example, Zhao et al.6 utilized WGCNA to investigate the relationships underlying the molecular and clinical characteristics of cholangiocarcinoma. However, a robust WGCNA network for cancer growth has not yet been established.
In the present study, the HCC transcriptome and tumor growth-related modules were explored by WGCNA. Focusing on rapid tumor growth, the integrative functional analysis was expanded to the levels of the growth signature, associated molecular events, and corresponding modulations.
Methods
Data preparation
The transcription profile of HCC was downloaded from the Gene Expression Omnibus (commonly referred to as GEO) with accession number GSE54236, which includes 78 primary HCC tumor samples representing the different speeds of tumor growth:7 Briefly, patients underwent two computed tomography scans 6 weeks apart to determine tumor volumes and HCC doubling time, which ranged from 30 to 621 days and were divided into the following quartiles: ≤53 days (n=19), 54–82 days (n=20), 83–110 days (n=20), and ≥111 days (n=19). Based on these quartiles, tumor growth was classified into slow, fast, faster, and fastest states, respectively. Low and non-expressed genes were removed by selecting probes with a mean expression in the top 50% of all probes. Next, genes with expression variance above average level were selected. Different probes targeting the same gene were collapsed. These steps finally resulted in 5511 genes to infer coexpression networks.
In addition, GSE14520, GSE25097, GSE62232, GSE36376 datasets were obtained from the GEO database. RNA-seq expression profiles from nine cancer types, including adrenocortical carcinoma, kidney renal clear cell carcinoma, kidney renal papillary cell carcinoma, brain lower grade glioma, liver hepatocellular carcinoma (LIHC), lung adenocarcinoma, mesothelioma, pancreatic adenocarcinoma, and sarcoma, were obtained from The Cancer Genome Atlas (commonly known as TCGA) database. Detailed information about the datasets is shown in Supplementary Table 1.
Table 1Hub genes in the HCC rapid growth module
| Top10 GS | GS.Fastest | p.GS.Fastest | Top10 K.in | K.in | Top10 MM | MM.green | p.MM.green |
|---|
| CCNA2 | 0.5795 | 7.97E-16 | CDCA5 | 51.7516 | CCNB1 | 0.9639 | 2.53E-93 |
| NDC80 | 0.5733 | 1.89E-15 | CCNB1 | 50.1223 | CDCA5 | 0.9611 | 8.49E-91 |
| CDCA8 | 0.5700 | 2.99E-15 | PRC1 | 50.1064 | PRC1 | 0.9503 | 1.69E-82 |
| CENPE | 0.5580 | 1.48E-14 | CDK1 | 48.4094 | CDK1 | 0.9492 | 9.73E-82 |
| KIF11 | 0.5576 | 1.55E-14 | EXO1 | 47.6809 | EXO1 | 0.9478 | 7.54E-81 |
| H2AFX | 0.5571 | 1.66E-14 | BIRC5 | 46.9120 | DLGAP5 | 0.9431 | 6.61E-78 |
| RAD54L | 0.5534 | 2.69E-14 | DLGAP5 | 46.5064 | CENPF | 0.9404 | 2.24E-76 |
| NCAPH | 0.5520 | 3.20E-14 | CENPF | 44.2346 | BIRC5 | 0.9354 | 1.11E-73 |
| RFC4 | 0.5495 | 4.42E-14 | PTTG1 | 44.0655 | PTTG1 | 0.9352 | 1.32E-73 |
| HJURP | 0.5491 | 4.61E-14 | CDCA8 | 43.8984 | NUF2 | 0.9350 | 1.83E-73 |
Coexpression network construction
We constructed the coexpression network using the WGCNA package.8 Briefly, the steps included: (1) defining the similarity matrix; (2) choosing the soft threshold, β, and inferring the adjacency matrix; (3) defining the topological overlap matrix; (4) performing hierarchical clustering; (5) performing the dynamic tree cut method to identify the modules; and (6) computing the module eigengene (ME) of each module. The ME can be considered as a representative of the gene expression profiles in a module. The average-linkage hierarchical clustering method was employed to cluster the MEs of all modules, and the modules with high similarity were merged to obtain the coexpression network.9 Another tool, the CEMITool package, was used to validate the gene modules, as described previously.10
The module preservation statistic Zsummary was used to assess the overlap between network modules,8 which takes into account the overlap in module membership (MM), the density (mean connectivity) and connectivity (sum of connections) patterns of modules. A module was considered not being preserved if preservation Zsummary < 2, moderately preserved if 2≤Zsummary<10, and highly preserved if Zsummary ≥10.
Identification of hub genes and growth signature
Hub genes (genes that are highly interconnected with the nodes of the module) are of functional importance. MM was defined as the correlation between the ME and gene expression values. The MM measure is highly related to intramodular connectivity (K.in).8 Highly connected intramodular hub genes tend to have high MM values to the respective module. In short, the larger the MM value of the gene, the higher the correlation between the gene and a given module. In addition, the gene significance (GS) was defined as mediated p-value of each gene (GS=lgP) in the linear regression between gene expression and the clinical traits.5
We used the network screening function based on GS (representing the correlation between the gene and a given clinical trait) and MM or K.in (representing the correlation between the gene and a given module) in the WGCNA package to directly identify the top hub genes in the fastest growth module, and further assembled them into the growth signature. The growth activity was quantified by applying the single-sample gene-set enrichment analysis (ssGSEA).6 We defined the growth signature as either high or low by using median cut-off.
Functional annotation
Functional annotations of the gene sets were performed using webgestalt11 or Enrichr.12
Pathway activity computation
Using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways as templates for signal simulation, a canonical circuit was defined as any possible route the signal can traverse to be transmitted from a particular input to a specific output node.13 The computation of the signal intensity across the different circuits of the entire pathways was performed by the Hipathia program.13 In addition, the growth signature associated cancer pathway activity in pan-cancer was performed using GSCALite.14
Modeling genetic-gene expression
Multi-layered profiles for DNA copy numbers, mRNA expression, and mutations in the LIHC data were obtained from the Xena portal (http://xena.ucsc.edu/ ). A linear modeling approach that measured the association of expression levels was used on a gene-by-gene basis with a number of potential predictors, including gene mutations or genomic alterations, as previously described.15 Somatically acquired mutations and genomic alterations were presented by using Maftools16 and encoded as being present/absent. Linear expression models were fit with the Limma package.17 For the expression of gene k in patient i, Yik is fitted by the following equation:
Yik=∑j=1Xijβjk+β0k+ε
Xij is the mutation matrix for patient i and mutation j, with entries Xij=1 denoting patient i has a mutation j and 0 otherwise. The coefficients βjk measure the expression change in gene k induced by the presence of a mutation j. The entry β0k implies the baseline expression level of gene k.15Immune and metabolism signatures
For each tumor sample, ESTIMATE was used to assess tumor purity.18 The cytolytic activity or the interferon score of the local immune infiltrate was calculated as previously described.19 Gene signatures of 28 tumor-infiltrating lymphocytes, including CD8 or CD4 T cells, B cells, natural killer cells, as well as markers from multiple types of oncoimmunology-containing genes were referenced in a prior study.20 The ssGSEA method was utilized to quantify the enrichment levels of metabolism-related signatures based on Reactome (https://reactome.org ) gene-sets.
Cellular experiments
The human liver cancer cell line HepG2 was purchased from the Cell Bank of Chinese Academy of Sciences (Shanghai, China). Cells were cultured in DMEM (Gibco, Life Technologies, Waltham, MA, USA) supplemented with 10% fetal bovine serum (Gibco), and used within 20 passages after culture.
Lentivirus production was performed as previously described.21 The targeting sequences of short hairpin (sh)RNAs were as follows: 5′-CACCGAACCAACCAACTTTAA-3′ for sh-NCAPH#1, 5′-ACTGACTCACCTCGCTTATTG-3′ for sh-NCAPH #2; and 5′-CTTTGTAATCATCCAGCTCTA-3′ for sh-RAD54L#1, 5′-TGGTCTGGGTGTAGCTCTTAG-3′ for sh-RAD54L#2. Knockdown efficiency was verified by quantitative PCR. The Cell Counting Kit-8 (Dojindo Molecular Technologies, Kumamoto, Japan) assays were used to measure cell proliferation.
Statistical analysis
Experimental data are represented as the average ± standard deviation. Unless otherwise indicated, the Student’s two-tailed unpaired t-test was used to determine statistical significance. The significance threshold was set at 0.05. For survival analysis, the LIHC data were analyzed with the GEPIA database.22 Samples were split into high and low groups based on the median value. Kaplan-Meier survival analysis was calculated using the log rank test, with a p value for significance of <0.05.
Results
Rapid growth module in HCC
In defining the gene clusters involved in HCC growth, 11 distinct gene modules were explored using WGCNA, as shown in Fig. 1A.
Module stability was verified by repeating network construction and module identification on expression data that consists of resampled sets of the original dataset or alternatively by another tool (the CEMITool) running the same dataset10 (Supplementary Fig. 1). The results proved the robustness of module assignments.
Next, we evaluated the relationship between each module and the growth status. Notably, the cirrhotic features of modules (red, blue, yellow, brown) and the aggressive proliferative HCC features of modules (magenta and green for the faster and fastest state respectively) were identified (Fig. 1B and C). Network features such as GS, MM and K.in23 of each module in different growth states were computed (Supplementary Table 2).
Moreover, we narrowed down our analyses to the fastest HCC growth. Module green (referred to as the growth module hereafter) was the best, as reflected by its strongly positive correlations to the fastest tumor growth (Fig. 1B–D).
Signal fluxes in the rapid growth module
To ask the question of whether the growth module was highly preserved across independent HCC datasets, external validation was performed. Using five different HCC cohorts, the growth module showed a higher preservation statistics summary than expected by random chance using bootstrapping procedures (Fig. 2A). Thus, the growth module was deemed to hold promise in independent tumor profiles from different patient cohorts. For function enrichment, the cut-off set with the false discover rate of <0.01, cell cycle, DNA replication, Fanconi anemia, etc. constituted the main KEGG pathways in the growth module (Fig. 2B).
To decompose the KEGG pathway into detail, a canonical circuit was defined as any possible route the signal can traverse to be transmitted from a particular input to a specific output node.13 Effector nodes at the end of the circuits trigger specific functions in the cell. Using gene expressions as proxies of node activation values, computation of the signal intensity across the different circuits of the pathways was performed by canonical circuit activity analysis to compute the transmission of the signal along the network.13 Thus, we estimated the level of activity of subpathways (signaling circuits) using the Hipathia program,13 and detected several pathways with perturbed activity in the growth module (Supplementary Fig. 2).
Focusing on the cell cycle pathway, five effector circuits were deemed ultimately responsible for the functions of DNA replication and cell cycling. These circuits were highlighted in the fastest state, one of them ending in the node including RB1, one including RAD21, one containing TFDP1 and E2F4, one ending in the node with protein genes for CDC45, MCM7, MCM6, MCM5, MCM4, MCM3 and MCM2, and the last one ending in the node with protein genes for ORC and MCM (Fig. 2C). Indeed, most nodes in the effector circuits have adverse outcome in the LIHC cohort (Supplementary Fig. 3). These results clearly suggested the signaling in the green module as providing multiple routes and broader activity to promote cell cycle progression, thus accelerating tumor growth.
Rapid growth program in HCC
Next, we examined the hub genes in the growth module. Hub genes,5 including the top 10 GS, K.in or MM genes, are shown in Fig. 3A and Table 1. Among them, 20 genes were identified.
There was step-wise activation of all these genes that accompanied increased speed of tumor growth (Fig. 3A and Supplementary Fig. 4). The trend was clearly consistent and coordinated. As expected, all these hub genes were involved in the advanced prognosis of HCC, as evidenced by the results from our survival analysis (Supplementary Fig. 4).
Since the top 20 hub genes were densely interacted by protein-protein interaction analysis,24 we categorized them among the rapid growth signature (Fig. 3B) and applied the ssGSEA algorithm to infer the growth activity for each sample.
Recent advances in molecular biomarkers of HCC have indicated various oncogenes and tumor suppressor genes.25 Indeed, growthhigh patients in our study were more likely to show higher expression of many known adverse prognostic biomarkers, such as AFP, DCP1A, GPC3, MDK, MCM6. In contrast, growthlow patients were likely to show higher expression of tumor suppressor genes, such as GPR155 and IFIT3 (Fig. 3C).
Furthermore, we found that rapid growth activity has a bad survival prognosis, both in LIHC and GSE14520 cohorts (Fig. 3D). These findings suggested a negative regulation relationship between the HCC growth program and HCC prognosis.
Next, we reasoned that genes in the signature would have a higher degree of association with cell proliferation. To test this hypothesis, we selected genes (NCAPH and RAD54L) for experimental validation. As expected, knockdown of NCAPH or RAD54L significantly suppressed the proliferation of HepG2 cells (Fig. 3E).
Somatic mutations and copy number alterations of the growth signature
We next investigated the growth signature at the genomic level. By focusing on the somatic non-silent mutation or copy number variation (commonly known as CNV) genes (Fig. 4A), the principal components analysis (commonly known as PCA) was computed to maximize the stability of the components. The first two principal components, respectively, account for 13.6% and 11.2% of the total variability in gene expression; the first 20 principal components cumulatively explain 67% of the variance (Fig. 4B). Notably, overlaying the status of the genetic alterations of hub genes (growthmut) on to the first two principal components demonstrated that mutations or CNV alterations correlated with general gene expression profiles (Fig. 4C).
The transcriptome was globally perturbed by growthmut, with expression levels of 6,565/15,569 (42%) genes significantly associated with at least one genetic change of the hub genes (r>0.3, false discovery rate of <0.001). For instance, genetic change of NUF2 co-occurred with other hub genes’ alteration, which led to 1069 genes’ differential expression (Fig. 4D–E). The observed variability can be largely explained by the presence of other hub genes’ alteration leading to strong up-regulation of NUF2 mRNA (Fig. 4D). The expression changes of the growthmut-related genes are summarized in Supplementary Table 3.
To understand the mechanism underlying growthmut, webgestalt11 analysis showed that growthmut-related genes were enriched in the KEGG pathways of cell cycle, oocyte meiosis and DNA replication, etc. (Fig. 4G). Given that oncogenic signaling accelerates cell cycle progression,26 these data indicated that growthmut cancers can amplify growth signaling to maintain cell proliferation.
Rapid growth is associated with tumor immune-metabolism
We then explored the correlation between growth activity and tumor-infiltrating lymphocytes (Fig. 5A). We observed a significant negative correlation of growth signature with Th1 cell types and positive correlation with Th2 cell types. Next, significant negative correlations between rapid growth and natural killer cells, plasmacytoid dendritic cells or macrophages were found (Fig. 5A and Supplementary Table 4).
Moreover, dysregulation of diverse immune signatures in HCC, including HLA expression, cytokines or chemokines, and interferon response were identified between growthhigh and growthlow groups (Fig. 5C). Notably, the growth signature showed strong correlation with neoantigens (r=0.37, p=1.52E-12 in LIHC; r=0.38, p=5.22E-09 in GSE14520) (Fig. 5B). In addition, tumor rapid growth showed no significant correlation with immunoinhibitors or immunostimulators, including well-known checkpoint genes (Table S4). Additionally, immune cytolytic activity or tumor burden was not significantly different between high versus low rapid growth tumors (Supplementary Fig. 5). Thus, the growth signature itself might not elicit the active immune response.
Next, we investigated the metabolic configuration according to the growth activity. Using curated metabolic gene sets from Reactome as indicators, we found rapid tumor growth was highly associated with diverse metabolism processes. For example, the rapid tumor growth was significantly positively correlated with synthesis of DNA (r=0.73, p=1.37E-58 in LIHC; r=0.85, p=2.52E-64 in GSE14520) and negatively correlated with bile acid metabolism (r=−0.43, p=1.70E-17 in LIHC; r=−0.46, p=1.92E-13 in GSE14520). Detailed growth-metabolism correlations are provided in Supplementary Table 5.
The representative metabolic activity between growthhigh and growthlow tumors is provided in Fig. 5D. For instance, dysregulation of cytochrome P450, xenobiotics, and biological oxidation are shown to be associated with poor prognosis. These results suggested that growthhigh tumors tend to present an immunetolerant and metabolism reconfigured microenvironment in HCC.
Prognostic role of growth signature in pan-cancer
When extending the growth signature to the pan-cancers, we found a significant hazard ratio between the growth genes and overall survival or recurrence-free survival in multiple cancers, including adrenocortical carcinoma, kidney renal clear cell carcinoma, kidney renal papillary cell carcinoma, brain lower grade glioma, LIHC, lung adenocarcinoma, mesothelioma, pancreatic adenocarcinoma, and sarcoma (Fig. 6A). Gene expression was significantly higher in the growthhigh group (Fig. 6D) and the growthhigh indicated advanced prognosis in these cancers (Fig. 6B–C). The pathway relation network also indicated that the signature was mostly involved in the cell cycle in nine cancers (Fig. 6E).
To investigate the clinical implications of the growth signature, we searched for targets of candidate drugs by using the L1000 project.27 The top 10 associations are presented in Fig. 6F. Palbociclib functions as a CDK4/6 inhibitor in multiple cancers.28 It also acts as a novel radiosensitizer, dovitinib, inducing apoptosis and overcoming sorafenib resistance in HCC through SHP-1-mediated inhibition of STAT3.29–32 As the topoisomerase inhibitor, mitoxantrone, was widely used in clinic.33 Canertinib is an irreversible EGFR inhibitor in cancers.34
Discussion
The theory that patients with a rapidly growing cancer have a poor prognostic outlook has remained persistent;2 thus, detailing the biological structure in the hidden layer of rapid growth is an interesting question.
As illustrated in this study, to understand the growth rate of HCC at the modular level and toward uncovering the critical drivers of the disease in a comprehensive manner, HCC transcriptomes were explored in the context of modular pattern, in which the green module was responsible for rapid tumor growth after internal and external validation. Briefly, the top enriched functional classes of this program are consistent with our existing knowledge of the cell cycle, as well as DNA repair, replication and cell proliferation in cancer.
Due to the presence of the highly interconnected top 20 hub genes, we assembled them into the growth signature. At the gene level, we demonstrated that the expression levels of all these genes were increased in coordination with the state of growth-rate. Accordingly, high expression of these genes predicted the adverse outcome of HCC. Indeed, various lines of evidence showed the involvement of these hub genes in previous studies. For example, CCNA2 is the leading gene according to the GS rank. A recent study revealed a new poor-prognosis HCC entity and a rearrangement signature related to replication stress, due to CCNA2 alterations.35 The top K.in gene, CDCA5, transcribed by E2F1, promotes oncogenesis by enhancing cell proliferation and inhibiting apoptosis via the AKT pathway in HCC.36 The top MM gene, CCNB1 was highly expressed in the samples of recurrent HCC, which was associated with significantly reduced recurrence-free survival.37 In addition, the vast majority of genes within the signature, such as RFC4, HJURP, ECT2, KIFC1, NUSAP1, CDK1, PRC1, KIF4A, etc., contribute to the pathogenesis of HCC, as reported previously.38–43 At the time of preparation of this paper, little is known about NCAPH and RAD54L in HCC. In this study, however, we found that knockdown of NCAPH and RAD54L expression is associated with growth inhibition.
Genetic background and in particular genomic alteration might contribute to gene expression in HCC related to tumor growth speed. Taking a genetic-centric view, we specified the set of gene expression changes that correlate with the alteration of the signature. The number of target genes whose expression are differentially affected varies widely across the different genetic change of growth signature. For example, BIRC5, playing dual roles in mitosis and cell survival of HCC,44 was independently correlated with expression levels of 920 genes. These data supported the notion that genetic change of these hub gene results in a rapid tumor growth.
Moreover, the growth signature exhibited a significant correlation with certain genomic features, such as tumor purity and neoantigens. Further, no obvious association between rapid growth and active antitumor immune signatures was found. In addition, a significant negative correlation was observed between growth activity and energy metabolism integration, biological oxidation, vitamins, fatty acid, and glucose metabolism, etc. (Supplementary Table 5). The above results supported the idea of a link between tumor growth, metabolic landscape reconfiguration, and clinical progression of cancer. Accordingly, aggressive malignant phenotypes of cancer cells obtained by accumulated mutations change metabolic phenotypes for proliferation represented by nucleotide and amino acid metabolism. Indeed, this hypothesis of highly coordinated growth program warrants further study.
Finally, candidate drugs have been inferred based on the growth signature. Numerous cell cycle inhibitors have been designed over the past decades.45 The current findings of growth-specific drugs for HCC would have potential implications in warranting future studies toward developing targeted and combinatorial therapeutics for HCC.
Supporting information
Supplementary Table 1
HCC data sets.
(DOCX)
Supplementary Table 2
Summary of WGCNA and CEMITool analysis.
(XLSX)
Supplementary Table 3
Growthmut-affected genes.
(XLSX)
Supplementary Table 4
Correlation between growth signature and diverse immune-related signatures.
(XLSX)
Supplementary Table 5
Correlation between growth signature and diverse metabolism-related signatures.
(XLSX)
Supplementary Fig. 1
Verification of module assignment.
(A) Resampling survey on the full data set, including the dendrogram of the genes together with module assignment in the full data as well as in the 50 resampled data sets. Identified modules are remarkably stable. (B) Dynamics of module significance of each module from cirrhotic, slow to fastest growth state of HCC. (C) Module assignment by using WGCNA or CEMITool respectively. (D) The left panel indicated the modules identified by the CEMITool using adjacent cirrhotic tissues and tumor tissues. Focusing on tumors with different growth rates, the corresponding modules are presented on the right panel. (E) Overlap in gene compositions between WGCNA and CEMITool-identified modules, e.g. module green (WGCNA) with M5 (CEMITool) and module magenta (WGCNA) with M9 (CEMITool). (F) Significance of pairwise module-module overlap between two methods using Jaccard index or odds ratio respectively.
(TIF)
Supplementary Fig. 2
Signal transduction route of KEGG pathways enriched for the rapid growth module, including cell cycle, Fanconi anemia pathway, oocyte meiosis pathway, human T lymphotropic virus-I infection pathway, TP53 signaling and the pathway in cancer.
Arrows in red represent activated circuits. Genes in red represent genes up-regulated in the cancer with respect to the adjacent normal tissues; genes in blue represent down-regulated genes and genes with no color were not differentially expressed. Squares at the end of the circuit represent the cell functions triggered by the circuits. The ending nodes of the effector circuits are highlighted.
(TIF)
Supplementary Fig. 3
Survival contribution of multiple ending nodes (genes) of the activated circuits in the cell cycle pathway, estimated using the Mantel-Cox test.
(TIF)
Supplementary Fig. 4
Top hub genes in the green module.
(A) Expression level of hub genes from cirrhotics, from the slow to the fastest state. (B) Kaplan-Meier plots and log-rank tests for the hub genes that are associated with the overall survival of HCC patients in LIHC cohort, generated by GEPIA.
(TIF)
Supplementary Fig. 5
Tumor burden or immune cytolytic activity was not significantly different between the growthhigh and growthlow groups.
(TIF)