SuperAPI AI Directory
Partner with us中文AI newsFavorites
Skillsalterlab-pydeseq2

alterlab-pydeseq2

Coding & development

View repository

Updated 10-16 05:00

About this skill

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