ChIP-seq maps where a protein sits on DNA. This page simulates a wildtype and an R175H p53 experiment on the real CDKN1A locus, calls peaks in both, and then checks the p53 motif in the real sequence against real p53 experiments.
Five Views showed several routes by which the TP53 pathway can fail. ChIP-seq (chromatin immunoprecipitation followed by sequencing) gives direct evidence for one of them: where a protein was bound on the genome when the cells were fixed. Does the R175H mutant of p53 still bind the CDKN1A promoter? The reads below are simulated, so at CDKN1A the answer is built in; what the page can show is how the evidence is read, where a caller draws its line, and what the real sequence and real experiments say about the sites.
The assay cross-links proteins to DNA, shears the DNA into fragments of a few hundred bases, pulls down the fragments bound by the target protein with an antibody, and sequences them. Mapping the reads back to the reference gives a pileup: how many fragments cover each position. Because the assay sees whole fragments, a bound site shows up as a hill about a fragment wide, not as the protein's 20-base footprint.
Any antibody pull-down also collects background: open chromatin breaks more easily during shearing, amplified regions have more copies, repeats confuse the mapping, and some sequences are simply sticky. The input control, sheared DNA from the same cells with no pull-down, measures that local background. A good transcription-factor peak rises above input, appears in biological replicates, and usually has the expected motif near its summit.
Figure 1. Simulated reads on the real locus: fragment pileup (fragments covering each 20-bp bin) over 14 kb of chr6 around CDKN1A, positions in kb from the TSS of its Ensembl canonical transcript ENST00000244741 (GRCh38). The gene models are real; the reads are simulated. Wildtype p53 is given the two classic response elements (dashed), the 5′ site at and the 3′ site at , plus a weaker site in intron 1 at , where real experiments also find one (Figure 5). R175H keeps 10% of wildtype occupancy at all three. All three samples share a 2-fold bump of open chromatin at the promoter. The triangle is the caller's summit (Figure 2). Hover, tap, or focus the tracks and use the arrow keys for the counts at one position.
The wildtype track rises to fragments at the 5′ site against about in the flanks, and the hill is bp wide at half height, about one fragment long, since the site is seen through the fragments that cover it. The two response elements are bp apart and read as two hills. They sit kb upstream of the TSS, and the coding sequence starts another kb downstream, in exon 2.
Input is flat except at the promoter, where every sample, pulled down or not, has extra fragments from open chromatin. The R175H track shows nothing a reader would call a peak. R175H is a structural hotspot: arginine 175 sits in the loops held by the zinc ion that keep the DNA-binding surface in shape, so the mutation can unfold the surface needed for sequence-specific binding. Contact hotspots such as R248 and R273 instead change residues that touch the DNA.
Motif searches, target lists and differential binding all need a discrete list of bound regions. Peak calling turns the continuous pileup into one. MACS, a widely used caller, compares the ChIP pileup with a background rate estimated locally from the input. The caller on this page follows that recipe in miniature. For each 20-bp bin it counts the ChIP fragments covering it, \(k\), and asks how surprising \(k\) would be from background alone:
\(p = P(X \ge k), \;\; X \sim \text{Poisson}(\lambda_{\text{local}})\)
\(\lambda_{\text{local}} = \max(\lambda_{\text{BG}}, \lambda_{1\,\text{kb}}, \lambda_{5\,\text{kb}}, \lambda_{10\,\text{kb}})\)
The rate is the largest of the region-wide average and the input's rate in 1, 5 and 10 kb around the bin, so a bin in open chromatin is compared with open chromatin. The tail is summed exactly, in logarithms, so even a summit with p near \(10^{-60}\) gets its own value. The p-values of all bins in the 40-kb simulated region go through Benjamini-Hochberg to give q-values, and runs of bins above the threshold, joined across gaps of up to 40 bp, become peak regions. The summit is the bin with the most fragments.
Figure 2. A toy caller in the spirit of MACS, on the simulated samples of Figure 1. Bars are −log10 q for each 20-bp bin (square-root axis). Bins at or above the dashed threshold are called and joined into regions (shaded), each with a summit dot. The default, q ≤ 0.05, is MACS2's default cutoff; the slider is logarithmic.
In practice the decision is not made at the caller's threshold. ENCODE's pipelines run the caller loose on purpose, so the weak end of the list still holds noise, then keep the peaks that independently prepared replicates rank the same way, a cut made on the irreproducible discovery rate (IDR). The same guidelines set the scale of a usable experiment: at least 10 million uniquely mapped reads per replicate, and a fraction of reads falling in peaks (FRiP) of 1% or more, which 787 of the 1,052 ENCODE data sets they reviewed reached.
Has R175H lost binding everywhere, or only moved it? Figure 3 asks the first half, running the same simulation and caller at ten canonical p53 targets, two replicates per sample, and calling each site from the result rather than from a label.
Figure 3. Simulated differential binding. Each row is its own simulated 40-kb region run through the Figure 2 caller, with two replicates per sample. Dots are fragments within 250 bp of the site in each replicate (log axis); a filled dot means that replicate's peak was called there at q ≤ 0.05, and the gray tick is the input. The line joins the wildtype mean to the R175H mean, blue for less and red for more. Verdicts: lost, called in both wildtype replicates and neither R175H; gained, the reverse; reduced or kept, called in both samples, with log2 fold change of fragments above input at most −1 or above it; unclear, the replicates disagree. The gene names are real HGNC symbols; the occupancies are chosen for the simulation. The lower three rows are promoters bound by gain-of-function mutants in Zhu et al. (2015); the simulation gives R175H binding there and wildtype none. Focus or tap a row for its numbers.
A table of canonical sites cannot see binding that moved, because it only looks where wildtype p53 was; that needs peaks called in the mutant sample on its own. Zhu et al. (2015) did this for the gain-of-function mutants R273H, R248Q and R249S, not R175H. The mutants bound the promoters of the chromatin regulators KMT2A (MLL1), KMT2D (MLL2) and KAT6A (MOZ), which wildtype p53 did not bind, while wildtype p53 bound CDKN1A, MDM2 and BBC3 as expected. The mutant peaks were enriched for ETS motifs rather than the p53 motif, and the mutant protein binds the transcription factor ETS2, which the authors propose brings it there. The lower rows of Figure 3 borrow that result, so in this simulation the answer is both, because both were put in; in the published experiments, binding moved to sites the p53 motif does not mark.
If p53 finds its sites by reading DNA, why not find them from the sequence? The p53 response element is two copies of the ten-letter pattern RRRCWWGYYY (R is A or G, W is A or T, Y is C or T), separated by 0 to 13 bases. Both CDKN1A elements are real instances, with no spacer, and each matches of the 20 letters. A consensus counts every letter alike. JASPAR's TP53 matrix (MA0106.3) instead counts each base at 18 positions in about sequences that p53 picked out in vitro, and scores a window by adding up the log-odds of each of its bases in those sequences against background.
Figure 4. Top: JASPAR MA0106.3 as a sequence logo. A stack's height is the position's information content, 2 − H bits against equal base frequencies, from the matrix's raw counts ( sequences per position); letters are stacked by frequency, the commonest on top. With that many sequences the small-sample correction is at most bits, so it is not applied. Below: the 20-letter consensus and the six best-scoring windows of Figure 5's scan, down to the 3′ element, each read on the strand that matches the matrix and extended one base at each end. Lower case misses the consensus. Shading is what a base costs against the best base at that position under the scanning matrix (log-odds with pseudocount 0.8 against the region's GC content); a row's costs add up to its distance from the best possible score, bits. After each row: its score in bits (on wide screens) and how many of the ReMap TP53 data sets have a peak on it, bold above 5% of them. All real data. Hover, tap or focus a row (arrow keys step rows and positions); the same window is ringed in Figure 5.
Tall stacks are where p53 is strict: the C and G that open and close each CATG core carry bits, the outer letters of each half . So two missed letters can cost little or a lot. The 5′ element's misses fall on loose positions and cost it bits in all. The 3′ element has a C where its second core wants G, which costs bits on its own. Figure 5 scans 20 kb of the real locus with the matrix and checks every strong match against ReMap 2022, which holds p53 ChIP-seq data sets with a peak in this region, mostly from cells with wildtype p53, many treated with DNA damage or the drug nutlin to activate it.
| Rank | Window | Score (bits) | ReMap data sets |
|---|
Figure 5. Real sequence, real binding: chr6 within 10 kb of the CDKN1A TSS (GRCh38). Top: the JASPAR MA0106.3 score of every 18-bp window reaching 75% of the matrix's range, the better of its two strands, against a background of the window's own GC content () with pseudocount 0.8; matches at or above the threshold are dark and ranked. Below it: how many ReMap 2022 TP53 data sets put a peak summit in each 100-bp bin. Table: the ranked matches and how many data sets have a peak over each (of ). Histogram: the number of matches at the threshold in the real sequence (line) against 100 dinucleotide shuffles of it, which keep its length and every two-letter count (Altschul-Erickson shuffle, seeds 1 to 100). Hover, tap or focus a table row to ring its match in the plot and, for the top six, its row in Figure 4.
ChIP-seq says where a protein was bound at fixation, and nothing about whether that binding changed transcription. A bound site counts as regulation only when the nearby gene's expression moves with it, which is why ChIP-seq is paired with RNA-seq from the same cells. The next article takes that second measurement, where CDKN1A turns up among the genes that truly change.