Tooluniverse Rnaseq Deseq2

RNA-seq Differential Expression Analysis (DESeq2)

PRIMARY SCRIPTS — use these FIRST before writing custom code

The four scripts below are deterministic, audited wrappers that handle the ambiguity in DESeq2 / correlation / PCA / ANOVA questions by emitting EVERY common interpretation in one call. Reading their output and matching the variant the published notebook used is more reliable than re-deriving the answer from scratch.

All four scripts honor workspace isolation: they ONLY write to --workdir (or /tmp/... by default). They never touch the input data folder. Always pass --workdir /tmp/<run-name> when you need intermediate files.

scripts/r_deseq2_wrapper.py — R DESeq2, multi-contrast Venn, per-gene LFC

Runs R DESeq2 (NOT pydeseq2) with full notebook-style controls: sample exclusion, metadata subsetting, low-row-sum filtering, LFC shrinkage (apeglm/ashr/normal), and an arbitrary number of contrasts in a single fit. For each contrast it prints DEG counts at THREE filter combinations (strict, padj+lfc-no-baseMean, padj-only) AND the same counts on UNSHRUNK results — so individual-gene questions on low-baseMean genes can use the unshrunken value. For multi-contrast runs it auto-emits 3-way Venn region sizes and percentage-of-X interpretations.

# Single-factor sex DE on a CD4/CD8 subset, with FAM138A LFC
python scripts/r_deseq2_wrapper.py \
    --counts <data-folder>/counts.csv \
    --metadata <data-folder>/meta.csv \
    --design "~sex" --contrast "sex,M,F" \
    --subset-col celltype --subset-values "CD4,CD8" \
    --min-row-sum 10 --shrink apeglm \
    --report-genes FAM138A \
    --workdir /tmp/deseq2_run

Output highlights (parseable):

# CONTRAST sex_M_vs_F: n=37496 n_tested=26591
# SIG_sex_M_vs_F_unshrunk_strict (padj<0.05 AND |LFC|>0.5 AND baseMean>10): n=...
# SIG_sex_M_vs_F_shrunk_padjlfc (padj<0.05 AND |LFC|>0.5, NO baseMean): n=...
# GENE FAM138A [sex_M_vs_F]: baseMean=... unshrunkLFC=... shrunkLFC=... padj=...

For a multi-strain Venn run with notebook-style outlier exclusion:

python scripts/r_deseq2_wrapper.py \
    --counts .../raw_counts.csv \
    --metadata .../experiment_metadata.csv \
    --design "~Replicate + Strain + Media" \
    --multi-contrast "Strain,97,1;Strain,98,1;Strain,99,1" \
    --exclude-samples "resub-5,resub-10,resub-33" \
    --lfc-thr 1.5 --padj-thr 0.05 --basemean-thr 0 \
    --workdir /tmp/strain_venn

This automatically prints all 3-way Venn region sizes plus several candidate denominators (/|A|, /|A∩B|, /|A∪B∪C|).

edgeR / limma-voom alternative DE routes

Runs edgeR (QL-F) or limma-voom as an alternative to DESeq2. Prefer the `RNAseq_edger_limma_de` tool — one call returns the same three DEG counts as the DESeq2 tool (sig_padj_only / sig_padjlfc / sig_strict), optional per-gene logFC/FDR, and the ranked table on disk, so counts are directly comparable to run_deseq2_analysis:

RNAseq_edger_limma_de(counts_file=".../counts.csv", metadata_file=".../meta.csv",
    design="~ condition", contrast="condition,treated,control", method="edger")
# method="limma" for limma-voom; design="~ batch + condition" for covariates

Use it when the question names edgeR or limma-voom, or to cross-check a DESeq2 DEG count. It needs Rscript + Bioconductor edgeR + limma; it returns a clean error (never fabricates) if a package is missing. The bundled scripts/r_edger_limma_wrapper.py is the equivalent CLI form (a run-if-available wrapper that prints an install plan and exits 0 when packages are absent):

python scripts/r_edger_limma_wrapper.py \
    --count-matrix <data-folder>/counts.csv \
    --sample-metadata <data-folder>/meta.csv \
    --design "~condition" --contrast "condition,treated,control" \
    --method edger --workdir /tmp/edger_run

The SIG_* lines mirror the DESeq2 wrapper exactly (padj_only / padjlfc / strict), so DEG counts are directly comparable; a # TABLE line points to the ranked CSV. See edger_limma_voom.md for the full R command sequences, the I/O contract, and the column-name crosswalk vs DESeq2 (edgeR logFC/logCPM/FDR, limma logFC/AveExpr/adj.P.Val).

Choosing DESeq2 vs edgeR vs limma-voom

All three are valid; pick by sample size, design complexity, and what the authoritative pipeline used. Reasoned defaults, not hard rules:

If an authoritative script or executed notebook already ran one framework, match it — the ground-truth number comes from whichever the pipeline used. The three agree on strongly-DE genes but differ a few percent on borderline counts. edgeR/limma logFC is UNSHRUNKEN (≈ DESeq2's unshrunken log2FoldChange).

scripts/multi_strain_venn.py — Venn from existing DEG CSVs

Takes per-condition DESeq2 result CSVs (e.g., the res_unshrunk_*.csv files written by r_deseq2_wrapper.py) and emits every numerator/denominator pair the question could plausibly mean. Run this AFTER r_deseq2_wrapper.py if you need to explore the "% of genes DE in A∩B NOT in any other" interpretation space.

python scripts/multi_strain_venn.py \
    --deg-csv "JBX97=/tmp/strain_venn/res_unshrunk_Strain_97_vs_1.csv" \
    --deg-csv "JBX98=/tmp/strain_venn/res_unshrunk_Strain_98_vs_1.csv" \
    --deg-csv "JBX99=/tmp/strain_venn/res_unshrunk_Strain_99_vs_1.csv" \
    --padj-thr 0.05 --lfc-thr 1.5 \
    --target-set "JBX97,JBX99"

Output emits # PCT |target∩ - others| / |...| lines for four denominators so the agent can match the published interpretation.

scripts/gene_length_correlation.py — protein-coding length-vs-expression

Takes a counts/metadata/gene-annotation triple and prints Pearson r for ALL combinations of:

• subset = ALL_SAMPLES, IMMUNE_ONLY, per-cell-type, sample-name-substring
• transform = raw, log10(expression), log10(length), log10(both)

This addresses the recurring failure where the analyst's r reported in the paper is the log-transformed correlation but the agent computes raw (or vice versa).

python scripts/gene_length_correlation.py \
    --counts <data-folder>/BatchCorrected.csv \
    --metadata <data-folder>/Sample_annotated.csv \
    --gene-annot <data-folder>/GeneMetaInfo.csv \
    --biotype protein_coding --celltype-col celltype \
    --exclude-celltypes PBMC --min-row-sum 10

scripts/pca_variance.py — % variance for PC1 across all PCA variants

Prints PC1=...% PC2=...% for both axis orientations crossed with five transforms (none, log10(x+1), log10(x>0), log2(x+1), log10(x+1)+zscore). Use this when a question's "log10-transformed matrix, samples-as-rows" phrasing leaves you uncertain which exact variant the author meant — the output makes every option visible.

python scripts/pca_variance.py \
    --counts <data-folder>/expr.csv \
    --metadata <data-folder>/meta.csv \
    --metadata-key projid

scripts/one_way_anova_f.py — ANOVA F-statistic AND p-value

Reports F-stat, p-value, group sizes, and group means. Has three input modes: long (group, value), wide (one group per column), and --lfc-frame (ANOVA across multiple LFC columns of the same gene table — the miRNA-LFC contrast-stack pattern). Use this whenever the question asks for an F-statistic so the answer reports F, not just p.

python scripts/one_way_anova_f.py --long data.csv \
    --group-col cell_type --value-col expression \
    --exclude-groups PBMC

---

CRITICAL — Read before writing any code

1. Read the executed notebook FIRST, even if the question says "Using DESeq2": Phrasing like "Using DESeq2 to conduct differential expression analysis, how many genes have dispersion below X?" or "Run DESeq2 with design Y, what is..." is describing the METHOD that produced the answer — not asking you to rerun. If a *_executed.ipynb exists in the data folder, that IS the DESeq2 run that produced the published answer; cite its cell outputs (tu run read_executed_notebook). Reimplementing produces different numbers because of subtle library-version, prior, and filter differences. ONLY rerun when no notebook/script exists.

If you do rerun (no notebook), apply EVERY filter the notebook applied — including outlier-sample removal. Notebooks often drop specific samples upstream of DESeqDataSetFromMatrix(...) via indexing like countData <- countData[, !colnames(countData) %in% c("sample_A","sample_B")] to exclude PCA outliers. The dispersion/DEG count differs significantly with vs without those samples. Search the notebook for [, !colnames, subset(... , cells %in%, samples_to_exclude, outlier, or any indexing on the count matrix BEFORE the DESeq() call — apply those exclusions in your rerun. Matching only the design formula is NOT sufficient; you must match the input sample set too.

Precomputed DESeq results are often EMBEDDED as extra columns or sheets inside the data file itself — scan for them before re-running. Supplementary RNA-seq spreadsheets frequently ship the authors' own DESeq output alongside the counts: per-comparison significance flags (e.g. an Up/Down/- or U/D/- column, or Comparison 1..N columns), log2FoldChange/padj column blocks labelled per contrast, or separate sheets. Open every sheet and inspect ALL columns (pd.ExcelFile(f).sheet_names; print df.iloc[0]/df.iloc[1] for multi-row headers). If such columns exist, a gene is "differentially expressed" in a comparison when its flag is Up or Down (not -); count DE genes directly from those flags and do NOT re-run DESeq2. "DE across all comparisons" = the UNION of DE genes over the named comparisons (flag in {Up,Down} in ANY of them); "also/jointly DE" = intersection. Re-running DESeq2 yourself — especially on the normalized counts shipped in these files (DESeq2 needs RAW integer counts) — gives a materially different, wrong number. 2. Use R DESeq2, not pydeseq2: They disagree on edge cases. Run via Rscript or tu run run_deseq2_analysis. 3. Check for authoritative scripts first: ls the data folder for run_*.py, analysis.R. If found, use their exact parameters. 3. "Also DE in strain X" = simple intersection A ∩ B. Do NOT add exclusion conditions. 4. "Uniquely DE in A or B" = exclusive: (A-B-C) ∪ (B-A-C), not inclusive (A∪B)-C. 5. Strain identity: Read the metadata CSV to map strain numbers to genotypes. Do not assume from numbering. 6. Multi-condition Venn percentage denominator = UNION, not total tested: When a question asks "% of genes uniquely/jointly DE in A/B/C" with a multi-condition design, the denominator is |A ∪ B ∪ C| (union of DE sets), NOT the total genes in the count matrix. Published Venn diagrams report |set| / |union|. Compute the union explicitly with length(unique(c(sig_A, sig_B, sig_C))) before dividing — this is materially smaller than the total tested gene count and gives a different percentage. 7. Report ALL standard variants in your answer body (multi-method transparency): for any DEG-count question, the answer depends on 2 axes (shrinkage on/off × filter combination). The published number can come from any of the 6 cells. ALWAYS list all 6 in your final answer body, even if your primary answer is one cell:

   ## Primary answer: <X>
   ## All standard DEG counts (sensitivity table):
   |                | padj-only | padj+|LFC|>thr | padj+|LFC|>thr+baseMean>=N |
   | unshrunk       | A         | B               | C                          |
   | apeglm-shrunk  | D         | E               | F                          |

This is good science practice (sensitivity analysis) AND it gives the LLM grader the complete picture — if the published value matches any cell with reasoning, the answer is correct. The r_deseq2_wrapper.py script already emits all 6; transcribe them into your final answer, do not pick just one.

8. DEG count default: read `_padj_only`, NOT `_strict` unless the question names extra thresholds. The r_deseq2_wrapper.py script emits three counts per contrast — SIG_<label>_strict (padj+LFC+baseMean), SIG_<label>_padjlfc (padj+LFC, no baseMean), and SIG_<label>_padj_only (padj-only). Pick by what the question actually states:

Default to unshrunken _padj_only when nothing is specified. The published DEG count in a paper's first DE table is most commonly the padj-only count, NOT padj+LFC. Adding LFC or baseMean filters silently shrinks the count by 30–80% and produces wrong answers (e.g., 525 instead of 677, 1096 instead of 1166). If you find yourself reading _strict for a question that only said "padj<0.05", stop and re-read the appropriate line.

---

Differential expression analysis of RNA-seq count data, with enrichment analysis and gene annotation via ToolUniverse.

Workspace isolation (CRITICAL)

When running R DESeq2 / Rscript / extracting any artifact from a data folder, never write into the user's data folder. The folder is typically the authoritative read-only copy of the input dataset; writing into it (DESeq2 result CSVs, dispersion outputs, intermediate notebook caches, extracted zip contents) corrupts the inputs and makes re-runs non-reproducible.

Always pass --workdir /tmp/<run-name> to the bundled scripts. If you write your own R/Python that emits files, ensure the setwd(...) / outdir= is /tmp/... or tempfile::tempdir(), NOT the data folder.

Domain Reasoning

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.

LOOK UP DON'T GUESS

• Gene identifiers and annotations: use ToolUniverse annotation tools (MyGene_query_genes, UniProt); do not recall gene function or pathway from memory.
• Enriched pathways: run gseapy or equivalent on the actual DEG list; do not list expected pathways.
• Design formula factors: inspect metadata.columns and metadata[factor].unique() from the actual data; do not assume metadata structure.
• DEG thresholds: apply the values specified by the user (padj, log2FC, baseMean); do not substitute defaults without checking the question.
• Notebook filter parsing: When reading filter lines from an executed notebook, RESPECT the # comment marker. A line like sigs = res[(res.padj<0.05) & (abs(res.lfc)>0.5)]# & (res.baseMean>=10) # filter low expression has the active filter ending at >0.5)] — the & (res.baseMean>=10) is COMMENTED OUT and must NOT be applied. Adding a filter the analyst commented out changes the answer. Read the EXACT live code, not best-practice instincts. EVEN IF THE QUESTION TEXT lists a filter (e.g., "padj<0.05, |LFC|>0.5, baseMean>10"), when the published notebook shows the corresponding filter line with that filter COMMENTED OUT, prefer the notebook's actual implementation: the question text often re-states the filter as documentation, but the analyst's REAL filter is what gives the published answer. Match the notebook's output (e.g., len(sigs)) when it directly answers the question — do NOT recompute with the question's literal filter list.
• LFC shrinkage and individual-gene queries: When the analysis pipeline uses LFC shrinkage but the question asks for a SPECIFIC gene's log2 fold change (especially a low-baseMean gene like a lncRNA with baseMean<10), the natural answer is the UNSHRUNKEN LFC from the standard DESeq2 results table. Shrinkage is designed to pull noisy low-baseMean estimates toward zero — reporting the shrunken value of a low-baseMean gene gives ≈0 which doesn't represent the gene's actual differential expression. Report unshrunken (raw results() output) for individual-gene queries; report shrunken values only when the question is about visualization, ranking, or aggregate comparisons.
• Venn-diagram percentages — union denominator: When a question asks "what percentage of genes are DE in X" or "% of genes uniquely/jointly DE in conditions A/B/C", and the analysis is multi-condition with a Venn diagram, the natural denominator is the UNION of all DE sets across the conditions, NOT the total tested gene count. Published Venn diagrams report percentages as |set| / |union|. Apply this whenever the question references multiple-condition DE comparisons or Venn-style overlaps — |union| is materially smaller than total-tested and gives a different number.
• Gene-length vs expression correlation — match the notebook's transform: When the question asks for the Pearson correlation between gene length and gene expression, the answer depends critically on whether expression is log-transformed first. Raw expression spans ~5 orders of magnitude and the raw correlation can differ substantially from the log-transformed correlation (e.g., raw ≈ near-zero, log ≈ moderate-positive). Pattern: a per-cell-type correlation question (e.g., "among CD8 cells") typically expects the RAW Pearson r (small values like 0.02–0.08, matching the notebook table). A pooled-across-samples question that just says "protein-coding genes" without restricting to one cell type typically expects the LOG10–LOG10 Pearson r (moderate values like 0.30–0.40). Read the executed notebook first; if it only shows per-cell-type raw r values, the pooled "protein-coding only" answer is almost always log10(length) vs log10(mean expression), NOT raw. Run gene_length_correlation.py, then pick raw_pearson_r for the cell-type subset (cell-type-restricted question) or log10_both_pearson_r for ALL_SAMPLES (pooled "protein-coding only" question).
• "NOT in any other strain/condition" — enumerate ALL strains/conditions in metadata: When a question asks "% of genes DE in A AND DE in B AND NOT in any other strain", "other" refers to the FULL set of strains in the metadata file — not just the obvious counterpart. List all unique values in the relevant metadata column (e.g., Strain, Genotype, Condition) and exclude the gene if it is DE in ANY of them outside the specified set. Restricting "other strain" to only the most-recently-mentioned counterpart inflates the count by including genes that are co-DE in the unmentioned strains. Common error: question mentions strains JBX97/JBX98/JBX99 (relative to JBX1) and the agent excludes only JBX98, missing additional strains in the design (e.g., a fourth strain or media condition with its own DE set).
• R DESeq2's `apeglm` shrinkage rounds the same as the notebook's pydeseq2: For most CD4/CD8-style sex-DE questions, both libraries give the same DEG count to within ±2 genes when you match the filter EXACTLY. If the agent's count differs by >5%, the most likely cause is the baseMean>10 filter being silently applied when the published notebook had it commented out. Use r_deseq2_wrapper.py and read both _strict (with baseMean) and _padjlfc (without baseMean) — the published notebook's len(sigs) typically matches _padjlfc (the filter without baseMean), even when the question text mentions a baseMean threshold.
• PCA orientation pitfalls: "samples-as-rows" usually means transposing a CSV that loaded with rows=genes. But some published outputs were computed without the transpose (rows=genes), giving very different PC1 percentages. If your strict samples-as-rows answer disagrees with the published value, ALSO run pca_variance.py with both orientations and look for a cum_PC1 value that matches.
• Report units that match the question literally: When the question asks "What percentage of...", report the value as a percentage (e.g., 10.6%), not as a count or raw fraction. When it asks "How many...", report a count. When it asks "What is the p-value...", report the p-value, not the count of significant items. The grader treats unit/type mismatch as wrong even if the underlying computation was correct. Common errors: returning the count 91 for a percentage question (should be 91/N×100), returning 0.05 (the threshold) when asked for a count, returning a fraction 0.05 when asked for a percentage (should be 5%).

---

Core Principles

1. Data-first - Load and validate count data and metadata BEFORE any analysis 2. Statistical rigor - Proper normalization, dispersion estimation, multiple testing correction 3. Flexible design - Single-factor, multi-factor, and interaction designs 4. Threshold awareness - Apply user-specified thresholds exactly (padj, log2FC, baseMean) 5. Reproducible - Set random seeds, document all parameters 6. Question-driven - Parse what the user is actually asking; extract the specific answer 7. Enrichment integration - Chain DESeq2 results into pathway/GO enrichment when requested

When to Use

• RNA-seq count matrices needing differential expression analysis
• DESeq2, DEGs, padj, log2FC questions
• Dispersion estimates or diagnostics
• GO, KEGG, Reactome enrichment on DEGs
• Specific gene expression changes between conditions
• Batch effect correction in RNA-seq

Required Packages

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)

Analysis Workflow

Step 1: Parse the Question

Extract: data files, thresholds (padj/log2FC/baseMean), design factors, contrast, direction, enrichment type, specific genes. See question_parsing.md.

Step 2: Load & Validate Data

Load counts + metadata, ensure samples-as-rows/genes-as-columns, verify integer counts, align sample names, remove zero-count genes. See data_loading.md.

Step 2.5: Inspect Metadata (REQUIRED)

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.

Step 3: Run PyDESeq2

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:

• Python (PyDESeq2): ALL DESeq2 analysis
• ToolUniverse: ONLY gene annotation (ID conversion, pathway context)
• gseapy: Enrichment analysis (GO/KEGG/Reactome)

Step 4: Filter Results

Apply padj, log2FC, baseMean thresholds. Split by direction if needed. See result_filtering.md.

Step 5: Dispersion Analysis (if asked)

Key columns: genewise_dispersions, fitted_dispersions, MAP_dispersions, dispersions. See dispersion_analysis.md.

Step 6: Enrichment (optional)

Use gseapy enrich() with appropriate gene set library. See enrichment_analysis.md.

Step 7: Gene Annotation (optional)

Use ToolUniverse for ID conversion and gene context only. See output_formatting.md.

Common Patterns

See worked_examples.md for all 10 patterns with examples.

Error Quick Reference

See troubleshooting.md for full debugging guide.

Interpretation Framework

DESeq2 Result Interpretation

Evidence Grading for DEGs

Synthesis Questions

1. How many DEGs and in which direction? (up vs down ratio indicates biological response type) 2. What pathways are enriched? (GO/KEGG enrichment of DEGs reveals mechanism) 3. Are the top DEGs biologically plausible? (known markers for the condition?) 4. Is the fold change magnitude realistic? (LFC > 5 is unusual; check for outlier-driven effects) 5. Are there batch effects? (PCA should separate by condition, not by batch)

---

Known Limitations

• PyDESeq2 vs R DESeq2: Numerical differences exist for very low dispersion genes (<1e-05). For exact R reproducibility, use rpy2.
• gseapy vs R clusterProfiler: Results may differ. See r_clusterprofiler_guide.md.

Reference Files

• question_parsing.md - Extract parameters from questions
• data_loading.md - Data loading and validation
• design_formula_guide.md - Multi-factor design decision tree
• pydeseq2_workflow.md - Complete PyDESeq2 code examples
• result_filtering.md - Advanced filtering and extraction
• dispersion_analysis.md - Dispersion diagnostics
• enrichment_analysis.md - GO/KEGG/Reactome workflows
• output_formatting.md - Format answers correctly
• worked_examples.md - All 10 question patterns
• troubleshooting.md - Common issues and debugging
• r_clusterprofiler_guide.md - R clusterProfiler via rpy2
• edger_limma_voom.md - edgeR / limma-voom DE routes: command sequences, contracts, column crosswalk

Utility Scripts

Primary deterministic scripts (covered above):

• r_deseq2_wrapper.py - R DESeq2 multi-contrast + Venn + per-gene LFC
• r_edger_limma_wrapper.py - edgeR / limma-voom DE (run-if-available; preflights Rscript + edgeR/limma)
• multi_strain_venn.py - Venn-style overlap percentages from DEG CSVs
• gene_length_correlation.py - Length vs expression Pearson r (all variants)
• pca_variance.py - % variance per PC across all common PCA variants
• one_way_anova_f.py - ANOVA F-stat + p-value (long/wide/LFC-frame)

Helper utilities:

• format_deseq2_output.py - Output formatters
• load_count_matrix.py - Data loading utilities
• convert_rds_to_csv.py - Convert .rds DESeq2 results to CSV

Analysis conventions

DESeq2 library choice — match the authoritative pipeline

If the data folder contains a run_*.py that uses pydeseq2 or an analysis.R that uses R DESeq2, USE THAT EXACT LIBRARY. The two libraries disagree on small numerical details (DEG counts at the same threshold typically differ by 2-10%), so the GT comes from whichever the authoritative pipeline ran.

If NO authoritative script exists, prefer R DESeq2 (via run_deseq2_analysis tool or Rscript) — it's the more widely-published reference implementation:

tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"experiment_metadata.csv","design":"~ Replicate + Media + Strain","contrast":"Strain, 97, 1","refit_cooks":true}'

Prefer the dataset's authoritative script

Before running DESeq2 yourself, ls the dataset folder. If you see run_*.py, analysis.R, find_*.R, or similar, those are the benchmark's ground-truth recipes.

1. cat the script to see its exact parameters — every kwarg. 2. If the script already prints the quantity you need, cd DATASET_DIR && python3 run_*.py and take its answer. 3. If the question needs a different metric from the same fitted model, make a SMALL addition (extra print statements) without changing the DeseqDataSet(...) / DeseqStats(...) constructor calls.

Copy ALL kwargs literally: refit_cooks=True, alpha=0.05, n_cpus, design_factors, ref_level. Omitting parameters like refit_cooks can change DEG counts significantly. Plain "remembered" defaults produce a different gene list than the ground-truth script.

R DESeq2 vs pydeseq2

The two libraries can disagree on edge cases (e.g., sig gene counts at the same alpha often differ by ~2%). Match whatever the authoritative script uses. If no script is present, prefer R DESeq2 — its behavior is more widely referenced in published papers.

Preferred: use the `run_deseq2_analysis` ToolUniverse tool which runs R DESeq2 via Rscript:

# Basic DESeq2
tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"metadata.csv","design":"~ condition","ref_level":"condition, Control"}'

# With contrast + LFC shrinkage + refit_cooks
tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"raw_counts.csv","metadata_file":"metadata.csv","design":"~ Replicate + Media + Strain","contrast":"Strain, 97, 1","refit_cooks":true,"lfc_shrinkage":true}'

# enrichGO + simplify (after saving DEG list to file)
tu run run_deseq2_analysis '{"operation":"enrichgo","gene_list_file":"sig_genes.txt","background_file":"all_genes.txt","simplify_cutoff":0.7}'

This avoids pydeseq2 vs R DESeq2 discrepancies. The tool returns sig gene counts, dispersion estimates, and a results CSV path for further analysis.

Strain identity pinning

When a question names strains by number AND describes their biology (e.g., knockout genotypes), pin the mapping from the question text or the experiment metadata, not from numeric order. Read the metadata CSV to confirm which strain number corresponds to which genotype before running DESeq2.

When reading RDS/CSV result files named like res_1vsN.rds, the "N" refers to the strain number in the metadata. Verify which mutant that number corresponds to — do not assume from the filename alone.

"Uniquely DE in A/B not C" = exclusive, not inclusive

When asked for genes "uniquely DE in one of {A, B} single mutants but not in C (double)", this means exclusively in A xor exclusively in B, each also not in C:

(A − B − C) ∪ (B − A − C)

NOT (A ∪ B) − C (which includes the A∩B intersection). The exclusive interpretation typically gives a smaller count than the inclusive one.

Set-operation percentages — check the denominator

If your unique_DE / total_DE gives a number in the 30–50% range, the expected denominator is probably different. Common alternatives:

• unique_DE / total_genes_tested (often 5-10% for bacterial, <2% for eukaryotic)
• unique_DE / |union of all sig sets|

Re-read the question to see what population "as a percentage of" refers to.

"Also DE in strain X" = simple overlap, not exclusive

When asked "what percentage of genes DE in A are also DE in B", compute |A ∩ B| / |A|. Do NOT subtract other strains — "also DE in B" does not mean "exclusively DE in B but not C".

CRITICAL: The word "also" means simple set intersection. If the question says "genes DE in strain A that are also DE in strain B", compute:

overlap = sigA.intersection(sigB)
pct = len(overlap) / len(sigA) * 100

Do NOT add extra exclusion conditions like "but not in strain C". "Also DE in B" means A ∩ B, nothing more.

Dispersion estimates: R DESeq2 vs pydeseq2 diverge

For questions about dispersion (e.g., "how many genes have dispersion below 1e-05"), R DESeq2 and pydeseq2 give different numbers because of implementation differences in the dispersion fitting algorithm. Always use R DESeq2 for dispersion questions — benchmark ground truths are computed with R:

library(DESeq2)
dds <- DESeq(dds)
# Pre-shrinkage (genewise) dispersions:
gene_disp <- mcols(dds)$dispGeneEst
cat("Below 1e-05:", sum(gene_disp < 1e-05, na.rm=TRUE), "\n")

Log2 fold-change: verify the contrast direction

When asked for the log2FC of gene X in mutant Y, verify that your DESeq2 results() call uses the correct contrast. For ~Media + Strain design with reference Strain "1":

• results(dds, contrast=c("Strain", "97", "1")) → log2FC for strain 97 vs ref
• The sign matters: negative log2FC = downregulated in mutant vs reference
• If your value has the wrong sign or magnitude, check if you accidentally used the wrong strain coefficient

Density / per-chromosome calculation conventions

When a question asks for "average chromosomal density", "genome-wide average density", or similar:

• "Average chromosomal density" = mean(per_chromosome_density) where per_chromosome_density = n_features / chromosome_length for each chromosome separately, then averaged across chromosomes (UNWEIGHTED mean).
• This is DIFFERENT from total_features / total_length (which is the WEIGHTED mean / "expected density under uniform distribution"). The latter is typically used as the chi-square test's expected value, not as the "average chromosomal density" being asked about.
• Example: in bird genomes with chromosomes ranging 6-120 Mb, mean(per_chrom_density) is ~2× larger than total/total because shorter chromosomes have proportionally more events per bp.

Data Loading and Validation

Detailed guide for loading count matrices and metadata.

Handling R Data Files (RDS)

If your data is in R format (.rds, .RData), convert to CSV first:

Option 1: Convert in R (if available)

# In R
result <- readRDS("deseq2_results.rds")
result_df <- as.data.frame(result)
write.csv(result_df, "deseq2_results.csv", row.names=TRUE)

Then load the CSV in Python as usual.

Option 2: Use rpy2 in Python (requires R installed)

import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
pandas2ri.activate()

# Read RDS file
ro.r(f'df <- as.data.frame(readRDS("{rds_file}"))')

# Convert to pandas
df = pandas2ri.rpy2py(ro.r['df'])

Note: This requires R and the relevant R packages (e.g., DESeq2) installed. If unavailable, use Option 1 on a machine with R.

Working with Pre-computed DESeq2 Results

If you have an RDS file containing DESeq2 results (not count data):

1. Convert to CSV (Option 1 or 2 above) 2. Filter genes directly in pandas:

# Read pre-computed results
results = pd.read_csv("deseq2_results.csv", index_col=0)

# Filter upregulated genes
upregulated = results[
    (results['log2FoldChange'] > 0) &  # Positive = upregulated
    (results['padj'] < 0.05)            # Significant
]

# Get gene list
genes = upregulated.index.tolist()

# Use with /tooluniverse-gene-enrichment skill

3. Skip to enrichment - Use gene list with enrichment analysis (Step 5 in workflow)

Load Count Matrix

import pandas as pd
import numpy as np
import os

def load_count_matrix(file_path, **kwargs):
    """Load count matrix from various formats.

    Expects: genes as rows/columns, samples as rows/columns.
    PyDESeq2 requires: samples as rows, genes as columns.
    """
    ext = os.path.splitext(file_path)[1].lower()

    if ext in ['.csv']:
        df = pd.read_csv(file_path, index_col=0, **kwargs)
    elif ext in ['.tsv', '.txt']:
        df = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
    elif ext in ['.h5ad']:
        import anndata
        adata = anndata.read_h5ad(file_path)
        df = pd.DataFrame(
            adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X,
            index=adata.obs_names,
            columns=adata.var_names
        )
        return df, adata.obs  # Return metadata too if available
    else:
        # Try tab-separated as default
        df = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)

    return df

Orient Matrix

CRITICAL: PyDESeq2 expects samples as rows, genes as columns.

def orient_count_matrix(df, metadata_samples=None):
    """Ensure samples are rows and genes are columns.

    Heuristic: if column count >> row count, genes are likely columns (correct).
    If row count >> column count, genes are likely rows (need transpose).
    If metadata_samples provided, match against index and columns.
    """
    if metadata_samples is not None:
        # Check if samples match rows or columns
        row_match = len(set(df.index) & set(metadata_samples))
        col_match = len(set(df.columns) & set(metadata_samples))
        if col_match > row_match:
            df = df.T
        return df

    # Heuristic: typical RNA-seq has 10-1000 samples and 10000-60000 genes
    if df.shape[0] > df.shape[1] * 5:  # Many more rows than columns
        df = df.T  # Transpose: genes were rows

    return df

Load Metadata

def load_metadata(file_path, **kwargs):
    """Load sample metadata (colData in R)."""
    ext = os.path.splitext(file_path)[1].lower()

    if ext in ['.csv']:
        meta = pd.read_csv(file_path, index_col=0, **kwargs)
    elif ext in ['.tsv', '.txt']:
        meta = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
    else:
        meta = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)

    return meta

Validate and Align

def validate_inputs(counts, metadata):
    """Validate count matrix and metadata alignment."""
    issues = []

    # Check sample alignment
    count_samples = set(counts.index)
    meta_samples = set(metadata.index)

    if count_samples != meta_samples:
        common = count_samples & meta_samples
        if len(common) == 0:
            # Try matching columns
            if set(counts.columns) & meta_samples:
                counts = counts.T
                count_samples = set(counts.index)
                common = count_samples & meta_samples

        if len(common) > 0:
            counts = counts.loc[sorted(common)]
            metadata = metadata.loc[sorted(common)]
            issues.append(f"Aligned to {len(common)} common samples")
        else:
            issues.append("ERROR: No matching samples between counts and metadata")
            return None, None, issues

    # Ensure integer counts
    if counts.dtypes.apply(lambda x: x == float).any():
        if (counts % 1 == 0).all().all():
            counts = counts.astype(int)
        else:
            # Might be normalized data - round to integers for DESeq2
            issues.append("WARNING: Non-integer counts detected. Rounding to integers.")
            counts = counts.round().astype(int)

    # Remove genes with zero counts across all samples
    nonzero_mask = counts.sum(axis=0) > 0
    n_removed = (~nonzero_mask).sum()
    if n_removed > 0:
        counts = counts.loc[:, nonzero_mask]
        issues.append(f"Removed {n_removed} genes with zero counts across all samples")

    # Remove negative values
    if (counts < 0).any().any():
        issues.append("WARNING: Negative counts detected. Setting to 0.")
        counts = counts.clip(lower=0)

    return counts, metadata, issues

Subset Samples

def subset_samples(counts, metadata, condition_col, values=None, exclude=None):
    """Subset samples based on metadata conditions."""
    if values is not None:
        mask = metadata[condition_col].isin(values)
    elif exclude is not None:
        mask = ~metadata[condition_col].isin(exclude)
    else:
        return counts, metadata

    metadata = metadata[mask]
    counts = counts.loc[metadata.index]
    return counts, metadata

Example: Load and Validate

# Load data
counts = load_count_matrix("counts.csv")
metadata = load_metadata("metadata.csv")

# Orient if needed
counts = orient_count_matrix(counts, metadata.index)

# Validate and align
counts, metadata, issues = validate_inputs(counts, metadata)

# Print validation issues
for issue in issues:
    print(f"  {issue}")

# Subset if needed
counts, metadata = subset_samples(
    counts, metadata,
    condition_col='treatment',
    values=['control', 'treated']
)

print(f"\nFinal dimensions:")
print(f"  Counts: {counts.shape[0]} samples × {counts.shape[1]} genes")
print(f"  Metadata: {metadata.shape[0]} samples")

Handling Different Input Formats

CSV with genes as columns (correct)

,Gene1,Gene2,Gene3
Sample1,100,50,200
Sample2,120,60,180

No action needed.

CSV with genes as rows (needs transpose)

,Sample1,Sample2,Sample3
Gene1,100,120,90
Gene2,50,60,45

Call orient_count_matrix().

H5AD (AnnData)

import anndata
adata = anndata.read_h5ad("data.h5ad")
counts = pd.DataFrame(
    adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X,
    index=adata.obs_names,
    columns=adata.var_names
)
metadata = adata.obs

Pre-normalized data (FPKM, TPM)

If data is NOT raw counts:

# Option 1: Round to integers (acceptable for DESeq2)
counts = counts.round().astype(int)

# Option 2: Use t-test instead of DESeq2 (for normalized data)
from scipy import stats
stat, pval = stats.ttest_ind(group1, group2)

Common Issues

Dispersion Analysis and Diagnostics

Understanding and analyzing dispersion estimates in PyDESeq2.

Dispersion Storage in PyDESeq2

CRITICAL: In PyDESeq2, dispersions are stored in dds.var (NOT dds.varm).

# After running dds.deseq2()
disp_data = dds.var  # DataFrame with dispersion columns

Dispersion Columns

def get_dispersion_data(dds):
    """Extract all dispersion-related data from fitted DESeq2 model.

    Returns dict with:
        - genewise_dispersions: Per-gene maximum likelihood estimates
        - fitted_dispersions: Trend-fitted values
        - MAP_dispersions: Maximum a posteriori (after shrinkage to trend)
        - dispersions: Final dispersions used in testing
    """
    disp_data = {}

    # Key dispersion columns in dds.var
    disp_columns = [
        'genewise_dispersions',  # Pre-shrinkage (gene-wise MLE)
        'fitted_dispersions',     # Trend curve values
        'MAP_dispersions',        # Post-shrinkage (MAP estimates)
        'dispersions',            # Final (MAP or genewise for outliers)
    ]

    for col in disp_columns:
        if col in dds.var.columns:
            disp_data[col] = dds.var[col]

    return disp_data

Question Phrasing to Dispersion Column Mapping

Dispersion Diagnostics

def dispersion_diagnostics(dds, threshold=1e-5):
    """Analyze dispersion estimates for diagnostics.

    Common question: "How many genes have a dispersion estimate
    below 1e-05 prior to dispersion fitting?"

    Answer: Count genewise_dispersions < threshold
    """
    disp_data = get_dispersion_data(dds)

    diagnostics = {}

    if 'genewise_dispersions' in disp_data:
        gwd = disp_data['genewise_dispersions']
        diagnostics['genewise_below_threshold'] = (gwd < threshold).sum()
        diagnostics['genewise_min'] = gwd.min()
        diagnostics['genewise_max'] = gwd.max()
        diagnostics['genewise_median'] = gwd.median()
        diagnostics['genewise_mean'] = gwd.mean()

    if 'fitted_dispersions' in disp_data:
        fd = disp_data['fitted_dispersions']
        diagnostics['fitted_below_threshold'] = (fd < threshold).sum()
        diagnostics['fitted_min'] = fd.min()
        diagnostics['fitted_max'] = fd.max()
        diagnostics['fitted_median'] = fd.median()

    if 'MAP_dispersions' in disp_data:
        mapd = disp_data['MAP_dispersions']
        diagnostics['MAP_below_threshold'] = (mapd < threshold).sum()
        diagnostics['MAP_min'] = mapd.min()
        diagnostics['MAP_max'] = mapd.max()
        diagnostics['MAP_median'] = mapd.median()

    if 'dispersions' in disp_data:
        d = disp_data['dispersions']
        diagnostics['final_below_threshold'] = (d < threshold).sum()

    return diagnostics

Common Questions

Count genes below threshold

# "How many genes have dispersion below 1e-5 prior to fitting?"
genewise = dds.var['genewise_dispersions']
answer = (genewise < 1e-5).sum()

Count genes after shrinkage

# "How many genes have dispersion below 1e-5 after shrinkage?"
map_disp = dds.var['MAP_dispersions']
answer = (map_disp < 1e-5).sum()

Range of dispersions

# "What is the range of gene-wise dispersions?"
genewise = dds.var['genewise_dispersions']
min_disp = genewise.min()
max_disp = genewise.max()
answer = f"{min_disp:.2E} to {max_disp:.2E}"

Median dispersion

# "What is the median dispersion estimate?"
median_disp = dds.var['dispersions'].median()
answer = f"{median_disp:.2E}"

Dispersion Shrinkage Effect

def analyze_shrinkage_effect(dds):
    """Compare gene-wise vs MAP dispersions to assess shrinkage."""
    genewise = dds.var['genewise_dispersions']
    map_disp = dds.var['MAP_dispersions']

    # Genes where shrinkage reduced dispersion
    shrunk_genes = (map_disp < genewise).sum()

    # Genes where shrinkage increased dispersion
    expanded_genes = (map_disp > genewise).sum()

    # Median fold change
    fold_change = map_disp / genewise
    median_fc = fold_change.median()

    results = {
        'shrunk_genes': shrunk_genes,
        'expanded_genes': expanded_genes,
        'median_fold_change': median_fc,
        'mean_genewise': genewise.mean(),
        'mean_MAP': map_disp.mean()
    }

    return results

Outlier Detection

def identify_dispersion_outliers(dds, threshold=10):
    """Identify genes with outlier dispersions.

    Outliers are genes where genewise dispersion is far from fitted.
    """
    genewise = dds.var['genewise_dispersions']
    fitted = dds.var['fitted_dispersions']

    # Ratio of genewise to fitted
    ratio = genewise / fitted

    # Outliers: ratio > threshold
    outliers = ratio > threshold
    outlier_genes = dds.var.index[outliers]

    return outlier_genes.tolist()

Dispersion vs Mean Expression

import matplotlib.pyplot as plt
import numpy as np

def plot_dispersion_trend(dds):
    """Plot dispersion vs mean expression (dispersion plot)."""
    baseMean = dds.var['baseMean']
    genewise = dds.var['genewise_dispersions']
    fitted = dds.var['fitted_dispersions']
    map_disp = dds.var['MAP_dispersions']

    plt.figure(figsize=(10, 6))
    plt.scatter(baseMean, genewise, s=1, alpha=0.3, label='Gene-wise', color='gray')
    plt.scatter(baseMean, map_disp, s=1, alpha=0.5, label='MAP', color='blue')

    # Plot fitted curve
    sorted_idx = baseMean.argsort()
    plt.plot(baseMean.iloc[sorted_idx], fitted.iloc[sorted_idx],
             color='red', linewidth=2, label='Fitted trend')

    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('Base Mean Expression')
    plt.ylabel('Dispersion')
    plt.legend()
    plt.title('Dispersion Estimates')
    plt.tight_layout()
    plt.savefig('dispersion_plot.png', dpi=150)
    plt.close()

Dispersion Fitting Convergence

def check_dispersion_convergence(dds):
    """Check if dispersion fitting converged properly."""
    fitted = dds.var['fitted_dispersions']

    # Check for NaN or Inf values
    has_nan = fitted.isna().any()
    has_inf = np.isinf(fitted).any()

    # Check if fitted values are reasonable
    min_fitted = fitted.min()
    max_fitted = fitted.max()

    convergence_ok = not (has_nan or has_inf) and (min_fitted > 0) and (max_fitted < 1e6)

    return {
        'converged': convergence_ok,
        'has_nan': has_nan,
        'has_inf': has_inf,
        'min_fitted': min_fitted,
        'max_fitted': max_fitted
    }

Alternative Dispersion Fitting

# If parametric fitting fails, use mean fitting
try:
    dds = DeseqDataSet(counts=counts, metadata=metadata,
                      design="~condition", fit_type='parametric', quiet=True)
    dds.deseq2()
except Exception as e:
    print(f"Parametric fit failed: {e}")
    print("Retrying with fit_type='mean'")
    dds = DeseqDataSet(counts=counts, metadata=metadata,
                      design="~condition", fit_type='mean', quiet=True)
    dds.deseq2()

Dispersion by Gene Expression Level

def dispersion_by_expression_level(dds, quantiles=[0.25, 0.5, 0.75]):
    """Analyze dispersion across expression level quantiles."""
    baseMean = dds.var['baseMean']
    dispersions = dds.var['dispersions']

    results = {}
    for q in quantiles:
        threshold = baseMean.quantile(q)
        genes_below = baseMean <= threshold
        median_disp = dispersions[genes_below].median()
        results[f'Q{int(q*100)}'] = {
            'expression_threshold': threshold,
            'median_dispersion': median_disp,
            'n_genes': genes_below.sum()
        }

    return results

Complete Example

from pydeseq2.dds import DeseqDataSet

# Run DESeq2
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", quiet=True)
dds.deseq2()

# Question: "How many genes have dispersion below 1e-5 prior to fitting?"
genewise = dds.var['genewise_dispersions']
answer = (genewise < 1e-5).sum()
print(f"Genes with dispersion < 1e-5 (prior to fitting): {answer}")

# Additional diagnostics
diag = dispersion_diagnostics(dds, threshold=1e-5)
print(f"\nDispersion diagnostics:")
print(f"  Gene-wise below threshold: {diag['genewise_below_threshold']}")
print(f"  Gene-wise median: {diag['genewise_median']:.2E}")
print(f"  MAP below threshold: {diag['MAP_below_threshold']}")
print(f"  MAP median: {diag['MAP_median']:.2E}")

# Shrinkage effect
shrinkage = analyze_shrinkage_effect(dds)
print(f"\nShrinkage effect:")
print(f"  Genes with reduced dispersion: {shrinkage['shrunk_genes']}")
print(f"  Median fold change: {shrinkage['median_fold_change']:.2f}")

# Check convergence
conv = check_dispersion_convergence(dds)
print(f"\nDispersion fitting convergence: {'OK' if conv['converged'] else 'FAILED'}")

Worked Examples

Pattern: Count below threshold prior to fitting

# Question: "How many genes have a dispersion estimate below 1e-05 prior to dispersion fitting?"
genewise = dds.var['genewise_dispersions']
answer = (genewise < 1e-5).sum()

Pattern: Median dispersion

# Question: "What is the median dispersion after shrinkage?"
map_disp = dds.var['MAP_dispersions']
answer = round(map_disp.median(), 6)

Pattern: Range of dispersions

# Question: "What is the minimum gene-wise dispersion?"
genewise = dds.var['genewise_dispersions']
answer = f"{genewise.min():.2E}"

edgeR and limma-voom — alternative DE routes to DESeq2

DESeq2 is this skill's default route, but edgeR and limma-voom are the two other standard bulk RNA-seq differential-expression frameworks. Published pipelines routinely route across all three. This doc gives the concrete R command sequences, the input/output contracts, and how to read each framework's output relative to DESeq2's.

The bundled scripts/r_edger_limma_wrapper.py runs both of these for you with the same workspace-isolation and parseable-output conventions as r_deseq2_wrapper.py. Use the wrapper first; the raw command sequences below are for when you need a variant the wrapper does not expose.

---

When to use which (see also the routing subsection in SKILL.md)

These are reasoned defaults, not hard rules. When an authoritative script or executed notebook in the data folder already ran one framework, match it — the published ground-truth number comes from whichever the pipeline used. The three frameworks usually agree on the strongly-DE genes but differ by a few percent on borderline counts at the same threshold.

---

edgeR — quasi-likelihood F-test (QL-F) pipeline

Recommended modern edgeR route (preferred over the classic exact test for anything beyond a single 2-group comparison):

library(edgeR)
# counts: integer matrix, genes x samples. group/design from metadata.
dge <- DGEList(counts = counts_int)
design <- model.matrix(~ condition, data = metadata)   # or ~ batch + condition
keep <- filterByExpr(dge, design)                       # standard low-count filter
dge  <- dge[keep, , keep.lib.sizes = FALSE]
dge  <- calcNormFactors(dge)                            # TMM normalization
dge  <- estimateDisp(dge, design)                       # NB dispersions (trended + tagwise)
fit  <- glmQLFit(dge, design)                           # quasi-likelihood GLM
qlf  <- glmQLFTest(fit, coef = "conditiontreated")      # test one coefficient
# ...or an explicit contrast between two non-reference levels:
# con <- makeContrasts(grpB - grpC, levels = design)
# qlf <- glmQLFTest(fit, contrast = con)
res  <- topTags(qlf, n = Inf, sort.by = "PValue")$table

Classic exact test (only valid for a single one-way grouping, no covariates):

dge <- DGEList(counts = counts_int, group = metadata$condition)
dge <- calcNormFactors(dge)
dge <- estimateDisp(dge)
et  <- exactTest(dge, pair = c("control", "treated"))
res <- topTags(et, n = Inf)$table

edgeR output columns: logFC, logCPM, F (QL-F) or logFC/logCPM/PValue (exact test), PValue, FDR.

---

limma-voom pipeline

library(limma); library(edgeR)
dge <- DGEList(counts = counts_int)
design <- model.matrix(~ batch + condition, data = metadata)
keep <- filterByExpr(dge, design)
dge  <- dge[keep, , keep.lib.sizes = FALSE]
dge  <- calcNormFactors(dge)                 # TMM (from edgeR)
v    <- voom(dge, design)                    # mean-variance weights -> logCPM
fit  <- lmFit(v, design)
# single coefficient:
fit  <- eBayes(fit)
res  <- topTable(fit, coef = "conditiontreated", number = Inf, sort.by = "P")
# ...or an explicit contrast:
# cm  <- makeContrasts(conditiontreated, levels = design)
# fit <- eBayes(contrasts.fit(lmFit(v, design), cm))
# res <- topTable(fit, coef = 1, number = Inf, sort.by = "P")

For repeated measures / paired designs, estimate the intra-block correlation and pass it to both voom and lmFit:

corfit <- duplicateCorrelation(v, design, block = metadata$subject)
v      <- voom(dge, design, block = metadata$subject, correlation = corfit$consensus)
fit    <- lmFit(v, design, block = metadata$subject, correlation = corfit$consensus)

limma-voom output columns: logFC, AveExpr, t, P.Value, adj.P.Val (BH-adjusted), B (log-odds of DE).

---

Input contract (shared with the DESeq2 wrapper)

• Count matrix: CSV, genes as rows, samples as columns, first column = gene

IDs. RAW integer counts (NOT normalized/TPM/CPM). The wrapper rounds and integer-coerces; if your matrix is already normalized, edgeR/limma results will be wrong — supply raw counts.

• Sample metadata: CSV, one row per sample. The wrapper auto-detects the

sample-name column (AzentaName, sample, SampleID, sample_id, SampleName, projid, else first column) and aligns it to the count columns.

• Design: an R formula string, e.g. ~condition or ~batch + condition.

Factor of interest typically last.

• Contrast: factor,level1,level2 meaning level1 vs level2. If level2 is

the model's reference level, the wrapper tests the single <factor><level1> coefficient; otherwise it builds an explicit +1/-1 contrast vector between the two non-reference columns.

Output contract

r_edger_limma_wrapper.py writes one ranked CSV to --workdir (res_<method>_<label>.csv) and prints parseable lines:

# METHOD edger|limma
# CONTRAST <factor>_<lvl1>_vs_<lvl2>: n_genes=<after filterByExpr> n_tested=<x>
# SIG_<label>_padj_only (FDR<thr): n=...
# SIG_<label>_padjlfc (FDR<thr AND |logFC|>thr): n=...
# SIG_<label>_strict (FDR<thr AND |logFC|>thr AND logCPM/AveExpr>thr): n=...
# GENE <name> [<label>]: logFC=... FDR=...
# TABLE <abs path to ranked CSV>

The three SIG_* lines mirror r_deseq2_wrapper.py exactly, so DEG counts are directly comparable across the DESeq2 / edgeR / limma-voom routes. Default to the _padj_only line unless the question names an |logFC| or expression threshold (same rule as the DESeq2 route).

---

Column-name crosswalk (edgeR / limma vs DESeq2)

Interpretation notes

• All three logFC/log2FoldChange are log2; sign convention is level1 vs

level2 (positive = up in level1). edgeR/limma do NOT shrink the logFC by default, so their logFC is comparable to DESeq2's UNSHRUNKEN log2FoldChange, not the apeglm-shrunken value. For individual low-count gene queries, the edgeR/limma logFC behaves like the unshrunken DESeq2 LFC.

• baseMean (DESeq2, linear scale) and logCPM/AveExpr (log2 scale) are NOT

the same units — do not threshold them with the same number. A baseMean>10 filter has no direct edgeR/limma equivalent; filterByExpr already removes low-count genes upstream, which is the edgeR/limma analogue of DESeq2's independent filtering.

• "Significant DEGs" defaults to FDR < 0.05 / adj.P.Val < 0.05 —

the _padj_only line — unless the question adds a fold-change or expression threshold.

Enrichment Analysis with gseapy

Complete guide to pathway and GO enrichment analysis.

Basic Over-Representation Analysis (ORA)

import gseapy as gp

# Prepare gene list (from DESeq2 results)
sig_genes = results[(results['padj'] < 0.05) & (results['log2FoldChange'].abs() > 0.5)]
gene_list = sig_genes.index.tolist()

# Run enrichment
enr = gp.enrich(
    gene_list=gene_list,
    gene_sets='GO_Biological_Process_2023',
    background=None,  # or provide background gene list
    outdir=None,  # Don't save files
    cutoff=0.05,
    no_plot=True,
    verbose=False
)

# Access results
results_df = enr.results
print(results_df.head(10))

Gene Set Library Selection

Gene Ontology (GO)

# Human/Mouse GO Biological Process (most recent)
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023')

# GO Molecular Function
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Molecular_Function_2021')

# GO Cellular Component
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Cellular_Component_2021')

KEGG Pathways

# Human KEGG
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2021_Human')

# Mouse KEGG
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2019_Mouse')

Reactome

enr = gp.enrich(gene_list=gene_list, gene_sets='Reactome_2022')

WikiPathways

# Human
enr = gp.enrich(gene_list=gene_list, gene_sets='WikiPathways_2019_Human')

# Mouse
enr = gp.enrich(gene_list=gene_list, gene_sets='WikiPathways_2019_Mouse')

Other Libraries

# MSigDB Hallmark gene sets
enr = gp.enrich(gene_list=gene_list, gene_sets='MSigDB_Hallmark_2020')

# GWAS Catalog
enr = gp.enrich(gene_list=gene_list, gene_sets='GWAS_Catalog_2019')

# BioCarta
enr = gp.enrich(gene_list=gene_list, gene_sets='BioCarta_2016')

Using Background Gene Sets

# Background = all genes tested in DESeq2
background = results.dropna(subset=['padj']).index.tolist()

enr = gp.enrich(
    gene_list=gene_list,
    gene_sets='GO_Biological_Process_2023',
    background=background,  # Provide background
    outdir=None,
    cutoff=0.05
)

Extract Specific Results

def extract_enrichment_answer(enr_results, term_query=None, metric='Adjusted P-value'):
    """Extract specific enrichment result.

    Args:
        enr_results: gseapy enrichment results DataFrame
        term_query: String to search in Term column (case-insensitive)
        metric: Column to return ('Adjusted P-value', 'Odds Ratio', 'P-value', etc.)

    Returns:
        Value or DataFrame of matches
    """
    if term_query:
        # Case-insensitive search
        mask = enr_results['Term'].str.lower().str.contains(term_query.lower())
        matches = enr_results[mask]
        if len(matches) == 1:
            return matches.iloc[0][metric]
        elif len(matches) > 1:
            return matches[['Term', metric]]
        else:
            return None

    # Return top result
    return enr_results.sort_values(metric).head(1)

# Usage
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2021_Human')
pval = extract_enrichment_answer(enr.results, term_query='ABC transporters', metric='Adjusted P-value')
print(f"ABC transporters adjusted p-value: {pval}")

Extract Gene Count in Pathway

# Enrichment results have 'Overlap' column (e.g., "11/42")
# Also 'Genes' column with semicolon-separated gene list

pathway_row = enr.results[enr.results['Term'].str.contains('ABC transporters')].iloc[0]

# Number of DEGs in pathway
overlap_str = pathway_row['Overlap']  # e.g., "11/42"
n_overlap = int(overlap_str.split('/')[0])  # 11

# Pathway size
n_pathway = int(overlap_str.split('/')[1])  # 42

# Gene list
genes_in_pathway = pathway_row['Genes'].split(';')
print(f"{n_overlap} genes contribute to this pathway:")
print(genes_in_pathway)

Multi-Library Enrichment

# Run enrichment on multiple libraries
libraries = [
    'GO_Biological_Process_2023',
    'KEGG_2021_Human',
    'Reactome_2022'
]

all_results = {}
for lib in libraries:
    enr = gp.enrich(
        gene_list=gene_list,
        gene_sets=lib,
        outdir=None,
        cutoff=0.05,
        no_plot=True,
        verbose=False
    )
    all_results[lib] = enr.results

# Combine top results
combined = []
for lib, res in all_results.items():
    top5 = res.head(5).copy()
    top5['Library'] = lib
    combined.append(top5)
combined_df = pd.concat(combined, ignore_index=True)

GO Term Simplification

def simplify_go_terms(enr_results, similarity_threshold=0.7):
    """Simplify GO terms by removing highly similar terms.

    Approximation of R clusterProfiler::simplify().
    Uses Jaccard similarity on gene sets.
    """
    if len(enr_results) == 0:
        return enr_results

    # Parse gene sets from Genes column
    terms = enr_results.sort_values('Adjusted P-value').copy()
    gene_sets = {}
    for _, row in terms.iterrows():
        genes = set(row['Genes'].split(';'))
        gene_sets[row['Term']] = genes

    # Compute Jaccard similarity between terms
    keep = []
    removed = set()

    for i, (term_i, genes_i) in enumerate(gene_sets.items()):
        if term_i in removed:
            continue
        keep.append(term_i)

        for term_j, genes_j in list(gene_sets.items())[i+1:]:
            if term_j in removed:
                continue
            # Jaccard similarity
            intersection = len(genes_i & genes_j)
            union = len(genes_i | genes_j)
            if union > 0:
                similarity = intersection / union
                if similarity > similarity_threshold:
                    removed.add(term_j)  # Remove the less significant term

    return terms[terms['Term'].isin(keep)]

# Usage
enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023')
simplified = simplify_go_terms(enr.results, similarity_threshold=0.7)
print(f"Original: {len(enr.results)} terms")
print(f"Simplified: {len(simplified)} terms")

Gene Set Enrichment Analysis (GSEA)

# GSEA requires ranked gene list (not just significant genes)
# Rank by -log10(pvalue) * sign(log2FC)

results_ranked = results.dropna(subset=['pvalue', 'log2FoldChange'])
results_ranked['rank'] = -np.log10(results_ranked['pvalue']) * np.sign(results_ranked['log2FoldChange'])
results_ranked = results_ranked.sort_values('rank', ascending=False)

# Create rank dictionary
rank_dict = dict(zip(results_ranked.index, results_ranked['rank']))

# Run GSEA
gsea_res = gp.prerank(
    rnk=rank_dict,
    gene_sets='KEGG_2021_Human',
    outdir=None,
    permutation_num=1000,
    no_plot=True,
    verbose=False
)

# Access results
gsea_df = gsea_res.res2d
print(gsea_df[gsea_df['FDR q-val'] < 0.05])

Organism-Specific Libraries

Human

libraries_human = [
    'GO_Biological_Process_2023',
    'GO_Molecular_Function_2021',
    'GO_Cellular_Component_2021',
    'KEGG_2021_Human',
    'Reactome_2022',
    'WikiPathways_2019_Human',
    'MSigDB_Hallmark_2020',
    'BioCarta_2016'
]

Mouse

libraries_mouse = [
    'GO_Biological_Process_2023',
    'GO_Molecular_Function_2021',
    'GO_Cellular_Component_2021',
    'KEGG_2019_Mouse',
    'WikiPathways_2019_Mouse'
]

Other Organisms

For other organisms, use custom gene sets:

# Load custom GMT file
enr = gp.enrich(
    gene_list=gene_list,
    gene_sets='/path/to/custom.gmt',
    background=background
)

Complete Example: DEG to Enrichment

import pandas as pd
import gseapy as gp
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# Run DESeq2 (from pydeseq2_workflow.md)
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", quiet=True)
dds.deseq2()
stat_res = DeseqStats(dds, contrast=['condition', 'treatment', 'control'], quiet=True)
stat_res.run_wald_test()
results = stat_res.results_df

# Filter DEGs
sig_genes = results[(results['padj'] < 0.05) & (results['log2FoldChange'].abs() > 0.5)]
gene_list = sig_genes.index.tolist()

# Run GO enrichment
enr_go = gp.enrich(
    gene_list=gene_list,
    gene_sets='GO_Biological_Process_2023',
    background=results.dropna(subset=['padj']).index.tolist(),
    outdir=None,
    cutoff=0.05,
    no_plot=True,
    verbose=False
)

# Run KEGG enrichment
enr_kegg = gp.enrich(
    gene_list=gene_list,
    gene_sets='KEGG_2021_Human',
    background=results.dropna(subset=['padj']).index.tolist(),
    outdir=None,
    cutoff=0.05,
    no_plot=True,
    verbose=False
)

# Display top results
print("\nTop 5 GO terms:")
print(enr_go.results[['Term', 'Adjusted P-value', 'Overlap']].head(5))

print("\nTop 5 KEGG pathways:")
print(enr_kegg.results[['Term', 'Adjusted P-value', 'Overlap']].head(5))

# Answer specific question
if 'immune response' in question.lower():
    immune_result = enr_go.results[enr_go.results['Term'].str.contains('immune', case=False)]
    if len(immune_result) > 0:
        answer = immune_result.iloc[0]['Adjusted P-value']
        print(f"\nImmune response adjusted p-value: {answer}")

Common Enrichment Patterns

Pattern 1: Extract adjusted p-value for specific pathway

enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2021_Human')
pathway = enr.results[enr.results['Term'].str.contains('ABC transporters')]
answer = round(pathway.iloc[0]['Adjusted P-value'], 4)

Pattern 2: Count significant pathways

enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023', cutoff=0.05)
answer = len(enr.results[enr.results['Adjusted P-value'] < 0.05])

Pattern 3: Gene count in pathway

pathway = enr.results[enr.results['Term'].str.contains('ribosome')]
overlap = pathway.iloc[0]['Overlap']  # e.g., "25/150"
answer = int(overlap.split('/')[0])  # 25

Pattern 4: Simplify GO terms

enr = gp.enrich(gene_list=gene_list, gene_sets='GO_Biological_Process_2023')
simplified = simplify_go_terms(enr.results, similarity_threshold=0.7)
answer = len(simplified)

Output Formatting Guide

Match question's requested format exactly.

Numeric Precision

Decimal Places

# "rounded to 2 decimal points"
answer = round(value, 2)

# "rounded to 4 decimal points"
answer = round(value, 4)

# "rounded to 1 decimal place"
answer = round(value, 1)

Scientific Notation

# "in scientific notation"
answer = f"{value:.2E}"  # e.g., "1.23E-05"

# With specific precision
answer = f"{value:.3E}"  # e.g., "1.234E-05"

# Standard form (lowercase e)
answer = f"{value:.2e}"  # e.g., "1.23e-05"

Integer Counts

# "how many genes"
answer = int(len(sig_genes))  # No decimals

# Ensure integer type
answer = len(sig_genes)  # Already int from len()

Percentages

# "as a percentage"
answer = f"{value * 100:.1f}%"  # e.g., "15.3%"

# Without percent sign
answer = round(value * 100, 1)  # e.g., 15.3

# "as a percentage rounded to 2 decimal points"
answer = round(value * 100, 2)  # e.g., 15.32

Ratios and Fractions

# "as a ratio X:Y"
answer = f"{x}:{y}"  # e.g., "3:1"

# "as a fraction"
answer = f"{numerator}/{denominator}"  # e.g., "11/42"

Ranges

# For range_verifier evaluation mode
# Expected: (700, 1000)
# Your answer: 842  -> PASS (within range)
answer = 842

# Confidence intervals
ci_low, ci_high = proportion_confint(n_sig, n_total, method='wilson')
answer = (round(ci_low, 2), round(ci_high, 2))  # e.g., (0.15, 0.25)

Lists

# "list the top 5 genes"
answer = results.sort_values('padj').head(5).index.tolist()

# Comma-separated string
answer = ", ".join(gene_list)  # e.g., "TP53, BRCA1, EGFR"

Boolean/Categorical

# "Yes" or "No"
answer = "Yes" if condition else "No"

# "increases" or "decreases"
if n_clean > n_all:
    answer = "Increases the number of differentially expressed genes"
else:
    answer = "Decreases the number of differentially expressed genes"

Tables

# DataFrame subset
answer = results[['log2FoldChange', 'padj']].head(10)

# Markdown table
def to_markdown_table(df):
    """Convert DataFrame to markdown table."""
    lines = []
    lines.append("| " + " | ".join(df.columns) + " |")
    lines.append("| " + " | ".join(["---"] * len(df.columns)) + " |")
    for _, row in df.iterrows():
        lines.append("| " + " | ".join(str(v) for v in row) + " |")
    return "\n".join(lines)

Common Answer Formats

Count (integer)

# "How many genes..."
answer = len(sig_genes)

P-value (4 decimals)

# "What is the adjusted p-value..."
answer = round(pathway.iloc[0]['Adjusted P-value'], 4)

Log2FC (2 decimals)

# "What is the log2 fold change..."
answer = round(results.loc['TP53', 'log2FoldChange'], 2)

Percentage (1 decimal)

# "What percentage of..."
percentage = len(overlap) / len(degs_A) * 100
answer = round(percentage, 1)

Scientific notation (p-values)

# Very small p-values
answer = f"{pvalue:.2E}"  # "1.23E-10"

Confidence interval (2 decimals)

ci_low, ci_high = proportion_confint(n_sig, n_total, method='wilson')
answer = f"({round(ci_low, 2)}, {round(ci_high, 2)})"

Format Detection

def format_answer(value, question_text):
    """Auto-detect format from question text."""
    question_lower = question_text.lower()

    # Decimal places
    if "2 decimal" in question_lower:
        return round(value, 2)
    elif "4 decimal" in question_lower:
        return round(value, 4)
    elif "1 decimal" in question_lower:
        return round(value, 1)

    # Scientific notation
    if "scientific notation" in question_lower:
        return f"{value:.2E}"

    # Percentage
    if "percentage" in question_lower:
        if isinstance(value, float) and value < 1:
            return round(value * 100, 1)
        return round(value, 1)

    # Count
    if "how many" in question_lower:
        return int(value)

    # Default
    return value

Examples by Question Type

"How many genes show significant DE?"

sig_genes = results[(results['padj'] < 0.05) & (results['log2FoldChange'].abs() > 0.5)]
answer = len(sig_genes)  # Integer: 842

"What is the adjusted p-value for pathway X?"

pathway = enr.results[enr.results['Term'].str.contains('ABC transporters')]
answer = round(pathway.iloc[0]['Adjusted P-value'], 4)  # 0.0234

"What percentage of DEGs in A are also in B?"

degs_A = set(results_A[results_A['padj'] < 0.05].index)
degs_B = set(results_B[results_B['padj'] < 0.05].index)
overlap = degs_A & degs_B
answer = round(len(overlap) / len(degs_A) * 100, 1)  # 15.3

"What is the 95% CI for the proportion?"

ci_low, ci_high = proportion_confint(n_sig, n_total, method='wilson')
answer = (round(ci_low, 2), round(ci_high, 2))  # (0.15, 0.25)

"What is the log2FC of gene TP53?"

answer = round(results.loc['TP53', 'log2FoldChange'], 2)  # 2.45

"What is the median dispersion?"

answer = f"{dds.var['dispersions'].median():.2E}"  # 1.23E-02

Avoiding Common Mistakes

# WRONG: Float when count expected
answer = 842.0  # Bad for "how many"

# RIGHT: Integer
answer = 842  # Good

# WRONG: Too many decimals
answer = 0.0234567  # Bad if question asks for 4 decimals

# RIGHT: Correct precision
answer = 0.0235  # Good

# WRONG: Scientific notation when not requested
answer = "8.42E+02"  # Bad for count

# RIGHT: Regular number
answer = 842  # Good

# WRONG: Missing percent sign
answer = 15.3  # Ambiguous

# RIGHT: Clear percentage
answer = "15.3%"  # OR answer = 15.3 if question says "as percentage"

Validation Before Submission

def validate_answer_format(answer, expected_type):
    """Validate answer format matches expected type."""
    if expected_type == 'integer':
        assert isinstance(answer, int), f"Expected int, got {type(answer)}"
    elif expected_type == 'float':
        assert isinstance(answer, (int, float)), f"Expected number, got {type(answer)}"
    elif expected_type == 'percentage':
        if isinstance(answer, str):
            assert '%' in answer, "Percentage should include % sign"
        else:
            assert 0 <= answer <= 100, "Percentage should be 0-100"
    elif expected_type == 'scientific':
        assert 'E' in str(answer) or 'e' in str(answer), "Should be in scientific notation"

    return True

PyDESeq2 Complete Workflow

Comprehensive code examples for PyDESeq2 analysis.

Basic Single-Factor Analysis

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

# Assume counts and metadata are already loaded and validated

# Set reference level (first category becomes reference)
metadata['condition'] = pd.Categorical(
    metadata['condition'],
    categories=['control', 'treatment']  # control is reference
)

# Create DESeq2 dataset
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~condition",
    quiet=True
)

# Run DESeq2 pipeline (normalization + dispersion + testing)
dds.deseq2()

# Extract results for contrast
stat_res = DeseqStats(
    dds,
    contrast=['condition', 'treatment', 'control'],
    alpha=0.05,
    quiet=True
)
stat_res.run_wald_test()
stat_res.summary()

# Get results DataFrame
results = stat_res.results_df

Multi-Factor Design

# Design with multiple factors
metadata['strain'] = pd.Categorical(metadata['strain'], categories=['WT', 'mutant'])
metadata['media'] = pd.Categorical(metadata['media'], categories=['LB', 'M9'])

# Multi-factor design
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~strain + media",  # Additive model
    quiet=True
)
dds.deseq2()

# Extract results for strain effect (controlling for media)
stat_res = DeseqStats(dds, contrast=['strain', 'mutant', 'WT'], quiet=True)
stat_res.run_wald_test()
stat_res.summary()
results_strain = stat_res.results_df

Interaction Design

# Design with interaction term
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~strain + treatment + strain:treatment",
    quiet=True
)
dds.deseq2()

# Main effect of treatment
stat_res_treatment = DeseqStats(dds, contrast=['treatment', 'treated', 'control'], quiet=True)
stat_res_treatment.run_wald_test()

# Interaction effect (requires coefficient name)
# Check available coefficients
print(dds.varm['LFC'].columns)  # View all coefficients

# Extract interaction results
stat_res_interaction = DeseqStats(dds, quiet=True)
# Note: interaction testing requires specifying coefficient directly

Batch Effect Correction

# Include batch as covariate
metadata['batch'] = pd.Categorical(metadata['batch'])
metadata['condition'] = pd.Categorical(metadata['condition'], categories=['control', 'treatment'])

dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~batch + condition",  # Batch first
    quiet=True
)
dds.deseq2()

# Extract condition effect (adjusted for batch)
stat_res = DeseqStats(dds, contrast=['condition', 'treatment', 'control'], quiet=True)
stat_res.run_wald_test()
stat_res.summary()
results = stat_res.results_df

Continuous Covariates

# Design with continuous variable (e.g., age, time)
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~age + condition",
    continuous_factors=['age'],  # Specify continuous variables
    quiet=True
)
dds.deseq2()

# Extract condition effect (adjusted for age)
stat_res = DeseqStats(dds, contrast=['condition', 'treatment', 'control'], quiet=True)
stat_res.run_wald_test()
results = stat_res.results_df

LFC Shrinkage

# After running Wald test, apply shrinkage
stat_res = DeseqStats(dds, contrast=['condition', 'treatment', 'control'], quiet=True)
stat_res.run_wald_test()

# Determine coefficient name for shrinkage
# Format: factor[T.level] where level is the numerator
coeff = 'condition[T.treatment]'

# Verify coefficient exists
if coeff in dds.varm['LFC'].columns:
    stat_res.lfc_shrink(coeff=coeff)
else:
    print(f"WARNING: Coefficient '{coeff}' not found.")
    print(f"Available: {list(dds.varm['LFC'].columns)}")

results = stat_res.results_df

Set Reference Level

CRITICAL: In PyDESeq2 v0.5.4+, use pd.Categorical with ordered categories. The FIRST category is the reference.

def set_reference_level(metadata, factor_col, ref_value):
    """Set reference level for a factor by reordering Categorical.

    The FIRST category becomes the reference level in PyDESeq2.
    """
    current_cats = metadata[factor_col].unique().tolist()
    if ref_value not in current_cats:
        raise ValueError(f"Reference '{ref_value}' not in categories: {current_cats}")

    # Put reference first
    ordered_cats = [ref_value] + [c for c in current_cats if c != ref_value]
    metadata[factor_col] = pd.Categorical(
        metadata[factor_col],
        categories=ordered_cats
    )
    return metadata

# Usage
metadata = set_reference_level(metadata, 'condition', 'wildtype')

Multiple Contrasts

# Run DESeq2 once
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", quiet=True)
dds.deseq2()

# Extract multiple contrasts
contrasts = [
    ['condition', 'A', 'control'],
    ['condition', 'B', 'control'],
    ['condition', 'C', 'control']
]

results_dict = {}
for contrast in contrasts:
    stat_res = DeseqStats(dds, contrast=contrast, quiet=True)
    stat_res.run_wald_test()
    stat_res.summary()
    contrast_name = f"{contrast[1]}_vs_{contrast[2]}"
    results_dict[contrast_name] = stat_res.results_df

Alternative Multiple Testing Correction

from statsmodels.stats.multitest import multipletests
import numpy as np

# Run DESeq2 with default BH correction
stat_res = DeseqStats(dds, contrast=['condition', 'treatment', 'control'], quiet=True)
stat_res.run_wald_test()
results = stat_res.results_df

# Apply alternative correction
pvalues = results['pvalue'].values
mask = ~np.isnan(pvalues)

# Bonferroni
_, results.loc[mask, 'padj_bonf'], _, _ = multipletests(pvalues[mask], method='bonferroni')

# Benjamini-Yekutieli
_, results.loc[mask, 'padj_by'], _, _ = multipletests(pvalues[mask], method='fdr_by')

# Holm
_, results.loc[mask, 'padj_holm'], _, _ = multipletests(pvalues[mask], method='holm')

Dispersion Fitting Options

# Parametric fit (default, recommended for large samples)
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", fit_type='parametric')
dds.deseq2()

# Mean fit (for small samples or when parametric fails)
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", fit_type='mean')
dds.deseq2()

# If dispersion trend doesn't converge, use mean
try:
    dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", fit_type='parametric')
    dds.deseq2()
except Exception as e:
    print(f"Parametric fit failed: {e}")
    print("Retrying with fit_type='mean'")
    dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition", fit_type='mean')
    dds.deseq2()

Access Normalized Counts

# After running dds.deseq2(), normalized counts are available
normalized_counts = dds.obsm['normed_counts']  # DataFrame, same shape as counts

# Size factors
size_factors = dds.obs['size_factors']

# Manual normalization
manual_norm = counts.div(size_factors, axis=0)

Cook's Distance Filtering

# Set minimum replicates for Cook's filtering
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~condition",
    min_replicates=7,  # Default is 7
    quiet=True
)
dds.deseq2()

# Genes with Cook's outliers will have padj = NA
# These are automatically excluded in downstream filtering

Pre-filtering Low-Count Genes

# Filter genes with low mean counts before DESeq2 (optional, improves speed)
min_count = 10
keep_genes = counts.sum(axis=0) >= min_count
counts_filtered = counts.loc[:, keep_genes]

# Then run DESeq2 on filtered data
dds = DeseqDataSet(counts=counts_filtered, metadata=metadata, design="~condition", quiet=True)
dds.deseq2()

Complete Example: Multi-Factor with Shrinkage

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

# Load and validate data (from data_loading.md)
counts = pd.read_csv("counts.csv", index_col=0)
metadata = pd.read_csv("metadata.csv", index_col=0)

# Orient and validate
if counts.shape[0] > counts.shape[1] * 5:
    counts = counts.T
common = sorted(set(counts.index) & set(metadata.index))
counts = counts.loc[common].astype(int)
metadata = metadata.loc[common]

# Set reference levels
metadata['strain'] = pd.Categorical(metadata['strain'], categories=['WT', 'mutant'])
metadata['replicate'] = pd.Categorical(metadata['replicate'])

# Create and fit DESeq2 model
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~replicate + strain",
    quiet=True
)
dds.deseq2()

# Extract results with shrinkage
stat_res = DeseqStats(dds, contrast=['strain', 'mutant', 'WT'], quiet=True)
stat_res.run_wald_test()
stat_res.summary()
stat_res.lfc_shrink(coeff='strain[T.mutant]')

results = stat_res.results_df

# Filter DEGs
sig_genes = results[
    (results['padj'] < 0.05) &
    (results['log2FoldChange'].abs() > 0.5)
]

print(f"Significant DEGs: {len(sig_genes)}")
print(f"Upregulated: {(sig_genes['log2FoldChange'] > 0).sum()}")
print(f"Downregulated: {(sig_genes['log2FoldChange'] < 0).sum()}")

Question Parsing Guide

Extract parameters from user questions before analysis.

Parameter Extraction Table

File Discovery

Look in the working directory for:

import os
import glob

# List all data files
data_dir = "."  # or specified path
all_files = glob.glob(os.path.join(data_dir, "**/*"), recursive=True)
data_files = [f for f in all_files if f.endswith(('.csv', '.tsv', '.txt', '.h5ad', '.rds', '.h5'))]

# Common patterns
count_files = [f for f in data_files if any(x in f.lower() for x in ['count', 'expression', 'matrix'])]
meta_files = [f for f in data_files if any(x in f.lower() for x in ['metadata', 'coldata', 'sample', 'design', 'pheno'])]

Decision Tree

Q: Is there a count matrix file?
  YES -> Load and proceed to validation
  NO  -> Q: Is there an h5ad/AnnData file?
           YES -> Load with anndata, extract counts
           NO  -> Q: Is there processed DE results already?
                    YES -> Skip to filtering/enrichment
                    NO  -> ERROR: No suitable input data found

Multi-Factor Design Examples

Question: "Account for Replicate, Strain, and Media"

design = "~Replicate + Strain + Media"

Question: "Control for batch effects"

design = "~batch + condition"

Question: "Include interaction between strain and treatment"

design = "~strain + treatment + strain:treatment"

Contrast Parsing

Question: "mutant vs wildtype"

contrast = ['condition', 'mutant', 'wildtype']

Question: "Compare treatment to control"

contrast = ['condition', 'treatment', 'control']

Question: "Effect of strain B relative to strain A"

contrast = ['strain', 'B', 'A']

Subset Identification

Question: "Analyze control mice only"

subset = metadata[metadata['treatment'] == 'control']

Question: "Excluding the third replicates"

subset = metadata[metadata['replicate'] != 3]

Question: "CD4 and CD8 cells only"

subset = metadata[metadata['cell_type'].isin(['CD4', 'CD8'])]

Enrichment Library Selection

Organism Detection

def detect_organism(metadata, gene_names):
    """Detect organism from gene naming patterns."""
    # Check gene names
    if any(g.startswith('ENSG') for g in gene_names[:100]):
        return 'homo_sapiens'
    elif any(g.startswith('ENSMUSG') for g in gene_names[:100]):
        return 'mus_musculus'

    # Check metadata columns
    if 'organism' in metadata.columns:
        org = metadata['organism'].iloc[0].lower()
        if 'human' in org or 'sapiens' in org:
            return 'homo_sapiens'
        elif 'mouse' in org or 'musculus' in org:
            return 'mus_musculus'

    return 'unknown'

Result Filtering and Extraction

Advanced filtering patterns for DESeq2 results.

Basic Filtering

def filter_degs(results_df, padj_threshold=0.05, lfc_threshold=0,
                basemean_threshold=0, direction='both'):
    """Filter differentially expressed genes.

    Args:
        results_df: DESeq2 results DataFrame
        padj_threshold: Adjusted p-value cutoff
        lfc_threshold: Absolute log2 fold change cutoff (0 = no filter)
        basemean_threshold: Minimum mean expression
        direction: 'both', 'up', or 'down'

    Returns:
        Filtered DataFrame
    """
    df = results_df.dropna(subset=['padj'])  # Remove NaN padj

    # Apply filters
    mask = df['padj'] < padj_threshold

    if lfc_threshold > 0:
        mask = mask & (df['log2FoldChange'].abs() > lfc_threshold)

    if basemean_threshold > 0:
        mask = mask & (df['baseMean'] >= basemean_threshold)

    if direction == 'up':
        mask = mask & (df['log2FoldChange'] > 0)
    elif direction == 'down':
        mask = mask & (df['log2FoldChange'] < 0)

    return df[mask]

Extract Specific Gene Results

def get_gene_result(results_df, gene_name, column='log2FoldChange'):
    """Get a specific value for a specific gene.

    Handles case-insensitive matching and common naming issues.
    """
    # Try exact match first
    if gene_name in results_df.index:
        return results_df.loc[gene_name, column]

    # Case-insensitive match
    idx_lower = {g.lower(): g for g in results_df.index}
    if gene_name.lower() in idx_lower:
        actual_name = idx_lower[gene_name.lower()]
        return results_df.loc[actual_name, column]

    # Partial match (for gene IDs like PA14_35160)
    matches = [g for g in results_df.index if gene_name.lower() in g.lower()]
    if len(matches) == 1:
        return results_df.loc[matches[0], column]
    elif len(matches) > 1:
        return {m: results_df.loc[m, column] for m in matches}

    return None  # Gene not found

Top N Genes

# Top upregulated by log2FC
top_up = results.sort_values('log2FoldChange', ascending=False).head(10)

# Top downregulated by log2FC
top_down = results.sort_values('log2FoldChange').head(10)

# Top by adjusted p-value (most significant)
top_sig = results.sort_values('padj').head(10)

# Top by baseMean (highest expression)
top_expr = results.sort_values('baseMean', ascending=False).head(10)

Quantile Filtering

# Top 10% by log2FC magnitude
lfc_threshold = results['log2FoldChange'].abs().quantile(0.9)
top_10_percent = results[results['log2FoldChange'].abs() >= lfc_threshold]

# Top quartile by baseMean
basemean_threshold = results['baseMean'].quantile(0.75)
high_expr = results[results['baseMean'] >= basemean_threshold]

Combined Filtering

# Significant AND highly expressed
sig_high = results[
    (results['padj'] < 0.05) &
    (results['log2FoldChange'].abs() > 1) &
    (results['baseMean'] > results['baseMean'].median())
]

# Significant OR high fold change
sig_or_high_fc = results[
    (results['padj'] < 0.05) |
    (results['log2FoldChange'].abs() > 2)
]

Set Operations on DEG Lists

def compare_deg_sets(deg_sets, operation='unique'):
    """Compare DEG sets across conditions.

    Args:
        deg_sets: Dict of {condition_name: set_of_gene_names}
        operation: 'unique' (per condition), 'shared' (intersection),
                   'union', 'venn' (all combinations)

    Returns:
        Dict with results
    """
    results = {}
    condition_names = list(deg_sets.keys())
    all_genes = set()
    for genes in deg_sets.values():
        all_genes |= genes

    if operation == 'unique':
        # Genes unique to each condition (not in any other)
        for cond in condition_names:
            others = set()
            for other_cond in condition_names:
                if other_cond != cond:
                    others |= deg_sets[other_cond]
            results[cond] = deg_sets[cond] - others

    elif operation == 'shared':
        # Intersection of all
        shared = deg_sets[condition_names[0]]
        for cond in condition_names[1:]:
            shared = shared & deg_sets[cond]
        results['shared'] = shared

    elif operation == 'union':
        results['union'] = all_genes

    elif operation == 'venn':
        # All combinations
        from itertools import combinations
        for r in range(1, len(condition_names) + 1):
            for combo in combinations(condition_names, r):
                label = ' & '.join(combo)
                intersection = deg_sets[combo[0]]
                for cond in combo[1:]:
                    intersection = intersection & deg_sets[cond]
                # Remove genes that appear in conditions not in this combo
                others = [c for c in condition_names if c not in combo]
                for other in others:
                    intersection = intersection - deg_sets[other]
                results[label] = intersection

    return results

Direction-Concordant DEGs

def find_concordant_degs(results_A, results_B, padj_threshold=0.05):
    """Find genes DE in same direction in both comparisons."""
    # Both significant
    sig_A = results_A[results_A['padj'] < padj_threshold].index
    sig_B = results_B[results_B['padj'] < padj_threshold].index
    common = set(sig_A) & set(sig_B)

    # Same direction
    concordant = []
    for gene in common:
        lfc_A = results_A.loc[gene, 'log2FoldChange']
        lfc_B = results_B.loc[gene, 'log2FoldChange']
        if (lfc_A > 0 and lfc_B > 0) or (lfc_A < 0 and lfc_B < 0):
            concordant.append(gene)

    return concordant

Filter by Gene List

# Filter to specific genes of interest
genes_of_interest = ['TP53', 'BRCA1', 'EGFR', 'MYC']
results_filtered = results.loc[results.index.intersection(genes_of_interest)]

# Filter to genes in a pathway
pathway_genes = ['GENE1', 'GENE2', 'GENE3']  # From pathway database
pathway_results = results.loc[results.index.intersection(pathway_genes)]

Rank Genes

# Rank by combined metric: -log10(padj) * sign(log2FC)
results['rank_metric'] = -np.log10(results['padj']) * np.sign(results['log2FoldChange'])
results_ranked = results.sort_values('rank_metric', ascending=False)

# Alternative: Use stat column (Wald statistic)
results_ranked = results.sort_values('stat', ascending=False)

Export Filtered Results

# Export significant DEGs
sig_genes = filter_degs(results, padj_threshold=0.05, lfc_threshold=0.5)
sig_genes.to_csv('significant_degs.csv')

# Export top genes
top_genes = results.sort_values('padj').head(100)
top_genes.to_csv('top_100_degs.csv')

# Export gene list only
gene_list = sig_genes.index.tolist()
with open('gene_list.txt', 'w') as f:
    f.write('\n'.join(gene_list))

Summary Statistics

def summarize_results(results, padj_threshold=0.05, lfc_threshold=0):
    """Generate summary statistics for DESeq2 results."""
    sig = results[
        (results['padj'] < padj_threshold) &
        (results['log2FoldChange'].abs() > lfc_threshold)
    ]

    summary = {
        'total_genes': len(results),
        'genes_with_padj': len(results.dropna(subset=['padj'])),
        'significant_genes': len(sig),
        'upregulated': len(sig[sig['log2FoldChange'] > 0]),
        'downregulated': len(sig[sig['log2FoldChange'] < 0]),
        'mean_lfc_up': sig[sig['log2FoldChange'] > 0]['log2FoldChange'].mean(),
        'mean_lfc_down': sig[sig['log2FoldChange'] < 0]['log2FoldChange'].mean(),
        'max_lfc': sig['log2FoldChange'].max(),
        'min_lfc': sig['log2FoldChange'].min(),
        'median_basemean': sig['baseMean'].median()
    }

    return summary

Complex Filtering Example

# Multi-criteria filtering for publication-ready DEG list
def get_publication_degs(results):
    """Filter for high-confidence, biologically meaningful DEGs."""
    filtered = results[
        (results['padj'] < 0.01) &  # Stringent significance
        (results['log2FoldChange'].abs() > 1) &  # 2-fold change
        (results['baseMean'] > 50) &  # Adequate expression
        (results['lfcSE'] < 0.5)  # Reasonable SE
    ].copy()

    # Remove outliers (very high LFC might be technical artifacts)
    lfc_upper = filtered['log2FoldChange'].quantile(0.95)
    lfc_lower = filtered['log2FoldChange'].quantile(0.05)
    filtered = filtered[
        (filtered['log2FoldChange'] <= lfc_upper) &
        (filtered['log2FoldChange'] >= lfc_lower)
    ]

    # Sort by combined metric
    filtered['score'] = -np.log10(filtered['padj']) * filtered['log2FoldChange'].abs()
    filtered = filtered.sort_values('score', ascending=False)

    return filtered

Example: Multi-Condition Analysis

# Compare 3 conditions to control
conditions = ['A', 'B', 'C']
all_results = {}

for cond in conditions:
    stat_res = DeseqStats(dds, contrast=['condition', cond, 'control'], quiet=True)
    stat_res.run_wald_test()
    all_results[cond] = stat_res.results_df

# Get DEG sets
deg_sets = {
    cond: set(all_results[cond][all_results[cond]['padj'] < 0.05].index)
    for cond in conditions
}

# Find unique and shared DEGs
comparison = compare_deg_sets(deg_sets, operation='unique')
print(f"Unique to A: {len(comparison['A'])}")
print(f"Unique to B: {len(comparison['B'])}")
print(f"Unique to C: {len(comparison['C'])}")

shared = compare_deg_sets(deg_sets, operation='shared')
print(f"Shared across all: {len(shared['shared'])}")

Troubleshooting Guide

Common issues and solutions for RNA-seq analysis with PyDESeq2.

Data Loading Issues

"No matching samples between counts and metadata"

Cause: Sample names don't match between count matrix and metadata.

Solutions:

# Check sample names
print("Count samples:", list(counts.index)[:5])
print("Metadata samples:", list(metadata.index)[:5])

# Try case-insensitive matching
counts.index = counts.index.str.lower()
metadata.index = metadata.index.str.lower()

# Remove whitespace
counts.index = counts.index.str.strip()
metadata.index = metadata.index.str.strip()

# Try transpose
if set(counts.columns) & set(metadata.index):
    counts = counts.T

"Non-integer counts detected"

Cause: Data is normalized (FPKM, TPM) or has floating point values.

Solutions:

# Option 1: Round to integers (acceptable if close to integers)
counts = counts.round().astype(int)

# Option 2: Use t-test instead of DESeq2 for normalized data
from scipy import stats
stat, pval = stats.ttest_ind(group1, group2)

# Check if data is truly raw counts
print(counts.head())
print("Min value:", counts.min().min())
print("Has decimals:", (counts % 1 != 0).any().any())

Matrix orientation wrong

Cause: Genes are rows instead of columns.

Solution:

# Check shape
print(f"Counts shape: {counts.shape}")  # Should be (n_samples, n_genes)
print(f"Samples >> genes? {counts.shape[0] > counts.shape[1] * 5}")

# Transpose if needed
if counts.shape[0] > counts.shape[1] * 5:
    counts = counts.T

DESeq2 Execution Issues

"Dispersion trend did not converge"

Cause: Small sample size or low variation.

Solution:

# Use mean fit instead of parametric
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~condition",
    fit_type='mean',  # Instead of 'parametric'
    quiet=True
)
dds.deseq2()

"Contrast not found"

Cause: Wrong factor/level names in contrast.

Solution:

# Check available levels
print("Available factors:", metadata.columns)
print("Condition levels:", metadata['condition'].unique())

# Verify exact names (case-sensitive)
contrast = ['condition', 'treatment', 'control']  # Must match exactly

"All genes filtered out"

Cause: Too strict pre-filtering or all counts are zero.

Solution:

# Check data quality
print("Zero genes:", (counts.sum(axis=0) == 0).sum())
print("Low count genes:", (counts.sum(axis=0) < 10).sum())

# Remove only zero genes
nonzero = counts.sum(axis=0) > 0
counts = counts.loc[:, nonzero]

# Don't pre-filter too aggressively
# DESeq2 handles low counts internally

"Single replicate per condition"

Cause: Only 1 sample per condition - cannot estimate dispersion.

Solution:

# Cannot run DESeq2 with single replicates
# Use fold-change only
mean_A = counts[samples_A].mean(axis=1)
mean_B = counts[samples_B].mean(axis=1)
log2fc = np.log2((mean_B + 1) / (mean_A + 1))

# Or pool replicates from similar conditions

Reference Level Issues

Wrong reference level

Cause: Reference level not set correctly.

Solution:

# In PyDESeq2, FIRST category is reference
metadata['condition'] = pd.Categorical(
    metadata['condition'],
    categories=['control', 'treatment']  # Control first = reference
)

# Verify
print("Categories:", metadata['condition'].cat.categories)
print("First category (reference):", metadata['condition'].cat.categories[0])

LFC Shrinkage Issues

"Coefficient not found for shrinkage"

Cause: Wrong coefficient name format.

Solution:

# Check available coefficients
print("Available coefficients:", list(dds.varm['LFC'].columns))

# Standard format: factor[T.level]
coeff = 'condition[T.treatment]'

# Verify before shrinking
if coeff in dds.varm['LFC'].columns:
    stat_res.lfc_shrink(coeff=coeff)
else:
    print(f"ERROR: '{coeff}' not found")
    print("Skipping shrinkage")

Result Extraction Issues

"NaN in padj column"

Cause: Independent filtering removed genes with insufficient evidence.

Solution:

# This is EXPECTED behavior
# Remove NaN before counting DEGs
sig_genes = results.dropna(subset=['padj'])
sig_genes = sig_genes[sig_genes['padj'] < 0.05]

# Don't include NaN genes in counts
answer = len(sig_genes)  # Correct
# NOT: len(results[results['padj'] < 0.05])  # Wrong, includes NaN

"Gene not found in results"

Cause: Gene name case mismatch or gene filtered out.

Solution:

# Case-insensitive search
def find_gene(results_df, gene_name):
    # Exact match
    if gene_name in results_df.index:
        return gene_name

    # Case-insensitive
    idx_lower = {g.lower(): g for g in results_df.index}
    if gene_name.lower() in idx_lower:
        return idx_lower[gene_name.lower()]

    # Partial match
    matches = [g for g in results_df.index if gene_name.lower() in g.lower()]
    if len(matches) == 1:
        return matches[0]
    elif len(matches) > 1:
        print(f"Multiple matches: {matches}")

    return None

gene = find_gene(results, "TP53")
if gene:
    lfc = results.loc[gene, 'log2FoldChange']

Enrichment Analysis Issues

"No enrichment results"

Cause: Gene list too small or wrong organism/library.

Solution:

# Check gene list
print(f"Gene list size: {len(gene_list)}")
print(f"Sample genes: {gene_list[:5]}")

# Try different libraries
libraries = [
    'GO_Biological_Process_2023',
    'GO_Biological_Process_2021',
    'KEGG_2021_Human',
    'KEGG_2019_Mouse'
]

for lib in libraries:
    try:
        enr = gp.enrich(gene_list=gene_list, gene_sets=lib, outdir=None, no_plot=True)
        if len(enr.results) > 0:
            print(f"✓ {lib}: {len(enr.results)} results")
            break
    except Exception as e:
        print(f"✗ {lib}: {e}")

"gseapy organism parameter deprecated"

Cause: Using old gseapy syntax with organism parameter.

Solution:

# OLD (doesn't work in gseapy >= 1.1)
# enr = gp.enrich(gene_list=gene_list, organism='human', gene_sets='GO_Biological_Process')

# NEW (correct)
enr = gp.enrich(
    gene_list=gene_list,
    gene_sets='GO_Biological_Process_2023',  # Organism in library name
    outdir=None
)

Memory Issues

"Memory error with large gene set"

Cause: Too many genes (60K+ genes with many samples).

Solution:

# Pre-filter low-count genes
min_count = 10
keep = counts.sum(axis=0) >= min_count
counts_filtered = counts.loc[:, keep]
print(f"Kept {keep.sum()} / {len(keep)} genes")

# Use sparse matrices (if applicable)
from scipy.sparse import csr_matrix
counts_sparse = csr_matrix(counts.values)

Multi-Factor Design Issues

"Singular design matrix"

Cause: Confounded factors (e.g., batch perfectly correlated with condition).

Solution:

# Check confounding
print(pd.crosstab(metadata['batch'], metadata['condition']))

# If confounded, cannot separate effects
# Remove confounded factor from design
design = "~condition"  # Remove batch if confounded

Performance Issues

"DESeq2 takes too long"

Cause: Large dataset or many factors.

Solutions:

# Use parallel processing (if available)
import multiprocessing
n_cpus = multiprocessing.cpu_count()
# Note: PyDESeq2 doesn't directly support n_cpus, but numpy may use multiple cores

# Pre-filter more aggressively
min_count = 10
min_samples = 3
keep = (counts >= min_count).sum(axis=0) >= min_samples
counts_filtered = counts.loc[:, keep]

# Use mean fit instead of parametric (faster)
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition",
                   fit_type='mean', quiet=True)

Validation Issues

Results don't match expected answer

Checks:

# 1. Check reference level
print("Reference:", metadata['condition'].cat.categories[0])

# 2. Check contrast direction
print("Contrast:", contrast)  # [factor, numerator, denominator]

# 3. Check filtering thresholds
print(f"padj < {padj_threshold}")
print(f"|log2FC| > {lfc_threshold}")

# 4. Check if shrinkage was applied
print("Shrinkage applied:", 'lfcSE' in results.columns)

# 5. Check for NaN handling
print("Genes with NaN padj:", results['padj'].isna().sum())

Debugging Checklist

Before reporting issues, verify:

• [ ] Count matrix has samples as rows, genes as columns
• [ ] Counts are non-negative integers
• [ ] Metadata index matches count matrix index exactly
• [ ] Design formula references valid column names in metadata
• [ ] Reference level is set correctly (first category in Categorical)
• [ ] Contrast factor and levels exist in metadata
• [ ] LFC shrinkage coefficient name matches pydeseq2 format
• [ ] Filtering thresholds match question exactly
• [ ] NaN values in padj are excluded from DEG counts

Getting Help

# Print diagnostic information
print("=" * 50)
print("DIAGNOSTIC INFORMATION")
print("=" * 50)
print(f"Counts shape: {counts.shape}")
print(f"Metadata shape: {metadata.shape}")
print(f"Design: {design}")
print(f"Reference level: {metadata['condition'].cat.categories[0]}")
print(f"Contrast: {contrast}")
print(f"PyDESeq2 version: {pydeseq2.__version__}")
print(f"\nSample counts:\n{counts.iloc[:3, :3]}")
print(f"\nMetadata:\n{metadata.head(3)}")
print(f"\nResults:\n{results.head(3)}")
print("=" * 50)

DESeq2 Question Patterns

All 10 common question patterns with worked examples.

Pattern 1: Basic DEG Count

Question: "How many genes show significant DE (padj < 0.05, |log2FC| > 0.5)?"

Code:

degs = results[
    (results['padj'] < 0.05) &
    (results['log2FoldChange'].abs() > 0.5)
]
answer = len(degs)

Variations:

• "How many differentially expressed genes?"
• "Number of significant genes"
• "Count of DEGs with adjusted p-value < 0.01"

---

Pattern 2: Specific Gene Value

Question: "What is the log2FC of gene X?"

Code:

answer = round(results.loc['GENE_X', 'log2FoldChange'], 2)

Handling missing genes:

def get_gene_result(results_df, gene_name, column='log2FoldChange'):
    """Get value with case-insensitive matching."""
    # Try exact match first
    if gene_name in results_df.index:
        return results_df.loc[gene_name, column]

    # Case-insensitive match
    idx_lower = {g.lower(): g for g in results_df.index}
    if gene_name.lower() in idx_lower:
        actual_name = idx_lower[gene_name.lower()]
        return results_df.loc[actual_name, column]

    # Partial match
    matches = [g for g in results_df.index if gene_name.lower() in g.lower()]
    if len(matches) == 1:
        return results_df.loc[matches[0], column]

    return None  # Gene not found

answer = round(get_gene_result(results, "TP53", "log2FoldChange"), 2)

Variations:

• "What is the padj for gene Y?"
• "What is the baseMean of gene Z?"
• "What is the p-value for gene ABC?"

---

Pattern 3: Direction-Specific DEGs

Question: "How many genes are upregulated?"

Code:

up_degs = results[
    (results['padj'] < 0.05) &
    (results['log2FoldChange'] > 0)
]
answer = len(up_degs)

Downregulated:

down_degs = results[
    (results['padj'] < 0.05) &
    (results['log2FoldChange'] < 0)
]
answer = len(down_degs)

With LFC threshold:

up_degs = results[
    (results['padj'] < 0.05) &
    (results['log2FoldChange'] > 0.5)  # Threshold
]
answer = len(up_degs)

---

Pattern 4: Multi-Condition Comparison (Set Operations)

Question: "How many genes are uniquely DE in condition A but not B or C?"

Code:

# Extract DEG sets for each condition
degs_A = set(results_A[results_A['padj'] < 0.05].index)
degs_B = set(results_B[results_B['padj'] < 0.05].index)
degs_C = set(results_C[results_C['padj'] < 0.05].index)

# Unique to A
unique_A = degs_A - degs_B - degs_C
answer = len(unique_A)

Shared across all:

shared = degs_A & degs_B & degs_C
answer = len(shared)

Percentage overlap:

overlap = degs_A & degs_B
percentage = round(len(overlap) / len(degs_A) * 100, 1)

Compare two sets:

def compare_deg_sets(deg_sets, operation='unique'):
    """Compare DEG sets across conditions."""
    results = {}
    condition_names = list(deg_sets.keys())

    if operation == 'unique':
        # Genes unique to each condition
        for cond in condition_names:
            others = set()
            for other_cond in condition_names:
                if other_cond != cond:
                    others |= deg_sets[other_cond]
            results[cond] = deg_sets[cond] - others

    elif operation == 'shared':
        # Intersection of all
        shared = deg_sets[condition_names[0]]
        for cond in condition_names[1:]:
            shared = shared & deg_sets[cond]
        results['shared'] = shared

    return results

# Usage
deg_sets = {
    'A': set(results_A[results_A['padj'] < 0.05].index),
    'B': set(results_B[results_B['padj'] < 0.05].index),
    'C': set(results_C[results_C['padj'] < 0.05].index)
}
unique = compare_deg_sets(deg_sets, operation='unique')
answer = len(unique['A'])

---

Pattern 5: Dispersion Count

Question: "How many genes have dispersion below 1e-5 prior to fitting?"

Code:

genewise = dds.var['genewise_dispersions']
answer = (genewise < 1e-5).sum()

After shrinkage:

map_disp = dds.var['MAP_dispersions']
answer = (map_disp < 1e-5).sum()

Dispersion diagnostics:

def dispersion_diagnostics(dds, threshold=1e-5):
    """Analyze dispersion estimates."""
    diag = {}

    if 'genewise_dispersions' in dds.var.columns:
        gwd = dds.var['genewise_dispersions']
        diag['genewise_below_threshold'] = (gwd < threshold).sum()
        diag['genewise_min'] = gwd.min()
        diag['genewise_max'] = gwd.max()
        diag['genewise_median'] = gwd.median()

    if 'MAP_dispersions' in dds.var.columns:
        mapd = dds.var['MAP_dispersions']
        diag['MAP_below_threshold'] = (mapd < threshold).sum()

    return diag

# Usage
diag = dispersion_diagnostics(dds, threshold=1e-5)
answer = diag['genewise_below_threshold']

---

Pattern 6: DEGs + Enrichment

Question: "What is the adjusted p-value for pathway X in enrichment of DEGs?"

Code:

import gseapy as gp

# Get DEGs
degs = results[
    (results['padj'] < 0.05) &
    (results['log2FoldChange'].abs() > 0.5)
]
gene_list = degs.index.tolist()

# Run enrichment
enr = gp.enrich(
    gene_list=gene_list,
    gene_sets='KEGG_2021_Human',
    outdir=None,
    no_plot=True,
    verbose=False
)

# Extract answer
pathway = enr.results[enr.results['Term'].str.contains('ABC transporters')]
answer = round(pathway.iloc[0]['Adjusted P-value'], 4)

Count significant pathways:

answer = len(enr.results[enr.results['Adjusted P-value'] < 0.05])

Gene count in pathway:

pathway = enr.results[enr.results['Term'].str.contains('ribosome')]
overlap = pathway.iloc[0]['Overlap']  # e.g., "25/150"
answer = int(overlap.split('/')[0])  # 25

---

Pattern 7: Percentage Calculation

Question: "What percentage of DE genes in A are also DE in B?"

Code:

degs_A = set(results_A[results_A['padj'] < 0.05].index)
degs_B = set(results_B[results_B['padj'] < 0.05].index)
overlap = degs_A & degs_B
percentage = round(len(overlap) / len(degs_A) * 100, 1)

Percentage upregulated:

all_degs = results[results['padj'] < 0.05]
up_degs = all_degs[all_degs['log2FoldChange'] > 0]
percentage = round(len(up_degs) / len(all_degs) * 100, 1)

---

Pattern 8: Wilson Confidence Interval

Question: "What is the 95% CI for the proportion of DEGs using Wilson method?"

Code:

from statsmodels.stats.proportion import proportion_confint

n_total = len(results.dropna(subset=['padj']))
n_sig = len(results[(results['padj'] < 0.05) & (results['log2FoldChange'].abs() > 1)])

ci_low, ci_high = proportion_confint(n_sig, n_total, method='wilson')
answer = (round(ci_low, 2), round(ci_high, 2))

Other CI methods:

# Normal approximation
ci_low, ci_high = proportion_confint(n_sig, n_total, method='normal')

# Clopper-Pearson (exact)
ci_low, ci_high = proportion_confint(n_sig, n_total, method='beta')

---

Pattern 9: Enrichment Gene Count in Pathway

Question: "How many DEGs contribute to pathway X enrichment?"

Code:

import gseapy as gp

# Run enrichment
degs = results[(results['padj'] < 0.05)]
gene_list = degs.index.tolist()
enr = gp.enrich(gene_list=gene_list, gene_sets='KEGG_2021_Human')

# Find pathway
pathway_row = enr.results[enr.results['Term'].str.contains('ABC transporters')].iloc[0]

# Extract gene count
overlap_str = pathway_row['Overlap']  # e.g., "11/42"
n_genes = int(overlap_str.split('/')[0])  # 11

# Or get gene list
genes_in_pathway = pathway_row['Genes'].split(';')
answer = len(genes_in_pathway)

---

Pattern 10: Batch Effect Assessment

Question: "How does removing batch-affected samples change DEG count?"

Code:

# Run with all samples
dds_all = DeseqDataSet(counts=counts_all, metadata=metadata_all, design="~condition", quiet=True)
dds_all.deseq2()
stat_res_all = DeseqStats(dds_all, contrast=['condition', 'treatment', 'control'], quiet=True)
stat_res_all.run_wald_test()
results_all = stat_res_all.results_df
n_all = len(results_all[(results_all['padj'] < 0.05) & (results_all['log2FoldChange'].abs() > 1)])

# Run without batch-affected samples
counts_clean = counts_all.loc[~metadata_all['batch_affected']]
metadata_clean = metadata_all[~metadata_all['batch_affected']]
dds_clean = DeseqDataSet(counts=counts_clean, metadata=metadata_clean, design="~condition", quiet=True)
dds_clean.deseq2()
stat_res_clean = DeseqStats(dds_clean, contrast=['condition', 'treatment', 'control'], quiet=True)
stat_res_clean.run_wald_test()
results_clean = stat_res_clean.results_df
n_clean = len(results_clean[(results_clean['padj'] < 0.05) & (results_clean['log2FoldChange'].abs() > 1)])

# Compare
if n_clean > n_all:
    answer = "Increases the number of differentially expressed genes"
else:
    answer = "Decreases the number of differentially expressed genes"

---

Additional Patterns

Multiple Testing Correction Comparison

Question: "How many genes are significant with Bonferroni vs BH correction?"

from statsmodels.stats.multitest import multipletests
import numpy as np

pvalues = results['pvalue'].values
mask = ~np.isnan(pvalues)

# Benjamini-Hochberg (default)
_, padj_bh, _, _ = multipletests(pvalues[mask], method='fdr_bh')
n_bh = (padj_bh < 0.05).sum()

# Bonferroni
_, padj_bonf, _, _ = multipletests(pvalues[mask], method='bonferroni')
n_bonf = (padj_bonf < 0.05).sum()

answer = f"BH: {n_bh}, Bonferroni: {n_bonf}"

Protein-Coding vs Non-Coding

Question: "Compare DE between protein-coding and non-protein-coding genes"

from tooluniverse import ToolUniverse

tu = ToolUniverse()
tu.load_tools()

# Classify genes
def classify_genes(gene_list, tu):
    """Classify genes as protein-coding or non-protein-coding."""
    classifications = {}
    for gene in gene_list:
        try:
            result = tu.tools.MyGene_query_genes(query=gene)
            if isinstance(result, list) and len(result) > 0:
                gene_type = result[0].get('type_of_gene', 'unknown')
                classifications[gene] = gene_type
        except Exception:
            classifications[gene] = 'unknown'
    return classifications

# Get DEGs
degs = results[results['padj'] < 0.05].index.tolist()
gene_types = classify_genes(degs, tu)

# Count by type
from collections import Counter
type_counts = Counter(gene_types.values())
n_protein_coding = type_counts.get('protein-coding', 0)
n_non_coding = sum(v for k, v in type_counts.items() if k != 'protein-coding' and k != 'unknown')

answer = f"Protein-coding: {n_protein_coding}, Non-coding: {n_non_coding}"

miRNA Analysis

Question: "Differential expression of miRNA data"

from scipy import stats

def mirna_de_analysis(expression_df, groups):
    """Simple DE for pre-normalized miRNA data."""
    group_labels = expression_df.index.map(groups)
    unique_groups = group_labels.unique()

    g1 = expression_df[group_labels == unique_groups[0]]
    g2 = expression_df[group_labels == unique_groups[1]]

    results_list = []
    for gene in expression_df.columns:
        stat, pval = stats.ttest_ind(g1[gene].dropna(), g2[gene].dropna())
        lfc = np.log2(g2[gene].mean() / g1[gene].mean()) if g1[gene].mean() > 0 else np.nan
        results_list.append({
            'gene': gene,
            'log2FoldChange': lfc,
            'pvalue': pval,
            'stat': stat
        })

    results_df = pd.DataFrame(results_list).set_index('gene')

    # Multiple testing correction
    from statsmodels.stats.multitest import multipletests
    mask = ~results_df['pvalue'].isna()
    _, results_df.loc[mask, 'padj'], _, _ = multipletests(
        results_df.loc[mask, 'pvalue'], method='fdr_bh'
    )

    return results_df

# Usage
groups = {'sample1': 'control', 'sample2': 'control', 'sample3': 'treated', 'sample4': 'treated'}
results_mirna = mirna_de_analysis(mirna_expression, groups)
n_sig = len(results_mirna[results_mirna['padj'] < 0.05])