PyDESeq2: Differential Expression in Python

PyDESeq2 is a Python package for finding genes whose expression differs between conditions in bulk RNA-seq. It takes raw counts, the number of reads assigned to each gene, and sample information describing the experimental conditions, or metadata. It models sequencing depth and biological variation, then calculates the magnitude of each condition effect and its statistical evidence.

For example, if every experimental subject has one Control and one Treatment sample, PyDESeq2 models the changes repeated across subjects rather than directly dividing two counts. Its output is a table with columns such as log2FoldChange, pvalue, and padj for every gene.

How is it different from DESeq2?

DESeq2 is an R/Bioconductor package for bulk RNA-seq differential expression analysis. PyDESeq2 is a reimplementation of that method for Python users. It is not a wrapper that calls R or DESeq2 internally.

The two tools share the broad workflow of normalization, dispersion estimation, negative-binomial regression, and Wald testing. The official PyDESeq2 documentation notes, however, that values may differ slightly from R DESeq2 because of differences in the reimplementation and supported features. Record the version as well as the package name to make an analysis reproducible.

Pin the version in a reproducible analysis project. The example below uses 0.5.4.

uv add "pydeseq2==0.5.4"

The inputs are two tables

PyDESeq2 requires a count matrix and metadata. Both can be passed as pandas DataFrames.

1. Count matrix

sample	gene_A	gene_B	gene_C
S1_Control	120	8	520
S1_Treatment	210	5	490
S2_Control	980	12	430
S2_Treatment	1,150	7	470

• Each row is one sample.
• Each column is one gene.
• Values are non-negative integer raw counts.

A table such as an HTSeq file, with genes in rows and samples in columns, must be transposed. Separate non-gene summary counters such as __no_feature before providing the input. HTSeq raw counts explains how raw counts are produced.

2. Metadata

sample	subject	condition
S1_Control	S1	Control
S1_Treatment	S1	Treatment
S2_Control	S2	Control
S2_Treatment	S2	Treatment

The row names in the metadata must exactly match the sample names in the count matrix. subject and condition are sample attributes used by the statistical model. Store subject IDs as categorical variables that distinguish subjects, not as numeric measurements whose magnitude is compared.

The analysis design defines the comparison question

design = "~ subject + condition"

This formula tells the model to estimate the condition effect after accounting for each subject’s expression baseline. Whether to include only condition or also include batch or subject depends on the experimental structure and analysis question. Design formulas explains how to read ~ and +.

Specify the direction of the difference with a contrast.

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

This contrast means Treatment - Control. A positive log2FoldChange therefore means expression is higher in Treatment, and a negative value means it is lower in Treatment.

What PyDESeq2 does internally

PyDESeq2 proceeds in the following order.

1. Estimate size factors: correct for differences in sequencing depth among samples.
2. Estimate dispersion: for every gene, estimate how much counts vary among biological replicates in the same condition.
3. Fit the model: model raw counts with a negative-binomial distribution and estimate each effect in the design.
4. Wald test: test whether the condition effect specified by the contrast is distinguishable from zero.
5. Multiple-testing correction: correct p-values across thousands of genes to calculate padj.

This is not a workflow that first turns raw counts into one normalized table and then repeats a simple t-test. Normalization, subject and condition effects, and gene-specific variation are handled together in one count model.

Two objects divide the work in code

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~ subject + condition",
    n_cpus=1,
)
dds.deseq2()

stats = DeseqStats(
    dds,
    contrast=["condition", "Treatment", "Control"],
)
stats.summary()

results = stats.results_df

Object	Role
DeseqDataSet	Takes counts, metadata, and a design, then fits size factors, dispersion, and fold changes
DeseqStats	Runs the Wald test and multiple-testing correction for a specified contrast and produces the result table

If you create DeseqStats from an object before dds.deseq2() finishes, the required model estimates are incomplete. Conversely, one fitted DeseqDataSet can produce several comparisons by specifying different contrasts.

Reading result-table columns

stats.results_df contains the following values for each gene.

Column	Meaning
baseMean	Mean count across all samples after size-factor normalization
log2FoldChange	Log2 change showing how much higher the tested level is than the reference level
lfcSE	Standard error of the estimated log2 fold change
stat	Wald test statistic
pvalue	P-value evaluating the result under a model with no condition effect
padj	P-value corrected for testing many genes at once

A log2FoldChange of 1 means the tested level is about twice as high, while -1 means it is about half as high. The direction changes with the order of the contrast, so always record the comparison name with the result.

Do not select candidates from padj alone. Examine the magnitude, log2FoldChange; uncertainty, lfcSE; mean expression, baseMean; and the per-sample pattern together. Reading DESeq2 Results explains how the six columns are calculated and connected, while Reading a Volcano Plot shows how to narrow candidates on two axes. Statistical Testing and Multiple Testing explains the general meaning of p-values and multiple testing.

Values for DEG testing differ from values for visualization

Differential expression testing takes integer raw counts. PCA and heatmaps use transformed values such as VST to moderate the variance of large counts.

Purpose	Value to use
Fit the DEG model	Integer raw counts
PCA, sample distances, and heatmaps	VST or another visualization transform
Interpret the result	log2FoldChange, padj, and per-sample expression patterns

Do not rerun DEG testing on VST values or provide TPM as PyDESeq2 count input. Statistical testing and visualization require different properties from values derived from the same data.

Checks before running

• Does the count matrix have samples in rows and genes in columns?
• Are all counts non-negative integers?
• Do the sample names and order match between counts and metadata?
• Does every column named in the design exist in the metadata?
• Are identifiers such as subject IDs and batch IDs categorical?
• Does each comparison group have biological replicates?
• Did you filter low-expression genes using a defined criterion?
• Is the order of the tested and reference levels in the contrast intentional?
• Did you record the PyDESeq2 version and filtering criteria?

PyDESeq2 cannot judge the biological meaning of incorrectly constructed metadata. A computation can succeed yet accurately answer the wrong question if its design does not match the study question.

Official resources

• Official PyDESeq2 0.5.4 documentation
• Official PyDESeq2 basic workflow
• Official PyDESeq2 GitHub repository
• PyDESeq2 paper