Pediatric liver cancers (PLC) comprise diverse diseases affecting infants, children, and adolescents. Despite overall good prognosis, PLCs display heterogeneous response to chemotherapy. Integrated genomic analysis of 126 pediatric liver tumors showed a continuum of driver mechanisms associated with patient age, including new targetable oncogenes. In 10% of patients with hepatoblastoma, all before three years old, we identified a mosaic premalignant clonal expansion of cells altered at the 11p15.5 locus. Analysis of spatial and longitudinal heterogeneity revealed an important plasticity between “hepatocytic,” “liver progenitor,” and “mesenchymal” molecular subgroups of hepatoblastoma. We showed that during chemotherapy, “liver progenitor” cells accumulated massive loads of cisplatin-induced mutations with a specific mutational signature, leading to the development of heavily mutated relapses and metastases. Drug screening in PLC cell lines identified promising targets for cisplatin-resistant progenitor cells, validated in mouse xenograft experiments. These data provide new insights into cisplatin resistance mechanisms in PLC and suggest alternative therapies.

Significance:

PLCs are deadly when they resist chemotherapy, with limited alternative treatment options. Using a multiomics approach, we identified PLC driver genes and the cellular phenotype at the origin of cisplatin resistance. We validated new treatments targeting these molecular features in cell lines and xenografts.

This article is highlighted in the In This Issue feature, p. 2355

Pediatric liver cancers (PLC) are rare tumors and therefore have not been molecularly characterized in large series. Hepatoblastomas represent about 67% to 80% of all PLCs worldwide, generally developing before five years of age on nonfibrotic liver (1, 2). Although some rare syndromes such as familial adenomatous polyposis, Beckwith–Wiedemann syndrome, or Simpson-Golabi Behmel syndrome predispose to hepatoblastoma development, the etiology of hepatoblastomas is poorly understood because most of these tumors are sporadic. Hepatoblastomas are characterized by their histologic heterogeneity, with three main histology patterns—fetal, embryonal, and mesenchymal—that often coexist within a single tumor. A handful of genomics studies (3–7) have established hepatoblastomas as genetically simple tumors, with the smallest mutation burden among 24 pediatric cancer types (8). Beyond activating CTNNB1 (β-catenin) alterations found in most hepatoblastomas (70%–90%), only few recurrent driver mutations have been described, including NFE2L2 (5%–10%) and TERT promoter mutations (2%–5%). Pediatric hepatocellular carcinomas (HCC; incidence = 0.24–0.65 per 1,000,000) resemble adult HCC in their histology and frequently develop on fibrotic/cirrhotic liver as a consequence of hepatitis B virus infection or rare congenital disorders (1, 9). Previous genomic studies identified recurrent alterations in WNT signaling and telomerase pathways, but no oncogenic driver in pediatric HCC with underlying liver disease (10, 11). Fibrolamellar carcinomas (FLC), a rare subtype of HCC, develop in adolescents and young adults on healthy liver and have been characterized by a recurrent DNAJB1PRKACA driver gene fusion (12). Finally, benign lesions such as hepatocellular adenomas (HCA) and focal nodular hyperplasias are usually related to congenital malformation associated with vascular abnormalities or metabolic genetic diseases, or occur after chemotherapy (2, 13).

Hepatoblastomas are usually treated with cisplatin-based neoadjuvant chemotherapy and subsequent surgical removal of the tumor, leading to >80% five-year survival (14, 15). However, some hepatoblastomas develop resistance to chemotherapy during the initial neoadjuvant chemotherapy or after tumor recurrence, and the molecular determinants of cisplatin resistance are yet to be discovered. In contrast to hepatoblastomas, pediatric HCCs respond poorly to chemotherapy, and as in adults, they have a poor prognosis if not completely removed by surgery.

We aimed to establish the detailed driver landscape of 126 pediatric liver tumors, then we analyzed the plasticity of hepatoblastoma tumors in relation to cisplatin resistance, and explored new therapeutic strategies to overcome this resistance.

The Driver Landscape of PLCs

We analyzed a cohort of 126 pediatric liver tumors comprising 104 hepatoblastomas that developed in 65 patients (87 primary tumors including 18 prechemotherapy and 17 relapses/metastases), 10 HCCs (9 patients), 7 FLCs and 5 HCAs sequenced by whole-genome sequencing (WGS; n = 65), whole-exome sequencing (WES; n = 57), RNA sequencing (RNA-seq; n = 120), and reduced representation bisulfite sequencing (RRBS; n = 92; Supplementary Fig. S1; Supplementary Table S1). Among the 65 patients with hepatoblastoma, 14 were older than 5 years of age at diagnosis. HCC developed in fibrotic or cirrhotic liver related to various constitutional liver diseases (tyrosinemia, mitochondrial cytopathy, progressive familial intrahepatic cholestasis, associated with germline mutations of FAH, NDUFA11, NDUFB9, TJP2, and ABCB11; Supplementary Table S2) whereas the other tumors arose on normal liver. Germline truncating mutations of BRCA1 and BRCA2 were also identified in two hepatoblastomas and one FLC. In WGS analyses, primary hepatoblastomas displayed a small number of somatic mutations (median = 886, 0.3 mutations/Mb), whereas hepatoblastoma metastases and relapses showed a massive mutation load (median = 12,824, 4.3 mutations/Mb) with a high proportion of doublet-base substitutions (8%; Fig. 1A). Overall, HCC displayed a higher mutation rate than hepatoblastomas (median 5,318, P = 0.002).

Figure 1.

Genomic landscape of PLCs. A, Number of somatic mutations and SV identified in 65 PLCs by WGS. Samples are ordered by diagnosis and mutation burden. Alteration types are indicated with a color code, and specific structural variant phenotypes are highlighted. B, Heat map representation of driver alterations across 122 PLCs analyzed by WGS or WES. Samples are ordered by diagnosis and patient age at sampling. C, Frequency of copy-number alterations along the genome for 122 PLCs analyzed by WGS or WES. The top axis indicates the frequency of low-amplitude changes (gains, losses, and LOH); the bottom axis indicates the frequency of high-amplitude changes (focal amplifications and homozygous deletions). Target genes of amplifications and homozygous deletions are indicated. Correlation between gene expression (variance-stabilized) and copy number is displayed for three selected genes. HB, hepatoblastoma.

Figure 1.

Genomic landscape of PLCs. A, Number of somatic mutations and SV identified in 65 PLCs by WGS. Samples are ordered by diagnosis and mutation burden. Alteration types are indicated with a color code, and specific structural variant phenotypes are highlighted. B, Heat map representation of driver alterations across 122 PLCs analyzed by WGS or WES. Samples are ordered by diagnosis and patient age at sampling. C, Frequency of copy-number alterations along the genome for 122 PLCs analyzed by WGS or WES. The top axis indicates the frequency of low-amplitude changes (gains, losses, and LOH); the bottom axis indicates the frequency of high-amplitude changes (focal amplifications and homozygous deletions). Target genes of amplifications and homozygous deletions are indicated. Correlation between gene expression (variance-stabilized) and copy number is displayed for three selected genes. HB, hepatoblastoma.

Close modal

The WNT/β-catenin pathway was the most frequently altered oncogenic pathway in PLCs (84.5%), with different activating mechanisms across diagnoses. In hepatoblastomas, CTNNB1 alterations activating β-catenin were identified in 92% of the tumors (Supplementary Tables S3 and S4), with missense mutations exclusively observed in young patients (<4 years), whereas exon 3 in-frame deletions were observed later in life (Fig. 1B; Supplementary Fig. S2). The remaining five patients with hepatoblastoma without CTNNB1 mutation showed germline truncating mutations of APC (n = 3) or AXIN1 (n = 1) with somatic inactivation of the second allele in the tumor, and only one hepatoblastoma remained without an identified alteration in the pathway. CTNNB1 was also altered in 40% of HCAs. In HCC, no CTNNB1 mutations were detected, but biallelic inactivation of AXIN1 and/or AMER1 occurred in six of nine patients (67%).

The 11p15.5 imprinted locus, containing the IGF2 oncogene, was the second most frequently altered locus in hepatoblastomas and HCC (84% and 89%, respectively), mostly through copy-neutral loss of heterozygosity (cn-LOH; 51%–56%; Fig. 1C). Adding to LOH, we found epimutations of the imprinting control regions IC1 [gain of methylation (GOM), 22% of hepatoblastomas and 33% of HCC] and IC2 [loss of methylation (LOM), 5% of hepatoblastomas], and recurrent somatic mutations of CDKN1C in four patients with hepatoblastoma. In six patients with hepatoblastoma (10%), imbalanced B allele frequency (BAF) profiles revealed that the cn-LOH found in tumor cells was present as a mosaic in the adjacent normal liver (Fig. 2A and B). These young patients (median = 8.4 months, P = 0.045) were not diagnosed with Beckwith–Wiedemann syndrome (BWS). Yet, the cn-LOH of 11p15.5 was detected in a significant fraction of normal liver cells (6%–58%) without the other driver alterations identified in the corresponding tumor, indicating premalignant clonal expansions of hepatocytes. Oncogenic transformation involved the acquisition of CTNNB1 missense activating mutation in the six cases (Fig. 2B). Accordingly, IGF2 was highly overexpressed in the nontumor liver samples of two patients harboring a mosaic 11p15.5 alteration, whereas WNT/β-catenin target genes were overexpressed only in tumor cells (Fig. 2C). In patient #3559 (30% mosaic cn-LOH in the nontumor liver), RNAscope in situ hybridization revealed a massive overexpression of IGF2 in both the tumor and adjacent nontumor tissue, whereas glutamine synthetase and β-catenin immunostainings demonstrated an oncogenic activation of the WNT/β-catenin pathway specifically in tumor cells (Fig. 2D). These data support the idea that premalignant clonal expansions of “normal” hepatocytes with 11p15.5 cn-LOH overexpress IGF2 and can lead to hepatoblastoma formation after oncogenic CTNNB1 mutation. Finally, one patient with BWS (#3180) displayed IC2 mosaic epimutation in both blood cells and the liver (27% of cells). Overall, the mechanism of 11p15.5 alteration in PLC was related to age, with cn-LOH occurring more frequently in young patients and gain of methylation of IC1 mostly observed in older patients (Fig. 1B; Supplementary Fig. S2).

Figure 2.

Premalignant clonal expansions with 11p15 alteration in patients with hepatoblastoma. A, Identification of a copy-neutral LOH (cn-LOH) at 11p15 locus in the nontumor (NT) liver of patient #3559. BAFs of heterozygous single-nucleotide polymorphisms (SNP) are represented along chromosome 11. SNPs with a BAF greater (respectively lower) than 0.5 in the tumor are colored in red (respectively blue). In the cn-LOH region, red (respectively blue) SNPs correspond to those for which the B allele was retained (respectively lost). The same BAF imbalance is identified in the NT sample, with the same boundaries, demonstrating the presence of the cn-LOH. The amplitude of BAF changes indicates that the cn-LOH is present in 30% of cells in the NT sample. B, Premalignant expansions with cn-LOH at 11p15 locus were identified in six patients with hepatoblastoma. The proportion of NT cells harboring the alteration is indicated below. Copy-neutral LOH became clonal in matched hepatoblastoma that had acquired in addition activating CTNNB1 mutations. C, Expression levels of IGF2 and the β-catenin targets GLUL, LGR5, and AXIN2 (2−ΔΔCt) in tumor and NT samples from patients with and without 11p15.5 cn-LOH premalignant expansions. D, Representative areas of FFPE slides from two patients with hepatoblastoma: patient #3559 with mosaic 11p15.5 cn-LOH and control patient #4217. Four types of staining were performed: hematoxylin and eosin staining (H&E), IGF2 RNAscope in situ hybridization, and β-catenin and glutamine synthetase immunostainings.

Figure 2.

Premalignant clonal expansions with 11p15 alteration in patients with hepatoblastoma. A, Identification of a copy-neutral LOH (cn-LOH) at 11p15 locus in the nontumor (NT) liver of patient #3559. BAFs of heterozygous single-nucleotide polymorphisms (SNP) are represented along chromosome 11. SNPs with a BAF greater (respectively lower) than 0.5 in the tumor are colored in red (respectively blue). In the cn-LOH region, red (respectively blue) SNPs correspond to those for which the B allele was retained (respectively lost). The same BAF imbalance is identified in the NT sample, with the same boundaries, demonstrating the presence of the cn-LOH. The amplitude of BAF changes indicates that the cn-LOH is present in 30% of cells in the NT sample. B, Premalignant expansions with cn-LOH at 11p15 locus were identified in six patients with hepatoblastoma. The proportion of NT cells harboring the alteration is indicated below. Copy-neutral LOH became clonal in matched hepatoblastoma that had acquired in addition activating CTNNB1 mutations. C, Expression levels of IGF2 and the β-catenin targets GLUL, LGR5, and AXIN2 (2−ΔΔCt) in tumor and NT samples from patients with and without 11p15.5 cn-LOH premalignant expansions. D, Representative areas of FFPE slides from two patients with hepatoblastoma: patient #3559 with mosaic 11p15.5 cn-LOH and control patient #4217. Four types of staining were performed: hematoxylin and eosin staining (H&E), IGF2 RNAscope in situ hybridization, and β-catenin and glutamine synthetase immunostainings.

Close modal

A paucity of structural variants (SV) was observed in primary hepatoblastomas (median = 6/tumor) except for six tumors with chromoplexy, which mostly developed in older patients (Fig. 1A; Supplementary Table S5). In HCC, we identified an unusual number of focal deletions (median = 47/tumor) of small size (median = 27 kb). This remarkable deletor phenotype led to recurrent cancer driver alterations, particularly on chromosome X with complete inactivation of AMER1 (44%), GPC3 (44%), RPS6KA3 (33%), SMARCA1 (22%), and BCORL1 (11%), but also on chromosome 16 with homozygous deletion or combined deletion and truncating mutation of AXIN1 (44%) and CREBBP (11%; Fig. 1B and C; Supplementary Fig. S3). Of note, GPC3, BCORL1, CREBBP, RPS6KA3, and AXIN1 were also inactivated in hepatoblastomas through truncating or damaging mutations, including one germline GPC3 mutation in a patient with Simpson-Golabi syndrome. Recurrent copy-number alterations were identified in two genes controlling p53 degradation: focal deletions at 4q35 pinpointed inactivation of the tumor suppressor IRF2 in 30% of hepatoblastomas and 56% of HCC, whereas focal amplifications at 1q32.1 led to a high expression of MDM4 in four patients with hepatoblastoma (Fig. 1C; Supplementary Fig. S3; Supplementary Table S6). Interestingly, gains of the entirety of 1q were also observed in 50% of hepatoblastomas and HCC (Fig. 1C). At the chromosome arm level, HCC and hepatoblastomas had a roughly similar profile of gains but HCC harbored more losses, including loss of the 13q arm (encompassing RB1) in 33% of cases (Supplementary Fig. S4). Finally, focal amplifications of CCND1/FGF19 were found in three patients with hepatoblastoma, leading to the overexpression of both genes. Other recurrent driver mutations in hepatoblastomas involved TERT promoter in nine patients (all older than 40 months, P = 1.5 × 10−5; Supplementary Fig. S2), NFE2L2 (four patients) and ARID1A (two patients). No recurrent gene fusion was identified in the cohort except the PRKACADNAJB1 fusion pathognomonic of FLC. We also identified recurrent mutations of ERBB4 in 2/7 FLCs (29%) and of HNF1A in 3 HCAs (60%), including one germline mutation. Overall, hepatoblastomas and HCC shared common pathways altered by diverse genes and mechanisms, while HCA and FLC had specific driver alterations (Fig. 1B).

Phenotypic Plasticity of Hepatoblastoma Cells across Three Differentiation States

Unsupervised transcriptomic classification of 100 hepatoblastoma samples from 64 patients revealed 4 robust molecular groups characterized by diverse differentiation states, cell proliferation, and immune infiltration levels (Fig. 3A). The differentiated “hepatocytic” hepatoblastoma group comprised 44 samples, 73% of which had only “fetal” histologic component in their corresponding mirror block (vs. 10%,P = 1.5 × 10−10). They strongly expressed transcription factors (TF) involved in hepatic differentiation (HNF1A and HNF4A; Fig. 3A). “Hepatocytic” hepatoblastomas were divided into two clusters: the “hepatocytic hot” subgroup was defined by a strong signature of polymorphic immune infiltration, including a mixture of T, B, and natural killer (NK) cells and macrophages, with 20/21 of these samples collected after chemotherapy; the “hepatocytic cold” subgroup was characterized by lower levels of immune infiltrates, and it was enriched in prechemotherapy samples (39%,P = 0.010; Fig. 3A; Supplementary Fig. S5A–S5C). The “liver progenitor” group comprised 44 samples enriched in highly proliferative, immune-cold tumors with embryonal compartments, 11p15.5 alterations, and MDM4 amplifications. They expressed TFs involved in hepatic differentiation but also TFs involved in self-renewal and pluripotency maintenance (MYCN, MIXL1, and LIN28B). The last group, “mesenchymal” hepatoblastomas, included all samples with a mesenchymal histology that were derived from mixed hepatoblastomas with both epithelial and mesenchymal components, and was found mostly in younger patients (Supplementary Fig. S2). These mesenchymal samples displayed a distinct differentiation program, with no expression of liver differentiation genes but a strong expression of mesenchymal stem cell TFs (TWIST1, TBX5, and MSX2), and variable levels of immune infiltration. We validated these four transcriptomic subgroups and their associations with histology and immune infiltration in an independent RNA-seq cohort of 34 hepatoblastoma samples (ref. 16; Supplementary Fig. S6A–S6C).

Figure 3.

Molecular plasticity of hepatoblastoma across three differentiation states. A, Gene expression–based classification of hepatoblastoma. Unsupervised hierarchical clustering of 100 hepatoblastoma (HB) samples from 64 patients and 4 nontumor liver samples revealed 4 molecular groups. Clinical and molecular annotations are depicted below the dendrogram, with P values indicating their association with molecular groups. A heat map shows the expression of key TF and marker genes representative of each group, as well as molecular scores of hepatic differentiation, cell proliferation, and immune infiltration. B, Top, proportion of transcriptomic groups and a selection of driver alteration frequencies at different steps of hepatoblastoma progression. Middle, transcriptomic group switches identified in 24 patients with multiple sampling, including patients with pre/post-chemotherapy samples, synchronous samples at distinct locations in the primary tumor, and/or paired primary and relapse/metastasis. The number of patients with/without molecular switch is indicated for each type of multiple sampling, and the transcriptomic switch is represented by a color code on the arrows. Bottom, examples of histologic heterogeneity matching transcriptomic group switches. H, hepatocytic; M, mesenchymal; LP, liver progenitor. C, Projection of hepatoblastoma and nontumor liver samples over two independent methylation components. Hepatoblastomas are colored by their transcriptomic group, and samples from a same patient are linked with black lines. D, Association of intrasample histologic heterogeneity with transcriptomic groups. E, Association of intratumor histologic heterogeneity with 11p15.5 locus alteration.

Figure 3.

Molecular plasticity of hepatoblastoma across three differentiation states. A, Gene expression–based classification of hepatoblastoma. Unsupervised hierarchical clustering of 100 hepatoblastoma (HB) samples from 64 patients and 4 nontumor liver samples revealed 4 molecular groups. Clinical and molecular annotations are depicted below the dendrogram, with P values indicating their association with molecular groups. A heat map shows the expression of key TF and marker genes representative of each group, as well as molecular scores of hepatic differentiation, cell proliferation, and immune infiltration. B, Top, proportion of transcriptomic groups and a selection of driver alteration frequencies at different steps of hepatoblastoma progression. Middle, transcriptomic group switches identified in 24 patients with multiple sampling, including patients with pre/post-chemotherapy samples, synchronous samples at distinct locations in the primary tumor, and/or paired primary and relapse/metastasis. The number of patients with/without molecular switch is indicated for each type of multiple sampling, and the transcriptomic switch is represented by a color code on the arrows. Bottom, examples of histologic heterogeneity matching transcriptomic group switches. H, hepatocytic; M, mesenchymal; LP, liver progenitor. C, Projection of hepatoblastoma and nontumor liver samples over two independent methylation components. Hepatoblastomas are colored by their transcriptomic group, and samples from a same patient are linked with black lines. D, Association of intrasample histologic heterogeneity with transcriptomic groups. E, Association of intratumor histologic heterogeneity with 11p15.5 locus alteration.

Close modal

Strikingly, 14/24 patients (58%) with multiple synchronous and/or metachronous samples displayed transcriptomic group switches (Fig. 3B; Supplementary Fig. S7). Spatial transcriptomic heterogeneity was identified in the primary tumors of seven patients and matched histologic heterogeneity between fetal and embryonal parts of the tumor (Fig. 3B). Longitudinal transcriptomic group changes were also identified between pre- and post-chemotherapy samples (4 out of 6 patients) and between primary tumors and relapses/metastases (5 out of 8 patients). Phylogenetic trees revealed private driver alterations in some cases, but no recurrent gene associated with specific transcriptional group changes. In contrast, DNA methylation profiles were closely associated with differentiation states, and transcriptomic group changes were associated with DNA methylation reprogramming (Fig. 3C; Supplementary Fig. S8A). In particular, mesenchymal samples displayed coordinated hypermethylation of HNF4 and PPAR binding sites, both in our cohort (Supplementary Fig. S8B and C) and in an independent data set (ref. 16; Supplementary Fig. S6B).

The molecular plasticity of hepatoblastomas recapitulated the heterogeneity observed at the histologic level (Fig. 3B). A systematic histologic review at the sample (mirror block) and whole tumor levels revealed that 80% of primary hepatoblastomas displayed spatial heterogeneity with a mixture of embryonal, fetal, or mesenchymal areas. The “liver progenitor” molecular group was associated with increased intrasample histologic heterogeneity (P = 1.0 × 10−6; Fig. 3D), with frequent coexistence of fetal and embryonal cells (Fig. 3A and B). Primary tumors with alterations in the 11p15.5 locus also displayed more histologic heterogeneity (89% vs. 46%,P = 0.0021; Fig. 3E). Thus, the phenotypic plasticity of hepatoblastomas may relate to the multipotency of progenitor cells with 11p15.5 alteration.

Consistent with this plasticity, the proportion of molecular groups evolved across disease stages (P = 0.0014; Fig. 3B). The “hepatocytic hot” group was enriched in post- versus prechemotherapy samples (30% vs. 6%), whereas the “liver progenitor” group was enriched in metastases and relapses (64% vs. 40%). WNT/β-catenin alterations were always trunk in phylogenetic trees and already present in 100% of prechemotherapy primary samples. In contrast, 11p15 alterations occurred late in 8 out of the 20 affected phylogenetic trees and their frequency increased in later disease stages, which was also the case for IRF2 deletions (Fig. 3B).

To validate this phenotypic plasticity at the single-cell level, we analyzed matched nontumor liver, primary hepatoblastoma, and lung metastasis from one patient (#3981) by single-nucleus RNA-seq (snRNA-seq). Virtual copy-number profiles were consistent with the copy-number changes observed in bulk WES data, with gains of 1q, 2q, 5q, and 15q identified in all tumor cells, and an extra gain of Xq in all metastatic cells (Fig. 4A). In contrast to this genetic homogeneity, we identified clusters of tumor cells with distinct transcriptomic profiles, corresponding to the “mesenchymal,” “liver progenitor,” and “hepatocytic” signatures (Fig. 4B). The three populations were present in different proportions in each sample, with a larger “mesenchymal” contingent in the primary tumor, consistent with the molecular classification of bulk samples (“mesenchymal” for the primary, “liver progenitor” for the metastasis). These data demonstrate the existence of the three differentiation states at the single-cell level, and the plasticity of hepatoblastoma cells across these differentiation states. Indeed, both the last common ancestor (LCA) of the primary tumor and the LCA of the metastasis (with additional Xq gain) were able to generate the three cell types repeatedly during tumor evolution in this patient (Fig. 4C).

Figure 4.

Single-nucleus RNA-seq reveals molecular plasticity along tumor progression in one patient. A, Virtual copy-number profiles discriminate tumor and normal cells. In agreement with WES data, all tumor cells display gains at 1q, 2q, 5q, and 15q, whereas the Xq gain is specific to metastatic cells. B, UMAP of tumor nuclei in the primary tumor (top) and metastasis (bottom). Mean expression of marker genes for each differentiation state are displayed with a color code. The last UMAP on the right shows the assignment of each cell cluster to “mesenchymal,” “liver progenitor,” or “hepatocytic” cell type based on the expression of marker genes. C, Schematic representation of genetic and nongenetic evolution in patient #3981. LCA, last common ancestor; M, mesenchymal; LP, liver progenitor; H, hepatocytic; HB, hepatoblastoma.

Figure 4.

Single-nucleus RNA-seq reveals molecular plasticity along tumor progression in one patient. A, Virtual copy-number profiles discriminate tumor and normal cells. In agreement with WES data, all tumor cells display gains at 1q, 2q, 5q, and 15q, whereas the Xq gain is specific to metastatic cells. B, UMAP of tumor nuclei in the primary tumor (top) and metastasis (bottom). Mean expression of marker genes for each differentiation state are displayed with a color code. The last UMAP on the right shows the assignment of each cell cluster to “mesenchymal,” “liver progenitor,” or “hepatocytic” cell type based on the expression of marker genes. C, Schematic representation of genetic and nongenetic evolution in patient #3981. LCA, last common ancestor; M, mesenchymal; LP, liver progenitor; H, hepatocytic; HB, hepatoblastoma.

Close modal

Thus, dynamic switches between three differentiation states operate in hepatoblastoma cells, resulting in high intrapatient heterogeneity in space and time. This phenotypic plasticity is more frequent in tumors with 11p15.5 alterations and involves transcriptional and epigenetic reprogramming of TF modules. Molecular phenotypes match well with histologic cell types and display drastically different immune infiltrates.

Cisplatin Resistance Results from the Expansion of Progenitor Cell Clones Acquiring Massive Mutation Load

Mutational signature analysis of 65 pediatric liver cancer genomes identified four single-base substitution (SBS), two doublet-base substitution (DBS), and five insertion/deletion (indel) signatures, most of which matched signatures previously described in pan-cancer studies (refs. 17, 18; Fig. 5A; Supplementary Fig. S9A and S9B; Supplementary Table S7). The majority of signatures, including the clock-like signatures SBS1 and SBS5, were found ubiquitously in the tumors, while the signature SBS18, commonly found in neuroblastoma (17) and related to oxidative DNA damage (19), was identified in a subset of hepatoblastomas. Finally, we identified a massive load of mutations due to the signatures SBS35, DBS5, and ID3 in 20/66 primary hepatoblastomas resected after neoadjuvant chemotherapy (vs. 0/17 hepatoblastomas sampled before cisplatin treatment, P = 0.009). Consistently, SBS35 and DBS5 are known to reflect the diverse mutation types caused by cisplatin adducts on DNA (18, 20). In primary hepatoblastomas after neoadjuvant chemotherapy, cisplatin mutations were subclonal and almost exclusively found in the “liver progenitor” tumor subgroup (18/27 vs. 2/38; P = 1.2 × 10−7; Fig. 5B and C). Furthermore, in four tumors with spatial heterogeneity, cisplatin signature was restricted to the “liver progenitor” component even though the non-“progenitor” samples were exposed to the same chemotherapy regimen (Supplementary Fig. S7).

Figure 5.

Massive load of cisplatin-induced mutations in chemoresistant hepatoblastomas (HB). A, SBS signatures identified in PLCs. Each signature is displayed according to the 96-substitution classification defined by substitution type and sequence context immediately 5′ and 3′ to the mutated base. B, Unsupervised classification of 65 PLC genomes based on their mutational signature exposures. Clinical and molecular annotations are depicted below the dendrogram. Bar graphs indicate the proportion of the four SBS, two DBS, and five indel (ID) signatures in each sample. C, Number of subclonal (left) and clonal (right) mutations attributed to signature SBS35 (cisplatin) in hepatoblastoma samples according to sample type (pre-chemotherapy biopsy, primary tumor, or relapse/metastasis) and molecular group (hepatocytic or mesenchymal vs. liver progenitor). D, Phylogenetic trees reconstructed for eight patients with hepatoblastoma with primary tumors and relapses/metastases analyzed by WGS or WES. The time and molecular group of each sample are shown above the trees, together with chemotherapy treatment. Driver alterations are indicated. Branch lengths are proportional to the number of mutations acquired, with a color code indicating the contribution of each mutational signature.

Figure 5.

Massive load of cisplatin-induced mutations in chemoresistant hepatoblastomas (HB). A, SBS signatures identified in PLCs. Each signature is displayed according to the 96-substitution classification defined by substitution type and sequence context immediately 5′ and 3′ to the mutated base. B, Unsupervised classification of 65 PLC genomes based on their mutational signature exposures. Clinical and molecular annotations are depicted below the dendrogram. Bar graphs indicate the proportion of the four SBS, two DBS, and five indel (ID) signatures in each sample. C, Number of subclonal (left) and clonal (right) mutations attributed to signature SBS35 (cisplatin) in hepatoblastoma samples according to sample type (pre-chemotherapy biopsy, primary tumor, or relapse/metastasis) and molecular group (hepatocytic or mesenchymal vs. liver progenitor). D, Phylogenetic trees reconstructed for eight patients with hepatoblastoma with primary tumors and relapses/metastases analyzed by WGS or WES. The time and molecular group of each sample are shown above the trees, together with chemotherapy treatment. Driver alterations are indicated. Branch lengths are proportional to the number of mutations acquired, with a color code indicating the contribution of each mutational signature.

Close modal

Interestingly, all 16 hepatoblastoma relapses and metastases that developed after chemotherapy displayed a deluge of clonal cisplatin-induced mutations (Fig. 5C), leading to high mutation loads (median = 4.3 mutations/Mb), comparable to those of adult tumors and 15 times higher than those of primary liver tumors of the same age (Supplementary Fig. S9C–S9E). We analyzed in detail eight patients with chemoresistant hepatoblastomas at different stages of the disease including primary tumors, relapses, and metastases operated on between 4 months and 12 years after chemotherapy. All relapses were clonally related to their matched primary tumors, even the lung metastasis of patient #3538 (germline APC), detected 12 years after initial surgery (Fig. 5D). Phylogenetic trees revealed that a sudden burst of mutations is acquired during cisplatin treatment. Every relapse/metastasis was derived from a single common ancestor cell that acquired between 5,000 and 13,000 mutations during cisplatin treatment, with various modes of metastatic seeding. In patients #3370 and #3981, SBS35 mutations were private to each sample, indicating that metastatic seeding involved independent resistant cells. In patient #3694, the two relapses displayed a mixture of shared and private SBS35 mutations, indicating successive rounds of cisplatin mutagenesis followed by clonal expansion of resistant cells. Finally, patients #3529 and #3949 displayed late branching of relapse/metastasis samples derived from a same cisplatin-resistant clone. The seeding abilities of these clones are clearly illustrated by patient #3529 whose primary tumor was treated with liver transplantation: the same cisplatin-resistant clone gave rise first to a metastasis in the spleen, followed by two subsequent metastases on the grafted liver (Fig. 5D).

Overall, these data suggest that “liver progenitor” cells are able to bypass cisplatin–DNA adducts and proliferate under chemotherapy while accumulating SBS35 mutations. Accordingly, we found an enrichment of genes related to cisplatin resistance (21) or DNA repair among genes overexpressed in “liver progenitor” hepatoblastomas (63 genes, P < 2.2 × 10−16). Among those, 20 genes were also upregulated in post- versus prechemotherapy “liver progenitor” hepatoblastomas (Supplementary Fig. S10A–S10D), including genes involved in inhibition of apoptosis (BIRC5, coding for Survivin) and DNA repair through homologous recombination (BRCA1, RAD54L, and EXO1), Fanconi anemia (FANCA, FANCB, FANCD2, and FANCI) or base-excision repair pathways (LIG1, LIG3, POLE, and POLE2). Of note, most of these genes were also overexpressed in fetal liver samples (13th–30th weeks of amenorrhea) versus postnatal liver.

The subclonal presence of the signature SBS35 in primary tumors after neoadjuvant chemotherapy is thus a marker of cisplatin-resistant cell proliferation and was associated with poor progression-free and overall survival (log-rank P = 0.012 and P = 0.032, respectively; Supplementary Fig. S11). Heavily mutated resistant cells later give rise to relapses and metastases. A median of 70 coding sequence mutations per relapse/metastasis occurred due to the extra mutation load attributed to cisplatin. Although some of these mutations affected known cancer genes including NF1, BRAFL485W, KMT2C, KMT2D, BCORL1, and NOTCH1, no recurrent driver gene associated with tumor progression was identified.

New Therapeutic Strategies Targeting Hepatoblastoma Based on Molecular Features

To investigate new therapeutic options, we characterized a panel of nine pediatric liver cancer cell lines (PL-CCL; eight hepatoblastomas and one HCC; Fig. 6A; Supplementary Table S8). Our panel of cell lines, which are mainly derived from older patients and of “liver progenitor” type with multiple alterations in cancer driver genes, reproduced the major driver events identified in PLC. In particular, three cell lines with a matched primary or relapsed hepatoblastoma displayed genomic alterations globally similar to the original tumors (Supplementary Table S9). All hepatoblastoma cell lines carried CTNNB1 alteration whereas Hep3B, derived from an HCC, carried TP53, AXIN1, and RPS6KA3 alterations with frequent homozygous deletions in agreement with the deletor phenotype that we observed in pediatric HCC. Furthermore, all cell lines exhibited 11p15.5 locus alterations through cn-LOH (5/9), GOM IC1 (3/9), or CDKN1C mutation (1/9).

Figure 6.

New therapeutic strategies targeting hepatoblastoma (HB) based on molecular features. A, Heat map representation of somatic alterations in nine pediatric liver cancer cell lines (PL-CCL) and corresponding patient age at diagnosis (top). B, Correlation between sensitivity to cisplatin assessed with AUC and “liver progenitor” 7-gene expression signature. Color gradient intensity reflects levels of “liver progenitor” signature. Pearson correlation test was performed. C, Accumulation of SBS35 signature in four out of five PL-CCL. Trapezium indicates cell passaging and triangle symbolizes cell line limit dilution and clonal expansion. D, Correlation between sensitivity to neutralizing antibody anti-IGF2 (AUC) and IGF2 expression in five PL-CCL (Pearson test). E, Trametinib sensitivity (AUC) in nine PL-CCL (top) and corresponding genetic alterations in the RAS–MAPK pathway (bottom). F, Volcano plot representing differentially expressed genes in hepatoblastomas belonging to “hepatocytic” and “mesenchymal” and “Liver progenitor” transcriptomic subgroups. G, PLK1, BIRC5, and CHEK1 interactions and roles in cell cycle, DNA damage response, and apoptosis. H, Drug response assessed with AUC in nine PL-CCL. Color gradient intensity reflects levels of “liver progenitor” signature. Pearson correlations are performed between AUC and levels of “liver progenitor” signature. I, Sensitivity to cisplatin, doxorubicin, and BI-2536 in vivo in HepG2 xenografts from 24 nude mice and (J) mice weight follow-up. Mann–Whitney–Wilcoxon tests were performed at day 37.

Figure 6.

New therapeutic strategies targeting hepatoblastoma (HB) based on molecular features. A, Heat map representation of somatic alterations in nine pediatric liver cancer cell lines (PL-CCL) and corresponding patient age at diagnosis (top). B, Correlation between sensitivity to cisplatin assessed with AUC and “liver progenitor” 7-gene expression signature. Color gradient intensity reflects levels of “liver progenitor” signature. Pearson correlation test was performed. C, Accumulation of SBS35 signature in four out of five PL-CCL. Trapezium indicates cell passaging and triangle symbolizes cell line limit dilution and clonal expansion. D, Correlation between sensitivity to neutralizing antibody anti-IGF2 (AUC) and IGF2 expression in five PL-CCL (Pearson test). E, Trametinib sensitivity (AUC) in nine PL-CCL (top) and corresponding genetic alterations in the RAS–MAPK pathway (bottom). F, Volcano plot representing differentially expressed genes in hepatoblastomas belonging to “hepatocytic” and “mesenchymal” and “Liver progenitor” transcriptomic subgroups. G, PLK1, BIRC5, and CHEK1 interactions and roles in cell cycle, DNA damage response, and apoptosis. H, Drug response assessed with AUC in nine PL-CCL. Color gradient intensity reflects levels of “liver progenitor” signature. Pearson correlations are performed between AUC and levels of “liver progenitor” signature. I, Sensitivity to cisplatin, doxorubicin, and BI-2536 in vivo in HepG2 xenografts from 24 nude mice and (J) mice weight follow-up. Mann–Whitney–Wilcoxon tests were performed at day 37.

Close modal

We tested our panel of cell lines with conventional therapeutic regimens (cisplatin, carboplatin, and doxorubicin treatments; Supplementary Table S10). In the nine cell lines analyzed, resistance to cisplatin and carboplatin correlated with a high expression of “liver progenitor” markers (Fig. 6B; Supplementary Fig. S12A). Additionally, we treated five cell lines with long-term exposure to low doses of cisplatin (0.5 μmol/L, 4 weeks; Supplementary Fig. S12B). Only B6-2 was completely sensitive after two weeks. In the four remaining cell lines, WES of resistant clones revealed an accumulation of the SBS35 cisplatin mutational signature. These results demonstrate that the SBS35 signature is directly induced by cisplatin exposure in resistant cells (Fig. 6C). Interestingly, doxorubicin, which is commonly used in resistant or high-risk patients with hepatoblastoma, exhibited a significant antitumor effect in all cell lines.

Next, we tested treatments targeting specific genomic alterations and their consequences. Three out of five tested cell lines were sensitive to anti-IGF2 antibodies, and sensitivity to this drug correlated with the level of IGF2 expression (P = 0.03, R = −0.91, Pearson correlation; Fig. 6D). We also tested the monoclonal therapeutic antibody xentuzumab, and we obtained between 14% and 47% growth inhibition (Supplementary Fig. S12C), consistent with a previous report in mice xenografts (22). Surprisingly, sensitivity to xentuzumab was not correlated with IGF2 expression, possibly because it neutralizes both IGF1 and IGF2. We then tested our cell lines with the MEK1/2 inhibitor trametinib (Fig. 6E). Remarkably, five of nine cell lines were sensitive to trametinib and harbored alterations within the MAP kinase pathway through RPS6KA3 homozygous deletions, NRAS-activating mutations, or FGF19 amplifications as previously shown in adult HCC cell lines (23). Of note, one of our patients with a metastatic hepatoblastoma (#3538) displayed a hotspotMAPK1E322K trunk mutation and acquired a BRAFL485W mutation in the metastasis that has previously been associated with response to an ERK inhibitor (24).

To identify new targets in cisplatin-resistant “progenitor” cells, we performed a differential expression analysis between the “liver progenitor” and other hepatoblastoma subgroups. Of the genes overexpressed in “liver progenitor” samples, we selected three key targetable genes: PLK1, BIRC5 (Survivin), and CHEK1, which are involved in cell-cycle, apoptosis, and DNA damage response, respectively (Fig. 6F and G). In patients with hepatoblastoma, high expression of PLK1, BIRC5, and CHEK1 was associated with (i) poor response to chemotherapy as assessed with AFP reduction in the sera (P = 4.1 × 10−4, 2.0 × 10−4, and 4.7 × 10−4, respectively), (ii) embryonal histology (P = 1.2 × 10−5, 5.1 × 10−6, and 7.1 × 10−6), (iii) the cisplatin-resistance signature SBS35 (P = 8.8 × 10−10, 1.3 × 10−9, 7.6 × 10−9), and (iv) a high expression of proliferation genes (P ≤ 2.2 × 10−16; Supplementary Fig. S13A–S13D). We thus tested inhibitors of PLK1 (BI-2536), Survivin (YM155), and CHEK1 (AZD-7762) in our collection of nine cell lines. Interestingly, the three drugs were more efficient than both cisplatin and carboplatin in all the cell lines except B6-2, which has the least “liver progenitor”-like phenotype. All cell lines were sensitive to BI-2536, including those that are resistant to cisplatin (Fig. 6H). Notably, the cell lines most resistant to cisplatin were also the most sensitive to YM155 and AZD-7762 (Supplementary Fig. S14A–S14C). Finally, we tested the most effective drug, BI-2536, in vivo in 24 nude mice xenografted with HepG2 cells randomized in four groups (cisplatin, cisplatin + doxorubicin, BI-2536, and control). HepG2 xenografts were resistant to cisplatin treatment, both in isolation and in combination with doxorubicin. In contrast, tumors in mice receiving BI-2536 responded to treatment (P = 0.003) without significant toxicity, unlike with cisplatin (Fig. 6I and J).

In conclusion, the PLK1 inhibitor BI-2536 showed efficacy in both in vitro and in vivo models and appears to be a good candidate drug for the treatment of PLC. Trametinib and anti-IGF2 antibodies show promise for the treatment of resistant hepatoblastomas with CCND1/FGF19 amplification or high IGF2 expression, respectively.

Altogether, our findings suggest a major role of cell plasticity in the spatiotemporal evolution of hepatoblastomas and their resistance to cisplatin-based chemotherapy. In this study, genomic analysis of multiple parts of primary tumors and of longitudinal samples at various time points allowed us to profile the clonal evolution of hepatoblastoma along treatment. Cisplatin induced hundreds of mutations in progenitor tumor compartments, a phenomenon reproduced in cell lines, and targeting specific oncogenic processes with different inhibitors could be used to target resistant cells.

Primary liver cancers showed a low mutation burden, as seen in other pediatric cancers (8). In PLC, β-catenin activation, mainly caused by CTNNB1 mutations, acted as the major early and common mechanism of liver tumorigenesis. However, we identified an earliest recurrent premalignant alteration with mosaic cn-LOH of 11p15 in the liver of very young patients with hepatoblastoma (10%). This finding is in line with similar premalignant clonal expansions recently identified as precursors of Wilms tumors, associated with hypermethylation of IC1 at the 11p15 locus (25). Our data reinforce the role of 11p15 alterations in premalignant stages that occur during pregnancy or early in life. Systematic searches of 11p15 locus alterations in the liver and other organs of patients with pediatric cancer could reveal additional mosaic cases.

The present molecular transcriptomic classification revealed four robust hepatoblastoma subgroups defined by both their differentiation states and immune infiltration levels. We verified that our clustering was not driven by patients with multiple samples (Supplementary Fig. S15), and we validated our classification in an independent RNA-seq cohort of 34 hepatoblastomas (Supplementary Fig. S6). Our transcriptomic classification overlaps the previously published subtypes C1/C2 (26), C1/C2A/C2B (27), and MRS1/MRS2/MRS3a/MRS3b (16), with a good match between “hepatocytic” and C1, “liver progenitor” and C2A, and “mesenchymal” and C2B subgroups (Fig. 3A; Supplementary Fig. S16). Our “liver progenitor” subgroup also overlaps with the very high-risk MRS-3b subgroup identified by Carrillo-Reixach and colleagues (ref. 16; P = 2.0 × 10−6). The most original findings of our classification are (i) the definition of clear transcriptomic signature of “mesenchymal” histologic cell type validated by careful histologic review of mirror blocks, and (ii) the identification of two subgroups of “hepatocytic” hepatoblastomas showing varying levels of immune infiltration, which is potentially enhanced by cisplatin-based chemotherapy. Hepatoblastoma subgroups are not defined by specific driver events, as opposed to medulloblastoma subgroups (28), for example. In turn, each state of hepatoblastoma differentiation is defined by specific TFs and a characteristic DNA methylation landscape. Our methylation-based classification revealed two main hepatoblastoma clusters reminiscent of the previously described Epi-CA/CB groups (16): a normal-like cluster matching the Epi-CA group comprising most “hepatocytic” samples, and a differentially methylated cluster matching the Epi-CB group, which we further divided into three subgroups associated with “liver progenitor,” “mesenchymal” phenotype, and older age. Importantly, hepatoblastoma cells are able to switch between differentiation states along the course of the disease, leading to striking spatial and temporal heterogeneity. Our single-cell data support the notion that hepatoblastoma cells have the inherent ability to switch across three differentiation states. Yet, clinical and molecular features associated with transcriptomic subgroups suggest that the cell of origin, patient age, chemotherapy, and driver alterations can favor the predominance of a given subtype in each tumor (Supplementary Table S11; Supplementary Fig. S17). This molecular plasticity, driven by specific TF modules and epigenetic landscapes, is reminiscent of neuroblastoma cells that can transdifferentiate between committed adrenergic cells and undifferentiated mesenchymal cells (29, 30), the latter being chemoresistant and enriched in relapse tumors. Similarly, we showed that hepatoblastoma “progenitor” cells proliferate under neoadjuvant chemotherapy and are enriched in relapses and metastases.

Moreover, PLC development follows evolutionary trajectories that vary according to age, involving similar pathways but diverse mechanisms of activation. In contrast to hepatoblastomas, pediatric HCCs occur in older children, on cirrhotic liver due to various causes. Like hepatoblastomas, HCC displayed frequent 11p15.5 alterations, but also a characteristic chromosome deletor phenotype leading to WNT/β-catenin activation through AXIN1 or AMER1 deletions, as well as frequent GPC3, RPS6KA3, and SMARCA1 alterations. Because a deletor phenotype is seen in only 3% of adult HCC (31), this oncogenic mechanism could be associated with young age. Interestingly, many HCC driver genes are located on the X chromosome (AMER1, RPS6KA3, GPC3, and SMARCA1), which may partly explain the enrichment of pediatric HCC in male patients described previously (9) and retrieved in the present cohort (7/9 males, 78%). Pediatric HCCs are commonly resistant to cisplatin, but their mechanism of resistance likely differs from that involved in hepatoblastomas because they do not exhibit similar progenitor and stem cell features.

“Progenitor” cells are able to proliferate under treatment, bypass cisplatin-induced adduct formation, acquire a huge number of cisplatin-induced mutations (median 11,142) with a characteristic SBS35 signature, and give rise to highly mutated relapse tumors. Various DNA polymerases and DNA repair genes involved in double-strand break repair, homologous recombination, and Fanconi pathway are already overexpressed both in “liver progenitor” hepatoblastomas before chemotherapy and in the fetal liver, showing that these mechanisms are constitutively active in liver progenitor cells. Relatedly, a previous work had shown that Fanconi anemia inhibitors were able to block the growth of hepatoblastoma cells in vitro and in vivo (27). Pich and colleagues estimated that ∼15 of every 1,000 cisplatin-induced mutations affect the sequence of coding genes, of which ∼0.7 are expected to affect the sequence of known cancer genes (32). Thus, an average hepatoblastoma relapse or metastasis is at risk of acquiring 167 coding mutations, including 7.8 in cancer genes due to cisplatin treatment. Although most cancer gene mutations are unlikely to play a functional role in hepatoblastoma cells, this increased mutational burden provides an opportunity for resistant cells to acquire additional oncogenic capabilities. These data highlight the necessity to identify alternative treatments for cisplatin-resistant hepatoblastomas, in which more aggressive chemotherapy regimens may just promote the selection of aggressive progenitor cells with a massive extra mutational burden.

Our findings revealed new potential therapies to combat hepatoblastoma resistant to cisplatin. Targeting the different mechanisms involved in cisplatin resistance (21) such as apoptosis (YM155 targeting Survivin), DNA repair (AZD-7762 targeting CHEK1), or cell-cycle control (BI-2536 as a PLK1 inhibitor) is very promising because an efficient antitumor effect was observed in PLC cell lines with a progenitor molecular signature. In the same line, other drugs targeting the proteasome (bortezomib) could also be efficient in resistant hepatoblastomas (27). Immunotherapy could be another appealing type of treatment for hepatoblastomas because we showed that cisplatin treatment can induce an intratumor polymorphic immune response leading to “hepatocytic hot” hepatoblastomas, consistent with our previous observation in hepatoblastomas with germline APC mutations (33). Reinforcement of an efficient intratumor immune response could be improved by the use of oxaliplatin to induce an immune cell death (34) or by treatments with immunomodulators. However, “liver progenitor” hepatoblastomas are immune- cold despite showing multiple cisplatin-induced mutations. In these cases, penetrating the progenitor compartments with immune cells will be challenging. Ongoing and planned clinical trials aim to test immunomodulators or targeted therapies in high-risk hepatoblastoma; however, these trials are limited by the small number of patients.

In conclusion, PLC showed various mechanisms of tumorigenesis related to age and cell of origin. Hepatoblastomas demonstrated a striking spatial and longitudinal phenotypic plasticity related to the progenitor compartment associated with cisplatin resistance and the mutational signature of DNA adduct bypass. And, finally, the identification of drugs targeting “progenitor” cells opens new avenues to treat children at high risk of resistance.

Clinical Samples

A series of 126 liver tumor samples and their nontumor counterparts were collected from 86 patients surgically treated in various French hospitals. The study was approved by the local Ethics Committee (CCPRB Paris Saint-Louis). Written informed consent was obtained in accordance with French legislation. All samples were immediately frozen in liquid nitrogen and stored at −80°C. Tumors included in this study comprised 87 primary hepatoblastoma (18 prior chemotherapy and 69 after chemotherapy) and 17 hepatoblastoma relapses and metastases (from 65 patients with hepatoblastoma including 14 older than five years old), five HCAs (five patients), seven FLCs (seven patients), and ten HCCs (nine patients). HCCs were developed in various etiologic contexts including tyrosinemia (two patients), mitochondrial cytopathy (three patients), and progressive familial intrahepatic cholestasis (three patients). Samples were analyzed by a combination of WGS (n = 65), WES (n = 57), RNA-seq (n = 120) and RRBS (n = 92). A summary of the cohort is provided in Supplementary Fig. S1, and detailed clinical characteristics of each sample are provided in Supplementary Table S1. Also, seven human fetal liver samples between the 13th and 30th weeks of amenorrhea were analyzed by RNA-seq.

Pathologic Reviewing

All tumors were reviewed by three expert pathologists specializing in pediatric liver tumors. For hepatoblastoma, fractions of histologic components (fetal, embryonal, mesenchymal, cholangioblastic, and small cell undifferentiated) were estimated according to the consensus classification (14) for the whole tumor as well as for mirror blocks corresponding to frozen samples when available. Thus, spatial heterogeneity was defined by the coexistence of at least two histologic components, either at the whole tumor level or at the intrasample level (mirror block).

WGS

We extracted DNA using the Maxwell DNA extraction kit (Promega) or the AllPrep DNA/RNA/miRNA Universal Kit (Qiagen). Sixty-five tumors and matched nontumor samples were sequenced at the Centre National de Génotypage (Evry, France) on an Illumina HiSeqX5 as paired-end 150 bp reads, with an average depth of 90× for tumors and 30× for nontumor liver samples. Sequence reads were aligned to the hg19 version of the human genome using BWA (35). We used Picard tools (http://broadinstitute.github.io/picard/) to remove PCR duplicates and GATK (36) for local indel realignment and base quality recalibration, as recommended in GATK best practices (37).

WES

Whole-exome data from 57 tumors and matched nontumor samples were analyzed in this study, as well as nine pediatric liver tumor cell lines. Sequence capture, enrichment, and elution of genomic DNA were performed by IntegraGen (Evry). Agilent in-solution enrichment was used with the manufacturer's biotinylated oligonucleotide probe library SureSelect Clinical Research Exome V2 or Twist Bioscience Human Core Exome Enrichment System, according to the manufacturer's instructions. The eluted enriched DNA samples were sequenced on an Illumina HiSeq4000 as paired-end 75 bp reads (n = 35) or Illumina NovaSeq as paired-end 100 bp reads (n = 31), with an average depth of 100× for tumors and 65× for nontumor liver samples. We used BWA to align reads on the hg38 version of the human genome and sambamba to remove duplicate reads.

Somatic Mutation Calling

We used MuTect2 to call somatic mutations from WES and WGS data by comparing each tumor sample with its matched nontumor counterpart. The three cell lines derived from patient tumors within our cohort were compared with the matched nontumor samples, while the six other cell lines were compared with a panel of normal samples. We excluded mutations belonging to the ENCODE Data Analysis Consortium blacklisted regions (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeDacMapabilityConsensusExcludable.bed.gz) and regions covered by < 6 reads in the tumor or normal sample. We then selected only single-nucleotide variants (SNV) with a MuTect2 flag among “PASS,” “clustered_events,” “t_lod_fstar,” “alt_allele_in_normal,” or “homologous_mapping_event” and small insertions and deletions (indels) with a MuTect2 flag among “PASS,” “clustered_events,” or “str_contraction.” To improve specificity in the calling of mutations with low variant allele frequency (VAF), we quantified the number of high-quality variant reads in the tumor (mapping quality ≥ 20, base quality ≥ 20) and the number of variant reads in the nontumor sample with no quality threshold using bamreadcount (https://github.com/genome/bam-readcount). Only variants matching the following criteria were finally retained: VAF ≥ 2% in the tumor with ≥ 3 variant reads, VAF ≤ 5% in the nontumor samples with ≤ 2 variant reads, and a VAF ratio ≥ 5 between the tumor and nontumor sample.

Germline Mutation Calling

Germline variant calling was performed independently for WGS and WES data according to GATK (version 4.0) best practices recommendations (38). For each data set, we performed SNV and indel discovery using HaplotypeCaller and joint genotyping across all non-tumor samples simultaneously using GenomicsDBImport followed by GenotypeGVCFs. We used hard filtering for WES data with different parameters for SNVs (ExcessHet > 54.69, FS > 10.0, MQ < 50.0, MQRankSum < −5.0, QD < 5.0, QUAL < 50.0, ReadPosRankSum < −5.0, SOR > 2.0, DP < 8.0, GQ < 20.0) and indels (ExcessHet > 54.69, FS > 200.0, QD < 2.0, QUAL < 30.0, ReadPosRankSum < −20.0, DP < 8.0, GQ < 20.0), and we used variant quality score recalibration for WGS (truth sensitivity level: SNVs = 99.6%, indels = 95.0%). Resulting high-quality variants were annotated using the Variant Effect Predictor toolset (39). We selected variants with an allele frequency lower than 0.01 or lacking in the gnomAD database (40). All candidate pathogenic variants were manually verified using the Integrative Genomics Viewer (41).

Copy-Number and Structural Rearrangement Analysis

We used MANTA (42) software to identify somatic structural rearrangements in WGS data. To keep only the most reliable events, we selected only rearrangements supported by ≥ 15 reads and with a variant allele fraction ≥ 10% in the tumor, and not more than 1 variant read in the nontumor counterpart. We used the cgpBattenberg (43) algorithm to reconstruct absolute copy-number profiles from WGS data and estimate tumor purity. We used the Genome Alteration Print method (44) to infer absolute copy-number profiles and tumor purity from WES data, and the circular binary segmentation algorithm implemented in the Bioconductor package DNAcopy (45) to identify focal homozygous deletions and high-level amplifications, as previously described (46).

Identification of Driver Genes

We used two different approaches to identify driver genes in PLCs, using both WES and WGS data. First, the MutSigCV tool (47) was used to identify genes with significantly recurrent mutations while taking into account gene size and genomic covariates. High-level amplifications and homozygous deletions were included in this analysis as additional mutation categories. Second, we used Oncodrive (48) to identify genes with a significant enrichment of mutations with functional impact. Finally, we defined as putative drivers genes with damaging alterations in ≥2 patients and either (i) a P value < 0.05 with both MutSigCV and Oncodrive tools or (ii) identified as driver genes in previous pediatric pan-cancer studies (25, 49). After removing three genes not expressed in normal liver nor in tumors (95th FPKM quantile < 0.1), we obtained 22 candidate driver genes (Supplementary Table S2). Only genes with damaging alterations in ≥3 patients were represented in Fig. 1B. The frequency of driver alterations per tumor type was estimated by patient and not by sample, in order to avoid biases due to patients with multiple samples. To this end, a patient was considered as altered for a gene as long as one of its samples harbored an alteration in this gene.

Identification of Mosaic Copy-Neutral LOH of 11p15.5 Locus in Nontumor Liver

For each tumor with a cn-LOH of 11p15.5 locus, we searched for the presence of the same cn-LOH in matched nontumor liver tissue. To that aim, we calculated the BAF of common single-nucleotide polymorphisms (SNP), obtained from WGS or WES data, in the tumor (BAFT) and in the nontumor (BAFNT) sample. We considered a cn-LOH to be present in the nontumor liver if there was a significant overlap between SNPs with the B allele retained in the tumor (BAFT > 0.5) and those with a BAFNT > 0.5 (binomial test). We then estimated the proportion of nontumor liver cells carrying the cn-LOH as 2*RAFNT-1 with RAFNT the median BAF of retained alleles (those with a BAFT > 0.5) in the nontumor sample.

Mutational Signature Analysis

We used Palimpsest (50) to extract signatures of SBS, DBS, and indels in WGS data of 65 PLCs, and to compare them with known signatures from the COSMIC database (v3). De novo analysis revealed four SBS signatures corresponding to the COSMIC SBS1, 5, 18, and 35 (cosine similarity scores > 0.8), two DBS signatures corresponding to the COSMIC DBS5 and a new signature (DBSnew), and five indel signatures corresponding to the COSMIC ID1, 2, 3, 5, and 8. To ensure comparability with other data sets, signature exposures in each tumor were recalculated using the COSMIC versions of signatures, except for signature SBS35 (we kept our own version as it was based on a larger number of mutations than the COSMIC version) and DBSnew (absent in COSMIC). SBS signature exposures were also calculated in WES data. Finally, we estimated the cancer cell fraction (CCF), i.e., the proportion of tumor cells carrying each mutation (see the “Clonality” section below). Mutations with an upper bound of their CCF confidence interval < 0.95 were considered subclonal. We then quantified SBS signature exposures among clonal and subclonal mutations separately.

Clonality Analysis and Phylogenetic Tree Reconstruction

We used Palimpsest to estimate the CCF of each mutation from its variant allele fraction, taking into account tumor purity and absolute copy number, as previously described (31). For 23 patients with multiple tumor samples, we used a Bayesian Dirichlet process in multiple dimensions (51) to identify clusters of mutations with a similar CCF distribution across all samples, hence belonging to the same branch of the phylogenetic tree. We then manually reconstructed the phylogenetic tree of each patient by organizing branches to fit the observed clonal composition of each sample. Damaging mutations affecting candidate driver genes were annotated on the trees, and we quantified the contribution of each mutational process on each branch with Palimpsest.

RNA-seq

We performed RNA-seq for 120 tumors, 4 nontumor liver samples, 7 fetal liver samples, and 9 pediatric liver tumor cell lines. RNA samples were enriched for polyadenylated RNA from 1 μg of total RNA, and the enriched samples were used to generate sequencing libraries with the Illumina TruSeq Stranded mRNA Kit or the NEBNext Ultra II Directional RNA Library Prep Kit and associated protocol as provided by the manufacturer. Libraries were sequenced by IntegraGen (Evry) on an Illumina HiSeq 4000 as paired-end 75 bp reads or Illumina NovaSeq as paired-end 100 bp reads. Full Fastq files were aligned to the reference human genome hg38 using TopHat2 (52). We removed reads mapping to multiple locations, and we used HTSeq (53) to obtain the number of reads associated with each gene in the Gencode database. We used DESeq2 (54) to import raw read counts into R statistical software and apply variance stabilizing transformation to the raw count matrix. FPKM scores (number of fragments per kilobase of exon model and millions of mapped reads) were calculated by normalizing the count matrix for the library size and the coding length of each gene.

Gene Fusion Detection

Fusions detected by TopHat2 (ref. 52; –fusion-search –fusion-min-dist 2000 –fusion-anchor-length 13 –fusion-ignore-chromosomes chrM) were filtered using the TopHatFusion-post algorithm and validated using FusionInspector (ref. 55; https://github.com/FusionInspector). We kept only fusions validated by BLAST, with at least 10 split-reads or read pairs spanning the fusion event, and an FFPM (fusion fragments per million reads) ≥ 0.1. We removed fusions identified recurrently in a cohort of 36 normal liver samples, involving genes with inconsistent orientations, noncoding genes, and putative read-through transcripts.

Detection of Large In-Frame CTNNB1 Deletions

Large in-frame CTNNB1 deletions were screened in WGS, WES, and RNA-seq data using dedicated approaches. We used MANTA software (42) to identify structural variations spanning CTNNB1 region in WGS and WES (with the –exome option) data, with the following parameters to increase sensitivity: minEdgeObservations = 2, minScoredVariantSize = 10. For RNA-seq, we analyzed the junctions.bed output files from TopHat2 (52) to identify abnormal junctions involving CTNNB1 exons 3 or 4. We also screened for paired reads with one mate in exon 2 and the other in exon 3/4, and an abnormally long insert size, greater than the mean insert size of the library + 1 standard deviation.

Gene Expression–Based Classification of Hepatoblastoma

We restricted the transcriptomic classification to hepatoblastoma samples, excluding other diagnoses (HCC, FLC, and HCA) with a limited number of samples, in order to identify robust subgroups. Hierarchical clustering was performed on 100 hepatoblastomas and 4 nontumor liver samples, based on the vst-normalized expression of the 3,000 autosomal genes with sufficient expression (95th FPKM percentile ≥ 0.1) and the highest standard deviation. The hclust R function was used with Euclidean distance and Ward.D2 linkage method, after median-centering the data. To verify that hepatoblastoma subgroups were not driven by patients with multiple samples, we reproduced the clustering with only one sample by patient and we obtained similar results (Supplementary Fig. S15). For external validation of hepatoblastoma clusters, we downloaded an RNA-seq data set from Carrillo-Reixach and colleagues (16) comprising 34 hepatoblastomas (Gene Expression Omnibus accession number GSE132219) and performed a hierarchical clustering based on the expression of the 19 markers of hepatic differentiation, progenitor, and mesenchymal subgroups represented in Fig. 3, with Cosine distance and Ward.D2 linkage method (Supplementary Fig. S6). We used the Bioconductor limma (56) package to identify differentially expressed genes between hepatoblastoma subgroups. We applied a q-value threshold of ≤ 0.05 to define differentially expressed genes.

Gene-Expression Signatures

We used previously established molecular signatures from the MSigDB database to quantify the level of hepatic differentiation (“Hsiao liver-specific genes”), cell proliferation (“Hallmark REACTOME cell cycle”), and inflammation (“Hallmark inflammatory response”) in each sample. We used MCPcounter (57) to estimate the infiltration by diverse immune cell types and classify hepatoblastoma samples accordingly. We also classified each hepatoblastoma sample according to the previously described C1/C2 (27), C1/C2A/C2B (28), and RMS-1/2/3a/3b (16) classifications. For C1/C2 (respectively C1/C2A/C2B) classification, we generated a consensus clustering (1,000 resampling iterations of hierarchical clustering with Euclidean distance and Ward.D2 linkage) of our cohort in two (respectively three) groups based on the 16 marker genes defined by Cairo and colleagues (ref. 27; respectively the four markers VIM, HSD17B6, TOP2A, and ITGA6 defined by Hooks and colleagues; ref. 28), and we assigned each consensus cluster to the relevant subgroup based on the expression of their respective markers. For the RMS classification, we first performed a hierarchical clustering based on the expression of 68 genes of the 14q32 imprinted locus to define two groups (moderate vs. strong 14q32-gene signature). We verified that our two main DNA methylation-based clusters were consistent with Carrillo-Reixach's Epi-CA/CB groups (16). Finally, we combined the 14q32 signature, Epi-CA/CB methylation groups, and C2A transcriptomic subgroup to define the RMS classification as in Carrillo-Reixach and colleagues (16). A signature of the “Liver progenitor” subgroup of hepatoblastomas comprising the seven representative genes indicated in Fig. 3A was used to quantify this phenotype in cell lines.

Single-Nucleus Isolation from Frozen Tissues

Single nuclei were isolated from four samples for snRNA-seq as previously described (58), using EZ Lysis buffer workflow with slight modifications. Briefly, tissue samples were thawed in PBS and cut into pieces < 0.5 cm. Approximately 35 mg of tissue was poured in a glass Dounce tissue grinder (Sigma, cat. no. D8938) and homogenized 25 times with pestle A and 25 times with pestle B in 1.5 mL of ice-cold nuclei EZ lysis buffer. Samples were then incubated on ice for 5 minutes with an additional 3 mL of cold EZ lysis buffer. Nuclei were centrifuged at 500 × g for 5 minutes at 4°C, washed with 5 mL ice-cold EZ lysis buffer, and incubated on ice for 5 minutes. After centrifugation, the nucleus pellet was washed with 5 mL of Nuclei Wash buffer containing 1 × PBS, 0.1%, nonacetylated BSA (Thermo AM2618), and 200 units/mL RNase inhibitor (NEB M0307L). Isolated nuclei were resuspended in 2 mL of Nuclei Suspension Buffer containing 1 × PBS, 1% nonacetylated BSA (Thermo AM2618), and 200 units/mL RNase inhibitor (NEB M0307L), filtered through a 70-μm and then a 30-μm MACS SmartStrainer (Miltenyibiotec 130-098-462 and 130-098-458), and counted under microscope using C-chip disposable hemocytometer. A final concentration of 1,000 nuclei per μL was used for loading on a 10× channel.

snRNA-seq

snRNA-seq was performed by Integragen SA (Evry) on matched nontumor liver, primary hepatoblastoma, and lung metastasis samples from one patient (#3981), following the Chromium Next GEL Single Cell 3′ V3.1 protocol. In short, about 8,800 single nuclei were loaded into each channel of a Chromium single-cell 3′ chip. Single nuclei were partitioned into droplets with gel beads in the Chromium Controller. After emulsions were formed, barcoded reverse transcription of RNA took place, followed by cDNA amplification, fragmentation, and adapter and sample index ligation, according to the manufacturer's recommendations. Libraries from the 10X channels were pooled together and sequenced as paired-end 100 bp reads on an Illumina NovaSeq. We used 10X Genomics Cell Ranger 5.0 (59) to align snRNA-seq reads to the human genome (GrCh38/hg38) and generate UMI counts for each sample, including intronic reads. We obtained, respectively, 9,142, 5,859, and 2,875 nuclei, with a median of 2,042, 2,334, and 3,446 genes per nuclei from samples #3981N (nontumor liver), #3984T (primary hepatoblastoma), and #3988T (lung metastasis). We then filtered the feature-barcode matrix to retain only good-quality nuclei and reliable genes. We removed nuclei with <1,000 read counts, <500 genes detected, or >5% of UMI counts in mitochondrial genes, and we removed genes detected in <3 nuclei as well as ERCC and mitochondrial genes. After QC we kept, respectively, 9,140, 5,685, and 2,826 nuclei, with a total of 16,470, 18,384, and 18,294 genes for samples #3981N, #3984T, and #3988T. All secondary analyses were performed using Seurat v3 (60). We normalized each data set using the SCTransform function with default parameters, performed principal component analysis on the 3,000 most variable genes and ran Louvain graph-based clustering on the 30 principal components with a resolution of 0.5. We used a Uniform Manifold Approximation and Projection (UMAP) with default settings to visualize the results. We used inferCNV v1.6 (Tickle T and colleagues, available from https://github.com/broadinstitute/inferCNV) to reconstruct virtual copy-number profiles, using healthy hepatocytes from the nontumor liver sample as reference, and genes with an average read count >0.1 in reference nuclei. We performed a clustering of copy-number profiles in each sample separately. We identified tumor cells using both inferCNV clusters and Seurat classification, and we redid the whole analysis (SCTransform normalization, dimensionality reduction, clustering and UMAP visualization) restricted to tumor cells. To characterize tumor cell clusters, we selected marker genes of each differentiation state (“hepatocytic,” “liver progenitor,” and “mesenchymal”). To that end, we selected the top 70 genes overexpressed in the corresponding group in bulk RNA-seq data, restricted to genes identified in snRNA-seq data. We obtained ∼50 marker genes for each subgroup, and we computed for each nucleus the mean log-normalized expression (NormalizeData pipeline) over each set of markers. We then assigned a status to each cluster based on the expression of these marker genes.

RRBS

RRBS was performed in two distinct projects (PJ17, n = 82 and PJ20, n = 23) by Integragen SA (Evry) as described by Gu and colleagues (61), with the Diagenode Premium RRBS kit. In brief, 100 ng of qualified genomic DNA were digested with MspI. After end-repair, A-tailing and ligation to methylated and indexed adapters, the size-selected library fragments were subjected to bisulfite conversion and PCR-amplified. Samples of the PJ17 project (respectively PJ20) were then sequenced on an Illumina HiSeq4000 (respectively NovaSeq) sequencer as paired-end 75 bp (respectively 100 bp) reads. Image analysis and base calling were performed using Illumina Real-Time Analysis with default parameters. Reads were aligned to the hg38 version of the human genome using BS_Seeker2 (62). Sorted bam files were converted into CGmap files using CGmaptools (63), and the methylation level of each CpG site was defined as the ratio between the number of effective cytosines after bisulfite conversion (= methylated cytosines) and the total number of cytosines and thymines after bisulfite conversion (= methylated + unmethylated cytosines). On average, approximately 7 million CpG sites were covered in each sample after discarding CpG sites located in ENCODE-blacklisted genomic regions (wgEncodeDacMapabilityConsensusExcludable tract from the UCSC genome browser). We next integrated methylation levels across 100-bp-long genomic regions (tiles) using the tileMethylCounts function from the Methylkit package (64). Tile coverage was defined as the coverage sum of all CpGs inside the tile, and methylation was defined as the ratio between the total number of methylated CpGs and tile coverage. On average, approximately 1.3 million tiles were covered in each sample. Due to the different read lengths in PJ17 and PJ20, some tiles displayed heterogeneous coverage leading to systematic biases between the two projects. We thus compared the nontumor liver samples from PJ20 (n = 4) and PJ17 (n = 9), and we removed 58,179 tiles (∼3%) with a methylation difference ≥ 0.05 or ≤ −0.05 between all PJ17 and all PJ20 nontumor samples.

To determine the methylation status of 11p15.5 locus, we computed for each sample the mean methylation over imprinted regions IC1 (chr11:1998745–2003509, hg38) and IC2 (chr11:2697587–2700983, hg38). We then used a K-means clustering to identify samples with GOM of IC1 and/or LOMof IC2.

DNA Methylation Changes in Hepatoblastoma

We generated a DNA methylation–based classification of 84 hepatoblastomas and 13 nontumor liver samples. To do so, we selected the 5,000 tiles with coverage >50× in all samples and the highest standard deviation, and we used the hclust R function to perform a hierarchical clustering with Pearson distance and Ward.D linkage method. We used Independent Component analysis as previously described (65) to characterize independent sources of DNA methylation changes in hepatoblastomas, based on ∼212,000 tiles with coverage >10× in all samples and standard deviation >0.035. We identified five methylation components (MC) including components related to the “liver progenitor” (MC1) and “mesenchymal” (MC2) subgroups, and to patient age (MC3). We explored various (epi)genomic features associated with differentially methylated tiles and MCs, including gene features (promoter/gene body), CpG island features (island/shore/shelf/open sea), chromatin states (66), and DNA methylation domains (ref. 67; highly methylated domains vs. partially methylated domains) in normal liver, and replication timing in the HepG2 cell line (68). We used ELMER v2 package (69) to identify transcription factor binding sites enriched within tiles hypermethylated in the “mesenchymal” component. We used an in-house adaptation of the gene set enrichment analysis method to identify overrepresented gene sets (MSigDB v6 database) among genes paired with tiles hypermethylated in “mesenchymal” samples. We compared the vst-normalized expression of TFs involved in liver differentiation (HNF1A, HNF4A, PPARG, and NR2F6) with the methylation of tiles containing their respective binding motifs and hypermethylated in the “mesenchymal” component. These correlations were validated in Carrillo-Reixach's data set comprising 26 hepatoblastomas analyzed with Illumina 850k methylation arrays (16).

Cell Lines

Nine pediatric liver tumor cell lines were collected from collaborations or obtained from commercial sources (Supplementary Table S8). Cells were grown in either DMEM or Advanced DMEM F-12 supplemented with 10% fetal bovine serum and using the usual conditions (100 U/mL penicillin/streptomycin, 1% glutamine at 37°C, 5% CO2, identity confirmed using WES, Mycoplasma-free verified with the MycoAlert Mycoplasma PLUS detection kit; Lonza).

Determination of Drug Sensitivity

Doxorubicin (S1208, Selleck Chemicals), YM155 (S1130, Selleck Chemicals), BI-2536 (S1109, Selleck Chemicals), AZD7762 (S1532, Selleck Chemicals), trametinib (S2673, Selleck Chemicals) were dissolved in DMSO at 10 mmol/L final concentration. Cisplatin (S1166, Selleck Chemicals) was dissolved in H2O, 0.9% NaCl, 0.3% Tween20 at a 0.5 g/L final concentration, carboplatin (S1215, Selleck Chemicals) was dissolved in H2O, 0.1% Triton. Cells were seeded at 2,500–4,500 cells/well (Supplementary Table S10). After overnight incubation at 37°C and 5% CO2, cells were treated with five different concentrations of drugs (0.001, 0.01, 0.1, 1, and 10 μmol/L) using HP D300 digital dispenser (Tecan). Growth inhibition was measured 72 hours after treatment with MTS diluted 1:6 in fresh culture medium. Cell viability was assessed by recording absorbance at 490 nm using a FLUOstar microplate reader. Dose–response curves were performed using GraphPad Prism 6 Software to determine two parameters reflecting drug sensitivity: GI50 and the area under the curve (AUC). When the GI50 was not reached, the values were set to the highest concentration tested (10 μmol/L). Each concentration was tested in duplicate. Polyclonal goat anti-IGF2 antibody (Ref AF-292-NA; Biotechne) and control anti-Igg goat antibody (Ref AB-108-C; Bio-techne) were tested using five concentrations (0.1, 0.5, 1, 5, and 10 μg/mL). Monoclonal human anti-IGF1/IGF2 therapeutic antibody (Xentuzumab, #TAB-475CQ, Creative Biolabs) was resuspended in PBS at 1 mg/mL final concentration and tested at 1 μmol/L.

In Vivo Xenograft Treatments

Mice were housed in a specific pathogen-free facility, and experiments were conducted using protocols and conditions approved by the institutional animal ethical committee (Authorization no. 2015082610113065.01, Ethics Committee Paris-Nord C2EA 121). At day 0, 5 × 106 HepG2 cells with 1:1 matrigel were inoculated subcutaneously on each flank of female BALB/cAnNRj-Foxn1 five-week-old mice. At day 20, when tumors started growing (volume >200 mm3), mice were randomized in four groups of six mice allocated to the following treatment arms: cisplatin, cisplatin/doxorubicin, BI-2536, and vehicle. Mice in the cisplatin group were injected intraperitoneally at days 20 and 27 with 5 mg/kg cisplatin. Mice from the cisplatin/doxorubicin group received one injection of cisplatin at day 20 (5 mg/kg) and one intraperitoneal doxorubicin at day 27 (2 mg/kg). Mice in the BI-2536 group were intravenously injected with 40 mg/kg BI-2536 at days 20 and 27, whereas mice in the vehicle control group received injections of H2O, 0.9% NaCl. Mice were weighted two times a week, and tumor volume was measured at the same time using a caliper and the following formula: (length*width²)/2 and was expressed as a percentage of initial tumor volume at day 20.

Generation of SBS35 Signature In Vitro

To test cisplatin ability to induce SBS35 mutational signature in hepatoblastoma, six pediatric liver tumor cell lines were treated with 0.5 μmol/L cisplatin for four weeks. At day 0, 60,000 cells were seeded in a 6-well plate and reseeded each week at the same concentration. When mortality was too high, cells were not split and fresh medium supplemented with 0.5 μmol/L cisplatin was added. At day 28, cells were grown in fresh DMEM and limit dilutions were performed in order to expand clonal cell lines. After amplification, samples were extracted using the AllPrep DNA/RNA/miRNA universal extraction kit (Qiagen). Finally, mutational signature analysis was performed on cells derived from limit dilution and compared with bulk nontreated baseline mutational profile.

Methylation-Specific Multiplex Ligation-Dependent Probe Amplification

Multiplex ligation-dependent probe amplification (MLPA) was used to determine the status of locus 11p15.5 in 27 samples without RRBS data available. MLPA was carried out with 50–100 ng of DNA diluted in Tris-HCl, 0.1 mmol/L EDTA. DNA samples were screened using ME030-C3 BWS/RSS kit (MRC-Holland) containing 42 (methylation-specific) MLPA probes: 26 probes specific of BWS 11p15 region, 2 probes targeting NSD1 gene, 13 reference probes targeting other chromosomes, and 1 digestion control probe. Two reference probes targeting, respectively, 2p25 and 2q24 were excluded because of the presence of very frequent copy-number variation in hepatoblastoma tumors. Among the 26 probes targeting the 11p15 region, 10 were methylation-specific, allowing for the detection of methylation abnormalities. To determine thresholds for IC1 GOM and IC2 LOM, we used the K-means clustering method as well. Samples with a gain of methylation in at least two of three probes targeting IC1 were annotated as GOM IC1, whereas samples with a loss of methylation in at least two of four probes were considered LOM IC2. One probe covering IC1 (H19.11.001.976583) was excluded from analysis because its distribution among samples was not discriminating samples with and without GOM IC1.

In Situ Hybridization IGF2 in FFPE Slides

In situ hybridization experiments were performed using RNAscope 2.5 HD Detection Reagent BROWN Kit (cat. #322300) and IGF2 probe (cat. #594361) according to the manufacturer's protocol. FFPE slides underwent target retrieval under standard conditions (15 minutes in target retrieval reagent >98°C) and a 20-minute protease digestion at 40°C in the HybEZ oven.

IHC

IHC analyses of anti–β-catenin (BD Transduction, clone 14, 610154, 1/300) and anti-GS (Bioscience, 1/500) were performed as previously described (33).

TERT Promoter Screening

TERT promoter mutations were identified with WGS data when available and completed with Sanger sequencing as previously described (70) for other samples.

Data Availability

The sequencing data reported in this paper have been deposited to the EGA (European Genome–Phenome Archive) database (accession numbers EGAS EGAS00001005108 and EGAS00001003536).

T.Z. Hirsch reports grants from Cancéropôle Ile-de-France, Fondation d'Enterprise Bristol-Myers Squibb pour la Recherche en Immuno-Oncologie, and CARPEM during the conduct of the study. J. Pilet reports grants and personal fees from Fondation pour la recherche médicale during the conduct of the study. G. Morcrette reports grants from Assistance Publique Hôpitaux de Paris, CARPEM, and Ligue Nationale Contre le Cancer during the conduct of the study. A. Roehrig reports grants from Fondation pour la Recherche médicale during the conduct of the study. L. Molina reports grants from Chateaubriand Fellowship, Embassy of France in the US, Office of Science and Technology during the conduct of the study. S. Caruso reports grants from Labex OncoImmunology and CARPEM during the conduct of the study. E. Jacquemin reports personal fees from Laboratoire CTRS, France, and Vivet Therapeutics, France, outside the submitted work. J. Zucman-Rossi reports grants from France Génomique, Institut National du Cancer, SFCE, association Etoile de Martin, Fédération Enfants et Santé, association Hubert Gouin “Enfance et Cancer,” INSERM Plan Cancer, Cancéropôle Ile-de-France, Inserm ITMO Cancer 3R, Ligue Nationale contre le Cancer, Labex OncoImmunology, IREB, Fondation Bettencourt-Shueller, SIRIC CARPEM, Fondation pour la recherche médicale (FRM), Ligue Contre le Cancer Comité de Paris, Fondation Mérieux, Fondation d'Entreprise Bristol-Myers Squibb, and INCA during the conduct of the study; personal fees from CAIXA Fundation and EASL JHEP Report outside the submitted work. No disclosures were reported by the other authors.

T.Z. Hirsch: Conceptualization, investigation, visualization, writing–original draft. J. Pilet: Conceptualization, investigation, visualization, writing–review and editing. G. Morcrette: Conceptualization, visualization, resources, and investigation. A. Roehrig: Investigation and visualization. B.J.E. Monteiro: Investigation and visualization. L. Molina: Investigation. Q. Bayard: Investigation. E. Trepo: Data curation and investigation. L. Meunier: Investigation and visualization. S. Caruso: Investigation. V. Renault: Resources and data curation. J.-F. Deleuze: Resources. B. Fresneau: Resources. C. Chardot:Resources. E. Gonzales: Resources. E. Jacquemin: Resources. F. Guerin: Resources. M. Fabre: Resources. I. Aerts: Resources. S. Taque: Resources. V. Laithier: Resources. S. Branchereau: Resources. C. Guettier: Resources. L. Brugieres: Resources. S. Rebouissou: Conceptualization, supervision, and investigation. E. Letouze: Conceptualization, supervision, investigation, visualization, writing–original draft. J. Zucman-Rossi: Conceptualization, supervision, investigation, writing–original draft.

We thank all the clinicians, radiologists, pathologists, laboratory technicians, and scientists who participated in this work: Jayendra Shinde, Sandrine Imbeaud, Samantha Schaeffer, Barkha Gupta (Inserm UMR 1138); Yves Lebouc (Inserm UMR 938), Isabelle Janoueix (Inserm UMR 830); Anne Boland, Robert Olaso (CEA, CNRGH); Pr. Frédéric Gauthier, Pr. Hélène Martelli, Cécile Jeannot, Dr. Danièle Pariente, Dr. Stéphanie Franchi-Abella, Isabelle Friteau, Dr. Charlotte Mussini, Dr. Jeanine Quillard, Katia Posseme, Olivier Trassard, Martine Prsle, Véronique Bruna, Marie-Josée Redon, Véronique Mamesier, Aurélie Jove-Garcia, Damien Schmidt, Ludivine Meni, Violette Gendre, Delphine Colmant, Christophe Avignon, Julie Tisserand (CHU Bicêtre, APHP, Le Kremlin Bicêtre); Dr. Dominique Debray, Gisèle Le Gall (CHU Necker, APHP, Paris); Pr. Michel Peuchmaur, Pr. Dominique Berrebi, Pascal Blain, Muriel Mateo-Gajardo, Elodie Dauvergne, Maxette Pierin, Louise Le Guen (CHU Robert Debré, APHP, Paris); Pr. Aurore Coulomb L'hermine, Dr. Sabah Boudjemaa, Dr. Hélène Dariane (CHU Trousseau, APHP, Paris); Dr. Julien Calderaro (CHU Henri Mondor, APHP, Créteil); Pr. Jean Yves Scoazec, Brenda Mallone (Institut Gustave Roussy, Villejuif); Dr. Paul Freneaux, Dr. Daniel Orbach (Institut Curie, Paris); Dr. Carole Coze, Pr. Dominique Figarella Branger, Dr. Andre Maues de Paula (CHU de la Timone, Marseille); Pr. Valérie Rigau, Dr. Yuri Musizzano (CHU Gui de Chauliac, Montpellier); Pr. Jannick Selves (CHU Toulouse); Dr. Estelle Thebaud, Dr. Louise Galmiche, Pr. Jean-François Mosnier (CHU Nantes); Dr. Cécile Dumesnil, Dr. Arnaud François (CHU Charles Nicolle, Rouen); Dr. Anne Lutin, Dr. Antoine Gourmel, Dr. Gilles Morin (CHU Amiens); Dr. Claire Berger (CHU Saint-Etienne); Dr. Frédérique Dijoud, Dr. Cécile Picard (CHU Lyon); Dr Clémence Yguel (CHU Nancy); Charles Marcaillou (Integragen); Josh Waterfall (Institut Curie, INSERM U830), Andrei Zinovyev (Institut Curie, INSERM U900); HEPATOBIO database. We warmly thank Xin Chen, USCF (USA) and Karl-Dimiter Bissig, Baylor College, for their kind gift of the Hep293TT and B6-2 cell lines. This work was supported by France Génomique (GEPELIN project), Institut National du Cancer (INCa, PELICAN.resist project), SFCE, association Etoile de Martin, Fédération Enfants et Santé (FES), association Hubert Gouin “Enfance et Cancer,” INSERM with the HTE (plan Cancer), Cancéropôle Ile-de-France “ThéBioCaPe” project, Inserm ITMO Cancer 3R “PEILICANS” project. The group is supported by the Ligue Nationale contre le Cancer (Equipe Labellisée), Labex OncoImmunology (investissement d'avenir), grant IREB, Coup d'Elan de la Fondation Bettencourt-Shueller, the SIRIC CARPEM, FRM prix Rosen, Ligue Contre le Cancer Comité de Paris (prix René et Andre Duquesne), and Fondation Mérieux. T.Z. Hirsch was supported by a fellowship from Cancéropôle Ile-de-France, Fondation d'Entreprise Bristol-Myers Squibb pour la Recherche en Immuno-Oncologie, and CARPEM; G. Morcrette was supported by CARPEM, the Ligue Nationale contre le Cancer and Assistance Publique Hôpitaux de Paris. Q. Bayard was funded by a fellowship from the Ministry of Education and Research. This work was supported by the Fondation pour la Recherche Médicale, grant number ECO201906008977 to A. Roehrig and grant number ECO20170637540 to J. Pilet. L. Molina was funded by the Chateaubriand Fellowship from the Embassy of France in the USA. S. Caruso was funded by Labex OncoImmunology and CARPEM. HEPATOBIO was supported by SFCE, Etoile de Martin, Enfant Sans Cancer, INCA/DGOS (Grant INCA-TRANSLA N° 2013-209).

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.

1.
Darbari
A
,
Sabin
KM
,
Shapiro
CN
,
Schwarz
KB
. 
Epidemiology of primary hepatic malignancies in U.S. children
.
Hepatology
2003
;
38
:
560
6
.
2.
Hadzic
N
,
Finegold
MJ
. 
Liver neoplasia in children
.
Clin Liver Dis
2011
;
15
:
443
62
, vii–x.
3.
Arai
Y
,
Honda
S
,
Haruta
M
,
Kasai
F
,
Fujiwara
Y
,
Ohshima
J
, et al
Genome-wide analysis of allelic imbalances reveals 4q deletions as a poor prognostic factor and MDM4 amplification at 1q32.1 in hepatoblastoma
.
Genes Chromosomes Cancer
2010
;
49
:
596
609
.
4.
Eichenmüller
M
,
Trippel
F
,
Kreuder
M
,
Beck
A
,
Schwarzmayr
T
,
Häberle
B
, et al
The genomic landscape of hepatoblastoma and their progenies with HCC-like features
.
J Hepatol
2014
;
61
:
1312
20
.
5.
Jia
D
,
Dong
R
,
Jing
Y
,
Xu
D
,
Wang
Q
,
Chen
L
, et al
Exome sequencing of hepatoblastoma reveals novel mutations and cancer genes in the Wnt pathway and ubiquitin ligase complex
.
Hepatology
2014
;
60
:
1686
96
.
6.
Sumazin
P
,
Chen
Y
,
Treviño
LR
,
Sarabia
SF
,
Hampton
OA
,
Patel
K
, et al
Genomic analysis of hepatoblastoma identifies distinct molecular and prognostic subgroups
.
Hepatology
2017
;
65
:
104
21
.
7.
Sekiguchi
M
,
Seki
M
,
Kawai
T
,
Yoshida
K
,
Yoshida
M
,
Isobe
T
, et al
Integrated multiomics analysis of hepatoblastoma unravels its heterogeneity and provides novel druggable targets.
NPJ Precis Oncol
2020
;
4
:
1
12
.
8.
Gröbner
SN
,
Worst
BC
,
Weischenfeldt
J
,
Buchhalter
I
,
Kleinheinz
K
,
Rudneva
VA
, et al
The landscape of genomic alterations across childhood cancers
.
Nature
2018
;
555
:
321
7
.
9.
Khanna
R
,
Verma
SK
. 
Pediatric hepatocellular carcinoma
.
World J Gastroenterol
2018
;
24
:
3980
99
.
10.
Iannelli
F
,
Collino
A
,
Sinha
S
,
Radaelli
E
,
Nicoli
P
,
D'Antiga
L
, et al
Massive gene amplification drives paediatric hepatocellular carcinoma caused by bile salt export pump deficiency
.
Nat Commun
2014
;
5
:
3850
.
11.
Haines
K
,
Sarabia
SF
,
Alvarez
KR
,
Tomlinson
G
,
Vasudevan
SA
,
Heczey
AA
, et al
Characterization of pediatric hepatocellular carcinoma reveals genomic heterogeneity and diverse signaling pathway activation
.
Pediatr Blood Cancer
2019
;
66
:
e27745
.
12.
Honeyman
JN
,
Simon
EP
,
Robine
N
,
Chiaroni-Clarke
R
,
Darcy
DG
,
Lim
IIP
, et al
Detection of a recurrent DNAJB1-PRKACA chimeric transcript in fibrolamellar hepatocellular carcinoma
.
Science
2014
;
343
:
1010
4
.
13.
Franchi-Abella
S
,
Branchereau
S
. 
Benign hepatocellular tumors in children: focal nodular hyperplasia and hepatocellular adenoma
.
Int J Hepatol
2013
;
2013
:
215064
.
14.
López-Terrada
D
,
Alaggio
R
,
de Dávila
MT
,
Czauderna
P
,
Hiyama
E
,
Katzenstein
H
, et al
Towards an international pediatric liver tumor consensus classification: proceedings of the Los Angeles COG liver tumors symposium
.
Mod Pathol
2014
;
27
:
472
91
.
15.
Feng
T-C
,
Zai
H-Y
,
Jiang
W
,
Zhu
Q
,
Jiang
B
,
Yao
L
, et al
Survival and analysis of prognostic factors for hepatoblastoma: based on SEER database
.
Ann Transl Med
2019
;
7
:
555
.
16.
Carrillo-Reixach
J
,
Torrens
L
,
Simon-Coma
M
,
Royo
L
,
Domingo-Sàbat
M
,
Abril-Fornaguera
J
, et al
Epigenetic footprint enables molecular risk stratification of hepatoblastoma with clinical implications
.
J Hepatol
2020
;
73
:
328
41
.
17.
Alexandrov
LB
,
Nik-Zainal
S
,
Wedge
DC
,
Aparicio
SAJR
,
Behjati
S
,
Biankin
AV
, et al
Signatures of mutational processes in human cancer
.
Nature
2013
;
500
:
415
21
.
18.
Alexandrov
LB
,
Kim
J
,
Haradhvala
NJ
,
Huang
MN
,
Tian Ng
AW
,
Wu
Y
, et al
The repertoire of mutational signatures in human cancer
.
Nature
2020
;
578
:
94
101
.
19.
Pilati
C
,
Shinde
J
,
Alexandrov
LB
,
Assié
G
,
André
T
,
Hélias-Rodzewicz
Z
, et al
Mutational signature analysis identifies MUTYH deficiency in colorectal cancers and adrenocortical carcinomas
.
J Pathol
2017
;
242
:
10
5
.
20.
Boot
A
,
Huang
MN
,
Ng
AWT
,
Ho
S-C
,
Lim
JQ
,
Kawakami
Y
, et al
In-depth characterization of the cisplatin mutational signature in human cell lines and in esophageal and liver tumors
.
Genome Res
2018
;
28
:
654
65
.
21.
Galluzzi
L
,
Senovilla
L
,
Vitale
I
,
Michels
J
,
Martins
I
,
Kepp
O
, et al
Molecular mechanisms of cisplatin resistance
.
Oncogene
2012
;
31
:
1869
83
.
22.
Martinez-Quetglas
I
,
Pinyol
R
,
Dauch
D
,
Torrecilla
S
,
Tovar
V
,
Moeini
A
, et al
IGF2 is up-regulated by epigenetic mechanisms in hepatocellular carcinomas and is an actionable oncogene product in experimental models
.
Gastroenterology
2016
;
151
:
1192
205
.
23.
Caruso
S
,
Calatayud
A-L
,
Pilet
J
,
La Bella
T
,
Rekik
S
,
Imbeaud
S
, et al
Analysis of liver cancer cell lines identifies agents with likely efficacy against hepatocellular carcinoma and markers of response
.
Gastroenterology
2019
;
157
:
760
76
.
24.
Chang
MT
,
Bhattarai
TS
,
Schram
AM
,
Bielski
CM
,
Donoghue
MTA
,
Jonsson
P
, et al
Accelerating discovery of functional mutant alleles in cancer
.
Cancer Discov
2018
;
8
:
174
83
.
25.
Coorens
THH
,
Treger
TD
,
Al-Saadi
R
,
Moore
L
,
Tran
MGB
,
Mitchell
TJ
, et al
Embryonal precursors of Wilms tumor
.
Science
2019
;
366
:
1247
51
.
26.
Cairo
S
,
Armengol
C
,
De Reyniès
A
,
Wei
Y
,
Thomas
E
,
Renard
C-A
, et al
Hepatic stem-like phenotype and interplay of Wnt/beta-catenin and Myc signaling in aggressive childhood liver cancer
.
Cancer Cell
2008
;
14
:
471
84
.
27.
Hooks
KB
,
Audoux
J
,
Fazli
H
,
Lesjean
S
,
Ernault
T
,
Dugot-Senant
N
, et al
New insights into diagnosis and therapeutic options for proliferative hepatoblastoma
.
Hepatology
2018
;
68
:
89
102
.
28.
Northcott
PA
,
Buchhalter
I
,
Morrissy
AS
,
Hovestadt
V
,
Weischenfeldt
J
,
Ehrenberger
T
, et al
The whole-genome landscape of medulloblastoma subtypes
.
Nature
2017
;
547
:
311
7
.
29.
van Groningen
T
,
Koster
J
,
Valentijn
LJ
,
Zwijnenburg
DA
,
Akogul
N
,
Hasselt
NE
, et al
Neuroblastoma is composed of two super-enhancer-associated differentiation states
.
Nat Genet
2017
;
49
:
1261
6
.
30.
Boeva
V
,
Louis-Brennetot
C
,
Peltier
A
,
Durand
S
,
Pierre-Eugène
C
,
Raynal
V
, et al
Heterogeneity of neuroblastoma cell identity defined by transcriptional circuitries
.
Nat Genet
2017
;
49
:
1408
13
.
31.
Letouzé
E
,
Shinde
J
,
Renault
V
,
Couchy
G
,
Blanc
J-F
,
Tubacher
E
, et al
Mutational signatures reveal the dynamic interplay of risk factors and cellular processes during liver tumorigenesis
.
Nat Commun
2017
;
8
:
1315
.
32.
Pich
O
,
Muiños
F
,
Lolkema
MP
,
Steeghs
N
,
Gonzalez-Perez
A
,
Lopez-Bigas
N
. 
The mutational footprints of cancer therapies
.
Nat Genet
2019
;
51
:
1732
40
.
33.
Morcrette
G
,
Hirsch
TZ
,
Badour
E
,
Pilet
J
,
Caruso
S
,
Calderaro
J
, et al
APC germline hepatoblastomas demonstrate cisplatin-induced intratumor tertiary lymphoid structures
.
Oncoimmunology
2019
;
8
:
e1583547
.
34.
Tesniere
A
,
Schlemmer
F
,
Boige
V
,
Kepp
O
,
Martins
I
,
Ghiringhelli
F
, et al
Immunogenic death of colon cancer cells treated with oxaliplatin
.
Oncogene
2010
;
29
:
482
91
.
35.
Li
H
,
Durbin
R
. 
Fast and accurate short read alignment with Burrows-Wheeler transform
.
Bioinformatics
2009
;
25
:
1754
60
.
36.
McKenna
A
,
Hanna
M
,
Banks
E
,
Sivachenko
A
,
Cibulskis
K
,
Kernytsky
A
, et al
The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data
.
Genome Res
2010
;
20
:
1297
303
.
37.
DePristo
MA
,
Banks
E
,
Poplin
R
,
Garimella
KV
,
Maguire
JR
,
Hartl
C
, et al
A framework for variation discovery and genotyping using next-generation DNA sequencing data
.
Nat Genet
2011
;
43
:
491
8
.
38.
Van der Auwera
GA
,
Carneiro
MO
,
Hartl
C
,
Poplin
R
,
Del Angel
G
,
Levy-Moonshine
A
, et al
From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline
.
Curr Protoc Bioinformatics
2013
;
43
:
11.10.1
11.10.33
.
39.
McLaren
W
,
Gil
L
,
Hunt
SE
,
Riat
HS
,
Ritchie
GRS
,
Thormann
A
, et al
The ensembl variant effect predictor
.
Genome Biol
2016
;
17
:
122
.
40.
Karczewski
KJ
,
Francioli
LC
,
Tiao
G
,
Cummings
BB
,
Alföldi
J
,
Wang
Q
, et al
The mutational constraint spectrum quantified from variation in 141,456 humans
.
Nature
2020
;
581
:
434
43
.
41.
Robinson
JT
,
Thorvaldsdóttir
H
,
Winckler
W
,
Guttman
M
,
Lander
ES
,
Getz
G
, et al
Integrative genomics viewer
.
Nat Biotechnol
2011
;
29
:
24
6
.
42.
Chen
X
,
Schulz-Trieglaff
O
,
Shaw
R
,
Barnes
B
,
Schlesinger
F
,
Källberg
M
, et al
Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications
.
Bioinformatics
2016
;
32
:
1220
2
.
43.
Nik-Zainal
S
,
Van Loo
P
,
Wedge
DC
,
Alexandrov
LB
,
Greenman
CD
,
Lau
KW
, et al
The life history of 21 breast cancers
.
Cell
2012
;
149
:
994
1007
.
44.
Popova
T
,
Manié
E
,
Stoppa-Lyonnet
D
,
Rigaill
G
,
Barillot
E
,
Stern
MH
. 
Genome Alteration Print (GAP): a tool to visualize and mine complex cancer genomic profiles obtained by SNP arrays
.
Genome Biol
2009
;
10
:
R128
.
45.
Olshen
AB
,
Venkatraman
ES
,
Lucito
R
,
Wigler
M
. 
Circular binary segmentation for the analysis of array-based DNA copy number data
.
Biostatistics
2004
;
5
:
557
72
.
46.
Schulze
K
,
Imbeaud
S
,
Letouzé
E
,
Alexandrov
LB
,
Calderaro
J
,
Rebouissou
S
, et al
Exome sequencing of hepatocellular carcinomas identifies new mutational signatures and potential therapeutic targets
.
Nat Genet
2015
;
47
:
505
11
.
47.
Lawrence
MS
,
Stojanov
P
,
Polak
P
,
Kryukov
GV
,
Cibulskis
K
,
Sivachenko
A
, et al
Mutational heterogeneity in cancer and the search for new cancer-associated genes
.
Nature
2013
;
499
:
214
8
.
48.
Mularoni
L
,
Sabarinathan
R
,
Deu-Pons
J
,
Gonzalez-Perez
A
,
López-Bigas
N
. 
OncodriveFML: a general framework to identify coding and non-coding regions with cancer driver mutations
.
Genome Biol
2016
;
17
:
128
.
49.
Ma
X
,
Liu
Y
,
Liu
Y
,
Alexandrov
LB
,
Edmonson
MN
,
Gawad
C
, et al
Pan-cancer genome and transcriptome analyses of 1,699 paediatric leukaemias and solid tumours
.
Nature
2018
;
555
:
371
6
.
50.
Shinde
J
,
Bayard
Q
,
Imbeaud
S
,
Hirsch
TZ
,
Liu
F
,
Renault
V
, et al
Palimpsest: an R package for studying mutational and structural variant signatures along clonal evolution in cancer
.
Bioinformatics
2018
;
34
:
3380
1
.
51.
Yates
LR
,
Gerstung
M
,
Knappskog
S
,
Desmedt
C
,
Gundem
G
,
Van Loo
P
, et al
Subclonal diversification of primary breast cancer revealed by multiregion sequencing
.
Nat Med
2015
;
21
:
751
9
.
52.
Kim
D
,
Pertea
G
,
Trapnell
C
,
Pimentel
H
,
Kelley
R
,
Salzberg
SL
. 
TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions
.
Genome Biol
2013
;
14
:
R36
.
53.
Anders
S
,
Pyl
PT
,
Huber
W
. 
HTSeq–a Python framework to work with high-throughput sequencing data
.
Bioinformatics
2015
;
31
:
166
9
.
54.
Love
MI
,
Huber
W
,
Anders
S
. 
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
550
.
55.
Haas
BJ
,
Dobin
A
,
Li
B
,
Stransky
N
,
Pochet
N
,
Regev
A
. 
Accuracy assessment of fusion transcript detection via read-mapping and de novo fusion transcript assembly-based methods
.
Genome Biol
2019
;
20
:
213
.
56.
Ritchie
ME
,
Phipson
B
,
Wu
D
,
Hu
Y
,
Law
CW
,
Shi
W
, et al
limma powers differential expression analyses for RNA-sequencing and microarray studies
.
Nucleic Acids Res
2015
;
43
:
e47
.
57.
Becht
E
,
Giraldo
NA
,
Lacroix
L
,
Buttard
B
,
Elarouci
N
,
Petitprez
F
, et al
Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression
.
Genome Biol
2016
;
17
:
218
.
58.
Slyper
M
,
Porter
CBM
,
Ashenberg
O
,
Waldman
J
,
Drokhlyansky
E
,
Wakiro
I
, et al
A single-cell and single-nucleus RNA-Seq toolbox for fresh and frozen human tumors
.
Nat Med
2020
;
26
:
792
802
.
59.
Zheng
GXY
,
Terry
JM
,
Belgrader
P
,
Ryvkin
P
,
Bent
ZW
,
Wilson
R
, et al
Massively parallel digital transcriptional profiling of single cells
.
Nat Commun
2017
;
8
:
14049
.
60.
Stuart
T
,
Butler
A
,
Hoffman
P
,
Hafemeister
C
,
Papalexi
E
,
Mauck
WM
, et al
Comprehensive integration of single-cell data
.
Cell
2019
;
177
:
1888
902
.
61.
Gu
H
,
Smith
ZD
,
Bock
C
,
Boyle
P
,
Gnirke
A
,
Meissner
A
. 
Preparation of reduced representation bisulfite sequencing libraries for genome-scale DNA methylation profiling. Nature Protocols
.
Nature
2011
;
6
:
468
81
.
62.
Guo
W
,
Fiziev
P
,
Yan
W
,
Cokus
S
,
Sun
X
,
Zhang
MQ
, et al
BS-Seeker2: a versatile aligning pipeline for bisulfite sequencing data
.
BMC Genomics
2013
;
14
:
774
.
63.
Guo
W
,
Zhu
P
,
Pellegrini
M
,
Zhang
MQ
,
Wang
X
,
Ni
Z
. 
CGmapTools improves the precision of heterozygous SNV calls and supports allele-specific methylation detection and visualization in bisulfite-sequencing data
.
Bioinformatics
2018
;
34
:
381
7
.
64.
Akalin
A
,
Kormaksson
M
,
Li
S
,
Garrett-Bakelman
FE
,
Figueroa
ME
,
Melnick
A
, et al
methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles
.
Genome Biol
2012
;
13
:
R87
.
65.
Meunier
L
,
Hirsch
TZ
,
Caruso
S
,
Imbeaud
S
,
Bayard
Q
,
Roehrig
A
, et al
DNA methylation signatures reveal the diversity of processes remodeling hepatocellular carcinoma methylomes
.
Hepatology
2021
;
doi: 10.1002/hep.31796
.
66.
Roadmap Epigenomics Consortium
,
Kundaje
A
,
Meuleman
W
,
Ernst
J
,
Bilenky
M
,
Yen
A
, et al
Integrative analysis of 111 reference human epigenomes
.
Nature
2015
;
518
:
317
30
.
67.
Salhab
A
,
Nordström
K
,
Gasparoni
G
,
Kattler
K
,
Ebert
P
,
Ramirez
F
, et al
A comprehensive analysis of 195 DNA methylomes reveals shared and cell-specific features of partially methylated domains
.
Genome Biol
2018
;
19
:
150
.
68.
ENCODE Project Consortium
. 
An integrated encyclopedia of DNA elements in the human genome
.
Nature
2012
;
489
:
57
74
.
69.
Silva
TC
,
Coetzee
SG
,
Gull
N
,
Yao
L
,
Hazelett
DJ
,
Noushmehr
H
, et al
ELMER v.2: an R/Bioconductor package to reconstruct gene regulatory networks from DNA methylation and transcriptome profiles
.
Bioinformatics
2019
;
35
:
1974
7
.
70.
Nault
JC
,
Mallet
M
,
Pilati
C
,
Calderaro
J
,
Bioulac-Sage
P
,
Laurent
C
, et al
High frequency of telomerase reverse-transcriptase promoter somatic mutations in hepatocellular carcinoma and preneoplastic lesions
.
Nat Commun
2013
;
4
:
2218
.