Reading DESeq2 Results

A differential expression result table from DESeq2 or PyDESeq2 is a statistical result that calculates the expression difference between two conditions and its uncertainty for each gene. From integer raw counts and sample information, it outputs one row per gene with six values: baseMean, log2FoldChange, lfcSE, stat, pvalue, and padj.

This document explains results from the basic Wald test. The examples use the following direction, comparing Treatment with Control.

contrast = ["condition", "Treatment", "Control"]

A positive log2FoldChange therefore means expression is higher in Treatment, while a negative value means it is higher in Control. Reversing the comparison reverses the sign.

Start by reading one result row

The following values are a hypothetical result created for explanation.

gene	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
gene_A	500	2.00	0.50	4.00	$6.3 \times 10^{-5}$	0.003

Read the row as follows:

After normalization, gene_A was observed at an average count of about 500 across all samples. Its expression is estimated to be about four times higher in Treatment than Control, and the standard error of the change estimate 2.00 is 0.50. Dividing the change by its standard error gives a Wald statistic of 4.00. The p-value for this individual gene is about 0.000063, and the p-value corrected for testing all genes together is 0.003.

The six values are not an unrelated list of names. Grouping their relationships makes the table easier to read.

Role	Column	Relationship to the next value
Summary of overall observation level	baseMean	Not the direct starting value used to calculate the other five columns
Magnitude of the condition effect	log2FoldChange	Used with lfcSE to calculate stat
Uncertainty in the effect estimate	lfcSE	Used with log2FoldChange to calculate stat
Ratio of effect to uncertainty	stat	Used to calculate pvalue
Test result for one gene	pvalue	Corrected together with values from all genes
Multiple-testing result	padj	Statistical evidence used in candidate selection

The key point is that baseMean is not the input from which log2FoldChange is calculated. Both come from the same count data but answer different questions.

baseMean: how much of this gene was observed overall?

baseMean is the mean normalized count across all samples. Because sequencing depth differs among samples, it does not average raw counts directly. It divides each value by that sample’s size factor first.

For raw count $K_{ij}$ of gene $i$ in sample $j$, sample size factor $s_j$, and total sample count $m$, it can be understood as:

$$\operatorname{baseMean}_i = \frac{1}{m} \sum_{j=1}^{m} \frac{K_{ij}}{s_j}$$

baseMean = 500 does not mean that one sample contained 500 reads. It means that after correcting sequencing depth across samples from all conditions, the average observation level was about 500.

baseMean alone does not reveal direction

Two genes can have the same baseMean and completely different condition effects. The table uses simple means for illustration.

Gene	Control mean	Treatment mean	Overall mean	Simple fold change
gene_A	50	50	50	1-fold
gene_B	10	90	50	9-fold

Both genes have an overall mean of 50. But gene_A does not differ between conditions, while gene_B is higher in Treatment. You cannot calculate log2FoldChange or infer direction from baseMean.

baseMean is useful for questions such as:

• Was this gene observed at a sufficient level overall?
• Did a large fold change arise from extremely low counts?
• Which expression range contains the gene on the horizontal axis of an MA plot?

Conversely, baseMean does not separately show the Treatment and Control means. Inspect normalized counts separately to see condition-level values and patient-level patterns.

log2FoldChange: in which direction and by what factor?

log2FoldChange expresses the condition effect of interest in log base 2 units. The sign gives the direction, and the absolute value gives the magnitude.

log2FoldChange	Treatment / Control	Interpretation
2	4	About four times higher in Treatment
1	2	About twice as high in Treatment
0	1	No difference
-1	1/2	Treatment is about half of Control
-2	1/4	Treatment is about one quarter of Control

Convert back to the original factor with:

$$\text{fold change}=2^{\text{log2FoldChange}}$$

The hypothetical log2FoldChange = 2 means expression is about four times higher in Treatment because $2^2=4$.

Why use log base 2?

Untransformed fold changes are asymmetric: a two-fold increase is 2, while a two-fold decrease is 0.5. In log base 2, equal increases and decreases become 1 and -1, symmetric around zero.

Logs also turn multiplication into addition.

$$\log_2(a \times b)=\log_2(a)+\log_2(b)$$

Within the statistical model, factors that act multiplicatively on expression, including sequencing depth, patient-specific baselines, and condition effects, can be handled as additive terms. Base 2 is not mathematically required, but it makes a value of 1 read directly as two-fold and is convenient in result tables.

It is not a direct ratio of two group means

With only two numbers, you could calculate log2 fold change as:

$$\log_2\left( \frac{\text{Treatment expression}} {\text{Control expression}} \right)$$

DESeq2 and PyDESeq2 do not use the following simple calculation:

log2((treatment_mean + 1) / (control_mean + 1))

For each gene, they put integer raw counts from all samples into a negative-binomial model and estimate the condition effect specified in the design formula. With this paired design:

design = "~ patient + condition"

the model allows a different baseline for every patient and estimates the condition effect repeated across patients. log2FoldChange is that model coefficient expressed in log base 2 units. baseMean is not an ingredient in that coefficient; it is a separately reported summary of observation level.

The negative-binomial model

Count $K_{ij}$ for gene $i$ in sample $j$ is modeled with a negative-binomial distribution having mean $\mu_{ij}$ and gene-specific dispersion $\alpha_i$.

$$K_{ij} \sim \operatorname{NB}(\mu_{ij}, \alpha_i)$$

The expected count can be separated into the sample size factor $s_j$ and the gene’s condition-dependent expression level $q_{ij}$.

$$\mu_{ij}=s_jq_{ij}$$

Taking a logarithm turns multiplication into addition.

$$\log(\mu_{ij}) = \log(s_j) + \beta_{i0} + x_{j1}\beta_{i1} + \cdots$$

$\log(s_j)$ accounts for sample-specific sequencing depth. Values such as $x_{j1}$ are columns of the design matrix created from metadata, and $\beta_{i1}$ is the effect of that column.

If condition has the two categories Control and Treatment, with Control as the reference category, the condition coefficient is the Treatment-versus-Control log fold change. Internal model calculations may use natural logarithms, but DESeq2 and PyDESeq2 report coefficients and standard errors in log base 2 units.

Check whether LFC shrinkage was applied

Genes with low counts or high dispersion can have unstable, very large log2 fold changes. LFC shrinkage is an optional step that moves change estimates for information-poor genes toward zero, stabilizing gene ranking and visualization.

In PyDESeq2 0.5.4, call lfc_shrink() separately.

stats.lfc_shrink(coeff="condition[T.Treatment]")

log2FoldChange and lfcSE may differ before and after this call. In the official example, the existing Wald-test stat, pvalue, and padj remain unchanged. Record whether shrinkage was applied and which coefficient was used.

lfcSE: how precisely was the log2 fold change estimated?

lfcSE is the log2 fold change standard error. If log2FoldChange is the estimated condition effect, lfcSE describes the uncertainty of that estimate.

Standard error is different from standard deviation.

Value	Question it answers
Standard deviation (SD)	How widely are individual sample values spread?
Standard error (SE)	How uncertain is the effect estimate calculated from those samples?

For a simple mean, $SE=SD/\sqrt{n}$, but lfcSE is not this formula applied directly to counts. It is the standard error of the condition coefficient calculated from the negative-binomial model, gene-specific dispersion, size factors, and design matrix together.

The following conditions can increase lfcSE:

• Few samples
• Large count variation within one condition
• Extremely low counts for the gene
• Many effects to estimate in the design
• Two effects that the data cannot separate well

Conversely, more replicates and consistent changes in the same direction tend to reduce the standard error.

Read it relative to the change

There is no universal standalone threshold for lfcSE. Do not interpret it as “good below 0.5”; compare it with log2FoldChange.

In the hypothetical result, log2FoldChange = 2.00 and lfcSE = 0.50. The change is four times its standard error.

$$2.00 \div 0.50=4.00$$

Under the Wald approximation, an approximate 95% interval is:

$$2.00 \pm 1.96 \times 0.50 = 1.02 \sim 2.98$$

Converted into fold change, the interval is about $2^{1.02}=2.0$-fold to $2^{2.98}=7.9$-fold. The exact magnitude is uncertain, but the entire interval is above zero, so the direction of increase in Treatment is relatively stable.

This interval illustrates the basic Wald approximation. If you used shrinkage or another test configuration, consult the inference procedure provided by that method.

stat: how many standard errors is the change from zero?

For the default Wald test, stat is calculated as:

$$\operatorname{stat} = \frac{\text{log2FoldChange}-0}{\text{lfcSE}}$$

The zero in the numerator is the null hypothesis of no condition effect, or log2FoldChange = 0. In the hypothetical result, 2.00 ÷ 0.50 = 4.00, so stat = 4.00.

• Large positive value: increased in Treatment and far from zero
• Large negative value: decreased in Treatment and far from zero
• Value near zero: estimated effect is small or uncertainty is large

A large absolute stat gives stronger statistical evidence that the data do not fit the no-difference hypothesis. It does not mean that the gene is biologically important. Even a small change can have a large stat when its standard error is very small.

The testing threshold also changes if lfc_null is not zero. For example, you can test |log2FoldChange| > 1 rather than merely testing for any difference. In that case, do not always read stat as log2FoldChange ÷ lfcSE.

stat means something different in a likelihood ratio test

The equation in this document applies to the default Wald test. With a likelihood ratio test (LRT), stat is the deviance difference between the reduced and full models, and its p-value is calculated from a chi-squared distribution.

pvalue: how inconsistent is the result with a no-difference model?

The null hypothesis of the default two-sided Wald test is:

$$H_0:\text{log2FoldChange}=0$$

pvalue is the probability, assuming this null hypothesis is true, of observing a result at least as extreme as the current stat. The Wald statistic is compared with a standard normal distribution, and the probabilities in both tails are added.

$$p = 2\left(1-\Phi\left(|\operatorname{stat}|\right)\right)$$

$\Phi$ is the cumulative distribution function of the standard normal distribution. The hypothetical stat = 4.00 gives a p-value of about $6.3 \times 10^{-5}$.

In a CSV file, 6.3E-05 is scientific notation for $6.3 \times 10^{-5}$. In decimal form, it is 0.000063.

What a p-value does not mean

pvalue = 0.01 does not mean any of the following:

• There is a 1% probability that the result occurred by chance.
• There is a 1% probability that the null hypothesis is true.
• There is a 99% probability that this gene is causal.
• The change is large or biologically important.

A p-value describes how extreme the data are under a particular statistical model with no condition effect. Separately verify that the model and design match the research question and that count quality and outliers are acceptable.

padj: correcting for testing thousands of genes

RNA-seq tests thousands of genes at once, not one gene. In a simple scenario where no gene truly differs, applying only pvalue < 0.05 to 20,000 genes would let about 1,000 pass by chance on average.

padj is an adjusted p-value that considers all p-values together to handle this accumulation of false positives. By default, the Benjamini-Hochberg procedure controls the false discovery rate (FDR). A false discovery in this setting is different from one FP observation in a classification confusion matrix.

FDR limits the expected proportion of false discoveries within the selected result set. padj < 0.05 does not mean that each gene has exactly a 5% chance of being wrong.

padj cannot be calculated from one row

pvalue is calculated from that gene’s stat, while padj depends on the ranks of p-values from all genes in the analysis. Even if one gene’s p-value stays the same, its padj can change when the tested gene set or filtering settings change.

Record the prefilter criterion and number of genes included in the analysis. Arbitrarily removing genes you dislike after seeing the result and recalculating padj lets the selection process itself affect the result.

Independent filtering may also be applied

DeseqStats in PyDESeq2 0.5.4 uses independent filtering by default. It excludes genes with such low baseMean that they have little testing power from multiple-testing correction, while finding a cutoff that can increase discoveries at the chosen alpha.

This can produce rows that have a pvalue but a NaN padj. That is not the same as padj = 1. It means the gene did not pass the independent filter, so no adjusted value was calculated.

When does NaN appear?

Missing values do not all have the same cause. Check PyDESeq2 settings and execution logs together.

Visible result	Possible cause	First check
Both pvalue and padj are NaN	Count outlier detected by Cook’s distance	Does one sample dominate the result?
pvalue exists but only padj is NaN	Low-baseMean row removed by independent filtering	Independent-filter settings and cutoff
All counts are zero	No information for estimating the effect or variance	Prefilter and original counts
Very large LFC and large SE	Unstable ratio caused by low counts in one condition	Counts by condition and sample

Before replacing all NaN values with zero or one, determine why they are missing. If a visualization treats missing padj as 1 for convenience, state that this is a display classification rather than a change to the statistical result.

An order for reading all six values

Do not select DEG candidates from one column alone. Use the following order to separate observation level, effect size, and statistical evidence.

1. Check comparison direction: What are the contrast’s tested and reference levels?
2. Check baseMean: Was the gene observed sufficiently, or is it extremely low?
3. Check log2FoldChange: In which direction and by what factor did it change?
4. Check lfcSE: How uncertain is the effect size?
5. Check stat: How large is the change relative to its standard error?
6. Check padj: Does statistical evidence remain after multiple testing?
7. Check per-sample counts: Is the result driven by a subset of samples or one patient?

Common combinations can be interpreted as follows.

Result combination	Interpretation
Large LFC, small SE, small padj	Large, stably estimated condition difference
Small LFC, small SE, small padj	Small but consistent condition difference across replicates
Large LFC, large SE, large padj	Change looks large but is uncertain
Low baseMean, large LFC	Check whether a difference of only a few reads produced a large ratio
Small padj, only a few extreme samples	Recheck outliers and model fit

Thresholds such as padj < 0.05 and |log2FoldChange| >= 1 can be used to select candidates, but they are not universal biological boundaries. Choose them for the research question, sample size, cost of follow-up validation, and required effect size.

What the result table cannot establish

The six values show the magnitude and statistical evidence of a condition difference. They do not directly prove that:

• The gene causes the disease.
• The expression change produces the same change in protein abundance or activity.
• A change in bulk tissue occurred within a particular cell type.
• The effect will reproduce in another dataset or population.
• The gene with the smallest padj is the most biologically important.

DEGs from bulk RNA-seq can reflect both regulatory changes within cells and changes in the tissue’s cellular composition. Recheck top candidates in per-patient normalized counts, raw-read quality, tissue information, and independent data.

What to record for reproducibility

To reproduce the same result table, save more than the CSV. Record:

• Exact DESeq2 or PyDESeq2 version
• Annotation and quantification-tool version used to produce raw counts
• Low-expression prefilter criterion
• Design formula
• Order of the contrast’s tested and reference levels
• Wald test or LRT selection
• alpha, lfc_null, and alt_hypothesis
• Independent-filtering and Cook’s-filtering settings
• Whether LFC shrinkage was applied and the coefficient used
• All samples used to produce the result and any excluded samples

PyDESeq2 explains the tool’s inputs and execution flow. Statistical testing and multiple testing explains the general meaning of p-values and FDR. Differential expression testing across several samples covers the full DEG workflow.

Official resources

• PyDESeq2 0.5.4 basic workflow
• PyDESeq2 0.5.4 DeseqStats API
• Official DESeq2 package manual
• Original DESeq2 paper