alterlab-pydeseq2

PyDESeq2

Overview

PyDESeq2 is a Python implementation of DESeq2 for differential expression analysis with bulk RNA-seq data. It supports complete workflows from data loading through result interpretation, including single-factor and multi-factor designs, Wald tests with multiple-testing correction, optional apeGLM shrinkage, and integration with pandas and AnnData.

When to Use This Skill

Use this skill when:

• Analyzing bulk RNA-seq count data for differential expression
• Comparing gene expression between experimental conditions (e.g., treated vs control)
• Performing multi-factor designs accounting for batch effects or covariates
• Converting R-based DESeq2 workflows to Python
• Integrating differential expression analysis into Python-based pipelines
• Users mention "DESeq2", "differential expression", "RNA-seq analysis", or "PyDESeq2"

Does NOT Trigger

Scenario	Use Instead
Single-cell differential expression or cluster marker genes	alterlab-scanpy (markers) or alterlab-scvi-tools (model-based DE)
Turning FASTQ into the count matrix (salmon/kallisto/STAR quantification, tximport)	alterlab-rnaseq-quant
Generic regression / GLM / mixed models on non-count data	alterlab-statsmodels
Somatic or germline variant calling from sequencing reads	alterlab-nf-core-sarek
Microbiome feature-table differential abundance	alterlab-qiime2-amplicon

Installation and Requirements

uv pip install "pydeseq2>=0.5,<0.6"

System requirements (pydeseq2 0.5.x, current 0.5.4): Python ≥3.11; numpy ≥2.0, pandas ≥2.2, scipy ≥1.12, scikit-learn ≥1.4, anndata ≥0.11, formulaic ≥1.0.2 and formulaic-contrasts ≥0.2 (parse the ~ design formula and build contrast vectors), matplotlib ≥3.9. These are pulled in automatically as dependencies.

API note (0.4+): parallelism is configured through an inference object, not a bare n_cpus= kwarg:

from pydeseq2.default_inference import DefaultInference
inference = DefaultInference(n_cpus=8)
dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition", inference=inference)
ds = DeseqStats(dds, contrast=["condition", "treated", "control"], inference=inference)

Core Workflow

1. Prepare data — load counts as samples × genes (transpose with .T if loaded genes × samples); filter low-count genes (e.g., total reads < 10); drop samples with missing metadata.
2. Specify the design — Wilkinson formula ("~condition", "~batch + condition"); put adjustment variables before the variable of interest.
3. Fit — DeseqDataSet(...).deseq2() runs the full pipeline (size factors → dispersions → LFCs → Cook's outliers).
4. Test — DeseqStats(dds, contrast=[var, test, ref]).summary(); read results_df.
5. (Optional) shrink — ds.lfc_shrink() for visualization/ranking only; p-values stay unshrunken.
6. Interpret/export — filter on padj < 0.05, plot volcano/MA, save CSV/pickle.

Minimal skeleton:

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

dds = DeseqDataSet(counts=counts_df, metadata=metadata, design="~condition")
dds.deseq2()
ds = DeseqStats(dds, contrast=["condition", "treated", "control"])
ds.summary()
significant = ds.results_df[ds.results_df.padj < 0.05]

Command-Line Script

This skill includes a complete standalone script for standard analyses:

python scripts/run_deseq2_analysis.py \
  --counts counts.csv \
  --metadata metadata.csv \
  --design "~batch + condition" \
  --contrast condition treated control \
  --output results/ \
  --min-counts 10 --alpha 0.05 --n-cpus 4 --plots

It handles data loading/validation, gene+sample filtering, the full DESeq2 pipeline, statistical testing with customizable parameters, result export (CSV, pickle), and optional volcano/MA plots. Refer users to scripts/run_deseq2_analysis.py for batch-processing multiple datasets.

Routing Guidance

• Running a standard analysis (load → fit → test → export), or any specific design (two-group, multi-comparison, batch, covariate) → references/pipeline_steps.md.
• Interpreting results, ranking genes, plotting volcano/MA, or quality metrics → references/interpretation_and_plots.md.
• Hitting an error (index mismatch, all-zero counts, "not full rank", no significant genes) → Troubleshooting in references/interpretation_and_plots.md.
• Need exact class/method parameters or object attributes → references/api_reference.md.
• Complex experimental designs or in-depth workflow → references/workflow_guide.md.

Key Reminders

1. Data orientation matters: counts usually load genes × samples but need samples × genes — transpose with .T if needed.
2. Sample filtering: remove samples with missing metadata before analysis.
3. Gene filtering: drop low-count genes (e.g., < 10 total reads) to improve power.
4. Design formula order: adjustment variables before the variable of interest ("~batch + condition").
5. LFC shrinkage timing: shrink after testing, for visualization/ranking only — p-values stay unshrunken.
6. Significance: use padj < 0.05 (Benjamini-Hochberg FDR), not raw p-values.
7. Contrast format: [variable, test_level, reference_level]. contrast also accepts a raw numpy contrast vector over the design matrix columns for comparisons a three-element list cannot express (e.g. interaction terms, averaging several levels).
8. Save intermediates: pickle the DeseqDataSet to avoid re-running the expensive fit.
9. Test against a fold-change threshold, not just zero: DeseqStats(..., lfc_null=1.0, alt_hypothesis="greaterAbs") asks "is |LFC| > 1?" inside the model. That is the statistically correct way to demand an effect size — filtering a lfc_null=0 result on `abs(log2FoldChange)

   1` afterwards does not control the FDR for that claim.

10. Zero-heavy or sparse counts: DeseqDataSet(..., size_factors_fit_type="poscounts") uses the positive-counts estimator instead of the median-of-ratios default, which fails when no gene is detected in every sample. control_genes= restricts size-factor estimation to spike-ins or housekeeping genes.

Reference Index

• references/pipeline_steps.md — Quick-start, the six pipeline steps with full code (data prep, design, fitting, testing, shrinkage, export), and four common experimental designs.
• references/interpretation_and_plots.md — Filtering/ranking significant genes, quality metrics, volcano and MA plots, and a troubleshooting guide.
• references/api_reference.md — Complete PyDESeq2 class/method/parameter and data-structure documentation.
• references/workflow_guide.md — In-depth complete workflows, data-loading patterns, multi-factor designs, and best practices.

Additional Resources

• Official Documentation: https://pydeseq2.readthedocs.io
• GitHub Repository: https://github.com/owkin/PyDESeq2
• Publication: Muzellec et al. (2023) Bioinformatics, DOI: 10.1093/bioinformatics/btad547
• Original DESeq2 (R): Love et al. (2014) Genome Biology, DOI: 10.1186/s13059-014-0550-8