Saturation analysis of ChIP-seq data for reproducible identification of binding peaks

Abstract

Chromatin immunoprecipitation coupled with next-generation sequencing (ChIP-seq) is a powerful technology to identify the genome-wide locations of transcription factors and other DNA binding proteins. Computational ChIP-seq peak calling infers the location of protein-DNA interactions based on various measures of enrichment of sequence reads. In this work, we introduce an algorithm, Q, that uses an assessment of the quadratic enrichment of reads to center candidate peaks followed by statistical analysis of saturation of candidate peaks by 5 ′ ends of reads. We show that our method not only is substantially faster than several competing methods but also demonstrates statistically significant advantages with respect to reproducibility of results and in its ability to identify peaks with reproducible binding site motifs. We show that Q has superior performance in the delineation of double RNAPII and H3K4me3 peaks surrounding transcription start sites related to a better ability to resolve individual peaks. The method is implemented in C++ and is freely available under an open source license. [Supplemental material is available for this article.] Chromatin immunoprecipitation (ChIP) followed by massively parallel sequencing (ChIP-seq) is designed to detect genome-wide protein-DNA interaction. ChIP-seq can identify both sharp peaks typically associated with sequence-specific transcription factors , as well as broad histone-modification signals (Park 2009; Peng and Zhao 2011), and has become a central technology for the investigation of gene regulation. The ChIP-seq procedure involves formaldehyde-mediated crosslinking of chromatin followed by fragmentation of protein-DNA complexes into short fragments, which are then subjected to immunoprecipitation using an anti-body directed against a protein of interest (e.g., a transcription factor or a modified histone), thereby enriching genomic segments that are bound by the protein of interest prior to sequencing (Laajala et al. 2009). A crucial challenge in the computational analysis of ChIP-seq data pertains to finding peaks in ChIP-seq data that correspond to protein-DNA binding sites. Numerous peak calling algorithms have been presented, most of which address the same basic analytical tasks with methods to estimate the mean DNA fragment length from the data, to shift or extend the reads toward the center of the binding peak, to identify candidate peak regions, and to evaluate the statistical significance of the read depth of the candidate peaks. The sequence reads represent only the 5 ′ ends of the co-precipitated DNA fragments, which are generally 100-to 500-bp in length. Around true binding sites of the target protein, this results in a characteristic bimodal distribution of reads on the forward and reverse strands, which depends on the distribution of fragment lengths in the library and can be exploited for signal detection and evaluation. Therefore, an initial step in many algorithms is the estimation of the actual fragment-length distribution. Following fragment-lengt...

Characterization of hundreds of regulatory landscapes in developing limbs reveals two regimes of chromatin folding

“…For figure display purposes, replicate ChIP-seq tracks were merged. RAD21 ChIP-seq tracks were produced using the coverage function of the Q-PeakCaller (Hansen et al 2015). Peaks of CTFC and Rad21 were called on merged replicates with the Q-PeakCaller.…”

Q-nexus: a comprehensive and efficient analysis pipeline designed for ChIP-nexus

“…3 a). Here, we present a new method for the estimation of the protected-region width in ChIP-nexus data that uses the distribution of qfrag [19] lengths. A qfrag is defined to be the genomic interval between any pair of 5’ read mapping positions on the forward and reverse strand.…”

“…Since PCR artifacts (IMIB) are already removed, IMUB reads are kept, assuming that such reads stem from different molecules because they have different random barcodes. Our algorithm (see Methods) implements the method of qfrag-length distribution with pseudo-control in order to estimate the protected-region width ℓ ′′′ , which is then used to combine pairs of 5’ ends on the forward and reverse strand by forming qfrags [19] with a minimal allowed distance q min = ℓ ′′′ −5 and a maximal allowed distance of q max = ℓ ′′′ +5. The qfrag-depth at any one position is the total number of qfrags that cover the position.…”

“…For ChIP-nexus, it is not the size of the original fragments that is important, but rather the segment of DNA that cannot be digested because of steric interference by the formaldehyde cross-linked protein (the “protected region”). We present a method for estimating the average width of the protected region that is based on the notion of qfrags [19] and show that, on the ChIP-nexus data, it yields unbiased signatures that are not affected by the so-called phantom peak [20], which is not the case for the cross-correlation method developed for ChIP-seq [2]. The estimates of the protected region width are in line with distances that were derived in a previous study of same data [11] using integrated 5’ end coverage plots around predicted and motif-centered sites.…”

“…We have previously developed a method using “qfrag-analysis” to identify candidate peaks in ChIP-seq analysis [19]. Here, we adapted that algorithm for ChIP-nexus analysis.…”