QSEA-modelling of genome-wide DNA methylation from sequencing enrichment experiments

Abstract

Genome-wide enrichment of methylated DNA followed by sequencing (MeDIP-seq) offers a reasonable compromise between experimental costs and genomic coverage. However, the computational analysis of these experiments is complex, and quantification of the enrichment signals in terms of absolute levels of methylation requires specific transformation. In this work, we present QSEA, Quantitative Sequence Enrichment Analysis, a comprehensive workflow for the modelling and subsequent quantification of MeDIP-seq data. As the central part of the workflow we have developed a Bayesian statistical model that transforms the enrichment read counts to absolute levels of methylation and, thus, enhances interpretability and facilitates comparison with other methylation assays. We suggest several calibration strategies for the critical parameters of the model, either using additional data or fairly general assumptions. By comparing the results with bisulfite sequencing (BS) validation data, we show the improvement of QSEA over existing methods. Additionally, we generated a clinically relevant benchmark data set consisting of methylation enrichment experiments (MeDIP-seq), BS-based validation experiments (Methyl-seq) as well as gene expression experiments (RNA-seq) derived from non-small cell lung cancer patients, and show that the workflow retrieves well-known lung tumour methylation markers that are causative for gene expression changes, demonstrating the applicability of QSEA for clinical studies. QSEA is implemented in R and available from the Bioconductor repository 3.4 (www.bioconductor.org/packages/qsea).

Figure 1.

Overview of the QSEA workflow. Green boxes represent data input, functions implemented in QSEA are depicted in blue, and red boxes describe the respective analysis step performed within these functions.

Figure 2.

Modelling MeDIP enrichment. (A) MeDIP read density (depicted exemplary for normal sample from patient 3) has linear relation to methylation level, and (B) gets saturated for high CpG density. (C) CpG dependent enrichment profile plot shows heat-color coded density of MeDIP reads by CpG density, the mean MeDIP coverage (red line), the observed MeDIP coverage at fully methylated regions (black line), and the sigmoidal function fitted to the coverage of fully methylated regions (dashed green line). (D) Exemplary illustration of methylation estimates depending on MeDIP read density, assuming four background reads and 25 reads at fully methylated windows.

Figure 3.

Benchmark and comparison to alternative methods. (A) Spearman correlation of QSEA, BayMeth and MeSiC estimates vs BS-seq for IMR-90 cell line. (B) High density scatterplots of MeDIP methylation estimates and 450k methylation levels for IMR 90 cell line. (C) Spearman correlation of QSEA, BayMeth and MeSiC estimates vs Methyl-Seq for chr1 of PDX samples. (D) Clustering of methylation differences between PDX and normal tissue for QSEA estimates.

Figure 4.

Validation and functional interpretation of differentially methylated Regions: (A) Scatterplot of methylation difference from QSEA methylation estimate and Methyl-Seq of CpG island promoter DMRs. Blue: directly covered by Methyl-Seq. Orange: neighbourhood covered by Methyl-Seq, red: not covered by Methyl-Seq. (B) Methylation difference of 20 most hypermethylated tumour suppressor genes form QSEA estimates (red) and Methyl-Seq(blue). (C) Histogram of correlation between CpG island promoter methylation and gene expression. (D) Mean methylation difference versus gene expression log2 ratio of differentially expressed genes with differentially methylated CpG island promoter. Red numbers are gene counts per quadrant. (E) Gene expression in four NSCLC cell lines after demethylation validation experiment.

Figure 5.

Enrichment of differentially methylated regions. (A) Sixteen most enriched transcription factor binding sites for hypermethylated regions and (B) hypomethylated regions respectively. (C) Tumour-Normal methylation differences chr22 show large hypomethylated blocks (LHB). Black dots represent mean methylation difference between PDX and normal, dashed lines are smoothed methylation differences for individual patients. Violet bars indicate two LHB, abundant in most of the tumour samples.