Production-ready RNA-seq differential expression analysis using PyDESeq2. Performs DESeq2 normalization, dispersion estimation, Wald testing, LFC shrinkage, and result filtering. Handles multi-factor designs, multiple contrasts, batch effects, and integrates with gene enrichment (gseapy) and ToolUniverse annotation tools (UniProt, Ensembl, OpenTargets). Supports CSV/TSV/H5AD input formats and any organism. Use when analyzing RNA-seq count matrices, identifying DEGs, performing differential expression with statistical rigor, or answering questions about gene expression changes.
Differential expression analysis of RNA-seq count data using PyDESeq2, with enrichment analysis (gseapy) and gene annotation via ToolUniverse.
BixBench Coverage: Validated on 53 BixBench questions across 15 computational biology projects.
DESeq2 assumes that most genes are NOT differentially expressed — this is its normalization assumption. If this assumption is violated (e.g., global transcriptional shutdown, where the majority of genes genuinely decrease), size factor normalization will inflate expression in the treatment group and produce artifactually upregulated genes. Always check the MA plot: the fold-change cloud should be centered on zero across all expression levels. A systematic upward or downward shift indicates a normalization problem, not biology.
MyGene_query_genes, UniProt); do not recall gene function or pathway from memory.metadata.columns and from the actual data; do not assume metadata structure.metadata[factor].unique()import pandas as pd, numpy as np
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats
import gseapy as gp # enrichment (optional)
from tooluniverse import ToolUniverse # annotation (optional)
Extract: data files, thresholds (padj/log2FC/baseMean), design factors, contrast, direction, enrichment type, specific genes. See question_parsing.md.
Load counts + metadata, ensure samples-as-rows/genes-as-columns, verify integer counts, align sample names, remove zero-count genes. See data_loading.md.
List ALL metadata columns and levels. Categorize as biological interest vs batch/block. Build design formula with covariates first, factor of interest last. See design_formula_guide.md.
Set reference level via pd.Categorical, create DeseqDataSet, call dds.deseq2(), extract DeseqStats with contrast, run Wald test, optionally apply LFC shrinkage. See pydeseq2_workflow.md.
Tool boundaries:
Apply padj, log2FC, baseMean thresholds. Split by direction if needed. See result_filtering.md.
Key columns: genewise_dispersions, fitted_dispersions, MAP_dispersions, dispersions. See dispersion_analysis.md.
Use gseapy enrich() with appropriate gene set library. See enrichment_analysis.md.
Use ToolUniverse for ID conversion and gene context only. See output_formatting.md.
| Pattern | Type | Key Operation |
|---|---|---|
| 1 | DEG count | len(results[(padj<0.05) & (abs(lfc)>0.5)]) |
| 2 | Gene value | results.loc['GENE', 'log2FoldChange'] |
| 3 | Direction | Filter log2FoldChange > 0 or < 0 |
| 4 | Set ops | degs_A - degs_B for unique DEGs |
| 5 | Dispersion | (dds.var['genewise_dispersions'] < thr).sum() |
See bixbench_examples.md for all 10 patterns with examples.
| Error | Fix |
|---|---|
| No matching samples | Transpose counts; strip whitespace |
| Dispersion trend no converge | fit_type='mean' |
| Contrast not found | Check metadata['factor'].unique() |
| Non-integer counts | Round to int OR use t-test |
| NaN in padj | Independent filtering removed genes |
See troubleshooting.md for full debugging guide.
| Metric | Threshold | Interpretation |
|---|---|---|
| padj | < 0.05 | Statistically significant after multiple testing correction |
| log2FoldChange | > 1 or < -1 | Biologically meaningful fold change (2x up or down) |
| baseMean | > 10 | Gene is expressed at detectable levels |
| lfcSE | < 1.0 | Fold change estimate is precise |
| Grade | Criteria | Action |
|---|---|---|
| Strong DEG | padj < 0.01, | LFC |
| Moderate DEG | padj < 0.05, | LFC |
| Weak DEG | padj < 0.1 or | LFC |
| Not significant | padj >= 0.1 | Do not report as differentially expressed |