MBD-seq as a cost-effective approach for methylome-wide association studies: demonstration in 1500 case--control samples

Abstract

Aim: We studied the use of methyl-CpG binding domain (MBD) protein-enriched genome sequencing (MBD-seq) as a cost-effective screening tool for methylome-wide association studies (MWAS).

Materials & methods: Because MBD-seq has not yet been applied on a large scale, we first developed and tested a pipeline for data processing using 1500 schizophrenia cases and controls plus 75 technical replicates with an average of 68 million reads per sample. This involved the use of technical replicates to optimize quality control for multi- and duplicate-reads, an in silico experiment to identify CpGs in loci with alignment problems, CpG coverage calculations based on multiparametric estimates of the fragment size distribution, a two-stage adaptive algorithm to combine data from correlated adjacent CpG sites, principal component analyses to control for confounders and new software tailored to handle the large data set.

Results: We replicated MWAS findings in independent samples using a different technology that provided single base resolution. In an MWAS of age-related methylation changes, one of our top findings was a previously reported robust association involving GRIA2. Our results also suggested that owing to the many confounding effects, a considerable challenge in MWAS is to identify those effects that are informative about disease processes.

Conclusion: This study showed the potential of MBD-seq as a cost-effective tool in large-scale disease studies.

Figure 1. Overview of the data processing pipeline and key results at each step

FDR: False-discovery rate; MBD-seq: Methyl-CpG binding domain sequencing; MWAS: Methylome-wide association studies; PC: Principal component; PCA: Principal component analysis; QC: Quality control; QV: Quality value.

Figure 2

Correlations calculated for 73 duplicates after different quality control procedures for duplicate- and multi-reads.

Figure 3. CpG density versus ‘raw’ coverage

For CpG density a simple ‘coupling’ factor was calculated by counting the number of CpGs within ±100 bp. The vertical lines indicate the relative frequencies of CpGs with that specific density.

Figure 4

Distribution of in silico coverage for CpGs in loci overlapping and not overlapping with RepeatMasker.

Figure 5. Size of stage 2 blocks across biological features

The mean, median, SD and 99th percentile for the stage 2 block size, in base pairs, is given for all blocks included and for blocks overlapping with CpG islands, CpG shores, regions marker by Rep. Mask, genes, exons, introns, 5′-UTRs, 3′-UTRs, regions within 8 kb upstream of transcriptional start sites corresponding to potential promoter regions (promoter) and conserved regions (cons). The percentage of blocks overlapping with the biological features are given in the legend on the x-axis. Note that a single block can overlap with multiple biological features. Rep. Mask: RepeatMasker; SD: Standard deviation.