# NEXT-peak: a normal-exponential two-peak model for peak-calling in ChIP-seq data

- Nak-Kyeong Kim
^{1}Email author, - Rasika V Jayatillake
^{1}and - John L Spouge
^{2}

**14**:349

https://doi.org/10.1186/1471-2164-14-349

© Kim et al.; licensee BioMed Central Ltd. 2013

**Received: **25 January 2013

**Accepted: **20 May 2013

**Published: **25 May 2013

## Abstract

### Background

Chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) can locate transcription factor binding sites on genomic scale. Although many models and programs are available to call peaks, none has dominated its competition in comparison studies.

### Results

We propose a rigorous statistical model, the normal-exponential two-peak (NEXT-peak) model, which parallels the physical processes generating the empirical data, and which can naturally incorporate mappability information. The model therefore estimates total strength of binding (even if some binding locations do not map uniquely into a reference genome, effectively censoring them); it also assigns an error to an estimated binding location. The comparison study with existing programs on real ChIP-seq datasets (STAT1, NRSF, and ZNF143) demonstrates that the NEXT-peak model performs well both in calling peaks and locating them. The model also provides a goodness-of-fit test, to screen out spurious peaks and to infer multiple binding events in a region.

### Conclusions

The NEXT-peak program calls peaks on any test dataset about as accurately as any other, but provides unusual accuracy in the estimated location of the peaks it calls. NEXT-peak is based on rigorous statistics, so its model also provides a principled foundation for a more elaborate statistical analysis of ChIP-seq data.

## Keywords

## Background

ChIP-seq experiments use chromatin immunoprecipitation and then high-throughput sequencing, primarily to locate transcription factor binding sites across entire genomes, and to better our understanding of biological control systems [1]. As a brief overview of the relevant experimental protocols, they begin by irreversibly cross-linking a transcription factor (TF) molecule to its binding site in genomic DNA. They then shear the DNA into millions of short sequence fragments. Usually, the ends of a fragment are near the corresponding cross-link, but the exact distance between the end of the fragment and the cross-link is random. Moreover, the fragments on the two DNA strands show different systematic biases in the positions of their ends relative to the cross-link. Antibodies to the TF then precipitate each TF molecule along with its attached fragment. Fragments are dissociated from the TF molecules, amplified by polymerase chain reaction (PCR). The fragments are then sequenced into short subsequences, called “tags”. Computational analysis then enters by mapping the tag sequences to a reference genome. If a tag sequence is long enough, the tag matches only one genomic coordinate. Sometimes, however, its sequence is short and maps to more than one coordinate, making the mapping ambiguous. The possibilities of ambiguous mapping, false positive tag-reads, and other experimental errors have motivated the development of programs to analyze ChIP-seq experiments.

Peak-calling programs locate potential binding sites as “peaks” where mapped tags concentrate. Peak-calling programs use widely differing approaches, none of which has yet emerged as dominant in reviews [2–4], because relative accuracy of programs varies with the dataset examined [2, 3]. Improvement is probably possible, however, because the models underlying existing programs do not consider mapping ambiguities directly, despite the existence of packages for enumerating the ambiguities, e.g., the PeakSeq suite [5]. Moreover, some programs ignore strand-specific information [6, 7].

Most programs use (sometimes implicitly) kernel smoothing to compensate for mapping ambiguities, the most popular kernel being the uniform density, which is equivalent to counting tags in sliding windows of fixed-length [8–12]. Programs also manipulate information about a tag’s strand in various ways: as mentioned, some ignore it [6, 7]; some use it explicitly (e.g. SISSRs [12], spp [11]); and some use it indirectly by transposing tag locations over to the other strand (e.g. MACS [8], QuEST [13], CisGenome [9]). Programs that combine strand information and windowing are essentially using two separate uniform densities as kernels for the forward and backward strands. For example, QuEST [13] (among other computer programs) estimates tag densities with a kernel smoother than the uniform density, to mimic the observed shape of tag peaks. QuEST did not dominate in comparisons, however, perhaps because it transposes tag locations, rather than estimating two separate tag densities, one for each strand.

The significance of peaks can be reported either as the number of tags in a window, a p-value, a q-value, or a posterior probability. Although p-values guide a naïve user better than tag numbers, they introduce problematic assumptions. To derive a p-value, some programs [7, 13] assume a globally uniform background intensity of tags, an assumption known to fail in ChIP experiments. To assign the significance of peaks, different programs use various model assumptions such as Poisson [7, 12, 13], local Poisson [8], binomial [9], and hidden Markov [10] models. Some programs [8, 9, 11–13] use control data to account for local variations in background tag intensity or to compensate for experimental artifacts like PCR over-amplification, which can cause a spurious concentration of tags in a few specific locations. The reproducibility of control data is suspect, however, because it varies across cell types and ChIP protocol [4]. Although control data mitigates some experimental artifacts, its unreliability ultimately undermines any inference based on the corresponding p-value.

Mapping ambiguities can also be problematic for a naïve p-value calculation. Consider, e.g., a large genomic region where exactly *A* locations are ambiguously mapped (where *A* is fixed). In the same region, consider a window containing a total of *L* locations, including the *A* ambiguous locations. The window length is an arbitrary parameter (within limits), and as it increases, *L* increases. The *A* ambiguous locations are essentially censored data, so the simplest maximum likelihood estimate of the tag count in the window is *L*/(*L-A*) times the observed tag count. Thus, if observed tag count is fixed, the estimated tag count decreases as the window length increases. If a p-value decreases with the estimated tag count, it then depends on the window length. False discovery rates (FDR) depend on p-values, so the use of FDRs does not remove the dependency. Under the circumstances described, therefore, the arbitrary choice of window length influences the number of peaks reported. A recent study on the number of binding sites in a genome [14] indicates that many real binding sites from ChIP-seq data go unreported, suggesting that the assumptions and approximations underlying current p-value estimates leave room for improvement.

Intuitively, the spatial resolution of a peak should also improve as more tags contribute to it. In principle, therefore, a program should also assign errors to its location estimates, but in fact, existing programs do not infer the accuracy of their estimated peak locations.

**Summary of programs used for comparison**

Profile | Strand specific | Statistical model | Peak criteria | Rank by | Ref | |
---|---|---|---|---|---|---|

HPeak | Sliding window | Indirect | Hidden Markov model | Tag counts | Tag counts | [10] |

WTD | Sliding window | Direct | Score | P-value | [11] | |

MTC | Sliding window | Direct | Score | P-value | [11] | |

CisGenome | Sliding window | Indirect | Binomial | Tag counts | Tag counts | [9] |

MACS | Sliding window | Indirect | Local Poisson | P-value | P-value | [8] |

QuEST | Gaussian kernel density | Indirect | Poisson | Height threshold | Q-value | [13] |

SISSRs | Sliding window | Direct | Tag counts with sign change | P-value | [12] | |

NEXT-peak | NEXT-peak kernel density | Direct | Poisson per base pair | Likelihood | Estimated binding intensity (2ν) |

## Results

### Fitting the NEXT-peak model to ChIP-seq data

Using ChIP-seq datasets for three TFs: STAT1 [15], NRSF, and ZNF143 (see Methods for details), results with and without mappability information were examined; for each dataset, only the better of the two are presented here. Mappability information improved results for STAT1, but degraded results for NRSF and ZNF143.

^{-6}for all datasets. See Methods for details on the p-value computation for finding motif sites. Figure 1a shows a density of the normal-exponential two-peak (NEXT-peak) model (see Methods). Figure 1b-d displays the tag number, normalized to a probability density, for each location around the position of the candidate sites. The observed tag density is superimposed on the estimated density (derived from model estimates of the expected tag numbers, ${\lambda}_{j}^{R}$ or ${\lambda}_{j}^{L}$ in the Methods). Maximum likelihood estimation on the NEXT-peak model produced parameter estimates underlying ${\lambda}_{j}^{R}$ and ${\lambda}_{j}^{L}$. Table 2 reports estimated parameter values for each dataset.

**Summary of ChIP-seq datasets**

Dataset | tag length | motif length | number of tags | $\widehat{\mathit{\beta}}$ | $\widehat{\mathit{\sigma}}$ | Genome |
---|---|---|---|---|---|---|

STAT1 | 27 | 15 | 15.1 million | 73.5 | 52.8 | Human build 36.1 |

NRSF | 36 | 21 | 33.1 million | 30.4 | 23.6 | Human build 36.1 |

ZNF143 | 36 | 20 | 25.2 million | 35.3 | 21.6 | Human build 36.1 |

For STAT1, the observed tag counts follow the density curve of the NEXT-peak model with a small difference in terms of average trend, except for two unexplained dips (Figure 1b). The dips (one for right tags, and one for left tags) display symmetry around the binding site. The tag counts for NRSF (Figure 1c) also follow the NEXT-peak density with a small trend difference. The tag counts for ZNF143 (Figure 1d) show a larger trend difference from the NEXT-peak density, perhaps because the ChIP experiment was noisier.

### Examples of regions with unmappable locations

### Comparison of the programs: top 2,000 peaks

To compare peak-calling programs, we run NEXT-peak along with other popular programs like HPeak [10], spp package (WTD and MTC) [11], CisGenome [9], MACS [8], QuEST [13], and SISSRs [12] (see Background for details). On running these programs including NEXT-peak, following standard practice, we used the default parameters to ensure reproducibility. As the first stage of comparison, we consider the 2,000 top peaks from each. The NEXT-peak program uses the estimated tags per binding event ($\widehat{\nu}$) to rank its peaks. We considered every peak called within 250 bp of a candidate site (as determined by position-specific scoring matrix search) to be a true positive (TP). Our primary performance measure was the number of TPs within the 2,000 highest peaks. (Precision provides a standard but equivalent performance measure: the precision at 2,000 positives is the number of TPs within the 2,000 top peaks divided by the constant, 2,000.) Our secondary performance measures considered placement of TP peaks: (1) the mean distance between a TP peak center and the nearest motif site, and (2) the mean bias between a TP peak center and the nearest motif site. A TP peak upstream of the nearest motif site contributes to a negative bias; downstream, a positive bias. Thus, small distances and biases are desirable.

**Program result summary for ChIP-seq datasets**

STAT1 | NRSF | ZNF143 | |||||||
---|---|---|---|---|---|---|---|---|---|

Program | #motif sites | mean distance | mean bias | #motif sites | mean distance | mean bias | #motif sites | mean distance | mean bias |

HPeak | 726 | 23.6 | 12.0 | 1489 | 7.5 | 0.2 | 698 | 26.0 | 0.5 |

WTD | 732 | 20.2 | 0.7 | 1487 | 17.4 | −16.5 | 698 | 26.1 | −13.4 |

MTC | 737 | 18.8 | 1.3 | 1498 | 16.7 | −15.8 | 697 | 25.7 | −12.2 |

CisGenome | 716 | 28.5 | 11.3 | 1471 | 8.0 | 0.3 | 681 | 27.6 | −15.1 |

MACS | 740 | 24.4 | −0.9 | 1491 | 18.5 | −17.4 | 709 | 27.2 | −13.8 |

QuEST | 715 | 24.9 | 11.2 | 1105 | 8.3 | −0.4 | 703 | 22.0 | 1.2 |

SISSRs | 659 | 19.6 | −0.2 | 1413 | 18.1 | −16.8 | 661 | 29.1 | −14.3 |

NEXT-peak | 781 | 16.7 | −1.6 | 1507 | 6.0 | −0.5 | 707 | 20.3 | 0.6 |

### Comparison of the programs: top peaks in general

*r*in the x-axis, the precision is computed for cumulative peaks between rank 1 and rank

*r*. As expected, precision generally decreased with the length of the list. For STAT1, NEXT-peak had the largest precision (of all programs examined) over the full range of lengths, up to 10,000 peaks. For NRSF, NEXT-peak had nearly the best precision up to 4,500 peaks (and in fact, it had the best precision at 2,000 peak; see the previous section.); NEXT-peak had the best precision between 4,500 and 10,000 peaks. For ZNF143, NEXT-peak had near the best precision up and 10,000 peaks. For ZNF143, MACS performed similarly to NEXT-peak between 1,500 and 10,000 peaks, but MACS had a significantly poorer performance between 0 and 1,500 peaks compared to other programs.

### Correlation of estimated standard deviation and distance to motif site

## Discussion

Alone among existing peak-calling programs, NEXT-peak analyzes data with a parametric statistical model. Of the existing programs, therefore, it alone provides a principled foundation for elaborating the statistical analyses of ChIP-seq data. One obvious elaboration is to model multiple binding events in a region. This work is currently underway and the results will be reported elsewhere.

NEXT-peak can estimate the average fragment length, even if the experiment does not measure the average fragment length. Let *d* denote the tag length, e.g., *d*=27 for STAT1. In the NEXT-peak model, the average distance from a fragment end to a cross-link is *β*. The average distance between the fragment ends is therefore 2*β* + *d* − 1 (because the “location” of the right end is the leftmost position of the corresponding tag). For the STAT1 dataset, $\widehat{\beta}$ =73.5 and *d* =27, so the NEXT-peak model estimate of the average fragment length is 173.0, consistent with a previous estimate of 174 [15]. For NRSF data, $\widehat{\beta}$ =30.4 and *d*=36, the estimate of the fragment length is 95.8, and for GABP data, $\widehat{\beta}$ =35.3 and *d*=36, the estimate is 105.6.

Existing programs simply discard ambiguously mapped reads. In contrast, NEXT-peak explicitly models the locations where reads do not map uniquely into the reference genome. NEXT-peak can therefore adjust for ambiguous mapping while estimating the total number of tags in a region, thereby sharpening its estimates of TF binding strength. Sharper estimates of binding strength can promote better physical interpretation of ChIP-seq results.

Existing peak-calling programs require tedious visual screening of up to tens of thousands of binding regions, to eliminate experimental artifacts like spikes in tag numbers caused by PCR anomalies. The goodness-of-fit tests in NEXT-peak can reduce the burden of visual screening. Moreover, the same tests can detect the presence of multiple TF binding sites, which are usually found in regions longer than those containing PCR anomalies. The long regions with small p-values can therefore be set aside for further, more intensive analyses, such as searching for multiple binding events or sequence motifs.

For STAT1, the observed tag distribution follows the NEXT-peak density closely, indicating that the NEXT-peak model captured the essence of the physical processes in the ChIP-seq experiment. Consequently, NEXT-peak outperformed its competitors, possibly because the NEXT-peak model successfully mimicked the true experimental kernel. On the other hand, for ZNF143, the observed tag distribution is somewhat deviated from the NEXT-peak density, possibly degrading NEXT-peak’s performance slightly. The observed tag density might reflect a mixture of multiple binding events, however, resulting from TF binding fluctuating between different protein complexes. Mass and structural differences between the protein complexes could cause binding locations or mean fragment lengths to fluctuate. Conventional motif analysis or a more elaborate model including multiple binding sites might expose the protein-protein interactions, however.

By adding mappability information, STAT1 increased true-positive binding sites by 4.0% on average. Unlike STAT1, mappability information for NRSF and ZNF143 actually degraded the performance of NEXT-peak: on average, it decreased true-positives by 0.5% and 0.6%, a surprising result given that both NRSF and ZNF143 had large numbers of mapped tags (33.1 millions and 25.2 millions). The truncated read length for NRSF and ZNF143 was 36, however, much larger than read length of 27 for STAT1. Thus, fewer genomic locations were mapped ambiguously for NRSF or ZNF143 (10.3%) than for STAT1 (16.2%), diminishing NEXT-peak’s ability to enhance its performance by adding mappability information.

This article examined three ChIP-seq datasets with a single dominant binding motif, permitting motif sites to serve as surrogates for the true binding sites. In general, however, even with an antibody specific to a protein, protein-protein interactions between TF molecules might cause multiple TFs (and hence, multiple motifs) to cross-link to an antibody. The two global parameters *σ* (the standard deviation for the cross-link locations) and *β* (the intensity of the Poisson process modeling shearing) then require delicate estimation. One could select a few hundred of the most tag-rich regions. One could screen the regions visually, choosing the ones with a good fit to the dual normal-exponential density and then estimate *σ* and *β*. Alternatively, one could perform a motif search on the tag-rich regions. The observed tag density for each motif then can be fit to the NEXT-peak model. Thus, NEXT-peak can analyze any ChIP-seq experiment, even without specific information on the protein interactions.

## Conclusions

We proposed a new statistical model for identifying binding sites from ChIP-seq data. The model successfully mimics the underlying data-generating process in ChIP-seq experiments by using the dual density of a normal-exponential two-peak model. The NEXT-peak program produced better prediction with more true positives and a better spatial resolution than any other program tested. The NEXT-peak program tests the validity of its underlying NEXT-peak model without depending (as many programs do) on an unrealistic assumption of a global uniform background tag distribution. The NEXT-peak program stands alone in quantifying errors by reporting a standard error for its estimates of binding intensity. Moreover, smooth scatterplots showed that its standard errors are informative about errors in motif location, as estimated from external standards. The NEXT-peak program also provides a goodness-of-fit test, automating screening of the spurious binding, and work is in progress to extend its model to locate multiple binding events in a region.

## Methods

### ChIP-seq datasets

Our analysis used ChIP-seq datasets corresponding to three TFs: STAT1, NRSF, and ZNF143. Because the three datasets correspond to known binding motifs, they provide a gold-standard for evaluating peak-calling programs [2]. Table 2 presents summary statistics for the three datasets.

The STAT1 dataset [15] was downloaded at http://www.bcgsc.ca/downloads/chiptf/human/STAT1/stimulated/july_23_2008/. The NRSF[SRA:SRR577995] and ZNF143[SRA:SRR243553] datasets were downloaded from the SRA database at http://www.ncbi.nlm.nih.gov/sra. Bowtie [17] mapped tags into a reference human genome (NCBI Build 36.1) for all three datasets. Mismatches of up to 2 bases were permitted, if they mapped uniquely within the genome.

For all datasets, the PeakSeq suite [5] then determined whether tag sequences in the reference genome were unique. PeakSeq requires tags of a uniform length. The tags for the downloaded STAT1 dataset, however, had varying length although most tags had length 27. We truncated the tags to length 27, if they were longer, or we discarded them, if they were shorter. Thus, we used the mappability information for 27 bp tags to approximate the complete STAT1 data. The downloaded NRSF had tags of length 50 and the downloaded ZNF143 had tags of length 40. We truncated tags from both datasets to length 36 to make it easy to investigate the effect of the mappability information (36 bp mappability information was used). Additional file 1 also reports on results from three additional datasets, MAX, GABP, and FoxA1.

### Notations for the ChIP-seq data

Let some preliminary method (see NEXT-peak algorithm for detail) flag possible cross-links in candidate genomic regions *R*_{
r
} (*r* = 1, …, *S*). Computational time is a consideration, because *S* can be on order of 10^{4} or more. Consider a specific genomic region *R*_{
s
}, where *s* ∈ {1, …, *S*}, and let *R*_{
s
} have width *w*_{
s
}, with the nucleotide bases having coordinates 1, …, *w*_{
s
}. Call the forward and backward DNA strands “left” and “right”, so the bases at each location *j* on the left and right strands are complementary (*j* = 1, …, *w*_{
s
}). Let ${\mathbb{X}}^{0}$ denote the set of locations within *R*_{
s
} where a tag sequence is not unique within the genome, so the corresponding tag maps ambiguously.

The superscripts “*L*” and “*R*” distinguish quantities pertaining to the left and right strands: note in particular that “*L*” and “*R*” are not exponents. Within *R*_{
s
}, let a total of *n*^{
L
} “left tags” be observed on the left strand; *n*^{
R
} “right tags”, on the right strand. Define the “location” of a tag as its leftmost position. Let the location of the left tags map to ${x}_{i}^{L}\in \left\{1,\dots ,{w}_{s}\right\}$ (*i* = 1, …, *n*^{
L
}); the location of the right tags, to ${x}_{i}^{R}\in \left\{1,\dots ,{w}_{s}\right\}$ (*i* = 1, …, *n*^{
R
}). Thus, the data are $\mathbf{X}=\left({x}_{1}^{L},\dots ,{x}_{{n}^{L}}^{L},{x}_{1}^{R},\dots ,{x}_{{n}^{R}}^{R},{\mathbb{X}}^{0}\right)$ where no location ${x}_{i}^{L}$ or ${x}_{i}^{R}$ is in ${\mathbb{X}}^{0}$.

An alternative representation is occasionally useful. Let ${y}_{j}^{L}\in \left\{0,1,2,\dots \right\}$ be the number of left tags observed at location *j* ∈ {1, …, *w*_{
s
}}; ${y}_{j}^{R}$, the number of right tags observed at position *j*. Tags cannot be observed at locations $j\in {\mathbb{X}}^{0}$. Thus, $\left({y}_{1}^{L},\dots ,{y}_{w}^{L},{y}_{1}^{R},\dots ,{y}_{w}^{R},{\mathbb{X}}^{0}\right)$ provides an equivalent representation of the data **X**. The model parameters for the data in *R*_{
s
} are (*μ*_{
s
}, *ν*_{
s
}, *ρ*_{
s
}) and (*σ*^{2}, *β*), defined below. Parameters specific to the region *R*_{
s
} are subscripted with “*s*”; the parameters common to all genomic regions lack subscripts.

### Dual density for a binding event

*ϕ*(•) and the cumulative distribution function

*Φ*(•). Each observed right tag location ${x}_{i}^{R}$ in

*R*

_{ s }corresponds to an underlying (and unobservable) random variate

*ξ*

_{ i }, the coordinate of the corresponding cross-link. For mathematical convenience, assume

*ξ*

_{ i }~

*N*(

*μ*

_{ s },

*σ*

^{2}), i.e.,

*ξ*

_{ i }is a normal variate with mean

*μ*

_{ s }(specific to

*R*

_{ s }) and variance

*σ*

^{2}(common to

*R*

_{ r }for

*r*= 1, …,

*S*). The density of the distribution of

*ξ*

_{ i }is

*β*of the exponential distribution is common to all regions

*R*

_{ r }(

*r*= 1, …,

*S*), so the density function corresponding to ${x}_{i}^{R}$ is

for ${x}_{i}^{R}>{\xi}_{i}$, and 0 otherwise.

*ξ*

_{ i }out, to derive the density function of ${x}_{i}^{R}$:

*σ*and

*β*in a ChIP-seq experiment. Figure 1a shows a normal-exponential density for both left and right tags (as discussed in Results).

### Poisson regression model for the observed tags

*R*

_{ s }: (1) let

*ν*

_{ s }be the expected number of right tags for each TF molecule that binds; (2) let

*ρ*

_{ s }be the uniform background intensity of right tags; and (3) let ${\lambda}_{j}^{R}$ be the expected number of right tags at location

*j*. Assume

Approximate the sum over all locations *j* with an integration, to produce a consistency check: *ν*_{
s
} = ∫ *ν*_{
s
} · *f*^{
R
}(*x*|**θ**)*dx*. The expected total number of right tags within *R*_{
s
} due to binding is therefore approximately *ν*_{
s
}, *ν*_{
s
} = 0 being equivalent to the absence of TF binding in *R*_{
s
}. On the other hand, the expected total number of right tags within *R*_{
s
} due to noise is *w*_{
s
}*ρ*_{
s
}.

*j*, and assume

*Y*

_{ j }

^{ R }has a Poisson distribution with mean ${\lambda}_{j}^{R}$, i.e.,

for *y* ϵ {0, 1, 2, …}.

*μ*

_{ s },

*σ*

^{2},

*β*,

*ν*

_{ s },

*ρ*

_{ s }) and differ only in mirroring the tag location densities around the link location

*μ*

_{ s }, i.e.,

The models for the left and right tags share *ν*_{
s
} and *ρ*_{
s
}, so as in the model for the right tags, *ν*_{
s
} is the expected number of left tags for each TF molecule that binds, and *ρ*_{
s
} is the uniform background intensity of left tags.

Figure 1a shows a normal-exponential two-peak (NEXT-peak) density. In this example, *μ* = 0, *β* = 60, and *σ* = 40. The two density curves mirror each other around the center location *μ* = 0. The curves are asymmetrical distributions skewed in opposite directions. Although both tails of each density rapidly approach zero, one tail approaches zero much faster than the other. Figure 1a motivates the model name: the “normal-exponential two-peak” (NEXT-peak) model.

**θ**

_{ s }= (

*μ*

_{ s },

*ν*

_{ s },

*ρ*

_{ s }). Under the NEXT-peak model in Figure 1a, the likelihood of the data $\left({y}_{1}^{L},\dots ,{y}_{w}^{L},{y}_{1}^{R},\dots ,{y}_{w}^{R},{\mathbb{X}}^{0}\right)$ in

*R*

_{ s }is

where **θ** = (**θ**_{
r
} : *r* = 1, …, *S*).

Because *σ* and *β* do not depend on the region *R*_{
s
}, training data can yield maximum likelihood estimates (MLEs) $\widehat{\sigma}$ and $\widehat{\beta}$. Fix the values $\sigma =\widehat{\sigma}$ and $\beta =\widehat{\beta}$, so the remaining parameters requiring estimation for the NEXT-peak model are **θ** = (**θ**_{1}, …, **θ**_{
S
}). Maximization of the profile likelihood $L\left({\mathbf{\theta}}_{s}\right)=L\left(\widehat{\sigma},\widehat{\beta},{\mu}_{s},{\nu}_{s},{\rho}_{s}\right)$ then yields the estimate ${\widehat{\theta}}_{s}=\left({\widehat{\mu}}_{s},{\widehat{\nu}}_{s},{\widehat{\rho}}_{s}\right)$ within each region *R*_{
s
}. The estimate’s components are: (1) the estimated mean location ${\widehat{\mu}}_{s}$ of a binding event, (2) the estimated mean total number ${\widehat{\nu}}_{s}$ of right (or left) tags within *R*_{
s
} due to the binding event, and (3) the estimated uniform background intensity ${\widehat{\rho}}_{s}$ of right (or left) tags within *R*_{
s
}. As usual, the inverse of the Fisher Information matrix (i.e., the inverse of the negative expectation of the Hessian of the log-likelihood) estimates the asymptotic variance-covariance matrix for **θ**_{
s
}.

**θ**

_{ s }= (

*μ*

_{ s },

*ν*

_{ s },

*ρ*

_{ s }) under the restriction

*ν*

_{ s }= 0. To test whether or not a binding event occurred in

*R*

_{ s }, under the null hypothesis

*H*

_{0}:

*ν*

_{ s }= 0, the likelihood ratio (LR) statistic

has an asymptotic chi-square distribution with 1 degree of freedom. Unlike the null hypotheses in many existing programs, the NEXT-peak model does not assume a globally uniform background intensity. Instead, its locally uniform background intensity is equivalent to assuming that the expected number of background tags per location varies slowly enough to be almost constant within each region *R*_{
r
} (*r* = 1, …, *S*).

### Goodness-of-fit test

The following can test whether observed tag data are consistent with the NEXT-peak model. PCR over-amplification and other experimental artifacts can cause spikes in the observed tag distribution. Likewise, multiple binding events in a region *R*_{
r
} cause deviations from the unimodal density of the model. Accordingly, consider Pr (**D**|*H*_{0}), the probability of the data **D** under the null hypothesis *H*_{0} of a NEXT-peak model. Under *H*_{0}, the counts ${y}_{j}^{R}$ of right tags at location *j* are generated with the Poisson intensity ${\lambda}_{j}^{R}$ given above, the counts of the left tags being generated with the mirror intensity ${\lambda}_{j}^{L}$. Consider also the probability Pr (**D**|*H*_{1}) of the data under an alternative model *H*_{1} where the underlying intensity ${\lambda}_{j}^{R}$ is unrestricted. The LR statistic 2 log[Pr(**D**|*H*_{1})/Pr(**D**|*H*_{0})] for the models follows a *χ*^{2} distribution, with degrees of freedom equal to the difference of the number of parameters in *H*_{0} and *H*_{1}. The LR test yields a p-value for each region *R*_{
s
}, with a small p-value indicating a poor fit within *R*_{
s
} to the NEXT-peak model underlying *H*_{0}.

### P-value computation for finding motif sites

From a position-specific score matrix, any segment receives a score by adding the corresponding columns scores from the position-specific score matrix. The probability of observing the score or higher is computed based on the null model that nucleotides (A, C, G, and T) can appear at random with equal probabilities. We used the Staden’s method [18] for the convolution computation of the score distribution, The principal of our p-value computation has been used for PSI-BLAST [19] and A-GLAM [20], among others, for their p-value computation, where p-values are used to compute E-values.

### Screening based on goodness-of-fit tests

Like many peak-calling programs, NEXT-peak masks locations having anomalously many tags, but it does so in a principled manner, with p-values based on goodness-of-fit tests. Long regions with small p-values suggest multiple binding sites, so for all datasets NEXT-peak masked regions less than a certain length long and with small p-values. When motif site locations are available, based on the area under the precision curve up to the top 10,000 peaks (e.g. Figure 3), the NEXT-peak program automatically reports a cut-off recommendation for each dataset. The Results section uses the length 400 and the p-value 10^{-8} as a cut-off for STAT1 and the length 300 and the p-value 10^{-4} as a cut-off for NRSF and the length 400 and the p-value 10^{-6} as a cut-off for ZNF143, all suggested by the NEXT-peak output.

### NEXT-peak algorithm

*σ*and

*β*by maximizing the likelihood using motif site locations. For a known TF, a publicly available motif pattern is used, e.g. from JASPAR [16]. For an unknown motif, run the NEXT-peak program with default values (

*σ*= 30 and

*β*= 50), and identify the strongest motif from a motif search. Alternatively, a user can estimate these parameters with selected regions. (4) For each region, estimate

*μ*and

*ν*by maximizing the likelihood. It computes the standard deviation estimates for these estimates. Then, perform a goodness-of-fit test for each region. (5) As a post-processing step, compute the region length and p-value cut-off recommendations to screen out potential spurious regions when the motif site locations are available.

### NEXT-peak software

The NEXT-peak program is implemented in C++. The code is publicly available at http://www.odu.edu/~nxkim/nextpeak/. A typical computation time is about 15 ~ 45 minutes, depending on the size of the input data.

## Declarations

### Acknowledgements

This work was supported by the Summer Research Fellowship Program of Old Dominion University Research Foundation and the Intramural Research Program of the National Library of Medicine at the National Institutes of Health.

## Authors’ Affiliations

## References

- Johnson DS, Mortazavi A, Myers RM, Wold B: Genome-wide mapping of in vivo protein-DNA interactions. Science. 2007, 316 (5830): 1497-1502. 10.1126/science.1141319.View ArticlePubMedGoogle Scholar
- Wilbanks EG, Facciotti MT: Evaluation of algorithm performance in ChIP-seq peak detection. PLoS One. 2010, 5 (7): e11471-10.1371/journal.pone.0011471.PubMed CentralView ArticlePubMedGoogle Scholar
- Laajala TD, Raghav S, Tuomela S, Lahesmaa R, Aittokallio T, Elo LL: A practical comparison of methods for detecting transcription factor binding sites in ChIP-seq experiments. BMC Genomics. 2009, 10: 618-10.1186/1471-2164-10-618.PubMed CentralView ArticlePubMedGoogle Scholar
- Pepke S, Wold B, Mortazavi A: Computation for ChIP-seq and RNA-seq studies. Nat Methods. 2009, 6 (11 Suppl): S22-S32.PubMed CentralView ArticlePubMedGoogle Scholar
- Rozowsky J, Euskirchen G, Auerbach RK, Zhang ZD, Gibson T, Bjornson R, Carriero N, Snyder M, Gerstein MB: PeakSeq enables systematic scoring of ChIP-seq experiments relative to controls. Nat Biotechnol. 2009, 27 (1): 66-75. 10.1038/nbt.1518.PubMed CentralView ArticlePubMedGoogle Scholar
- Fejes AP, Robertson G, Bilenky M, Varhol R, Bainbridge M, Jones SJ: FindPeaks 3.1: a tool for identifying areas of enrichment from massively parallel short-read sequencing technology. Bioinformatics. 2008, 24 (15): 1729-1730. 10.1093/bioinformatics/btn305.PubMed CentralView ArticlePubMedGoogle Scholar
- Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5 (7): 621-628. 10.1038/nmeth.1226.View ArticlePubMedGoogle Scholar
- Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nussbaum C, Myers RM, Brown M, Li W: Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008, 9 (9): R137-10.1186/gb-2008-9-9-r137.PubMed CentralView ArticlePubMedGoogle Scholar
- Ji H, Jiang H, Ma W, Johnson DS, Myers RM, Wong WH: An integrated software system for analyzing ChIP-chip and ChIP-seq data. Nat Biotechnol. 2008, 26 (11): 1293-1300. 10.1038/nbt.1505.PubMed CentralView ArticlePubMedGoogle Scholar
- Qin ZS, Yu J, Shen J, Maher CA, Hu M, Kalyana-Sundaram S, Chinnaiyan AM: HPeak: an HMM-based algorithm for defining read-enriched regions in ChIP-Seq data. BMC Bioinforma. 2010, 11: 369-10.1186/1471-2105-11-369.View ArticleGoogle Scholar
- Kharchenko PV, Tolstorukov MY, Park PJ: Design and analysis of ChIP-seq experiments for DNA-binding proteins. Nat Biotechnol. 2008, 26 (12): 1351-1359. 10.1038/nbt.1508.PubMed CentralView ArticlePubMedGoogle Scholar
- Jothi R, Cuddapah S, Barski A, Cui K, Zhao K: Genome-wide identification of in vivo protein-DNA binding sites from ChIP-Seq data. Nucleic Acids Res. 2008, 36 (16): 5221-5231. 10.1093/nar/gkn488.PubMed CentralView ArticlePubMedGoogle Scholar
- Valouev A, Johnson DS, Sundquist A, Medina C, Anton E, Batzoglou S, Myers RM, Sidow A: Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data. Nat Methods. 2008, 5 (9): 829-834. 10.1038/nmeth.1246.PubMed CentralView ArticlePubMedGoogle Scholar
- Kuznetsov VA, Singh O, Jenjaroenpun P: Statistics of protein-DNA binding and the total number of binding sites for a transcription factor in the mammalian genome. BMC Genomics. 2010, 11 (Suppl 1): S12-10.1186/1471-2164-11-S1-S12.PubMed CentralView ArticlePubMedGoogle Scholar
- Robertson G, Hirst M, Bainbridge M, Bilenky M, Zhao Y, Zeng T, Euskirchen G, Bernier B, Varhol R, Delaney A: Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing. Nat Methods. 2007, 4 (8): 651-657. 10.1038/nmeth1068.View ArticlePubMedGoogle Scholar
- Sandelin A, Alkema W, Engstrom P, Wasserman WW, Lenhard B: JASPAR: an open-access database for eukaryotic transcription factor binding profiles. Nucleic Acids Res. 2004, 32 (Database issue): D91-D94.PubMed CentralView ArticlePubMedGoogle Scholar
- Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10 (3): R25-10.1186/gb-2009-10-3-r25.PubMed CentralView ArticlePubMedGoogle Scholar
- Staden R: Methods for calculating the probabilities of finding patterns in sequences. Comput Appl Biosci. 1989, 5 (2): 89-96.PubMedGoogle Scholar
- Altschul SF, Koonin EV: Iterated profile searches with PSI-BLAST–a tool for discovery in protein databases. Trends Biochem Sci. 1998, 23 (11): 444-447. 10.1016/S0968-0004(98)01298-5.View ArticlePubMedGoogle Scholar
- Tharakaraman K, Marino-Ramirez L, Sheetlin S, Landsman D, Spouge JL: Alignments anchored on genomic landmarks can aid in the identification of regulatory elements. Bioinformatics. 2005, 21: I440-I448. 10.1093/bioinformatics/bti1028.PubMed CentralView ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.