RNA-seq differential expression analysis with DESeq2, edgeR, and limma-voom — DEG lists, fold changes, dispersion estimation, design formulas including covariates, multi-condition contrasts, and Venn-set operations across groups. Routes across DESeq2 (default), edgeR (QL-F / exact test for small rep

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. Output highlights (parseable): For a multi strain Venn run with notebook style outlier exclusion: 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 : 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): 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](references/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: Situation Prefer Why Standard 2 group, modest n, default ask DESeq2 Most widely published reference; shrinkage + independent filtering tuned for small n. This skill's default. Very small replicate counts (n=2 3/group), simple 2 group edgeR (exact test / QL F) Empirical Bayes dispersion moderation is robust at tiny n; QL F controls FDR well. Large n, complex/multi factor designs, many contrasts, or speed matters limma voom Fast, flexible per gene linear model; voom weights handle heteroscedasticity; duplicateCorrelation for repeated measures; extends to interactions. 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. 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). 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. 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. 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: 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: Question phrasing Read which line "significant DEGs", "padj < 0.05", "DEGs at p.adj<0.05" (alone) SIG padj only "DEGs with \ LFC\ X" or "fold change Y" SIG padjlfc with matching lfc thr "DEGs with baseMean N" or "expressed DEGs" SIG strict (need all three thresholds) "shrunk" / "apeglm" / "ashr" in question The shrunk variant of the matching line "before shrinkage" / "unshrunk" / nothing said about shrinkage The unshrunk variant 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 produc

tooluniverse-rnaseq-deseq2