20 Epigenomics and Chromatin: Concepts

20.1 The problem

Every cell in an organism carries essentially the same genome, and they are not the same cell. What differs is which parts of that genome are reachable, which proteins are sitting on them, and which chemical marks have been placed on the DNA itself. Epigenomics is the attempt to measure those differences with the same sequencing machine that reads the sequence.

The measurements all work by the same trick: turn a physical property of chromatin into a difference in how many reads land somewhere. An accessible region yields more fragments than a closed one. A region bound by a protein yields more fragments than an unbound one. That conversion is what makes the question answerable, and it is also the source of every difficulty in the rest of this chapter, because read density is affected by a great many things that are not the property you meant to measure.

Nothing here measures chromatin directly. Every assay measures read density as a proxy, and the analysis is the argument that a local excess of reads is attributable to the biology rather than to copy number, mappability, chromatin being generally more soluble where it is open, or the library having been amplified from very little material. The tools cannot make that argument. They can only make it checkable.

20.2 Three assays, one file format

The assays that dominate this field produce BAM files that look alike and answer questions that are not alike at all.

ChIP-seq immunoprecipitates chromatin with an antibody against a protein or a histone modification, and sequences what came down. Its inference problem is enrichment: reads accumulate where the antibody found its target, above a background of chromatin that was pulled down for other reasons. The quality of the antibody is a parameter of the experiment that no amount of computation can recover.

ATAC-seq uses a transposase that inserts sequencing adapters preferentially into accessible DNA, so accessible regions are fragmented and sequenced while closed ones are not (Buenrostro et al. 2013). Its inference problem is also enrichment, but there is no antibody and no target — the readout is a property of the chromatin itself, which makes it easier to run and harder to interpret, because “accessible” is a description of many different underlying situations.

Bisulfite sequencing and its relatives do something categorically different: they change the DNA chemically so that a modified base and an unmodified base sequence as different letters (Frommer et al. 1992). The readout is not a count of fragments at a location but a ratio at a single base, and this is why it does not belong to the same analysis family despite arriving in the same format.

A newer group of targeted methods — cutting or tagging chromatin in place rather than fragmenting and pulling it down (Skene and Henikoff 2017; Kaya-Okur et al. 2019) — sit with the first two in kind but not in statistics. They produce far less background, which sounds like an unambiguous improvement and is mostly one, but it invalidates a modelling assumption the older methods rely on. A method that estimates local background from the data needs there to be background.

Enrichment assays ask where, and their answer is a set of intervals. Modification assays ask how much, here, and their answer is a fraction at each base. The file format hides this distinction and the analysis cannot afford to: almost every mistake in this field comes from applying the reasoning of one to data from the other.

20.3 What a peak is

A peak is not an observation. It is a claim that the number of reads in a window exceeds what a background model predicts for that window, at some threshold. Both halves of that sentence are choices, and the background model is where the method actually lives.

A genome-wide average makes a poor background. Read density varies across the genome for reasons that have nothing to do with the protein you precipitated: copy number differs between regions and between cell lines, mappability differs because repetitive sequence cannot be uniquely placed, chromatin that is open is more accessible to everything including the reagents, and amplification is uneven. A region with twice the average depth for any of these reasons will clear a genome-wide threshold everywhere along its length.

The standard response is to estimate the expected count locally, from a window surrounding the candidate rather than from the genome. This absorbs whatever varies slowly along the chromosome, which is most of the confounders above, and it is the single design decision that separates a usable caller from a coverage threshold.

The second half — the threshold — determines how many peaks you get, and the number of peaks is not a property of the biology. It is a property of the threshold, the sequencing depth and the signal-to-noise ratio of that particular library. Two experiments on the same tissue at different depths will report different peak counts, and neither number means anything on its own.

A peak count is not a result

“We identified 12,483 peaks” is a statement about a threshold applied to a library, not a finding. Sequence the same sample twice as deeply and the number goes up, because weak enrichment that previously fell below the cut now clears it. If a comparison between conditions rests on how many peaks each one had, it is measuring depth and library quality alongside whatever else it found.

20.4 Signal, noise and regions that lie

Before any peak is called there is a prior question: did the experiment enrich for anything at all? A failed immunoprecipitation still produces reads, still produces a BAM, and still produces peaks when a peak caller is pointed at it.

The useful diagnostics all ask some version of “how concentrated is the signal”. If a small fraction of the genome holds a large fraction of the reads, the assay worked; if reads are spread evenly, it did not, whatever the peak caller subsequently reports. The consortia that had to do this at scale wrote their practices down, and those guidelines remain the reference for what a defensible enrichment experiment looks like (Landt et al. 2012).

Two structural traps sit underneath this.

Control libraries. An enrichment experiment without a matched control — input chromatin, or a non-specific pull-down — has no way to distinguish enrichment from the background structure described above. The control is not a formality; it is the term that makes the model identifiable.

Problematic regions. Some parts of every assembled genome accumulate high read counts in essentially every experiment, regardless of the assay or the sample: unresolved repeats, collapsed segmental duplications, regions where the assembly disagrees with the biology. They produce enormous, extremely significant, entirely artifactual peaks. Curated exclusion lists exist for the common assemblies (Amemiya et al. 2019), and filtering against one is a routine step whose omission is visible in the results as a set of implausibly strong peaks that recur across unrelated papers.

The exclusion list is assembly-specific

A list built for one reference does not transfer to another, and applying the wrong one silently removes the wrong intervals. This is a coordinate problem of exactly the kind Section 2.3 describes, with the added hazard that the failure produces a plausible result rather than an error.

20.5 A peak set is not a measurement

The most common analysis in this field compares conditions, and the most common way of doing it is wrong in an instructive way: call peaks separately in each condition, then intersect the sets and treat what does not overlap as differential.

The problem is that a peak set is a thresholded object. A region with genuinely equal signal in two conditions, sitting near the threshold, will appear in one set and not the other about as often as not. The resulting “condition-specific” peaks are dominated by regions where the evidence was weakest, which is the opposite of what the analysis was meant to find. Nothing about the output distinguishes this from a real result, and the effect gets stronger as the two libraries differ in depth.

The correct shape is the one the expression chapters use (Section 16.6): define a set of regions once, count every sample over that same set, and test the counts with a model that knows about replicate variability. What the regions are can be decided by peak calling, by fixed windows, or by an annotation — but they must be the same regions for every sample, and the testing must happen on counts rather than on set membership.

There is a circularity to watch for even in the correct version. If the regions are chosen by calling peaks on the same data you then test, the choice has already used the signal. This is usually tolerable, because the selection is made without reference to the group labels, but it stops being tolerable the moment peaks are called separately per condition and then pooled. Choose regions in a way that is blind to the comparison.

20.6 What a methylation measurement is

Bisulfite treatment deaminates unmethylated cytosine to uracil, which is read as thymine, while methylated cytosine is protected and still reads as cytosine. The measurement at a given base is therefore the proportion of reads carrying C rather than T — a ratio with a denominator, not a count of fragments.

Three consequences follow, and all three are structural rather than implementation details.

The genome loses a letter. After conversion, most cytosines in the reads are thymines, so aligning against the ordinary reference fails. The standard approach converts both reads and reference into a three-letter alphabet, aligns in that reduced space, and restores the original bases afterwards. Reducing the alphabet reduces the information available to place a read, so bisulfite alignment is intrinsically less unique than ordinary alignment, and mapping rates are lower for reasons that are not a bug.

A C-to-T SNP is indistinguishable from an unmethylated cytosine. Both produce a T where the reference has a C. Nothing in the reads separates them, so a study that cares about individual-level methylation needs genotypes, and one that does not will silently report a fixed genetic difference as complete demethylation.

Incomplete conversion inflates methylation. Any cytosine that escapes the chemistry reads as methylated. The conversion rate is therefore a parameter of every experiment, and it is estimated from sequence known to be unmethylated — a spiked-in control, or the mitochondrial genome, or non-CpG cytosines in a tissue where those are not methylated.

Sequencing technologies that detect base modifications directly, from the raw signal rather than from a chemical conversion (Simpson et al. 2017), avoid all three problems at once: the read retains its own sequence, so genotype and methylation come from the same alignment and no conversion efficiency is involved. What they introduce instead is a per-base probability that has to be thresholded, which moves the uncertainty rather than removing it.

20.7 Sites, regions, and what is being tested

Methylation is measured per cytosine and is almost never interesting per cytosine. Neighboring sites are strongly correlated, individual sites are noisy at realistic coverage, and the biology people care about — a promoter, an enhancer, an imprinted locus — is a region.

So the analysis aggregates, and the aggregation is where the statistics are usually lost. The tempting procedure is to test each cytosine, join adjacent significant ones into regions, and report those regions with the best p-value found inside each. That p-value describes a single site chosen for being extreme. It says nothing valid about the region, and correcting it for the number of sites tested does not fix the problem, because the object being reported is not a site.

Getting this right means treating the region as the unit of inference from the start, and obtaining a null distribution for whole regions — typically by permuting the sample labels and seeing what region-level statistics arise by chance. This is more work and it is the difference between a false discovery rate that means what it says and one that does not.

Coverage is part of the measurement

A site where three of four reads are methylated and a site where 300 of 400 are both give a ratio of 0.75, and they are not equally informative. Any method that works from ratios has discarded the denominator, and the sites it is most confident about will be the low-coverage ones, where a ratio can reach 0 or 1 by chance. Models that keep the counts — methylated out of total — do not have this failure mode.

20.8 Reproducibility is the assay’s own question

In most of this book, replicates exist so that a statistical test has an error estimate. Here they do something additional: they are the primary evidence that the assay worked at all.

An enrichment experiment has no ground truth and no internal control for antibody quality or transposition efficiency. What it has is the observation that two independent preparations of the same material should find the same regions, and should rank them similarly. When they do not, the disagreement is about the experiment rather than about the biology, and no downstream analysis can repair it.

This is why the reproducibility of a ranking, rather than the overlap of two thresholded sets, is the quantity worth measuring — it uses the whole signal rather than only the part above a cut, and it degrades gracefully instead of collapsing when one replicate is slightly shallower. The same logic argues against merging replicates before calling: pooling produces a deeper library and hides exactly the disagreement you needed to see.

Two replicates that agree strongly tell you the assay is repeatable. They do not tell you the antibody bound what it was supposed to bind, and a consistently wrong experiment is consistent. Reproducibility is necessary and it is not sufficient, which is why the enrichment diagnostics of Section 20.4 are a separate check rather than the same one.