A novel statistical method for quantitative comparison of multiple ChIP-seq datasets

Abstract

Motivation: ChIP-seq is a powerful technology to measure the protein binding or histone modification strength in the whole genome scale. Although there are a number of methods available for single ChIP-seq data analysis (e.g. 'peak detection'), rigorous statistical method for quantitative comparison of multiple ChIP-seq datasets with the considerations of data from control experiment, signal to noise ratios, biological variations and multiple-factor experimental designs is under-developed.

Results: In this work, we develop a statistical method to perform quantitative comparison of multiple ChIP-seq datasets and detect genomic regions showing differential protein binding or histone modification. We first detect peaks from all datasets and then union them to form a single set of candidate regions. The read counts from IP experiment at the candidate regions are assumed to follow Poisson distribution. The underlying Poisson rates are modeled as an experiment-specific function of artifacts and biological signals. We then obtain the estimated biological signals and compare them through the hypothesis testing procedure in a linear model framework. Simulations and real data analyses demonstrate that the proposed method provides more accurate and robust results compared with existing ones.

Availability and implementation: An R software package ChIPComp is freely available at http://web1.sph.emory.edu/users/hwu30/software/ChIPComp.html.

Fig. 1.

Scatterplots of IP counts versus estimated background signals from the peak regions, in logarithm scale. The red dashed line is the result from cubic smoothing spline fitting (Color version of this figure is available at Bioinformatics online.)

Fig. 2.

Comparison of differential peak detection accuracies from simulations. The proportions of true discovery among top-ranked candidate regions is plotted against the number of top-ranked regions. (a) and (b) data are generated based on the proposed model. (c) and (d) data are generated based on the additive model. (a) and (c) are based on H3K27ac. (b) and (d) are based on PolII (Color version of this figure is available at Bioinformatics online.)

Fig. 3.

Histogram of P-values reported from different methods, based on null model that there’s no differential regions. The data are generated from the proposed model

Fig. 4.

Comparison of FDR estimations from different methods, based on simulation. X-axis shows the FDR reported from different methods, and y-axis shows the observed FDR (Color version of this figure is available at Bioinformatics online.)

Fig. 5.

Comparison of differential peak accuracies from real datasets. All results are for two-condition comparisons on different histone modification or protein binding, as marked in figure titles (Color version of this figure is available at Bioinformatics online.)