Performing bulk RNA-seq differential expression analysis allows researchers to uncover transcriptomic alterations across biological conditions, disease states, or perturbations. While R's DESeq2 has historically served as the standard framework for negative binomial modeling of read counts, Python users frequently required bridge packages or language switching to complete downstream workflows. PyDESeq2 offers a native Python re-implementation of the DESeq2 pipeline, fitting generalized linear models directly within the scientific Python stack.

In this guide, you will learn how to prepare read count matrices, configure experimental design formulas, execute differential testing using PyDESeq2, and validate output files for reproducibility. Because this tutorial is designed as an agent-compatible planning guide rather than a benchmarked execution report, you should treat each step as a structured protocol to adapt to your specific dataset and computational environment.

---

Upstream Prerequisites and Infrastructure

Differential expression analysis does not operate directly on raw FASTQ sequencing reads. Before running statistical tests, your pipeline must generate an unnormalized gene-level count matrix. Upstream workflows handled by tools referenced in the bulk-rnaseq skill typically manage read alignment with STAR or pseudo-alignment with Salmon and Kallisto. Additionally, inspecting genome coverage tracks or signal distributions using deeptools helps verify library balance and track consistency across genomic regions.

Reproducibility is paramount in computational biology. Variations in underlying shared libraries can alter statistical outcomes or numerical convergence. When setting up your runtime environment, apply these two operational principles:

  1. Use containerized environments: Docker packages software and its operational dependencies into portable, isolated images. On shared high-performance computing (HPC) clusters where users lack administrative rights, Singularity/Apptainer serves as the standard alternative because it runs without root privileges.
  2. Pin code states precisely: Recording the exact git commit hash via git rev-parse HEAD pins the analytical code version far more precisely than relying on mutable branch names or arbitrary release tags.

If you are orchestrating this workflow with autonomous assistants, note that Agent Skills follows an open format in which each skill is a folder containing a SKILL.md file with YAML frontmatter containing name and description fields, along with optional execution scripts and reference files. AI agents discover analytical skills by reading the SKILL.md file and loading referenced resources on demand.

---

Step-by-Step Bulk RNA-seq Differential Expression with PyDESeq2

The pydeseq2 package provides clean object-oriented classes—principally DeseqDataSet and DeseqStats—to estimate size factors, fit gene-wise dispersions, and compute hypothesis tests.

Step 1: Format Count Matrices and Metadata

PyDESeq2 expects two distinct tabular structures:

  • Counts table: A non-negative integer matrix where rows represent individual samples and columns represent genes or transcript identifiers. Do not supply normalized counts (such as TPM, FPKM, or CPM), because the negative binomial model explicitly models integer sampling variability.
  • Metadata table: A table where rows match the count matrix samples, and columns define experimental factors, such as genotype, treatment, batch, or sequencing run.
import numpy as np
import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# Load raw count matrix (samples as rows, genes as columns)
counts_df = pd.read_csv("gene_counts.csv", index_col=0)

# Load sample metadata (samples as rows, covariates as columns)
metadata = pd.read_csv("sample_metadata.csv", index_col=0)

# Verify that sample identifiers align precisely
assert (counts_df.index == metadata.index).all(), "Sample IDs do not match."

# Filter unexpressed or exceedingly low-count genes to reduce testing burden
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_filtered = counts_df[genes_to_keep]

Step 2: Initialize and Fit the DeseqDataSet

Instantiate a DeseqDataSet by passing the filtered count matrix, the sample metadata, and the target design factor. Calling .deseq2() triggers dispersion parameter fitting, size factor estimation, and generalized linear model regression.

# Initialize DeseqDataSet with the primary experimental factor
dds = DeseqDataSet(
    counts=counts_filtered,
    metadata=metadata,
    design_factors="condition",
    refit_cooks=True,
    n_cpus=4
)

# Run estimation of size factors, dispersions, and dispersion shrinkage
dds.deseq2()

Step 3: Run Statistical Tests and Extract Contrast Results

Once the base model has converged, evaluate differential expression between experimental levels using the DeseqStats class. Explicitly defining contrasts prevents baseline ambiguities.

# Construct statistical contrast: factor name, test condition, reference condition
stat_res = DeseqStats(
    dds, 
    contrast=["condition", "treated", "control"], 
    n_cpus=4
)

# Run Wald testing and independent filtering
stat_res.summary()

# Extract the final results DataFrame
results_df = stat_res.results_df
results_df.to_csv("differential_expression_results.csv")

The resulting table reports base mean counts, log2 fold changes, standard errors, Wald test statistics, raw p-values, and Benjamini-Hochberg adjusted p-values (padj).

Step 4: Transition to Downstream Biological Interpretation

After identifying candidate differentially expressed genes, export prioritized gene lists filtered by significance thresholds (such as adjusted p-values below defined cutoffs) into downstream functional annotation workflows. Integrating the pathway-enrichment skill enables over-representation analysis (ORA) or gene set enrichment analysis (GSEA) to map gene expression changes to functional pathways.

---

Checks Before You Trust the Result

Never evaluate differential expression results purely on raw numerical output. Before adopting statistical findings, verify the following diagnostic indicators:

  • Independent Filtering and Dispersion Convergence: Check whether the gene-wise dispersion estimates properly shrink toward the fitted dispersion-mean curve. High-variance outliers should be flagged or refit via Cook's distance diagnostics.
  • P-value Uniformity Under the Null: Plot a histogram of raw p-values. A healthy distribution typically features an elevated peak near zero representing true differential signal, superimposed on a flat, uniform distribution across the interval from zero to one. A severely skewed or U-shaped distribution suggests unmodeled batch effects, sample swap errors, or model misspecification.
  • Reference Factor Verification: Verify that the sign of log2 fold changes aligns with the biological question. If the contrast was inverted (for example, comparing control against treated instead of treated against control), the reported directions will be reversed.
  • Sample Clustering and Outlier Inspection: Compute principal component analysis (PCA) plots on variance-stabilized transformations of the count data. Replicates should cluster predominantly by biological condition rather than by sequencing run, extraction date, or library preparation batch.
  • Execution Provenance: Confirm that your run logged the container runtime specification and recorded the Git commit hash using git rev-parse HEAD. This ensures that downstream collaborators can re-execute the identical software environment.

---

Running statistical differential testing is an iterative process that relies heavily on proper data hygiene and upstream quality checks. Before executing these steps in automated pipelines or deployment clusters, review the respective source documentation and parameter schemas linked across our pydeseq2, bulk-rnaseq, deeptools, and pathway-enrichment skill entries on ScienceSkills.