What should a good p-value histogram look like for differential ATAC-seq?

Flat across 0 to 1 with a spike in the first bin or two. The flat part is the peaks that are not changing, whose p-values are uniform under the null. The spike is the peaks that truly differ between conditions.

Why is my DESeq2 p-value histogram U-shaped?

The most common cause in ATAC-seq is a pile of low-count consensus peaks whose discrete p-values collect near one. Filter before testing, for example `keep = 5) >= ceiling(ncol(dds)*0.1)`, and replot. If the right-hand pile remains, look at your contrast and design.

Is a hill-shaped p-value histogram a problem?

Yes, it costs you power. It means variance is overestimated for the null peaks, so p-values run higher than they should. Check whether within-group variance differs between conditions, and whether an unmodeled batch is inflating dispersion.

Can I filter peaks after running DESeq2 or edgeR?

Not by anything correlated with the test statistic. The csaw book says filtering is valid only when it is independent of the test statistic under the null, for example filtering by overall mean abundance. Filtering on p-values or fold changes after testing breaks error control.

Does a flat histogram with no spike mean my ATAC-seq experiment failed?

Not necessarily. It means the test found no detectable signal, which can be real if the conditions barely differ. Check fragment-size distribution, TSS enrichment and fraction of reads in peaks before you decide the library is the problem.

Sanity check · ATAC-seq

How to Read a P-Value Histogram in ATAC-seq

Your differential accessibility table is only as good as the shape of its p-values, and one hist() call tells you whether to trust it before you look at a single peak.

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

You counted reads in a consensus peak set, ran DESeq2 or edgeR, and got a table of peaks with adjusted p-values. Before you pick a cutoff and send the list to motif enrichment, plot the raw p-values. It takes one line, and it tells you whether the test was calibrated.

What is at stake: a miscalibrated test gives you a peak list that looks like biology and is not. A U shape means low-count peaks are polluting the result. A hill means the model overestimates variance and you are losing real signal. A pile-up near one means something was filtered after the test. Each has a different fix, and applying the wrong one makes things worse.

This page shows you how to read the five shapes you will meet, which checks to run in what order, and how to filter consensus peaks correctly. You can finish all of it in the next hour on the results object you already have.

What it looks like when it's happening

`hist(res$pvalue)` shows a tall bar at 0 and a second tall bar at 1, with a trough in the middle (a U shape).
The histogram is a hill: low near 0, highest in the middle, falling toward 1, and almost no peaks pass FDR despite clear PCA separation.
The histogram is flat with no spike at 0, and the adjusted p-value column is nearly all 1 or NA.
The bar nearest 1 towers over everything else, while the rest of the histogram is flat or sloped.
Thousands of peaks pass FDR 0.05, yet the PCA shows samples separating by batch or library prep date rather than by condition.
The p-value histogram has visible gaps or stripes in specific bins, because many peaks have tiny counts and discrete test statistics.
`results()` returns many NA p-values and you cannot tell whether that is independent filtering or zero counts.

Why it happens

A p-value is uniform on [0,1] only when the null is true and the model is right. In a differential accessibility test, most peaks do not change between your conditions, so they contribute a flat floor. Peaks that truly change pile up near zero. That is the healthy picture: flat floor, spike at zero. Any other shape means some assumption failed for a group of peaks.

The ATAC-seq-specific problem is the consensus peak set. Merging peaks across all libraries gives you many regions that are open in one sample or one rare population and barely covered elsewhere. Their counts are small integers, so the test statistic is discrete and the p-values cannot spread evenly. They collect near one because they cannot provide evidence against the null. That is the U shape, and filtering low-count features before testing removes it.

A hill is a different failure. It means the null peaks get p-values higher than uniform, which happens when variance is overestimated. DESeq2 estimates one dispersion per peak across treatments, so if variance differs between your groups, or if an unmodeled batch adds noise inside groups, the dispersion is inflated and real differences are hidden. In ATAC-seq, batch is not a small worry: Tn5 is a single-turnover enzyme, so the nuclei-to-Tn5 ratio varies between preps and shifts fragment yield per sample.

A pile-up near one with an otherwise flat histogram points to filtering after testing, or a test that is conservative in its null. And one caveat: after you filter peaks sensibly, the histogram may still not be perfectly uniform even when the null is true, because peak calling itself enriches for sites that differ between conditions or vary within conditions. Read the shape as a diagnostic, not a pass/fail test.

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

01Plot the raw p-value histogram
Take the unadjusted p-values from the results object, drop NAs, and plot with enough bins to see the edges. Use the full result table, not the subset that passed a cutoff. For edgeR, get the full table with topTags(..., n = Inf, sort.by = 'PValue') and plot its PValue column.
r
```
dds <- DESeq(dds)
res <- results(dds)
hist(res$pvalue, col = 'grey', border = 'white', xlab = '', ylab = '', main = 'frequencies of p-values')

# ggplot2 version
res.df <- as.data.frame(res)
ggplot(res.df[!is.na(res.df$pvalue), ], aes(x = pvalue)) +
  geom_histogram(alpha = .5, position = 'identity', bins = 50) +
  labs(title = 'Histogram of unadjusted p-values') +
  xlab('Unadjusted p-values') + xlim(c(0, 1.0005))
```
Healthy
A roughly flat floor across the bins with a spike in the first bin or two. The height of the spike relative to the floor tells you how many peaks changed.
Red flag
A U shape, a hill, a bar at 1 that towers over the floor, or a histogram with no spike at all. Name the shape before you do anything else.
02Count NA p-values and where they come from
Count how many peaks have NA p-values. In DESeq2, NA means zero counts, an outlier flagged by Cook's distance, or removal by independent filtering. Sort out which, because independent filtering NAs are by design and the others are not.
r
```
sum(is.na(res$pvalue))
sum(is.na(res$padj) & !is.na(res$pvalue))   # removed by independent filtering
sum(rowSums(counts(dds)) == 0)               # all-zero peaks
```
Healthy
A modest set of NAs explained by zero-count peaks and independent filtering.
Red flag
A large fraction of peaks with NA p-values and no clear reason, which suggests the consensus set is dominated by barely covered regions.
03Look at the count distribution of the peaks piling up near 1
Split the p-values at a high threshold such as 0.9 and compare the mean normalized counts of peaks above and below it. Plot baseMean on a log scale for the two groups. If the peaks near 1 are the low-count ones, you have the U shape from discrete p-values.
Healthy
Peaks near 1 have a baseMean distribution similar to the rest of the table.
Red flag
Peaks with p-values near 1 are concentrated at the lowest baseMean values. That is the low-count signature, and it points to filtering before testing.
04Apply a pre-filter and replot
Filter out low-count consensus peaks before running DESeq(), rebuild the object, and plot again. The filter in the Bioconductor support thread keeps features with at least 5 counts in at least 10% of samples. Treat that as a starting point and check what it does to your peak count and to the histogram.
r
```
keep <- rowSums(counts(dds) >= 5) >= ceiling(ncol(dds)*0.1)
dds <- dds[keep, ]
dds <- DESeq(dds); res <- results(dds); hist(res$pvalue)
```
Healthy
The pile-up at 1 shrinks or disappears, and the floor flattens.
Red flag
The pile-up at 1 is unchanged after filtering. Then low counts are not the explanation, and you should check design and contrast. With few samples, 10% of samples rounds up to a single library, which may be too lenient: look at how many peaks you keep.
05Check PCA and sample metadata for batch
Run PCA on variance-stabilized counts and color points by condition, then by prep date, sequencing lane, operator and Tn5 batch. Compare how much of PC1 and PC2 is explained by each variable. If a technical variable lines up with condition, the batch effect overlaps the experimental groups.
Healthy
Replicates cluster by condition, or technical variables are spread evenly across conditions.
Red flag
Samples separate by a technical variable, especially one confounded with condition. This can give a hill, an inflated spike, or both.
06Compare within-group variability between conditions
Compute per-peak variance or look at the dispersion plot separately for each condition, using normalized counts. DESeq2 fits one dispersion per peak across treatments, so a condition with much higher variance can lift the null p-values and make a hill.
Healthy
Similar spread within each condition.
Red flag
One condition has visibly more spread, for example because one group contains a poor library. Check the fragment-size histogram, TSS enrichment and fraction of reads in peaks for that library.
07Confirm the design and contrast test what you think
Print resultsNames(dds) and the design formula, and confirm the contrast matches the comparison you intended. A composite null, such as testing a contrast that spans several groups, can produce a U shape even on well-filtered data. Rerun with the simplest two-group comparison and see if the shape changes.
Healthy
The shape is stable between the full design and a simple pairwise comparison.
Red flag
The shape changes sharply between the two. The contrast or the design is the problem, not the peaks.
08Audit the order of filtering and testing
Read your script from top to bottom and confirm that nothing removes peaks after testing on the basis of p-value, fold change or significance. Any filter must use a statistic independent of the test under the null, such as overall mean abundance or average log-CPM. Per the csaw book, filtering belongs at the peak-calling stage by pooling libraries with equal contributions, not after count-based testing.
Healthy
All filtering happens before the model is fit, using overall abundance.
Red flag
Peaks removed by condition-specific counts, p-value or fold change after testing. That gives a conservative pile-up near one and breaks FDR control.

What to do about it

Filter low-count consensus peaks before testing

When: The histogram has a U shape and the peaks near 1 are the low-baseMean ones.

Remove peaks with too little coverage before fitting the model. One documented starting point: keep peaks with at least 5 counts in at least 10% of samples, using rowSums(counts(dds) >= 5) >= ceiling(ncol(dds)*0.1). Pick the threshold on overall abundance only, never on condition labels or results, then refit and replot.

Caveat: With 2 to 4 replicates per condition, a 10% sample rule may keep a peak supported by one library. Check how many peaks survive and consider requiring support in at least the size of your smallest group.

Model the batch in the design

When: PCA or metadata shows a technical variable that is not fully confounded with condition, and the histogram is a hill or over-inflated.

Add the batch term to the design, for example ~ batch + condition, and rerun. Do not pre-correct the matrix and then test on the corrected values: put the term in the model so degrees of freedom are accounted for.

Caveat: If batch is perfectly confounded with condition, no design can separate them. You need to say so and treat the result as unreliable.

Switch to edgeR quasi-likelihood for small replicate numbers

When: You have few replicates and a complex design, and the DESeq2 histogram stays hill-shaped after batch is modeled.

Fit with edgeR using TMM normalization and the quasi-likelihood F-test, then extract the full table and replot the PValue column. Compare the shape to the DESeq2 histogram on the same filtered peak set.

Caveat: Both tools assume a negative binomial model. If the problem is a bad library or confounded batch, switching tools just gives you a second wrong answer.

Move filtering to the peak-calling stage

When: You filter by sample-specific signal or after testing, or the histogram shows a conservative pile-up near one.

Follow the csaw approach: define the peak set from pooled libraries with equal contributions, filter by overall average abundance before testing, and do no further p-value based removal afterward. Refit from the filtered counts.

Caveat: The csaw book covers ChIP-seq more than ATAC-seq, and guidance on filtering order for ATAC-seq is inconsistent across sources. Document your choice.

Remove or rerun a failed library

When: One condition has much higher within-group variance and QC on a library (fragment sizes, TSS enrichment, fraction of reads in peaks) is poor.

Check that library against the others. If it fails, drop it and state why, then refit. Do this on QC evidence, not because removing it improves the histogram.

Caveat: With 2 to 3 replicates, dropping one leaves very little. Sometimes the honest outcome is that the condition needs re-sequencing.

When not to "fix" it

Do not chase a flat histogram with no spike if your two conditions are close biologically. Two sorted populations that barely differ, or a mild treatment, can truly have almost no differential peaks, and the flat floor is the correct answer. Forcing a spike with looser filters, extra covariates or a different tool invents signal. Likewise, do not tune filters until the histogram looks perfect: filtered peak sets are not guaranteed uniform even under the null, because peak calling enriches for variable sites. Stop when the floor is roughly flat and the cause of any remaining shape is understood. And if a very large sample or a tiny p-value pile makes everything significant, the histogram is fine and the problem is effect size, so filter on log2 fold change.

Five things experienced analysts do here

Plot the p-value histogram before you look at the volcano plot or the top peaks. Once you have seen exciting hits, it is hard to believe the histogram.
Save the histogram before filtering and after filtering in the same report. The change between them documents what the filter did.
Decide the filter threshold from overall counts and write it down before you see results. A threshold tuned to the p-values is filtering after testing in disguise.
Plot baseMean against p-value or use a coloring by baseMean bins. It shows at once whether a shape comes from low-count peaks or from the whole table.
Read QC alongside the histogram: fragment-size distribution, TSS enrichment and fraction of reads in peaks. A weird histogram plus one weak library usually means the library, not the statistics.

Questions people ask

What should a good p-value histogram look like for differential ATAC-seq?: Flat across 0 to 1 with a spike in the first bin or two. The flat part is the peaks that are not changing, whose p-values are uniform under the null. The spike is the peaks that truly differ between conditions.
Why is my DESeq2 p-value histogram U-shaped?: The most common cause in ATAC-seq is a pile of low-count consensus peaks whose discrete p-values collect near one. Filter before testing, for example keep <- rowSums(counts(dds) >= 5) >= ceiling(ncol(dds)*0.1), and replot. If the right-hand pile remains, look at your contrast and design.
Is a hill-shaped p-value histogram a problem?: Yes, it costs you power. It means variance is overestimated for the null peaks, so p-values run higher than they should. Check whether within-group variance differs between conditions, and whether an unmodeled batch is inflating dispersion.
Can I filter peaks after running DESeq2 or edgeR?: Not by anything correlated with the test statistic. The csaw book says filtering is valid only when it is independent of the test statistic under the null, for example filtering by overall mean abundance. Filtering on p-values or fold changes after testing breaks error control.
Does a flat histogram with no spike mean my ATAC-seq experiment failed?: Not necessarily. It means the test found no detectable signal, which can be real if the conditions barely differ. Check fragment-size distribution, TSS enrichment and fraction of reads in peaks before you decide the library is the problem.

Guide · How to Choose a Normalization Method in ATAC-seq
Guide · How to Handle Multiple Testing and FDR in ATAC-seq
Guide · How to Read a P-Value Histogram in Bulk RNA-seq
Guide · How to Read a P-Value Histogram in ChIP-seq
Compare · DESeq2 vs edgeR: Which One Should You Use?
Glossary · ChIP-seq
Glossary · Principal component analysis (PCA)
Glossary · Effect size

Sources

Understanding p value, multiple comparisons, FDR and q value — Flat histogram with a spike near zero; conservative pile-up near one.
U-shape of p-value histograms is removed by filtering for genes that have >5 counts in >10% of samples before applying DESeq2 analysis — Low-count features cause the U shape; the filter code and threshold.
"U" shaped p values distribution in DESeq2 RNA-seq DE analysis using the contrast argument — Hill-shaped histograms and one dispersion per gene across treatments.
Chapter 4: Filtering out uninteresting windows (csaw book) — Independent filtering by average abundance; filtering must be independent of the test statistic.
Filtering for ATAC-seq — Filtered peak sets may be non-uniform under the null.
A generic reference defined by consensus peaks for single-cell ATAC-seq data analysis — Consensus peak counting yields many low-count peaks.
Reducing batch effects in single cell chromatin accessibility measurements by pooled transposition with MULTI-ATAC — Batch effects and Tn5 stoichiometry in ATAC-seq.
regionReport: regionReport package documentation — ggplot2 p-value histogram code for DESeq2 results.
edgeR: topTags function documentation — Extracting the full edgeR result table with p-value ranking.

Part of the P-value histograms series.