A linear time algorithm for detecting long genomic regions enriched with a specific combination of epigenetic states

Calculating a similarity vector

We need to generate a modification vector at each position from epigenomic signal data. For example, to create benchmark datasets in this study, we binarized the modification signal level at each position using BinarizeBed in ChromHMM [15], which classified the signal at each position into 0 or 1 according to a Poisson background model. Subsequently, we defined a modification vector as the vector with binary scores of modifications at each position.

Definition 1. Let w_l,w₂...,w _n be non-overlapping windows of the same length (e.g., 200 bp in this study) in a DNA sequence. Let $s_i^1,\dots ,s_i^k$ be binary or real-valued signals of k modifications in window i. The modification vector of w _i is defined as $M_i=\left(s_i^1,\dots ,s_i^k\right)$. Let F denote the feature vector of a focal modification pattern with k elements. The similarity score of M _i and F is defined as their inner product minus a given threshold τ.

Example. Suppose that k = 3, τ = 1.3, F = (1,1,0), M₁ = (1,1,0), M₂ = (1,0,1) and M₃ = (0,0,1). Similarity scores of F and M _i are 0.7, -0.3, and -1.3 for i = 1,2,3

When the inner product of M _i and F is positive for all i = 1,...,n, the optimal set of regions that maximizes the sum of similarity scores in the regions becomes the entire region, [1,n], which may not be informative. If we want to select a set of regions whose modification vectors are closer to the feature vector F, we can set the threshold τ to an appropriate positive value to yield a negative similarity score for the inner product that is lower than τ. Positions with negative similarity scores are less likely to be included in the optimal set of regions. A higher threshold is likely to divide the entire genome into smaller regions with a higher precision, while a lower threshold yields an opposite trend. In this manner, for a series of windows w_l,w₂...,w _n in a DNA sequence, we generate a series of similarity scores.

Detecting an optimal set of disjoint regions

To detect regions of focal epigenetic states such as K27HMD, we present an algorithm for calculating an optimal set of disjoint regions in a sequence that maximizes the sum of similarity scores for all regions. In addition, to identify sufficiently long regions, we define a minimum length threshold of regions such that each region is longer than or equal to the minimum length. The problem can be defined as follows.

Readers may find the last condition strange because it appears to disallow the situation that the first region starts at position a₁ >m₀. We used the condition to simplify the presentation of our linear-time algorithm, which is described later. To obtain such an optimal series of regions that the first region starts at a₁ >m₀, for example, you can temporarily add m₀ negative weights in front of L, calculate the optimal series, and restore the coordinate.

To calculate a C that maximizes the sum of weights in C, $\sum _{i\in I\in C}{L}_i$, we used a dynamic programming algorithm developed by Csuros [18]. Here, we outline the algorithm.

Definition 3. We assume that all series meet the conditions given in Definition 2. Let w(C) denote the sum of weights in C, $\sum _{i\in I\in C}{L}_i$. We consider two cases: that in which the last region of C ends at I and that in which it does not. When the last region does not end at i, let $C_{i,m}^0$ denote a series of regions that maximizes w(C) among all series, such that the last region ends at position b _k ≤ i - m, where m ≥ 1. When the last region ends at i, let $C_{i,m}^1$ denote a series of regions that maximizes w(C) among all series, such that the last region is of length ≥ m (≥1); specifically, a _k + m - 1 ≤ i (= b _k ).

According to this definition, C maximizing w(C) is either $C_{n,1}^0$ or $C_{n,m_1}^1$. For calculating these two series, we define here $w\left(C_{i,m}^0\right)$ and $w\left(C_{i,m}^1\right)$ recursively for i = 1,...,n and m ≥ 1.

Recall that C maximizing w(C) is either $C_{n,1}^0$ or $C_{n,m_1}^1$. From $W_{long}^1\left(n\right)$, we can build the series of regions, $C_{n,m_1}^1$, by tracing back the process of calculating $W_{long}^1\left(n\right)$. Similarly, from $W_{short}^0\left(n\right)$, we can obtain $C_{n,1}^0$.

We implemented the above idea. We call the program CSMinfinder.

Data sets

To compare CSMinfinder with other available methods for detecting large K27HMD, we used real biological datasets from the medaka-fish and human genomes, each of which was a set of vectors of DNA methylation levels at CpG sites, determined by bisulfite sequencing, and H3K4me3 and H3K27me3 histone modification Chip-seq data [8]. We set the window size to 200 bp, normalized the data using a Poisson distribution model, and set the binarized score of a window to 1 if its probability was < 0.0001 and to 0 otherwise.

Detecting large K27HMD in medaka epigenomic data

We compared CSMinfinder with ChromHMM [15] and Nakamura's method [8].

• Using ChromHMM, we estimated six chromatin states and divided the given DNA sequence into these six states. Specifically, ChromHMM asks users to input the number of epigenetic states beforehand. Thus, we started with inputting a small number into ChromHMM, increased the number gradually one by one until we found a state similar to K27HMD, hypo-methylated DNA modification and H3K27me3 histone modification, and called the number sufficient. Inputting a value larger than the sufficient number into ChromHMM did not make much sense because it just output a state similar to K27HMD. The sufficient number was six. Among the six states, one represented hypomethylated DNA modifications and the H3K27me3 histone modification. We therefore treated the state as K27HMD.

• Nakamura's method detects a hypomethylated domain on a DNA sequence that has more than nine contiguous CpG sites with low methylation (methylation level <0.4) and no more than four contiguous highly methylated CpG sites. Parameters are selected heuristically. A hypomethylated domain is treated as a K27HMD if it contains H3K27me3 peaks detected by QuEST [19], such that each peak is more than three times larger than the average.

• In CSMinfinder, we used two types of minimum length thresholds, 4 kbp and 8 kbp, to evaluate the effect of this constraint. We set the minimum length of any interval between regions to 600 bp.

Comparing the performance in detection of large K27HMD around genes in the medaka genome

Large K27HMD regions of length >4 kbp suppress the expression of many developmental genes [8]. Thus, we verified the effectiveness of CSMinfinder for detecting large K27HMD regions surrounding genes in the medaka genome. Nakamura's method could detect 246 large K27HMD regions containing the promoter regions of developmental genes (e.g., hox clusters) that were relevant to transcriptional regulation and the developmental process. CSMinfinder detected 911 K27HMD regions, and of these, 386 regions contained promoter regions of >4 kbp in size and contained 242 of the 246 regions identified using Nakamura's method. Indeed, our regions covered 91% of bases in the entire regions detected by Nakamura's method. Specifically, although the exact boundaries of individual regions estimated by the two methods were unlikely to be consistent, these regions largely overlapped each other. These results demonstrate the high concordance between CSMinfinder and Nakamura's methods as well as the ability of CSMinfinder to identify more K27HMD regions than did Nakamura's method.

We assessed the quality of each K27HMD region in terms of their low average DNA methylation level because this property is considered to be essential in maintaining the suppression of developmental gene expression in embryonic cells [8]. Indeed, Figure 2 shows the tendency of the average methylation level in the vertical axis to become lower for a longer K27HDM region, the length of which is displayed in the horizontal axis. This trend was also observed with all three methods.

We then compared the performance of the three methods by examining the length distributions of K27HMD regions in the medaka genome. Figure 3A shows the length distributions of large K27HMD regions (>4 kb in size) estimated by each of the three methods. Setting the minimum length threshold to 4 kbp in CSMinfinder detected more regions of length > 6 kbp but fewer regions of length >7 kbp compared with Nakamura's method. CSMinfinder allows us to output longer regions by setting the minimum length threshold to a higher value. For example, setting the minimum length to 8 kbp, CSMinfinder found more regions than did Nakamura's method (Figure 3C).

Analysis of large K27HMD regions in human epigenomic data

We also compared CSMinfinder with the other two for processing human epigenomic data. For ChromHMM, we calculated the sufficient number for the human data according to the procedure described before, and we classified epigenetic modification data into seven states rather than six so as to identify a state similar to K27HMD. The sufficient numbers of epigenetic states in the human and medaka data differed due to the difference in data quality. The sufficient number in the medaka data was smaller than that in the human data presumably because epigenetic state signals in the medaka data were clearer.

In CSMinfinder, we set the minimum length threshold to 8 kbp and the interval between regions to 600 bp. We also searched an ideal value of threshold τ by repeated trials to detect large continual regions, and we set τ to 1.3 and 1.6 in the respective medaka and human data.

Because the human genome is longer than the medaka genome, we focused on large K27HMD regions of length > 8 kbp. Nakamura's method detected 314 regions, and CSMinfinder identified 542 regions, including 291 of those found using Nakamura's method. Again, there was high concordance between the results obtained by the two methods. Figure 4 shows examples of large K27HMD regions detected around developmental genes. Although CSMinfinder and Nakamura's method yielded slightly different regions with distinct boundaries in the output, each created regions of similar sizes. In contrast, ChromHMM yielded shorter regions than the other two did. Specifically we compared the length distribution of large K27HMDregions estimated by each of the three methods (Figure 5). We found that CSMinfinder and Nakamura's method were comparable. Precisely, although the number of extremely large regions longer than 12 kbp is slightly smaller than the number found by Nakamura's method, CSMinfinder could detect similar numbers of large regions between 8 kbp to 12 kbp. Later we will discuss the reason why ChromHMM were inferior to the other two methods.

Computational performance and software availability

We observed the computational performance of CSMinfinder using Intel Xeon CPU E5-2670 processor with a clock rate of 2.60 GHz and 66GB of main memory. The computation time to calculate the optimal series of regions was negligible. Figure 6 shows that the average elapsed time was less than 2 seconds when we processed the epigenetic data of any of human and medaka chromosomes. Furthermore, Figure 6 also illustrates that the elapsed time is almost proportional to the size of each chromosome, thereby confirming experimentally that the worst-case time complexity of the algorithm is linear in the input size. CSMinfinder does not consume a large amount of main memory. CSMinfinder is made available at the following site: URL: http://mlab.cb.k.u-tokyo.ac.jp/~ichikawa/Segmentation/