Should I use DESeq2 median-of-ratios or edgeR TMM for ATAC-seq?

For a standard design with most peaks unchanged, either works. Both estimate scale factors from the unchanged majority and fail the same way under global shifts. Choose one, run it on raw consensus-peak counts, and confirm with an MA plot. edgeR's quasi-likelihood framework is the usual pick for complex designs with few replicates.

Can I use CPM or TPM for differential accessibility?

No. CPM corrects for depth only and is meant for comparing values across samples, and TPM adds a length term designed for transcripts. Count-based tests like DESeq2 and edgeR need raw counts and apply their own factors. Use CPM for tracks and display.

Why do my CPM tracks differ from the DESeq2 results?

The scale factors come from different places. CPM uses total reads, or reads in peaks if you scale that way, and when FRiP differs between samples these give different tracks. DESeq2 uses the unchanged majority of peaks. Compare size factors across methods to see which one is moved by library quality.

Does SCTransform or LogNormalize work for ATAC-seq?

The sources I have do not cover SCTransform on ATAC peak matrices, so I would not assume it transfers. Single-cell ATAC workflows commonly use TF-IDF followed by LSI, and drop the first component when it tracks sequencing depth. Validate any choice by checking that depth does not drive your embedding.

What if I expect a global change in accessibility?

Then median-of-ratios and TMM are the wrong default. Use a non-linear method such as loess, or scale with spike-in, and check the MA plot afterward. Confirm the result against an independent assay such as RNA-seq on the same samples.

Sanity check · ATAC-seq

How to Choose a Normalization Method in ATAC-seq

The method you pick can turn 33 significant regions into 24,450 on the same data, so decide it on purpose and check it with an MA plot.

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

You have a consensus peak set, a count matrix, and a bigWig per sample. Someone asks for a heatmap, a differential accessibility table and a motif result, and each of those can take a different normalization. You pick whatever the last tutorial used, and it works until a reviewer asks why the treated sample has no signal anywhere.

The stakes are large. In a published comparison of ATAC-seq normalization strategies, the number of significant regions at FDR < 0.10 ranged from 33 to 24,450 depending on the method alone. Same reads, same peaks, same contrast. If you do not know which assumption your method makes, you cannot tell whether your hits are biology.

This page gives you a decision you can make in the next hour: which normalization fits a track, a heatmap or a statistical test, which ones must never feed each other, and the diagnostic plot that tells you whether the assumption held in your data.

What it looks like when it's happening

An MA plot of treated vs control peaks is tilted: the cloud of unchanged peaks sits above or below M = 0 instead of centered on it.
Almost every significant peak goes in one direction, for example 90% or more of hits are 'gained' in one condition.
Two normalization methods run on the same count matrix return hit lists with little overlap, or hit counts that differ by orders of magnitude.
Two samples with different FRiP look different in CPM bigWigs in the same loci where the count-based test says nothing changed.
The aggregate TSS profile from computeMatrix shows one condition flat at every site, or shows a gain at all sites together.
A heatmap clusters by sequencing depth or library quality instead of by condition.
PCA PC1 tracks total reads or FRiP, not the biology you sorted the cells for.

Why it happens

Every library-size method turns raw counts into a comparison by dividing each sample by a scale factor. The only question is what the scale factor is computed from. CPM uses total reads. Median-of-ratios (DESeq2) and TMM (edgeR) use the bulk of regions that are assumed not to change, so they estimate the factor from the unchanged majority. All of them rest on the same premise: most regions are equally accessible in both conditions. That is what lets one number per sample absorb depth differences.

ATAC-seq breaks that premise more often than RNA-seq does. A treatment, a differentiation step or a knockout of a chromatin remodeler can open or close chromatin broadly. A sorted population can also differ in composition, so a subpopulation with different accessibility shifts a large share of peaks. When the unchanged-majority assumption is false, median-of-ratios and TMM force the bulk to zero and push the real global change into noise or into the opposite direction. The MA plot shows it as a tilt or a skew.

The denominator also matters in a way RNA-seq users do not expect. Total reads include background reads outside peaks, and the share of reads in peaks (FRiP) varies by library. The packet's guidance is a minimum of 0.2, with 0.3 considered good. CPM computed from all reads and CPM computed from reads in peaks are two different scale factors. When FRiP differs between samples they give different tracks, and both mislead if accessibility changes globally. This is a known source of confusion in the deepTools tracker for fragment-filtered CPM tracks.

TPM and RPKM add a length term built for transcripts. Peak widths come from the caller or from your consensus step, not from biology, so dividing by them mostly injects an artifact. The lesson on RNA-seq normalization is blunt about CPM and TPM: they are for comparing values across samples or genes, not for differential testing. Feeding CPM or TPM into a count-based test, or normalizing a matrix twice, breaks the model, because DESeq2 and edgeR model raw counts with their own offsets.

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

01Decide the question before the method
Write one sentence for each output: a browser track, a heatmap or aggregate plot, a differential accessibility test. Then ask whether you expect most peaks to stay constant between conditions. If the biology predicts a global opening or closing, or a large shift in cell composition, flag the project as a global-shift design before any normalization runs.
bash
```
# Tn5 shift first, then normalize tracks. Shift BAMs before bamCoverage.
alignmentSieve --ATACshift --input input.bam --output output.bam
```
Healthy
You can state, per output, which normalization it uses and why. A mainstream design (sorted populations, a modest perturbation) is marked as 'most peaks unchanged'.
Red flag
You are using one normalized matrix for tracks, heatmaps and the statistical test, or you cannot say whether a global shift is plausible.
02Look at FRiP and fragment size per sample before comparing anything
Compute FRiP for every sample against the consensus peak set and plot it next to total reads. Check the fragment-size histogram and TSS enrichment too, since a failed library can look fine by read count. Mark any sample where FRiP differs strongly from its replicates or from the other condition.
Healthy
FRiP above 0.2 for every sample, with above 0.3 considered optimal, and similar values within a condition.
Red flag
FRiP below 0.2, or FRiP that differs systematically between conditions. Then CPM over all reads and CPM over peaks diverge, and a depth-based factor reflects library quality, not accessibility.
03Compare size factors from different methods
Compute total-read library sizes, DESeq2 size factors and edgeR TMM factors on the same raw peak count matrix. Plot each against the others as a scatter, colored by condition. Large disagreement means the composition of the data is driving the factors.
r
```
library(DESeq2); library(edgeR)
dds <- DESeqDataSetFromMatrix(countData = raw_counts, colData = coldata, design = ~ condition)
dds <- estimateSizeFactors(dds)
sf_rle <- sizeFactors(dds)

y <- DGEList(counts = raw_counts)
y <- calcNormFactors(y, method = "TMM")
sf_tmm <- y$samples$norm.factors * y$samples$lib.size

lib <- colSums(raw_counts)
pairs(log2(cbind(lib = lib, rle = sf_rle, tmm = sf_tmm)),
      col = as.integer(factor(coldata$condition)), pch = 19)
```
Healthy
The three sets of factors track each other closely and the points do not separate by condition.
Red flag
The DESeq2 or TMM factors split by condition while the library sizes do not. The methods are seeing a composition or global shift.
04MA plot after normalization
For each contrast, plot M (log2 fold change) against A (mean normalized log intensity) for all consensus peaks, after the normalization you plan to use. Draw a horizontal line at M = 0 and overlay a smoother. The published ATAC-seq comparison uses this plot to expose systematic bias and asymmetry.
r
```
res <- results(DESeq(dds))
plotMA(res, ylim = c(-3, 3))
abline(h = 0, col = "red")

# edgeR equivalent
plotMD(y, column = 1); abline(h = 0, col = "red")
```
Healthy
The bulk of peaks forms a cloud centered on M = 0 across the range of A, with the significant peaks scattered on both sides.
Red flag
The cloud is shifted off zero, curved along A, or one-sided. Either the method assumption failed or the data need a non-linear method such as loess.
05Check the direction balance of significant peaks
Count significant gained and lost peaks at your FDR cutoff. Then rerun with a second method (for example TMM if you used median-of-ratios, or loess) and compare both the counts and the overlap of the hit lists.
Healthy
Gains and losses are in a plausible ratio for the biology, and the two methods return overlapping lists with similar counts.
Red flag
Hit counts that differ by large factors between methods, or nearly all hits in one direction with no biological reason. That is the 33 versus 24,450 problem showing up in your data.
06Compare CPM tracks computed two ways
Make bigWigs with bamCoverage using --normalizeUsing CPM, and compute a second version where the scale factor comes from reads in peaks only (scale the track by 1e6 divided by peak reads). Open both at loci where the count-based test says nothing changed and at loci it calls differential. If you filter fragment lengths, remember the fragment filter affects what gets counted in the denominator, as the deepTools issue on mono-nucleosome CPM discusses.
bash
```
bamCoverage --normalizeUsing CPM -i input.bam -o output.cpm.bw
bamCoverage --normalizeUsing RPGC --effectiveGenomeSize <size> -i input.bam -o output.rpgc.bw
```
Healthy
Tracks agree at unchanged loci and agree on direction at the differential ones, with only modest scale differences.
Red flag
Tracks disagree at unchanged loci, or the all-reads and in-peak versions flip which condition is higher. FRiP differences are driving the display.
07Test for a global shift with an external anchor
If you have spike-in or a set of control regions you trust to be constant, compute factors from them and compare to the factors from step 3. Spike-in factors are calculated by dividing the spike-in count of the sample with the lowest spike-in counts by each sample's spike-in count, so the lowest sample gets 1. Also check an aggregate TSS profile from computeMatrix and plotHeatmap under the candidate normalization.
Healthy
Spike-in or control-region factors agree with the data-derived factors, and the aggregate profile differences look like modest changes, not a uniform offset.
Red flag
Spike-in factors disagree with TMM or median-of-ratios factors in the same direction as the unexplained tilt. The unchanged-majority assumption is wrong for this experiment.
08Orthogonal sanity check against RNA-seq
If the same samples have RNA-seq, take genes with clear expression changes and compare the accessibility change at their promoters or linked peaks under each normalization. Prefer the method whose output agrees with the independent assay across the board, and report the comparison.
Healthy
Direction of accessibility change at promoters of strongly changed genes agrees with expression in most cases, under the method you chose.
Red flag
Agreement appears only under one method, or the sign flips for most genes under your chosen method. Revisit the normalization before reporting biology.

What to do about it

Median-of-ratios or TMM on raw consensus-peak counts

When: Most peaks should be unchanged: sorted populations, a modest perturbation, and an MA plot centered on zero. This is the default for differential accessibility.

Count reads in the consensus peak set to get an integer matrix. For DESeq2 build the dataset from raw counts, call estimateSizeFactors, and test with DESeq. For edgeR use DGEList, calcNormFactors with TMM, and the quasi-likelihood F-test, which suits complex designs and small replicate numbers. Pick one of the two and do not feed one tool's normalized output into the other.

Caveat: Both fail when global accessibility shifts. The packet does not give a rule for choosing between RLE and TMM, so treat them as near-equivalents and verify with the MA plot.

Loess-based normalization for suspected global change

When: The MA plot is curved or tilted, or the biology predicts broad opening or closing of chromatin.

Use a non-linear method such as the csaw loess option, which the published comparison prefers when global shifts are expected. Rerun the MA plot after normalization and confirm the cloud is flat at M = 0 across A.

Caveat: Loess can remove a real global signal too, so confirm with an external anchor such as spike-in or RNA-seq before you claim 'no global change'.

Spike-in or control-based factors

When: You designed the experiment with spike-in, or you know a global shift is the question.

Count the spike-in reads per sample, set the lowest-spike-in sample to a factor of 1 and scale the others by dividing that lowest count by their own, as described in the RNA-seq normalization lesson. Apply the factors as size factors in your count model, not as a second normalization on top.

Caveat: Only as good as the spike-in addition. Uneven spike-in amounts become a technical batch effect that you then call biology.

Pick the right deepTools option for tracks and heatmaps

When: You need bigWigs for viewing or computeMatrix, not a statistical test.

Use bamCoverage --normalizeUsing CPM for library-size scaling, or RPGC with --effectiveGenomeSize for 1x depth normalization. Run alignmentSieve --ATACshift first. Keep bigWig scaling and the DE normalization separate in your methods text. Do not use RPKM or BPM for ATAC peaks, since the length term has no biological meaning for peak widths.

Caveat: CPM over all reads and over peaks differ when FRiP differs. Tracks are for display, so do not read a differential call off them.

Single-cell ATAC: TF-IDF and LSI

When: The data are scATAC, not bulk replicates.

Use TF-IDF followed by LSI, and inspect the first LSI component. It often captures sequencing depth, so drop it if it correlates with total counts.

Caveat: The packet has no source on SCTransform for ATAC. It is a single-cell RNA method, so do not assume it carries over to peak matrices.

When not to "fix" it

If the MA plot is centered and the factors agree, do not switch methods to chase more hits. Trying several normalizations until the hit list grows is how you land on 24,450 regions that mean nothing. Pick the method on assumptions before looking at hit counts, and report the others as sensitivity checks.

Also leave a real global shift alone. If a knockout closes chromatin broadly and spike-in or RNA-seq agree, a tilted MA plot is biology, and forcing it flat with median-of-ratios or loess erases the finding. In that case the right move is to keep the shift and explain it, using an external anchor for scaling.

Five things experienced analysts do here

Keep raw integer counts in the object forever. DESeq2 and edgeR need them, and every normalized matrix (CPM, TPM, size-factor scaled) should be derived and labeled, never saved over the original.
Name every matrix and bigWig by its normalization, such as `peaks_raw`, `peaks_rle`, `tracks_cpm_allreads`. Mixing a normalized matrix into a count-based test is the error that never throws a warning.
Run the MA plot for every contrast before reading any hit list, and save it with the results. A reviewer who sees the tilt gets the answer to your normalization choice in one glance.
Write down the global-shift question at the design stage, before the sequencing. If the answer is 'maybe', add spike-in then, because you cannot add it later.
Treat hit counts as a sensitivity result: run a second method and report the overlap. If two sensible methods disagree sharply, your data are telling you the assumption is shaky, not that one tool is buggy.

Questions people ask

Should I use DESeq2 median-of-ratios or edgeR TMM for ATAC-seq?: For a standard design with most peaks unchanged, either works. Both estimate scale factors from the unchanged majority and fail the same way under global shifts. Choose one, run it on raw consensus-peak counts, and confirm with an MA plot. edgeR's quasi-likelihood framework is the usual pick for complex designs with few replicates.
Can I use CPM or TPM for differential accessibility?: No. CPM corrects for depth only and is meant for comparing values across samples, and TPM adds a length term designed for transcripts. Count-based tests like DESeq2 and edgeR need raw counts and apply their own factors. Use CPM for tracks and display.
Why do my CPM tracks differ from the DESeq2 results?: The scale factors come from different places. CPM uses total reads, or reads in peaks if you scale that way, and when FRiP differs between samples these give different tracks. DESeq2 uses the unchanged majority of peaks. Compare size factors across methods to see which one is moved by library quality.
Does SCTransform or LogNormalize work for ATAC-seq?: The sources I have do not cover SCTransform on ATAC peak matrices, so I would not assume it transfers. Single-cell ATAC workflows commonly use TF-IDF followed by LSI, and drop the first component when it tracks sequencing depth. Validate any choice by checking that depth does not drive your embedding.
What if I expect a global change in accessibility?: Then median-of-ratios and TMM are the wrong default. Use a non-linear method such as loess, or scale with spike-in, and check the MA plot afterward. Confirm the result against an independent assay such as RNA-seq on the same samples.

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

Sources

ATAC-seq normalization method can significantly affect differential accessibility analysis and interpretation — Source for 33 to 24,450 significant regions across eight strategies, loess for global shifts, MA plot inspection, and orthogonal validation.
bamCoverage, deepTools Documentation — Source for the RPKM, CPM, BPM, RPGC and none options, the RPGC effective genome size requirement, and the ATACshift preprocessing step.
Normalizations, deepTools Wiki — Overview of deepTools normalization and GC bias correction with correctGCbias.
CPM normalized bigWig for mono-nucleosome fragments in ATAC-seq (GitHub Issue #841) — Shows how CPM depends on which reads define library size when FRiP differs.

Part of the Normalization choice series.