Why do I have zero significant peaks after FDR correction in ATAC-seq?

Most often because you tested 100,000+ peaks instead of ~20,000 genes, so the Benjamini-Hochberg threshold is far stricter than what you're used to from RNA-seq. Check the raw p-value histogram first: if it's flat, there's genuinely little signal at your current power; if it spikes near 0 but padj is still all NA or high, look at independent filtering and how many peaks you actually tested.

What's the difference between padj and pvalue in DESeq2 for ATAC-seq?

pvalue is the per-peak test result before any correction; padj (the BH-adjusted value, DESeq2's estimate of the FDR) accounts for the total number of peaks tested and for independent filtering. Always report padj for significance calls, but look at pvalue and log2FoldChange together to judge whether a near-miss peak is worth following up.

Should I use Benjamini-Hochberg or Bonferroni for ATAC-seq peaks?

Use Benjamini-Hochberg (or DESeq2/edgeR's built-in BH adjustment). Bonferroni controls the family-wise error rate, which is appropriate for something like GWAS where a single false positive is very costly, but it is far too conservative for exploratory differential accessibility work across 100,000+ peaks.

What is independent filtering in DESeq2 and why does it matter for ATAC-seq?

Independent filtering removes low-count peaks before FDR adjustment using a criterion (mean normalized count) that is independent of the test statistic, which increases power without inflating false positives. In ATAC-seq it matters more than in RNA-seq because a much larger fraction of called peaks are low-signal MACS2 calls that would otherwise dilute the BH correction for every real peak.

Is FDR < 0.05 the right cutoff for ATAC-seq differential peaks?

It's a reasonable default, not a law. With only 2-3 replicates per condition and modest accessibility changes, FDR < 0.05 combined with a fold-change filter can leave you with nothing to work with even when real biology is present. Decide your cutoff before looking at the results, and treat effect size (log2FC) as at least as important as the adjusted p-value.

Sanity check · ATAC-seq

How to Handle Multiple Testing and FDR in ATAC-seq

Q: Is FDR < 0.05 the right cutoff for ATAC-seq differential peaks?

It's a reasonable default, not a law. With only 2-3 replicates per condition and modest accessibility changes, FDR < 0.05 combined with a fold-change filter can leave you with nothing to work with even when real biology is present. Decide your cutoff before looking at the results, and treat effect size (log2FC) as at least as important as the adjusted p-value.

Zero significant peaks after BH correction is usually the statistics doing exactly what you asked, not a broken pipeline.

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

You called peaks with MACS2, built a consensus set, counted fragments, ran DESeq2 or edgeR, and got a results table where padj is NA or above 0.9 for every single peak. RNA-seq on the same samples gave you 200 DE genes. Something feels wrong, so the instinct is to lower the FDR threshold to 0.1, drop the "bad" replicate, or quietly switch to raw p-values in the paper draft.

That instinct is the expensive mistake. A consensus ATAC-seq peak set is often 100,000 to 300,000 regions, five to fifteen times the number of genes in a typical RNA-seq comparison. Benjamini-Hochberg correction gets stricter as the number of tests grows, so the same effect size and the same p-value that would survive correction in RNA-seq can vanish in ATAC-seq purely because you tested more things. Separately, if your peaks were called per-condition instead of on a pooled, condition-blind library, you've already violated the independence assumption FDR control relies on, and no threshold will fix that.

This page gives you a short sequence of checks to tell which of those two problems you have, in under an hour, before you touch a single cutoff.

What it looks like when it's happening

results(dds) or topTags() shows padj/FDR equal to NA or above 0.9 for every single peak in the table
summary(res) reports a large fraction of peaks removed by independent filtering, sometimes more than half the consensus set
a raw p-value histogram (hist(res$pvalue)) is flat across 0 to 1 instead of spiking near zero
the consensus peak set has 100,000 to 300,000 regions, far more than the ~20,000 genes tested in a parallel RNA-seq comparison on the same samples
MA plot shows log2 fold changes shrunk tightly around zero even for peaks with visibly different read pileups in a genome browser
edgeR's topTags() and DESeq2's results() disagree sharply on the number of peaks below FDR 0.05 from the same count matrix
peaks were called separately per sample or per condition, then combined into a consensus set by union or overlap

Why it happens

Benjamini-Hochberg adjusts each p-value based on its rank among all p-values tested, so the adjusted threshold gets stricter as the total number of tests (m) grows. A gene at rank 50 out of 20,000 RNA-seq genes clears a much looser bar than a peak at rank 50 out of 200,000 ATAC-seq peaks, even with an identical raw p-value and identical effect size. Since ATAC-seq consensus peak sets routinely run five to fifteen times larger than a gene-level RNA-seq matrix, the same biology produces far fewer FDR-significant hits by construction, not because the accessibility changes are smaller.

A second, more damaging mechanism is peak-calling circularity. MACS2 calls a large fraction of ATAC-seq peaks near the noise floor because Tn5 cuts both nucleosome linker DNA and genuinely accessible regulatory DNA, producing a continuum of signal strength rather than a clean on/off call. If peaks are called per-sample or per-condition and then unioned into a consensus set, peaks that happen to look stronger in one condition are more likely to survive into the consensus set for that condition, which means the filtering statistic (was this region called a peak) is no longer independent of the test statistic (does this region differ between conditions) used downstream. That correlation inflates the false discovery rate silently: BH correction assumes the set of hypotheses was fixed before you looked at group differences, and pooling peaks this way breaks that assumption regardless of which tool runs the test afterward.

A third contributor is that most ATAC-seq peaks in a consensus set are real Tn5 insertion events but biologically uninteresting: promoters and repeat-adjacent open regions that don't move between conditions. Counting these low-variance, low-signal peaks alongside real differential regions dilutes the test set and lowers power for everything, which is exactly what independent filtering (DESeq2's built-in filtering, or filterByExpr in edgeR) is designed to remove before BH correction is applied, not after.

Finally, accessibility changes are frequently partial rather than binary. A transcription factor losing half its binding events at a locus produces a 30-50% change in Tn5 insertion frequency, not the 4-8 fold changes common in strongly induced RNA-seq genes. With 2-4 replicates, the negative binomial dispersion estimate for a peak with that modest a fold change often isn't precise enough to push the p-value low enough to survive correction against 100,000+ competitors, even when the effect is real.

The checks

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

0/8 checked · saved in this browser

01Count how many peaks you actually tested
Check the number of rows going into DESeq2 or edgeR before any filtering: nrow(counts) or dim(y$counts). Compare it to a typical RNA-seq gene count (roughly 20,000).
r
```
nrow(counts)
# compare to ~20,000 protein-coding genes in RNA-seq
```
Healthy
A few tens of thousands after reasonable pre-filtering; the raw consensus set can be much larger before filtering.
Red flag
You're running BH correction on 150,000+ unfiltered peaks, most of which have near-zero counts across all samples.
02Read the independent filtering summary
Run summary(res) on the DESeq2 results object and look at the line reporting how many peaks were removed by independent filtering (low counts) versus how many were tested for differential accessibility.
r
```
res <- results(dds, contrast=c('condition','treated','control'))
summary(res)
```
Healthy
A moderate fraction filtered (commonly 20-50%), with a clearly stated number remaining that FDR was computed over.
Red flag
Almost nothing filtered (most junk peaks retained, diluting the correction) or almost everything filtered (real peaks discarded, leaving too few tests to detect anything).
03Plot the raw p-value histogram before looking at padj
Histogram the unadjusted p-values from the results table. This is the fastest diagnostic for whether there is any signal at all versus a model or filtering problem.
r
```
hist(res$pvalue, breaks=50, main='raw p-values')
```
Healthy
A spike near 0 (real differential peaks) sitting on top of a roughly uniform background (the null peaks).
Red flag
Perfectly flat across 0 to 1 (little to no detectable signal at current power) or a spike near 1 (a covariate or batch term is misspecified in the design formula).
04Verify peak calling was blind to condition
Check the peak-calling command in your pipeline log. Peaks should come from a single pooled BAM merging all samples with equal contribution, not from calling MACS2 separately per sample or per condition and then taking a union.
bash
```
samtools merge -h sample1.bam sample2.bam sample3.bam ... > pooled.bam
macs2 callpeak -t pooled.bam -f BAM --nomodel --shift -100 --extsize 200 -q 0.05 --name pooled_peaks
```
Healthy
One consensus peak set called on the pooled, condition-blind BAM, then used to count fragments per sample.
Red flag
Separate MACS2 runs per sample or per condition, combined by union or overlap into the consensus set used for counting.
05Confirm your low-count filter isn't condition-aware
Check whether your pre-filtering step uses an overall abundance measure (filterByExpr with the full design, or mean/aveLogCPM across all samples) rather than a per-condition or per-library threshold that could preferentially keep peaks that already look different between groups.
r
```
filter <- filterByExpr(y, design=design)
y_filtered <- y[filter,]
```
Healthy
Filtering criterion is symmetric across conditions and computed before the differential test.
Red flag
Peaks were kept because they passed a threshold in one condition's samples but not the other's, which reintroduces the same circularity as condition-aware peak calling.
06Look at effect sizes among the near-misses
Regardless of padj, sort by raw p-value and inspect log2FoldChange for peaks with pvalue < 0.05 but padj > 0.05. This tells you whether you have real, sizeable effects that correction is (correctly) discounting for lack of power, or nothing but noise-level fold changes.
r
```
subset(res[order(res$pvalue),], pvalue < 0.05 & abs(log2FoldChange) > 1)
```
Healthy
A handful of peaks with |log2FC| > 1 and small raw p-values clustering just above the FDR line, biologically plausible as underpowered true positives.
Red flag
Near-miss peaks all have |log2FC| well under 0.5, meaning there's no meaningful effect size hiding behind the correction.
07Cross-check DESeq2 against edgeR on the same count matrix
Run both tools on the identical filtered peak-by-sample count matrix and design, then compare the number of peaks at FDR < 0.05 and the overlap between the two significant peak lists.
r
```
y <- DGEList(counts=counts, group=groups)
y <- calcNormFactors(y, method='TMM')
design <- model.matrix(~groups)
y <- estimateDisp(y, design=design)
fit <- glmFit(y, design=design)
lrt <- glmLRT(fit, contrast=c(-1,1))
topTags(lrt, n=Inf, adjust.method='BH', p.value=0.05)
```
Healthy
Broadly similar peak counts and substantial overlap in the significant lists between DESeq2 and edgeR.
Red flag
One tool returns hundreds of hits and the other returns zero from the identical input matrix, which points to a dispersion estimation or design formula problem rather than a real biological absence of signal.
08Rule out a failed library masking real signal
Before concluding the effect is genuinely small, check each sample's fragment-size histogram, TSS enrichment score, and fraction of reads in peaks (FRiP). A failed library can suppress true differential signal by adding noise that a downstream DE test can't distinguish from biology.
Healthy
Fragment-size histogram shows the expected nucleosome-periodic pattern, TSS enrichment and FRiP within your lab's usual range across all samples in the comparison.
Red flag
One or more samples in the comparison have a flat fragment-size distribution or low TSS enrichment, meaning the zero-hit result may reflect a bad library, not the true condition effect.

What to do about it

Apply independent filtering before correction, not after

When: A large share of your consensus peaks have near-zero counts across all samples, diluting the BH correction for the peaks that matter.

Use DESeq2's built-in independent filtering (on by default in results()) or edgeR's filterByExpr(y, design=design) on the full design matrix, applied once before any per-condition comparison, then run the test only on peaks that pass.

Caveat: The filtering criterion must depend only on overall abundance, never on which condition a peak looks stronger in, or you reintroduce the same circularity you're trying to fix.

Recall peaks on a pooled, condition-blind BAM

When: Your consensus peak set was built from separate per-sample or per-condition MACS2 runs and then unioned.

Merge all sample BAMs with equal contribution into one pooled BAM, call peaks once on the pool (macs2 callpeak -t pooled.bam -f BAM --nomodel --shift -100 --extsize 200 -q 0.05), then recount fragments per sample against that single consensus set.

Caveat: Requires re-running peak calling and re-counting for every sample; budget the compute and re-run downstream DE analysis from scratch, since the old counts are built on a biased peak set.

Report by effect size when FDR-significant hits are near zero

When: Raw p-values spike near 0 for a small set of peaks with real fold changes, but few or none survive BH correction given the size of the peak set.

Rank peaks by log2FoldChange among those with raw p < 0.05, treat this as an exploratory hit list for motif enrichment or follow-up validation, and state the relaxed criteria explicitly in the methods rather than silently lowering the FDR cutoff after seeing the result.

Caveat: This list carries a higher false discovery rate than the label 'FDR < 0.05' implies; anything drawn from it needs orthogonal validation (ChIP-seq, qPCR, or replication in a new cohort) before being treated as confirmed.

Use a quasi-likelihood test and shrunken fold changes with small n

When: You have only 2-3 replicates per condition and need to avoid inflated false positives from unstable dispersion estimates.

In edgeR, prefer glmQLFit/glmQLFTest over the plain likelihood ratio test; in DESeq2, use lfcShrink() (apeglm) to get more reliable fold-change estimates for ranking and downstream motif work.

Caveat: These methods produce more conservative, more trustworthy estimates, but they do not create statistical power that more replicates would provide; a real small effect can still fail to clear FDR with n=2.

Reduce test count by merging redundant overlapping peaks

When: Your consensus set contains many fragmented or highly overlapping peak calls that inflate m without adding independent information.

Merge peaks within a small distance (a few hundred bp, matched to your fragment size distribution) into single consensus regions before counting, so each biological regulatory element is tested once.

Caveat: Overly aggressive merging can collapse genuinely distinct regulatory elements (e.g., adjacent enhancer and promoter) into one region, losing resolution you may need for interpretation.

When not to "fix" it

If peaks were called on a pooled, condition-blind BAM, filtering was done on overall abundance rather than per-condition counts, and the raw p-value histogram is genuinely flat, zero significant peaks at FDR < 0.05 is the correct answer, not a failure. Accessibility differences between two well-matched cell populations, especially with only 2-3 replicates, are often real but small, and BH correction across 100,000+ peaks is supposed to be conservative. Lowering the threshold, re-picking a smaller peak set to shrink m, or switching from padj to raw p-values until something turns significant is p-hacking dressed up as troubleshooting, and it will not replicate.

Five things experienced analysts do here

Plot the raw p-value histogram before you ever look at padj; a flat histogram tells you immediately that lowering the FDR cutoff won't help because there's no real signal to recover.
Always call ATAC-seq peaks on one pooled, equal-contribution BAM across all conditions, never per-sample then merged, so the peak set is fixed before you look at group differences.
Report log2FoldChange next to padj in every results table you hand off; a biologist evaluating a hit list needs the effect size to judge relevance, not just a pass/fail flag.
Never compare 'number of significant peaks' across two ATAC-seq experiments without checking they tested a similar number of total peaks first, since BH correction is calibrated to m.
Decide your filtering rule (which peaks to test, what FDR cutoff) before looking at the differential results, and write it down; deciding it after seeing the result is circular and inflates your real false discovery rate.

Questions people ask

Why do I have zero significant peaks after FDR correction in ATAC-seq?: Most often because you tested 100,000+ peaks instead of ~20,000 genes, so the Benjamini-Hochberg threshold is far stricter than what you're used to from RNA-seq. Check the raw p-value histogram first: if it's flat, there's genuinely little signal at your current power; if it spikes near 0 but padj is still all NA or high, look at independent filtering and how many peaks you actually tested.
What's the difference between padj and pvalue in DESeq2 for ATAC-seq?: pvalue is the per-peak test result before any correction; padj (the BH-adjusted value, DESeq2's estimate of the FDR) accounts for the total number of peaks tested and for independent filtering. Always report padj for significance calls, but look at pvalue and log2FoldChange together to judge whether a near-miss peak is worth following up.
Should I use Benjamini-Hochberg or Bonferroni for ATAC-seq peaks?: Use Benjamini-Hochberg (or DESeq2/edgeR's built-in BH adjustment). Bonferroni controls the family-wise error rate, which is appropriate for something like GWAS where a single false positive is very costly, but it is far too conservative for exploratory differential accessibility work across 100,000+ peaks.
What is independent filtering in DESeq2 and why does it matter for ATAC-seq?: Independent filtering removes low-count peaks before FDR adjustment using a criterion (mean normalized count) that is independent of the test statistic, which increases power without inflating false positives. In ATAC-seq it matters more than in RNA-seq because a much larger fraction of called peaks are low-signal MACS2 calls that would otherwise dilute the BH correction for every real peak.
Is FDR < 0.05 the right cutoff for ATAC-seq differential peaks?: It's a reasonable default, not a law. With only 2-3 replicates per condition and modest accessibility changes, FDR < 0.05 combined with a fold-change filter can leave you with nothing to work with even when real biology is present. Decide your cutoff before looking at the results, and treat effect size (log2FC) as at least as important as the adjusted p-value.

Guide · How to Handle Multiple Testing and FDR in Bulk RNA-seq
Compare · DESeq2 vs edgeR: Which One Should You Use?
Guide · How to Avoid Pseudoreplication in Single-Cell RNA-seq
Guide · How to Detect Batch Effects in Bulk RNA-seq
Guide · How to Handle Outlier Samples in PCA in Bulk RNA-seq
Glossary · Count matrix
Glossary · Effect size
Glossary · False discovery rate (FDR)

Sources

Understanding p value, multiple comparisons, FDR and q value — Definition of FDR and q-value used throughout the checks section
Filtering for ATAC-seq (Bioconductor support thread) — Peak calling must be blind to condition; per-sample filtering introduces circularity
Differential ATAC-seq and ChIP-seq peak detection using ROTS — edgeR/DESeq2 fold-change and FDR thresholds used in practice for ATAC-seq
Comparison of differential accessibility analysis strategies for ATAC-seq data — Published comparison showing FDR behavior differs across DA methods

Part of the Multiple testing and FDR series.