Pydeseq2 Differential Expression

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

bashCopy code
pip install pydeseq2 matplotlib seaborn

Pre-flight Interview

Settle these with the user before writing any analysis code.

yamlCopy 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.

pythonCopy code
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.

pythonCopy code
# 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).

pythonCopy code
from pydeseq2.dds import DeseqDataSet

dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~condition",   # Wilkinson-style formula
    refit_cooks=True,      # Refit after Cook's outlier removal
    n_cpus=4               # Parallel threads
)

# Run: size factors → dispersions → trend → MAP shrinkage → LFC fitting
dds.deseq2()

# Inspect normalization
print(f"Size factors (first 5): {dds.obsm['size_factors'][:5]}")
print(f"Size factor range: {dds.obsm['size_factors'].min():.2f} - {dds.obsm['size_factors'].max():.2f}")

Step 4: Statistical Testing (Wald Test)

Perform Wald tests to identify differentially expressed genes. Specify the contrast as [variable, test_level, reference_level].

pythonCopy code
from pydeseq2.ds import DeseqStats

ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],
    alpha=0.05,              # FDR threshold
    cooks_filter=True,       # Filter Cook's outliers
    independent_filter=True  # Independent filtering for power
)

ds.summary()

# Access full results
results = ds.results_df
print(f"Total genes tested: {len(results)}")
print(f"Significant (padj < 0.05): {(results.padj < 0.05).sum()}")

Step 5: LFC Shrinkage (Optional)

Apply apeGLM shrinkage to reduce noise in log2 fold change estimates. Use shrunk values for visualization and ranking, not for significance calls.

pythonCopy code
# Apply shrinkage — modifies results_df.log2FoldChange in place
ds.lfc_shrink()

# Compare pre/post shrinkage effect
print(f"Max |LFC| after shrinkage: {results.log2FoldChange.abs().max():.2f}")
print(f"Genes with |LFC| > 2: {(results.log2FoldChange.abs() > 2).sum()}")

Step 6: Result Filtering and Export

Filter significant genes and export results for downstream analysis.

pythonCopy code
import numpy as np

# Significance + effect size filter
significant = results[
    (results.padj < 0.05) &
    (results.log2FoldChange.abs() > 1.0)
].copy()

# Separate up/down-regulated
up = significant[significant.log2FoldChange > 0].sort_values("padj")
down = significant[significant.log2FoldChange < 0].sort_values("padj")

print(f"Upregulated: {len(up)}, Downregulated: {len(down)}")

# Export
results.to_csv("deseq2_all_results.csv")
significant.to_csv("deseq2_significant.csv")
up.to_csv("deseq2_upregulated.csv")
down.to_csv("deseq2_downregulated.csv")
print("Results exported to CSV files")

Step 7: Visualization — Volcano Plot

Prepare the DEG table, then read skills/data-visualization/omics-plotting/SKILL.md and follow its "Volcano" recipe on the exported CSV (→ figures/volcano_plot.png).

pythonCopy code
# Volcano input = log2FoldChange vs padj; saved in Step 6 as deseq2_all_results.csv
volcano_input = results[["log2FoldChange", "padj"]].dropna()
print(f"Volcano input: {len(volcano_input)} genes -> deseq2_all_results.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Volcano" recipe.

Step 8: Visualization — MA Plot

Prepare the DEG table, then read skills/data-visualization/omics-plotting/SKILL.md and follow its "MA plot" recipe on the exported CSV (→ figures/ma_plot.png).

pythonCopy code
# MA input = log2FoldChange vs mean normalized count (baseMean); from deseq2_all_results.csv
ma_input = results[["baseMean", "log2FoldChange", "padj"]].dropna(subset=["padj"])
print(f"MA input: {len(ma_input)} genes -> deseq2_all_results.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "MA plot" recipe.

Key Parameters

Common Recipes

Recipe: Multiple Contrasts from One Model

When comparing multiple treatment groups against a shared control.

pythonCopy code
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()

contrasts = {
    "A_vs_ctrl": ["condition", "treatment_A", "control"],
    "B_vs_ctrl": ["condition", "treatment_B", "control"],
    "C_vs_ctrl": ["condition", "treatment_C", "control"],
}

for name, contrast in contrasts.items():
    ds = DeseqStats(dds, contrast=contrast, alpha=0.05)
    ds.summary()
    n_sig = (ds.results_df.padj < 0.05).sum()
    ds.results_df.to_csv(f"results_{name}.csv")
    print(f"{name}: {n_sig} significant genes")

Recipe: Batch-Corrected Analysis

When samples come from multiple batches or sequencing runs.

pythonCopy code
# Verify batch is not confounded with condition
print(pd.crosstab(metadata["batch"], metadata["condition"]))

dds = DeseqDataSet(
    counts=counts_df, metadata=metadata,
    design="~batch + condition"  # batch first, then variable of interest
)
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()
print(f"Significant genes (batch-corrected): {(ds.results_df.padj < 0.05).sum()}")

Recipe: P-value Distribution QC

Diagnostic check — a healthy analysis shows a flat histogram with a spike near 0.

pythonCopy code
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# P-value histogram
axes[0].hist(ds.results_df.pvalue.dropna(), bins=50, edgecolor="black")
axes[0].set_xlabel("P-value")
axes[0].set_ylabel("Frequency")
axes[0].set_title("P-value Distribution")

# Dispersion plot
axes[1].scatter(
    np.log10(dds.varm["_normed_means"] + 1),
    np.log10(dds.varm["dispersions"]),
    alpha=0.3, s=5
)
axes[1].set_xlabel("Log10(Mean Expression + 1)")
axes[1].set_ylabel("Log10(Dispersion)")
axes[1].set_title("Dispersion vs Mean")

plt.tight_layout()
plt.savefig("qc_diagnostics.png", dpi=300)
print("Saved qc_diagnostics.png")

Recipe: Save and Reload DESeq Dataset

For resuming analysis without re-running the expensive fitting step.

pythonCopy code
import pickle

# Save
with open("dds_fitted.pkl", "wb") as f:
    pickle.dump(dds.to_picklable_anndata(), f)
print("Saved fitted DESeqDataSet")

# Reload
with open("dds_fitted.pkl", "rb") as f:
    adata = pickle.load(f)
# Note: reload requires re-constructing DeseqDataSet from the AnnData

Expected Outputs

• deseq2_all_results.csv — Full results table (baseMean, log2FoldChange, lfcSE, stat, pvalue, padj) for all tested genes
• deseq2_significant.csv — Filtered results (padj < 0.05 and |LFC| > 1)
• deseq2_upregulated.csv — Significant upregulated genes sorted by padj
• deseq2_downregulated.csv — Significant downregulated genes sorted by padj
• figures/volcano_plot.png — Volcano plot with significance and fold change thresholds
• figures/ma_plot.png — MA plot showing fold change vs mean expression
• figures/qc_diagnostics.png — P-value distribution and dispersion plot

Troubleshooting

References

• PyDESeq2 Documentation — Official API reference and tutorials
• Muzellec et al. (2023) — PyDESeq2: a Python package for bulk differential expression analysis with DESeq2. Bioinformatics.
• Love et al. (2014) — Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology.
• DESeq2 Vignette (Bioconductor) — Comprehensive R DESeq2 guide with design formula examples