/pydeseq2-differential-expression

name: "pydeseq2-differential-expression"
description: "Bulk RNA-seq DE with PyDESeq2: load counts, normalize, fit negative binomial models, Wald test (BH-FDR), LFC shrinkage, volcano/MA plots. Use for two-group comparisons, multi-factor designs with batch correction, multiple contrasts."
license: "CC-BY-4.0"

PyDESeq2 Differential Expression Analysis

Overview

PyDESeq2 is a Python reimplementation of the R DESeq2 package for differential gene expression analysis from bulk RNA-seq count data. It fits negative binomial generalized linear models per gene, estimates dispersion with empirical Bayes shrinkage, and performs Wald tests with Benjamini-Hochberg FDR correction. This skill covers the full pipeline from raw counts to publication-ready result tables and visualizations.

When to Use

• Identifying differentially expressed genes between two or more experimental conditions from bulk RNA-seq
• Performing two-group comparisons (e.g., treated vs control) with proper statistical testing
• Running multi-factor designs that account for batch effects or covariates (e.g., `~batch + condition`)
• Applying log2 fold change shrinkage (apeGLM) for ranking and visualization
• Use `omics-plotting` SKILL after DE for publication-quality plots of differential expression results
• Converting R-based DESeq2 workflows to a pure Python environment
• Integrating DE analysis into larger Python bioinformatics pipelines (e.g., with scanpy, pandas)
• Use **DESeq2** (R/Bioconductor) or **edgeR** instead for the reference R implementations with the broadest method support and community validation

Prerequisites

• **Python packages**: `pydeseq2>=0.4`, `pandas>=1.4`, `numpy>=1.23`, `scipy>=1.11`, `scikit-learn>=1.1`, `anndata>=0.8`
• **Data requirements**: Raw (unnormalized) integer count matrix (samples x genes) + sample metadata DataFrame
• **Environment**: Python 3.10+; optional `matplotlib`, `seaborn` for visualization

pip install pydeseq2 matplotlib seaborn

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: design
    kind: required
    source: data
    ask: "Which column of the sample sheet separates the groups you want to compare?"
    default: null

  - id: D2
    param: contrast
    kind: required
    source: user
    depends_on: [D1]
    ask: "Within that column, which group is the baseline the others are measured against?"
    default: null

  - id: D3
    param: design
    kind: required
    source: data
    depends_on: [D1]
    ask: "Is there a nuisance variable to adjust for (batch, donor, sex, sequencing run)?"
    default: "none - compare groups without covariates"
    skip_if: "sample sheet carries no column besides the grouping variable"

  - id: D4
    param: alpha
    kind: optional
    source: user
    ask: "How strong must the evidence be before a gene counts as changed (false-discovery rate)?"
    default: 0.05

  - id: D5
    param: lfc_cutoff
    kind: optional
    source: user
    ask: "How large a fold change is worth reporting, on top of the statistical call?"
    default: "|log2FC| > 0.58 (1.5x change)"

  - id: D6
    param: min_total_counts
    kind: optional_conditional
    source: data
    ask: "Genes with very few reads cannot be tested - where should the floor sit?"
    default: 10

  - id: D7
    param: lfc_shrink
    kind: optional
    source: user
    ask: "Shrink fold-change estimates so low-count genes stop dominating the rankings and plots?"
    default: "applied for ranking and visualization, not for significance calls"

  - id: D8
    param: refit_cooks, cooks_filter, independent_filter
    kind: optional
    source: user
    ask: "Should DESeq2 retain its default outlier replacement and adaptive filtering behavior?"
    default: true

  - id: D9
    param: n_cpus
    kind: never_ask
    source: data
    reason: "Affects runtime only, not the result"
    default: 4

`design` appears twice because the formula is assembled from both D1 and D3; D2 then selects which coefficient of that fitted model is tested.

Workflow

Step 1: Data Loading and Validation

Load the count matrix and metadata. PyDESeq2 expects counts as a samples x genes DataFrame with non-negative integers, and metadata as a samples x variables DataFrame with matching indices.

import pandas as pd

# Load data — typical CSV has genes as rows, samples as columns
counts_raw = pd.read_csv("counts.csv", index_col=0)
metadata = pd.read_csv("metadata.csv", index_col=0)

# Transpose if needed: PyDESeq2 requires samples x genes
if counts_raw.shape[0] > counts_raw.shape[1]:
    counts_df = counts_raw.T  # genes x samples → samples x genes
else:
    counts_df = counts_raw

# Validate alignment
common_samples = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common_samples]
metadata = metadata.loc[common_samples]

print(f"Samples: {counts_df.shape[0]}, Genes: {counts_df.shape[1]}")
print(f"Metadata columns: {list(metadata.columns)}")
print(f"Condition counts:\n{metadata['condition'].value_counts()}")

Step 2: Gene Filtering

Remove lowly expressed genes to improve statistical power and reduce multiple testing burden.

# Filter genes with total counts below threshold
min_total_counts = 10
gene_counts = counts_df.sum(axis=0)
genes_to_keep = gene_counts[gene_counts >= min_total_counts].index
counts_df = counts_df[genes_to_keep]

# Optional: require minimum counts in a minimum number of samples
min_count_per_sample = 5
min_samples = 3
genes_expressed = (counts_df >= min_count_per_sample).sum(axis=0) >= min_samples
counts_df = counts_df.loc[:, genes_expressed]

print(f"Genes after filtering: {counts_df.shape[1]}")

Step 3: DeseqDataSet Initialization and Fitting

Create the DESeq dataset object, specify the design formula, and run the full pipeline (size factor estimation, dispersion estimation, model fitting).

from pydeseq2.dds import D