Chatomics Field GuideWhat They Don't Teach You

Sanity check · Bulk RNA-seq

How to Handle Multiple Testing and FDR in Bulk RNA-seq

Zero genes at padj < 0.05 with a clean PCA is a design or filtering problem to diagnose, not a cutoff to loosen.

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

You ran DESeq2 or edgeR on 3-6 replicates per condition, the PCA plot separates your groups cleanly on PC1, and you're confident the treatment did something. Then results() comes back with zero genes at padj < 0.05, or a padj column that's mostly NA. The temptation is to lower the cutoff, drop the "outlier" sample, or quietly switch the write-up to raw p-values until the list looks respectable.

That temptation is exactly what independent filtering, Benjamini-Hochberg, and Cook's-distance outlier detection exist to prevent. Each of those steps changes the effective number of tests (m) the correction divides by, and a mismatch between your design formula, your alpha, and how you're interpreting NA values is what turns real signal into an empty results table, or turns an empty table into a false claim of "no effect."

This page gives you an ordered set of checks, cheapest first, to find out whether zero significant genes means a broken design, a filtering mismatch, or an honest null result, so you can decide what to actually fix in the next hour instead of re-running the same test with a different cutoff.

What it looks like when it's happening

  • `results(dds)` returns 0 rows at padj < 0.05, but the raw `pvalue` column has hundreds of values below 0.05
  • `plotPCA(vst(dds))` shows condition cleanly separating samples on PC1 or PC2
  • Most of the `padj` column is `NA` instead of a number
  • The gene count barely moves between padj < 0.1 and padj < 0.05, or collapses to zero when you tighten it
  • A volcano plot shows a cloud of points below the raw p = 0.05 line but none colored significant after adjustment
  • `summary(res)` reports "0 genes with padj < 0.05" alongside a large count of low-count or outlier-filtered genes
  • Re-running the identical contrast in edgeR or limma-voom gives a very different number of hits
  • Toggling `cooksCutoff` or `independentFiltering` in `results()` swings the significant gene count a lot

Why it happens

Benjamini-Hochberg ranks all m tested genes by p-value and multiplies each by m/rank. m is not 20,000 by default, it's whatever DESeq2 (or edgeR/limma) decided actually got tested after independent filtering drops genes with low mean normalized counts and after Cook's-distance outlier detection sets padj to NA for genes with one influential sample. Skip pre-filtering and thousands of near-zero-count genes with unstable, near-uniform p-values inflate m, so every real gene's adjusted p-value grows because it's competing against a denominator full of noise. This is also why filtering low-count genes before testing can paradoxically produce more significant hits, not fewer: it shrinks m.

On the biology side, replicate variability and sequencing-depth imbalance widen the fitted negative-binomial dispersion for every gene. With 3 replicates per group, one noisy or shallow sample pulls the dispersion curve (plotDispEsts) upward, and every gene's test statistic shrinks accordingly, even when PC1 shows a clean split. A clean PCA is driven by the top few hundred most variable genes; it says nothing about whether the other 15,000 genes have enough replicates to survive correction.

Design and labeling problems produce the identical symptom through a different mechanism. If a batch variable is confounded with your condition, fully or partially, the model attributes real variance to the wrong term, standard errors balloon, and the coefficient you're testing never separates from zero. If someone fed TPM or FPKM into DESeq2 instead of raw or expected counts, the negative-binomial variance model is fit against numbers it was never built for, and dispersion estimation degrades across the whole matrix. Both look identical downstream: an empty gene list dressed up as "no significant genes."

Depth imbalance adds a second, subtler failure mode: samples sequenced deeper detect more weakly-expressed genes, and a naive comparison calls those genes "up-regulated" purely from depth, not biology. When depth correlates with condition (a common artifact of batch processing samples on different flow cells), this inflates false discoveries at low expression and simultaneously distorts the dispersion curve that FDR correction relies on.

The checks

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

0/9 checked · saved in this browser

  1. Histogram the unadjusted res$pvalue (or qlf$table$PValue for edgeR) before you touch padj at all. BH correction only has something to promote if this distribution has a spike near zero.

    r
    hist(res$pvalue, breaks = 50, col = "grey",
         main = "raw p-value distribution")
    Healthy
    A spike of low p-values near 0 sitting on top of a roughly flat/uniform tail out to 1, that spike is the real signal for BH to work with.
    Red flag
    A flat, uniform histogram across [0,1] means there's no real signal for FDR correction to find; padj will legitimately have few or zero hits. A histogram spiking near 1 instead usually means a swapped contrast or coding error.
  2. Read the independent-filtering and outlier lines straight out of summary(res), then plot the filtering diagnostic to see where the rejection curve peaks.

    r
    summary(res)
    metadata(res)$filterThreshold
    plot(metadata(res)$filterNumRej, type = "b",
         xlab = "quantiles of filter", ylab = "number of rejections")
    Healthy
    `filterNumRej` rises then falls, peaking at some intermediate quantile, that's the threshold DESeq2 picked, rejecting a modest fraction of genes (the truly low-count ones), not most of the genes tested.
    Red flag
    Thousands of genes with padj = NA from independent filtering, or a rejection curve that never peaks, m is inflated because low-count genes weren't pre-filtered before `DESeqDataSet` was built.
  3. Independent filtering is optimized for the FDR target set by alpha in results(), default 0.1. If you report padj < 0.05, rerun with alpha=0.05 explicitly rather than filtering the default table yourself.

    r
    res05 <- results(dds, alpha = 0.05)
    summary(res05)
    Healthy
    The gene count is stable or slightly larger than what you'd get by filtering the default-alpha table down to padj < 0.05 afterward.
    Red flag
    You've been filtering a default alpha=0.1 results table for padj < 0.05 downstream, the independent-filtering threshold was optimized for the wrong FDR and is quietly discarding power.
  4. Separate the three legitimate reasons DESeq2 sets padj to NA: zero counts across all samples, Cook's-distance outliers, and independent filtering.

    r
    sum(res$baseMean == 0)
    sum(!is.na(res$pvalue) & is.na(res$padj))
    table(is.na(res$pvalue), is.na(res$padj))
    Healthy
    NAs concentrate in `baseMean == 0` genes (never expressed) plus a small independent-filtering exclusion.
    Red flag
    A large block of NA padj with nonzero baseMean and no Cook's flag, low-count genes weren't pre-filtered before building `dds`, inflating m and burying real hits.
  5. Pull the per-gene, per-sample Cook's distance matrix and check whether a single sample drives outlier calls across many genes.

    r
    cooks <- assays(dds)[["cooks"]]
    boxplot(log10(cooks), range = 0, main = "Cook's distance per sample")
    Healthy
    Cook's distances are similar in spread across all samples, no single sample stands out.
    Red flag
    One sample's boxplot sits visibly above the rest, it's driving NA padj (and possibly the confounding pattern); decide whether to exclude it rather than just setting `cooksCutoff = FALSE`.
  6. Compare sizeFactors(dds) across samples and look at whether the dispersion plot's fitted curve tracks the point cloud.

    r
    sizeFactors(dds)
    plotDispEsts(dds)
    Healthy
    Size factors within roughly 2-fold of each other across samples; gene-wise dispersions (black) shrinking cleanly toward the fitted curve (red), with final estimates (blue) close to it.
    Red flag
    Size factors differing 5-10x between samples, or gene-wise dispersions scattered far above the fitted curve, depth imbalance or high inter-replicate variance is widening every gene's test statistic.
  7. Print the model formula and the exact coefficient name, then cross-tabulate condition against any batch/covariate term.

    r
    design(dds)
    resultsNames(dds)
    table(colData(dds)$condition, colData(dds)$batch)
    Healthy
    The coefficient you extracted matches your condition of interest, and the condition-by-batch table has both batch levels represented in both conditions.
    Red flag
    The condition-by-batch table is diagonal (each batch maps to only one condition), that's full confounding, and no formula change fixes it after the fact.
  8. Check that the matrix dds was built from is non-negative integers (or near-integer expected counts from tximport), not normalized abundances.

    r
    head(counts(dds))
    any(counts(dds) %% 1 > 0.01)
    Healthy
    Whole numbers (or near-whole for tximport expected counts) ranging from 0 to the tens of thousands, varying freely in total per sample.
    Red flag
    Small decimals bounded roughly 0-1000 summing to a fixed total per sample, that's TPM/FPKM, and the negative-binomial model DESeq2/edgeR fits is invalid on it.
  9. Re-run the identical two-group contrast in edgeR or limma-voom with explicit BH adjustment and compare hit counts and top genes.

    r
    qlf <- glmQLFTest(fit, coef = 2)
    topTags(qlf, n = Inf, adjust.method = "BH")
    Healthy
    Roughly concordant top genes and a similar order-of-magnitude hit count between tools; small differences are normal.
    Red flag
    One tool returns hundreds of hits and another returns zero on the identical contrast and samples, that's a tool-specific design or filtering setting, not biology.

What to do about it

Pre-filter low-count genes before running DESeq2/edgeR, not after

When: The independent-filtering diagnostic (checks 2 and 4) shows m inflated by thousands of near-zero-count genes.

Filter the count matrix before constructing dds, e.g. keep genes with rowSums(counts(dds) >= 10) >= smallest_group_size, then rebuild the DESeqDataSet and rerun dispersion estimation and testing from scratch, dispersion fitting depends on which genes are included.

Caveat: Filtering too aggressively drops lowly-expressed but biologically important genes; pick the threshold from your design's smallest group size, not a fixed number copied from a tutorial.

Match alpha to your reporting cutoff

When: You report padj < 0.05 but ran `results()` with the default alpha = 0.1.

Call results(dds, alpha = 0.05) so independent filtering optimizes for the same FDR you're reporting, and use that object for both the gene count and downstream tables.

Caveat: This aligns the filtering threshold with your stated cutoff, it doesn't add power, don't expect a dramatically larger gene list from this alone.

Add the batch term to the design instead of pre-correcting counts

When: The condition-by-batch table (check 7) shows partial overlap, not full confounding.

Refit with design = ~ batch + condition and re-extract the condition coefficient. Do not run limma::removeBatchEffect() on the matrix and feed the corrected values into DESeq2/edgeR/limma's test.

Caveat: If batch and condition are fully confounded (each batch maps to exactly one condition), no formula fixes it, you need new samples that break the confound; modeling around full confounding invents statistical power that isn't there.

Re-import from raw counts if TPM/FPKM went into the model

When: Check 8 shows decimal, capped-scale values instead of integers.

Re-run tximport requesting counts for your quantifier (salmon/kallisto), or re-pull raw gene counts from your aligner's count table, then rebuild the DESeqDataSet/DGEList from scratch.

Caveat: This is a full redo of the DE step, not a patch on existing results; anything downstream (volcano plots, GSEA rankings) needs to be regenerated too.

Inspect and, if justified, exclude a Cook's-outlier sample rather than disabling cooksCutoff globally

When: Check 5 shows one sample dominating Cook's distances across many genes, and that sample also looks off on PCA or sample correlation.

Confirm the sample is genuinely an outlier via an independent signal (PCA position, degradation metrics, a suspected label swap), remove it, and refit dispersion and testing on the remaining replicates. Reserve results(dds, cooksCutoff = FALSE) for cases with few replicates where you're confident no single sample is driving false positives.

Caveat: Dropping a replicate reduces degrees of freedom and can make everything less significant if you're already at 3 per group; disabling cooksCutoff instead lets outlier-driven false positives back into the list.

Use shrunken effect sizes or IHW when the raw-vs-adjusted gap is real, not a filtering bug

When: The p-value histogram (check 1) shows a genuine spike near zero and every other check comes back clean, but the effect is simply underpowered at your replicate count.

Run lfcShrink(dds, coef = ..., type = 'ashr') for reliable effect sizes independent of the FDR test, and consider results(dds, filterFun = ihw) to weight genes by mean expression instead of a hard filter, which can recover a few extra genes without inflating FDR.

Caveat: Neither shrinkage nor IHW manufactures power your experiment doesn't have; if you're still at zero significant genes with 3 replicates per group after this, the honest fix is more replicates, not a different filtering function.

When not to "fix" it

If the raw p-value histogram (check 1) is flat and uniform across [0,1] instead of spiking near zero, and the condition-by-batch table and count matrix check out clean, zero significant genes is the correct answer: the experiment doesn't have a detectable transcriptional effect at this depth and replicate count. Don't chase significance by lowering alpha, dropping samples until the model "improves," or quietly reporting raw p-values instead of padj in the write-up, that's p-hacking. A small log2 fold change reported honestly at p = 0.06 is more useful to your collaborators than a manufactured padj < 0.05 list that won't replicate.

Five things experienced analysts do here

  1. Plot the raw p-value histogram before you ever open the padj column, its shape tells you in ten seconds whether BH correction has any real signal to work with.
  2. Report log2 fold change (ideally after `lfcShrink`) next to padj in every table and plot, a long list of tiny padj values with fold changes under 1.2x is usually a depth or replicate-count artifact, not biology.
  3. Set alpha in `results()` to whatever padj cutoff you actually report; independent filtering is optimized for that number, and leaving it at the default 0.1 while filtering the table yourself for 0.05 throws away power silently.
  4. Print `resultsNames(dds)` and `design(dds)` before trusting any gene count, testing the wrong coefficient produces the exact same empty-list symptom as a real FDR problem.
  5. Never run `limma::removeBatchEffect()` on the count or expression matrix before feeding it to DESeq2/edgeR/limma's test, correct for batch inside the model formula, or you'll bias your p-values and manufacture significance instead of controlling it.

Questions people ask

Why do I have zero significant genes after FDR correction in DESeq2?

Usually one of three things: too many low-count genes inflated the number of tests BH divides by, dispersion is high because of few replicates or depth imbalance, or your design has batch confounded with condition so the coefficient you're testing never separates from zero. Check the raw p-value histogram first, if it's flat, there may genuinely be no detectable signal at your current depth and replicate count.

What's the difference between padj and p-value in DESeq2 results?

The p-value is the per-gene test result before any correction for testing thousands of genes at once. padj applies Benjamini-Hochberg to control the false discovery rate across all genes tested, so at padj < 0.05 you expect about 5% of your significant list to be false positives, not 5% of all 20,000 genes tested.

Should I use padj < 0.05 or raw p-value < 0.05 to call genes significant?

Use padj. Testing 20,000 genes at raw p < 0.05 gives roughly 1,000 false positives from chance alone. padj < 0.05 controls the expected fraction of false positives within your significant gene list, which is the number you actually want when reporting hits.

What does independent filtering do in DESeq2 and why does it change my gene list?

It removes genes with low mean normalized counts before multiple testing correction, because those genes have too little power to ever reach significance and only add noise to the correction. Removing them shrinks the total number of tests (m) in the BH formula, which can paradoxically increase the number of significant genes among those that remain.

Why does changing the alpha parameter in DESeq2's results() function change my number of significant genes?

alpha sets the target FDR that independent filtering is optimized for, with a default of 0.1. If you're reporting padj < 0.05 but left alpha at the default, the filtering threshold wasn't tuned for your actual cutoff, and you likely have fewer significant genes than you should. Set alpha=0.05 in results() to match.

Related pages

Related reading on the blog

Sources

  1. Gene-level differential expression analysis with DESeq2, HBC Training — BH formula, independent filtering mechanics, and how padj cutoffs translate into expected false positives
  2. DESeq2: Differential gene expression analysis based on the negative binomial distribution, Bioconductor — alpha parameter, independent filtering, and IHW integration via results(filterFun=ihw)
  3. Outliers and Filtering with DESeq2 to overcome NA values in padj, Bioconductor Support — The three sources of NA padj: zero counts, Cook's distance outliers, and independent filtering
  4. Are published RNA seq data analyses often wrong in calculating p-values and FDR?, Bioconductor Support — Troubleshooting zero or unexpected significant gene counts and pointing to design/data-quality causes
  5. Group Heteroscedasticity - A Silent Saboteur of Power and False Discovery in RNA-Seq Differential Expression, bioRxiv — How unequal group variances distort FDR control and inflate type I/II errors
  6. Understanding p value, multiple comparisons, FDR and q value, Tommy Tang's Blog — Bonferroni vs Benjamini-Hochberg vs Storey's q-value explained conceptually

Part of the Multiple testing and FDR series.