Chatomics Field GuideWhat They Don't Teach You

Sanity check · Bulk RNA-seq

How to Read a P-Value Histogram in Bulk RNA-seq

A misshapen histogram is the cheapest bug report your DE test will ever hand you, and almost nobody opens it before trusting padj.

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

You ran DESeq2 or edgeR, called results(), and jumped straight to summary(res) or a padj < 0.05 filter to count how many genes changed. That's the moment this page is for. Before you trust that gene list, look at one plot most people skip entirely: the raw p-value histogram.

The shape of that histogram is the cheapest diagnostic in your whole workflow, one line of code, and it tells you whether the statistics underneath your gene list are calibrated at all. A U shape, a hump in the middle, or a suspiciously flat line each point to a specific, fixable problem: low-count genes flooding the test, an unmodeled batch effect, or a dispersion fit that doesn't match your data. If you skip this check, you inherit whatever is wrong straight into your FDR correction, your gene list, your enrichment analysis, and your figure.

In the next hour: run hist(res$pvalue) on your own results, match what you see against the shapes below, and work through the ordered checks to find which of a handful of causes is actually responsible. Most of the time the fix lives in your design formula or your filtering step, not in a fancier statistical test.

What it looks like when it's happening

  • hist(res$pvalue) shows two spikes, one near 0 and one near 1, with a valley between them (U shape)
  • The histogram has a visible bulge or hump centered around 0.4 to 0.6 instead of being flat outside the zero spike
  • summary(res) reports almost no genes at padj < 0.05 despite an obvious biological perturbation like a knockout or drug treatment
  • Nearly every tested gene has a p-value below 1e-10, giving you thousands of 'significant' genes for a modest perturbation
  • The histogram is close to flat everywhere with no spike near zero at all
  • plotDispEsts(dds) shows many gene-wise dispersion points that refuse to shrink toward the fitted trend line
  • A PCA plot shows PC1 or PC2 separating samples by processing day, replicate number, or batch rather than by condition

Why it happens

Under a correctly specified null hypothesis, a gene's p-value is uniformly distributed on [0,1] by construction. Stack the null p-values from thousands of non-DE genes and you get a flat histogram; overlay the minority of truly differentially expressed genes and you get a spike near zero sitting on top of that flat baseline. Any other shape means one of the assumptions behind the test, the mean-variance relationship, the dispersion estimate, or the design formula, doesn't match your data.

The U shape shows up after multiple-testing correction when a large number of genes are barely above detection. A negative binomial fit for a gene with near-zero counts across most samples is unstable: there's essentially no signal relative to noise, so the test can't reject the null and the p-value piles up near 1 instead of spreading uniformly. Stack enough of those low-count genes and you get a second spike at 1, with a valley in between. This interacts with DESeq2's independent filtering, which is a separate mechanism from removing genes outright: independent filtering sets padj to NA for low-mean genes rather than deleting their raw p-values, so the U shape is visible in res$pvalue even when independent filtering is doing its job downstream.

The hill or hump shape comes from a different failure: an unmodeled batch effect or an overestimated dispersion inflates the residual variance the model attributes to "within-group noise." That miscalibration pushes p-values toward the middle of the range rather than to either tail; the test is neither confidently rejecting nor confidently failing to reject. Because the miscalibration lives in how the GLM estimates variance, the fix has to happen inside the model. Including batch in the design formula (~ batch + condition) lets DESeq2 estimate and subtract its contribution to variance correctly. Pre-correcting the count matrix with removeBatchEffect() before DESeq2 throws that information away before the GLM ever sees it, which is why it isn't the recommended fix.

Uneven sequencing depth that correlates with your biological variable produces a look-alike of the batch problem, not a separate mechanism. Size-factor and TMM normalization in DESeq2, edgeR, and limma-voom handle moderate depth differences adequately, but weakly expressed genes stay sensitive to depth-driven sampling noise. If depth happens to track condition, the resulting artifact shows up in the same PCA and histogram checks you'd use to catch a batch effect.

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

  1. Right after results(), before summary(), before a volcano plot, run a plain histogram of the unadjusted p-values with enough bins to see structure (50 is a reasonable default).

    r
    hist(res$pvalue, breaks = 50, col = "grey", main = "p-value distribution")
    Healthy
    A sharp spike in the leftmost bin (genes with real DE), roughly flat across the rest of [0,1], and for count data a small extra bump right at 1 is normal and not a red flag on its own.
    Red flag
    Two spikes with a valley between them (U shape), a bulge in the middle (hill shape), or a histogram that's flat everywhere with no spike near zero at all.
  2. Check the results() call that produced res: was lfcThreshold set to anything above 0?

    Healthy
    `lfcThreshold = 0` (the default) was used, so any deviation from a flat-plus-spike shape is a real diagnostic signal.
    Red flag
    `lfcThreshold` was set >0 and the histogram is U-shaped. That U shape is the expected output of DESeq2's composite null test in this mode, not a sign your data is broken. Stop here if this is your case.
  3. Print the design formula and compare it against every technical variable you recorded: extraction date, sequencing batch, RNA prep person, lane.

    r
    design(dds)
    Healthy
    Every known technical covariate that isn't the thing you're testing is included in the formula, e.g. `~ batch + condition`.
    Red flag
    `design(dds)` is `~ condition` only, despite a documented batch, prep-date, or lane variable sitting unused in colData.
  4. Tally the number of samples per condition using your sample sheet or colData(dds), and compare it against the common 3-replicate default.

    Healthy
    More than the bare minimum of 3 per group where the budget allows it; small n makes dispersion estimates noisy and destabilizes the whole histogram shape.
    Red flag
    n ≤ 3 per group combined with a hill-shaped or bimodal-looking histogram. Treat the p-values as underpowered before treating them as wrong.
  5. Look at the distribution of baseMean (or mean counts) for the genes DESeq2 actually tested.

    r
    hist(log10(res$baseMean + 1), breaks = 50)
    sum(res$baseMean < 10, na.rm = TRUE)
    Healthy
    A minority of tested genes sit in the very-low-count tail; most genes have enough counts for the negative binomial fit to be stable.
    Red flag
    A large fraction of tested genes have baseMean near zero. Their p-values pile up near 1, and after multiple-testing correction that pile-up produces the U shape.
  6. Run a variance-stabilized PCA and color by both condition and every technical covariate you have.

    r
    vsd <- vst(dds, blind = TRUE)
    plotPCA(vsd, intgroup = c("condition", "batch"))
    Healthy
    The axis that explains the most variance (PC1) tracks your biological condition, or if it tracks a technical variable, that variable is already in the design formula.
    Red flag
    PC1 or PC2 cleanly separates samples by processing day, replicate number, or batch instead of condition, and that variable is not in `design(dds)`.
  7. Plot gene-wise dispersion estimates against the fitted trend.

    r
    plotDispEsts(dds)
    Healthy
    Gene-wise (black) points shrink cleanly toward the red fitted trend line; blue points (post-shrinkage) sit close to the curve.
    Red flag
    A large scatter of points capped far above the trend, or the trend itself looks wrong for your data, and your p-value histogram has a hump instead of a flat baseline. This points to overestimated dispersion or a group-to-group variance mismatch the NB model can't absorb.
  8. Refit the same counts and the same design in edgeR (glmQLFit/glmQLFTest) or limma-voom and compare histogram shapes.

    Healthy
    Broadly similar shapes: spike near zero, flat elsewhere. Minor differences in spike height are normal since the tools use different tests.
    Red flag
    DESeq2 gives a hill shape but edgeR or limma-voom on the identical counts and design gives a normal spike-plus-flat shape (or vice versa). That disagreement means the problem is in one tool's fit, not in your samples.

What to do about it

Pre-filter low-count genes before testing, not after

When: The histogram is U-shaped and you confirmed a large fraction of tested genes sit at very low baseMean.

Filter the count matrix (or the dds object) to keep genes with, for example, counts ≥5 in at least 10% of samples, then run DESeq() fresh on the filtered object so dispersion estimation and testing both see the cleaner gene set.

Caveat: This changes the number of tests performed, which changes every padj value, even for genes you didn't touch. Filter only on count magnitude, never on padj or fold change, or you're shaping the histogram to look the way you want instead of fixing the model.

Add the batch term to the design formula

When: PCA shows a technical variable (batch, prep date, lane) separating samples, and the histogram is hill-shaped.

Rebuild dds with the covariate in the formula, e.g. design = ~ batch + condition, and refit dispersions and the GLM from that object. Let DESeq2 estimate and account for the batch term inside the model.

Caveat: This only works if batch and condition are not fully confounded. If every batch corresponds to exactly one condition, no formula separates their effects; you need new samples with batch and condition crossed, not a statistical fix.

Use a dispersion estimator more robust to small designs

When: `plotDispEsts()` shows poor shrinkage or clearly unequal variance between groups, and n per group is small.

Switch to edgeR's quasi-likelihood pipeline (glmQLFit + glmQLFTest), which is designed to handle dispersion uncertainty in small-replicate designs better than a plain likelihood-ratio test.

Caveat: A different test handles moderate dispersion misestimation better, it does not fix a genuinely broken design like confounded batch or a mislabeled sample. Diagnose the cause first.

Trace and, if justified, drop an outlier sample

When: One sample stands apart on PCA or has unusually low correlation with its replicates, and you suspect a lab-side error (swap, degraded RNA, failed prep).

Check raw count correlations between replicates and cross-reference the sample sheet for a documented lab issue tied to that specific sample. Only drop it if you can point to a concrete cause.

Caveat: Dropping a sample lowers your n and your power. Never remove a sample just because the histogram looks better without it; that's p-hacking the diagnostic itself.

When not to "fix" it

If you set lfcThreshold above 0 in results(), a U-shaped histogram is the expected output of DESeq2's composite null test, not a bug. Don't chase a fix for it.

If the histogram shows a sharp spike near zero and a genuinely flat tail, and the perturbation is strong (a knockout of a hub gene, a drug at a saturating dose, a developmental time-course endpoint), a large fraction of low p-values is real biology. This is the field's "anti-conservative but expected" category, not a calibration problem. Don't dilute it by adding covariates that aren't actually in your design just to make the histogram look more textbook.

A histogram that's uniformly flat everywhere, no spike near zero at all, can be a correct result if your treatment genuinely produced no detectable differential expression at your depth and replicate number. That's a disappointing result, not a broken one. Don't keep relaxing filters or switching tools until you manufacture a spike.

Five things experienced analysts do here

  1. Run `hist(res$pvalue)` before summary() or a volcano plot every single time. It's one line and it tells you whether the numbers downstream are worth trusting at all.
  2. Look at the PCA before you ever call DESeq() or glmQLFit(). A batch axis is far cheaper to spot in a PCA plot than to diagnose after the fact from a warped p-value histogram.
  3. Never filter genes by padj or fold change to make the histogram look better. Filter on raw count magnitude, and do it before the test runs, or you're curve-fitting your own statistics.
  4. Keep `design(dds)` in sync with your full sample sheet, including covariates you didn't manipulate. The histogram diagnoses the model, and the fix almost always lives in the formula, not in a pre-corrected count matrix.
  5. If two tools give you visibly different histogram shapes on the same counts and the same design, trust neither until you understand why. That disagreement means the model is unstable, not that one tool is simply 'better'.

Questions people ask

What does a U-shaped p-value histogram mean in DESeq2?

After multiple-testing correction, a U shape almost always comes from a large number of low-count genes whose p-values pile up near 1 because the negative binomial fit is unstable at low expression. Filtering genes with counts below roughly 5 in at least 10% of samples before running DESeq2, not after, typically removes it. If you deliberately set lfcThreshold above 0, a U shape is the expected result of the composite null test and isn't a problem.

Why does my p-value histogram have a hump around 0.5 instead of being flat?

A hill or hump shape means your p-values aren't well calibrated. The usual causes are an unmodeled batch effect, overestimated dispersion, or unequal within-group variability between conditions. Check your PCA for a batch axis and your design formula for missing covariates before assuming the data itself is bad.

Is a spike at p = 0 in edgeR always a sign of real differential expression?

Not necessarily. edgeR often shows a spike at exactly p = 0 even when there's no true DE signal, because very small p-values underflow to zero in floating-point arithmetic rather than representing literally impossible probabilities. Look at the shape of the rest of the distribution, not just the height of the zero bin.

What does a healthy p-value histogram look like for bulk RNA-seq?

A sharp spike in the lowest bin from genuinely differentially expressed genes, sitting on top of an approximately flat distribution across the rest of [0,1] from non-DE genes. For count-based tests a small additional bump right at p = 1 is normal and not a red flag by itself.

Should I filter low-count genes before or after running DESeq2?

Before. Filter on raw count magnitude, for example counts ≥ 5 in at least 10% of samples, then run DESeq2 from scratch. Filtering by padj or fold change after the fact doesn't fix the histogram shape, it just hides the genes that were dragging it down, and DESeq2's own independent filtering (which sets padj to NA for low-mean genes) is a separate mechanism from pre-filtering the count matrix.

Related pages

Related reading on the blog

Sources

  1. DESeq2 dispersion and p-value histograms — Hill-shaped histograms from overestimated dispersion or unequal within-group variability
  2. U-shaped p-value distributions in DESeq2 RNA-seq analysis — U shape is expected, not a bug, when lfcThreshold is set
  3. U-shape removal by filtering for genes with >5 counts in >10% of samples — Concrete pre-filtering fix for U-shaped histograms
  4. Hill-shaped p-value histogram in DESeq2 — Hill shape as a sign of miscalibrated p-values
  5. A field-wide assessment of differential expression profiling by high-throughput sequencing reveals widespread bias — 75% of 557 experiments had malformed histograms; category breakdown and low-count filtering improvement
  6. Batch Effect: To Correct or Not for Bulk RNA-seq Data — Design formula fix for batch-driven hill-shaped histograms
  7. Understanding p value, multiple comparisons, FDR and q value — Background on interpreting p-values, FDR, and q-values
  8. DESeq2 package page, Bioconductor — Tool reference
  9. edgeR package page, Bioconductor — Tool reference

Part of the P-value histograms series.