domains / science / rna-seq-differential-expression / pydeseq2

pydeseq2

Skill card

combined

The K-Dense-AI copy is a rigorous, version-pinned (PyDESeq2 0.5.4) skill that includes a CLI driver, strict data validation rules (rejecting silent intersections), and specific output artifacts. The other three copies (Vibe-Skills, medical-research-skills, OpenClaw-Medical-Skills) are nearly identical to each other, providing a standard Python API workflow with inline code examples and a manual volcano plot script. The K-Dense copy lacks the specific Python API code patterns found in the others, while the others lack the strict validation logic, CLI usage, and version constraints of the K-Dense copy. Combining them provides both the robust procedural guardrails and the practical code implementation.

Repositories

Other copies

Combined skill

---
name: pydeseq2
description: Performs bulk RNA-seq differential expression analysis with PyDESeq2, including count validation, formula designs, explicit contrasts, Wald tests, FDR correction, coefficient-matched LFC shrinkage, and result visualization.
---

# PyDESeq2

Use this skill for bulk RNA-seq negative-binomial GLMs with explicit biological replicates. Do not use for single-cell per-cell testing.

## Environment

- **Python:** >= 3.11 (Recommended 3.13 for current stack)
- **PyDESeq2:** 0.5.4
- **AnnData:** 0.13.4 (if using h5ad)

Install in an isolated environment:
```bash
uv venv --python 3.13 .venv-pydeseq2
uv pip install --python .venv-pydeseq2/bin/python 'pydeseq2==0.5.4' 'anndata==0.13.4'
```

## Analysis Contract

1. **Experimental Unit:** Ensure samples are independent biological replicates.
2. **Data Orientation:** Counts must be **samples × genes** with non-negative integer counts. If input is genes × samples, transpose it.
3. **Validation:** 
   - **Never silently intersect** sample IDs. Resolve unmatched IDs and missing design annotations explicitly.
   - Reject duplicate headers, unequal sample sets, and invalid count values.
4. **Prefiltering:** Apply a prespecified gene prefilter (e.g., total counts >= 10) after sample exclusions.
5. **Design:** Build a full-rank formula. Put adjustment variables (batch) before the variable of interest (condition).

## Workflow

### Option A: CLI Driver (Recommended for reproducibility)

Use the bundled driver for strict validation and standardized outputs.

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

*Note: The driver expects counts CSV as genes × samples by default. Use `--no-transpose` if already samples × genes.*

**Outputs:**
- `deseq2_results.csv`: Shrunken estimates (if requested) with original Wald stats.
- `deseq2_results_unshrunken.csv`: Original Wald table.
- `significant_genes.csv`: Genes with `padj < alpha`.
- `deseq_dataset.h5ad`: Fitted AnnData snapshot.
- `volcano_plot.png` / `ma_plot.png`: Visualizations.

### Option B: Python API (For custom pipelines)

```python
import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# 1. Load and Orient Data
counts_df = pd.read_csv("counts.csv", index_col=0).T  # Ensure samples x genes
metadata = pd.read_csv("metadata.csv", index_col=0)

# 2. Validate Alignment (Explicit check, do not silently intersect)
if not counts_df.index.equals(metadata.index):
    raise ValueError("Sample IDs in counts and metadata do not match exactly.")

# 3. Filter Low-Count Genes
genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10]
counts_df = counts_df[genes_to_keep]

# 4. Fit Model
dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design="~batch + condition",
    refit_cooks=True
)
dds.deseq2()

# 5. Statistical Testing
ds = DeseqStats(
    dds,
    contrast=["condition", "treated", "control"],
    alpha=0.05,
    cooks_filter=True,
    independent_filter=True
)
ds.summary()

# 6. Results and Shrinkage
results = ds.results_df.copy()

# Optional: Shrink LFC for visualization/ranking (p-values remain unchanged)
ds.lfc_shrink()
results_shrunk = ds.results_df.copy()

# 7. Export
results.to_csv("deseq2_results_unshrunken.csv")
results_shrunk.to_csv("deseq2_results_shrunk.csv")
significant = results[results.padj < 0.05]
significant.to_csv("significant_genes.csv")
```

## Statistical Interpretation

- **Contrast:** `["condition", "treated", "control"]` tests **treated / control**.
- **Significance:** Use `padj < alpha` (e.g., 0.05). Missing `padj` is never significant.
- **Shrinkage:** Use shrunk `log2FoldChange` for volcano plots and ranking. Report unshrunken `pvalue` and `padj` for statistical significance.
- **Volcano Plot:** Plot `log2FoldChange` vs `-log10(padj)`. Cap zero adjusted p-values for plotting only.