Why does FindMarkers give me a p-value of 1e-200 in scATAC-seq?

Because the default per-cell test treats every cell as an independent observation, when in reality cells from the same donor share that donor's chromatin state, tagmentation batch, and library prep. With thousands of cells but only 2-4 donors, the test's within-group variance collapses and even small or coincidental differences get scored with absurd confidence.

How many biological replicates do I need for scATAC-seq differential accessibility?

Treat it the same as bulk RNA-seq: 3 or more independent biological samples per condition is the practical floor for a DESeq2-style pseudobulk test. Cell count does not substitute for donor count; a study with 10,000 cells from 2 donors per condition still has an effective n of 2.

Does pseudobulk aggregation lose information compared to single-cell differential accessibility?

Yes, it collapses cell-to-cell heterogeneity within a sample into one summed count per peak. In exchange, published comparisons show pseudobulk gives better statistical power and tighter type I error control than naive single-cell tests, and it runs much faster than mixed models.

Can I use a mixed model instead of pseudobulk for scATAC-seq?

Yes. Generalized estimating equations or GLMMs that include donor as a random effect model the within-sample correlation explicitly and keep single-cell resolution, which matters if you need to test interactions between cell state and condition. They cost more compute and tend to show slightly worse type I error control than pseudobulk in head-to-head comparisons.

Is pseudoreplication in scATAC-seq worse than in scRNA-seq?

It compounds rather than replaces the scRNA-seq problem. The peak-by-cell matrix is sparser and more binary than a gene-by-cell matrix, so each cell carries less true information, and if peaks were called on pooled data instead of per cluster, rare cell types are underrepresented in the feature set feeding the test in the first place.

Sanity check · Single-Cell ATAC-seq

How to Avoid Pseudoreplication in Single-Cell ATAC-seq

Five thousand cells from two donors are not five thousand replicates, and your p-values are lying to you about it.

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

You extended an scRNA-seq differential accessibility script to your 10x Multiome or scATAC-seq data. Cluster the cells, run FindMarkers or getMarkerFeatures between two conditions, and the peak list comes back with p-values like 3e-180 across hundreds of peaks. It looks like a slam dunk result worth writing up immediately.

Those p-values are almost entirely an artifact of treating each of your cells as an independent observation. Your real n is the number of donors, mice, or biological samples, often 2 to 4, not the thousands of cells you sequenced from them. A peak "significant" at 1e-180 under this framing frequently fails to replicate when a third donor is added, and a collaborator or reviewer can burn months chasing a peak that was never real to begin with.

This page walks through how to tell whether your differential accessibility result is pseudoreplicated, how to rebuild it correctly as a pseudobulk test in DESeq2 or a mixed model that respects the sample as the unit of replication, and why the sparsity of the scATAC-seq peak-by-cell matrix makes this worse, not better, than the scRNA-seq case you may already be used to.

What it looks like when it's happening

FindMarkers or getMarkerFeatures returns p-values in the 1e-100 to 1e-300 range for peaks with a log2FC under 0.25.
The list of 'significant' peaks barely changes whether you filter at adjusted p < 0.05 or p < 1e-50; the p-value distribution is saturated near zero.
Your design has only 2 to 4 biological donors or mice per condition but thousands of cells, and the DA test's group.by, ident.1, or ident.2 argument was set at the cell or cluster level with no sample term anywhere.
A peak called significant with donors A and B vs donors C and D flips direction or drops out entirely when you add a third donor per condition.
The volcano plot shows a wall of points hugging the p ≈ 0 axis paired with only modest fold changes.
Pseudobulk fragment totals per sample vary by an order of magnitude because minCells/maxCells constraints were never applied, so one or two large donors dominate every consensus estimate.

Why it happens

Standard statistical tests, whether it is the Wilcoxon rank-sum or likelihood-ratio test inside FindMarkers, or a negative binomial GLM applied per cell, assume each row is an independent draw from the population you want to generalize to. Cells from the same donor are not independent draws: they share that donor's genotype, disease stage, cell-state composition, tagmentation batch, and library prep. Testing 5,000 cells from 2 donors as though they were 5,000 independent replicates is really testing n = 2 with 2,500-fold pseudoreplication. The test's within-group variance collapses toward zero because cell-to-cell noise within one donor is far smaller than the true biological variance between donors, so almost any separation between conditions, real or coincidental, gets scored with absurd confidence.

scATAC-seq compounds this rather than escaping it. The peak-by-cell matrix carries only a few thousand Tn5 insertion events per cell, spread thin and near-binary across tens of thousands of peaks, so each cell carries less true information than the UMI counts in an scRNA-seq gene-by-cell matrix. The apparent "n" you get from counting cells is even more inflated relative to the actual information content behind it. On top of that, if peaks were called on the pooled dataset instead of per cluster, rare cell types are systematically underrepresented in the feature set feeding the test, which can manufacture peak-level differences that track cell-type composition shifts rather than a true accessibility change within a shared population.

The biological cause underneath all of this: donors differ in genotype, age, disease stage, cell-state mix, and how their samples were handled at the bench, from tagmentation batch to fragment size distribution. Any of these donor-level differences looks exactly like a condition effect when you only have 2 donors per condition and 5,000 cells standing in for them. No amount of sequencing depth or cell count lets you tell "this donor's biology" apart from "the condition" when your donor count sits at n = 2.

The checks

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

0/6 checked · saved in this browser

01Count biological replicates, not cells
Before looking at any p-value, tabulate your metadata by sample and condition. Count unique donors/mice/patients per condition, not the number of cells contributed by each.
r
```
table(atac$sample, atac$condition)
```
Healthy
3 or more independent biological samples per condition. This is the same bar DESeq2 documentation uses for bulk RNA-seq, and scATAC-seq is not exempt just because you have thousands of cells.
Red flag
2 samples per condition, or samples that are technical replicates of the same donor (re-tagmented, re-sequenced) rather than distinct biological units.
02Inspect the actual test call for a sample term
Grep your analysis script for the differential accessibility call (Signac FindMarkers or ArchR getMarkerFeatures) and check whether sample/donor appears anywhere in the arguments.
bash
```
grep -n "FindMarkers\|getMarkerFeatures" analysis.R
```
Healthy
The call either operates on a pre-aggregated pseudobulk matrix, or passes a sample/donor grouping into a batch-aware argument.
Red flag
group.by is set to cluster or condition only, ident.1/ident.2 select individual cells directly, and 'sample' or 'donor' never appears in the call.
03Check p-value magnitude against effect size
Plot a histogram of adjusted p-values and a volcano plot of log2FC against -log10(p) for the differential accessibility results.
r
```
hist(da_peaks$p_val_adj, breaks = 50)
plot(da_peaks$avg_log2FC, -log10(da_peaks$p_val_adj))
```
Healthy
P-values spread across a range, and the strongest significance corresponds to the largest fold changes.
Red flag
A wall of peaks sitting at p = 1e-100 to 1e-300 with log2FC under 0.25-0.5. Hundreds of peaks 'significant' regardless of whether you use p < 0.05 or p < 1e-50 is a saturation artifact, not a strong signal.
04Rebuild the test as pseudobulk and compare
Aggregate raw fragment counts per peak, summed within each sample x cell-type group, then run DESeq2 on the resulting sample-level matrix instead of the per-cell matrix.
r
```
pb <- aggregate.Matrix(t(counts(sce)), groupings = paste(sce$sample, sce$celltype), fun = "sum")
dds <- DESeqDataSetFromMatrix(countData = pb, colData = metadata, design = ~ condition)
dds <- DESeq(dds)
res <- results(dds)
```
Healthy
The significant peak count drops substantially, often by one to two orders of magnitude, and the surviving hits show larger, more biologically plausible fold changes.
Red flag
The pseudobulk result returns nearly the same hit list and the same implausibly small p-values as the per-cell test. That usually means the aggregation grouped by cluster only and never collapsed to one row per sample.
05Leave-one-donor-out sensitivity check
Rerun the pseudobulk DESeq2 test N times, each time excluding one biological sample, and compare which peaks stay significant across all runs.
r
```
for (s in unique(metadata$sample)) {
  keep <- metadata$sample != s
  dds_sub <- DESeqDataSetFromMatrix(pb[, keep], metadata[keep, ], design = ~ condition)
  dds_sub <- DESeq(dds_sub)
  # record which peaks are significant
}
```
Healthy
A core set of peaks stays significant regardless of which single donor is dropped.
Red flag
Significance depends heavily on including one specific donor. The result is being driven by an outlier sample, not the condition.
06Check zero-inflation at the pseudobulk level
Even after aggregation, scATAC-seq pseudobulk counts carry excess zeros. Check the fraction of zero counts per peak across samples and whether DESeq2's negative binomial fit looks reasonable.
r
```
zero_frac <- rowMeans(pb == 0)
summary(zero_frac)
plotDispEsts(dds)
```
Healthy
A moderate, roughly consistent zero fraction across peaks, with dispersion estimates that shrink smoothly toward the fitted trend.
Red flag
Many peaks with zero counts in most samples but a large outlier count in one sample, or dispersion estimates that scatter wildly. This is a sign plain DESeq2's NB model is not enough and a zero-inflated method (e.g., scaDA) is worth trying.

What to do about it

Pseudobulk aggregation + DESeq2 or edgeR

When: You have 3 or more biological samples per condition and want the standard, well-supported answer.

Sum raw fragment counts per peak within each sample x cell-type group with aggregate.Matrix(), build a DESeqDataSet with sample-level metadata (design = ~ condition), run DESeq(), and extract shrunken fold changes with lfcShrink().

Caveat: You lose cell-to-cell resolution and cannot test interactions between cell state and condition within a sample. If you need that, use a mixed model instead.

ArchR pseudo-bulk replicates

When: You are already in an ArchR workflow with per-cluster peaks.

Call addGroupCoverages(ArchRProj, groupBy = "CellType", minCells = 40, maxCells = 500) to build sample-aware pseudo-bulk replicates, then run getMarkerFeatures, which uses these replicates rather than raw cells for the test.

Caveat: minCells and maxCells are data-dependent and there is no universal threshold; too low and you get noisy pseudo-replicates, too high and small samples get dropped entirely.

Mixed model or GEE with sample as a random effect

When: You need to keep single-cell resolution, or you want to test whether a condition effect differs across cell states.

Fit a GLMM or generalized estimating equation with donor/sample as a random effect and a Poisson or negative binomial link, modeling the within-sample cell-to-cell correlation directly instead of aggregating it away.

Caveat: Heavier compute, more convergence issues on sparse scATAC-seq data, and published comparisons show somewhat worse type I error control than pseudobulk-DESeq2 in most tested scenarios.

Zero-inflated model at the pseudobulk level (scaDA)

When: Your pseudobulk counts still show excess zeros beyond what a negative binomial model expects, i.e. peaks accessible in only a handful of samples.

Use a method like scaDA that jointly models mean, zero-inflation prevalence, and dispersion with empirical Bayes shrinkage, built specifically for the sparsity and overdispersion of scATAC-seq pseudobulk data.

Caveat: Newer and less battle-tested than DESeq2, with a smaller user community and fewer worked examples to debug against.

Add biological replicates instead of more cells

When: Your pilot dataset only has 2 donors per condition and the differential accessibility result matters for a downstream decision.

Sequence additional donors/animals before drawing a conclusion. Statistical power against a condition effect comes from the number of independent biological units, not the number of cells sampled from each.

Caveat: Costs time and money you may not have. If you can't add samples, report the result as hypothesis-generating and say explicitly that it has not been replicated across donors.

When not to "fix" it

If your goal is to describe within-sample heterogeneity, not to make a between-condition or between-donor claim, per-cell testing is the right tool and pseudobulk correction would throw away exactly the resolution you need. Characterizing which cell states differ within a single tumor's chromatin landscape, or ranking peaks by accessibility across cell types within one Multiome dataset, does not require a population-level p-value defended against pseudoreplication. Also skip the fix when your "replicates" really are intentional technical replicates from the same biological sample, profiled to estimate technical variance rather than to generalize to a population; just label the claim as technical, not biological, and do not present its p-value as evidence of a donor-level effect.

Five things experienced analysts do here

Ask 'what is my real n' before trusting any p-value: count donors, not cells, and write that number at the top of your analysis notebook.
Keep raw, unnormalized fragment counts around through clustering and cell typing so you can aggregate.Matrix() them into a pseudobulk matrix later; DESeq2 needs integer counts, not TF-IDF-transformed values.
Run the naive per-cell test and the pseudobulk test side by side and compare hit lists. If the pseudobulk result shrinks by one or two orders of magnitude, that gap is the pseudoreplication you almost got fooled by.
Call peaks per cluster or cell type before you aggregate, not on the pooled data. A peak set biased toward abundant cell types will bias every pseudobulk count matrix built on top of it.
Budget donors the way you would replicates in bulk RNA-seq: 3 minimum, more if the expected effect size is subtle. No amount of extra cells per donor buys you power at the donor level.

Questions people ask

Why does FindMarkers give me a p-value of 1e-200 in scATAC-seq?: Because the default per-cell test treats every cell as an independent observation, when in reality cells from the same donor share that donor's chromatin state, tagmentation batch, and library prep. With thousands of cells but only 2-4 donors, the test's within-group variance collapses and even small or coincidental differences get scored with absurd confidence.
How many biological replicates do I need for scATAC-seq differential accessibility?: Treat it the same as bulk RNA-seq: 3 or more independent biological samples per condition is the practical floor for a DESeq2-style pseudobulk test. Cell count does not substitute for donor count; a study with 10,000 cells from 2 donors per condition still has an effective n of 2.
Does pseudobulk aggregation lose information compared to single-cell differential accessibility?: Yes, it collapses cell-to-cell heterogeneity within a sample into one summed count per peak. In exchange, published comparisons show pseudobulk gives better statistical power and tighter type I error control than naive single-cell tests, and it runs much faster than mixed models.
Can I use a mixed model instead of pseudobulk for scATAC-seq?: Yes. Generalized estimating equations or GLMMs that include donor as a random effect model the within-sample correlation explicitly and keep single-cell resolution, which matters if you need to test interactions between cell state and condition. They cost more compute and tend to show slightly worse type I error control than pseudobulk in head-to-head comparisons.
Is pseudoreplication in scATAC-seq worse than in scRNA-seq?: It compounds rather than replaces the scRNA-seq problem. The peak-by-cell matrix is sparser and more binary than a gene-by-cell matrix, so each cell carries less true information, and if peaks were called on pooled data instead of per cluster, rare cell types are underrepresented in the feature set feeding the test in the first place.

Guide · How to Call Peaks You Can Trust in Single-Cell ATAC-seq
Guide · How to Detect Batch Effects in Single-Cell ATAC-seq
Guide · How to Tell If You Overclustered in Single-Cell ATAC-seq
Guide · How to Avoid Pseudoreplication in Single-Cell RNA-seq
Guide · How to Avoid Pseudoreplication in Single-Nucleus RNA-seq
Glossary · Effect size
Glossary · P-value
Glossary · Peak calling

Sources

How to create pseudobulk from single-cell RNAseq data — Foundational pseudobulk aggregation workflow this page adapts for scATAC-seq
Strategies for addressing pseudoreplication in multi-patient scRNA-seq data (bioRxiv 2024) — Comparison of pseudobulk vs mixed models vs naive single-cell tests, including scATAC-seq applications and zero-inflation findings
scaDA: A Novel Statistical Method for Differential Analysis of Single-Cell Chromatin Accessibility Sequencing Data — Zero-inflated negative binomial method for pseudobulk-level scATAC-seq differential accessibility
Pseudobulk differential expression analysis with DESeq2 (HBC Training) — Step-by-step aggregate.Matrix() and DESeq2 workflow used in the checks and fixes

Part of the Pseudoreplication series.