Integrative Transcriptome Analysis Reveals Common Molecular Subclasses of Human Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) affects approximately half a million patients worldwide and is the most rapidly increasing cause of cancer death in the United States owing to the lack of effective treatment options for advanced disease (1). Numerous lines of clinical and histopathologic evidence suggest that HCC is a heterogeneous disease, but a coherent molecular explanation for this heterogeneity has yet to be reported. Genomic approaches to the classification of HCC therefore hold promise for a molecular taxonomy of the disease.

Mutations in the WNT signaling pathway have been found to be common in HCC, but other DNA level classification approaches have proven challenging. This relates to the enormous complexity of the genomic alterations observed in HCC, likely attributable to the accumulation of chromosomal rearrangements resulting from decades of chronic viral hepatitis and cirrhosis. This complexity makes it difficult to identify the causal genetic events promoting HCC development and progression (2, 3). An alternate approach to HCC classification has been to study tumors at the level of their gene expression profiles. Although several such profiling efforts have been reported (4–11), a cohesive view of expression-based subclasses of HCC has yet to emerge. In part, this is because each of the reported studies analyzed different patient populations (most of them small) on a different microarray platform, with a different primary biological or clinical question in mind. Perhaps not surprisingly then, each study reported a somewhat different view of the heterogeneity of HCC, and it has been therefore impossible to see whether there exists a common biological thread that links these disparate studies.

We believe that any biologically or clinically meaningful classification system should be informative across multiple patient populations and should be independent of any particular microarray platform. In the present study, we therefore set out to define molecular subclasses of HCC that existed across all available HCC data sets, including eight previously reported studies and one new one reported here, totaling 603 patients. We report that indeed there exist three distinct molecular subclasses of HCC that are present in all nine data sets examined, regardless of technical differences between the microarray platforms used to generate the profiles. We show that these subclasses are correlated with histologic, molecular, and clinical features of HCC, and we highlight the important role of transforming growth factor-β (TGF-β) signaling in one of the HCC subclasses. These findings thus create a solid foundation for HCC classification on which to build informed clinical trials for patients with HCC and also suggest new opportunities for therapeutic intervention.

Identification of common HCC subclasses. To define and validate a gene expression–based model of common molecular subclasses of HCC, we collected publicly available gene expression data sets from eight independent cohorts profiled on a wide variety of microarray platforms (HCC-A, HCC-B, HCC-C, HCC-D, HCC-E, HCC-F, HCC-G, and HCC-H; see Supplementary Tables S1 and S2 for details; refs. 4–11). Between the training data sets (HCC-A, HCC-B, and HCC-C) chosen as larger data sets covering major etiologies of HCC to avoid overfitting a model to any particular cohort or microarray platform, corresponding subgroups of the samples were defined by subclass mapping (SubMap) method (12) based on subclasses that identified three unsupervised clustering methods: hierarchical clustering, k-means clustering, and nonnegative matrix factorization, which finds clusters of samples after collapsing the data set into representative “meta-genes” (13).

For each subclass defined by SubMap, meta-analysis marker genes were selected as overexpressed genes compared with the rest of the subclasses (HCC subclass signature) in all the three clustering methods to avoid defining gene expression–based subclasses that were unique to a particular clustering algorithm. Prediction of the subclass was performed for each sample using nearest template prediction method (14, 15) to accommodate the diverse microarray platforms (see Supplementary Data for details).

Molecular annotation of HCC subclasses. Functional characterization of the HCC subclasses was performed using gene set enrichment analysis (GSEA; ref. 16). Two categories of gene sets in Molecular Signature Database (MSigDB)¹¹

Clinical data analysis. The hazard of tumor recurrence was calculated to estimate the pattern of HCC recurrence over time after the surgery as previously reported (17, 18). Continuous and proportional data were tested by Wilcoxon rank sum test and Fisher's exact test, respectively. All clinical data analyses were performed using the R statistical package version 2.4.0.¹²

We generated a ninth data set of formalin-fixed, paraffin-embedded (FFPE) tumors (HCC-I), reasoning that any meaningful classification system should be applicable to routinely collected fixed (as opposed to frozen) tissues. We analyzed FFPE tissue blocks from 118 HCC patients who consecutively underwent surgical resection during 1990 to 2001 at Toranomon Hospital. Ethical approval for use of the FFPE tissues, obtained and archived as part of routine clinical care, was acquired from the institutional review board granted on condition that all samples be made anonymous. Total RNA was extracted from macrodissected 10-μm tissue slices (three to four slices for each sample) using High Pure RNA Paraffin kit (Roche). Expression of transcriptionally informative 6,000 genes, selected to capture global state of transcriptome (14), was profiled using DNA-mediated annealing, selection, extension, and ligation (DASL) assay (Illumina; see Supplementary Data; ref. 19).

Total RNA was isolated using Trizol reagent (Invitrogen) according to the manufacturer's instruction. Microarray experiment was performed using HT_HG-U133A High-Throughput Arrays (Affymetrix). The raw data were normalized using BioConductor affy package.¹³

All microarray data sets are available through National Center for Biotechnology Information Gene Expression Omnibus¹⁴

Immunohistochemical staining was performed on 10-μm FFPE sections using antibodies for β-catenin (Becton Dickinson), phospho-AKT (Cell Signaling), and p53 (Immunotech) followed by detection using the Envision+ DAB system (Dako). The stains were evaluated by a pathologist blinded to the results of the gene expression profiling, and the results were scored in a binary system. For immunofluorescence staining, cells grown on multiwell chamber slides were fixed by 4% paraformaldehyde and stained for β-catenin (see Supplementary Data for details).

Human HCC cell lines, Huh-7 (Riken Bioresource Center) and SNU-387 (American Type Culture Collection), were grown in DMEM supplemented with 10% heat-inactivated fetal bovine serum at 37°C in a 5% CO₂ atmosphere. No β-catenin mutation has been found in these cell lines.

Cell lysates were separated on NuPAGE 4% to 12% gels (Invitrogen) and transferred to polyvinylidene difluoride membranes (Bio-Rad) and blotted for α-fetoprotein (AFP; Santa Cruz Biotechnology), β-catenin, phospho-SMAD3 (Cell Signaling), and proliferating cell nuclear antigen (Santa Cruz Biotechnology) antibodies.

Cells were infected with the indicated short hairpin RNA (shRNA) vectors (construct 1, TRCN0000003843; construct 2, TRCN0000003844; The RNAi Consortium¹⁶

Cells were transfected using Lipofectamine 2000 (Invitrogen) with either TOP-flash or FOP-flash constructs and stimulated with 100 pmol/L TGF-β (R&D Systems) for 48 h. Luciferase activity was measured using Dual-Glo kit (Promega). All transfections were performed in triplicate and measurements were normalized to a SV40 promoter-driven Renilla luciferase construct.

Three common molecular subclasses of HCC. The SubMap method identified three robust subclasses in each of the three initial data sets analyzed. We refer to these subclasses as S1, S2, and S3 (Fig. 1), and the complete list of genes that correlated with each of the subclasses is available in Supplementary Table S3. As with any unsupervised clustering–based definition of cancer subclasses, it is essential to establish the validity of the newfound classification system. In the following sections, we describe three independent validations of the three-class structure of HCC. First, we show that the three subclasses are detected with statistical significance in each of the six remaining HCC data sets (totaling 371 patients). Second, we show that the subclasses are associated with clinical parameters. Third, we show that the subclasses are associated with biological mechanism known to be operative in the pathogenesis of HCC.

Statistical validation of subclasses across nine HCC cohorts. As a statistical measure of the validity of the three subclasses, we determined the confidence with which HCC samples could be classified into one of the three subclasses using the HCC subclass signature–based classifier, including 619 genes. As expected, subclass labels were assigned with high confidence [false discovery rate (FDR) < 0.05] to 94% of the training samples (HCC-A, HCC-B, and HCC-C), which were used to define the subclasses in the first place. More importantly, high confidence subclass labels were assigned to 84% of the 371 samples in the validation set (HCC-D, HCC-E, HCC-F, HCC-G, HCC-H, and HCC-I). In contrast, a classifier based on the same number of genes chosen randomly yielded high confidence predictions in <1% of the samples. In addition, our classification system was superior to those reported previously (10, 11, 20, 21) when those classifiers were tested across all of the validation data sets (high confidence predictions using reported signatures were 27–75%; Supplementary Fig. S1). These results, taken together, indicate the statistical significance of our three subclasses and point to the limitation of defining subclasses based on only a single cohort, where overfitting often leads to failure of the classifier to validate on new samples, particularly when profiled on a different microarray platform.

Clinical relevance of HCC subclasses. Having established the statistical validity of the HCC subclasses, we next asked whether any of the subclasses were associated with clinical parameters to add the validity of the subclasses. Clinical data were available for 197 patients, as summarized in Table 1.

Our first observation was that tumors in subclass S2 were larger than the others, whereas tumors in S3 were smaller compared with the rest (P = 0.003). In addition, subclass S3 included the majority of well-differentiated tumors (37 of 49; P < 0.001), whereas there was no histologic difference between S1 and S2 (P = 0.73). We also examined the serum levels of the one clinically used serum biomarker of HCC—AFP. Serum AFP was the highest in S2 (0.001), further supporting the notion that our subclasses are clinically relevant.

Next, we sought to determine whether the HCC subclasses were associated with clinical outcome following surgical resection. We were careful to analyze the two major patterns of HCC recurrence: early recurrence, which is related to residual dissemination of primary tumor cells within the liver, and late recurrence, which is attributable to new primary tumors arising in a hypercarcinogenic state of a cirrhotic liver (17, 18, 22). “Early” recurrence is associated with malignant characteristics of the primary tumor itself, and we reported that it has less effect on patient survival in earlier-stage HCC, in which “late” recurrence is the major determinant of survival (14, 23). We found that subclass S1 was associated with a significantly greater risk of earlier recurrence (P = 0.03 within 1 year; Supplementary Fig. S2). This association remained to be significant even after correcting for tumor size in multivariate Cox regression modeling (Supplementary Table S4). Consistent with this observation, S1 tumors exhibited more vascular invasion and satellite lesions (both known risk factors for early recurrence; Supplementary Table S5). These results may suggest that the S1 subclass is associated with a more invasive/disseminative phenotype. Interestingly, we found that a recently reported signature of poor survival defined in patients with more advanced HCC, where early recurrence is the major determinant of survival (4), was associated with S1 and S2, whereas the good survival signature defined in that study was enriched in S3 tumors (Supplementary Table S6), lending further credence to our subclass model. Importantly, our HCC subclasses were not associated with late recurrence, consistent with our recent study indicating that late recurrence is determined not by the characteristics of the tumor itself but rather by the biological state of the surrounding liver at risk (14). Furthermore, consistent with our previous observation, the HCC subclasses showed no association with survival (P = 0.12) in our data set (HCC-H) consisting mostly of early-stage tumors.

Molecular pathways associated with HCC subclasses. We next asked whether we could ascribe biological meaning to our validated HCC subclasses. The GSEA results indicated that indeed our HCC subclasses were associated with distinct biological processes, several of which have been implicated in HCC pathogenesis (Supplementary Tables S7 and S8). For example, S2 were tumors associated with a relative suppression of IFN target genes (7), of interest because of the use of IFN as an experimental chemopreventive strategy for HCC. S2 tumors were also enriched in MYC target genes, suggesting that MYC activation is a feature of S2 tumors. Consistent with this observation, we found that a recently reported mouse model of HCC based on MYC overexpression (24) exhibited the S2 subclass signature (Supplementary Fig. S3). S2 tumors were also strongly enriched in a signature of EpCAM positivity (Supplementary Table S6; ref. 25), and in addition, we found that S2 tumors overexpressed AFP (consistent with S2 patients having elevated serum AFP levels; Table 1). Lastly, S2 tumors were enriched in a signature of AKT activation (10), and validation experiment indicated a trend toward elevated phosphorylation of AKT as determined by immunohistochemistry (P = 0.07; Fig. 2A). An AFP-AKT association has been previously observed (10, 26), and we see here that this association is being driven primarily by S2 tumors. The mechanism of AKT activation in these tumors is not known but likely reflects upstream signaling of the phosphatidylinositol 3-kinase (PI3K) pathway (27). As PI3K inhibitors are now entering clinical development, it may be of value to examine their role in S2 tumors in particular.

GSEA also identified differential activation of p53 and p21 target gene sets, with these genes being more abundantly expressed in S3 tumors compared with S1 and S2, consistent with our observation that S3 tumors tend to be lower grade (Table 1). To further validate this result, we performed immunohistochemistry for p53, wherein nuclear accumulation of p53 protein is well known to reflect inactivating p53 mutation (28). As predicted by the GSEA analysis, S1 and S2 tumors exhibited significantly greater nuclear p53 staining compared with S3 (P = 0.001; Fig. 2B). The more well-differentiated nature of S3 tumors was also reflected in the S3 gene expression profile, with S3 tumors exhibiting relatively higher levels of expression of hepatocyte function-related genes involved in glycogen/lipid/alcohol metabolism (APO/ALDH/ADH family genes), detoxification (CYP family genes), coagulation, and oxygen radical scavenging (CAT, SOD1; Supplementary Tables S3, S7, and S8).

WNT pathway activation in S1. The WNT signaling pathway is perhaps the best characterized oncogenic pathway in HCC, with pathway activation occurring through β-catenin mutation (specifically via mutation in exon 3 in up to 44% of cases) and less frequently in AXIN1 (<10% of cases; refs. 2, 3). We addressed WNT status with regard to HCC subclass in two ways. First, we performed GSEA analysis using an experimentally defined WNT activation signature. We found strong enrichment of the WNT signature in subclass S1 compared with S2 or S3 (FDR = 0.03; Supplementary Table S7), suggesting preferential WNT activation in S1 tumors. We validated this result via immunohistochemistry for β-catenin (the principal downstream effector of WNT in HCC) and found that S1 tumors indeed had higher levels of cytoplasmic β-catenin protein expression compared with the other subclasses (0.001), again indicating preferential activation of the canonical WNT pathway in S1 (Fig. 2C; ref. 29). In addition, we found that shRNAs targeting β-catenin resulted in growth inhibition when introduced into the SNU-387 cell line (predicted to be subclass S1), thereby further supporting the hypothesis that WNT activation is functionally important in S1 tumors (Fig. 2D).

Mechanisms of WNT activation in S1 tumors. Having determined that S1 tumors exhibit preferential activation of the WNT/β-catenin pathway, we next addressed potential mechanisms for this activity. We first asked whether the S1 tumors were associated with β-catenin mutation in HCC-H data set, for which we previously sequenced exon 3 of β-catenin gene (11). Surprisingly, β-catenin mutations were preferentially found in S3 tumors, consistent with previously reported “CTNNB1” class representing a subset of S3 (11, 30). This result is also consistent with recent evidence indicating that β-catenin mutations do not regulate the canonical WNT target genes (e.g., cyclin D1 and MYC) that characterize our S1-associated WNT activation signature (31). These results further suggest that the WNT pathway is activated in S1 tumors by a mechanism other than β-catenin mutation.

To explore alternate explanations for WNT pathway activation, we again turned to GSEA, asking whether other gene sets (signatures) enriched in S1 tumors might provide insight into WNT activation in these tumors. Strikingly, we observed strong overexpression of TGF-β target gene sets (i.e., genes expressed as a result of experimental activation of TGF-β) in S1 tumors (Supplementary Table S7). We similarly observed enrichment of a gene set associated with epithelial-to-mesenchymal transition, a phenomenon implicated in tumor invasion and metastasis (32) and known to be regulated by TGF-β signaling in HCC (Supplementary Table S7; refs. 33, 34). Furthermore, a previously reported TGF-β activation signature associated with an invasive phenotype (35) showed strong enrichment in S1 (Supplementary Table S6). We observed no genomic copy number change associated with S1 subclass in TGFB1 locus, suggesting that chromosomal aberration is not the causative mechanism of the activation (Supplementary Fig. S4). These results indicate that TGF-β and WNT signaling co-occur in the same HCC subclass (subclass S1), and suggest the hypothesis that TGF-β might in some way lead to WNT activity that defines the S1 molecular phenotype.

TGF-β-WNT interactions. We next explored the hypothesis that TGF-β regulates WNT pathway activity in HCC cells. First, we treated the HCC cell line Huh-7 with intact WNT pathway components (classified as subclass S2 and with no activation of S1 and WNT activation signature) with TGF-β and monitored the genome-wide expression consequence. As predicted, TGF-β stimulation induced expression of WNT target genes (FDR < 0.001; Fig. 3A) and induced the expression of genes characteristic of subclass S1 (FDR = 0.04; Fig. 3B; Supplementary Data) characterized by WNT/TGF-β activity while suppressing expression of AFP protein, one of the top markers for S2 (Fig. 3C).

Second, we asked whether TGF-β could regulate the activity of a T-cell factor-lymphoid enhancer factor (TCF-LEF) reporter, further reflecting WNT/β-catenin activity. TGF-β stimulation of Huh-7 cells resulted in activation of a wild-type (TOP-flash) but not mutant (FOP-flash) TCF-LEF luciferase reporter (Fig. 4A). Interestingly, the superactivation of TCF-LEF activity was also observed in the presence of cotransfected mutant β-catenin. These results validate the hypothesis that TGF-β enhances WNT activity in HCC, consistent with the subclass S1 molecular profile.

We next explored the mechanism by which TGF-β augments WNT/β-catenin activity. A simple explanation would be that TGF-β induced expression of β-catenin RNA or protein levels, but we found no evidence for this (Fig. 4B). Strikingly, however, TGF-β treatment resulted in a marked change in β-catenin subcellular localization. Specifically, TGF-β treatment induced a shift from membranous β-catenin staining to a cytoplasmic distribution with focal perinuclear aggregation (Fig. 4C). This suggests that TGF-β enhances WNT signaling by modulating the intracellular pool of free β-catenin.

Taken together, these results validate the observation that TGF-β and WNT activity together typify the S1 subclass of HCC, and further suggest that TGF-β augments WNT activity via alteration of the subcellular localization of β-catenin, consistent with the cross-talk between these pathways observed in other biological contexts (34, 36). This implies that therapeutic cotargeting TGF-β and β-catenin in S1 tumors might be explored as a strategy for the treatment of S1 subclass HCC.

Advances in genome technologies are now supporting a breadth of cancer genome characterization studies, including those focusing on HCC. Along with this proliferation of studies has come, however, a certain confusion in the field—different studies often report different results relating to the same set of underlying questions. For example, ∼10 articles on the gene expression–based classification of HCC have been published in recent years, but a consensus molecular taxonomy of the disease has yet to emerge. This might lead some to believe that either expression technologies are insufficiently stable or HCC is so hopelessly heterogeneous and complex that regular, reproducible patterns in the data are nonexistent. We report here that in fact a highly reproducible molecular architecture of HCC is identifiable and is detected across all available HCC data sets.

Our analysis of nine HCC data sets totaling 603 patients indicated that there exist three major subclasses of HCC, which we refer to as subclasses S1, S2, and S3 (Fig. 5). Importantly, although the proportion of each subclass varied slightly from study to study, the subclasses were identifiable regardless of the geographic location of the study patients (Asia versus Europe versus United States) or the technology platform used (cDNA versus oligonucleotide arrays, and frozen versus FFPE tissues). Notably, the new data set generated in the present study used FFPE tissues, thereby showing that the three-class structure is readily detectable in specimens collected and stored in the routine clinical setting. This is relevant because the future deployment of diagnostic tests aimed at cancer classification should ideally be applicable to the standard FFPE specimens that are obtained in clinical practice.

Several biological insights can be made from the observed three-class structure of HCC. Class S1 is particularly notable for the prominence of a WNT activation gene expression signature. This is notable because such WNT activation is not simply explained by the presence of activated β-catenin mutations, suggesting that additional mechanisms of WNT activation seem to be at play, including TGF-β activation. This may be particularly important in the setting of clinical trials testing β-catenin inhibitors in HCC. Our data suggest that such inhibitors may be worth exploring in HCC beyond those patients harboring β-catenin mutation. Although additional mechanistic studies are clearly required, our data support the existence of an interaction between WNT activation and TGF-β activation in S1 tumors, an interaction that has been recently proposed in HCC (34).

Class S2 tumors were notable for their high level of expression of AFP, associated with elevated plasma levels of AFP protein compared with non-S2 patients. S2 tumors also tended to be enriched in MYC tumors harboring a MYC activation signature. This is of relevance because it suggests that genetically engineered mouse models of HCC based on MYC activation may be used to interrogate biological basis of the S2 subclass of human HCC. In addition, the finding of an AKT activation signature in S2 tumors suggests that AKT or PI3K inhibitors might be particularly worthy of exploration in this subclass. Further studies are required to establish the mechanism by which AKT is activated in these tumors.

S3 tumors were notable for their relative histologic evidence of differentiation, and the S3 gene expression program was accordingly suggestive of a molecular program of differentiated hepatocyte function. It is tempting to speculate that these tumors might be particularly well suited to differentiation therapy with agents such as retinoids, as has been previously suggested (37). Whether S3 tumors have distinct mechanisms of transformation or rather simply allow for more complete cellular differentiation remains to be determined. The preserved p53 function in S3 suggests that the abrogation of p53 is associated with stepwise malignant transformation of well-differentiated tumors rather than initiation of carcinogenesis. The less frequent β-catenin mutations in S1 and S2 may suggest that these tumors arose through different carcinogenic mechanisms compared with S3.

Clearly, much remains to be learned about the biological basis of our observed HCC subclasses. But the fact that they are observed in all studies of HCC examined to date suggests that they represent a reproducible classification framework for the disease. We therefore propose that it will be important to know the subclass of HCC patients entering clinical trials for the treatment of HCC because the response to targeted agents (e.g., β-catenin and PI3K) is likely to be different across the subsets (38). Early observations of differential sensitivity of these distinct tumor types may help guide the design of future clinical trials aimed at targeting agents to distinct patient populations.

No potential conflicts of interest were disclosed.

Grant support: NIH/National Cancer Institute grant 5U54 CA112962-03 (T.R. Golub), NIH/National Institutes of Diabetes, Digestive and Kidney Diseases grant 1R01DK076986 (J.M. Llovet), NIH (Spain) grant I+D Program SAF-2007-61898 (J.M. Llovet), and Samuel Waxman Cancer Research Foundation. Y. Hoshida is supported by Charles A. King Trust fellowship. S.M.B. Nijman was supported by the Netherlands Organisation for Scientific Research (Rubicon) and Dutch Cancer Society fellowships.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

We thank Menno Creyghton for reagents and helpful suggestion; So Young Kim, Ron Firestein, William Hahn, and David Root for the shRNA constructs; Megan Hanna for technical help; Weijia Zhang for critical reading of the manuscript; and Jadwiga Grabarek and Mariko Kobayashi for general support.