Chatomics Field GuideWhat They Don't Teach You

Sanity check · Single-Cell ATAC-seq

How to Call Peaks You Can Trust in Single-Cell ATAC-seq

Your top-scoring peaks might be centromeric repeats, not regulatory elements, and the pooled peak set your pipeline handed you is quietly erasing your rarest cell type.

By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 5 min read

You've extended your scRNA-seq project with 10x ATAC or Multiome, run Cell Ranger ATAC or your own MACS2 call, and pulled up the narrowPeak file. It looks fine at a glance: thousands of peaks, a q-value column, a signalValue column you can sort on. Then you sort by score and the top hits sit in satellite repeats and centromeric annotation, not promoters. There's no matched input sample anywhere in the pipeline, so nothing subtracted them out.

This is not a cosmetic problem. Every downstream step, the peak-by-cell matrix, TF-IDF and LSI, clustering, motif enrichment, inherits whatever is in that peak set. Artifact peaks in repeats and low-mappability regions contribute noisy, high-magnitude signal that clustering algorithms treat as real biological variation. You end up validating "cell states" that are actually mappability artifacts, or worse, missing the rare population you came here to find because it never crossed the peak-calling threshold in the pooled data.

In the next hour you can settle three questions: whether you need an input control at all (usually no), whether ENCODE blacklist filtering and a five-minute browser check catch what a control would have caught, and whether pooled peak calling is hiding a rare cell type you already suspect is there. The checks below are ordered so you can do the cheap ones first and only pay for a per-cluster MACS2 rerun once you know you need it.

What it looks like when it's happening

  • Sorting the narrowPeak file by score or -log10(q-value) puts satellite repeats, centromeric regions, or segmental duplications at the top instead of promoters.
  • A small number of peaks show enormous, flat, block-like pileups in a browser rather than the sharp, bimodal Tn5 insertion signal you'd expect at a real regulatory element.
  • FRiP (fraction of fragments in peaks) sits well below your expectation for the assay even though TSS enrichment score looks normal.
  • Two biologically implausible clusters share identical accessibility at the same small set of peaks, and those peaks turn out to be blacklist regions.
  • A rare population you can already see in the matched scRNA-seq or FACS data shows no distinguishing peaks when you cluster on the pooled, all-cell peak set.
  • Motif enrichment on a peak set returns low-complexity or repeat-associated motifs instead of the transcription factor motifs you expected for that cell type.
  • Peak widths in the raw MACS2 output range from under 100bp to several kilobases, with the wide ones sitting in the same repeat-rich regions flagged above.

Why it happens

Tn5 is not a clean reporter of transcription factor binding. It cuts nucleosome linker DNA and open chromatin indiscriminately, so the raw signal is already a mix of nucleosome positioning and regulatory information before you even call a peak. On top of that, Tn5 has real sequence bias: it prefers GC-rich sequence, CAG motifs, and open CpG islands, and its 9bp staggered cut leaves a reverse-complement palindrome pattern that can look like a transcription factor footprint when it is really an enzyme preference. None of this is corrected by a peak caller unless something explicitly models it.

The bigger problem for peak calling specifically is that most scATAC-seq experiments run with no input control. Input controls are the standard fix for exactly this kind of artifact in ChIP-seq, but they're rarely used in ATAC-seq because of the added sequencing cost. Without a control to subtract, regions with copy-number gain, generic open-chromatin bias, or poor mappability all get called as peaks with no counterbalance. Centromeric and satellite repeats are the classic offenders: multi-mapping and ambiguous alignment there produce large, flat pileups that MACS2 reads as a strong, broad signal, even though nothing about the underlying chromatin is unusually open.

The second failure mode is statistical, not artifactual. Cell Ranger ATAC and most default pipelines call peaks on all cells pooled together. A peak caller works on aggregate signal, so accessibility that is real but confined to a rare subpopulation, say a regulatory element open in 3% of cells, gets diluted into the pseudobulk profile of the whole sample and never crosses the significance threshold. The peak set you get back is biased toward whatever the most abundant cell types have open. This isn't a quality problem you can fix by tuning MACS2 parameters; it's a consequence of calling peaks before you know the population structure of the sample.

Both problems compound because scATAC-seq QC has a chicken-and-egg structure: you can't calculate FRiP, and you can't build the LSI-ready peak-by-cell matrix, until you have a peak set, but a trustworthy peak set often requires a first-pass clustering you can only get after building some matrix. Most workflows resolve this by starting from an intentionally provisional pooled peak set, cleaning it with blacklist filtering, and only refining to per-cluster peaks once a first clustering exists.

The checks

Run them in order. Each one tells you what healthy looks like and what the problem looks like.

0/7 checked · saved in this browser

  1. Download the ENCODE blacklist bed file for your exact genome build (the list differs between hg19, hg38, mm10, etc.) and remove any peak that overlaps it, before this peak set touches your peak-by-cell matrix.

    bash
    bedtools intersect -v -a peaks.narrowPeak -b hg38-blacklist.v2.bed.gz > peaks.blacklist_filtered.narrowPeak
    Healthy
    A small fraction of peaks drop out, typically well under a few percent, and the removed peaks are not your top-scoring ones.
    Red flag
    More than a handful of your highest-scoring peaks by q-value disappear, or the blacklist overlap is large enough that it visibly changes the shape of your peak count distribution.
  2. Reread the exact MACS2 command that produced your peaks (control used or not, --nomodel/--shift/--extsize for Tn5 cut-site centering, and the q-value cutoff) and plot a histogram of peak widths from the narrowPeak file.

    Healthy
    A documented, deliberate q-value threshold was used, and peak widths cluster in a consistent range (roughly 150 to 500bp before any fixed-width standardization, or a uniform width if you're using ArchR's iterative peak set).
    Red flag
    No explicit q-value threshold in the command (defaults silently pass borderline calls), or peak widths span from under 100bp to multiple kilobases with the wide ones overlapping repeat annotation.
  3. Sort the narrowPeak file by signalValue or -log10(q-value) descending, load the top peaks alongside a RepeatMasker or mappability track in IGV, and look at the shape of the pileup, not just the score.

    Healthy
    Top peaks sit at promoters or enhancers of genes you'd expect to be expressed, with a sharp, defined summit and mappability near 1.
    Red flag
    Top peaks land in satellite repeat, centromeric, or segmental duplication annotation, with flat, blocky coverage instead of a defined summit.
  4. Pull up promoters for a few housekeeping genes and, if you have matched RNA or FACS data, marker gene loci for the cell types you already know are in the sample, and confirm a peak is called there.

    Healthy
    Sharp, well-defined peaks at these loci regardless of how the rest of QC looks.
    Red flag
    No peak called at an obviously open promoter, which points to a broken pipeline step or an overly strict threshold, not a blacklist problem.
  5. For each cell, calculate the fraction of its fragments that fall inside the peak set (using Signac's FRiP function or ArchR's per-cell fragment-in-peak calculation), and look at the distribution across all cells, not just the mean.

    Healthy
    A unimodal distribution with a median in the range typical for your kit and library, consistent with a healthy TSS enrichment score for the same cells.
    Red flag
    Most cells sit well below the expected range, or the distribution is bimodal with a large low-FRiP mode that doesn't correspond to a known low-quality batch.
  6. Compute the ENCODE-standard TSS enrichment score (signal around a reference TSS set versus flanking background) for the same cells and compare it to FRiP.

    Healthy
    TSS enrichment and FRiP move together: cells with good TSS enrichment also have reasonable FRiP.
    Red flag
    TSS enrichment looks healthy but FRiP is low across the board, which points specifically at the peak set (artifacts, missing regions) rather than at library or cell quality.
  7. After a first-pass clustering (built on depth-independent LSI, with the first component dropped), call peaks separately per cluster and merge into a unified set, then check whether any cluster you suspect is a rare population gains peaks it didn't have in the pooled call.

    r
    peaks <- CallPeaks(atac, group.by = "seurat_clusters")
    # or, for the iterative overlap + fixed-width approach:
    archr_project <- addReproduciblePeakSet(archr_project)
    Healthy
    Rare clusters gain distinct peaks not present, or not significant, in the pooled peak set.
    Red flag
    The per-cluster and pooled peak sets are nearly identical, meaning rare-population accessibility is still being diluted away even after clustering.

What to do about it

Apply ENCODE blacklist filtering as a mandatory last step

When: Every scATAC-seq peak set, regardless of whether you have an input control.

Run bedtools intersect -v between your peak set and the ENCODE blacklist matched to your exact genome build immediately after peak calling, before building the peak-by-cell matrix.

Caveat: The blacklist only removes known recurrent artifacts; it won't rescue subtler mappability problems, and using the wrong genome build's list either misses artifacts or strips real peaks.

Skip the input control, but use Tn5-aware MACS2 settings and a documented q-value cutoff

When: Standard scATAC-seq where a matched input isn't affordable, which is most experiments.

Call peaks with the Tn5 cut-site shift model instead of the ChIP fragment model, and set an explicit q-value threshold rather than relying on defaults, so borderline calls don't silently make it into your peak set.

Caveat: Without a control you can't formally estimate a false discovery rate the way ChIP-seq can; blacklist filtering and manual spot checks only partially substitute, and residual artifacts persist especially in samples with copy-number changes.

Call peaks per cluster and merge into a unified peak set

When: Your sample has rare populations (roughly under 5 to 10 percent of cells) or heterogeneous tissue where pooled calling is plausibly collapsing signal.

Cluster first on depth-corrected LSI (drop the first component), then call peaks per cluster with Signac's CallPeaks(group.by=...) or ArchR's addReproduciblePeakSet, which does iterative overlap removal and produces fixed-width, non-overlapping peaks.

Caveat: You need a first-pass clustering before this peak set exists, so use a coarse scaffold (bins or gene activity scores) and expect to recluster once the refined peaks are ready; this also multiplies compute cost, one MACS2 run per cluster.

Recompute FRiP and TSS enrichment after every peak set change

When: Any time you apply blacklist filtering or move from pooled to per-cluster peaks.

Rerun per-cell FRiP and TSS enrichment against the new peak set and use both together, not either alone, as your cell-quality gate.

Caveat: Reusing thresholds calibrated on the old, artifact-inflated peak set will misclassify real low-depth cells, especially rare-population cells, as low quality.

Cross-reference suspiciously hot regions against a CNV estimate

When: Tumor samples, cell lines, or any sample where large-scale copy-number gain is plausible.

If you have matched WGS, bulk data, or an inferCNV-style estimate from RNA, check whether your strongest peaks sit in amplified regions; accessibility there scales with copy number, not chromatin state.

Caveat: Without an independent CNV estimate you can only flag this by pattern (implausibly broad, block-like signal), not confirm it.

When not to "fix" it

Not every unusually strong or broad peak is an artifact. Genuine super-enhancers and the promoters of extremely highly transcribed genes can produce very strong, broad accessibility signal that is biologically real, and blindly filtering anything "too big" will throw those out along with the repeats. Check gene annotation before you filter: a strong peak co-located with a real, ubiquitously expressed promoter is not the same finding as a strong peak sitting in RepeatMasker satellite annotation with no gene nearby.

If you're studying a tumor sample where clonal copy-number gain is the biology of interest, don't strip out accessible regions in amplified segments; that signal is the point of the experiment, not noise to be removed. And never apply a blacklist built for the wrong species or genome build. A human ENCODE blacklist run against a mouse alignment, or an hg19 list applied to hg38 coordinates, will misfire and remove real peaks instead of artifacts.

Five things experienced analysts do here

  1. Run the blacklist intersection as the literal last step before building any matrix, and rerun it every time you change genome build or add a new sample, don't assume last quarter's filtered file still applies.
  2. Sort every new peak set by score once and open the top 30 in IGV before trusting a QC summary; this five-minute manual check catches artifact peak sets that summary statistics alone miss.
  3. Never reuse FRiP or TSS enrichment thresholds from a previous project or a paper without recomputing them on your own peak set; an artifact-inflated peak set shifts the whole distribution and makes borrowed thresholds meaningless.
  4. Treat your first pooled MACS2 call as a scaffold for clustering, not a final answer; budget time for a second, per-cluster peak-calling pass before you generate any figure meant for publication.
  5. An antibody-free assay doesn't mean the biology check disappears: validate that peaks land at loci you already expect to be open, the same way you'd sanity-check a ChIP-seq peak against a known target gene.

Questions people ask

Do I need an input control for scATAC-seq peak calling?

No, most scATAC-seq experiments run without a matched input control, mainly because of the added sequencing cost. Compensate instead with ENCODE blacklist filtering, an explicit MACS2 q-value threshold, and a manual check of your top-scoring peaks in a browser.

Why are my strongest MACS2 peaks sitting in centromeres or repeats?

Repetitive, low-mappability DNA produces ambiguous or multi-mapping alignments that pile up regardless of true chromatin state, and with no input control to subtract that signal, MACS2 reads it as a strong peak. Filtering against the ENCODE blacklist and spot-checking your top peaks in a browser catches most of these.

What FRiP score is good enough for scATAC-seq?

There's no single universal cutoff; judge your own distribution rather than a borrowed number. If most cells sit well below your library's typical range, or FRiP is low while TSS enrichment looks fine, that points at a peak-set problem rather than a library-quality problem.

Should I call peaks on all cells pooled, or per cluster?

Pooled peak calling, which is what Cell Ranger ATAC does by default, favors accessibility from the most abundant cell types and can miss rare-population-specific regulatory elements entirely. Calling peaks per cluster with Signac's CallPeaks or ArchR's addReproduciblePeakSet and merging the result recovers those sites.

Narrow or broad peaks for ATAC-seq?

There isn't a settled, universal answer in the literature for scATAC-seq specifically. Default to narrow peak calling to match the sharp Tn5 cut-site signal and the ENCODE ATAC pipeline's MACS2 defaults, and validate against known marker loci regardless of which mode you choose.

Related pages

Related reading on the blog

Sources

  1. From reads to insight: a hitchhiker's guide to ATAC-seq data analysis — Input controls, repeat masking, and MACS2 usage in ATAC-seq peak calling
  2. The ENCODE Blacklist: Identification of Problematic Regions of the Genome — Origin and rationale of the ENCODE blacklist regions applied here
  3. Calculate scATACseq TSS enrichment score — ENCODE-standard TSS enrichment calculation used alongside FRiP
  4. Quality control of single-cell ATAC-seq data without peak calling using Chromap — FRiP's dependence on a pre-defined peak set and its limits as a QC metric

Part of the Peak calling controls and blacklists series.