Skill: RNA-seq Differential Expression

Use When

User has a count matrix or transcript quantification and wants to find differentially expressed genes
User wants to perform statistical DE testing with proper normalization
User needs volcano plots, MA plots, PCA, and heatmaps
User wants to compare DESeq2 vs edgeR results

Required:
- Count matrix (TSV with gene IDs as rows, samples as columns) OR Salmon/kallisto quantification directories plus a tx2gene mapping file
- Sample metadata (TSV with columns: sample_id, condition, and optionally batch)
Optional:
- Tool choice: deseq2 or edger (default: deseq2)
- Contrast: condition levels to compare (e.g., treated,control)
- FDR threshold (default: 0.05)
- log2 fold-change threshold (default: 1)
- Output directory (default: ./de_output)

Read count matrix or import transcript-level estimates via tximport.
Read sample metadata and validate that sample IDs match between counts and metadata.
Build design formula (~ condition, or ~ batch + condition if batch column is provided).
If DESeq2: create DESeqDataSetFromMatrix (or DESeqDataSetFromTximport), run DESeq(), extract results().
If edgeR: create DGEList, calcNormFactors, estimateDisp, glmQLFit, glmQLFTest.
Extract results: gene ID, log2 fold-change, p-value, adjusted p-value, baseMean (DESeq2) or logCPM (edgeR).
Generate plots: PCA of samples, volcano plot, MA plot, heatmap of top DE genes.
Write results table sorted by adjusted p-value.
Report: total genes tested, number significant (up-regulated and down-regulated), summary of top hits.

Requires R with DESeq2 or edgeR installed (Bioconductor packages).
Minimum 3 biological replicates per condition recommended for reliable statistics.
Not suitable for single-cell RNA-seq without modifications.
Python alternative: PyDESeq2 for DESeq2-equivalent analysis.
Common failure cases: sample ID mismatch between counts and metadata, too few replicates, incorrect contrast specification.