Whole-Genome and Epigenomic Landscapes of Etiologically Distinct Subtypes of Cholangiocarcinoma

Cholangiocarcinoma (CCA) is the second most common hepatobiliary malignancy, accounting for 10% to 20% of primary liver cancers (1). The highest rates of CCA are in South East Asia (Northeast Thailand, Cambodia, and Laos), where >8,000 cases are diagnosed annually due to infection by liver flukes such as Opisthorchis viverrini and Clonorchis sinensis (2). CCA is considered relatively rare in Western countries (<6 per 100,000 population) where risk factors such as primary sclerosing cholangitis, hepatolithiasis, and choledochal cysts predominate (3). However, the incidence of intrahepatic CCA in the United States appears to be increasing (1).

Current 5-year survival rates for CCA after surgery and chemotherapy remain poor (<20%, ref. 4), and clinical trials evaluating targeted therapies in unselected CCA populations have shown minimal benefits (5). Existing CCA classification systems are primarily based on either anatomic location (intrahepatic, perihilar, and distal) or pathologic features (cirrhosis, viral hepatitis, and primary sclerosing cholangitis), which do not provide insights into mechanisms of CCA tumorigenesis, nor potential targets for therapy. Although previous exome-sequencing studies by our group and others have revealed a complex CCA mutational landscape (6–8), no study to date has compared fluke-positive and fluke-negative CCAs at the whole-genome level, nor fully explored the extent and contribution of structural variants and noncoding regulatory mutations to CCA pathogenesis. Little is also known about epigenetic differences between fluke-positive and fluke-negative CCAs.

Here, on behalf of the International Cancer Genome Consortium, we report an integrated genomic, epigenomic, and transcriptomic analysis of CCA involving nearly 500 CCAs from 10 countries. Comprehensive integrative clustering revealed 4 CCA clusters likely driven by distinct etiologies, with separate genetic, epigenetic, and clinical features. Our analysis uncovered new driver genes (RASA1, STK11, MAP2K4, and SF3B1) and structural variants [FGFR2 3′ untranslated region (UTR) deletion]. Within the noncoding genome, we observed a significant enrichment of promoter mutations in genes regulated by epigenetic modulation, uncovered through a novel analysis framework incorporating experimentally derived protein–DNA binding affinities and pathway information. Strikingly, we found that 2 of the CCA clusters displayed distinct patterns of DNA hypermethylation enriched at different genomic regions (CpG islands vs. shores), demonstrating for the first time the existence of distinct DNA methylation subgroups of CCA. We propose that tumors from these two subtypes may have arisen through distinct mechanisms of carcinogenesis, driven by either extrinsic carcinogenic agents or intrinsic genetic insults.

Whole-genome sequencing (WGS) was performed on 71 CCA tumors and nonmalignant matched tissues, including both fluke-associated (Fluke-Pos, 22 O. viverrini, and 1 C. cinensis samples) and non–fluke-associated cases (Fluke-Neg, 48 samples). Of these, 27 samples have been previously analyzed at the exome level (Supplementary Table S1A and S1B; refs. 6–8). Sequencing was performed to an average depth of 64.2× (median 65.1×; Supplementary Table S1C). We called somatic mutations [single-nucleotide variants (sSNV) and short insertion-deletions (indel)] using both the Genome Analysis Toolkit and MuTect (Supplementary Methods). Orthogonal validation resequencing on 97 randomly selected sSNVs and 85 indels using either Ion Torrent or Sanger technologies determined accuracy rates to be 99% for sSNVs and 87% for indels. In total, we detected 1,309,932 somatic mutations across the 71 tumors, with 4,541 nonsilent sSNVs and 1,251 nonsilent indels in protein-coding genes. On average, each CCA had 82 nonsilent somatic mutations per tumor (median 47), consisting of 64 nonsilent somatic sSNVs (median 41) and 18 indels (median 6). These mutation counts are comparable with those previously observed in genomes from pancreatic cancer (74 sSNVs and 5 indels per tumor; ref. 9), liver cancer with biliary phenotype (79 sSNVs and 24 indels per tumor; ref. 10), and hepatocellular carcinoma (70 sSNVs and 6 indels per tumor; ref. 11). Three CCAs exhibited exceptionally high mutation levels (average 5.91 sSNVs/Mb and 24.17 indels/Mb, vs. 1.39 and 3.72 for other CCAs)—these tumors exhibited mutational signatures of microsatellite instability (MSI), and 2 of these cases exhibited PolE mutations. Excluding these 3 hypermutated cases, Fluke-Pos CCAs exhibited significantly more somatic mutations compared with Fluke-Neg CCAs (median 4,700 vs. 3,143 per tumor, P < 0.05, Wilcoxon rank-sum test).

Previous studies by our group have suggested that genomic alterations in CCA may differ according to causative etiology (6–8). However, definitive exploration of these differences has been missing, due to limitations in sample sizes and genes analyzed, reliance on a single genomic platform (exome), and lack of whole-genome information. To address these limitations, we assembled a cohort of 489 CCAs (Supplementary Table S1A and S1B), including 133 Fluke-Pos (132 O. viverrini and 1 C. cinensis samples) and 356 Fluke-Neg cases. Samples were analyzed using four different genomic platforms based on sample availability. Besides WGS (71 cases), these included exome sequencing of 200 cases (previously published; ref. 8), high-depth targeted sequencing of 188 cases, SNP array copy-number profiling of 175 cases, array-based DNA methylation profiling of 138 cases, and array-based expression profiling of 118 cases. To confirm the applicability of merging different data types, we used statistical models to confirm that our mutation calls were not biased by differences in sequencing platforms, and directly confirmed mutation concordance by analyzing those samples sequenced on overlapping platforms (Supplementary Methods). Using iClusterPlus, we performed integrative clustering combining data from somatic mutations, somatic copy-number alterations (sCNA), mRNA expression, and DNA methylation on 94 CCAs where all four data types were available. Randomized subsampling clustering confirmed the robustness of these integrative clusters (Supplementary Fig. S1A). To support the reliability of our conclusions, reanalysis using an expanded set of integrative clustered samples (121 samples), including samples with one or more missing platforms while retaining a cluster prediction accuracy of 90%, yielded similar results and associations (Supplementary Fig. S1B; Supplementary Methods). We also found that although the most discriminatory clustering was achieved by using all available genomic information, clustering by individual modalities recapitulated some of the integrative clusters, for sCNA, expression, and methylation data (Supplementary Fig. S1C; see later section on DNA methylation). Importantly, integrative clustering performed on samples stratified by anatomic location reproduced the original clusters within each anatomic site (≥88% concordance; Supplementary Fig. S1D), demonstrating that the molecular clusters are not simply recapitulating anatomic variation.

Integrative clustering revealed 4 distinct clusters characterized by different clinical features and genomic alterations (Fig. 1A). Cluster 1 comprised mostly Fluke-Pos tumors, with hypermethylation of promoter CpG islands (aberrantly methylated above normal level; Supplementary Methods), enrichment of ARID1A and BRCA1/2 mutations (P < 0.01 and P < 0.05, respectively, Fisher exact test), high levels of nonsynonymous mutations (P < 0.001, Wilcoxon rank-sum test; Supplementary Fig. S2A), and high levels of mutations in gene promoters with histone 3 lysine 27 trimethylation (H3K27me3) predicted to alter transcription factor binding (see later section on somatic promoter mutations). Cluster 2 was characterized by a mix of Fluke-Pos and Fluke-Neg tumors, with upregulated CTNNB1, WNT5B, and AKT1 expression (P < 0.05, Wilcoxon rank-sum test, Supplementary Fig. S2B) and downregulation of genes involving EIF translation initiation factors (Supplementary Table S2A). Most notably, clusters 1 and 2 were also significantly enriched in TP53 mutations and ERBB2 amplifications (P < 0.001 and P < 0.01, respectively, Fisher exact test) and elevated ERBB2 gene expression (P < 0.05, Wilcoxon rank-sum test; Fig. 1A and Supplementary Fig. S2C).

In contrast to clusters 1 and 2, clusters 3 and 4 comprised mostly Fluke-Neg tumors. Cluster 3 displayed the highest level of sCNAs, including enrichment of amplifications at chromosome arms 2p and 2q (q < 0.05, Fisher exact test; Supplementary Table S2B). Analysis of immune populations by ESTIMATE revealed that both clusters 2 and 3 displayed immune cell infiltration (Supplementary Fig. S2D), but only cluster 3 exhibited specific upregulation of immune checkpoint genes (PD-1, PD-L2, and BTLA; Fig. 1B and Supplementary Fig. S2E) and pathways related to antigen cross-presentation, CD28 co-stimulation, and T-cell signal transduction (Supplementary Table S2A). Cluster 4 was characterized by BAP1, IDH1/2 mutations, FGFR alterations (all P < 0.01, Fisher exact test), and upregulated FGFR family and PI3K pathway signatures (Fig. 1A and Supplementary Table S2A). Similar to cluster 1, cluster 4 tumors also exhibited DNA hypermethylation—however, rather than hypermethylation at CpG islands, cluster 4 hypermethylation was at CpG promoter shores (see later section on DNA methylation).

We sought to relate the CCA clusters to anatomic and clinical features. Clusters 1 and 2 were enriched in extrahepatic (consisting of perihilar and distal) tumors, whereas clusters 3 and 4 were composed almost entirely of intrahepatic tumors (P < 0.001, Fisher exact test). This was observed in both Fluke-Neg and Fluke-Pos tumors, and persisted after adjusting for fluke status (P < 0.001, multivariate regression; also confirmed in expanded clusters). Other CCA risk factors, such as hepatitis B virus (HBV), hepatitis C virus (HCV), and primary sclerosing cholangitis (PSC), were present in our cohort at frequencies of 10.4%, 2.9%, and 1.0%, respectively. Both HBV and PSC were associated with intrahepatic CCA (P < 0.05, Fisher exact test; ref. 12).

We further investigated the prevalence of driver genes according to anatomic location (among 459 samples with sequencing and anatomic location information). BAP1 and KRAS were more frequently mutated in intrahepatic cases (q < 0.1, Fisher exact test). This was observed in both Fluke-Neg and Fluke-Pos tumors and persisted even after adjusting for fluke status (q < 0.1, multivariate regression). No additional genes were identified as differentially mutated when the extrahepatic CCAs were analyzed as perihilar and distal CCAs.

Clinically, patients in clusters 3 and 4 had significantly better overall survival relative to the other 2 clusters (P < 0.001, log-rank test; Fig. 1C). As fluke infection was also associated with poorer survival (P < 0.001, log-rank test; Supplementary Fig. S2F), we performed multivariate analysis and confirmed that this cluster-associated survival difference persisted even after accounting for fluke association, anatomic location, and clinical staging (P < 0.05, Cox proportional hazards model; Supplementary Table S2C). To validate this finding, we assembled a separate validation cohort comprising newly classified samples from the expanded integrative clustering and a recently published set of CCA samples (Supplementary Methods; ref. 13). A survival analysis on this independent validation cohort reaffirmed the same survival trends, on both univariate (P < 0.05, log-rank test; Supplementary Fig. S2F) and multivariate analysis (P < 0.05, Cox proportional hazards model; Supplementary Table S2C). This result demonstrates that molecular clusters can provide additional prognostic information in a manner independent of fluke status and anatomic location. Figure 1D summarizes the salient features of each CCA cluster.

Driver gene mutation analysis across 459 CCAs (130 Fluke-Pos and 329 Fluke-Neg cases) revealed 32 significantly mutated genes (SMG; q < 0.1 by both MutSigCV and IntOGen; Supplementary Table S3A–S3D and Fig. 2A). Our analysis revealed 4 potentially new CCA driver genes not highlighted in previous CCA publications (6–8, 13–17): RASA1, STK11, MAP2K4, and SF3B1. Of these, RASA1, STK11, and MAP2K4 are related to RAS/MAPK signaling (Fig. 2B and C and Supplementary Table S3B). RASA1, encoding a p120 RAS GTPase-activating protein, was predicted to be inactivated in 4.1% of cases (10 frameshift, 4 nonsense; Fig. 2D). These inactivating mutations, along with observed focal RASA1 copy-number losses, were associated with decreased RASA1 expression (Fig. 2E). In CCA cell lines, shRNA-mediated knockdown of RASA1 resulted in significantly enhanced migration and invasion, supporting a tumor-suppressor role for RASA1 in CCA (Fig. 2F). STK11, a serine/threonine protein kinase, was mutated in 5% of cases, with most STK11 mutations also predicted to be inactivating (7 nonsense, 9 frameshift; Fig. 2G). SF3B1, an RNA splicing factor, was mutated in 4.6% of cases, at mutation hotspots (23% at codon 625 and 14% at codon 700) previously observed in uveal melanoma and breast cancer (Fig. 2H; refs. 18, 19). This latter finding may implicate a role for RNA splicing dysregulation in CCA tumorigenesis. MAP2K4 is a member of the mitogen-activated protein kinase and has been shown to activate p38 mitogen-activated protein kinase and JUN N-terminal kinase. We observed MAP2K4 focal homozygous deletions in 2 Fluke-Pos cases (Supplementary Fig. S3A–S3C), and MAP2K4 mutations in another 10 cases (2.2%). Half of these were predicted to be inactivating (frameshift, nonsense, splice-site mutations), consistent with a tumor suppressor role for MAP2K4 in CCA.

ERBB2 was amplified in 3.9% to 8.5% of CCAs (Supplementary Table S3E). ERBB2 amplifications were more frequent in Fluke-Pos cases (10.4% in Fluke-Pos vs. 2.7% in Fluke-Neg CCA, P < 0.01, Fisher exact test) with elevated ERBB2 gene expression in Fluke-Pos compared with Fluke-Neg tumors in these cases (Supplementary Fig. S2C; Supplementary Table S3E for validation samples). On average, ERBB2-amplified samples exhibited 14 ERBB2 copies [copy numbers determined by ASCAT (SNParray) or Quandico (sequencing data)], and gene set enrichment analysis (GSEA) confirmed upregulation of ERBB2-related gene sets among ERBB2-amplified samples (q < 0.2; Supplementary Table S3F). We independently validated the presence of ERBB2 amplification in selected cases by FISH (Supplementary Fig. S3D). Other pathways upregulated in ERBB2-amplified samples included biological oxidation, metabolism cytochrome P450, peroxisome proliferator-activated receptor signaling, and checkpoint signaling (q < 0.2). In addition to ERBB2 amplifications, we also detected activating ERBB2 mutations (S310F/Y, G292R, T862A, D769H, L869R, V842I, and G660D) in 9 cases (2%). Notably, previous studies in cell lines have shown that high ERBB2-expressing CCAs may be more sensitive to ERBB2-inhibitor treatment compared with low ERBB2-expressing cases, suggesting that tumors with high ERBB2 expression may be candidates for anti-ERBB2/HER2 therapy (20). Amplification of other selected oncogenes included MYC (n = 12), MDM2 (n = 9), EGFR (n = 11), and CCND1 (n = 7), whereas deletions included CDKN2A (n = 17), UTY (n = 17), and KDM5D (n = 16; Supplementary Table S3G).

Availability of WGS data also allowed us to investigate the role of structural variations (SV) in CCA. Using CREST, we identified approximately 93 somatic SVs per tumor (median 69; range 0–395), with a 91% (61/67) true positive rate by PCR. Most of the SVs were intrachromosomal (65%) and associated with cancer-related genes (ARID1A and CDKN2A/B), retrotransposon-associated genes (TTC28), and fragile sites (1q21.3; Supplementary Table S4). SV burden varied significantly across the 4 CCA clusters (P < 0.05, Kruskal–Wallis test), with Fluke-Neg tumors in cluster 4 associated with low burden (P < 0.05, 1-sided Wilcoxon rank-sum test). TP53, FBXW7, and SMAD4 were significantly associated with increased SV burden (q < 0.1, Wilcoxon rank-sum test; Supplementary Fig. S4A).

FGFR2 fusion genes, previously reported in CCA (21), are thought to deregulate FGFR2 signaling through the hijacking of 3′ fusion partners with dimerization motifs (21, 22). However, whether rearrangements affecting other FGFR members (besides FGFR2) exist in CCA, or if other categories of FGFR2 rearrangements (besides in-frame gene fusions) can contribute to CCA development, remains unclear. Analysis of the WGS data, followed by subsequent validation at the gene transcript level, revealed 5 in-frame gene fusions with intact tyrosine kinase domains—4 involving FGFR2 (FGFR2–STK26, FGFR2–TBC1D1, FGFR2–WAC, and FGFR2–BICC1; Fig. 3A) and 1 involving FGFR3 (FGFR3–TACC3; Fig. 3B). To our knowledge, this is the first report of FGFR3 fusions in CCA. FGFR3–TACC3 fusions have been previously reported in bladder cancer, glioblastoma, and lung cancer (23, 24) and shown to be oncogenic.

Besides FGFR in-frame fusions, we also identified recurrent truncating events translocating FGFR2, without its 3′UTR, to intergenic regions (Fig. 3C). FGFR2-truncated CCAs exhibited high expression levels compared with FGFR2 transcripts with intact 3′UTRs (P < 0.01, Wilcoxon rank-sum test). In vitro, luciferase reporter experiments confirmed diminished expression in constructs containing FGFR2 3′UTRs compared with control reporters (Fig. 3D). FGFR2 3′UTR loss may thus represent a new and additional mechanism for enhancing FGFR2 expression in CCA, similar to mechanisms reported for PD-L1 (25).

FGFR2 rearrangements were observed exclusively in cluster 4 (P < 0.001, Fisher exact test; Fig. 1A). In CCAs lacking FGFR rearrangements, we further identified other genomic alterations involving FGFR2, including indels (n = 3), SNVs (n = 10), and copy gains (n = 1)—several of these alterations have been previously shown to be activating (26, 27). Collectively, CCAs with altered FGFR2 genes (point mutations, indels, copy gain, and rearrangements) were significantly enriched in cluster 4 (P < 0.01, Fisher exact test) and also expressed significantly higher FGFR2 levels than FGFR2–wild-type tumors (Figs. 1A, 2B, and 3E).

In addition to FGFR fusions, WGS analysis also identified rearrangements affecting the catalytic subunit B of cAMP-dependent protein kinase A (PRKACB), including ATP1B1–PRKACB and LINC00261–PRKACB (Fig. 3F). Both PRKACB rearrangements retained the PRKACB pseudokinase domain, and thus may increase PKA activity and activate downstream MAPK signaling (8).

Long interspersed nuclear element-1 (LINE-1 or L1) repeats are autonomous retrotransposons collectively occupying approximately 17% of the human genome. Recent studies have shown that certain L1 elements are active in cancer, displaying somatic retrotransposition and potentially contributing to genomic instability and tumorigenesis (28–30). In the CCA WGS samples, we observed frequent somatic L1 retrotranspositions, particularly originating from an L1 element in intron 1 of the TTC28 gene (52 events in 20/71 tumors, 28.2%; Supplementary Fig. S4B). PCR testing validated 98% (46/47) of these retrotransposition events (Supplementary Table S4). CCAs with L1 retrotransposition were associated with Fluke-Pos tumors (P < 0.01, Fisher exact test) and increased SV burden (P < 0.05, Wilcoxon rank-sum test), suggesting a relationship between genomic instability and L1 endonuclease activity. Further analysis revealed that these intragenic insertions were not overtly associated with cancer-related genes, such as tumor suppressors.

Somatic mutations in noncoding regulatory regions have been proposed to play crucial roles in carcinogenesis (31). However, systematic approaches for identifying such mutations are lacking and their effects in CCA remain poorly understood. To date, only TERT promoter mutations have been observed in CCA (8), and in our WGS cohort, only 2 CCAs (2.8%) harbored TERT-promoter mutations (chr5:1295228). Besides TERT, no other recurrent noncoding promoter-region mutations were observed in the WGS cohort.

Even when integrating mutations over regulatory regions or gene promoters, the low recurrence rate of most noncoding mutations may lead to a lack of statistical power to identify potential drivers in the promoters of individual genes. We hypothesized that the effects of noncoding promoter mutations might instead be detectable at the gene-set level. To test this hypothesis, we developed a novel method, FIREFLY (FInding Regulatory mutations in gEne sets with FunctionaL dYsregulation; Fig. 4A), which identifies gene sets dysregulated by somatic promoter mutations that alter transcription factor (TF) binding. Compared with approaches used in previous cancer genome studies (32, 33), FIREFLY differs in three important respects. First, it uses experimentally determined high-throughput TF-DNA binding data for 486 TFs representing a broad range of TF families (34, 35), as opposed to position weight matrices (PWM; ref. 36), to predict mutation-associated changes in TF binding affinity. Second, FIREFLY condenses the large numbers of highly nonrecurrent noncoding mutations into biologically meaningful gene sets, shortlisting those sets with an overrepresentation of mutations, as assessed by multiple statistical tests. Third, it orthogonally validates the transcriptional consequences of the binding-change predictions using expression data from primary tumors.

To identify sets of genes that had dysregulated transcription in the aggregate due to promoter mutations, we applied FIREFLY to 70 WGS samples (1 hypermutated sample was excluded), representing 6,639 somatically mutated gene promoters. Based on binding-change predictions, FIREFLY identified 138 sets of genes that were enriched for binding-change mutations in promoters (q < 0.1, Fisher exact test with Benjamini–Hochberg false discovery rate). FIREFLY's second statistical test then compared the number of binding-change mutations in these sets with an expected null distribution determined from synthetically mutated data, created using the tumor-specific mutation rates for each type of mutation in its trinucleotide context. Nineteen sets passed the second test (q < 0.1). Finally, FIREFLY orthogonally assessed the transcriptional impact of the binding-change predictions, by testing whether tumors with increasing numbers of binding-change mutations for a given gene set also exhibit greater transcriptional dysregulation in that set (q < 0.1, by Gene Set Analysis; ref. 37). Four of the 19 sets were validated by this test (Fig. 4B and C and Supplementary Fig. S5A for example nonsignificant gene sets). This was significantly greater than expected for a null distribution based on randomly selected gene sets of similar sizes (P < 0.01; Supplementary Fig. S5B). We note that FIREFLY results do not imply that every somatic promoter mutation in the 4 gene sets affects gene expression. Instead, the results indicate that the identified sets have excesses of promoter mutations that likely affect gene expression. We also note that FIREFLY does not use predicted gain or loss of binding to infer directionality of change in expression levels. This is because, in general, TFs can act as either activators or repressors depending on the regulatory context. Consequently, one cannot in general predict whether a gain or loss of a binding site will result in upregulation or downregulation.

To validate a fundamental assumption in the FIREFLY pipeline, that is, that mutations predicted to alter TF binding also affect transcription, we selected 3 mutations and tested them in luciferase reporter assays. We confirmed altered regulatory activity for 2 of the 3 mutations, providing direct evidence that mutations predicted to change TF binding can in fact alter gene expression (Supplementary Fig. S5C).

FIREFLY identified 4 gene sets that were likely dysregulated by altered TF-DNA binding. Interestingly, 2 of these (MIKKELSEN_MCV6_HCP_WITH_H3K27ME3 and MIKKELSEN_MEF_ICP_WITH_H3K27ME3) are subsets of PRC2 target genes that have promoters with histone modification H3K27me3 in certain contexts (Fig. 4B). Mutations in all 4 enriched gene sets occurred across all 4 of the CCA clusters. It was noteworthy, however, that cluster 1 was enriched for mutations in 3 of the 4 gene sets, including 2 associated with H3K27me3 (Supplementary Fig. S5D). This, together with the observed hypermethylation of polycomb repressive complex 2 (PRC2) target genes in cluster 1 (see later section on DNA methylation), provides additional evidence of the importance of alterations in PRC2 regulation in cluster 1.

As shown in Fig. 1, iCluster analysis revealed 2 distinct hypermethylated CCA subgroups (clusters 1 and 4). To validate these differences across a larger CCA series, we performed unsupervised DNA-methylation clustering on an expanded panel of 138 CCAs. Clustering based solely on DNA methylation recapitulated 2 hypermethylated clusters highly concordant with clusters 1 and 4 (96.3% and 86.1% concordance, respectively; Fig. 5A), and a third group of low-methylation tumors representing a mix of clusters 2 and 3. Cluster 1, enriched in Fluke-Pos CCAs, was dominated by hypermethylation in promoter CpG islands, whereas cluster 4, enriched in Fluke-Neg CCAs, was dominated by hypermethylation in promoter CpG island shores (Fig. 5B). Different sets of gene promoters were targeted for hypermethylation in clusters 1 and 4. However, GSEA revealed that both affected common pathways, including PRC2 targets. We observed significant inverse correlations between transcript levels and promoter methylation in both clusters 1 and 4 (q < 0.05; Supplementary Fig. S6A), consistent with these epigenetic alterations exerting transcriptional impact.

DNA hypermethylation in the 2 hypermethylated clusters may be driven by distinct epigenetic mechanisms. In cluster 1, we observed downregulation of the DNA demethylation enzyme TET1 and upregulation of the histone methyltransferase EZH2 (Supplementary Fig. S6B), suggesting a possible role for these genes in establishing the hypermethylation phenotype (38). In contrast, cluster 4 CCAs were significantly enriched in IDH1/2 mutations, which are known to be associated with CCA hypermethylation (31.6% in cluster 4 vs. 1.0% in other clusters, q < 0.001, multivariate regression; Fig. 5A; refs. 7, 13, 39). Among cluster 4 CCAs lacking IDH1/2 mutations (68.4%), BAP1 mutations were enriched (q < 0.001 and 0.05, respectively for inactivating point mutations and regional deletions; Fig. 5A). BAP1 mutated cases were also associated with increased CpG hypermethylation relative to BAP1–wild-type cases (Supplementary Fig. S6C). Notably, BAP1 mutations have been associated with DNA hypermethylation in CCA and renal cell carcinoma (13, 40).

To explore mutation patterns between these 2 clusters, we identified 10 established mutation signatures in the WGS cohort. These included Catalogue of Somatic Mutations in Cancer (COSMIC) Signatures 1, 5, 8, 16, and 17, and signatures associated with activated APOBECs (signatures 2 and 13), mismatch-repair deficiency (MMR; Signatures 6 and 20), and aristolochic acid exposure (signature 22, Supplementary Fig. S6D). Signature 5 burdens were correlated with patient age (Spearman correlation 0.25, P < 0.05), as previously reported for other cancer types (41). Fluke-Pos CCAs were enriched for activated APOBEC mutation burden (P < 0.001, multivariate regression). Signatures of MMR and signatures 8, 16, and 17 have not been previously reported in CCA.

We observed elevated levels of signature 1 (CpG>TpG mutations) in cluster 1, even after adjusting for patient age (P < 0.001, multivariate regression; Fig. 5C). Importantly, this elevation is signature 1–specific, as it was not observed for signature 5. We note that CpG dinucleotides are known mutation hotspots, due to spontaneous deamination of 5-methylcytosine (5mC) to thymine (CpG>TpG mutation; ref. 42). To investigate if hypermethylated CpGs in cluster 1 might provide susceptible genomic substrates for deamination and subsequent signature 1 mutations, we integrated the locations of the CpG>TpG mutations with regions of hypermethylation. In cluster 1, CpG>TpG mutations were indeed located preferentially near hypermethylated regions (P < 0.001, Fisher exact test; Fig. 5D and Supplementary Fig. S6E), whereas in cluster 4, these mutations showed no such regional preferences. These results support a significantly increased level of DNA hypermethylation–related deamination events in cluster 1 compared with cluster 4.

We further investigated if the differences in genome-wide mutation patterns between clusters are accompanied by distinct clonal structures harboring these mutations. Distribution analysis of variant allele frequencies (VAF) for point mutations (in copy-neutral regions and adjusted for tumor purity) revealed a wide spread of VAFs in cluster 1 compared with cluster 4 (Fig. 5E), indicating the presence of heterogeneous subclones in cluster 1 tumors, but not in cluster 4 tumors. Together, these distinct patterns of hypermethylation-related deamination and tumor heterogeneity suggest disparate somatic-evolution processes during tumorigenesis in these clusters (see Discussion).

Surgery is the only proven treatment modality for CCA (43), and all formal evaluations of targeted therapies to date, performed in unselected CCA cohorts, have proved unsuccessful (44). In this study, we analyzed a cohort of nearly 500 CCAs from distinct geographic regions, including 94 CCAs covered by four genomic platforms. Integrative clustering of mutation, copy number, gene expression, and epigenetic data revealed 4 subtypes of CCA, each exhibiting distinct molecular and clinicopathologic features. Four lines of evidence highlight that clustering based on molecular profiles provides additional information beyond anatomic site. First, anatomic site does not drive molecular subtypes, as evidenced by the reproducibility of the molecular subtypes within each anatomic site separately (Supplementary Fig. S1D). Second, tumors in different anatomic sites may exhibit similarities at the molecular level, whereas tumors located in the same anatomic site can display profound differences in their molecular profiles. This is exemplified by clusters 1 and 2 comprising mixtures of intrahepatic and extrahepatic tumors, while intrahepatic tumors are split among all 4 clusters. Third, from a disease prognosis standpoint, CCAs in different anatomic sites do not differ in their survival trends, whereas the molecular clusters show significant differences in survival, in both the original and an independent validation cohort (Fig. 1C and Supplementary Fig. S2F). Finally, we note that current medical oncology guidelines do not discriminate CCA treatments based on their anatomic site (1) as classifications based on anatomy do not provide information regarding potential therapeutic opportunities. In contrast, the molecular profiles highlight several potential cluster-specific therapies (see next paragraph). Taken collectively, our results indicate that molecular CCA subtypes based on integrative molecular clustering are likely to offer enhanced information regarding CCA biology and clinical behavior, beyond that provided by anatomic location alone.

Examination of signature genomic alterations in each subtype highlights potential therapeutic opportunities, although we emphasize that such findings require further clinical validation. For example, it is possible that CCAs in clusters 1 and 2 with ERBB2 amplification may be appropriate for therapies targeting ERBB2/HER2 signaling (45). The elevated expression of immune related genes and pathways in cluster 3 CCAs suggests a therapeutic opportunity for immunotherapy; however, this conclusion should be treated with caution due to the small sample size of this cluster. The specific mechanisms driving the elevated expression of immune-related genes in cluster 3 remain unclear; however, increased levels of immunogenicity and upregulation of MHC protein and antigen processing complexes have been independently observed in aneuploid tumors (46, 47), consistent with the elevated level of chromosomal aberrations in cluster 3. Cluster 4 CCAs, which are associated with IDH1/2 mutations and FGFR2 and PRKA-related gene rearrangements, might also be tested with recently described IDH inhibitors (ClinicalTrials.gov identifier: NCT02073994) or FGFR-targeting agents. Notably, our study suggests that loss of FGFR2 3′UTRs may represent an alternative mechanism of FGFR2 activation, beyond FGFR2 in-frame gene fusions and activating mutations. This, along with our additional observation of FGFR3 rearrangements, expands the proportion of potential FGFR-targetable cases in CCA (48).

Previous cancer studies have discovered recurrent regulatory mutations in individual promoters such as TERT (33, 49). In complement to these studies, we developed in this work an alternative analysis framework (FIREFLY), which uses protein binding microarray (PBM)–based analysis to examine the effects of noncoding promoter mutations on gene-expression pathways in cancer. FIREFLY's use of PBM data with a k-mer—based approach provides advantages for estimating the effects of mutations on TF-binding affinities compared with traditional PWM models. For example, PWMs implicitly assume that each base pair contributes independently to binding affinity; however, this is not always true. Moreover, PWMs do not capture the multiple modes of binding associated with many TFs, and PWM scores are not easily comparable between TFs. FIREFLY addresses these shortcomings by using experimentally determined PBM data, which provides actual binding affinities of every possible 8-mer. Using FIREFLY, we identified several pathways with strong statistical evidence for recurrent systematic dysregulation in CCA. Interestingly, 2 of these pathways reflected epigenetically modulated cellular differentiation processes, which have been shown to be dysregulated in cancer in general by other means such as DNA hypermethylation or histone modification.

Of the 4 clusters, 2 clusters (clusters 1 and 4) are most clearly distinguished by their highly distinctive patterns of genome-wide DNA hypermethylation, targeting either promoter CpG islands or promoter CpG shores. To our knowledge, such epigenetically distinct tumor subtypes have not been previously reported in the literature, particularly for tumors from the same tissue type. Further analysis demonstrated that cluster 1 CCAs are Fluke-Pos with elevated mutation rates, Mutation Signature 1 enrichment, and increased point-mutation subclonality, whereas cluster 4 CCAs are Fluke-Neg and by comparison relatively clonal. We propose that these differences are consistent with a model where early in tumorigenesis, cluster 1 CCAs are likely driven by external carcinogenic agents and early epigenetic deregulation (“epimutations”), whereas cluster 4 CCAs are likely driven by pioneer genetic events such as IDH1/2 or BAP1 mutations, with epigenetic aberrations arising as a downstream consequence (Fig. 6). In this model, fluke infection, by inducing chronic inflammation, metabolic disruption of host bile homeostasis (50), or secretion of fluke growth factors and excretory vesicles for modulating host–pathogen interactions (51), induces genome-wide epigenetic deregulation. Cytosine residues experiencing aberrant methylation are then at higher risk of spontaneous deamination and mutation, consistent with the enrichment of signature 1 mutations in this subtype. CpG island hypermethylation in this subtype may also silence tumor suppressor genes, further enhancing cancer development. Moreover, because the processes of carcinogen-induced methylation, deamination, and mutation are inherently stochastic from cell to cell, such events would inevitably lead to increased levels of intratumor heterogeneity. In contrast, in cluster 4 CCAs which are Fluke-Neg, somatic mutations in critical chromatin modifier genes (e.g., IDH1/2 and BAP1) may occur as a primary event preceding epigenetic deregulation, driving both rapid clonal outgrowth and directly inducing DNA hypermethylation. Specifically, IDH1/2 mutations have been shown to increase 2-hydroxyglutarate oncometabolite production, leading to DNA hypermethylation.

We acknowledge that other models may also explain the striking molecular differences between clusters 1 and 4. These include differential vulnerabilities in distinct cells of origin, as the biliary system is known to contain multipotent stem/progenitor cells. Liver fluke infection primarily affects large intrahepatic and extrahepatic bile ducts, giving rise to intrahepatic and/or extrahepatic CCA. Conversely, parenchymal liver diseases exclusively affect canals of Hering and bile ductules, and are primarily associated with intrahepatic CCA (52).

A recent study from The Cancer Genome Atlas (TCGA) consortium reported a multi-omic analysis of a smaller and more homogenous CCA series (38 samples, exclusively Fluke-Neg, mostly intrahepatic and North American; ref. 13). Comparison of this study with our own data revealed an obvious difference—specifically our inclusion of Fluke-Pos samples allowed the discovery of another major CCA subtype (cluster 1) with distinct clinical and molecular features. Besides cluster 1, the other study's “IDH” (IDH mutants) and “METH3” (BAP1 mutants and FGFR rearrangements) groups are likely matches to our cluster 4 (characterized by IDH and BAP1 mutants, and FGFR rearrangements), whereas their “ECC” group (extrahepatic) matched cluster 2 (containing Fluke-Neg extrahepatic tumors). On the other hand, the TCGA “METH2” group (CCND1 amplifications) and our cluster 3 were not obviously matched. Taken collectively, these results suggest that most of the TCGA study's clusters are largely concordant with our own, and neither classification strictly precludes the other.

We acknowledge that limitations of sample resources (DNA, RNA, and paraffin-embedded tissues) were a major constraint in this study. We were unable to generate data using all platforms on all samples, which reduced the sample size in the integrative clustering analysis. To overcome this constraint, we sequenced an extended and separate sample cohort, to validate findings emerging from the integrative clustering. Center-specific differences in presample processing steps, including collection site, biopsy site, and sample processing protocols, may also result in sequencing biases. We attempted to control for these variations by reviewing the histology of all cases using standardized American Joint Committee on Cancer (AJCC) 7th criteria to confirm and harmonize the histology and anatomic subtype of samples. We also attempted to reduce biopsy bias and normal contamination originating from the biopsy sites by estimating tumor cell content through histopathologic review or SNP arrays. Lastly, our study merged data from different DNA sequencing platforms [WGS, whole-exome sequencing (WES), and targeted sequencing], thus limiting our analysis across the entire cohort to genomic regions common to these platforms. However, we overcome biases from different data processing centers and merging different data types by using one analysis pipeline, ensuring uniform analysis on each data type, and confirming good concordance between samples sequenced on multiple platforms (Supplementary Methods).

In summary, integrative analysis of a large CCA cohort has revealed a novel molecular taxonomy, with discovery of new potential CCA driver genes and gene rearrangements. This taxonomy may be clinically relevant; however, more functional data are required to test the validity of the genomic data. Analysis of noncoding promoter mutations, made possible by whole-genome analysis, revealed that they play a significant role in CCA targeting genes involved in PRC2 biology, and we identified 2 highly distinct CCA subtypes demonstrating distinct DNA hypermethylation patterns. We conclude by noting that these last two findings may carry conceptual relevance beyond biliary tract cancers. Specifically, although elevated DNA methylation levels have been observed in numerous cancers, these cases to date been largely consigned to a general CpG island methylator phenotype (CIMP) class—our data suggest that a more detailed examination of these cases may reveal potential for epigenetic heterogeneity. Finally, as exemplified by the TERT promoter, previous analyses of noncoding promoter mutations have largely focused on identifying individual regions of recurrent genomic aberration. Our FIREFLY analysis suggests that such cases are likely rare and that, similar to protein coding genes, mutations in noncoding regulatory regions may target genes in specific pathways rather than individually. Extending this concept to promoter mutation catalogues of other cancer types will certainly test the general applicability of this proposition.

Primary tumor and matched normal samples (nonneoplastic liver or whole blood) were obtained from the SingHealth Tissue Repository (Singapore), Fundeni Clinical Institute (Romania), Khon Kaen University (Thailand), ARC-Net Biobank (Italy), Centre de Resources Biologiques Paris-Sud (France), Department of Pathology, Yonsei University College of Medicine (South Korea), Hospital do Cancer de Barretos (Brazil), Linkou Chang Gung Memorial Hospital (Taiwan), and Department of Pancreatobiliary Surgery, The First Affiliated Hospital, Sun Yat-Sen University (China) with signed informed consent. WGS data from 10 Japanese cases were also included (National Cancer Center, Japan). The study was approved by the SingHealth Centralized Institutional Review Board (2006/449/B), Ethics Committee of the Clinical Institute of Digestive Diseases and Liver Transplantation, Fundeni (215/18.01.2010), Khon Kaen University (HE471214), Centre de Resources Biologiques Paris-Sud, National Cancer Center, Japan (G20-03), Severance Hospital, Yonsei University Health System (4-2014-0829), Hospital do Cancer de Barretos (716/2013), Linkou Chang Gung Memorial Hospital (100-2030B), The First Affiliated Hospital of Sun Yat-Sen University (2014/C_006), and ARC-Net Biobank at Verona University Hospital (n. prog. 1959). Clinicopathologic information for subjects, including age, sex, histology, tumor subtype, stage, and overall survival, were reviewed retrospectively. Cases were staged according to the AJCC Staging System 7th Edition. All patients had not received prior treatment. In total, 489 tumors with associated clinicopathologic data were obtained (133 Fluke-Pos: 132 O. viverrini, 1 C. sinensis; 39 HBV/HCV-positive; 5 PSC-positive). These were assayed on at least one profiling platform, which included: (i) WGS (71 CCAs); (ii) targeted sequencing surveying 404 genes (188 CCAs); (iii) published exome sequencing (ref. 8; 200 CCAs); (iv) HumanOmniExpress BeadChip arrays (SNP arrays; 175 CCAs); (v) DNA methylation 450k BeadChip arrays (138 CCAs); and (vi) HumanHT-12 Expression BeadChip arrays (gene expression arrays; 118 CCAs; Supplementary Tables S1 and S5). A detailed list of clinical data is included in Supplementary Table S1A.

H69 nonmalignant immortalized cholangiocyte cells were obtained in 2011 from D. Jefferson (New England Medical Center, Tufts University) and cultured as previously described (53). HEK293T cells were obtained from the ATCC (CRL-3216) in 2015 and cultured with DMEM, 10% FBS, 2 mmol/L L-glutamine, and 1% penicillin/streptomycin. EGI-1 was purchased from DSMZ (ACC 385) and maintained in DMEM supplemented with 10% FBS (Sigma-Aldrich). HUCCT1 (JCRB0425) cells were purchased from the Health Sciences Research Resources Back (HSRRB) in 2009 and maintained in RPMI-1640 medium with 10% FBS. M213 (JCRB1557) cells were obtained from the Liver Fluke and Cholangiocarcinoma Research Center in 2015 and cultured with Ham's F12 media (Gibco). Cells were cultured at 37°C in a 5% CO₂ humidified chamber. All cell lines were authenticated by short-tandem repeat profiling and found to be negative for Mycoplasma as assessed by the MycoSensor qPCR Assay Kit (Agilent Technologies).

Genomic DNA was extracted using the QIamp DNA mini kit (Qiagen). DNA yield and quality were determined using Picogreen (Invitrogen) and further visually inspected by agarose gel electrophoresis. RNA was extracted using an RNeasy mini kit (Qiagen), quantified by measuring Abs260 with a UV spectrophotometer, and quality assessed with the Agilent 2100 Bioanalyzer (Agilent Technologies). Sequencing libraries were prepared from DNA extracted from tumor and normal samples using the SureSelect XT2 Target Enrichment System for the Illumina Multiplexed Sequencing platform (Illumina) according to the manufacturer's instructions. WGS was performed using Illumina HiSeq X10, Illumina HiSeq2500, and Illumina HiSeq2000 instruments. To survey the frequency and distribution of somatic mutations in the validation cohort (188 CCAs), targeted sequencing of 404 genes was performed after capture with a custom SureSelect capture reagent designed using the SureDesign tool (Agilent Technologies). Target-enriched libraries were sequenced on the Illumina HiSeq 4000 sequencing platform. Coverage of coding regions based on the amplicon design was 99.6% (Supplementary Table S5).

For 3′ UTR reporter assays, control reporter plasmids (LUC) were generated by cloning the SV40 promoter into the pGL4.10 promoterless luciferase reporter vector (Promega). The LUC-FGFR2_3′UTR test plasmid was engineered by inserting the 1,666 bp 3′UTR region of FGFR2 (starting from stop codon) into the immediate 3′ end of the luciferase gene of the LUC plasmid. HEK293T or H69 cells were cotransfected with a Renilla-containing plasmid with either LUC or LUC-FGFR2_3′UTR. Transfected cells were incubated for 24 hours.

For promoter mutation reporter assays, luciferase constructs were generated by ligating fragments of promoter regions (2 kb upstream and 500 bp downstream of the transcription start site) prepared from genomic PCR using gene-specific primers. Mutant constructs were generated using the QuikChange II XL site-directed mutagenesis kit (Agilent Technologies) according to the manufacturer's instructions. H69 or EGI-1 cells were transfected and incubated for 48 hours, after which Dual-Luciferase Reporter Assays (Promega) were performed. Relative luciferase activity was calculated as the firefly luciferase activity normalized to Renilla (Promega) luciferase activity. Assays were conducted in triplicate and each experiment was repeated 3 independent times.

Construction of hairpin-pLKO.1 vectors (carrying a puromycin antibiotic resistance gene) containing shRNAs targeting RASA1 coding sequences were performed as follows: shRNA#1 (CCGGGCTGCAAGAACACTGATATTACTCGAGTAATATCAGTGTTCTTGCAGCTTTTTG; TRCN0000356449, Sigma) and shRNA#2 (CCGGCCTGGCGATTATTCA CTTTATCTCGAGATAAAGTGAATAATCGCCAGGTTTTT; TRCN0000005998, Sigma). Lentivirus particles were produced in HEK-293T cells transfected with psPAX2, PMD2G (Addgene), and pLKO.1-shRNA-containing plasmids. M213 and HUCCT1 cells were infected with the lentivirus, and stable cells were established by puromycin selection (Sigma). Migration assays were performed in Transwells (Corning Inc., 8.0-μm pore size). For migration, 2.5 × 10⁴ cells of RASA1 stably silenced M213 and HUCCT1 cells in serum-free medium were added to 24-well insert plates. Media containing 10% FBS were added to the lower wells and incubated for 5 hours. For cell invasion assays, the filters were precoated with Matrigel. RASA1 stably silenced M213 and HUCCT1 cells (2.5 × 10⁴) in serum-free medium were added into 24-well insert plates. Media containing 10% FBS were added to the lower well of the chambers and incubated for 5 hours. After incubation, the cells on the upper surface of the filter were completely removed and migrated cells were trypsinized and counted.

We aligned sequence data to the human reference genome (hs37d5) using BWA-MEM v0.7.9a (54). We removed PCR duplicates using SAMTools (55) and performed indel realignments using Genome Analysis Toolkit v1.0 (GATK; ref. 56). We used realigned data as input to both GATK Unified Genotyper and MuTect to call sSNVs. We used GATK IndelGenotyperV2 to identify indels. We applied filters and manual inspection to retain only high confidence sSNVs and indels. Further details are provided in the Supplementary Methods.

FIREFLY is a method for identifying gene sets dysregulated by somatic promoter mutations through modulation of TF binding. We extracted somatic promoter mutations by selecting noncoding sSNVs within ±2 kb of transcription start sites of GENCODE genes. We identified those mutations predicted to change TF-binding based on PBM data (details below). We then evaluated gene sets for enrichment in promoter binding-change mutations. For each gene set, we calculated the test statistic M = number of genes in the gene set with binding-change mutations, summed across all tumors. To identify gene sets enriched in binding-change mutations, we performed two statistical tests in sequence, where gene sets found significant in each test (q < 0.1) were then tested in the next: (i) Fisher exact test; (ii) a synthetic mutations test. For gene sets passing these tests, we then performed gene expression tests using the Gene Set Analysis (GSA) v1.03 R package (37) to assess their transcriptional dysregulation. Further details are provided in the Supplementary Methods.

We generated synthetic mutations in promoter regions of interest based on the genome-wide frequency of mutations at each trinucleotide observed for each tumor/normal pair. Synthetic mutations were generated separately for each tumor, based on the frequencies and mutation types observed in the actual tumor. Further details are provided in the Supplementary Methods.

To identify binding-change mutations, we used TF-DNA binding specificity PBM data from the cis-BP database (34). PBM experiments measure TF binding to all possible 8-mer sequences with a DNA binding enrichment score (E-score; ref. 35). Typically, E-scores > 0.35 correspond to specific TF-DNA binding (57). To call binding sites for a particular TF, we required that such sites contain at least 2 consecutive 8-mers with E-scores > 0.4. We also used PBM data to call “nonbinding sites,” defined as genomic regions containing only 8-mers with E-score < 0.3. For each somatic mutation and each TF with available PBM data, we analyzed the 15-bp genomic region centered at the mutation. If the region contained a TF binding site in the normal sample but not in the corresponding tumor sample, the mutation was called a “loss-of-binding” mutation for that TF. If the region contained a TF binding site in the tumor sample but not in the normal sample, then the mutation was called “gain-of-binding” for that TF. We describe a mutation as “binding change” if it is either loss-of-binding or gain-of-binding for any of the interrogated TFs. Further details are provided in the Supplementary Methods.

We integrated (i) 71 WGS CCAs, (ii) 188 targeted-sequenced CCAs (Supplementary Table S5), and (iii) 200 exome-sequenced CCAs (8), and performed gene significance analyses using MutSigCV (58) and IntOGen (59). For input, we used the list of all coding sSNVs and indels found in the 404 targeted genes across 459 samples (including both silent and nonsilent mutations). Both tools were run with default parameters and we retained genes found significant by both tools with q values < 0.1. Further details are provided in the Supplementary Methods.

BWA-MEM alignments from each tumor–normal pair were analyzed by CREST (Clipping REveals STructure; ref. 60) and PTRfinder (61). For most tumors we required ≥3 uniquely mapped split-read alignments at each SV breakpoint; for the shallower Japanese WGS data we required only ≥5 such alignments over both SV breakpoints combined. We considered a tumor SV to be somatic if no SV in the normal sample occurred within one-half of a read length from the tumor SV. Further details are provided in the Supplementary Methods.

We searched for sources of somatic L1 insertions by looking for highly recurrent SVs: ≥10 SVs in a 1 Mb region. We then selected the subset of these regions that contained a mobile L1 element in a database of retrotransposon insertion polymorphisms (dbRIP; ref. 62). Only SVs with ≥2 reads with poly-A tails at the putative L1 insertion site were retained for further analysis.

For genomic DNA, 100 ng of whole-genome amplified DNA of the tumor and normal matched cases were used as PCR templates. For cDNA, total cDNAs of tumor and normal matched control were synthesized using SuperScript III System according to the manufacturer's instructions (Invitrogen) and 40 ng of cDNA were used as PCR template. PCR was performed using fusion-specific primers with Platinum Taq DNA Polymerase system (Invitrogen). PCR products were cleaned up by the Exo/Sap enzyme system (Invitrogen) and bidirectionally sequenced using the BigDye Terminator v.3.1 Kit (Applied Biosystems) and an ABI PRISM 3730 Genetic Analyzer (Applied Biosystems). Sequencing traces were aligned to reference sequences using Lasergene 10.1 (DNASTAR) and analyzed by visual inspection.

Somatic L1 insertions were validated by PCRs using primers flanking predicted sites of insertion. PCRs were performed using AccuPrime Pfx DNA Polymerase (Invitrogen) with 200 ng of WGA DNA.

Raw SNP array data was processed using Illumina Genome Studio. We used ASCAT v2.0 (63) to estimate allele-specific copy-number profiles. We determined regions of copy-number alteration based on their relative copy number using the “copy-number” R package (64). For tumor–normal pairs without SNP array data, we estimated copy-number profiles based on sequencing data using Control-FREEC (65), Quandico (66), and Sequenza (67). We used GISTIC v2.0.22 (68) to determine regions of significant focal copy-number alterations, using ASCAT/Sequenza's inferred copy-number segments, and associated copy-number values were defined as log₂ of the segment's relative copy number. Full details are provided in the Supplementary Methods.

Gene expression microarray data were preprocessed using the “lumi” R package (69). Batch effects were removed using ComBat (70). We used GSEA v2.2.2 (71) with a classic weighting scheme to determine pathways upregulated or downregulated in each integrative CCA cluster relative to the others, employing canonical pathways in the MSigDB C2 catalogue of annotated gene sets. Full details are provided in the Supplementary Methods.

ESTIMATE (72) was used to determine the presence of infiltrating immune cells, using the ImmuneSignature geneset. A total of 126 genes were used to determine the immune score for each tumor.

DNA methylation profiles were obtained for 138 tumors and 4 normal samples. Data were preprocessed using the “minfi” (73) and “wateRmelon” (74) R packages. We selected 4,520 probes with the highest standard deviations in β-values across the tumors, and mean β < 0.5 in the normal samples, for clustering using the “RPMM” R package (75). In the hypermethylated methylation clusters (1 and 4), we considered a CpG site to be hypermethylated if the following conditions held: (1) β < 0.5 in normal samples; (2) M values were significantly different in the (i) hypermethylated cluster versus (ii) the combined normal samples and the low-methylation tumors—those not in methylation cluster 1 or 4 (q < 0.05, two-sided t test); and (3) its mean β in the hypermethylated cluster minus the mean β across the normal samples and low-methylation tumors was >0.2.

To explore associations between mutation signatures and hypermethylated CpGs, we considered only mutations located within 50 bp of CpG probes that had mean β < 0.5 in normal samples. In each tumor, the nearest CpG probe to a mutation was considered to be hypermethylated if (i) it was hypermethylated in that tumor's methylation cluster; (ii) its individual β was > 0.5; and (iii) its individual β minus the mean β across the normal samples and low-methylation tumors was >0.2. Other analyses were similar to community-standard analyses. Further details are provided in the Supplementary Methods.

The “iClusterPlus” R package (76) was used to perform integrative unsupervised clustering of 94 CCAs based on 4 genomic data types: (i) somatic point mutations in 404 targeted genes (gene by sample matrix of binary values), (ii) sCNAs defined as copy-number segments identified by ASCAT v2.0, (iii) the most variable expression probes (coefficient of variation > 0.1), and (iv) the most variable methylation probes (top 1% standard deviation in β value). We ran iClusterPlus.tune with different numbers of possible clusters (n = 2–7), choosing the number of clusters at which the percentage of explained variation leveled off (n = 4), and the clustering with the lowest Bayesian information criterion. Further details are provided in the Supplementary Methods.

The “survival” R package was used to perform survival analysis using Kaplan–Meier statistics, with P values computed by log-rank tests. Multivariate survival analysis was performed using the Cox proportional hazards method. To validate our survival analysis results, we also analyzed a separate validation cohort of 58 samples by combining two sources: (i) 25 samples (with survival data) newly classified into CCA clusters under the expanded integrative clustering; and (ii) 33 recently published CCA samples with survival data from Farshidfar and colleagues (2017; ref. 13). Further details are provided in the Supplementary Methods.

Nonnegative matrix factorization (NMF) was applied to the trinucleotide-context mutation spectra of CCAs to extract mutation signatures. Six stable and reproducible extracted mutational signatures were compared with the 30 signatures from COSMIC (http://cancer.sanger.ac.uk/cosmic/signatures) based on cosine similarities. We used supervised NMF to evaluate the contributions of the COSMIC signatures to mutations in CCA. Additional signatures were considered via visual inspection. We ignored signatures that contributed <5% of the total mutations in a particular tumor, and removal of these signatures did not substantially increase reconstruction errors. The MSI status in the prevalence set was determined by the indel counts in simple repeat sequences. Further details are provided in the Supplementary Methods.

The whole genome and targeted sequencing data done in this article have been deposited at the European Genome-phenome Archive (EGA; http://www.ebi.ac.uk/ega) under accession numbers EGAS00001001653 and at the International Cancer Genome Consortium Data Portal database (https://dcc.icgc.org/). The WES data (8) has been deposited at EGA under accession number EGA00001000950. Gene expression and methylation data have been deposited at the Gene Expression Omnibus (GEO; www.ncbi.nlm.nih.gov/geo) with GEO accession numbers GSE89749 and GSE89803, respectively.

S.P. Choo reports receiving a commercial research grant from Bristol-Myers Squibb and is a consultant/advisory board member for Novartis, Bristol-Myers Squibb, and Celgene. P. Tan is a consultant/advisory board member for Aslan Therapeutics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A. Jusakul, I. Cutcutache, J.Q. Lim, N. Padmanabhan, V. Nellore, D.G. Duda, S.G. Rozen, T. Shibata, C. Pairojkul, B.T. Teh, P. Tan

Development of methodology: A. Jusakul, I. Cutcutache, C.H. Yong, J.Q. Lim, V. Nellore, S. Kongpetch, C.C.Y. Ng, R. Gordan, S.G. Rozen, P. Tan

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Jusakul, N. Padmanabhan, L.M. Ng, S.P. Choo, C.K. Ong, S. Lie, A.S.T. Lim, T.H. Lim, J. Tan, N. Khuntikeo, V. Bhudhisawasdi, P. Yongvanit, S. Wongkham, Y. Arai, S. Yamasaki, P.K.H. Chow, A.Y.F. Chung, L.L.P.J. Ooi, K.H. Lim, S. Dima, I. Popescu, P. Broet, S.-Y. Hsieh, M.-C. Yu, A. Scarpa, J. Lai, D.-X. Luo, A.L. Carvalho, A.L. Vettore, H. Rhee, Y.N. Park, S.G. Rozen, T. Shibata, B.T. Teh

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A. Jusakul, I. Cutcutache, C.H. Yong, J.Q. Lim, M.N. Huang, N. Padmanabhan, V. Nellore, S. Kongpetch, A.W.T. Ng, L.M. Ng, S.P. Choo, S. Nagarajan, W.K. Lim, A. Boot, M. Liu, J.R. McPherson, Y. Totoki, H. Nakamura, K.H. Lim, S.-Y. Hsieh, L. Alexandrov, R. Gordan, S.G. Rozen, T. Shibata, P. Tan

Writing, review, and/or revision of the manuscript: A. Jusakul, I. Cutcutache, C.H. Yong, J.Q. Lim, M.N. Huang, V. Nellore, A.W.T. Ng, S.P. Choo, M. Liu, J.R. McPherson, P.K.H. Chow, L.L.P.J. Ooi, S. Dima, D.G. Duda, P. Broet, A. Scarpa, A.L. Carvalho, Y.N. Park, R. Gordan, S.G. Rozen, T. Shibata, C. Pairojkul, B.T. Teh, P. Tan

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Jusakul, C.H. Yong, J.Q. Lim, L.M. Ng, S.S. Myint, R. Thanan, C.C.Y. Ng, V. Rajasegaran, J. Tan, J.L. Loh, J.R. McPherson, P. Yongvanit, Y. Totoki, K.H. Lim, J. Lai, D.-X. Luo, Y.N. Park, S.G. Rozen, C. Pairojkul, P. Tan

Study supervision: V. Bhudhisawasdi, P.K.H. Chow, S.-Y. Hsieh, J. Lai, R. Gordan, S.G. Rozen, C. Pairojkul, B.T. Teh, P. Tan

Other (conception and design of methodology): R. Gordan

This study was performed as part of the International Cancer Genome Consortium (ICGC). We thank the Duke-NUS Genome Biology Facility for methylation and gene expression assays, the Genome Institute of Singapore Population Genetics facility for SNP arrays, and the Cytogenetics Laboratory, Department of Molecular Pathology, Singapore General Hospital for FISH analysis. We thank D. Jefferson (New England Medical Center, Tufts University) for H69 cells. We also thank the SingHealth Tissue Repository, Dr. Catherine Guettier, and Dr. Jean-Charles Duclos-Vallée (DHU Hepatinov, Hôpital Paul Brousse, AP-HP, Villejuif, France) for tissue samples.

We thank Sir Lambert Cornelias Bronsveld and Lady Bronsveld-Ngo Kim Lian for philanthropic support (B.T. Teh). This work was supported by the Singapore National Medical Research Council [NMRC/STaR/0006/2009 (B.T. Teh), NMRC/STaR/0024/2014 (B.T. Teh), NMRC/CG/012/2013 (B.T. Teh), NMRC/CIRG/1422/2015 (S.G. Rozen), and NMRC/STaR/0026/2015 (P. Tan)], Genome Institute of Singapore (P. Tan), Duke-NUS Medical School (P. Tan, S.G. Rozen), National University of Singapore, the National Research Foundation Singapore, Singapore Ministry of Education under the Research Centres of Excellence initiative (P. Tan), National University Cancer Institute, Singapore (NCIS) Centre Grant NR 13NMR 111OM (P. Tan), the Japan Agency for Medical Research and Development (Practical Research for Innovative Cancer Control, 15ck0106094h0002; T. Shibata), National Cancer Center Research and Development Funds (26-A-5; T. Shibata), Italian Cancer Genome Project (FIRB RBAP10AHJB; A. Scarpa), Associazione Italiana Ricerca sul Cancro (n. 12182; A. Scarpa), Italian Ministry of University Health (FIMP CUP_J33G13000210001; A. Scarpa), the NCI of the NIH (P01CA142538; R. Gordan), and the National Natural Science Foundation of China (NSFC81372825; J. Lai).