Lecture notes. Written from the material of MIT 7.091J, licensed CC BY-NC-SA 4.0. These are notes, not a transcript: the material has been reorganised and rewritten. This adaptation carries the same licence, and the original is linked above.

7. ChIP-seq Peak Calling and IDR¶

What this covers¶

ChIP-seq tells you, genome-wide, where a transcription factor is likely bound — but only indirectly, through a cloud of short sequencing reads scattered around the true binding site. This chapter covers how to turn that cloud into a trustworthy list of binding events: a probabilistic model of where a single event's reads land, the EM algorithm used to fit a sparse mixture of such events across an entire genome (GPS), a significance test that uses a matched control experiment to reject noise, and the Irreproducible Discovery Rate (IDR), which decides how many events to keep by checking consistency across replicates rather than any absolute score. It assumes the previous two lectures — aligning high-throughput reads to a genome, and assembling a genome from reads — but builds the probability model, the EM derivation and the hypothesis tests from scratch.

What a ChIP-seq experiment measures¶

Transcription factors are proteins — roughly 2,000 of them in humans — that bind specific DNA sequences and act as molecular switches: a DNA-binding domain recognizes the sequence, and an activation or repression domain recruits (or blocks) RNA polymerase, sometimes by folding the genome to bring distant regulatory sequences together. Finding out which genes a given regulator controls, and precisely where on the genome it sits while doing so, is one of the central techniques for understanding gene regulation. Chromatin immunoprecipitation followed by sequencing (ChIP-seq) is the standard way to do this, one protein at a time.

The protocol:

Crosslink proteins to DNA in living cells — a chemical "flash photograph" that freezes every protein, not just the one of interest, onto the genome.
Harvest and fragment the chromatin, typically by sonication, which breaks the DNA at essentially random positions (fragments of roughly 200–300 bp).
Immunoprecipitate with an antibody against the protein of interest, or against an epitope tag engineered onto it, to enrich for the fragments that protein was sitting on.
Reverse the crosslinks, amplify, and sequence — typically short reads (around 35 bp) from one end of each fragment.
In parallel, sequence the whole-cell extract (WCE) — the same chromatin, fragmented the same way, but without immunoprecipitation — as a control.

A sequencer always reads 5′ to 3′, so what is actually recorded for each fragment is the position of its 5′ end, on whichever strand it came from — not full fragment coverage, and not a count of how many times a base is sequenced. A read is counted once, at the position where it starts. Reads generated by a single binding event therefore form two offset distributions: forward-strand reads pile up just upstream of the site, reverse-strand reads just downstream, with a dip at the site itself, where the crosslinked protein blocked fragmentation and sequencing. With on the order of 10 million cells and a good antibody, this resolves binding locations to roughly 10 bp genome-wide — a spatial resolution that was not achievable even a decade before this lecture.

Why peak-finding alone is not enough¶

The naive approach is peak-finding: smooth the read density along the genome and report each local maximum as a binding event. This works when binding events are far apart, but as two events move closer together their read distributions overlap, and the two peaks merge into a single broad peak. The reads observed at such a locus are, to good approximation, the sum of the two individual per-event distributions — this additivity holds because the experiment pools millions of cells, not all of which have both sites occupied simultaneously, so the measured signal is an average over many independent realizations. Recovering two (or more) component events from their sum — so-called homotypic binding events, closely spaced copies of the same motif — is exactly the problem a peak-finder cannot solve and a generative model can.

Two binding events close together each produce their own read distribution (dashed); the experiment only shows their sum (shaded) — a single broad peak that peak-calling alone cannot split back into two events.

The spatial distribution of reads around one event¶

GPS (the algorithm this lecture builds up — "genome positioning system") starts by characterizing the shape a single binding event produces. For a read at position \(r_i\) caused by a binding event at position \(b_j\):

\[p(r_i \mid b_j) = p(r_i \mid z_{ij}=1) = \text{emp}\big((-1)^{S_i}(r_i - b_j)\big)\]

where \(S_i = 0\) for a forward-strand read and \(S_i=1\) for a reverse-strand read, and \(\text{emp}(d)\) is the empirical spatial distribution of reads at offset \(d\) from a binding site — a histogram built directly from data, not a fitted analytic density, over a window of a few hundred base pairs. For simplicity the model assumes the forward- and reverse-strand distributions are mirror images of each other, so only one needs to be learned and then flipped for the other strand. Because this distribution varies from protein to protein, the plan is to learn it from the data itself and re-learn it as better event calls become available: run the algorithm once with a generic distribution, use the resulting confident events to build a better empirical distribution, and re-run.

GPS: a genome-wide mixture model at single-base resolution¶

The key move is to model every base in the genome as a candidate center of a binding event, and let the data decide which candidates actually have mass. Let \(\pi\) be a vector with one entry per genome position \(m = 1, \dots, M\), where \(\pi_m\) is the probability (the "strength") of a binding event at position \(m\), subject to \(\sum_m \pi_m = 1\). If there are really only two true binding events, the ideal solution has \(\pi_{m_1} = \pi_{m_2} = 0.5\) and zero everywhere else — the constraint doesn't say only one event happened, it says there is exactly one unit of probability mass to be distributed over all the positions that, together, explain the reads.

Given \(\pi\), the probability of a single observed read \(r_n\) is a mixture over every possible originating position:

\[p(r_n \mid \pi) = \sum_{m=1}^M \pi_m\, p(r_n \mid m)\]

and the likelihood of the whole read set \(R = \{r_1, \dots, r_N\}\), assuming reads are generated independently, is the product over reads:

\[p(R \mid \pi) = \prod_{n=1}^N \sum_{m=1}^M \pi_m\, p(r_n \mid m), \qquad \sum_{m=1}^M \pi_m = 1\]

The number of reads \(N\) is typically \(10\)–\(20\) million; the number of positions \(M\) is the size of the genome (billions); the true number of binding events is typically \(5\)–\(30{,}000\) — many orders of magnitude smaller than \(M\), which is exactly why a method that favors a sparse \(\pi\) is the right kind of model, well before any prior is written down. Fitting \(\pi\) by maximum likelihood, \(\hat\pi = \arg\max_\pi p(R \mid \pi)\), would in principle solve the whole problem — but the sum inside the product makes direct maximization intractable, which is where EM comes in.

The EM derivation: if we knew the assignment, this would be trivial¶

Introduce a latent indicator \(g(z_n = m)\): 1 if read \(n\) was generated by the event at \(m\), 0 otherwise. If \(g\) were known, counting binding strength would be immediate: the number of reads attributable to the event at \(m\) is just

\[N_m = \sum_{n=1}^N g(z_n = m),\]

and the maximum-likelihood mixing weight is simply the fraction of all reads assigned to \(m\), \(\hat\pi_m = N_m / \sum_{m'} N_{m'}\). The entire problem is choosing where each read "belongs" — and that choice is latent. EM exploits the fact that a good guess for \(g\) gives a good estimate of \(\pi\), and a good estimate of \(\pi\) in turn gives a better guess for \(g\), and alternates between the two:

E step — given the current \(\pi\), compute \(\gamma(z_n=m)\), the soft, fractional assignment of read \(n\) to event \(m\):

\[\gamma(z_n = m) = \frac{\pi_m\, p(r_n \mid m)}{\sum_{m'=1}^M \pi_{m'}\, p(r_n \mid m')}\]

M step — given \(\gamma\), recompute the effective read counts and update \(\pi\) exactly as if \(\gamma\) were the true (soft) assignment:

\[N_m = \sum_{n=1}^N \gamma(z_n = m), \qquad \hat\pi_m^{(i)} = \frac{N_m}{\sum_{m'=1}^M N_{m'}}\]

\(\pi\) is initialized uniformly, \(\pi_m = 1/M\), and the two steps alternate until convergence; the final \(N_m\) — the effective number of reads assigned to position \(m\) — is reported as the strength of the binding event there. (Because each event's strength depends only on reads near it, the genome can be chopped into independent pieces and processed in parallel.) Like any EM run, the result is a local optimum that depends on the starting point, not a guaranteed global one.

Without a prior, the result is a smear¶

Run as just described, on synthetic data with two planted events 50 bp apart, the algorithm distributes \(\pi\)'s probability mass across a wide, fuzzy range of neighboring positions rather than collapsing it onto the two true locations — and the same happens on real Oct4 ChIP-seq data, where mass leaks out even into positions with only background noise. The model is not wrong, but nothing in it prefers a solution where most positions get exactly zero.

A sparsity prior forces component elimination¶

The fix is to add a prior that actively penalizes nonzero \(\pi_m\) — a negative Dirichlet-type prior:

\[p(\pi) \propto \prod_{m=1}^M \frac{1}{(\pi_m)^\alpha}, \qquad \alpha > 0\]

which grows rapidly as \(\pi_m \to 0\), so the posterior favors driving most components exactly to zero. Maximizing the posterior (now a MAP rather than ML estimate) leaves the E step unchanged, but changes the M step to:

\[\hat\pi_m^{(i)} = \frac{\max(0,\, N_m - \alpha)}{\sum_{m'=1}^M \max(0,\, N_{m'} - \alpha)}\]

Subtracting \(\alpha\) from every component's read count means: if a candidate event doesn't have at least \(\alpha\) reads of support, its count is clipped to zero and it is permanently eliminated — component elimination. The algorithm doesn't eliminate everything below threshold immediately; weak components are allowed to accumulate support over a few iterations before being zeroed, so that real but initially under-supported events aren't killed prematurely.

With this prior, the same synthetic two-event example converges cleanly onto exactly the two planted positions, with the rest of \(\pi\) at zero. On the real Oct4 data, the smear collapses into a handful of sharp components around the true peaks. And on closely-spaced CTCF homotypic events — where a single broad pile of reads conceals two adjacent binding sites — the sparse-prior EM recovers both locations separately, where methods that assume "one pile of reads = one binding event" cannot.

Is an event real? A binomial test against the control channel¶

Not every pile of reads is a real binding event. Genome assembly (an earlier lecture's subject) cannot always resolve the true copy number of repeat elements; when repeats are collapsed, reads originating from many genomic copies all map to the one assembled copy, producing a "tower" of reads that looks exactly like strong, narrow enrichment but is purely an assembly artifact.

This is why the whole-cell extract control exists. The approach is a two-step process:

Run the GPS discovery pipeline on the IP channel alone, ignoring the control, to get a candidate list of events and their strengths \(N_m\).
Filter that list using the control channel: at each candidate event, ask whether the IP read count is significantly larger than the control read count would allow under a null model of pure chance.

The null model treats reads at a given locus, pooled across the IP and control channels, as the outcome of a fair coin flip: each read is equally likely (under the null of no real binding) to have landed in either channel. This equal-likelihood assumption is only valid once the channels' overall read budgets are matched — since the IP channel spends many of its reads on real enrichment while the control does not, the control counts are first rescaled using regions of the genome known to be free of events, matching IP and control read density there before comparing at a candidate event. Given the (scaled) control count \(k\) and total pooled count \(n\) at a candidate locus, the one-sided binomial tail probability is:

\[F(k, n, P) = \sum_{l=0}^{\lfloor k \rfloor} \binom{n}{l} P^l (1-P)^{n-l}, \qquad P = 0.5\]

— the probability, under the null, of seeing as few as \(k\) (scaled) control reads out of \(n\) total. A small \(F\) rejects the null and supports calling the locus a real event.

Multiple testing: Benjamini–Hochberg at a desired FDR¶

A genome-wide scan produces p-values for tens or hundreds of thousands of candidate events simultaneously, so a fixed p-value cutoff is not principled on its own — some fraction of candidates will pass by chance no matter how small the cutoff, purely because so many are tested. The Benjamini–Hochberg procedure controls this directly: sort all candidate events by p-value, smallest to largest (rank \(1\) = most significant), and assign each a \(Q\)-value,

\[Q\text{-value} = P\text{-value} \times \frac{\text{Count}}{\text{Rank}}\]

where \(\text{Count}\) is the total number of events tested. Accept events at ranks \(1, \dots, k\) up to the largest \(k\) for which the \(Q\)-value is still below the desired false discovery rate (FDR) — commonly \(0.05\). The FDR has a direct interpretation: if \(1{,}000\) events are called at FDR \(0.1\), roughly \(100\) of them are expected to be false positives — it is the fraction of the accepted list expected to be wrong, not a per-event significance threshold.

Are two replicates consistent? The Irreproducible Discovery Rate¶

Every real ChIP-seq experiment should be done at least twice. Given two ranked lists of events from two replicate experiments — ranked from most to least significant, with matching events identified by shared genomic location — the question is how far down the two lists agreement can be trusted.

Spearman rank correlation: necessary, not sufficient¶

A natural first idea is to correlate the two experiments' rankings. The Spearman rank correlation is the ordinary (Pearson) correlation computed on ranks rather than raw scores:

\[\rho = \frac{\sum_i (x_i - \bar x)(y_i - \bar y)}{\sqrt{\sum_i (x_i - \bar x)^2 \sum_i (y_i - \bar y)^2}}\]

where \(x_i, y_i\) are the ranks of matched event \(i\) in lists \(X\) and \(Y\). This is useful precisely because it does not depend on the specific numeric scores each experiment happens to produce — only on their order — so two experiments whose raw scores correlate only loosely (say, a Pearson correlation of \(0.88\)) can still be in perfect rank agreement (\(\rho = 1\)). But a single correlation coefficient summarizes the whole list at once; it cannot say where, within the list, agreement breaks down — and that breakdown point is exactly what is needed to decide how many events to keep.

The correspondence curve¶

Define \(\Psi_n(t)\) as the fraction of the top \(n \cdot t\) events (by rank) that appear in both lists' top \(n \cdot t\) — i.e., how well the two replicates agree once you restrict to the top fraction \(t\) of the ranking. In the idealized case where the top half of all \(n\) events are genuine signal, perfectly concordant between replicates, and the bottom half is pure noise with no concordance at all, \(\Psi_n(t)\) rises with roughly unit slope up to \(t = 0.5\) (every additional event admitted is a true positive, matched in both replicates) and then flattens — additional entries past that point are noise, matched only by chance. The derivative \(\Psi_n'(t)\) makes the transition sharper to see: high while still in the reproducible region, dropping once the ranking crosses into noise.

In the idealized case, half the ranked events are true signal and perfectly concordant between replicates, half are noise with no concordance — the curve's break point is what IDR estimates from real, non-idealized data.

Fitting the break point: a two-component mixture¶

Real data isn't a clean two-piece idealization, so IDR models the paired ranks directly as a mixture of two components — a "reproducible" component and an "irreproducible" component — and fits the mixture's parameters from the observed pairs. Given a desired irreproducible discovery rate \(\alpha\), IDR selects the top \(\ell\) pairs such that the (fitted) probability that a selected pair actually belongs to the irreproducible component is \(\alpha\). This is how the ENCODE pipeline uses it in practice: each replicate is processed independently through event discovery, both resulting event lists are ranked by significance, and IDR is applied to the pair of rankings to decide how far down each list to trust. Different peak-calling methods reach IDR \(= 0.05\) after retaining very different numbers of peaks — some run out of reproducible signal quickly, others can go much further down their own ranked list before hitting the bound.

Consistency is not correctness¶

A caution worth holding onto: an algorithm that is consistently wrong — for example, one that always calls the same assembly-repeat "towers" as top events, in the same order, every time it is run — will score perfectly by IDR, because its answers are reproducible, not because they are correct. Reproducibility is necessary but not sufficient for validity. A useful sanity check is to run IDR on a pairing that should not agree — an experiment against a mismatched or unrelated dataset — and confirm that it correctly fails to find reproducible structure; a method that reports high agreement even there is not to be trusted on the real comparison either.

Pushing resolution further: motifs as a positional prior (GEM)¶

The lecture closes by pointing toward an extension built on exactly the same machinery: once you can resolve binding to within roughly 20 bp genome-wide, you can ask whether candidate binding events sit near known or discovered sequence motifs, and use that as evidence. GEM couples binding-event finding to motif discovery in both directions — biasing discovered motifs toward predicted binding sites, and biasing binding-event predictions toward motif positions — by replacing the uniform sparsity prior with a position-specific prior:

\[p(\pi) \propto \prod_{m=1}^M (\pi_m)^{-\alpha_s + \alpha_m}\]

where \(\alpha_s > 0\) is the same uniform sparseness parameter as before, and \(\alpha_m\) adds extra weight at positions where a motif has been found, making the model more tolerant of a weak read-count event there. The payoff the lecture gestures at — resolving binding events to the point of computing, genome-wide, the significant spacings between pairs of profiled transcription factors — is this same analysis: once events are localized this precisely, you can start reading off the regulatory grammar the genome actually uses (for instance, a substantial fraction of Oct4 sites sitting within 100 bp of a Sox2 site, suggesting cooperative or co-occurring binding).

Sources¶

Slides: lectures/07-slides/01-lecture-7.md (ChIP-seq background, protocol, and the per-event spatial distribution \(p(r_i \mid b_j)\)) and 02-gps-probabilistically-models-chip-seq-read-spatial-distribut.md (the GPS mixture model and EM derivation, the sparse prior and component elimination, the binomial significance test, Benjamini–Hochberg, Spearman correlation, the IDR correspondence-curve slides, and the GEM positional prior) — both from MIT OCW 7.91J Spring 2014, Lecture 7, CC BY-NC-SA 4.0. These slides are a model's reconstruction of a PDF with no text layer; treat the equations as a pointer into the original, not as verified transcriptions.
Transcript: recordings/lectures/07.md (captions, timestamped) — the protocol walkthrough and its rationale (00:00–14:40), the spatial-distribution and sparsity motivation and class discussion (24:00–50:00), the EM derivation given live on the board (41:00–48:00), the significance test and control-channel normalization discussion (59:00–1:08:00), Benjamini–Hochberg and FDR (1:03:00–1:05:00), and the IDR discussion including the Spearman example, the correspondence curve, and the "consistency is not correctness" caution (1:08:00–end).
The slide deck's closing section on GEM (positions ~210–290 of the second slide file) is only lightly gestured at in the recording — the professor explicitly runs out of time ("I will not have time to present all the rest of this") right after describing the spatial-constraint payoff this section ends with — so it is covered here in outline only, from the slides.
Two further slide sections — "Concept of a Transcriptional Regulatory Code" (Harbison et al., Nature 431, 2004) and the cross-species conservation of binding events (Odom et al., Nature Genetics 39, 2007) — appear at the end of the deck but are not reached or discussed in the recorded lecture at all, and are not covered in this chapter.
Named but not supplied here: the GPS paper itself (referenced in the transcript as "the GPS paper, which I posted," for the exact setting of \(\alpha\)); Figueiredo and Jain (2002), the source of the sparse Dirichlet-type prior; Kharchenko, Tolstorukov and Park (2008), on the shape of ChIP-seq read distributions; Wang and Zhang (2011), SeqSite; Li, Brown et al. (2011), the IDR paper; and Guo, Mahony et al. (2012), the GEM paper.

← 6. OLC and De Bruijn Assembly · Contents · 8. RNA-seq, Isoforms, and Expression Statistics →