Bulk RNAseq differential expression (DeSeq2)

DESeq2 Differential Expression Analysis

Core DESeq2 workflow for RNA-seq differential expression analysis with count data.

When to Use This Skill

Use DESeq2 when you have:

• ✅ Raw integer count data (not normalized TPM/FPKM)
• ✅ Biological replicates (≥2 per condition, ≥4 recommended)
• ✅ Need for log fold change shrinkage (ranking/visualization)
• ✅ Medium to large sample sizes (DESeq2's strength)

Don't use DESeq2 for:

• ❌ Normalized data (TPM/FPKM) → use limma-voom instead
• ❌ Very small samples (n=2-3) → consider edgeR quasi-likelihood

Quick Start (Example Data)

Test this skill with real RNA-seq data in ~2 minutes:

source("scripts/load_example_data.R")
data <- load_pasilla_data()  # Auto-installs pasilla package if needed (~2 min, ~50MB)
counts <- data$counts        # 14,599 genes × 7 samples
coldata <- data$coldata      # Metadata: treated vs untreated

# Run complete workflow
source("scripts/basic_workflow.R")  # Creates dds, res, resLFC objects + prints summary

What you get:

• Dataset: Drosophila pasilla gene RNAi knockdown (Brooks et al. 2011)
• Comparison: 3 treated vs 4 untreated samples
• Expected results: ~1,000 significant genes at padj < 0.1

For your own data: Replace data loading with your count matrix and metadata (see Inputs section).

Installation

Core packages (required):

# Set CRAN mirror first (required for installation)
options(repos = c(CRAN = "https://cloud.r-project.org"))

if (!require('BiocManager', quietly = TRUE))
    install.packages('BiocManager')
BiocManager::install(c('DESeq2', 'apeglm'))

Example data packages (optional - for testing/learning):

BiocManager::install(c('pasilla', 'airway'))  # ~70MB total, ~2-3 min

Visualization packages (required for QC plots):

# For publication-quality plots (required - generates PNG)
install.packages(c('ggplot2', 'ggprism', 'ggrepel'))

# For SVG export (optional - generates both PNG + SVG)
install.packages('svglite')

License: LGPL (>= 3) (commercial use permitted)

Inputs

Required:

• Count matrix: Raw integer counts (genes × samples)
  • Rows = genes (any identifier: Ensembl, symbols, etc.)
  • Columns = samples
  • Values = non-negative integers
• Sample metadata: Data frame with sample information
  • Row names must match count matrix column names
  • Required column: condition (factor with 2+ levels)
  • Optional: batch, covariates for complex designs

Alternative inputs:

• Salmon/Kallisto output (via tximport)
• SummarizedExperiment object
• featureCounts/HTSeq output
• Bioconductor data packages (pasilla, airway)

See references/deseq2-reference.md for loading examples.

Outputs

Primary results:

• deseq2_results.csv - Full differential expression table (baseMean, log2FC, lfcSE, pvalue, padj)
• deseq2_results_shrunk.csv - Shrunken LFC for visualization/ranking
• dds_object.rds - DESeqDataSet for further analysis

Normalized data:

• normalized_counts.csv - Size-factor normalized counts
• vst_transformed.csv / rlog_transformed.csv - Variance-stabilized values

QC plots (PNG always, SVG strongly preferred, 300 DPI):

• dispersion_plot.png / .svg - Dispersion estimates vs mean
• pca_plot.png / .svg - Principal component analysis
• ma_plot.png / .svg - Mean-average plot
• volcano_plot.png / .svg - Volcano plot (log2FC vs -log10 padj)
• ⚠️ SVG requires svglite package: install.packages('svglite') (falls back to PNG-only if unavailable)

Clarification Questions

⚠️ CRITICAL: Always ask question #1 first to check if user has provided input files before proceeding with analysis.

Before starting, gather:

1. Input Files (ASK THIS FIRST):

   • Do you have specific count matrix file(s) to analyze?
     • If uploaded: Is this the count matrix (genes × samples, raw integer counts)?
     • Expected formats: CSV/TSV, RDS (SummarizedExperiment), Salmon/Kallisto output
   • Or use example data for testing?
     • Use source("scripts/load_example_data.R"); data <- load_pasilla_data()
     • Requires installing pasilla package (~2 min, ~50MB)
   • ⚠️ If data is normalized (TPM/FPKM): Use limma-voom skill instead

2. Sample Metadata (if using own data):

   • What is the primary comparison (e.g., treated vs control)?
   • Which group is the reference/control?
   • Any covariates to adjust for (batch, sex, sequencing run)?
   • Validation: Confirm sample IDs match between count matrix columns and metadata rows

3. Experimental Design:

   • Simple: ~ condition | Multi-factor: ~ batch + condition | Paired: ~ individual + condition | Interaction: ~ genotype * treatment
   • See references/decision-guide.md#design-formulas

4. Sample Size Check:

   • n ≥ 4 per group (recommended) | n = 2-3 (consider edgeR) | n < 2 (insufficient)

5. Significance Thresholds:

   • Standard: padj < 0.05, |log2FC| ≥ 1 | Relaxed: padj < 0.1 | Stringent: padj < 0.01, |log2FC| ≥ 2

6. Analysis Goals:

   • Single pairwise comparison or multiple comparisons?
   • Need visualizations (volcano, heatmap)? → Use de-results-to-plots skill after
   • Need gene annotations? → Use de-results-to-gene-lists skill after

Typical Complete Workflow

This skill performs core differential expression analysis with QC plots. For a complete RNA-seq workflow:

1. This skill: Run DESeq2 → get dds, res, normalized counts, QC plots (PCA, MA, volcano, dispersion)
2. de-results-to-gene-lists: Filter significant genes → add annotations → export
3. de-results-to-plots (optional): Advanced visualizations (heatmaps, custom plots)

Quick start: "Run DESeq2 analysis and filter significant genes with annotations"

Why separate skills? Modular design works across DE methods (DESeq2, edgeR, limma). See Suggested Next Steps for details.

Standard Workflow

Note: Run from the OmicsClaw root directory and add the workflow scripts to sys.path:

import sys; import os; sys.path.insert(0, os.path.abspath('knowledge_base/scripts/bulk-rnaseq-counts-to-de-deseq2'))

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE 🚨

This skill uses low-freedom script execution. You must:

• ✅ Source the scripts using the exact commands below
• ✅ Wait for confirmation messages after each step
• ❌ NOT write inline DESeq2 code
• ❌ NOT rewrite plotting code
• ❌ NOT modify commands unless explicitly adapting for user-specific data

WHY USE SCRIPTS: They handle package installation, data validation, sample ID fixes, and error checking automatically. Writing inline code wastes time, introduces errors, and violates the skill design.

Step 1 - Load example data:

source("scripts/load_example_data.R")
data <- load_pasilla_data()
counts <- data$counts
coldata <- data$coldata

Step 2 - Run DESeq2 analysis:

source("scripts/basic_workflow.R")

DO NOT expand this into inline code. DO NOT write the DESeq2 steps manually. Just source the script.

Step 3 - Generate QC plots:

source("scripts/qc_plots.R")
run_all_qc(dds, res, output_dir = "results")

🚨 DO NOT write inline plotting code (ggsave, plotMA, etc.). Just source the script. 🚨

The script handles PNG + SVG export with graceful fallback for SVG dependencies.

Step 4 - Export results:

source("scripts/export_results.R")
export_all(dds, res, resLFC, output_dir = "results")

DO NOT write custom export code. Use export_all() to save all standard outputs including RDS and transformed counts.

✅ VERIFICATION - You should see these messages:

• After Step 1: "✓ Pasilla dataset loaded successfully" with dimensions
• After Step 2: "✓ Basic workflow completed successfully!" with summary statistics
• After Step 3: "✓ All QC plots generated successfully!" with file names
• After Step 4: "=== Export Complete ===" with list of 6-7 files saved

❌ IF YOU DON'T SEE THESE MESSAGES: You wrote inline code instead of using source(). Stop and use the commands above.

⚠️ CRITICAL - DO NOT:

• ❌ Write inline data loading code → STOP: This violates the skill design. Use source("scripts/load_example_data.R") instead. Inline loading causes sample ID mismatches and missing validations.
• ❌ Write inline DESeq2 workflow code → STOP: This violates the skill design. Use source("scripts/basic_workflow.R") instead. Inline workflow wastes time and introduces bugs.
• ❌ Write inline plotting code (ggsave, plotMA, etc.) → STOP: This violates the skill design. Use source("scripts/qc_plots.R") and run_all_qc() instead. If scripts fail, fix the script, don't rewrite inline.
• ❌ Write custom export code → STOP: This violates the skill design. Use source("scripts/export_results.R") and export_all() instead. Custom export code misses RDS objects and transformed counts needed downstream.
• ❌ Try to install svglite → script handles SVG fallback automatically
• ❌ Use absolute paths for scripts → Always use relative paths scripts/file.R
  • ❌ WRONG: source("/mnt/knowhow/workflows/bulk-rnaseq-counts-to-de-deseq2/scripts/load_example_data.R")
  • ❌ WRONG: setwd("/absolute/path/to/skill")
  • ✅ CORRECT: source("scripts/load_example_data.R") (skill should already be working directory)

⚠️ IF SCRIPTS FAIL - Script Failure Hierarchy:

1. Fix and Retry (90%) - Install missing package, re-run script
2. Modify Script (5%) - Edit the script file itself, document changes
3. Use as Reference (4%) - Read script, adapt approach, cite source
4. Write from Scratch (1%) - Only if genuinely impossible, explain why

NEVER skip directly to writing inline code without trying the script first.

📁 Output Directory Paths:

• ✅ Reco

---

Content truncated.