Differential gene expression analysis for bulk RNA-seq with PyDESeq2, including formulaic designs, Wald tests, FDR correction, LFC shrinkage, and result visualization.
SKILL.md
PyDESeq2
Overview
PyDESeq2 is a Python implementation of DESeq2 for differential expression analysis with bulk RNA-seq data. Design and execute complete workflows from data loading through result interpretation, including formulaic 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
This skill should be used 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"
Quick Start Workflow
For users who want to perform a standard differential expression analysis:
Data preparation — raw integer counts with genes as columns and samples as rows,
and matching metadata. Never feed normalized or transformed values to DESeq2.
Design specification — the design factors and the reference level for each.
DESeq2 fitting — size factors, dispersions, and the GLM fit.
Statistical testing — Wald tests for a named contrast.
Optional LFC shrinkage — for ranking and visualization.
Result export — the results table with adjusted p-values.
Issue: "Index mismatch between counts and metadata"
Solution: Ensure sample names match exactly
print("Counts samples:", counts_df.index.tolist())
print("Metadata samples:", metadata.index.tolist())
# Take intersection if needed
common = counts_df.index.intersection(metadata.index)
counts_df = counts_df.loc[common]
metadata = metadata.loc[common]
Issue: "All genes have zero counts"
Solution: Check if data needs transposition
print(f"Counts shape: {counts_df.shape}")
# If genes > samples, transpose is needed
if counts_df.shape[1] < counts_df.shape[0]:
counts_df = counts_df.T
Design Matrix Issues
Issue: "Design matrix is not full rank"
Cause: Confounded variables (e.g., all treated samples in one batch)
Solution: Remove confounded variable or add interaction term
# Check confounding
print(pd.crosstab(metadata.condition, metadata.batch))
# Either simplify design or add interaction
design = "~condition" # Remove batch
# OR
design = "~condition + batch + condition:batch" # Model interaction
No Significant Genes
Diagnostics:
# Check dispersion distribution
plt.hist(dds.var["dispersions"], bins=50)
plt.show()
# Check size factors
print(dds.obs["size_factors"])
# Look at top genes by raw p-value
print(ds.results_df.nsmallest(20, "pvalue"))
Possible causes:
Small effect sizes
High biological variability
Insufficient sample size
Technical issues (batch effects, outliers)
Reference Documentation
For comprehensive details beyond this workflow-oriented guide:
API Reference (references/api_reference.md): Complete documentation of PyDESeq2 classes, methods, and data structures. Use when needing detailed parameter information or understanding object attributes.
Workflow Guide (references/workflow_guide.md): In-depth guide covering complete analysis workflows, data loading patterns, multi-factor designs, troubleshooting, and best practices. Use when handling complex experimental designs or encountering issues.
Load these references into context when users need:
Detailed API documentation: Read references/api_reference.md
Troubleshooting guidance: Read references/workflow_guide.md (see Troubleshooting section)
Key Reminders
Data orientation matters: Count matrices typically load as genes × samples but need to be samples × genes. Always transpose with .T if needed.
Sample filtering: Remove samples with missing metadata before analysis to avoid errors.
Gene filtering: Filter low-count genes (e.g., < 10 total reads) to improve power and reduce computational time.
Design formula order: Put adjustment variables before the variable of interest (e.g., "~batch + condition" not "~condition + batch").
LFC shrinkage timing: Apply shrinkage after statistical testing and only for visualization/ranking purposes. P-values remain based on unshrunken estimates.
Result interpretation: Use padj < 0.05 for significance, not raw p-values. The Benjamini-Hochberg procedure controls false discovery rate.
Contrast specification: The format is [variable, test_level, reference_level] where test_level is compared against reference_level.
Save intermediate objects: Prefer dds.to_picklable_anndata().write_h5ad("dds_result.h5ad") for portable outputs. Only load pickle files that you created yourself and trust.
