/bio-de-deseq2-basics

name: bio-de-deseq2-basics
description: Perform differential expression analysis using DESeq2 in R/Bioconductor. Use for analyzing RNA-seq count data, creating DESeqDataSet objects, running the DESeq workflow, and extracting results with log fold change shrinkage. Use when performing DE analysis with DESeq2.
tool_type: r
primary_tool: DESeq2

Version Compatibility

Reference examples tested with: DESeq2 1.42+, Salmon 1.10+, edgeR 4.0+, scanpy 1.10+

Before using code patterns, verify installed versions match. If versions differ:

• R: `packageVersion('<pkg>')` then `?function_name` to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

DESeq2 Basics

Differential expression analysis using DESeq2 for RNA-seq count data.

Required Libraries

library(DESeq2)
library(apeglm)  # For lfcShrink with type='apeglm'

Installation

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

Creating DESeqDataSet

**Goal:** Construct a DESeqDataSet object from various input formats for DE analysis.

**Approach:** Wrap count data and sample metadata into the DESeq2 container, specifying the experimental design formula.

**"Load my RNA-seq counts into DESeq2"** → Create a DESeqDataSet from a count matrix, SummarizedExperiment, or tximport object with sample metadata and a design formula.

From Count Matrix

# counts: matrix with genes as rows, samples as columns
# coldata: data frame with sample metadata (rownames must match colnames of counts)
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData = coldata,
                               design = ~ condition)

From SummarizedExperiment

library(SummarizedExperiment)
dds <- DESeqDataSet(se, design = ~ condition)

From tximport (Salmon/Kallisto)

library(tximport)
txi <- tximport(files, type = 'salmon', tx2gene = tx2gene)
dds <- DESeqDataSetFromTximport(txi, colData = coldata, design = ~ condition)

Standard DESeq2 Workflow

**Goal:** Run the complete DESeq2 pipeline from raw counts to shrunken log fold change estimates.

**Approach:** Create dataset, pre-filter low-count genes, set reference level, run size factor estimation + dispersion estimation + Wald test, then apply LFC shrinkage.

**"Find differentially expressed genes between treated and control"** → Test for significant expression changes between conditions using negative binomial models with empirical Bayes shrinkage.

# Create DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData = coldata,
                               design = ~ condition)

# Pre-filter low count genes (recommended)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# Set reference level for condition
dds$condition <- relevel(dds$condition, ref = 'control')

# Run DESeq2 pipeline (estimateSizeFactors, estimateDispersions, nbinomWaldTest)
dds <- DESeq(dds)

# Get results
res <- results(dds)

# Apply log fold change shrinkage (recommended for visualization/ranking)
resLFC <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')

Design Formulas

**Goal:** Specify the experimental design to model biological and nuisance variables.

**Approach:** Build R formula objects that encode condition, batch, and interaction terms for the GLM.

# Simple two-group comparison
design = ~ condition

# Controlling for batch effects
design = ~ batch + condition

# Interaction model
design = ~ genotype + treatment + genotype:treatment

# Multi-factor without interaction
design = ~ genotype + treatment

Specifying Contrasts

**Goal:** Extract results for specific pairwise or complex comparisons from a fitted DESeq2 model.

**Approach:** Use coefficient names or contrast vectors to define which groups to compare.

# See available coefficients
resultsNames(dds)

# Results by coefficient name
res <- results(dds, name = 'condition_treated_vs_control')

# Results by contrast (compare specific levels)
res <- results(dds, contrast = c('condition', 'treated', 'control'))

# Contrast with list format (for complex designs)
res <- results(dds, contrast = list('conditionB', 'conditionA'))

Log Fold Change Shrinkage

**Goal:** Reduce noisy fold change estimates for low-count genes to improve ranking and visualization.

**Approach:** Apply empirical Bayes shrinkage (apeglm, ashr, or normal) to moderate log fold changes toward zero.

# apeglm method (default, recommended)
resLFC <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')

# ashr method (alternative)
resLFC <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'ashr')

# normal method (original, less recommended)
resLFC <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'normal')

Setting Significance Thresholds

**Goal:** Control the stringency of differential expression calls using adjusted p-value and fold change cutoffs.

**Approach:** Set alpha for multiple testing correction and optionally apply a minimum log fold change threshold.

# Default: padj < 0.1
res <- results(dds)

# Custom alpha threshold
res <- results(dds, alpha = 0.05)

# With log fold change threshold
res <- results(dds, lfcThreshold = 1)  # |log2FC| > 1

Accessing DESeq2 Results

**Goal:** Retrieve, filter, and sort DE results for downstream use.

**Approach:** Extract results as a data frame, subset by significance, and order by p-value or fold change.

# Summary of results
summary(res)

# Get significant genes
sig <- subset(res, padj < 0.05)

# Order by adjusted p-value
resOrdered <- res[order(res$padj),]

# Order by log fold change
resOrdered <- res[order(abs(res$log2FoldChange), decreasing = TRUE),]

# Convert to data frame
res_