iSeg: an efficient algorithm for segmentation of genomic and epigenomic data

High throughput genomic assays, such as microarrays and next-generation sequencing, are powerful tools for studying genetic and epigenetic functional elements at a genome scale [1]. A large number of approaches have been developed to exploit these technologies to identify and characterize the distribution of genomic and epigenomic features, such as nucleosome occupancy, chromatin accessibility, histone modifications, transcription-factor binding, replication timing, and DNA copy-number variations (CNVs). These approaches are often applied to multiple samples to identify differences in such features among different biological contexts. When detecting changes for such features, one needs to consider a very large number of segments that may undergo changes, and robustly calculating statistics for all possible segments is usually not feasible. As a result, heuristic algorithms are often needed to find the optimal solution for the objective function adopted by an approach. This problem is often called segmentation problem in the field of genomics, and change-point problem in other scientific disciplines.

Solving the segmentation problem typically involves dividing a sequence of measurements along the genome such that adjacent segments are different for a predefined criterion. For example, if segments without changes have a mean value of zero, then the goal could be to identify those segments of the genome whose means are significantly above or below zero. A large number of such methods have been developed for different types of genomic and epigenomic data [2–17]. Many methods are designed for specific data types or structures, but it is challenging to find versatile programs for data with different properties and different underlying statistical assumptions. The previous methods fall into several categories including change-point detection [2, 3, 9, 10, 12, 14, 18–25], Hidden Markov models [5, 15, 26–28], Dynamic Bayesian Network (DBN) models [29, 30], signal smoothing [31–34], and variational models [35, 36]. For review and comprehensive comparison, please refer to [16, 37–40].

Many currently available segmentation tools have poor performance, run slowly on large datasets, or are not straightforward to use. To address these challenges, we developed a general method for segmentation of sequence-indexed genome-wide data. Our method, called iSeg, is based on a simple formulation of the optimization problem [7]. Assuming the significance (i.e. p-value) of segments can be computed based on certain parametric or non-parametric models, iSeg identifies the most significant segments, those with smallest p-values. Once the segment with the smallest p-value is found, it will be removed from the dataset and segment with the second smallest p-value will be searched in the remaining of the data. The procedure repeats until no segments whose significance levels pass a predefined threshold. This simple objective function is intuitive from a biological perspective since the most statistically significant segments often biologically significant.

iSeg has several noteworthy features. First, the simple formulation allows it to serve as a general framework to be combined with different assumptions of underlying probability distributions of the data, such as Gaussian, Poisson, negative Binomial, or non-parametric models. As long as p-values can be calculated for the segments, the corresponding statistical model can be incorporated into the framework. Second, iSeg is a general segmentation method, able to deal with both positive and negative signals. Many of the existing methods, cannot deal with negative values in the data, because they are designed specifically for certain data types, such as genome-wide read densities, assuming data values with only zero or non-negative values. Negative values in genomic datasets can occur when analyzing data pair relationships, such as difference values or log2 ratios commonly used with fold-change analysis. Currently, most methods segment profiles separately and compare the resulting segmentations. The drawback of such treatment is that peaks with different starting and ending positions from different segmentations cannot be conveniently compared, and peaks with different magnitudes may not always be distinguished. Taking differences from two profiles to generate a single profile overcomes these drawbacks. Third, to deal with cases where segments are statistically significant, but the biological significance may be weak, we apply biological significance threshold to allow practitioners the flexibility to incorporate their domain knowledge when calling “significant” segments. Fourth, iSeg is implemented in C++ with careful design of data structures to minimize computational time to accommodate, for example, multiple or very long whole genome profiles. Fifth, iSeg is relatively easy to use with few parameters to be tuned by the users.

Here we describe the method in detail, followed by performance analysis using multiple data types including simulated, benchmark, and our own. The data types include DNA copy number variations (CNVs) for microarray-based comparative genomic hybridization (aCGH) data, copy number variations from next generation sequencing (NGS) data, and nucleosome occupancy data from NGS. From these tests, we found that iSeg performs at least comparably with popular, contemporary methods.

We compared our method with several previous methods for which we were able to obtain executable programs: HMMSeg [5], CGHSeg [10], DNAcopy [9, 24], fastseg [3], cghFLasso [34], BioHMM-snapCGH [27], mBPCR [12], SICER [47], PePr [48] and MACS [49]. Among them, CGHSeg, DNAcopy, BioHMM-snapCGH, mBPCR and cghFLasso are specifically designed for DNA copy number variation data; MACS, SICER and PePr are designed for ChIP-seq data; and HMMSeg is a general method for segmentation of genomic data. Each method has some parameters that can be tuned by a user to achieve better performance. In our comparative study, we carefully selected parameters on the basis of the recommendations provided by the authors of the methods. For each method including iSeg, a single set of parameters is used for all data sets except where specified. Post-processing is required by some of the methods to identify significant segments.

The methods CGHSeg, DNAcopy, and fastseg depend on random seeds given by a user (or at run-time automatically), and the F₁-score at different runs are very similar but not the same. These methods were run using three different random seeds. The averages of the F₁-score were used to measure their performance.

Performance on simulated data

The simulated profiles were generated under varying noise conditions, with signal to noise ratios (SNR) of 0.5, 1.0 and 2.0, which correspond to poor, realistic and best case scenarios, respectively. Ten different profiles of length 5000 were simulated.

For each profile, five different segments of varying lengths were predefined at different locations. Data points outside of these segments were generated from normal distribution with mean zero. The five segments were simulated with non-zero means and varying amplitudes, or ease of detection, in order to assess the robustness of the methods. Because this set of simulated data resembles more of the DNA copy-number variation data, we used it to compare iSeg to methods designed for DNA copy number data. Figure 2 shows an example of the simulated data and the segments identified by iSeg and other existing methods. Figure 3a shows the performance of iSeg and other methods on simulated data with SNR = 1.0. We can see that iSeg, DNACopy and CGHSeg perform similarly well, with HMMseg and CGHFLasso performing a little worse while fastseg did not perform as well as the other methods. iSeg is also tested using a set of 10 longer simulated profiles, each with length 100,000. Seven segments are introduced at varying locations along the profiles. iSeg still performs quite well in these very long profiles. The performance of these methods on long sequences is shown in Fig. 3b. In Fig. 4, we plotted Venn Diagrams of the overlapping called locations for several methods with better performance in terms of F₁-scores, for both the simulated short profiles (Fig. 4a) and long profiles (Fig. 4b). We can see that in the simulated profiles, there is a substantial overlap among most methods, while iSeg can detect more segments for relatively short profiles. In simulated long profiles, all methods, except HMMSeg, have similar segmentation results.

Performance on experimental data

DNA copy-number variation (CNV) data

To assess the performance of iSeg on experimental data, we use three different datasets. They were the Coriell dataset [50] with 11 profiles, the BACarray dataset [51] with three profiles, and the dataset from The Cancer Genome Atlas (TCGA) with two profiles. The 11 profiles in Coriell datasets correspond to 11 cell lines: GM03563, GM05296, GM01750, GM03134, GM13330, GM01535, GM07081, GM13031, GM01524, S0034 and S1514. We constructed “gold standard” annotations using a consensus approach. We first ran all the methods using several different parameter settings for each method. The resulting segments from all the different parameter settings of all the methods were combined to give an initial set of potential segments. The test statistics and p-values are then calculated for all the segments using the same probability distribution assumption described in Method. Benjamini-Hochberg procedure was then used to correct for multiple comparison. The 0.05 adjusted p-value cutoff was then used to select the set of segments as the gold standard. The annotations derived using the consensus approach are provided as Additional file 1.

The 11 profiles from the Coriell dataset were segmented using iSeg and the other methods. Segmentation result for one of the profiles is shown in Fig. 5, and the F₁-scores are shown in Fig. 6a. The performance of iSeg is robust with accuracy above 0.75 for all the profiles from this dataset and it was found to be comparable to or better than other methods. For HMMSeg, both no-smoothing and smoothing were used. The best smoothing scale for HMMSeg was found to be 2 for the Coriell dataset. In Fig. 5, we found that iSeg identified most of the segments. DNAcopy, fastseg, HMMSeg and cghseg missed single-point peaks, whereas cghFLasso, mBPCR and snapCGH missed some of the longer segments. The segmentation results for other profiles in Coriell dataset can also be found in the Additional file 1. We generated annotations using the consensus method for BACarray dataset similar to the Coriell dataset. The comparison of segmentation results for one profile of the BACarray dataset is shown in Fig. 7, and the comparison of F₁-scores is shown in Fig. 6b. iSeg returned better F1-scores than the other methods, consistent with the conclusions based on visual inspection.

For the TCGA datasets, since the profiles are rather long, we did not generate annotations using the consensus approach. Instead, we applied some of the methods to this dataset and compared their segmentation results visually (Fig. 8). Again, we found that iSeg identified most of the significant peaks. In this test, DNACopy performed well overall, but tended to miss some of the single-point peaks, whereas other methods performed even less well.

We compared the computational time of iSeg, shown in Table 1, with those of the other methods, and found that iSeg is the fastest for the three test datasets. Notably, iSeg took much less time than the other methods for very long profiles (length 100,000). This speed is achieved in part through dynamic programming and a power factor that provides rapid initial scanning of the profiles. The long profiles contain similar amount of data points that are signals (as opposed to background or noise) as the shorter profiles. The time spent on dealing with potentially significant segments is roughly the same between the two types of profiles. As a result, the overall running time of iSeg for the long profiles did not increase as much as that for the other methods. In summary, we observed that iSeg ran faster than the other methods, especially for profiles with sparse signals.

Method	Simulation (SNR ≃ 1.0, n = 5000)	Simulation (SNR ≃ 1.0, n = 100 K)	Coriell
iSeg (C++)	0.164	1.223	0.294
DNAcopy (R)	2.267	60.343	3.098
Fastseg (R)	0.647	48.139	0.630
CGHSeg (R)	54.480	157.626	24.36
HMMSeg (Java)	0.543	160.790	0.552

Differential nuclease sensitivity profiling (DNS-seq) data

We then tested our method on next generation sequencing data, for which discrete probability distribution models have been used in most of the previous methods. The dataset profiles were genome-wide reads from light or heavy digests (with zero or positive values) or difference profiles (light minus heavy, with positive or negative values) [51]. The difference plots are also referred to as sensitivity or differential nuclease sensitivity (DNS) profiles [51].

Segmentation of single nuclease sensitivity profiles

Figure 9 shows the significant segments (peaks) called by iSeg together with those of two other methods, MACS and SICER. Visual inspection revealed that iSeg successfully segmented clear peaks and their boundaries. The performance of iSeg was at least comparable to MACS and SICER. Follow-up analyses on the segmentation results also demonstrated its capability in identifying biological interesting functional regions [51, 52]. iSeg results with different biological significance cutoffs (BC) are displayed as genome browser tracks beside the input profile data to guide inspection of the segmentation results. The default is to output three BCs: 1.0, 2.0, and 3.0, which are sufficient for most applications.

Segmentation of difference profiles with both positive and negative values

Difference profiles between two conditions can be generated by subtracting one profile from the other at each genomic location. Pairwise comparisons are of great interest in genomics, as they allow for tests of differences within replicates, or across treatments, tissues, or genotypes. Analyzing such profiles can preserve the range or magnitude of differences, adding power to detect subtle differences between two profiles, compared to approaches that rely on calling peaks in the two files separately. Figure 10 shows the segmentation of iSeg on a typical DNS-seq profile. When analyzing these sets of DNS data, we merge segments only when they are consecutive, meaning the gap between the two segments is zero. The length of gaps between adjacent segments that can be merged is a parameter tunable by users. iSeg successfully identified both positive (peak, positive peak) and negative (valley, negative peak) segments, as shown in Fig. 11. Most existing ChIP-seq data analysis methods do not accommodate this type of data as input. To run MACS, SICER, and PePr, we assigned the light and heavy digestion read profiles as the treatment and control files respectively, as a fair way of comparison. Since the true biological significant segments are unknown, we compared the methods through careful visual inspections by domain experts. We found that iSeg performed satisfactorily and select it as the method of choice for analyzing the data from our own labs. The choice of parameters of competing methods is somehow subjective. We always start from the default values, and then do some tuning of each parameter while making the others fixed. Each time we will do a careful visualization until the final set of calls look reasonably satisfactory.

In this study, we designed an efficient method, iSeg, for the segmentation of large-scale genomic and epigenomic profiles. When compared with existing methods using both simulated and experimental data, iSeg showed comparable or improved accuracy and speed. iSeg performed equally well when tested on very long profiles, making it suitable for deployment in real-time, including online webservers able to handle large-scale genomic datasets.

In this study, we have assumed that the data follow a Gaussian (normal) distribution. The algorithm is not limited, however, to this distribution assumption. Other hypothesis tests, such as Poisson, negative binomial, and non-parametric tests, can be used to compute p-values for the segments. Data generated by next-generation sequencing (NGS) are often assumed to follow either a Poisson or a negative binomial distribution [53, 54]. However, a recent study by us (to be published) showed that normality-based tests such as t-tests or Welch’s t-test can perform equally well for NGS data after certain normalization, indicating Normal approximation of NGS data is likely adequate when analyzing NGS data, at least for some applications. For segmentation problems, when the segments tend to have relatively large sizes, the averages of signals within segments can be well approximated by normal distributions. This likely explains why iSeg performed well on NGS data. However, p-values of the segments can be computed based on either Poisson or negative binomial distributions, which can be directly used in the rest of the segmentation procedures. Such flexibility and general applicability make iSeg especially useful for segmenting a variety of genomic and epigenomic data.

The use of dynamic programming and a power factor makes iSeg computationally more efficient when analyzing very large data sets, especially with sparse signals as is typical for many types of NGS data sets. The refinement step identifies the exact boundaries of segments found by the scanning step. Merging allows iSeg to detect segments of any length. Together, these steps make iSeg an accurate and efficient method for segmentation of sequential data.

The statistic used in our method is very similar to that described for optimal sparse segment identification [7]. However in that study, the segments are identified using an exhaustive approach, which will not be efficient for segmenting very large profiles. To speed up computation, the optimal sparse segment method [7] employs the assumption that the segments have relatively short length, which is not true for some datasets. In contrast, the algorithm designed in this study allowed us to detect segments of any length with demonstrably greater efficiency.

The gold standard approach is desirable but the derivation of a true gold standard is complex and possibly subjective. A gold standard generated using a consensus approach does not guarantee that the true optimal segments will be identified. In addition, the F₁-scores may favor iSeg more as the test statistic used to generate the gold standard is not employed by the other methods. However, the statistic we used is based on model assumptions used by many existing methods. Visual inspection of segmentation results clearly shows how iSeg performs well in direct comparisons with other popular, contemporary methods. We expect that future research will benefit from the application of iSeg to compare multiple profiles simultaneously.

We have designed the method to make it flexible and versatile. This resulted in a number of parameters that users can tune. However, the default values work well for all the simulated and experimental datasets. In practice, to obtain satisfactory results, users are not expected to modify any parameters. The speed of iSeg would allow us and fellow researchers to implement it as an online tool to deliver segmentation results in real-time.

In this work, we developed a new method, iSeg, for the segmentation of large-scale genomic and epigenomic profiles. The performance of iSeg was demonstrated by comparing with existing segmentation methods using both simulated and real datasets. The computation framework utilized in iSeg can readily accommodate different assumptions on the underlying probability distribution of the data. The computational algorithms designed in iSeg make it able to search for, grow, and refine segments comprehensively for very long profiles without compromising too much of the computational speed. iSeg also provides a biological cutoff parameter, allowing researchers to select segments which are more likely to be biological significant by setting a proper biological cutoff value based on the problems at hand and/or their prior knowledge. We believe that iSeg will be a very useful tool for segmentation problems in biology.