7 Peak calling
A peak is a region with more reads than you would expect by chance. Peaks are where the antibody found its target, or where chromatin is open. Peak calling turns millions of reads into a list of regions. This chapter explains how the three peak callers in EpiNexus work, how to choose between them, and how to check that peaks are reproducible.
7.1 The question
Where in the genome is the signal real? Reads are never spread evenly. Some regions have more reads because of the target. Others have more reads because DNA there is easier to break, to sequence or to align. A peak caller must tell these apart.
7.2 How it works
Every peak caller compares the reads in a region with a background: how many reads you would expect there without real signal. The background can come from a control sample (input or IgG), from the rest of the genome, or from both. If a region has many more reads than its background, it becomes a peak.
Each peak gets a score. The p-value is the chance of seeing at least this many reads if the region had only background. The q-value corrects the p-value for the millions of regions tested. A q-value of 0.05 means that up to about 5 % of the peaks at this threshold are expected to be false.
7.3 The peak callers
MACS3: ChIP-seq and ATAC-seq
MACS is the most widely used peak caller for ChIP-seq (Zhang et al. 2008). EpiNexus uses its current version, MACS3. MACS does two things that suit ChIP-seq well.
It finds the centre of each site. With single-end reads, you only see one end of each fragment. Reads on the forward strand gather just before a binding site, and reads on the reverse strand just after it. MACS measures the distance between the two groups, called d. Then it “shifts all the tags by d/2 toward the 3′ ends to the most likely protein-DNA interaction sites” (Zhang et al. 2008). With paired-end reads, the fragments are known, so no shift is needed.
It uses a local background. Some regions have more background reads than others. MACS measures the background in windows of 1, 5 and 10 kb around each candidate peak, in the control sample. It uses the highest of these values, or the genome average if that is higher (Zhang et al. 2008). This makes MACS careful in regions where the input also has many reads.
MACS has two modes:
- Narrow mode finds sharp peaks. It suits transcription factors, sharp histone marks and ATAC-seq.
- Broad mode joins nearby enriched regions into wide domains. It suits broad marks such as H3K27me3.
SEACR: CUT&RUN and CUT&Tag
CUT&RUN and CUT&Tag have very little background. Peak callers made for ChIP-seq can then become too sensitive, because even a few reads look significant (Meers et al. 2019). SEACR was made for this kind of data (Meers et al. 2019).
SEACR splits the signal into blocks: stretches of the genome with continuous coverage. It adds up the signal in each block. Then it compares the blocks in the target with the blocks in the IgG control and chooses a threshold:
- stringent mode uses “the threshold at which the maximum percentage of target versus IgG signal blocks are retained” (Meers et al. 2019);
- relaxed mode uses a lower threshold, “halfway between the maximum and the ‘knee’” of the curve (Meers et al. 2019).
Without an IgG control, SEACR keeps a fixed top fraction of the blocks. For example, a threshold of 0.01 keeps the top 1 % of blocks by total signal.
GoPeaks: CUT&Tag for histone marks
GoPeaks was made for histone marks profiled with CUT&Tag (Yashar et al. 2022). It splits the genome into small bins, and uses a binomial test to find bins with more reads than the genome average. Then it joins significant bins that lie close together into peaks, and drops peaks narrower than 150 bp. The paper describes joining bins within 150 bp (Yashar et al. 2022); GoPeaks 1.0.0, the version EpiNexus runs, joins bins within 1,000 bp by default. For broad marks, the --broad option uses bins of 5,000 bp, moved 1,000 bp at a time. GoPeaks works on paired-end reads.
7.4 Choosing a peak caller
| Assay | Default in EpiNexus | Alternatives | Notes |
|---|---|---|---|
| ChIP-seq | MACS3 | — | Narrow or broad mode chosen from the mark |
| ATAC-seq | MACS3 | — | Narrow mode, with each summit reported (a peak can have several) |
| CUT&RUN | SEACR (stringent) | MACS3, GoPeaks | GoPeaks needs paired-end reads |
| CUT&Tag | SEACR (stringent) | MACS3, GoPeaks | GoPeaks needs paired-end reads |
Choose the peak caller in the run wizard, on the Parameters step (Peak Caller, Figure 7.1). The other settings are:
- MACS3. The q-value is 0.05 by default (Q-value Threshold). Broad marks use broad mode. ATAC-seq uses
--call-summits, which reports each summit inside a peak. EpiNexus gives MACS3 an effective genome size that matches the read length of each sample. If MACS3 cannot build its shift model, because there are too few strong peaks, EpiNexus runs it again without the model and with a fixed fragment size of 200 bp. - SEACR. Stringent mode. With an IgG control, SEACR compares the target with it. Without IgG, it keeps the top 1 % of signal blocks.
- GoPeaks. Default settings, and
--broadfor broad marks. - Min Peak Length. This field is empty by default. Each peak caller then uses its own minimum. For MACS3, that is the fragment size (MACS3 developers 2026). If you enter a length, MACS3 and GoPeaks do not report shorter peaks. SEACR has no such setting. (EpiNexus 1.2.2 and later. Before that, the field was saved but not used.)
After calling, EpiNexus removes peaks on unplaced or extra sequences and peaks that overlap the blacklist. Controls (input and IgG) are never called for peaks. They are only used as background.
Report the threshold your peaks were called at. MACS3’s own default q-value is 0.05 (MACS3 developers 2026). nf-core/cutandrun uses 0.01 when it calls CUT&Tag and CUT&RUN peaks with MACS (nf-core community 2026). EpiNexus uses the same values: 0.05 for ChIP-seq and ATAC-seq, and 0.01 for CUT&Tag and CUT&RUN when you choose MACS3. Since EpiNexus 1.2.2, runs started without the wizard use the same defaults.
7.5 Are the peaks reproducible?
A peak found in one replicate may be noise. A peak found in all replicates is more likely to be real. There are two ways to measure this.
IDR (irreproducible discovery rate) compares two replicates (Li et al. 2011). It ranks the peaks of each replicate by strength. Then it checks if strong peaks in one replicate are also strong in the other. Each peak gets an IDR value: the chance that it is not reproducible. ENCODE keeps peaks with an IDR below 0.05 (ENCODE Data Coordination Center 2026). IDR needs many weak peaks to model the noise. So peaks are first called with a loose threshold (p-value 0.01).
Overlap keeps regions where peaks from at least two replicates overlap. It is simpler and works for any number of replicates. It is the usual choice for histone marks.
EpiNexus checks reproducibility for each group of replicates (same mark and condition):
- IDR for transcription factors and ATAC-seq with exactly two replicates, called with MACS3. EpiNexus calls peaks again at p 0.01, keeps the top 300,000, and keeps peaks with IDR below 0.05.
- Overlap (at least two replicates) for histone marks, CUT&RUN, CUT&Tag and groups with more than two replicates.
By default, the differential analysis uses each sample’s own peaks. To use the reproducible peaks instead, choose Replicate-reproducible peaks (IDR / overlap) under Consensus peaks for differential. Check the size of the IDR set first: Chapter 16 shows an ATAC-seq dataset where it was much smaller than ENCODE expects.
7.6 Reading the result: GSE155100
In our ChIP-seq example, MACS3 found 54,372 to 66,954 peaks per H3K27ac library, with 36–45 % of the reads inside peaks (Table 6.2). H3K27ac is a narrow mark, so MACS3 used narrow mode. Each H3K27ac sample was compared with its own input.
Are these numbers reasonable? There is no correct number of peaks. It depends on the target, the tissue, the depth and the threshold. But you can compare with other studies of the same mark. For this dataset, ChIP-Atlas, a public database that reprocesses published data (Zou et al. 2024), called 21,000 to 27,000 peaks per sample with a much stricter threshold (q below 10−5). EpiNexus found 96–97 % of those peaks. The extra peaks come from the looser threshold of 0.05.
7.7 Pitfalls
Different peak callers and thresholds give very different numbers of peaks for the same data. Compare peak numbers only between samples analysed with the same caller and settings.
Without an input, open chromatin and regions with extra copies of DNA can look like peaks. If you have no input, treat the peaks with caution, and say so in your report.
More peaks in one condition does not mean more binding. Peak counts depend on depth and on thresholds. To compare conditions, use differential analysis (Chapter 9).