Implementationstooling reviewed 2026-09
21 Epigenomics and Chromatin: Implementations
Everything on this page is generated from tools/epigenomics.json. If something here is wrong or out of date, the fix belongs in that file, not in this prose.
For why a peak is a claim rather than an observation, and why a methylation ratio is a different kind of measurement from a read count, see Section 20.1.
This page is organized by the two halves of Section 20.2. The first three sections are the enrichment assays — ChIP-seq, ATAC-seq and the targeted methods — which share tools almost completely. The last two are DNA modification, which shares none of them.
21.1 Calling peaks
One caller is the default for both ChIP-seq and ATAC-seq, and the alternatives below exist because its model makes an assumption that some experiments break. The assumption is the one Section 20.3 identifies as the whole method: that enrichment is punctate, and that the surrounding genome carries enough background to estimate an expectation from.
Broad histone marks break the first half. The targeted low-background assays break the second.
The default for both ChIP-seq and ATAC-seq, and the reason is the local background model described in Section 20.3 rather than the peak-finding itself: it estimates the expected count from the surrounding window instead of from a genome-wide average, which is what stops copy-number and mappability variation from reading as enrichment. Everything downstream in this field expects its output, so choosing it also buys the annotations, the QC metrics and the pipelines. (Zhang et al. 2008)
The narrow-peak model assumes punctate enrichment. Broad marks such as H3K27me3 or H3K36me3 need the broad mode at minimum, and the alternatives below exist because even that is a compromise.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| epic2 | activelast release 2025-12 | pyC++ |
A reimplementation of SICER, whose model is built for domains rather than for points: it scores islands of enrichment separated by gaps, which is a better description of a spreading histone mark than a peak with a summit. Choose it when the mark genuinely forms domains and MACS3’s broad mode is fragmenting them. (Stovner and Sætrom 2019) |
| Genrich | activelast release 2026-08 | C |
Takes replicates jointly rather than calling each sample and intersecting afterwards, which is the right shape for the problem Section 20.5 describes: a peak just under threshold in one replicate should not vanish from the result. It also removes duplicates and filters mitochondrial reads in the same pass, which for ATAC-seq is most of the preprocessing. |
| SEACR | dormantlast release 2020-01 | shR |
Written for the sparse, near-zero-background data CUT&RUN produces, where the assumption a ChIP-seq caller makes — that background is dense enough to model locally — is simply false. It thresholds on the empirical distribution of contiguous signal blocks instead, which is why it finds peaks in libraries too shallow for MACS3 to model at all. (Meers et al. 2019) |
NoteMACS2 is not a separate entry
MACS2 and MACS3 are the same project — one repository, one lineage — so the registry gives them one record, at the current version. You will meet MACS2 constantly in published methods and in pipelines that pinned it, and it is not broken. MACS3 is a Python 3 port with the same model and the same sub-commands, so moving is usually a matter of the command name. Pipelines still running MACS2 are behind rather than wrong.
21.2 Signal and quality control
Both entries here answer the question Section 20.4 says to ask first: did this library enrich for anything. That question is worth answering before peak calling rather than after, because a peak caller will produce output either way and the output looks the same.
The first is the general toolkit — coverage tracks, matrices over region sets, and the enrichment plots. The second computes the two metrics that are specific to ATAC and that nothing general will give you.
The step between an alignment and anything you can look at: normalized coverage tracks, matrices over sets of regions, and the diagnostic plots that answer whether an experiment worked before you spend time on its peaks. Its fingerprint and correlation plots are the practical form of the enrichment question Section 20.4 poses, and they are cheaper to run than any of the analyses they would save you from. (Ramírez et al. 2016)
A 4.0.0 rewrite with a Rust backend has been in development on a branch through 2026 but is unreleased, and the default branch has not moved since the last release. Nothing is broken; the version you install is simply older than the work you can see happening.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| ataqv | dormantlast release 2023-02 | C++ |
Reports the two metrics that actually diagnose an ATAC library — TSS enrichment and the fragment-length distribution — and puts them in an interactive comparison across samples rather than one report per file. The nucleosomal periodicity in that distribution is the clearest single indicator that the transposition worked, and it is not something a generic read-QC tool computes. (Orchard et al. 2020) |
WarningNormalize the track before you compare tracks
A coverage track made without depth normalization is a picture of how deeply that sample was sequenced. Comparing two of them in a browser and concluding that one condition has more signal is one of the easiest mistakes in this field to make and one of the hardest to see in a figure. Normalize to a stated quantity, say which one, and be aware that every option assumes something — a normalization to total depth assumes the two libraries have comparable signal-to-noise, which is exactly what a failed experiment violates.
21.3 Testing for differences
Section 20.5 is the reason this is a section and not a bedtools intersect. Two conditions, peaks called in each, non-overlapping intervals declared condition-specific: the result is dominated by regions that sat near the threshold, and it gets worse as the libraries differ in depth.
The three entries take different routes past that. One tests windows across the genome and never needs a peak set. One uses a peak set but recounts every sample over a shared consensus before testing. One changes the question entirely, aggregating over motifs instead of testing regions.
Counts reads in fixed windows across the genome and tests those, which sidesteps the circularity Section 20.5 names: a peak set called from the same data you are testing has already used the signal you are about to measure. It hands the counts to edgeR, so the negative binomial machinery and the designs from Section 16.6 apply unchanged. (Lun and Smyth 2016)
Windows mean many more tests and a harder multiple-testing problem, which it handles by aggregating windows into regions before correction. That clustering step is a modelling choice and it deserves reading about rather than accepting.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| chromVAR | active | RC++ |
Asks a different question from the two above: not which regions differ, but which transcription factor motifs are associated with the variation in accessibility across samples. Aggregating over all peaks carrying a motif is what makes the question answerable in sparse data where no individual peak has enough counts. (Schep et al. 2017) |
| DiffBind | active | RC++ |
The peak-set route done carefully: it builds a consensus set across samples, recounts every sample over it, and then tests. That recounting is the part people skip when they do this by hand, and skipping it is what turns a difference in peak-calling sensitivity into a false differential result. (Ross-Innes et al. 2012) |
NoteThis is the RNA-seq problem with different regions
Once you are counting reads in a fixed set of intervals across replicated samples, the statistics are the ones Section 16.6 describes and the same packages are underneath. That is a convenience and a warning: the models expect raw counts and a design, so the input is a count matrix, not a set of normalized coverage values, and the same rule about not running several tests and choosing the friendliest applies here.
21.4 Measuring methylation
Two cards, because there are two genuinely different measurements here and not two implementations of one. The first is the bisulfite route of Section 20.6: convert the DNA, align in a three-letter alphabet, count C against T. The second involves no conversion at all — a nanopore read carries per-base modification probabilities alongside its sequence, so genotype and methylation come from one alignment.
The choice between them is made when the library is prepared, not when the analysis starts. The table holds two alternatives inside the bisulfite route.
Does the three-letter alignment of Section 20.6 and the methylation call in one run, and reports per-cytosine counts split by CpG, CHG and CHH context rather than collapsing them. Keeping the contexts separate matters because non-CpG methylation is a real signal in plants and in stem cells and an alignment artifact almost everywhere else, so the split is a diagnostic and not just a detail. (Krueger and Andrews 2011)
It drives an external aligner, so the aligner's own index and memory costs are yours as well. The bisulfite index is built once per genome and is roughly twice the size of a standard one, since both converted strands are indexed.
The tool for the other measurement entirely: nanopore reads carry per-base modification probabilities in the BAM’s MM and ML tags, so no chemical conversion and no separate library are involved, and the read still carries its own sequence. This turns those per-read probabilities into per-site counts, and it is the only maintained tool doing that job on the standard tags. (Oxford Nanopore Technologies 2026)
The license is Oxford Nanopore's own rather than an OSI-approved one; read it before redistributing. Its defaults also apply a probability filter whose threshold materially changes the calls, and the right value is dataset-dependent.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| BISCUIT | activelast release 2026-08 | C |
Calls genetic variants and methylation from the same bisulfite library rather than treating the genetics as a nuisance. That matters for the specific failure Section 20.6 describes: a C-to-T SNP is indistinguishable from an unmethylated cytosine unless something models both, so a study of allele-specific methylation needs this and a converted-genome aligner alone cannot give it. (Zhou et al. 2024) |
| MethylDackel | dormantlast release 2021-09 | C |
Extracts per-base methylation from a bisulfite BAM produced by any converted-genome aligner, which decouples the call from the alignment and is the reason to use it: you keep whichever aligner your project already validated. Its bias-plotting mode is the useful part, showing the methylation-by-read-position curve that reveals the end-repair artifact you would otherwise average into every result. |
WarningEstimate the conversion rate, every run
Incomplete bisulfite conversion reads as methylation, uniformly, in a direction that always inflates. None of the tools above will warn you: an unconverted cytosine is a methylated cytosine as far as the counting is concerned. Estimate the rate from something known to be unmethylated — a spike-in, the mitochondrial genome, or non-CpG context in a tissue where that is safe — and report it. A dataset whose conversion rate is unknown has an unknown bias of unknown size.
21.5 Differential methylation
The measurement is per cytosine and the finding is a region, and Section 20.7 is about how much of this field loses its error control in that gap. The card is the entry that treats the region as the unit of inference from the start. The table holds the two that test sites and then assemble, which is a defensible thing to do as long as the region-level numbers are read as descriptive.
All three are R packages, which is not an accident — this is a modelling problem rather than a throughput problem, and the counts are small by the time they arrive.
Tests regions as units and controls the false discovery rate at the region level, which is the specific thing Section 20.7 says most of this field gets wrong: a region assembled from per-cytosine tests and given the smallest p-value inside it has no valid error rate at all. It gets there by permuting the design to build a null distribution for whole regions, so the correction matches the object being reported. (Korthauer et al. 2019)
The permutation is the cost, in runtime and in needing a design with enough samples to permute. Last commit May 2025; the package is current in Bioconductor and the method is settled, but it is being kept rather than developed.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| DSS | active | RC |
Models the count nature of the measurement — methylated reads out of total reads at a cytosine — with a beta-binomial and shrinks the dispersion across sites, so a cytosine covered by eight reads is not treated as equally certain as one covered by eighty. That coverage-awareness is what a ratio-based method throws away, and it is the most common way a differential methylation result goes wrong. (Feng et al. 2014) |
| methylKit | active | RC++ |
Covers the whole path from an aligner’s methylation output to annotated differential results in one package, including the tiling, filtering and correlation steps people otherwise write themselves. It is the pragmatic choice for targeted designs such as RRBS, where the sites are few enough that per-site testing is the honest unit anyway. (Akalin et al. 2012) |
21.6 What is not here
HOMER is the conspicuous absence among the peak callers, and it is absent for the reason ANNOVAR is absent from the variant annotation part rather than for any judgment about the software. It has no public source repository; it is installed by running a Perl script that downloads from a university web page, and its versions are announced on that page rather than published to any feed. Every field a record would need is therefore unavailable, and both honest options are false — latest_version: null renders as “unrefreshed”, which claims nobody looked, and a hand-typed version is exactly the practice the refresh pass exists to prevent. So the position goes here: its motif-finding is still widely used and its peak caller is a reasonable one, an existing pipeline built on it is not broken, and for a new one the callers above are locatable, versioned and citable.
bwa-meth is missing for a related but distinct reason, and it is worth stating precisely because the record looks obtainable. It is a well-regarded and still-maintained bisulfite aligner — commits through 2026 — and its current version, 0.2.9, exists as a git tag. The problem is that the project restarted its version numbering: the tags run v0.01 to v0.10 through 2013 and 2014, then begin again at v0.2.0 in 2016. Under any ordering rule, 0.10 is greater than 0.2.9, so the newest tag by version is one from June 2014 and the newest tag by date carries a number that has been superseded. No automatic rule can resolve that, because the ambiguity is real. The book’s view is that the tool is a sound choice — pairing it with the extractor listed above is a standard and well-tested route — but it cannot carry a record whose version and date are both true.
Motif databases and motif enrichment. The substance of a motif analysis is which collection of position weight matrices you searched, and those are versioned datasets on their own release cycles. This is the same shape as the annotation databases in the variant part: ranking the search programs would put a review date on the wrong object, and the entry above that works on motifs takes the database as an argument rather than supplying one.
Footprinting — inferring bound factors from the fine structure of accessibility within a peak — is absent because the methods disagree with each other substantially and the quantity they infer depends strongly on the transposase’s own sequence bias. A survey with a “start here” verdict would overstate how settled it is.
Hi-C and chromosome conformation. Chromatin architecture is measured as a contact matrix between pairs of loci, not as a signal along the genome. Almost nothing on this page applies: the file formats are different, the normalization problem is different, and the unit of a finding is a pair of intervals rather than one. It is a subject, not a section.
Single-cell ATAC and multiome. Excluded on the same grounds part 9 excludes the other single-cell modalities. They share the peak concept and almost nothing else — the matrix is near-binary, the QC is different, and the questions asked of the result are different.
Consortium pipelines. The ENCODE processing pipelines are the reference implementation of much of what this page describes, and they are workflows rather than tools: they are written in workflow languages and run by the engines surveyed in Section 23.1. Reading one is an excellent way to see these tools composed correctly.
21.7 Language coverage: the most R-dependent part of the book
Seven of the sixteen entries here are R packages, which is more than any other part — more than bulk RNA-seq and more than single-cell — and the split is clean rather than incidental. Everything that touches every read is compiled: the peak callers are C, the methylation callers are C and Rust, and the one Python entry is Python around compiled I/O. Everything that fits a model to a count matrix is R and Bioconductor, without a competitive Python alternative in any of the three modelling sections.
That asymmetry is worth planning around, because it is stronger here than in the expression parts, where the single-cell chapter could offer two complete stacks. There is no Python route through differential binding or differential methylation that this book would recommend. A project doing chromatin analysis in Python will be calling peaks and building matrices in Python and then moving to R for the part where the conclusions are made.
The other thing the chips show is Rust arriving in earnest. Two of this part’s recommended tools are Rust, and one of them got there by being rewritten: Bismark’s Perl implementation was replaced in 2026 by a Rust suite reported byte-identical on its default path, which is the strongest form the argument can take — a rewrite that changes performance and nothing else. The general signal tool has a Rust backend under development on a branch as well. This is the first part of the book where an established tool has migrated language rather than a new tool arriving in a new one.
NoteWhat this means in practice
Budget for the handoff. The natural pipeline here crosses languages once, at the point where intervals and counts become a model, and the crossing is cheapest if you plan it: write the count matrix and the sample table to disk as a deliberate boundary rather than reaching for an in-process bridge. Everything before that boundary is a command line, and everything after it is R.