Reading RNA-seq differential expression results: padj, fold-change cutoffs, NA values and batch

Quick answer. In a DESeq2, edgeR or limma table, call genes differentially expressed by the adjusted p-value (padj or FDR, Benjamini–Hochberg) rather than the raw p-value: in our example, raw p < 0.05 called 344 genes including 75 that do not change at all; padj < 0.05 called 238 including 5. padj answers “did this gene change at all?”, so many called genes change only a little. Filtering by |log2 fold change| > 1 afterwards gives a useful list of genes estimated to change more than 2-fold, but it does not test that claim (in our simulation with 4 samples per group, about 6% of them truly changed less than 2-fold or not at all). If the size matters, test against the threshold (DESeq2 lfcThreshold, limma TREAT), which needs more samples: with 4 per group it found about a quarter of the truly 2-fold genes. NA in padj usually means the gene was filtered for low counts (or had all-zero counts or an extreme count outlier); leave such genes out of the number tested.

Terms used below. Differential expression (DE): a gene’s expression differs between conditions. Log2 fold change (LFC): log2 of the ratio between conditions; 1 means 2-fold up, −1 means 2-fold down. Raw p-value: the test result for one gene, without accounting for testing thousands. padj (adjusted p-value, FDR): the Benjamini–Hochberg adjusted p-value; calling genes with padj < 0.05 keeps the expected share of false calls among them at or below 5% when the tests are independent or positively related, the usual working assumption for genes. Threshold test (TREAT, lfcThreshold): tests whether |LFC| is greater than a chosen value, rather than greater than 0. Independent filtering: DESeq2 removes genes with low average counts before adjusting p-values, which leaves their padj as NA. LFC shrinkage: pulls noisy fold-change estimates of low-count genes towards 0, for ranking and plots. Volcano plot: −log10 p (or adjusted p) against log2 fold change for all genes. Wald test: DESeq2’s default test of a fold change against 0 (or against a threshold). apeglm: a shrinkage method for fold changes. Design formula: the model DESeq2 fits, such as ~ batch + condition (condition last). vst: DESeq2’s variance-stabilising transformation, giving log2-like values. coldata: the table describing each sample (condition, batch), in the same order as the count columns; resultsNames(dds) lists the coefficient names. Batch: samples processed together (extraction, library preparation, sequencing run). Power: the chance of detecting an effect that is really there. False positive: a “significant” result when there is no real effect.

Which analysis for your question

Question	Analysis	Report
Which genes changed at all?	padj (BH) < 0.05 from DESeq2, edgeR or limma	Number of genes, padj threshold, n per group
Which genes changed by more than 2-fold?	Threshold test: DESeq2 results(lfcThreshold = 1), limma treat, edgeR glmTreat; needs more samples than the test for any change	Genes with padj < 0.05 for |LFC| > 1
A shortlist of strongly changed genes for follow-up	padj < 0.05, then sort or filter by shrunken LFC; describe it as a filter	“padj < 0.05 and estimated |LFC| > 1”
Samples processed in batches	Batch in the design (~ batch + condition), with both conditions in every batch; if each batch holds one condition, no analysis can separate them	As above, adjusted for batch
Ranking and plots	Shrunken LFC (apeglm) for volcano plots and rankings	—

Love and colleagues describe DESeq2’s padj: “the Wald test P values from the subset of genes that pass an independent filtering step … are adjusted for multiple testing using the procedure of Benjamini and Hochberg”, and warn about cutoffs: “A common procedure is to disregard genes whose estimated LFC … is below some threshold … However, this approach loses the benefit of an easily interpretable FDR, as the reported P value and adjusted P value still correspond to the test of zero LFC” (Love et al. 2014). McCarthy and Smyth make the same point and test “the thresholded null hypothesis” instead (McCarthy & Smyth 2009). The DESeq2 vignette explains that if a gene “is filtered by automatic independent filtering, for having a low mean normalized count, then only the adjusted p value will be set to NA”, that shrinkage “is useful for visualization and ranking of genes”, and that “if the condition samples are balanced across experimental batches, by including the batch factor to the design, one can increase the sensitivity” (DESeq2 vignette). How many genes you find depends on replicates: “with three biological replicates, nine of the 11 tools evaluated found only 20%–40% of the significantly differentially expressed (SDE) genes identified with the full set of 42 clean replicates” (Schurch et al. 2016).

Worked example: 4 vs 4 samples, 2,000 genes

Constructed data: normalised log2 expression for 2,000 genes in 4 control and 4 treated samples (fewer than the six per condition the site recommends for RNA-seq, but a common design). By design, 1,700 genes stay the same and 300 change, with true |log2 fold change| evenly spread from 0.25 to 2 (1.2-fold to 4-fold), half up and half down; 171 of them change by more than 2-fold. Genes vary with a standard deviation of about 0.25 on the log2 scale between samples (a quiet assay; noisier data lower every detection rate below). Each gene is tested with a pooled-variance t-test on log2 expression, a simple stand-in for DESeq2, edgeR or limma, which share information across genes and have more power; for your own data use those tools (or Welch’s t-test if you must test genes one at a time). The random seed was the first one tried; the true fold changes were revised once, after a first run with only two effect sizes, so that some lie near the 2-fold threshold. Stored code and data: log2_expression.csv, truth.csv, model.py, make_data.py, analysis.py, analysis.R and the analysis output; analysis.py and analysis.R print the same output.

Rule	Genes called	Truly unchanged	True change below 2-fold (including unchanged)
Raw p < 0.05	344	75	173
padj < 0.05	238	5	69
padj < 0.05 and estimated |LFC| > 1, filtered after testing	176	0	10
Threshold test |LFC| > 1 (TREAT-style, on the t-test), adjusted p < 0.05	32	0	0

Raw p-values let 75 unchanged genes through; padj keeps that to 5, as intended. Of the 238 padj genes, 69 change by less than 2-fold, which is not an error: padj tests “any change”. The fold-change filter leaves 176 genes, of which 10 truly change less than 2-fold, so a sentence like “176 genes changed more than 2-fold” would be wrong for some of them. The threshold test controls the false discovery rate for that claim but, with 4 samples per group, called only 32 of the 171 genes that truly change more than 2-fold. In 2,000 simulated studies of this design (simulation.py, the simulation output):

Samples per group	Rule	Genes called (mean)	Share with a true change below 2-fold (including unchanged)	Genes truly above 2-fold found
4	padj < 0.05	234	29.5%	96.4%
4	padj < 0.05 and |LFC| > 1 after testing	167	6.4%	91.5%
4	Threshold test, adjusted p < 0.05	42	0.0%	24.8%
8	padj < 0.05	287	40.4%	100.0%
8	padj < 0.05 and |LFC| > 1 after testing	171	4.9%	95.3%
8	Threshold test, adjusted p < 0.05	104	0.0%	60.8%

The filter is a reasonable way to pick genes for follow-up, and in this simulation 6.4% (4 per group) or 4.9% (8 per group) of its genes truly changed by less than the threshold or not at all. The threshold test, which keeps such genes at or below 5% of its calls, let none through here but needs more samples; doubling the group size from 4 to 8 raised the share of truly 2-fold genes it found from 24.8% to 60.8%.

The code expects one row per gene with a gene column and one column of normalised log2 expression per sample (for example DESeq2’s vst output or log2 counts per million from edgeR), named control_… and treated_…; truth.csv exists only for the constructed example. Rename your columns to match, and run the code in the folder that holds the CSV files. Python needs NumPy, pandas and SciPy 1.11 or newer (for false_discovery_control); the R code uses base R only.

# Condensed: analysis.py also compares the rules with the known truth
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("log2_expression.csv")
ctl, trt = d.filter(like="control_").to_numpy(), d.filter(like="treated_").to_numpy()
lfc = trt.mean(1) - ctl.mean(1); n1, n2 = ctl.shape[1], trt.shape[1]; df = n1 + n2 - 2
se = np.sqrt(((n1 - 1) * ctl.var(1, ddof=1) + (n2 - 1) * trt.var(1, ddof=1)) / df * (1 / n1 + 1 / n2))
padj = stats.false_discovery_control(2 * stats.t.sf(np.abs(lfc) / se, df))          # any change
p_tr = stats.t.sf((np.abs(lfc) - 1) / se, df) + stats.t.sf((np.abs(lfc) + 1) / se, df)
padj_tr = stats.false_discovery_control(p_tr)                                       # change > 2-fold
print((padj < 0.05).sum(), ((padj < 0.05) & (np.abs(lfc) > 1)).sum(), (padj_tr < 0.05).sum())

# With DESeq2 on raw counts (the usual route; not run on this page's log2 table):
# dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ batch + condition); dds <- DESeq(dds)
# res <- results(dds, alpha = 0.05)                                   # padj for any change
# res2 <- results(dds, lfcThreshold = 1, altHypothesis = "greaterAbs", alpha = 0.05)   # padj for |LFC| > 1
# shr <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")   # for plots and ranking
d <- read.csv("log2_expression.csv")
ctl <- as.matrix(d[, grep("^control_", names(d))]); trt <- as.matrix(d[, grep("^treated_", names(d))])
lfc <- rowMeans(trt) - rowMeans(ctl); n1 <- ncol(ctl); n2 <- ncol(trt); df <- n1 + n2 - 2
se <- sqrt(((n1 - 1) * apply(ctl, 1, var) + (n2 - 1) * apply(trt, 1, var)) / df * (1 / n1 + 1 / n2))
padj <- p.adjust(2 * pt(abs(lfc) / se, df, lower.tail = FALSE), "BH")
padj_tr <- p.adjust(pt((abs(lfc) - 1) / se, df, lower.tail = FALSE) + pt((abs(lfc) + 1) / se, df, lower.tail = FALSE), "BH")
c(any = sum(padj < 0.05), filtered = sum(padj < 0.05 & abs(lfc) > 1), threshold = sum(padj_tr < 0.05))

Before you analyse

1. Decide the claim in advance: “changed” (padj) or “changed by more than X-fold” (threshold test), and the padj level.
2. Use biological replicates and plan their number (how many biological replicates); adding replicates raises power far more than deeper sequencing, which gives “diminishing returns on power” after about 10 million reads in one study (Liu et al. 2014).
3. Record and balance batches (both conditions in every batch), and put batch in the design (Leek et al. 2010 describe how batch effects distort comparisons).
4. Keep NA padj values out of the count of tested genes and report how many genes were tested after filtering; DESeq2 also sets NA for genes with all-zero counts or an extreme count outlier.
5. Use shrunken fold changes for volcano plots and rankings, and draw the padj line, not a raw p-value line.

Check what an AI assistant ran

1. The assistant called genes by raw p < 0.05. Ask for padj (BH) < 0.05.
2. The assistant reported “genes changed more than 2-fold” after filtering padj genes by fold change. Ask for a threshold test, or rephrase as “estimated fold change above 2”.
3. The assistant treated NA padj as non-significant or dropped them silently. Ask how many genes were filtered and why.
4. The assistant ignored batch. Ask for batch in the design.

Report it

Differential expression was tested with DESeq2 ([version]) using the design ~ [batch] + condition ([4] samples per condition). Genes with Benjamini–Hochberg adjusted p < 0.05 were called differentially expressed; genes changing by more than 2-fold were identified with a Wald test against |log2 fold change| > 1 (lfcThreshold = 1). Log2 fold changes were shrunk with apeglm for plots.

(Example numbers from this page’s constructed data.) [238] genes were differentially expressed (padj < 0.05), of [2,000] genes tested after filtering; [32] changed by more than 2-fold (padj < 0.05 for |log2 fold change| > 1).

If you filtered by fold change instead of testing against the threshold, write “[176] of these had an estimated |log2 fold change| above 1”. Replace the bracketed values with your own.

References

• Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550. doi:10.1186/s13059-014-0550-8
• McCarthy DJ, Smyth GK (2009). Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics 25:765–771. doi:10.1093/bioinformatics/btp053
• Love MI, Anders S, Huber W. Analyzing RNA-seq data with DESeq2. Bioconductor vignette (checked 12 November 2026). bioconductor.org
• Schurch NJ, Schofield P, Gierliński M, et al. (2016). How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22:839–851. doi:10.1261/rna.053959.115
• Liu Y, Zhou J, White KP (2014). RNA-seq differential expression studies: more sequence or more replication? Bioinformatics 30:301–304. doi:10.1093/bioinformatics/btt688
• Leek JT, Scharpf RB, Bravo HC, et al. (2010). Tackling the widespread and critical impact of batch effects in high-throughput data. Nature Reviews Genetics 11:733–739. doi:10.1038/nrg2825