A novel statistical method for modeling covariate effects in bisulfite sequencing derived measures of DNA methylation

Abstract

Identifying disease-associated changes in DNA methylation can help us gain a better understanding of disease etiology. Bisulfite sequencing allows the generation of high-throughput methylation profiles at single-base resolution of DNA. However, optimally modeling and analyzing these sparse and discrete sequencing data is still very challenging due to variable read depth, missing data patterns, long-range correlations, data errors, and confounding from cell type mixtures. We propose a regression-based hierarchical model that allows covariate effects to vary smoothly along genomic positions and we have built a specialized EM algorithm, which explicitly allows for experimental errors and cell type mixtures, to make inference about smooth covariate effects in the model. Simulations show that the proposed method provides accurate estimates of covariate effects and captures the major underlying methylation patterns with excellent power. We also apply our method to analyze data from rheumatoid arthritis patients and controls. The method has been implemented in R package SOMNiBUS.

Keywords: EM algorithm; differentially methylated region; generalized additive model; next-generation sequencing; penalized regression splines.

FIGURE 1

(A), The estimates (solid lines) and 95% pointwise confidence intervals (dashed lines) of the intercept, the smooth effect of RA and cell type (T cells versus monocytes) on methylation levels. (B), The predicted methylation levels in the logit scale (left) and proportion scale (right) for the four groups of samples with different disease and cell type status. The region‐based P‐values for the effect of RA status and T cell type are calculated as $1.11E-16$ and $6.37E-218$, respectively

FIGURE 2

The 14 simulation settings of methylation parameters $\pi (t)$ in Scenario 2. Methylation parameters for samples with $Z=1$ (dotted‐dashed black curve) are fixed across settings, whereas the methylation parameters for samples from group $Z=0$ (solid gray lines) vary across simulations corresponding to different degrees of closeness between methylation patterns in the two groups

FIGURE 3

Estimates of smooth covariate effects (gray) over the 100 simulations in Scenario 1, using SOMNiBUS. The black curves are the true functional parameters used to generate the data. Data with sample size $N=40$ were generated with error

FIGURE 4

Coverage probability of confidence intervals over 1000 simulations under different sample sizes ($N=40,100,150,400$). Data were generated with error, under simulation Scenario 1

FIGURE 5

Quantile‐Quantile (Q‐Q) plots of the region‐based P‐values for the null covariate Z ₃, obtained from the six methods, over 1000 simulations. Data were generated without error with a range of sample sizes ($N=40,100,150,400$), under simulation Scenario 1. Here, the Expected P‐values are uniformly distributed numbers, equal to $=(1/1001,2/1001,\dots ,1000/1001)$.

FIGURE 6

Powers to detect DMRs using the six methods for the 14 simulation settings in Scenario 2 under different levels of maximum deviation between ${\pi }_0(t)$ and ${\pi }_1(t)$, calculated over 100 simulations. (Sample size $N=100$).