Is log2(x+1) ever acceptable for RNA-seq data?

Yes, for quick visual QC on a single sample's count distribution or for a heatmap restricted to highly expressed genes, where the pseudocount's contribution is negligible relative to the actual counts. It's never acceptable as input to DESeq2, edgeR, or limma's statistical tests, and it's risky for fold changes on lowly expressed genes.

Should I use VST or rlog for my DESeq2 PCA?

Use VST for medium-to-large datasets (roughly n > 30 samples) where speed matters and the dispersion trend is well-behaved. Use rlog for smaller datasets, especially when sequencing depth varies widely across samples, since rlog handles that variability more robustly at the cost of speed.

Can I feed TPM or CPM values into DESeq2 for differential expression?

No. DESeq2 and edgeR expect raw, un-normalized counts because their internal normalization (median-of-ratios or TMM) and negative binomial dispersion estimation are calibrated for count data. Pre-normalized values like TPM or CPM break that model and will give you wrong p-values even though the functions won't throw an error.

Why does changing the pseudocount change my fold change so much?

Because the pseudocount is an additive offset, and for genes with very low counts the offset is a large fraction of the value being log-transformed. A gene with counts of 1 and 4 gives a very different log2 ratio depending on whether you add 1 or 0.1 to each value, even though nothing about the biology changed.

What should I use instead of a raw log2FoldChange for low-count genes?

Run DESeq2's lfcShrink() with the apeglm method. It uses an adaptive prior to pull fold-change estimates for lowly expressed or high-variance genes toward zero, giving you a more stable, less noise-driven number than the raw maximum-likelihood estimate from results().

Sanity check · Bulk RNA-seq

How to Log-Transform Counts Without Fooling Yourself in Bulk RNA-seq

The pseudocount you pick for low-count genes can move a fold change more than the biology does, and DESeq2's own statistical model never wants a log at all.

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

You've got a count matrix from STAR or salmon, a DESeq2 or edgeR object built, and before you run the actual test you want to eyeball a PCA or a heatmap. Log of zero is undefined, so you type log2(counts+1) without thinking about it, make your plot, and move on. Then, days later, you or a collaborator computes a fold change for a lowly expressed gene by hand, gets a number that looks dramatic, and it makes the top of the results table.

That pseudocount is not a neutral bookkeeping trick. For a gene with counts of 2 and 8 across two conditions, log2(x+1) gives you a very different answer than log2(x+0.1) would, and neither number has anything to do with the underlying biology. The distortion is worst exactly where you're least equipped to sanity-check it by eye: genes near the detection limit, where noise and signal look identical on a scatter plot.

This page gives you the ordered checks to run before you trust any log-transformed value, or before you let one anywhere near a p-value. In the next hour you can confirm whether your DESeq2/edgeR input is actually raw counts, whether your PCA is being driven by a pseudocount artifact instead of condition, and whether the "top hit" in your low-count genes survives shrinkage.

What it looks like when it's happening

A PCA built from log2(counts+1) separates samples by total library size rather than by condition, even though colSums() shows depth differs by 2x or more between groups.
The log2FoldChange for a low-count gene (single digit to low double digit counts) flips sign or changes by several fold when you rerun the same comparison with a different pseudocount.
A heatmap of your top DE genes is dominated by a handful of near-zero-count genes that look like solid blocks of color in one condition and near-black in the other.
DESeqDataSetFromMatrix() and DESeq() run without any error even though the input matrix contains non-integer or already-normalized values.
plotDispEsts(dds) shows dispersion estimates that don't shrink toward the fitted red trend line, or a trend line that looks flat or badly fit.
Two people computing 'the same' fold change for a lowly expressed gene get numbers 3-5x apart with no difference in the underlying counts, only in pseudocount choice.
An MA plot shows a cluster of enormous fold changes (over 100x) sitting at very low mean expression, with nothing comparable at higher expression levels.
Raw MLE log2FoldChange values collapse toward zero after lfcShrink() with apeglm, and the collapse is concentrated entirely in low-count genes.

Why it happens

A pseudocount is an arbitrary offset you add so log(0) doesn't blow up: log2(x+1), log2(x+0.1), log1p(x). For a gene sitting at hundreds or thousands of counts, the choice of offset barely matters, the ratio is dominated by the real counts. For a gene sitting at 0-10 counts, the offset dominates the ratio. Add 1 to counts of 1 and 4 and you get a log2 ratio of log2(5/2) = 1.32. Add 0.1 instead and you get log2(4.1/1.1) = 1.9. Same biology, same counts, a 44% swing in the reported fold change, purely from a number you picked without thinking about it.

DESeq2 and edgeR model raw counts with a negative binomial GLM: they estimate a mean-variance relationship (the dispersion trend) directly from the counts, then use median-of-ratios (DESeq2) or TMM (edgeR) normalization internally during the actual statistical test. Log-transforming counts before handing them to these tools breaks the count-based likelihood the model depends on. The dispersion estimation, the size-factor calculation, and the resulting p-values are all calibrated for count data, not for log-ratios with an arbitrary offset baked in. Feeding transformed values in doesn't produce an error, it produces wrong answers that look plausible.

Sequencing depth makes this worse, not just coincidentally. A sample sequenced to 40M reads detects more low-expression transcripts than a sample sequenced to 20M reads simply by sampling more molecules, independent of biology. If depth correlates with condition, even loosely, genes near the detection floor look "up-regulated" in the deeper samples for no reason except better sampling. Pseudocount-based fold changes amplify exactly this artifact, because the genes most sensitive to pseudocount choice are the same genes most sensitive to depth-driven detection noise.

VST and rlog exist to solve the problem correctly instead of papering over it. Both use the dispersion trend DESeq2 already estimated from your data to shrink low-count genes toward the gene's own average across samples, weighting by how much noise is expected at that count level. That is a principled variance-stabilizing correction, not a flat offset applied uniformly regardless of expression level. It is why VST/rlog output is safe for PCA, clustering, and heatmaps, and why a naive log2(x+1) is not.

The checks

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

0/7 checked · saved in this browser

01Grep your script for a transform applied before the model is built
Search your analysis script for log2, log1p, or cpm() calls that happen before DESeqDataSetFromMatrix() or DGEList(), or that feed into either.
bash
```
grep -n "log2\|log1p\|cpm(" analysis.R | grep -B2 -A2 "DESeqDataSetFromMatrix\|DGEList"
```
Healthy
No log/cpm transform appears upstream of DESeqDataSetFromMatrix() or DGEList(); those functions take the raw integer count matrix directly.
Red flag
A log2(counts+1) or cpm(log=TRUE) matrix is assigned to the same variable name later passed as countData or counts= into the model constructor.
02Confirm the model actually received raw integer counts
Check that the count matrix stored inside the DESeq2 object is whole-number valued, the way STAR/salmon output should be, not already normalized.
r
```
all(counts(dds) == round(counts(dds)))
```
Healthy
Returns TRUE, confirming every value in the matrix is a whole number consistent with raw sequencing counts.
Red flag
Returns FALSE, or throws an error, meaning the matrix contains fractional values, i.e. it was already normalized (TPM, CPM) or log-transformed before reaching DESeq2.
03Check whether sequencing depth is balanced across conditions
Sum the raw counts per sample and compare the totals across your condition groups to spot a depth imbalance that correlates with treatment.
r
```
colSums(counts(dds))
```
Healthy
Library sizes vary by less than roughly 2x across all samples, with no systematic difference between condition groups.
Red flag
One condition group is consistently 2-3x deeper than the other. That imbalance alone will inflate apparent up-regulation of low-count genes, independent of any pseudocount issue.
04Test pseudocount sensitivity directly on a handful of low-count genes
Pick a gene near your detection floor and recompute its log2 fold change with two different pseudocounts, then repeat for a highly expressed gene as a control.
r
```
x <- counts(dds)["LOW_COUNT_GENE", cond1_samples]
y <- counts(dds)["LOW_COUNT_GENE", cond2_samples]
log2((mean(y)+1)/(mean(x)+1))
log2((mean(y)+0.1)/(mean(x)+0.1))
```
Healthy
For the highly expressed control gene, the two pseudocounts give nearly identical log2FC.
Red flag
For the low-count gene, the two pseudocounts give values that differ by several-fold or flip sign.
05Inspect the dispersion fit
Plot per-gene dispersion estimates against the fitted trend after running DESeq().
r
```
plotDispEsts(dds)
```
Healthy
Gene-wise estimates (black) scatter around a smooth, monotonically decreasing red trend line and shrink toward it (blue).
Red flag
Estimates scatter wildly with no visible trend, or the trend line is flat or fails to fit. This is a strong sign the input to DESeq2 was not raw counts.
06Compare VST and rlog on your own data before picking one for plots
Run both transforms and compare the resulting PCA, especially if you have fewer than 30 samples and depth varies noticeably across them.
r
```
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
plotPCA(vsd)
plotPCA(rld)
```
Healthy
Sample clustering is broadly consistent between VST and rlog.
Red flag
VST's PCA looks noticeably different from rlog's, driven by a few low-count genes, in a dataset with small n and uneven depth. That divergence is a signal to prefer rlog for this dataset, per the DESeq2 vignette's own guidance.
07Compare raw MLE fold changes against shrunken estimates
Pull DESeq2's default results() alongside lfcShrink() with apeglm, and plot one against the other.
r
```
res <- results(dds, name="condition_treated_vs_control")
res_shrunk <- lfcShrink(dds, coef="condition_treated_vs_control", type="apeglm")
plot(res$log2FoldChange, res_shrunk$log2FoldChange)
```
Healthy
Points for well-expressed genes fall close to the y=x line; low-count genes shrink modestly toward zero.
Red flag
A handful of low-count genes have raw log2FC beyond ±5 that collapse toward zero after shrinkage. Those are pseudocount/noise artifacts, not biology, and they should not lead your results table.

What to do about it

Feed the model raw counts, always

When: Any time you're about to run DESeq(), edgeR's exactTest/glmQLFit, or build a DGEList.

Pass the untransformed integer count matrix into DESeqDataSetFromMatrix() or DGEList(). Let DESeq2's median-of-ratios or edgeR's TMM normalization happen internally during the call, not on a matrix you've already transformed.

Caveat: None functionally, this is just correct usage, but it means keeping a clean, untouched copy of the raw matrix separate from anything you transform for plotting.

Use VST or rlog for any visualization or distance-based step

When: Building a PCA, a sample-distance heatmap, or a clustering dendrogram, never for the statistical test itself.

Run vst(dds, blind=FALSE) for n > 30 samples where speed matters, or rlog(dds, blind=FALSE) for smaller datasets (n < 30) where sequencing depth varies widely across samples. Use the resulting matrix only for plots.

Caveat: rlog is noticeably slower on large sample counts; VST assumes a reasonably well-behaved dispersion trend and can behave oddly if that trend is poorly fit.

Shrink fold changes for low-count genes before reporting them

When: Any time you're presenting or ranking fold changes that include lowly expressed genes.

Run lfcShrink(dds, coef=..., type='apeglm') and report the shrunken log2FoldChange instead of the raw MLE estimate from results().

Caveat: Shrinkage biases the estimate toward zero by design. If you need an unbiased effect-size estimate for downstream meta-analysis or power calculations, use the raw MLE value and say so explicitly.

Use voom's log2-CPM with precision weights for the limma path

When: You're running limma instead of DESeq2/edgeR for the statistical test.

Build a DGEList, run calcNormFactors(), then voom(dge, design) to get log2-CPM values with precision weights that limma's linear model machinery expects. Don't substitute a manual log2(cpm(x)+1) for this.

Caveat: voom's weights are computed against the design matrix you supply; get the design wrong and the weights are wrong too, silently.

For pure visualization CPM checks in edgeR, use the built-in log option

When: You want a quick log2-CPM view for a QC plot in an edgeR workflow.

Call cpm(x, log=TRUE) on the DGEList, which already accounts for library size and any TMM factors from calcNormFactors().

Caveat: Still not meant as input to edgeR's statistical tests, only for plots and eyeballing.

When not to "fix" it

If you're just eyeballing a single sample's raw count distribution before you've built any DESeq2 object, checking for a failed library or an obviously bimodal sample, log2(x+1) on raw counts is fine. You're not comparing conditions or reporting a fold change, so pseudocount bias has nothing to distort. The same applies to a quick heatmap restricted to your top 50 highly expressed genes: at counts in the thousands, log2(x+1) and VST return nearly identical values, so the "wrong" transform doesn't move the picture. The problem only matters when you (a) feed transformed values into a count-based statistical test, or (b) compute or report fold changes for low-count genes and treat the pseudocount's contribution as if it were signal.

Five things experienced analysts do here

Keep two matrices with clearly different variable names in every script: one raw-counts object that only ever goes into DESeq()/DGEList(), and one vst/rlog/voom object that only ever goes into plots. Naming collisions between the two are how pseudocount bugs sneak into results tables.
Before calling DESeq(), run `all(counts(dds) == round(counts(dds)))` as a reflex, the same way you'd check a file exists before reading it. It catches TPM or pre-logged input in one line.
Never trust a raw log2FoldChange for a gene near your detection floor until you've compared it against lfcShrink's apeglm estimate. If the two disagree by a lot, the raw number is noise dressed up as an effect size.
Pick VST vs rlog based on your actual sample size and depth spread, not habit. Check both once per project on your own PCA before deciding, rather than defaulting to whichever one a tutorial used.
When someone hands you a fold change for a low-count gene without saying which pseudocount or which tool produced it, ask. The number is meaningless without that context, and it's usually the first thing that turns out to be wrong.

Questions people ask

Is log2(x+1) ever acceptable for RNA-seq data?: Yes, for quick visual QC on a single sample's count distribution or for a heatmap restricted to highly expressed genes, where the pseudocount's contribution is negligible relative to the actual counts. It's never acceptable as input to DESeq2, edgeR, or limma's statistical tests, and it's risky for fold changes on lowly expressed genes.
Should I use VST or rlog for my DESeq2 PCA?: Use VST for medium-to-large datasets (roughly n > 30 samples) where speed matters and the dispersion trend is well-behaved. Use rlog for smaller datasets, especially when sequencing depth varies widely across samples, since rlog handles that variability more robustly at the cost of speed.
Can I feed TPM or CPM values into DESeq2 for differential expression?: No. DESeq2 and edgeR expect raw, un-normalized counts because their internal normalization (median-of-ratios or TMM) and negative binomial dispersion estimation are calibrated for count data. Pre-normalized values like TPM or CPM break that model and will give you wrong p-values even though the functions won't throw an error.
Why does changing the pseudocount change my fold change so much?: Because the pseudocount is an additive offset, and for genes with very low counts the offset is a large fraction of the value being log-transformed. A gene with counts of 1 and 4 gives a very different log2 ratio depending on whether you add 1 or 0.1 to each value, even though nothing about the biology changed.
What should I use instead of a raw log2FoldChange for low-count genes?: Run DESeq2's lfcShrink() with the apeglm method. It uses an adaptive prior to pull fold-change estimates for lowly expressed or high-variance genes toward zero, giving you a more stable, less noise-driven number than the raw maximum-likelihood estimate from results().

Guide · How to Detect Batch Effects in Bulk RNA-seq
Compare · DESeq2 vs edgeR: Which One Should You Use?
Glossary · Log fold change shrinkage
Glossary · Log fold change (log2FC)

Sources

Analyzing RNA-seq data with DESeq2 — VST vs rlog recommendations, pseudocount and log-transform guidance, lfcShrink/apeglm usage
RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR — voom() transformation and cpm(log=TRUE) usage for visualization vs testing
Do you really understand log2Fold change in single-cell RNAseq data? — Pseudocount bias demonstration and why different pseudocounts give incomparable fold changes
Benchmarking differential expression analysis tools for RNA-Seq — Systematic comparison of normalization-based vs log-ratio transformation-based DE methods

Part of the Log transformation of counts series.