What log2 Fold Change Cutoff to Use in RNA-seq

If you want a single number to take away, test against a threshold of log2(1.5) = 0.585 rather than filtering afterward at |log2FC| ≥ 1. In DESeq2 that call is results(dds, lfcThreshold = log2(1.5), altHypothesis = "greaterAbs", alpha = 0.05), and in edgeR it is glmTreat(fit, coef = 2, lfc = log2(1.5)). Both compute a p-value for the hypothesis that the true effect exceeds the threshold, so the adjusted p-value you get out means what you think it means. The more common practice of taking padj < 0.05 from a test against zero and then discarding everything with |log2FC| < 1 gives you a gene list with no controlled error rate for the fold-change claim, and it systematically enriches for noisy, low-expression genes.

The sections below explain what the quantity is and why the post hoc filter fails. They then cover how to move the threshold inside the test and how to choose the threshold for your own design.

What log2 fold change is, and how to convert it

Before discussing thresholds, it helps to be precise about the quantity being thresholded. Log2 fold change is log₂(expression in condition B / expression in condition A), computed after normalization for library size and composition. The log scale makes the measure symmetric: a doubling is +1, a halving is −1, and no change is 0. On the raw ratio scale those same three cases would be 2, 0.5, and 1. That is awkward both for plotting and for the linear models underneath DESeq2 and edgeR, since both fit a negative binomial generalized linear model on the log scale.

Converting back to a fold change means raising 2 to the power of the value. A log2FC of 1.5 is 2^1.5 = 2.83-fold up, and a log2FC of −0.585 is 2^−0.585 = 0.67, a 33% decrease. Going the other way, a 1.5-fold increase is log₂(1.5) = 0.585, and a 20% increase is log₂(1.2) = 0.263. Negative values are ordinary and simply mean the gene is lower in the numerator condition, which is set by the contrast you specify rather than by any preference of the software.

Why a post hoc fold-change filter breaks the statistics

The trouble with filtering after the fact comes from a mismatch between the hypothesis you tested and the claim you end up making. The standard Wald or likelihood-ratio test asks whether the true log2FC differs from zero. A gene can pass that test at any effect size if it is measured precisely enough, and a highly expressed, tightly replicated gene may reach padj = 1e-12 with a log2FC of only 0.3. Conversely, a gene with 15 counts in one sample and 3 in another has a huge estimated log2FC with almost no information behind it.

When you intersect padj < 0.05 with |log2FC| > 1, you are applying a second selection step whose error rate was never modeled. The Benjamini-Hochberg adjustment was computed over the full set of hypotheses about zero, not over the subset that happened to clear a magnitude bar, so the resulting list does not control the false discovery rate for the statement “this gene changed by at least twofold.”

In practice the filter also inverts the precision ordering. The genes it admits most readily are those with the largest sampling variance in their log ratio, which are the low-count genes. Anyone who has compared one platform to another sees the same effect, where a fixed fold-change cutoff shifts which genes are called depending on dynamic range and detection floor rather than on biology ¹.

Testing against a threshold instead

The remedy is to put the magnitude you care about inside the hypothesis, and both major packages support this. In DESeq2, lfcThreshold changes the null hypothesis from |β| = 0 to |β| ≤ τ, and the Wald statistic becomes (|β̂| − τ) / SE(β̂) under altHypothesis = "greaterAbs". In edgeR, glmTreat implements the TREAT approach. It integrates the likelihood over the interval [−τ, τ] rather than testing a point null.

A DESeq2 workflow with the threshold built in looks like this:

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData, ~ subject + condition)
keep <- rowSums(counts(dds) >= 10) >= 3      # smallest group size
dds <- dds[keep, ]
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "post", "pre"),
               lfcThreshold = log2(1.5), altHypothesis = "greaterAbs",
               alpha = 0.05)
summary(res)

The edgeR equivalent, using quasi-likelihood fitting followed by TREAT, is:

library(edgeR)
y   <- DGEList(counts); y <- y[filterByExpr(y, design), , keep.lib.sizes = FALSE]
y   <- calcNormFactors(y); y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
tr  <- glmTreat(fit, coef = "conditionpost", lfc = log2(1.2))
topTags(tr, n = 50)

Note that the two snippets use different thresholds, and the difference is deliberate. TREAT is conservative by construction, and the edgeR authors recommend a small lfc of typically log2(1.2) to log2(1.5). The reason is that the effective cutoff on the observed fold change ends up well above the nominal one. Setting lfc = 1 in glmTreat usually empties your gene list. DESeq2’s threshold behaves closer to its nominal value, so log2(1.5) is a reasonable default there.

Shrink effect sizes before you rank anything

Testing is one half of the problem; the effect sizes you report and rank are the other. Whatever threshold you use for testing, avoid ranking or plotting the raw maximum-likelihood log2FC. Use lfcShrink(dds, coef = 2, type = "apeglm"), which applies a heavy-tailed prior fit from the data and pulls in estimates from low-count and high-dispersion genes while leaving well-measured large effects nearly untouched. The shrunken values are what you want on an MA plot, in a ranked list for GSEA, and in any downstream integration with proteins or metabolites.

apeglm can also take the threshold directly, in which case it returns s-values rather than p-values:

resS <- lfcShrink(dds, coef = "condition_post_vs_pre", type = "apeglm",
                  lfcThreshold = log2(1.5))
sum(resS$svalue < 0.005, na.rm = TRUE)

An s-value is the estimated false sign-or-magnitude rate of the set you keep, which is a cleaner statement than a filtered FDR.

Choosing the threshold from your own design

There is no universal τ. The right value depends on replicate structure and on what the gene list will be used for, and it should be fixed before you look at results. Three considerations drive the choice.

Replication comes first. With three biological replicates per group, a negative binomial dispersion of 0.1 to 0.3 in human tissue leaves you powered for roughly twofold changes on well-expressed genes and nothing subtler, so a τ near log2(1.5) is a fair reflection of what the experiment can resolve. Longitudinal sampling of a single person behaves differently, because within-person biological variance is much smaller than between-person variance. With eight or more timepoints modeled as ~ subject + timepoint you can resolve changes in the 20 to 40% range on genes above a few hundred counts, and that is where τ = log2(1.2) becomes defensible.

The second consideration is the purpose of the list. If the output feeds pathway enrichment, a magnitude threshold is often counterproductive, since coordinated pathway-level signal is frequently built from many genes moving 20 to 50%. Threshold-free approaches that use the whole ranked vector, including time-course methods that weight by both magnitude and significance, recover biology that a hard cutoff discards ². If instead the output is a shortlist of five genes for a targeted assay, a large τ is appropriate, because there you are optimizing for precision.

The third is consistency across assay types. If you are aligning transcript changes to protein changes from mass spectrometry, keep in mind that the two measurements have different variance structures and different sensible thresholds, and that proteomics fold changes are typically compressed relative to RNA ³. Applying |log2FC| > 1 uniformly to both layers will make the RNA list look rich and the protein list look empty.

For regulatory or otherwise audited work, adopt a published, pre-specified pipeline rather than inventing one. The R-ODAF framework fixes the filtering rules, normalization, and statistical thresholds in advance precisely so that results are reproducible across laboratories ⁴.

Filtering low counts, and where it interacts with the cutoff

Independent filtering of low-count genes matters more than the fold-change cutoff for the shape of your results, so it deserves attention in its own right. DESeq2 applies it automatically, choosing a mean-count quantile that maximizes the number of genes passing at your alpha. It is still worth doing a coarse pre-filter of your own, as in the keep line above, because that stabilizes dispersion estimation and speeds up fitting. The one thing to avoid is filtering on anything that depends on the group labels, such as removing genes with small observed fold changes, since that biases the null distribution and invalidates the FDR.

Putting the pieces together, a reasonable full recipe runs as follows. Quantify with salmon using --validateMappings --gcBias --seqBias and import with tximport at gene level. Filter to genes with ≥10 counts in at least as many samples as your smallest group. Fit with a paired or subject-aware design, test with lfcThreshold = log2(1.5), and report apeglm-shrunken effect sizes. Report the threshold in the methods, because a gene list is uninterpretable without it.

One interpretive caution applies to all of the above: none of this converts an expression change into a statement about your health. Differential expression in whole blood is dominated by cell-type composition, and shifts in neutrophil or lymphocyte fraction will move thousands of genes without anything changing inside any given cell. Deconvolve before you interpret, and take anything clinically consequential to a physician.

Questions people also ask

How do I convert log2 fold change to fold change? Compute 2^log2FC. A log2FC of 1.5 is a 2.83-fold increase, and −2 is a fourfold decrease (2^−2 = 0.25).

What does a fold change of 0.5 mean? Expression is half what it was in the reference condition, equivalent to log2FC = −1.

What log2 fold change is significant? Significance is a property of the test rather than of the magnitude on its own. Set a threshold you can defend from your replicate structure, such as log2(1.5) for small designs and log2(1.2) for dense longitudinal sampling. Put it inside the test with lfcThreshold or glmTreat, then report adjusted p-values or s-values against that threshold.

Is DESeq2 better than edgeR? They agree closely on most datasets. We default to DESeq2 for its apeglm shrinkage and its threshold-aware s-values, and reach for edgeR’s glmQLFit plus glmTreat when the design has many coefficients and quasi-likelihood error rate control matters. Use one, pre-specified.

Can you have a negative log2 fold change? Yes. The sign reflects the direction of your contrast, so always state which condition is the numerator when reporting results.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Rance Nault, Kelly A Fader, Tim Zacharewski. RNA-Seq versus oligonucleotide array assessment of dose-dependent TCDD-elicited hepatic gene expression in mice. BMC Genomics, 2015. https://doi.org/10.1186/s12864-015-1527-z ↩

2. Massimo Bionaz, Kathiravan Periasamy, Sandra L. Rodriguez-Zas, et al. A Novel Dynamic Impact Approach (DIA) for Functional Analysis of Time-Course Omics Studies: Validation Using the Bovine Mammary Transcriptome. PLoS ONE, 2012. https://doi.org/10.1371/journal.pone.0032455 ↩

3. Jennifer T Aguilan, Katarzyna Kulej, Simone Sidoli. Guide for protein fold change and p -value calculation for non-experts in proteomics. Molecular Omics, 2020. https://doi.org/10.1039/d0mo00087f ↩

4. M.C. Verheijen, T.W. Gant, W. Tong, et al. R-ODAF: Omics data analysis framework for regulatory application. Toxicology Letters, 2021. https://doi.org/10.1016/s0378-4274(21)00539-7 ↩