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.