Recent advances in ChIP-seq analysis: from quality management to whole-genome annotation.

Chromatin immunoprecipitation followed by sequencing (ChIP-seq) analysis can detect protein/DNA-binding and histone-modification sites across an entire genome. Recent advances in sequencing technologies and analyses enable us to compare hundreds of samples simultaneously; such large-scale analysis has potential to reveal the high-dimensional interrelationship level for regulatory elements and annotate novel functional genomic regions de novo. Because many experimental considerations are relevant to the choice of a method in a ChIP-seq analysis, the overall design and quality management of the experiment are of critical importance. This review offers guiding principles of computation and sample preparation for ChIP-seq analyses, highlighting the validity and limitations of the state-of-the-art procedures at each step. We also discuss the latest challenges of single-cell analysis that will encourage a new era in this field.

Introduction

Genome-wide investigations into cooperative interactions among genomic functions, e.g. DNA replication, segregation, translation, repair and rearrangement, are vital for systematically elucidating all biological activities on the genome. To this end, chromatin immunoprecipitation followed by sequencing (ChIP-seq) analysis was developed to understand the cooperation and interactions that occur in a wide variety of organisms using next-generation sequencing (NGS) [1-3]. ChIP-seq analysis is a mainstream method in genomics and epigenomics, and has led to important discoveries related to disease-associated transcriptional regulation [4-7], tissue-specificity of epigenetic regulation [8, 9] and chromatin organization [10-13].

ChlP-seq protocols have many steps involving sample preparation and computational analysis (Figure 1). In brief, cross-linked chromatin is sonicated, and purified with and without immunoprecipitation (ChIP and corresponding input DNA fragments, respectively). DNA fragments are sequenced as reads, which are then mapped onto the reference genome, and the genomic regions that are significantly enriched for ChIP reads, compared with input reads, are detected as peaks. Other genomic regions are regarded as non-specific background. Called peaks, which represent candidates of targeted protein/DNA-binding and histone modification sites, can be used to identify associated functional annotations, including binding motifs [14, 15] and gene ontology [16, 17]. ChIP-seq results can also be integrated with other types of genomic assays, including gene expression, DNA methylation and chromatin conformation, to understand mechanisms of genomic functions from multiple aspects [18-20].

Figure 1.

ChIP-seq analysis workflow. Boxes indicate the steps involved in ChIP-seq analyses for various aims discussed in this review. The considerations for each step are itemized. (A) Sample preparation, sequencing and mapping. This procedure is common to both (B) and (C). (B) Small-scale analysis (single or a few samples). In this case, adjusting peak-calling strategy and parameters to each sample’s property is possible. (C) Large-scale analysis (many samples). Left rectangles indicate the different experiments (e.g. same analysis for different cell types). Because integrative analysis is sensitive to the quality of input samples and one-by-one adjusting is difficult, objective quality metrics for multilateral quantitative assessment is necessary to filter poor-quality data automatically.

The shapes of the peaks vary among proteins, and are classified into three modes [1]: ‘sharp mode’, located at specific positions in the genome; ‘broad mode’, associated with large genomic domains; and ‘mixed mode’, which involves both peak modes. As most point-source transcription factors (TFs) and localized chromatin markers (e.g. H3K4me3) have sharp modes, a large majority of peak-calling algorithms have been designed for this mode, even though other proteins (e.g. heterochromatin protein HP1 [21]) and some histone modifications (e.g. H3K9me3) have broad modes. The mixed mode is observed for RNA polymerase II (Pol II) and transcription elongation factors [22]. The read distribution of Pol II-related NGS analyses, such as Nascent RNA-seq [23] and NET-seq [24], are also classified as mixed mode. Different peak-calling strategies are required for each shape.

Recent advances in sequencing technologies and analyses enable us to handle hundreds of ChIP samples simultaneously; such large-scale analyses revealed the high-dimensional interrelationship for regulatory elements [25-27] and annotate novel functional genomic regions de novo [28, 29]. Because a large-scale analysis is sensitive to the quality of input samples and adjusting the protocols for each sample’s quality is difficult, samples which have insufficient quality should be rejected automatically. As there are various factors (including antibody quality) during sample preparation that affect the quality of the obtained results [30, 31], multilateral quality assessments during the computational procedures are essential. Despite great efforts to streamline this process, no single workflow that is optimal under all circumstances exists, and there are many experimental considerations that are relevant to the method choice for a ChIP-seq analysis. Consequently, to obtain high-quality, unbiased and reasonable data, the overall protocol design and quality management, which are adjusted to the studies’ properties, are of critical importance.

In this review, we describe the computational protocol and sample preparation for a ChIP-seq analysis, and discuss the validity and limitations of emerging programs and quality measures currently available for specific analytical tasks by providing concrete examples. There are ChIP-seq-extended methods that can detect DNA-protein binding sites with base-pair resolution [32, 33]. For simplicity, we have limited the scope of this review to ChIP-seq protocols, for which the computational protocols are similar. The key issues described here also underpin protocols for other NGS-based analyses [34-37]. Finally, we discuss the latest challenges of single-cell analysis that will encourage a new era in this field.

Sample preparation and sequencing

Fragmented DNAs (150–500 bp) from ChIP-seq samples are sequenced as reads (36–100 bp). Single-end reads are often used for typical ChIP-seq analyses, while paired-end ones improve the library complexity and increase mapping efficiency at repetitive regions [38]. When research focuses on repetitive regions, longer and/or paired-end reads are preferred. While paired-end reads can be used to obtain the fragment size distribution, several methods exist that estimate it from single-end mapped data [12, 39, 40].

The chromatin accessibility during fragmentation is not uniform across the genome. In some open-chromatin regions (e.g. actively transcribed promoter regions), DNA is amenable to fragmentation and thereby preferentially represented in the fragmented sample, which causes false-positive read enrichment [41]. Tightly packed regions, e.g. heterochromatin, are sheared to a lesser extent by DNA fragmentation, thereby confounding weak enrichment of true binding sites for heterochromatin markers [38]. These fragmentation biases in a genome-wide read distribution profile should be taken into account when using null model analysis to obtain meaningful conclusions. One way to mitigate this fragmentation bias is to shear longer DNA fragments (350∼800 bp) further using ultrasonication after the immunoprecipitation step [4]. Although including longer fragments widens the obtained peaks, peak-summit resolution is not strongly affected [38].

Read mapping

Sequenced reads are mapped onto the genome using mapping tools [42, 43]. Most ChIP-seq experiments do not require gapped alignments that consider insertions and deletions (indels) because the sequenced reads do not contain them, unlike exon junctions in RNA-seq analyses. The exception is cross-species analysis, which maps reads onto different species’ genomes. If the information about heterozygous variants (e.g. single nucleotide polymorphisms and indels) in the reference genome is available, the allele-specific regulation analysis (personalized genome analysis) can be applied [3, 44, 45].

An important issue concerns the inclusion of multiple mapped reads (reads mapped to multiple loci on the reference genome). Allowing for multiple mapped reads increases the number of usable reads and the sensitivity of peak detection; however, the number of false positives may also increase [46]. In general, uniquely mapped reads are sufficient to analyze typical TFs, except for in-repeat analyses [47]. Considering the percentage of mapped reads (mapping ratio) is important, and desirable rate depends on the species and the read lengths.

Mappability

Recent central ChIP-seq studies [29, 45] used uniquely mapped reads. Instead of including multiple mapped reads, they considered the mappability [48] to correct for the loss of true signals in low-mappable regions. Theoretically, the mappability of a reference genome depends on the read length, read type (color or nucleotide space) and the mapping tool and parameters used [49], but calculating the genome-wide mappability for each is often time-consuming. Moreover, it is difficult to calculate the mappability of paired-end and gapped alignment data, although there has been an effort to calculate the former [50]. Consequently, it is practical to use the mappability data publicly available for similar parameter sets. When low-mappable regions (e.g. a ratio  <  0.25) are of interest, it might be better to include multiple mapped reads or use paired-end reads.

Library complexity

Library complexity is measured by the non-redundant fraction (NRF), the fraction N, where N and N are the numbers of non-redundant reads and the total number of mapped reads, respectively. The non-redundant reads are defined as reads mapped on the same genomic positons T times or less, where T is the threshold for redundant reads. The redundant reads (N − N) should be filtered from further analysis. For human, T is typically set 1 because the expected number of mapped reads per base pair (sequencing depth) ≪ 1. When it becomes  >1, due to the small genome under investigation (e.g. yeast), or in enriched regions for the very high signal-to-noise ratio (S/N) of the antibody (e.g. Pol II), it may be appropriate to relax the threshold T because stringent filtering has small effect on the sensitivity of the peak detection [38]. Moreover, when studying highly repetitive regions (e.g. rDNA regions in Saccharomyces cerevisiae), filtering redundant reads should be omitted.

Because the NRF score depends on the total number of mapped reads, read sampling is necessary when comparing NRF scores among samples. The ENCODE consortium endorses an NRF  > 0.8 for 10 million reads (T  =  1) [30]. A low library complexity often occurs when samples are prepared from a small amount of starting materials. Even if the number of sequenced reads is sufficient after polymerase chain reaction (PCR) amplification, the substantial read number may be small, resulting in a poor significance power.

Sequencing depth

The number of called peaks increases with the sequencing depth, because weaker sites become statistically significant with a greater number of reads [51]. Although early ChIP-seq analyses produced  <10 million reads per sample, it was reported that some highly active regions displayed a modest ChIP enrichment [30] and that weak protein binding by other factors may have important subfunctions [52]. Therefore, deep sequencing is required to include all functional sites. Broad markers that cover large genomic areas require more depth than those of sharp-mode peaks.

A sufficient sequencing depth depends on the S/N of the antibody. To determine sufficient sequencing depths, a ‘saturation analysis’ can be used, which subsamples the original read set in a stepwise fashion and calculates the proportion of identified peaks that overlap the original ones for each depth [1]. The point where the proportion is saturated is defined as sufficient depth. Although this approach is useful, there is no clear saturation point for most histone modifications [53]. Therefore, an agreeable depth has been determined empirically. For human samples, the ENCODE consortium suggested that at least two biological replicates with 10 million uniquely mapped reads for each replicate (providing at least 20 million reads per factor) is a minimum for typical TFs (sharp mode) [30]. Chen et al. [38] suggested that up to 60 million reads may be required for broad histone markers, and Jung et al. [53] suggested 40–50 million reads as a practical minimum for most broad histone markers. If the saturation point has not been detected at the available depth, it is still possible to apply tools for a sufficient depth estimation using a power analysis [54] or for predicting the benefits of additional sequencing [55].

Signal-to-noise ratio (S/N)

The S/N is evaluated by the number and strength of peaks obtained for each ChIP sample. This measure can also be used to assess the degree of noises in the input sample. The ENCODE consortium proposed two metrics, fraction of reads in peaks (FRiP) and cross-correlation profiles (CCPs) to measure the S/Ns [30]. The FRiP value is calculated as FRiP  =  N/N, where N is the number of reads falling within peak regions. This value correlates positively with the number and intensity of the identified peaks. ChIP samples that have too few peaks can be filtered using a cutoff for this score. However, because the FRiP score obviously depends on the sequencing depth and the parameters set for peak calling, it is not a perfectly objective metric. Conversely, CCPs assess the read-clustering levels without calling peaks beforehand [30]. This analysis plots the Pearson cross-correlations (CCs) between mapped read densities of positive and negative strands (y-axis) with shifting one strand (x-axis). Samples with large and small S/Ns typically have high CCs at the shift points corresponding to the fragment length (C) and the read length (C), respectively. Based on this observation, two quantitative measures are scored, the normalized strand coefficient NSC  =  C and relative strand correlation RSC = (C − C)/(C − C), where C is the minimum CC observed. The ENCODE consortium recommends an NSC  ≥ 1.05 and an RSC  ≥ 0.8 for typical TFs (sharp mode). Conversely, input and negative control samples should have low scores. Using this criteria, Marinov et al. [56] reported that a substantial minority (20%) of vertebrate ChIP-seq data sets for TFs in the Gene Expression Omnibus (GEO) were of insufficient quality, suggesting the necessity of quality check even for published data. Hansen et al. [12] has proposed the Hamming distance plot analogous to CCPs. CCPs are helpful; however, they have mainly been tested using only a few species and a more extensive investigation is necessary to understand the applicability of this approach to many other species. Moreover, a large S/N does not guarantee that the identified peaks are genuine binding sites—a large score merely means that there are many read-enriched regions in the genome. Samples that have many false-positive sites (e.g. non-specific binding sites) also have large S/Ns.

Example of the visualization and quality check

In Figure 2 we show an example of human ChIP-seq data obtained from the ENCODE project [45]. This example includes sharp mode (H3K4me3), broad mode (H3K36me3, H3K27me3 and H3K9me3) and mixed mode (RNA Pol II) samples with the control. Figure 2A summarizes the statistics of the quantitative quality metrics. The differences in average peak widths reflect the binding modes of each sample. The larger mapping ratio of H3K9me3 compared with H3K36me3 is mainly a consequence of the greater read quality. Conversely, the low affinity of the H3K9me3 antibody is indicated by the small RSC score and a comparable FRiP score against a much larger total peak width. These metric scores are used to automatically assess sample quality, leading to the rejection of poor-quality samples from further analyses.

Figure 2.

Statistics and visualization of ChIP-seq analysis for human K562 cells. A representative data set of ENCODE consortium [45]. The sequenced read files (fastq) and the reference peak lists (detected by Scripture [57] under the assumption of uniform background signal) were downloaded from GEO under accession number GSE29611. The fastq files were mapped onto the human genome (UCSC hg19) using Bowtie version 1.1.0 [42], allowing uniquely mapped reads only. (A) Summary statistics for each sample. The averaged read quality was obtained using fastqc version 0.11.4 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). The number of non-redundant reads, library complexity for 10 million mapped reads and FRiP scores were calculated using DROMPA3 version 3.0.0 [58]. Normalized strand coefficient (NSC) and relative strand correlation (RSC) scores were obtained using phantompeakqualtools version 1.1 (https://code.google.com/p/phantompeakqualtools). (B) The non-redundant read distribution for each sample with a RefSeq gene annotation (chromosome 1, 244.5–245.1 Mb). For the gene line (yellow box), genes in the upper and lower halves are on forward and reverse strands, respectively. The green and blue histograms represent the read distribution of ChIP and the input samples for 100-bp bins, respectively. The reference peak regions are highlighted in red. Note that the y axis indicates the read number normalized for the number of non-redundant reads, whereas the reference peak lists were identified based on raw read numbers. The gene reference was obtained from the UCSC genome browser [59]. (C) Visualization of the ChIP/Control enrichment distribution for 100-kb bins (chromosome 10). Bins with ChIP/control >1 are highlighted in red, and those with ChIP/control ≤1 are in gray. The GC contents and gene numbers for 500 kb windows are also plotted. The figures (B) and (C) were generated by DROMPA3.

Figure 2B shows the read distribution of each sample normalized for total non-redundant reads. The visualization of the read distribution is useful at the first step for judging the validity of the experimental results. In this figure, the active regions occupied by H3K36me3 and the silent regions occupied by H3K27me3 are distinguishable. On the other hand, H3K9me3 also has peaks, whereas their reads are not highly enriched. Considering H3K9me3 is a heterochromatin marker and generally not enriched in active gene regions [29], these peaks might be false positives. This is possibly because the mapped read number of H3K9me3 is more than twice that of the other samples (Figure 2A), resulting in a peak-calling threshold that is not stringent enough. In fact, because the S/N of the H3K9me3 antibody is smaller than those of H3K36me3 and H3K27me3, a high sequencing depth may be required to identify enriched regions. However, the chromosome-wide visualization of ChIP/input enrichment with a 100-kb bin clearly shows the exclusivity between H3K9me3 and H3K27me3 (Figure 2C). This large bin size enables us to directly use the ChIP/control distribution, even for large genomes, and provides an important insight.

Peak-calling

Establishing a definitive algorithm for peak detection has been the central topic in ChIP-seq analysis, resulting in the development of a plethora of programs. As the space is limited, we will just introduce here 20 representative programs.

Progress in the development of peak-calling algorithms

Peaks detection of a ChIP sample generally uses a corresponding input sample to estimate the background distribution at any genomic locus. Naked genomic DNA is less appropriate as a control because the input sample reflects the GC bias and chromatin structure rather than naked DNA [38, 41]. It was also reported that histone H3 ChIP-seq data can be used as a control [60]. ChIP samples for non-DNA-binding proteins, e.g. IgG, are often used to detect nonspecific binding sites. See [30] for a detailed discussion regarding control samples.

Early programs adopted the Poisson model, which assumes that the background reads are uniformly distributed along the genome (e.g. SICER [61] and CCAT [62]). However, a greater variation in the read distribution than allowed by the Poisson model is typically experimentally observed, and the negative binomial model that is an extension of the Poisson model was adopted to approximate such an overdispersion (e.g. CisGenome [63] and BayesPeak [64]). This model was extended to a zero-inflated negative binomial model to account for the zero-inflated read distribution caused by a lack of sequencing depth and low-mappable regions (e.g. MOSAiCS [65] and ZINBA [66]). For other strategies, MACS [39] uses the local Poisson model that estimates the parameter λ for each local genomic position. GPS [67] and PICS [68] predict protein-binding events using an EM algorithm. SISSRs [69], Peakzilla [70] and Q [12] focus on the equivalence between the read numbers of positive and negative strands to improve peak resolution. PePr [71] and JAMM [72] integrate information from multiple replicates to identify consistent or differential binding sites. The multiple hypothesis correction is performed to calculate false-discovery rates (FDRs) using the Benjamini–Hochberg procedure or the empirical method that calculates the peak number of the input sample compared with that of the ChIP sample.

Broad mode and mixed mode peaks depict weak and widespread enriched regions compared with the sharp mode, and there are no clear peak summits and sequence specificity. Although several peak-calling programs for the broad mode have been developed [61, 66, 73–75], and some peak-calling programs also have parameter settings for the broad mode, detection of such enrichment is still challenging. For proteins that are expected to be distributed within genic regions (e.g. Pol II and H3K36me3), a gene-annotation-based method is also useful, e.g. an aggregation plot around active genes and methods for differential gene expression analyses [76, 77]. For the characterization of Pol II occupancy, a travelling ratio (or pausing index) has been proposed [78]. The traveling ratio for gene i is defined as TR = , where d and d are the Pol II density in the promoter-proximal region and the gene body of gene i, respectively. This score indicates whether the promoter-proximal Pol II stalled at the gene. MUSIC discriminates between the binding modes, and between stalled and elongating forms of Pol II [74]. When investigating broad markers distributed in intergenic regions (e.g. H3K9me3), the gene-annotation-based method cannot be used. In such cases, it is still possible to use genome-wide visualization (Figure 2C) and compare the results with other public annotations, e.g. genome-wide maps of histone modifications [29].

While large genomes (e.g. human) require a statistical framework for peak calling owing to the low density and high variance of the reads, it is effective for small genomes (e.g. yeast [13, 79–81]) to inspect a genome-wide ChIP/input enrichment distribution itself (Figure 3).

Figure 3.

ChIP/input enrichment distribution of S. cerevisiae (chromosome I, 136–162 kb). Data from [11]. Smc6, Nse4 and ‘No tag (negative control)’ ChIP-seq data for a 100-bp bin with gene annotation obtained from the Saccharomyces Genome Database (http://www.yeastgenome.org). The reads were mapped onto the genome, allowing multiple mapped reads. For the yeast genome, inspecting the genome-wide ChIP/input enrichment distribution is effective because a read depth is large enough (>10-fold) and the division with the input sample can minimize the technical and biological biases of the conditions. The enriched regions of Smc6 and Nse4 that overlap those of the ‘No tag’ sample (black arrows) suggest false positives (e.g. hyper-ChIPable regions).

Which program is best for our analysis?

A performance comparison with a large set of programs in various aspects is of interest. However, several issues make such large-scale comparative study difficult. First, installing dozens of programs with multiple prerequired data (e.g. mappability) and assembling their results are often problematic owing to a large variety of existing file formats, and the availability and version-control of required tools in different computing environments. Second, because the programs have different underlying assumptions (e.g. using mappability and/or multiple mapped reads, subdividing peaks and necessity of replicates), the simultaneous comparison of them is difficult to be fair. Finally, the substantial evaluation of obtained peaks is difficult due to the lack of annotation for ‘true’ binding sites (see the next subsection). Although there are early studies of performance evaluations for peak-calling programs, based on the number, width and the distribution of peaks [82, 83], the reproducibility across replicates [38] and the accuracy against the known binding motif sites [84] or the manually curated benchmark data sets [85], the number of programs compared and quality metrics used in each study were limited. Moreover, the comparison with latest programs is lacked.

The appropriate method depends on the species, sample conditions and target proteins. Even though there is no clear consensus on which is best, the latest and widely used programs may be satisfactory for our needs.

How reliable are the obtained results?

Owing to the lack of annotation for true binding sites, the development of computational methods to evaluate the identified peaks has been limited. Motif-based evaluations are not applicable for proteins that do not have sequence specificity (e.g. histone modifications). Even for proteins with canonical motifs, there can be many tissue-specific binding sites recruited by other factors that do not involve the motif sequence [52, 86]. Another way of assessing the validity of the identified peaks is to focus on reproducibility. Cross-correlation coefficients of peak regions or whole-genomes measure the global similarity between two biological replicates [87]. The irreproducible discovery rate (IDR) assesses the rank consistency of common peaks between two replicates [88]. Based on a copula mixture model, IDR estimates the reproducibility of each peak pair, and reports the expected rate of irreproducible discoveries in the obtained peaks in a similar way to the FDR. In contrast to the CC coefficient, the IDR can assess each peak separately and therefore, can be used as a threshold robust for the technical variance owing to the analysis protocols and the choice of a peak-calling program. However, the poorest quality samples can become the bottleneck. When many true peaks do not appear in a replicate with a low S/N, they will be rejected as non-reproducible. Furthermore, several genomic regions tend to show artificially high enrichment levels, resulting in ‘reproducible’ false-positive peaks that cannot be filtered by the IDR. The ENCODE consortium summarizes empirically identified ‘blacklist regions’ for several species, which include repetitive and low-mappable regions [45]. Moreover, there are ‘hyper-ChIPable’ regions in the genome, in which peaks of unrelated proteins, including negative controls, overlap [79] (see Figure 3). These regions are positively correlated with promoters of well-expressed genes and are unchanged in cells containing the mutant protein of interest [89]. This result indicates that the peaks obtained by current methods may contain some (or a large) amount of unrecognizable false positives. Therefore, when investigating a protein with an unknown DNA-binding pattern, the obtained peaks should be validated carefully with a negative control (e.g. IgG), especially around transcription start sites (TSSs) of expressed genes.

How to treat low-quality samples?

There are various factors that can affect the data quality of the sample preparation step, including the quality of the antibody—e.g. its affinity and specificity, over-crosslinking, DNA fragmentation and overamplification by PCR and ChIP conditions. Different antibodies, even for the same protein (and even biological replicates of the same antibody), often produce completely different peak distributions. Despite an investigation of sequencing biases present in NGS data [90], it is still difficult to ascertain the exact sources of each bias in a sample preparation protocol [31]. In our experience, when a sample has a poor score for a quality measure, it often has other problems (e.g. low complexity causes a strong GC bias). As it is difficult to rescue poor-quality samples and include them in the analysis pipeline even with read normalization, fine-tuning sample preparation protocols to produce high-quality samples may be necessary for each project.

An efficient way to allow the use of ChIP-seq data of modest quality, while suppressing the noise, is by limiting the genomic regions to be investigated to a few candidate regions that satisfy the working hypothesis, and then validating them using other biological experiments, as suggested in [30]. This is a practical approach, rather than seeking the most accurate method by minutely tuning the parameters for each sample. Adding samples from related proteins and biological replicates increases reliability.

Normalization for a differential analysis

Relative-level difference (de novo normalization)

A typical procedure after obtaining peak sets is to summarize peak similarities and differences among samples in a binary (common or unique) or quantitative (variety of peak intensity) manner, which necessitates read normalization. The simplest normalization approach is to scale reads using the total read number (N) within the whole genome or background region, which assumes that the differences in the mapped reads among samples are small enough compared with the total read number. Instead of scaling reads using a constant factor, Taslim et al. [91] proposed the nonlinear method using a locally weighted regression (LOESS) to remove the effects of bias and systematic errors. However, the underlying assumption that the genome-wide distribution of read counts has an equal mean and variance across samples may not be valid in most cases (e.g. the different S/Ns between samples). Maehara et al. [92] proposed the co-localization score that measures the global similarity between two samples based on the common peaks, while it does not aim to identify individual differential peaks. Methods for differential gene expression analysis [76, 77] can be used to directly compare more than two groups, and these methods do not consider the different S/Ns among samples either. In contrast, MAnorm [93] and ChIPcomp [94] are designed to consider different S/Ns. MAnorm scales the reads of peaks common to two samples using a robust linear regression based on the MA plot. ChIPcomp performs quantitative comparison of multiple ChIP samples, which measures genomic background using control data and considers multiple-factor experimental designs. These tools assume that most of the common peaks should have similar binding intensities among the samples (e.g. same antibody for different conditions). See reference [94] for a detailed discussion regarding these de novo normalization programs.

Absolute-level difference (spike-in analysis)

In cases where the genome-wide peak distribution changes drastically (e.g. knock down analysis or stimulated versus non-stimulated) the aforementioned assumptions do not hold. Moreover, those normalizations are essentially limited to investigating relative differences in protein binding among samples [31]. For example, when the reads are relatively enriched in euchromatic regions and depressed in heterochromatic regions in a sample, it is difficult to discern whether protein binding was increased in the euchromatin or decreased in the heterochromatin in vivo. Recently, to investigate the absolute differences, several studies adopted spike-in analysis for human [95-97] and yeast [98, 99]. This analysis adds same quantities of chromatin DNA to all samples compared before or after immunoprecipitation. Because the number of reads derived from this reference chromatin should be same across samples, this number can be used as an internal control for read normalization. Thus, the spike-in analysis can detect global differences that cannot be identified by aforementioned de novo normalization methods. To discern reference reads from sample ones, the reference chromatin should be derived from a different genome.

Herein, we show an example of the spike-in data of H3K79me2 for EPZ5676-treated Jurkat cells from Orlando et al. (Figure 4). As shown in the original paper, spike-in normalization revealed a decrease in the H3K79me2 enrichment (Figure 4A and B). Next, the relative enrichment of H3K79me2 against that of 0% cells was visualized (Figure 4C). In the total read normalization for 100%/0% (fifth line), the read depth in the background (blue bars) was about four times greater owing to the genome-wide substantial decrease in H3K79me2-enriched regions (green bars). Using spike-in normalization (bottom line), the relative enrichment of the background declined to one, which indicates there is little difference within the background at an absolute level. A similar tendency can be observed for three other sample pairs, but they were less significant.

Figure 4.

Spike-in analysis of H3K79me2 ChIP-seq data for 0%, 25%, 50%, 75% and 100% EPZ5676-treated Jurkat cells. Data from [96] (GEO under accession number GSE60104). Spike-in normalization was implemented using the number of reads uniquely mapped onto the fly genome (UCSC dm3). (A) Read distribution near the RPL13A gene locus for 100-bp bins. Left: total read normalization, right: spike-in normalization. (B) Aggregation plots of total read normalization (left) and spike-in normalization (right) from 5-kb upstream to 10-kb downstream of the TSSs of the RefSeq genes. Shaded regions indicate a 95% confidence interval. (A) and (B) are identical visualizations of Figure 3C and E in reference [96], respectively. (C) Log-scale relative enrichment of H3K79me2 for 25%, 50%, 75% and 100% treated cells against 0% treated cells near the RPL13A gene locus (chromosome 19,49.88–50.13 Mb), with a 100-kb bin and 20-kb smoothing window. The top green line displays a H3K79me2 read distribution for 0% treated cells to roughly identify H3K79me2-enriched (green bars) and background regions (blue bars). Regions in which the enrichment (y-axis) is > 1 and <1 indicate a relative increase and decrease, respectively.

Although a spike-in analysis should be powerful and useful, its applicability and limitations, e.g. balancing the amount of spike-in relative to the chromatin of interest, are still not clear [31]. Moreover, we have occasionally observed a decrease in the read density owing to protein knockdowns in both the peak regions and in the background. It is challenging to determine whether a decrease in the read density of the background indicates that the protein of interest is also distributed in the background, and if we can conclude that the peak decreases owing to the knockdown when the read depth in the background decreases as much as in the peak regions. To answer these questions, extensive testing may be necessary, e.g. using a mock knockdown analysis. In addition, spiking genomes of multiple species may make the analysis more robust. The optimal normalization method should be chosen based on prior knowledge of the system and the statistics of the sequenced samples. The obtained results of a differential analysis should be evaluated using another method, e.g. quantitative PCR.

Integrative analysis for a de novo genome annotation

Although it is now possible to produce many ChIP-seq data sets at reasonable costs, their comparison and integration are not trivial. For example, when examining the differential ChIP-seq analysis of four proteins obtained under two conditions, and the knockdown effects on these data sets are of interest, it is necessary to investigate the differences among four proteins under the two conditions and/or between wild-type and knockdown cells simultaneously. As the results strongly depend on the peak-calling result of each sample [87], of which the protocol should be selected individually, the extensive work needed for tuning the protocol and integrating all results will be challenging. Consequently, there is a great demand for tools that jointly analyze all samples simultaneously.

Context-specific co-association identification

In the ENCODE consortium, Gerstein et al. [25] applied a machine-learning framework and examined the genome-wide co-association of 119 TFs that contained over 450 ChIP-seq samples. In this framework, peak calling is used with relaxed thresholds merely to obtain candidate regions for investigation. The identified sample peak sets are integrated into co-binding maps, and then ‘context-specific’ co-associations are identified, which are subsets of peaks binding to different TF sets in other genomic regions. These combinatorial patterns provide biological information on the high-dimensional interrelationship level for regulatory elements, which is difficult to ascertain using typical genome-wide pairwise comparisons. This concept has also been used for cross-species analyses of regulatory information [100-102]. Because this approach targeted the binding sites of point-source TFs, the results did not contain regions enriched in broad markers and background regions. The targeted regions may differ among experiments, and therefore it is difficult to directly compare the results across multiple experiments.

Joint analysis for a de novo genome annotation

Recently, integrative methods were developed to segment, classify and annotate a whole-genome sequence de novo, based on unsupervised machine-learning methods. These methods directly receive all ChIP sample data and analyze them simultaneously, instead of calling peaks and comparing them individually.

ChromHMM [103] and Segway [104] were developed to systematically identify the specific combination patterns of histone modifications as a chromatin state, which can detect large-scale variations of histone marks across the genome. ChromHMM is the most widely used tool, which models binary vectors (1 or 0) for each 200-bp bin converted from raw read counts using a sample-specific threshold as an independent Bernoulli random variable. Segway, by contrast, transforms the counts into real values and uses a dynamic Bayesian network at a 1-bp resolution. It can incorporate more complex relationships among samples in each region, although it requires a magnitude larger computational cost. Using these tools, high-quality chromatin-state maps for many cell lines are available [28, 29]. Several methods also exist that expand the methodologies of these tools to improve accuracy, computational cost or the interactive navigation [105-108]. Although this chromatin-state segmentation is not suitable for quantitative analysis among multiple experiments, several tools were developed to integrate and compare chromatin state sets from different experiments [109, 110]. hiHMM jointly infers chromatin state maps across multiple genomes and cell types [111].

There are also various joint analysis tools designed for TFs, based on several probabilistic models, such as the generalized EM algorithm and the Markov random field model [112-115]. These tools jointly model the dependencies among ChIP samples to identify global and local combinatorial enrichment patterns in whole-genome or specific functional regions. The number of classified enrichment patterns is not necessarily optimal for distinct functions of each protein, but these powerful approaches do efficiently reveal context-specific co-associations, and it may be possible to find out the false positives derived from hyper-ChIPable regions from total peak sets. Although the machine-learning approach for a large number of TFs is computationally daunting, these data-driven approaches may positively impact large-scale analyses in the near future.

Future potential: single-cell analysis

A limitation of ChIP-seq analysis is the requirement for large amounts of starting material (∼105 cells); accordingly, ChIP-seq analysis focuses on ensemble (averaged) features across large number of cells. To elucidate the internal heterogeneity within complex tissues and cell populations, the development of single-cell methodology is desired. In spite of rapid development of single-cell technologies in other genomic fields (reviewed in [116, 117]), there has been so far no report for single-cell ChIP-seq analysis. Recently, the first method for collection of chromatin data at single-cell resolution was published [118]. This ‘Drop-ChIP’ method adopts a droplet-based microfluidics system [119] for labeling chromatin from single cells before immunoprecipitation, which is also employed for single-cell RNA-seq analysis [120-122]. Labeled chromatins from all cells are sequenced, mapped and partitioned into single-cell reads by their barcode sequences. The single-cell profiles are classified by an unsupervised hierarchical clustering, and aggregated into profiles for subpopulations. The experiment using H3K4me2 antibody demonstrated that, although just a few hundred peaks were identified per cell owing to low sequencing depth (∼10 000 reads), this method could distinguish individual three cell types (mouse embryonic stem (ES) cells, embryonic fibroblasts and hematopoietic progenitors cells) with nearly 100% accuracy, and identify subpopulations in ES cells using the differences in chromatin signatures of pluripotency and differentiation priming.

This method involves several important considerations. First, because of low sensitivity and high technical read variance owing to low depth, this method can accept only antibodies for sharp mode peaks with high S/N so far. It is difficult to directly apply exiting techniques and quality metrics for traditional ChIP-seq to single-cell ChIP-seq data. Therefore, it is better to combine the result with typical ChIP-seq analysis, if possible. Second, this single-cell analysis aims to unveil internal heterogeneity by classifying patterns of subpopulations, rather than examine an individual cell (e.g. single-cell-level comparison with gene expression data). In this respect, this method is not interchangeable with low-input analyses [123-127], which aims to reduce the cell number required by a typical ChIP-seq workflow for precious biological samples (e.g. primary cells and clinical samples). Finally, identification of small subpopulation requires large number of sample cells [118]. Significance power of this method depends on the number of input cells and the performance of systems. Anyway, this first innovative challenge will encourage a new era in this field.

Concluding remarks

In this review, we discuss the computational aspects of ChIP-seq analysis and highlight key points associated with each step. We emphasize that the design of a ChIP-seq experiment is of critical importance and that a quality check of the data at each step is important, even when using published ChIP-seq data. By addressing a wide range of relevant subtopics within ChIP-seq analysis, which include integrative large-scale analysis and single-cell analysis, this review expands the topic and enhances the value of previous reviews. We believe that this review benefits researchers in related fields.

There are high-quality genome-wide maps of TFs and histone modifications provided by many consortia, and the open web servers (e.g. UCSC genome browser [59] and WashU Epigenome Browser [128]) enable researchers to effectively refer and track these resources for their own projects. Growth of publicly available resources will be a superb driving force for research progresses in the biological and bioinformatics field.

Remaining challenges for the future is to classify direct and indirect binding, capture temporary and non-site-specific (drifting) TF binding and investigate highly repetitive regions, e.g. centromeres. Finally, integration with other technologies, e.g. human genetic variation analyses [129-131], genome editing [132] and de novo assembly [133, 134], will make ChIP-seq analysis more fruitful and provide us with knowledge regarding the underlying mechanisms of genome function and evolution. Such comprehensive analyses facilitate the systematic elucidation of diverse biological activities intricately cooperating in the genome.