9 Differential binding and accessibility
Peaks tell you where a signal is. A differential analysis tells you where the signal changes between conditions. This chapter explains the five steps of a differential analysis, the choices at each step, and how to read the result. The example is the H3K27ac dataset GSE155100 (Chapter 14).
9.1 The question
Which regions have more, or less, signal in one condition than in another? And how sure can we be, with only a few replicates per condition?
A peak that is present in one condition and absent in the other is not enough. Peak calling depends on depth and on thresholds (Chapter 7). A differential analysis instead counts the reads in the same regions in every sample, and asks whether the difference between conditions is larger than the differences between replicates.
9.2 How it works
A differential analysis has five steps:
- Choose the regions to test.
- Count the reads of each sample in each region.
- Normalise: correct each sample for how deeply it was sequenced.
- Test each region: is the difference between conditions larger than expected from the replicates?
- Call the regions that pass the thresholds, and estimate how large each change is.
1. Which regions to test
The regions must be the same for every sample. Three kinds are common:
| Region type | Used for | Why |
|---|---|---|
| Peaks called on all libraries pooled, merged | Sharp histone marks such as H3K27ac and H3K4me3 | These marks cover a few kilobases. Whole peaks hold more reads than short windows. |
| Fixed windows around peak summits (501 bp) | Transcription factors and ATAC-seq | The sites are short. “Most DNA regulatory elements are less than 500-bp” (Grandi et al. 2022). |
| Sliding windows, merged (csaw) | Broad marks such as H3K27me3 | Broad domains have no clear summits (Lun and Smyth 2016). |
Lun and Smyth compared ways to choose regions for a differential test. They found that “peak-based methods are more powerful for detecting diffuse DB regions”, and that “peaks should be called from pooled libraries” (Lun and Smyth 2014). Peaks called on the pooled libraries do not depend on the condition of any one sample. This keeps the error rate of the later test close to the level you asked for.
2. Counting reads
For each sample and region, the software counts the reads that fall in the region. The result is a table with one row per region and one column per sample. Inputs are not counted: they only helped to call the peaks.
3. Normalisation
Samples are sequenced to different depths. Before comparing them, each sample’s counts are divided by a size factor, a number for how deeply that sample was sequenced. There are several ways to compute it:
- RLE (median-of-ratios), DESeq2’s default: for each region, the ratio of the sample’s count to a typical sample’s count, then the median of these ratios (Love et al. 2014).
- TMM (trimmed mean of M-values): the average log ratio between samples, after leaving out the regions with the largest changes and the most extreme counts (Robinson and Oshlack 2010).
- Library size (reads in peaks): each sample’s total number of reads in the tested regions.
- Spike-in: the reads from a fixed amount of DNA from another species, added to each sample (Chapter 10).
The first three methods assume that most regions do not change. TMM, for example, “assumes that most regions are not truly DA [differentially accessible]” (Reske et al. 2020). If many regions change in one direction, this assumption breaks. Chapter 16 shows an example, and Chapter 10 shows how spike-ins help.
4. The test
Read counts vary between replicates more than chance alone would explain. DESeq2 models this extra variation with a dispersion value for each region (Love et al. 2014). With two or three replicates, one region’s own dispersion is a poor estimate. So DESeq2 borrows from the other regions: “we assume that genes of similar average expression strength have similar dispersion” (Love et al. 2014). It pulls each region’s dispersion towards the trend of regions with similar counts.
Then DESeq2 tests each region (a Wald test) and gives it a p-value. With tens of thousands of regions, some p-values are small by chance. So the p-values are corrected into a false discovery rate (FDR) with the Benjamini–Hochberg method (Love et al. 2014). An FDR threshold of 0.1 means that up to about 10 % of the called regions are expected to be false.
Before the correction, DESeq2 sets aside regions with very few reads, which have little chance of reaching significance. This “independent filtering” makes the correction less strict for the other regions (Love et al. 2014). Regions set aside get no FDR.
5. Calling regions and measuring the change
A region is called gained or lost when it passes an FDR threshold and, often, a minimum fold change. The fold change is reported as log2: 1 means twice as many reads, −1 half as many.
Fold changes of regions with few reads are noisy: a count of 3 against 1 is a three-fold change. Shrinkage pulls such noisy fold changes towards zero, and leaves well-measured ones almost unchanged. EpiNexus uses the apeglm method (Zhu et al. 2019). Shrunk fold changes are better for ranking and plotting regions (Love et al. 2014).
Report the FDR and any fold-change cut-off with your results. DESeq2’s own default FDR is 0.1: “By default the argument alpha is set to 0.1”. To require a minimum change, DESeq2 can test against a fold-change threshold, instead of filtering the list afterwards (Love et al. 2026). No publication sets these values for ChIP-seq regions. EpiNexus calls a region significant at FDR < 0.1 and an unshrunk |log2 fold change| above 0.5; the second rule keeps tiny changes out of the list. Shrunk fold changes are used for ranking and for the MA plot.
9.3 In EpiNexus
On the Parameters step of the run wizard:
- Differential comparisons lists every pair of conditions for each mark, ticked. Untick the pairs you do not need. A label such as control, WT or untreated is used as the reference. Use the swap button to change the direction.
- Normalization: RLE (DESeq2 default), TMM or Library size (reads in peaks). The defaults are TMM for ChIP-seq, library size for ATAC-seq, and RLE with an E. coli spike-in for CUT&Tag and CUT&RUN.
- Samples with the same replicate number were processed together (paired design): tick this only if each replicate pair really was processed together, for example on the same day. The model then removes the differences between pairs.
EpiNexus chooses the regions from the mark (Table 9.1) and tests them with DESeq2 in R. Broad marks are tested with csaw (Lun and Smyth 2016). The Differential page shows the volcano plot, the MA plot, the sample PCA and a heatmap of the top regions. It states the regions, the normalisation and the thresholds that the run used. Download all N peaks (CSV) gives every tested region with its fold changes, p-value and FDR.
Figure 9.1 shows the top of the Differential page for our example.
9.4 Reading the result: GSE155100
The regions and the counts
For H3K27ac, EpiNexus called peaks on the four ChIP libraries pooled, against the four inputs pooled. After merging, and splitting the regions longer than 5 kb, this gave 92,052 regions. Their median width is 710 bp; 17 % are longer than 2 kb.
TMM gave size factors between 0.89 and 1.12. So the four libraries had similar depths, and no sample needed a large correction.
What the thresholds did
| Step | Regions |
|---|---|
| Tested | 92,052 |
| Set aside by independent filtering (too few reads) | 39,263 |
| FDR < 0.1 | 2,485 |
| FDR < 0.1 and |log2 fold change| > 0.5 (called) | 2,461 (2,267 lost, 194 gained) |
| For comparison: FDR < 0.05, no fold-change rule | 1,963 |
Two things stand out (Table 9.2):
- 43 % of the regions got no FDR. DESeq2’s independent filtering set them aside before the correction, because with so few reads they had little chance of passing it. They still have a p-value.
- The FDR threshold matters more than the fold-change rule. The fold-change rule removed only 24 regions. A stricter FDR of 0.05 would have removed about 500.
The called regions also hold more reads than average: their median mean count is 82, against 17 for all regions. More reads give more power to see a change.
The MA plot
The MA plot (Figure 9.2) shows each region’s mean count (A, across) against its fold change (M, up). It answers two questions.
- Is the normalisation sound? The grey cloud of unchanged regions should be centred on zero, at all counts. Here it is. If the cloud bent up or down, the normalisation would be suspect.
- Where are the changes? The lost regions span the whole range of counts. So the losses are not limited to weak regions.
The plot also shows shrinkage at work. At the left, regions with few reads sit close to zero. In this run, 777 regions have an unshrunk |log2 fold change| above 2 but a shrunk one below 0.5. Most have few reads: their median mean count is 6.6, against 17 for all regions. One has a mean count of 3.2: its unshrunk log2 fold change is 2.5, but its shrunk value is 0.03 and its p-value is 0.14. Without shrinkage, such regions would top a list ranked by fold change.
The sample PCA
The sample PCA (Figure 9.3) places the samples by their counts. Samples that are similar lie close together.
- The first component (63.9 %) separates the two conditions. The genotype is the largest source of difference.
- The second component (20.9 %) separates the replicates. In each condition, one sample is at the top and one at the bottom. In EpiNexus’s normalised counts, the two first replicates lie on one side and the two second replicates on the other. So the replicate pairs differ in a shared way. This can happen when each pair was prepared or sequenced together.
If you know that the replicate pairs were processed together, tick the paired design in the run wizard (see the box above). If you do not know, leave it unticked: with two replicates per condition, the paired model leaves very little to estimate the noise.
9.5 Pitfalls
Regions with few reads can show huge fold changes by chance. Rank by FDR, or by shrunk fold change.
A cloud of unchanged regions that is not centred on zero means that the normalisation does not fit the data. Look at the MA plot before you read the list of regions.
An FDR of 0.1 says that up to about one in ten regions in the list is expected to be false. It does not tell you which ones. Check the strongest regions in the genome browser (Chapter 8).
With two replicates, outliers cannot be detected, and the noise estimate leans on the trend across regions. EpiNexus says so on the Differential page. The strongest changes are the most likely to be real; changes near the threshold are less certain (Chapter 2).
9.6 Further reading
- Love et al. (2014): DESeq2, its dispersion model, independent filtering and shrinkage.
- Lun and Smyth (2014) and Lun and Smyth (2016): choosing regions for differential binding, and csaw.
- Robinson and Oshlack (2010): TMM normalisation.
- Reske et al. (2020): how normalisation changes differential ATAC-seq results.