BOREALIS: an R/Bioconductor package to detect outlier methylation from bisulfite sequencing data

Introduction

Multiple inborn diseases resulting from disorders of genomic methylation are well characterized. A growing body of literature reports associations between DNA methylation and conditions including Parkinson’s Disease (Chuang, et al., 2017) and methylation-based studies have also suggested the causative or correlative role of aberrant methylation in diverse rare inherited conditions (Guastafierro, et al., 2017; Sharp, et al., 2017; Sobreira, et al., 2017).

Methods for profiling DNA methylation traditionally focused on wide genomic regions, particularly CpG islands. While microarray methylation profiling continues to be utilized, next-generation sequencing has enabled analysis at read count-level resolution and greatly increased numbers of genomic CpG sites (Wu, et al., 2015). High-resolution methodologies have enabled complex pictures of methylation to emerge, including association of multiple specific-site or even single-site methylation with developmental processes or disease (Bui, et al., 2012; Choi, et al., 2018; Claus, et al., 2012; Fürst, et al., 2012; Hashimoto, et al., 2013; Nile, et al., 2008; Pogribny, et al., 2000; Scantamburlo, et al., 2017; Sohn, et al., 2010; Takahashi, et al., 2017).

DNA sequencing results in the diagnosis of up to 40% of genetic disease cases previously unsolved using standard clinical testing (Sawyer, et al., 2016). RNA-Seq has been investigated and has shown benefit in complementing DNA testing (Cummings, et al., 2017). With the growing evidence linking methylation to disease, the application of methylation-based sequencing for undiagnosed cases of rare inherited disease is a logical progression from current testing paradigms. Expanded methylation profiling will offer the ability to detect diagnostic signals unique to the epigenome, undetectable in DNA and RNA due to lack of measurable manifestation in those materials, or due to shortcomings in current technologies or analytical approaches.

An ideal method to detect deviant methylation should offer the ability to profile at single CpG sites while enabling flexibility to consolidate calls across regions. In the context of genetic disease, we can assume that disease-associated methylation aberrations in an individual will show significant differences from individuals with unrelated phenotypes. A similar rationale was successfully used by us and others in outlier-based RNA analysis (Jenkinson, et al., 2020). However, existing solutions for the detection of differentially methylated CpG sites from bisulfite sequencing focus on traditional group vs group or multi-group experimental designs (Wreczycka, et al., 2017) and are therefore not suited to genetic disease or other outlier-based analyses (Figure 1).

Figure 1. Conceptual differences between traditional case vs control analysis and the BOREALIS approach.

While traditional approaches to differential methylation analysis are based on group case vs. control analysis, BOREALIS utilizes a one vs. many outlier approach whereby individual(s) are compared to a cohort. This approach is especially useful when multiple similar cases are difficult to identify, as is the case in rare genetic disease studies. By comparing every affected individual to a cohort of heterogeneous individuals, outlier methylation can be identified at individual CpG sites for all members of the cohort, without the requirement for multiple similar cases.

Methods

At a given CpG site, we assume we have data in the form of methylated counts x_i and total read counts n_i for individuals i = 1, …, I in a cohort of size I. If every individual in the population had the exact same probability p of methylation at this site, i.e., p₁ = p₂ = … = p_i = p where p_i is the (true-but-unknown) probability of methylation for the i^th individual at this site, then the methylated counts x_i would be binomially distributed with parameters p and n_i. However, we expect varying degrees of sample-to-sample variability in the probability of methylation at a given site even in a healthy cohort. Therefore, it is more biologically accurate to assume that p_i for i=1, …, I have been sampled from a distribution over the unit interval. A common and mathematically convenient choice for this distribution is a beta distribution with parameters α and β. Thus the observed number of methylated reads x_i for the i^th individual can be viewed as being generated from a two-step process whereby the probability of methylation is selected p_i ~ Beta (α,β) and given this probability the number of methylated reads is binomially distributed x_i|p_i ~ Binomial (p_i,n_i). Viewed in this way, we say the methylated reads are beta-binomially distributed x_i ~ Beta-Binomial (n_i,α,β) with parameters α and β, and popular packages for differential DNA-methylation detection such as Dispersion Shrinkage for Sequencing (DSS) (Feng and Wu, 2019) use this same distribution for methylated counts.

What BOREALIS does differently from traditional tools such as DSS (Park and Wu, 2016) is that it builds its statistical model explicitly for the purpose of outlier detection compared to a cohort, which requires alternative statistical framing and considerations as compared to group-versus-group analyses (Jenkinson, et al., 2020). Specifically, at each CpG site, BOREALIS takes the input data {(x_i,n_i): i = 1, …, I} and estimates the population-level parameters α and β using gamlss library (https://www.gamlss.com) (Feng and Wu, 2019). We implement Laplace Smoothing on the counts (i.e. we use as counts $\widetilde{x_i}$= x_i+1 and $\widetilde{n_i}$= n_i+2) as a regularization step to help deal with any samples with low counts. From these estimated α and β parameters, we (for the i^th sample) compute the left-sided p-value by looking at the probability that a value of x_i or fewer methylated reads were generated from a Beta-Binomial (n_i,α,β), and likewise a right-sided p-value would evaluate the probability that a value of x_i or greater methylated reads came from this distribution. We implement this probability calculation using the pBB function of gamlss. The two-sided p-value is computed as two times the lesser of these one-sided p-values.

To validate performance of the BOREALIS method, we performed Monte Carlo simulations of an outlier sample in cohorts of varying sizes sequenced at varying depths of coverage. Namely, after selecting a cohort size I and an average depth of coverage D we conducted 10,000 Monte Carlo simulations wherein each sample i in the cohort has sequencing depth d_i drawn from a Poisson distribution with mean D. The number of methylated reads in cohort sample i is drawn from a Beta-Binomial distribution with mean 0.8 and dispersion 0.1, and then these simulated cohort data are fit using the BOREALIS model. An outlier sample is then simulated with sequencing depth d drawn from a Poisson distribution with mean D and number of methylated reads given by a Binomial distribution with mean 0.3. BOREALIS is then used to compute a p-value thresholded at level 0.05, and the power is given by the fraction of the 10,000 simulations that correctly reject the null hypothesis. Full code to replicate the power analysis is provided as described in the Data Availability section.

Results

The results of the power analysis are plotted in Figure 2 demonstrating that BOREALIS can accurately detect outlier methylation events in modest cohort sizes and sequencing depths. BOREALIS can successfully identify outlier methylation with high statistical power utilizing a wide range of sample numbers (3-100+) and read depths (< 10 – 100s). BOREALIS is also tolerant of multiple identical outliers provided sufficient sequencing depth and cohort size as shown in Figure 3. The method supports multithreaded computation, as well as splitting across chromosomes to facilitate parallelism across compute nodes in a cluster environment. To provide users with the ability to visually review the methylation distributions underlying any call made by BOREALIS, we provide a built-in plotting function whose outputs are illustrated in Figure 4.

Figure 2. BOREALIS power analysis and single site methylation profile output by BOREALIS.

Graphical summarization of Monte Carlo simulations of an outlier sample in a cohort of varying sizes and depths of sequencing coverage. Ten thousand simulations were performed for each set of experimental conditions whereby random sampling of sequencing depth for each sample and the proportion methylated reads was performed. BOREALIS built its beta-binomial model for each simulated cohort and an outlier sample was simulated with BOREALIS used to compute a p-value. Power estimation is based on the simulations correctly rejecting the null hypothesis at a level p < = 0.05. The parameters mu (mean methylation fraction at a given site), sigma (variability in methylation fraction at a given site) and muAb (deviation from the mean methylation level in the outlier sample) are fixed for the purposes of the simulation shown.

Figure 3. BOREALIS power analysis containing multiple individuals with identical outlier methylation events.

Graphical summarization of Monte Carlo simulations of (A) two outlier samples and (B) three outlier samples in a cohort of varying sizes and depths of sequencing coverage. While rare disease studies generally aim to identify single outlier individuals in a cohort, this demonstrates the ability of the BOREALIS approach to identify multiple, identical outliers in the presence of sufficient read coverage and cohort size. One thousand simulations were performed for each set of experimental conditions whereby random sampling of sequencing depth for each sample and the proportion methylated reads

was performed. BOREALIS built its beta-binomial model for each simulated cohort and an outlier sample was simulated with BOREALIS used to compute a p-value.

Figure 4. BOREALIS package functionality visually represents the methylation profile at any given CpG site, as compared to the cohort.

Here a single site within the LTB4R gene promoter is shown for Patient 72, from the BOREALIS Bioconductor package’s included test data. Full details on how to perform BOREALIS analysis and generate similar figures are detailed in the BOREALIS vignette, included with the package.

Conclusions

BOREALIS is a novel R/Bioconductor package that addresses an unmet need in bisulfite sequencing-based genetic disease studies. The method is suited for single or multi-site downstream analysis and can successfully identify outlier methylation with high statistical power, providing a new avenue of exploration in the quest for increased diagnostic rates in genetic disease patients. It is readily available and easily implemented, enabling seamless integration with other common pipelines and tools.

Author contributions

GRO’s contributions include Data Curation, Formal Analysis, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing – Original Draft Preparation, Writing – Review & Editing. WGJ was involved with Conceptualization, Data Curation, Formal Analysis, Methodology, Project Administration, Software, Supervision, Validation, Visualization, Writing – Review & Editing. RJO and LESR participated in study design, helped with analysis and edited the manuscript. EWK participated in study design, secured funding and edited the manuscript.

Software availability