Epigenetic alterations, such as promoter hypermethylation, may drive cancer through tumor suppressor gene inactivation. However, we have limited ability to differentiate driver DNA methylation (DNAme) changes from passenger events. We developed DNAme driver inference–MethSig–accounting for the varying stochastic hypermethylation rate across the genome and between samples. We applied MethSig to bisulfite sequencing data of chronic lymphocytic leukemia (CLL), multiple myeloma, ductal carcinoma in situ, glioblastoma, and to methylation array data across 18 tumor types in TCGA. MethSig resulted in well-calibrated quantile–quantile plots and reproducible inference of likely DNAme drivers with increased sensitivity/specificity compared with benchmarked methods. CRISPR/Cas9 knockout of selected candidate CLL DNAme drivers provided a fitness advantage with and without therapeutic intervention. Notably, DNAme driver risk score was closely associated with adverse outcome in independent CLL cohorts. Collectively, MethSig represents a novel inference framework for DNAme driver discovery to chart the role of aberrant DNAme in cancer.
MethSig provides a novel statistical framework for the analysis of DNA methylation changes in cancer, to specifically identify candidate DNA methylation driver genes of cancer progression and relapse, empowering the discovery of epigenetic mechanisms that enhance cancer cell fitness.
This article is highlighted in the In This Issue feature, p. 2113
DNA methylation (DNAme) is a central epigenetic modification of the human genome (1, 2). DNAme is also thought to be an important disease-defining feature in many cancers (3–6), pointing to the cancer's cell of origin and predictive of the outcome. Indeed, several tumor types harbor frequent mutations in genes that encode components of the methylation machinery (2).
DNAme changes in cancer have been described along two principal axes: global hypomethylation affecting retroviral elements and genome stability, and focal hypermethylation at promoters of tumor suppressor genes (TSG; refs. 1, 2). Promoter hypermethylation of TSGs has been surveyed across cancer in The Cancer Genome Atlas (TCGA) as well as other studies (1–3), and revealed that a plethora of cancer-related cellular pathways are disrupted by hypermethylation of TSG promoters, such as DNA repair (MLH1, RBBP8), cell cycle (CDKN2A, CDKN2B), p53 network (CDKN2A, TP73), apoptosis (WIF1, SFRP1), Ras signaling (RASSF1), Wnt signaling (SOX17), and tyrosine kinase cascades (SOCS3; refs. 7–9).
Although it is tempting to assume that all observed DNAme changes occur deterministically and drive the cancer phenotype, in vitro models and human cancers have shown that DNAme changes overwhelmingly follow a stochastic process (5, 6, 10, 11). While these changes are stochastic, they occur at different rates across the genome, correlated with features such as low gene expression and late DNA replication (5). Thus, stochastic DNAme changes in the growing malignant population result in a cancer methylome that displays locally disordered methylation and high intratumoral heterogeneity (5, 6). These data underscore the challenge of identifying candidate DNAme changes that are likely to be linked to the cancer phenotype among the highly abundant stochastic DNAme events across the genome, reminiscent of the challenge of distinguishing driver from passenger mutations in cancer.
However, unlike the field of cancer genomics where increasingly sophisticated tools have been developed to distinguish between driver and passenger mutations, accounting for confounding covariates (12, 13), inference tools in cancer epigenomics largely rely on uniform background models. Thus, widely used statistical methods produce hundreds or thousands of candidate promoter hypermethylation sites, likely overshadowing a much smaller number of DNAme changes that affect oncogenesis (referred to here as DNAme drivers).
To address this challenge, we developed a statistical inference framework accounting for varying stochastic hypermethylation rate across the genome and between patients: MethSig, analogous to leading approaches for cancer driver gene inference (12). MethSig estimates expected tumor promoter hypermethylation with an inference model that includes biological features known to affect the stochastic rate of DNAme changes (5). We applied MethSig to reduced representation bisulfite sequencing (RRBS; ref. 14) data across blood and solid tumor malignancies. Compared with benchmarked methods, MethSig delivers well-calibrated quantile–quantile (Q-Q) plots and more reproducible identification of DNAme drivers in independent cohorts. Importantly, MethSig achieved higher sensitivity and specificity in the inference of likely DNAme drivers compared with extant methods. Finally, the MethSig framework was extended to methylation array data (MethSig-array) and demonstrated the ability to identify candidate DNAme drivers, enriched in TSGs and associated with gene silencing, as well as disease outcomes. Thus, MethSig represents a novel statistical framework to infer DNAme drivers of cancer genesis and relapse, paving the way toward enhanced understanding of the role of epigenetic changes in cancer.
MethSig Infers Putative DNAme Drivers through the Application of an Optimized Background Model for Stochastic Hypermethylation
Promoter hypermethylation was measured using differentially hypermethylated cytosine ratio (DHcR), defined as the ratio of hypermethylated cytosines (HC) to the total number of CpGs profiled in promoters (Fig. 1A; Methods). We reasoned that a large number of high DHcR promoters may result from passenger hypermethylation due to the nonuniform distribution of stochastic hypermethylation rate across the genome (5). Extant inference tools relying on uniform background models will thus lead to a spuriously high number of significantly affected promoters. To illustrate this point, we compared simplified scenarios of constant versus varying HC rate across the genome, both sharing the same average promoter HC rate (Fig. 1B). This analysis demonstrated that when the HC rate varies across the genome, the uniform background assumption leads to many of the highly hypermethylated genes being falsely determined as significantly altered.
To overcome this challenge, we devised a model to estimate promoter-specific background hypermethylation rate (Fig. 1C; Supplementary Fig. S1A–S1E). We included covariates known to affect hypermethylation rate such as gene expression and replication time (Fig. 1D; Supplementary Fig. S1D and S1E), as well as promoter proportion of discordant reads (PDR; Fig. 1A; Supplementary Fig. S1B and S1C), a metric developed to characterize stochastic DNAme changes (5). Intuitively, loci and samples with high PDR (Fig. 1A, locus B) suggest lower reliability in DNAme driver identification compared with those with low PDR (Fig. 1A, locus A), akin to the role of background mutation rates in cancer driver gene inference.
To determine whether a single gene promoter is significantly hypermethylated in each sample, MethSig first generates the expected promoter hypermethylation (expected DHcR) based on a beta regression model and relevant covariates (Fig. 1C, step 1). Second, MethSig tests observed promoter hypermethylation (tumor DHcR) against the expected DHcR (Fig. 1C, step 2). Of note, the beta regression model was found to deliver good fits to tumor DHcR in most genes and patients (Supplementary Fig. S1F and S1G). Moreover, candidate DNAme drivers were not preferentially nominated by a small subgroup of patients (Supplementary Fig. S1H). Third, the hypermethylation signal is aggregated across the cohort as a stronger candidate DNAme driver is likely to affect a larger number of patients. This cross-patient aggregation procedure enables the estimation of hypermethylation enrichment at the cohort level (Fig. 1C, step 3).
To compare MethSig's performance to currently used methods, we applied three widely used methods to identify hypermethylation in cancer: t test, methylKit (15), and globalTest (ref. 16; see Methods and Supplementary Data for details). These methods were applied to prospective RRBS profiling of the CLL8 cohort (refs. 17, 18; Fig. 1E; Supplementary Table S1). Notably, a comparison of top candidates across methods showed relatively limited overlap, reinforcing the need to develop better statistical models to nominate candidate DNAme drivers (Supplementary Fig. S1I; Supplementary Tables S2 and S3).
In a well-calibrated statistical model, P values are uniformly distributed when the null hypothesis is true and all other assumptions are met. We thus evaluated the performance of MethSig and benchmarked methods through Q-Q plots (19), an established method to assess the uniformity of P value distribution in statistical genetics. As we anticipate only a small number of DNAme drivers (as compared with the much larger number of genome-wide stochastic changes), a well-calibrated Q-Q plot will mostly adhere to the diagonal, with few outliers with extreme P values. Benchmarked methods showed inflated Q-Q plots when applied to the CLL8 data set, which deviated from the expected line (dashed gray line) across the range of P values (Fig. 2A, first row). Nearly half of gene promoters were identified as candidate DNAme drivers, likely reflecting an underlying global phenomenon, such as an elevated rate of passenger DNAme alteration in chronic lymphocytic leukemia (CLL) compared with normal B cells, rendering the task of pinpointing candidate DNAme drivers, with biological and clinical significance, highly challenging. In contrast, MethSig exhibited a well-calibrated Q-Q plot, with a deviation factor that more closely approximated 1, and only a few candidate DNAme driver P values deviated from expected (Fig. 2A, first row).
Next, to test whether candidate DNAme driver nomination with MethSig is robust across data sets, we applied MethSig to an independent, previously published CLL RRBS data set (CLL-DFCI; Fig. 1E). Similarly, MethSig resulted in a well-calibrated Q-Q plot and a deviation factor closer to 1, compared with benchmarked methods (Fig. 2A, second row). To further test the generalizability of MethSig, we applied MethSig to available RRBS data sets of three additional tumor types (Fig. 1E; Supplementary Table S1). The performance of MethSig was maintained in a multiple myeloma cohort (MM-CNRS) and two solid tumor data sets (ductal carcinoma in situ, DCIS-MDACC, ref. 20; glioblastoma, GBM-MUV, ref. 21), resulting in well-calibrated Q-Q plots (Fig. 2A, third to fifth row). Here, too, benchmarked methods showed inflated Q-Q plots, suggesting that these methods are challenged in distinguishing oncogenic DNA hypermethylation from global DNAme changes.
We performed extensive model optimization to ensure the robustness of DNAme driver inference by MethSig, confirming that MethSig included informative covariates, parameters, and methodology (Supplementary Figs. S2A–S2H and S3A–S3G; see Supplementary Data for details). Notably, model optimization was also performed for benchmarked methods (e.g., using overdispersion correction option in methylKit); however, improvements of Q-Q plots were subtle (Supplementary Fig. S3H).
MethSig Provides Reproducible and Transcription-relevant Candidate DNAme drivers, Enriched in Genes Dysregulated across Cancer Types
Unlike passenger changes, DNAme drivers are anticipated to affect a large proportion of tumors and to be associated with silenced gene expression. Thus, we hypothesized that accurate inference of DNAme drivers can be assessed through reproducibility across independent patient cohorts and association with gene silencing.
Considering varied numbers of candidate DNAme drivers identified by different methods using an identical P value cutoff, we compared an equal number of top-ranking DNAme drivers to test the reproducibility of DNAme drivers nominated by different methods across the two CLL cohorts. MethSig resulted in a significantly higher overlap across the two cohorts, compared with benchmarked methods (Fig. 2B; Supplementary Fig. S4A).
Next, we tested whether DNAme drivers nominated by MethSig are more frequently linked to gene silencing compared with other methods (see Supplementary Data for details). Area under the receiver operating characteristic (AUROC) showed that MethSig achieved higher performance compared with benchmarked methods in identifying DNAme drivers associated with gene silencing (Fig. 2C; Supplementary Fig. S4B). Indeed, candidate DNAme drivers identified by MethSig were significantly more enriched in silenced genes compared with benchmarked methods or randomly selected genes (Fig. 2C). Similar findings were observed in the DCIS and GBM cohorts, where matched DNAme and RNA-sequencing (RNA-seq) data are available (Supplementary Fig. S4C and S4D).
Integrating data across the two CLL cohorts, MethSig nominated 189 candidate DNAme drivers out of 9,661 promoters captured by RRBS and with available input covariates (Supplementary Table S2; see Supplementary Fig. S5A and S5B; Supplementary Table S4; and Supplementary Data for additional analyses to rule out confounders including CpG density, B-cell subtype-specific epigenetic profiles, copy-number changes, and driver mutations). Samples where candidate DNAme drivers were found to be hypermethylated have a higher fraction of highly methylated promoters (DHcR > 0.75) compared with samples without hypermethylation (Supplementary Fig. S5C), suggesting high clonality level consistent with positive selection (see Supplementary Fig. S5D–S5F and Supplementary Data for further characterization of DNAme drivers). While known transcription factor binding motifs did not show enrichment in DNAme drivers, we observed significantly higher H3K27me3 signal at putative driver loci compared with nondriver loci, suggesting that MethSig candidates may in part conform to the model of promoting cancer development due to locking-in of repression by H3K27me3 (Supplementary Fig. S5G; Supplementary Data).
To interrogate their biological significance, we performed a pathway enrichment analysis of candidate DNAme drivers. Candidate CLL, multiple myeloma, DCIS, and GBM DNAme drivers were enriched in genes hypermethylated or silenced across tumor types, and associated with poor clinical outcome (ref. 22; Supplementary Table S5; Benjamini–Hochberg false discovery rate, BH-FDR, Q < 0.25). CLL and DCIS DNAme drivers were also enriched in genes downregulated by MYC and genes upregulated by p53 (22). Specifically, DCIS DNAme drivers were enriched in genes silenced in breast ductal carcinoma versus normal ductal breast cells (22), consistent with DNAme driver–mediated repression of corresponding genes.
Candidate CLL DNAme Drivers Include Established TSGs and Were Functionally Validated to Enhance Cancer Cell Fitness
Candidate CLL DNAme drivers included a well-established TSG, DUSP22, whose function as a TSG is silenced through promoter hypermethylation, as demonstrated previously in CLL (23). In addition to DUSP22, MethSig also identified other TSGs as putative CLL DNAme drivers such as RPRM and SASH1. RPRM is known to cooperate with p53 leading to cell-cycle arrest at G2 phase and has been reported to be hypermethylated or inactivated in carcinomas (24, 25). SASH1 encodes a scaffold protein involved in the TLR4 signaling pathway (26), which has been demonstrated to be a key signaling pathway in CLL (27).
To functionally validate candidate DNAme drivers identified with MethSig, given the limitations of demethylation agents and current dCas9-guided DNAme modification (Supplementary Data), we generated CRISPR/Cas9-mediated knockout (KO) to mimic gene silencing via promoter hypermethylation. We generated KO of three candidate DNAme drivers—DUSP22, RPRM, and SASH1 (see Supplementary Data for selection criteria). Of note, these candidates were suitable for functional validation given baseline gene expression and minimal promoter methylation in the HG3 cell line (Supplementary Fig. S6A–S6F).
After transduction with Cas9 and locus-specific targeting single guide RNA (sgRNA), HG3 cells were cultured with three leading CLL therapeutic agents: ibrutinib (a targeted BTK inhibitor), fludarabine (a key chemotherapy backbone in CLL chemoimmunotherapy regimens), and venetoclax (a BH3 mimetic; refs. 17, 28, 29). HG3 cells transduced with a nontargeting sgRNA (HG3-mock) were used as control. After 11 days (∼7 doubling times) of ibrutinib treatment, we observed higher fitness in cells with sgRNAs targeting all three candidate DNAme drivers (Fig. 3A). In contrast, only the DUSP22 KO cells showed higher proliferation after fludarabine treatment (Fig. 3B), and none of the KO led to greater proliferation with venetoclax, suggesting that DNAme drivers may have context-specific effects (Fig. 3C).
HG3 cells are known to show clonal diversity (30), which may affect bulk CRISPR/Cas9 KO. Furthermore, the DUSP22 locus is present in only one copy in HG3 cells due to a partial loss of the chromosome 6p (30), which may contribute to the greater effect in DUSP22 KO with fludarabine compared with RPRM (Fig. 3B). We therefore further generated stable KO HG3 clones of RPRM and DUSP22 through single-cell cloning (Fig. 3D; Methods; Supplementary Data). For RPRM, we identified a clone with biallelic frameshift-inducing indels and a second clone with monoallelic frameshift deletion (Fig. 3E). As DUSP22 locus is present in only one copy in the HG3 cell line, we generated two separate clones with complete gene KO by introducing frameshift indels in the remaining allele (Fig. 3F). A single cell–derived clone with a nontargeting sgRNA was used as a control (mock cell line). After culturing all clones without treatment for 7 days, we observed faster growth for the RPRM KO clones with a gene dose effect, compared with controls (Fig. 3G). Similarly, a KO clone for DUSP22 showed a significantly higher proliferation (Fig. 3G). These data are consistent with a fitness advantage in the absence of treatment, and the enrichment of these DNAme drivers in the previously untreated CLL8 cohort.
In agreement with our above results showing that candidate DNAme driver disruption confers resistance to treatment with leading CLL agents, we observed greater survival for the RPRM KO clones under ibrutinib and fludarabine with a gene dose effect (Fig. 3H; Supplementary Fig. S6G). Supporting the role of DUSP22 as a candidate DNAme driver, both KO clones for DUSP22 showed improved survival under ibrutinib and fludarabine treatment (Fig. 3I; Supplementary Fig. S6H). However, similar to the bulk transduction experiments, no fitness differences were observed with venetoclax treatment across KO clones (Supplementary Fig. S6G and S6H).
MethSig-nominated CLL DNAme Drivers Provide Independent Prognostic Information and Are Associated with Adverse Outcomes
We next sought to test the clinical significance of candidate DNAme drivers in the well-annotated CLL cohorts. Promoters whose hypermethylation is associated with failure-free survival (FFS) were defined as true positives (see Supplementary Data for details) while other promoters were defined as true negatives. In CLL8, MethSig resulted in the highest AUROC compared with benchmarked methods, and DNAme drivers identified by MethSig were enriched in genes associated with outcome compared with the benchmarked methods or randomly selected genes (Fig. 4A). We further validated this association with clinical outcome in the independent CLL-DFCI cohort (Supplementary Fig. S7A). Notably, MethSig also achieved higher AUROC compared with other methods when we combined two key features that were likely to be associated with DNAme drivers (i.e., either silenced by promoter hypermethylation or associated with FFS; Supplementary Fig. S7B). These results confirm that MethSig provides a nonincremental advance in the ability to identify likely DNAme drivers with high sensitivity and specificity.
Taking advantage of the large sample size in the CLL8 cohort, we sought to further triage the list of DNAme drivers by evaluating the clustering of methylated CpG positions. Intuitively, promoters with a nonrandom distribution of methylated CpGs are more likely to exert repression on corresponding genes and result in a substantial phenotypic impact (Fig. 4B). We used the maximum number of consecutive methylated CpGs to quantify the clustering degree of methylated CpGs (see Supplementary Fig. S7C and Supplementary Data for details). Indeed, we observed higher clustering of methylated positions in samples where the gene was predicted to be hypermethylated compared with other samples (Fig. 4B), and therefore applied an additional criterion of higher-level of clustering, decreasing the number of nominated CLL DNAme drivers to 122 (Supplementary Table S2).
To examine the prognostic value of DNAme drivers, we developed a clinical prediction score based on candidate DNAme drivers (n = 122). Elastic net regression (31) with a Cox proportional hazards model was used (see Supplementary Fig. S7D and Supplementary Data for model selection) to assign weights (coefficients) in terms of their contribution to the prediction of FFS to each candidate DNAme driver. To safeguard against overfitting and poor generalizability, CLL8 was used as the training set to select candidate DNAme driver coefficients, while CLL-DFCI was designated as an independent, test cohort not used in the training process.
Candidate DNAme drivers selected by the regression model included all three functionally validated TSGs, whose hypermethylation defined the subset with the least favorable prognosis (Fig. 4C). Higher risk score (greater than median) was significantly associated with shorter FFS in the training set [Fig. 4D, median FFS was 41.2 months in patients with high risk, and not reached in patients with low risk, HR 2.9; 95% confidence interval (CI), 2.1–4.0]. Notably, the model was also highly significant in distinguishing patients with high versus low risk of FFS in the test set (Fig. 4E; Supplementary Fig. S7E), and a regression model including established CLL risk indicators demonstrated that DNAme drivers contribute to adverse clinical outcome independently of previously established risk factors (Fig. 4F; Supplementary Fig. S7F–S7H; Supplementary Data).
MethSig Identifies Relapse-specific DNAme Drivers in CLL
Our data demonstrated that MethSig is an effective tool to nominate cancer DNAme drivers through the comparison of primary malignant (T1) versus controls. We sought to extend the application of MethSig to identify DNAme drivers of relapsed disease after fludarabine-based chemotherapy through the comparison of relapse (T2) versus control samples (Fig. 1E; CLL8). Notably, the application of MethSig within the context of relapsed CLL resulted in an equally well-calibrated Q-Q plot (Fig. 4G), consistent with its ability to identify the infrequent DNAme changes that likely contribute to the relapse phenotype.
We identified T2-specific (n = 32), T1 and T2 shared (n = 88), and T1-specific DNAme drivers (n = 101; Fig. 4H and I; Supplementary Table S2). In addition to previously observed DNAme drivers (e.g., T1 and T2 shared, DUSP22, SASH1), T2-specific DNAme drivers involve additional genes with potential tumor suppressor function, such as G0S2. G0S2 can promote apoptosis through BCL2, the therapeutic target of the BH3 mimetic venetoclax in CLL (29, 32). A pathway enrichment analysis of T2-specific DNAme drivers revealed enrichment in TP53 targets and DNA damage pathway (22), whereas T1- and T2-shared or T1-specific DNAme drivers were not enriched in these pathways (Fig. 4I; Supplementary Table S5). The enrichment in TP53 targets and DNA damage pathway of T2-specific DNAme drivers indicates that CLL relapse after chemotherapy may follow an alternative path compared with CLL progression in the absence of therapy, offering novel insights for therapeutic strategies to address drug-resistant or relapsed cancer.
MethSig-array Infers Candidate DNAme Drivers with Methylation Arrays
Considering the wide availability of methylation array data, we designed MethSig-array under the same statistical framework proposed by MethSig. Of note, promoter PDR cannot be estimated by array data, which does not provide read-level methylation information, and, as shown in the covariate analysis, promoter PDR provides an important contribution to the model (Supplementary Fig. S3B).
Nonetheless, we applied MethSig-array to Infinium HumanMethylation450 arrays of 18 tumor types in the TCGA Pan-Cancer analysis project (ref. 33; Supplementary Table S6). As anticipated, the deviation factors of Q-Q plots derived from MethSig-array were closer to 1, compared with higher deviation factors of benchmarked methods (Fig. 5A). To further evaluate the performance of MethSig-array, AUROC was used to assess the sensitivity and specificity in the inference of likely DNAme drivers, which were defined following three key readouts: association with gene silencing, association with disease outcome, and enrichment with TSGs using different published catalogs (Supplementary Data). MethSig-array achieved higher AUROC compared with benchmarked methods in the inference of likely DNAme drivers associated with gene silencing (Fig. 5B) and clinical outcome (Fig. 5C). MethSig-array also resulted in highest AUROC compared with benchmarked methods in the inference of TSGs (Supplementary Fig. S7I; OncoKB, ref. 34; or the TCGA cancer driver study, ref. 35). For example, SOX17 was identified as a DNAme driver in 13 different tumor types (Supplementary Data), which encodes a transcription factor involved in embryonic development and cell fate (9). Hypermethylation and downregulation of SOX17 have been described in multiple cancer types, which implies the broad tumor suppression function of SOX17 gene (9). Another important TSG is RASSF1, which was identified as a DNAme driver in six tumor types (Supplementary Data). RASSF1 is a microtubule-associated and multitasking scaffold protein communicating with the RAS pathway, estrogen receptor signaling, and Hippo pathway (36). RASSF1 methylation is proposed as a candidate maker in many cancer types (36). Other identified important TSGs included genes in a plethora of cancer-related cellular pathways, such as DNA repair (MLH1, RBBP8), apoptosis (WIF1, SFRP1), and tyrosine kinase cascades (SOCS3; refs. 7, 8). Collectively, these data confirm that MethSig can accurately infer likely DNAme drivers across cancer with both array and next-generation sequencing (NGS)–based methylation assays.
Aberrant gene function due to acquired epigenetic abnormalities has been highlighted as a key feature of cancer over the last decade, implicated in cancer initiation, progression, and treatment resistance (1, 2). Although the causal role of DNAme in cancer remains to be conclusively determined (37), genome-wide DNAme analyses have provided comprehensive surveys of the cancer epigenome and tumor-associated DNAme changes, and have proposed that these changes fuel the malignant process through TSG silencing and other mechanisms (1, 2).
However, in the context of steadily growing DNAme-sequencing data sets, a major challenge remains: to identify the DNAme changes involved in tumor progression among the abundant stochastic DNAme changes that occur in cancer cells. This challenge is reminiscent of the challenge of distinguishing driver from passenger mutation in cancer exome or genome data. While for the latter challenge, progress has been achieved through increasingly sophisticated inference tools that model the varying background mutation rate across the genome (12), inference tools in cancer epigenomics largely rely on uniform background models (38). Given the recent observation that stochastic DNAme varies widely in different genomic regions (5), statistical models relying on a uniform background assumption are anticipated to lead to spuriously high numbers of significantly affected regions.
Drawing on lessons learned in cancer genomics, we posited that robust nomination of oncogenic DNAme changes requires a rethinking of the statistical inference process to enable the differentiation of driver promoter hypermethylation changes (DNAme drivers) from the far larger number of stochastic DNAme changes without biological consequences (passenger DNAme changes). To address this challenge, we developed a statistical inference framework accounting for varying stochastic hypermethylation rate across the genome and between samples—MethSig. The model provides the expected promoter DNAme changes between tumor and control samples, allowing the identification of loci where the observed hypermethylation significantly exceeds expectation, potentially reflecting positive selection of fitness-enhancing candidate DNAme drivers.
We applied MethSig to methylation-sequencing data of two CLL cohorts, including 304 CLLs from a prospective clinical trial, as well as to other malignancies with available DNAme data (multiple myeloma, DCIS, and GBM), and benchmarked against state-of-the-art methods. Compared with benchmarked methods, MethSig resulted in well-calibrated Q-Q plots, higher reproducibility in DNAme driver inference across independent cohorts, and increased sensitivity/specificity in the inference of likely DNAme drivers. These observations confirm that MethSig allows separation of specific cancer-related DNAme drivers, which cause gene downregulation and phenotypic changes associated with tumoral progression, from stochastic passenger DNAme changes. Notably, the performance of MethSig was maintained across both hematologic malignancies (CLL and multiple myeloma) and solid tumors (DCIS and GBM), suggesting broad applicability for creating catalogs of candidate DNAme drivers of cancer genesis and relapse. Moreover, while MethSig was extended to array data and provided a nonincremental improvement in DNAme driver inference, we anticipate that future shift toward NGS data will leverage the even higher performance of MethSig with read-level data.
Our data showed that MethSig can also account for broad phenomena that alter DNAme profiles in identifying gene-specific DNAme drivers. For example, DNAme of the cell-of-origin represents one of the strongest sources of variation in the cancer epigenome (4). Indeed, in CLL, DNAme has been shown to strongly encode normal B-cell epigenetic reprogramming during differentiation, allowing high-resolution inference of the differentiation state of the initially transformed B cell (4). It is therefore notable that MethSig candidate DNAme drivers were not significantly enriched in the most variable methylated regions identified between naïve and class-switched memory B cells (Supplementary Data).
CLL candidate DNAme drivers included TSGs inactivated through hypermethylation, such as DUSP22, RPRM, and SASH1, which may play important roles in the initiation and relapse of CLL. To functionally validate these candidate DNAme drivers in CLL cells, we generated single or double allele frameshift KO. While transformed cells showed superior fitness in the absence of drug selection and with therapy, DNAme drivers showed context-specific effects, which may underlie some of the heterogeneity in CLL clinical course (5). These data also show that venetoclax therapy may uniquely overcome these mechanisms, providing rationale for future DNAme driver–guided trials. Of note, our data demonstrated improved risk stratifications based on candidate DNAme drivers, independent of known prognostic factors in CLL, suggesting that DNAme drivers contribute to adverse clinical outcome.
While we believe that this work presents a transformative advance in identifying DNAme drivers with higher sensitivity and specificity, further work will be needed to improve performance and reduce false-positive candidates. This may be achieved through expanded data sets and future discovery of additional informative covariates, following the example of genetic driver inference where successive versions resulted in a continuous improvement in performance and reduction in false positives (12, 13). Given the limitations of RRBS, including poor coverage of distal regulatory elements or other regions that are CpG poor, the future exploration of DNAme changes in other genomic areas will be greatly empowered by larger whole-genome DNAme-sequencing data sets. Finally, hypomethylation may also play important roles in cancer genesis by affecting genome stability. Further efforts will be needed to enable statistical inference of those epigenetic events.
Collectively, our data support a novel framework for the analysis of DNAme changes in cancer to specifically identify DNAme drivers of disease progression and relapse, empowering the discovery of candidate epigenetic mechanisms that may enhance cancer cell fitness. This work addresses a central gap between cancer epigenetics and genetics, where such tools have had a transformative impact in precision oncology and cancer gene discovery. We envision that inference tools such as MethSig, coupled with novel DNAme-sequencing modalities (39) and emerging tools for epigenetic editing (40), may herald a new era in cancer epigenomics in which large cohort studies will provide precision identification of oncogenic DNAme drivers for improved patient stratification and therapeutic targeting.
For the CLL8 cohort, blood was obtained from previously untreated patients enrolled in a prospective, randomized, open-label CLL8 trial (17, 18) before the first cycle of treatment. Written informed consent for genomic sequencing of patient samples was obtained prior to the initiation of sequencing studies. For the MM-CNRS cohort, bone marrow of patients presenting with previously untreated multiple myeloma (n = 24) or at relapse (n = 20) was obtained after patients' written informed consent in accordance with the Institutional Review Board and the Montpellier University Hospital Centre for Biological Resources (DC-2008–417). Genomic DNA was extracted from CLL and multiple myeloma cells.
RRBS libraries were generated by digesting genomic DNA with MspI to enrich for CpG-rich fragments, and then ligated to barcoded TruSeq adapters (Illumina) to allow immediate subsequent pooling. It was followed by bisulfite conversion and PCR, as described previously (14). Libraries were sequenced and aligned to the bisulfite-converted hg19 reference genome using Bismark v0.15.0 (RRID: SCR_005604; ref. 41).
Promoter (defined as ± 2 kb windows centered on RefSeq transcription start site) hypermethylation was measured using DHcR, defined as the ratio of HCs to the total number of CpGs profiled in promoter. HCs of each sample were defined as CpGs at which DNAme is statistically higher than the control (FDR = 20%, χ2 test; ref. 6). Only CpGs with read depth greater than 10 were included in the analysis. DHcR of each normal sample was calculated in the same way as for the tumor samples, testing against all normal sample controls. This was followed by averaging DHcR of all the normal samples as the normal DHcR.
If all the CpGs on a specific read are methylated or unmethylated, the read is classified as concordant. Otherwise, it is classified as discordant. At each CpG, the PDR is equal to the number of discordant reads divided by the total number of reads that cover that location. Promoter PDR is given by averaging the values of individual CpGs, as calculated for all CpGs within the promoter of interest that are covered by a minimum of 10 reads that contain at least 4 CpGs. The normal PDR was calculated by averaging PDR of all the normal samples.
The superscripts n, t, e, are shorthand for normal, tumor, and expected. The subscripts i and j represent gene i and sample j. The model was processed as following steps:
(i) Estimate expected hypermethylation of tumor samples (DHcRe): The independent variable matrix (X) has number of genes times the number of patients rows (Eq. A). The beta regression model was implemented by R package betareg (ref. 42; Eq. B). Next, predicted distribution of DHcRt (used as distribution of DHcRe in the following analysis) was estimated using and derived from the above beta regression model (Eq. C).
(ii) Evaluate whether DHcRt (observed) is significantly higher than DHcRe (expected): DHcRt was tested against the distribution of DHcRe. The patient-specific P value indicates the probability that observed promoter hypermethylation is significantly higher than expected (Eq. D). Only genes whose promoter hypermethylation significantly exceeds expectation will be assigned as patient-specific DNAme drivers (P < 0.05).
(iii) Determine whether promoter hypermethylation is overrepresented in patients (DNAme driver): Wilkinson P value combination method was used to combine P values from different patients to identify those frequently recurring DNAme drivers (43). To eliminate the effect of cohort size on P value combination results, MethSig randomly sampled equal number of patients (K = 10) iteratively (S = 100) and used lower quartile of combined P values to identify DNAme drivers. Wilkinson P value combination was performed by R package metap (https://cran.r-project.org/web/packages/metap/).
In the evaluation of benchmarked methods, t.test of R 3.3.2, methylKit 1.0.0 (RRID: SCR_005177), and globalTest 5.28.0 (RRID: SCR_001256) were used.
Pathway Enrichment Analysis
Pathway enrichment analysis was limited to the chemical and genetic perturbations of the C2 gene set collection (22), which includes gene sets that are more specific to cancer processes in different cancer types with or without perturbation. DNAme drivers were tested against all MethSig-inferred nondrivers and only pathways with at least 10 inferred genes were included. A hypergeometric test was used to measure the enrichment of DNAme drivers in each gene set, followed by a BH-FDR procedure.
HG3 (DSMZ #ACC-765, RRID: CVCL_Y547), PGA1 (DSMZ #ACC-766, RRID: CVCL_Y545), and MEC1 (DSMZ #ACC-497, RRID: CVCL_1870) cells were provided by Leibniz Institute DSMZ in August 2018. HEK293T cells (ATCC #CRL-3216, RRID: CVCL_0063) were provided by the ATCC in January 2017. Cells were routinely tested for Mycoplasma using the MycoAlert Mycoplasma Detection Kit (Lonza #LT07-318). All cell lines used for described experiments came from early frozen batches between 3 and 8 passages after cell receipt. HG3, PGA1, and MEC1 cells were used for qRT-PCR (Supplementary Data). HEK293T cells were used for the production of lentivirus and transduction (Supplementary Data). HG3 cells were used for all the other described experiments.
CRISPR/Cas9 Design and Cloning
sgRNAs were designed using CHOPCHOP (RRID: SCR_015723; ref. 44), to minimize in silico predicted off-target activity (≤1), target the first exon of the genes of interest, and have a good predicted efficiency (>56) based on the PAM sequence and the 3′ nucleotide (45). Two sgRNAs have been designed for each targeted locus and the empirically defined most efficient sgRNA was then used for the CRISPR/Cas9 experiments. The sgRNAs were cloned into lentiCRISPRv2 puro (Addgene #98290) as described previously (46).
Generation of Cas9-expressing KO HG3 Clones
To create single and double allele KO HG3 clones, we have transduced cells with CRISPR/Cas9 and a gene-specific sgRNA, performed a 10-day puromycin (0.5 μg/mL) selection, confirmed Cas9 efficiency using T7E1 – EnGen Mutation Detection Kit (NEB), created single-cell colony by cell sorting, and selected single-cell clones showing single or double allele KO (indels) by targeting the predicted CRISPR cutting site through NGS.
Risk Model of Patients with CLL
The regression model was implemented by R package glmnet (RRID: SCR_015505; refs. 47, 48). When evaluating the performance of the model in the training and test set, CLL cases were divided into two subgroups based on their predicted risk scores (patients with high vs. low risk, median risk score of the cohort was used as the cutoff) and the FFS difference between groups was evaluated using log-rank test.
Promoter DHcR and PDR cannot be estimated by array data, which does not provide read-level methylation information. MethSig-array was designed by using average promoter methylation instead of DHcR in the beta regression model and leaving out promoter PDR, under the same statistical framework proposed by MethSig.
Statistical analysis was performed with R version 3.3.2 (https://www.R-project.org/). All P values were two-sided and considered significant at the 0.05 level unless otherwise noted.
RRBS data of CLL8 are available via GEO accession number GSE143673. RRBS data of CLL-DFCI are available via dbGap accession number phs000435.v3.p1. RRBS data of DCIS-MDACC are available via GEO accession number GSE69994. RRBS data of GBM-MUV are available via EGA accession number EGAS00001002538.
The MethSig pipeline is available on GitHub at https://github.com/HengPan2007/MethSig.
E. Tausch reports grants and personal fees from Roche; grants, personal fees, and non-financial support from AbbVie; and personal fees and non-financial support from Janssen outside the submitted work. A.M. Fink reports grants from Celgene; personal fees from AbbVie, and personal fees from Janssen outside the submitted work. K. Fischer reports other support from Roche and personal fees from Roche during the conduct of the study; and personal fees from Abbvie outside the submitted work. C.J. Wu reports other support from Pharmacyclics outside the submitted work. O. Elemento reports personal fees and other support from Volastra Therapeutics, OneThree Biotech, and Champions Oncology, and other support from Freenome and Owkin during the conduct of the study. D.A. Landau reports grants from NIH and grants from LLS during the conduct of the study; non-financial support from Illumina, personal fees from Mission Bio, personal fees from C2i Genomics, and grants from BMS outside the submitted work; and has a patent related to this work pending. No disclosures were reported by the other authors.
H. Pan: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. L. Renaud: Validation, writing–original draft, writing–review and editing. R. Chaligne: Validation, writing–review and editing. J. Bloehdorn: Resources, data curation, writing–original draft, writing–review and editing. E. Tausch: Resources, data curation. D. Mertens: Resources, data curation. A.M. Fink: Resources, data curation. K. Fischer: Resources, data curation.C. Zhang: Methodology. D. Betel: Methodology. A. Gnirke: Resources, data curation. M. Imielinski: Methodology, writing–original draft, writing–review and editing. J. Moreaux: Resources, data curation, writing–review and editing. M. Hallek: Resources, data curation, writing–original draft, writing–review and editing. A. Meissner: Writing–original draft, writing–review and editing.S. Stilgenbauer: Resources, data curation, writing–original draft, writing–review and editing. C.J. Wu: Resources, data curation, writing–original draft, writing–review and editing. O. Elemento: Conceptualization, software, formal analysis, supervision, methodology, writing–original draft, writing–review and editing. D.A. Landau: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.
R. Chaligne is supported by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 750345. D.A. Landau is supported by the Burroughs Wellcome Fund Career Award for Medical Scientists, Vallee Scholar Award, Pershing Square Sohn Prize for Young Investigators in Cancer Research, the Sontag Foundation Distinguished Scientist award, and the NIH Director's New Innovator Award (DP2-CA239065). J. Moreaux is supported by Institut Universitaire de France (IUF), by grants from INCa (Institut National du Cancer; PLBIO2018-160 PIT-MM), ANR (PLASMADIFF-3D), ANR (TIE-Skip; 2017-CE15-0024-01), and SIRIC Montpellier Cancer INCa_Inserm_DGOS_12553. S. Stilgenbauer is supported by the DFG, SFB1074 subprojects B1 and B2. This work was supported by the NCI (1R01CA229902), the Leukemia & Lymphoma Society, Quest for Cures program, and LLS-SCOR 7012-16. The authors thank all patients and their physicians for trial participation and donation of samples, the DCLLSG; Sabrina Schrell and Christina Galler for their excellent technical assistance on CLL8 sample workup.