DNA Methylation-Based Prognosis and Epidrivers in Hepatocellular Carcinoma

Human Samples and Molecular Profiling.

Initially, the study included samples from 331 surgically resected HCC and 19 nontumor tissues, including 9 cirrhosis and 10 normal livers. The training set (Heptromic data set, n = 248; flow chart in Supporting Fig. 1) were samples obtained from two institutions of the HCC Genomic Consortium: IRCCS Istituto Nazionale Tumori (Milan, Italy; n = 217) and Hospital Clínic (Barcelona, Spain; n = 31). All samples included in this study were fresh-frozen. For RNA and DNA extraction, we used the Qiagen RNeasy Mini (500 ng of total RNA at a concentration of 100 ng/μL; Qiagen, Hilden, Germany) and Invitrogen ChargeSwitch genomic DNA Mini Tissue (1 μg of total DNA at a concentration of 100 ng/μL; Invitrogen, Carlsbad, CA) kits, respectively. Median sample storage time from collection to DNA/RNA extraction was 7 years. RNA profiling was conducted on 228 HCC and 168 nontumor liver adjacent cirrhotic tissues using the Affymetrix Human Genome U219 Array Plate (Affymetrix, Inc., Santa Clara, CA), which is able to interrogate more than 20,000 mapped genes. After hybridization, only one tumor sample was discarded for transcriptome analysis owing to poor quality.

Methylome profiling was performed on 248 samples with the Illumina Infinium HumanMethylation450 BeadChip array (Illumina, Inc., San Diego, CA) that interrogates more than 485,000 cytosine-phosphate-guanine (CpG) sites covering 96% of known CpG islands.¹² Of the 248 HCC samples from the training set, 221 HCC qualified for final analyses after quality filtering (i.e., less than 5% of the probes incorrectly interrogated; P < 0.01). In 205 patients, there was information on both transcriptome and methylation profiling. For validation of the methylation-based signature, we analyzed data from a cohort of 83 HCC patients treated with resection in two French institutions (Bordeaux and Créteil hospitals), profiled with the same technology as the training set. The institutional review boards of the participating centers approved the study.

Pyrosequencing Validation.

DNA methylation was evaluated with a pyrosequencing assay in a subset of samples previously analyzed by the Illumina Infinium Human Methylation450 array. A minimum of 500 ng of DNA was converted using the EZ DNA Methylation-Gold (Zymo Research Corporation, Irvine, CA) bisulfite conversion kit, following the manufacturer’s recommendations. Specific sets of primers for polymerase chain reaction (PCR) amplification and sequencing were designed using specific software (PyroMark assay design, version 2.0.01.15). These are summarized in Supporting Table 1. Primer sequences were designed, when possible, to hybridize with CpG-free sites to ensure methylation-independent amplification. PCR was performed under standard conditions with biotinylated primers, and the PyroMark Vacuum Prep Tool (Biotage AB, Uppsala, Sweden) was used to prepare single-stranded PCR products, according to the manufacturer’s instructions. PCR products were observed at 2% agarose gels before pyrosequencing. Reactions were performed in a PyroMark Q96 System (version 2.0.6; Qiagen) using appropriate reagents and protocols.

Data Analysis.

To study differential methylation between HCC and normal liver tissue, probes containing single-nucleotide polymorphisms (SNPs) or located on sex chromosomes were eliminated, leaving 434,728 probes for analysis. Only those with a high signal quality (P < 0.01) were considered. Probes hypomethylated (B < 0.33) in at least 90% of normal liver samples and hypermethylated (B > 0.5) in at least 5% of tumors (31,052 CpG sites) and probes hypermethylated (B > 0.5) in at least 90% of normal liver and hypomethylated (B < 0.33) in at least 5% of tumors were selected (68,086 CpG sites). The beta value (B) is used to estimate the methylation level of the CpG locus using the ratio of intensities between methylated and unmethylated alleles. For prognosis purposes, analyses focused only on CpG islands located in TSS1500, TSS200, 5′ UTR (untranslated region), and 1stExon (n = 84,448) shown to be differentially methylated in tumors versus normal liver according to predefined criteria (n = 11,307).¹³ This enabled us to select for those methylation markers primarily deregulated in cancer tissues.

In order to discover potential candidate HCC epidrivers, an F score directly proportional to intergroup variability (normal liver versus HCC) and inversely proportional to normal sample intragroup variability was calculated for all probes located in promoters and CpG islands. The first 500 ranked probes were grouped per gene. Genes with more than five hits were further studied as epidrivers. Previously reported epigenetically deregulated genes were also analyzed, considering the mean of all probes located in the TSS200 region, which shows the best correlation with expression. For those samples without TSS200 probes in the array, TSS1500 were used instead, as the closest marker available in the array.

We generate a methylome signature able to predict risk of death for each individual measured by a mortality index (MI). The prediction signature was generated using the random survival forest (RSF) method developed for variable selection.¹³ RSF is a method for prediction and variable selection for right-censored survival and competing risk data by growing survival trees to estimate a cumulative hazard function (CHF), which derives from each tree of the RSF. Input variables were the filtered 11,307 CpG sites. First, probes were randomly split in six sets of 2,000 or less probes to fit a model each using the RSF method by growing 1,000 survival trees. Variable importance (VIMP) scores were computed for all the probes used to grow the trees, and the 350 most informative from each model were selected to fit a new model using the same RSF method. Once the model was built, we removed redundant and “noisy” variables by generating incremental RSF models adding one more variable, ranked by the VIMP calculated in the previous model, at each step, getting an estimate of the error as a measure of how much misclassification increases, or decreases, for a new test case if a given variable was not available for that case. The final set of selected variables (36 probes) was used to build the model having the lowest error rate using only the two most significant decimals of the obtained error rate. Thus, the methylation signature was based upon those 36 probes. The predictive accuracy of the 36-signature was examined by bootstrap cross-validation.¹³ This allows estimating the prediction error curves of different survival models on 100 bootstrap samples and comparing the performance of the different models (i.e., Cox’s and RSF) with a reference (estimation for the whole group without splitting based on the prediction variables; Supporting Fig. 2).

From the last RSF model, we generated a DNA methylation-based mortality risk score (i.e., MI; range, 36–360) for every individual computed as a sum of the CHF for each patient evaluated at a set of distinct time points weighed by the number of individuals at risk at the different time points. High MI correlates with higher risk of death. Using this index, we were able to split our training set into two risk groups based on percentiles, selecting the patients with the 20% highest levels (high risk [MI >230] and nonhigh risk [MI <230]) and test the variables selected as a prediction signature by confronting the high- and low-risk groups. To test the significance of the MI index, the same MI value used to split the training set was also used for the validation set. In other words, for the validation, we did not use percentiles and split patients using the actual same value of the MI generated in the training set.

Processing of transcriptome data (i.e., normalization, background correction, and filtering) was conducted as previously reported.⁵ Prediction of liver cancer mRNA-based signatures in the training set was performed using the nearest template prediction method, as implemented in the specific module of Gene Pattern software.¹⁴ All mRNA signatures analyzed were already reported, being deposited in the Molecular Signature Database (www.broadinstitute.org/gsea/msigdb). To provide biological insight on samples with high mortality risk score (>230; percentile, 80), we used the Gene Set Enrichment Analysis (GSEA), as implemented in Gene Pattern. Gene ontology by PANTHER, INTERPRO, and KEGG pathways enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID; v6.7).

Outcome analyses were framed within current guidelines from the Progress Partnership for prognostic factor research.¹⁵ Kaplan-Meier’s plots and Cox’s regression were used to assess association with outcome of the methylation signature and clinical variables. Overall survival was defined as the time between surgical resection and death of any cause or last follow-up, whereas cancer-related survival (validation set) was the time between resection to HCC-related death, as previously defined.¹⁶ Time to recurrence is the time between resection and the radiological evidence of first tumor recurrence. Several clinical variables previously reported as outcome predictors in HCC were included in the analysis: gender; age; etiology of liver disease; tumor size; serum albumin; serum bilirubin; platelet count; serum alpha-fetoprotein (AFP); Barcelona Clinic Liver Cancer (BCLC) stage¹⁷; microvascular invasion; satellites; and degree of differentiation. Variables with a P value less than 0.05 in the univariate analysis (log-rank test) were separately evaluated in a multivariate Cox’s model aimed to identify independent predictors of survival and recurrence. Missing values were less than 10% for each variable included in the Cox’s model. Analyses were performed using the GenePattern analytical toolkit (www.broad.mit.edu/cancer/software/genepattern/) and the R statistical package (www.r-project.org).

Methylation Signature Predicts Survival in HCC.

RSFs enabled us to generate a methylation-based signature that accurately discriminated patients based on their survival (training set; n = 221). Table 1 summarizes the clinical characteristics of this cohort, that mostly included males (172; 78%), with a median of 66 years of age, viral-related HCC (hepatitis C virus [HCV]: 101 [47%]; hepatitis B virus [HBV]: 44 [20%]), and a median tumor size of 3.5 cm. Patients had predominantly uninodular disease (166; 75%) and no microvascular invasion (142; 65%) or satellites (158; 71%), being the majority at early clinical stages (BCLC O/A: 191 [87%]) with a median follow-up of 48.5 months. Using the array data, we generated a methylation-based signature formed by 36 unique probes (Supporting Table 2) that enable us to provide a DNA-methylation–based risk score of mortality (MI, ranging from 36 to 360; Fig. 1A). As predicted, these 36 probes were hypermethylated in tumor, compared to nontumor, tissue (Supporting Fig. 3). To ensure reproducibility of the methylation status of these probes, we performed technical validation and determine methylation levels by pyrosequencing in a subset of 10 probes in 20 HCCs. There was a significant and high concordance on the methylation status of these probes using both methods (Supporting Fig. 4).

The MI is a metric based on the weighed hazard of death for an individual across all trees (see Materials and Methods for details). This DNA-methylome–based score of mortality risk significantly correlated with known clinical and pathological predictors of survival, such as satellites (P = 0.02), multinodularity (P = 0.002), vascular invasion (P < 0.001), BCLC staging (P = 0.001), AFP (P = 0.001), bilirubin (P = 0.01), platelet count (P = 0.009), and albumin (P = 0.03; Fig. 1B-I). Concordantly, univariate analysis showed that DNA-methylation MI significantly correlated with risk of death (P < 0.001; Fig. 2A). Median survival was significantly lower for patients with DNA-methylation MI >230—which represents 20% of the cohort—than for those with MI <230 (13 vs. 79.7 months, respectively; P < 0.001; Fig. 2B). Multivariate Cox’s modeling confirmed DNA-methylation MI as an independent predictor of survival (hazard ratio [HR]: 13.35; 95% confidence interval [CI]: 7.94–22.42; P < 0.001) along with multinodularity (HR, 1.65; 95% CI: 1.03–2.62; P < 0.001) and platelet count (HR, 1.58; 95% CI: 1.01–2.41; P = 0.03; Table 2). Similarly, DNA-methylation MI score predicted tumor recurrence (Fig. 2C) and Cox’s modeling showed that the signature was also an independent predictor of overall recurrence (HR, 5.8; 95% CI: 3.1–11; P < 0.001) along with multinodularity (HR, 1.8; 95% CI: 1.1–3.0; P =0.007; Table 2).

We next sought to validate the prognostic performance of our methylation-based MI in an independent data set. The validation cohort also included 83 HCC patients treated by surgical resection from different French institutions. Unlike the training set, the main etiology of these patients was related to alcohol intake (47%; 39 of 83), presented with more aggressive disease (microvascular invasion 54%, satellites 42%, and tumors >3 cm in 83%) and had shorter median follow-up time (35 months). Patients with a high DNA-methylation–based MI had a significantly lower cancer-related survival (P = 0.01; Fig. 2D). Even when arbitrary categorizing MI by percentiles, there is a direct relation between MI and poor outcome, both in training and validation sets (Supporting Fig. 5). Overall, these data demonstrate how epigenetic deregulation correlates with patient outcome in surgically resected HCC and introduces a first-in-class DNA-validated, methylation-based signature.

Integrated Prognostic Classification With Transcriptome-Based Predictors.

In order to further characterize the subset of patients identified by the DNA-methylation signature, we integrate whole-genome mRNA and methylation data of 205 HCCs. DNA-methylation–based high MI was significantly enriched in tumors that harbored mRNA signatures capturing progenitor cell features, such as epithelial cell adhesion molecule (EpCAM)¹⁸ (P = 0.009) and S2 (a reported molecular class of potential progenitor cell origin¹⁹; P = 0.006; Fig. 3). This was further confirmed when we performed GSEA on the samples with the highest MI (Supporting Table 3). In addition, we found a significant enrichment of Gene Ontologies related to chromosomal dynamics and RNA processing in these samples (false discovery rate: <0.05). Conversely, tumors with DNA-methylation–based low MI were enriched in gene sets associated with metabolism and homeostasis (Supporting Table 3).

Noncoding mRNA expression profile and DNA methylation data were tested in 205 samples. Patients with a high DNA-methylation MI had a differential noncoding mRNA expression pattern comprising five microRNAs (miRNAs) and one small nucleolar RNA: miR-3193 (P < 0.001); miR-27b-star (P = 0.002); miR-422a (P = 0.02); miR-378c (P = 0.03); miR-938 (P = 0.03); and SNORD126 (annotated as miR-1201 in the array; P = 0.03). All six noncoding RNAs were significantly down-regulated in patients with DNA-methylation high MI score (Supporting Fig. 6).

Aberrant Methylation in Human HCC.

Similarly to other malignancies, there was a remarkable predominance of hypomethylated probes in HCC, compared to normal liver (68%; Fig. 4A). Hypomethylated probes were mainly located in the intergenic (39.9%) and body regions (34.5%), whereas hypermethylated probes were predominantly located in promoter areas (50.5%). Regarding CpG island relation, hypomethylated probes were mainly located in open sea regions (63.55%), whereas hypermethylated probes tended to locate in CpG islands (63.9%) and shores (24.8%; Fig. 4A).

To identify potential methylation markers of malignant transformation, we searched for probes that could accurately differentiate normal liver tissue and HCC, among the 11,307 above-described probes. Analyses focused on probes with low intragroup variability (tumor and nontumor), but high intergroup variability, as defined by an F score (mean-β differences between groups/mean-β differences within groups). The top-100–ranked probes using this score were able to accurately differentiate normal liver and HCC tissues (Fig. 4B; Supporting Table 4). These 100 probes were all hypermethylated in HCC and corresponded to 70 different genes. Gene ontology analyses showed enrichment in entries related to gene expression regulation (such as transcription factors [P = 2.9E⁻⁶]) and different signaling cascades, such as insulin growth factor (IGF; P = 2.1E⁻³), phosphatidylinositol 3-kinase (PI3K; P = 3E⁻³), transforming growth factor beta (TGF-β; P = 5.8E⁻³), and cadherin (P = 3.5E⁻²), among others. When compared to a recently published analysis in 24 HCCs,²⁰ we found 43 common probes (Supporting Fig. 7A; Supporting Table 5) and validated their ability to distinguish HCC from nontumor tissue (Supporting Fig. 7B). We next sought to explore the landscape of known and novel genes candidates as potential epidrivers in HCC. First, we validated other previously reported aberrantly DNA methylated genes in HCC such as hypermethylation of RASSF1 (82%), APC (78%) or NEFH (43%)²¹ and hypomethylation of IGF2 (51%), among others (Table 3). Then, we identified that known drivers in other tumors show aberrant promoter methylation in HCC, such as NOTCH3,²² nuclear receptor binding SET domain protein 1 (NSD1),²² and zinc finger of the cerebellum 1 (ZIC1). Finally, in order to describe novel candidate potential epidrivers not previously described in HCC, we selected the top 500 F-scored probes that were clustered in their target genes (Supporting Table 6) and identified genes involved in TGF-β, or fibroblast growth factor (FGF) signaling. Those genes with the highest number of deregulated hits between HCC and normal liver were highlighted as the most appealing candidates (Table 3). Some of these candidates were also validated with pyrosequencing (Supporting Fig. 8). Among them SEPT9, a tumor suppressor described in colon and ovarian cancer,²³ and ephrin-B2 ligand (EFNB2), of which hypermethylation was reported in patients with acute leukemia.²⁴ These new candidate epidrivers were also found significantly deregulated in the validation cohort of 83 HCCs (Fig. 5).