A unified hypothesis-free feature extraction framework for diverse epigenomic data

Abstract

Motivation: Epigenetic assays using next-generation sequencing have furthered our understanding of the functional genomic regions and the mechanisms of gene regulation. However, a single assay produces billions of data points, with limited information about the biological process due to numerous sources of technical and biological noise. To draw biological conclusions, numerous specialized algorithms have been proposed to summarize the data into higher-order patterns, such as peak calling and the discovery of differentially methylated regions. The key principle underlying these approaches is the search for locally consistent patterns.

Results: We propose $L_0$ segmentation as a universal framework for extracting locally coherent signals for diverse epigenetic sources. $L_0$ serves to compress the input signal by approximating it as a piecewise constant. We implement a highly scalable $L_0$ segmentation with additional loss functions designed for sequencing epigenetic data types including Poisson loss for single tracks and binomial loss for methylation/coverage data. We show that the $L_0$ segmentation approach retains the salient features of the data yet can identify subtle features, such as transcription end sites, missed by other analytic approaches.

Availability and implementation: Our approach is implemented as an R package "l01segmentation" with a C++ backend. Available at https://github.com/boooooogey/l01segmentation.

Figure 1.

(A) Binomial, Gaussian, and Poisson $L_0$ segmentation execution times are plotted in log-space (both x and y) with respect to various segment lengths. The data are simulated using the corresponding probability distributions. Even though the theoretical worst-case complexity of the algorithm is $O(n^2)$, in practice, the execution time scales linearly with the length of the input. (B and C) Comparison of the running time of our proposed method against similar approaches. CROCS is used for epigenetic data. PELT is the algorithm used in the Yokoyama et al. approach for segmenting methylation data. (D) To estimate how long it would take for the methods to segment the entire Hg38 we benchmarked the methods on a 300KB segment and multiplied the results with 10 000. This likely underestimates both CROCS’ and PELT’s actual execution time.

Figure 7.

Comparing $L_0$ to $L_1$ with binomial error on real WGBS data with different fold compression computed on the chromosome level. The reported differentially methylated region is highlighted. (A) POU5F1 expressing cell line, H1. (B) The non-expressing contrasting cell-line, IMR90. Note that unlike Fig. 6, fold compression is fixed at the chromosome level, and consequently the number of segments in this region is not constant across $L_0$ and $L_1$.

Figure 8.

Comparison of $L_0$ and $L_1$ segmentations with binomial error in preserving the characteristics of genomic regions, namely “TssA” (TSS) and “Tx” (Transcript) regions identified by ChromHMM, in two cell lines: H1 and IMR90.

Figure 2.

Poisson $L_0$ segmentation is applied to ChIP-seq assays of CTCF, RNA Polymerase II, and three H3 histone modifications. The hyperparameter ($\lambda$) is chosen automatically using cross-validation.

Figure 3.

Comparison of $L_0$ Poisson (circle), L1 Poisson (triangle), and fixed size binning (square) for the various compression ratios on the DNase track from the H1 and IMR90 cell lines. The methods are applied on a 10M base-pair segment of Chromosome 1. We evaluated the methods using 937 peaks for H1 and 1400 peaks for IMR90 discovered by MACS within this segment. (A and B) The ratio of the mean signal within peak regions to the mean background signal after segmentation. (C and D) Maximum Jaccard Index distributions for the methods at different compression ratios. Jaccard indices are calculated between peaks and the segments discovered by the methods. (E and F) Median of the distributions shown in C and D.

Figure 4.

The ratio of mean signal within ChromHMM clusters to the mean background signal after compressing the signal 10K fold using $L_0$ Poisson (L0), L1 Poisson (L1), and fixed size binning (binned). The results are shown for DNase, H3K27ac, and H3K27me3 ChIP-seq tracks from H1 and IMR90 cell lines. (A) For all tracks and cell lines, $L_0$ Poisson retains the information content of the signal, while L1 Poisson and binning either completely lose the structure or diminish the signal intensity. (B) A closer inspection highlights the efficiency of $L_0$ Poisson in preserving the integrity of epigenetic signals after compression.

Figure 5.

(A) Comparison of the performance of different data reduction approaches on the transcript end site (TES) discovery task. We plot the median absolute distance to TES based on known gene models. The degree of segmentation is adjusted using lambda, while for peaks, we adjust the FDR cutoff to get a comparable number of segments. (B) Examples of segmentation tracks for selected transcripts. Highlighted regions on the gene annotations indicate MACS peaks. While the Pol II signal is often noisy outside of the promoter region, a subtle drop-off in the signal can often be observed toward the end of transcripts. In many cases, this change in signal is correctly identified by L0 Poisson segmentation.

Figure 6.

Binomial $L_0$ segmentation accurately accounts for read coverage and discovers short regions with distinct methylation rates. We compare segmentation formulations on simulated data. Coverage and methylation are shown in the first and second rows respectively. (A) Gaussian (bottom row) and binomial error segmentation (the third row) with the same number of segments. Binomial error merges low-coverage regions with their neighbors. (B) Comparing $L_0$ (the third row) and $L_1$ (the bottom row). The $L_0$ penalty is more sensitive to local structure discovering the 2 CpGs simulated with different $\beta$.