Differential Expression with DESeq2

After quantification you have a gene-by-sample matrix of raw counts, and the next question is which genes change between conditions. DESeq2 models raw counts with a negative binomial distribution, estimates a size factor per sample (via the median-of-ratios method) to normalize for differences in sequencing depth and library composition, and shares information across genes to stabilize dispersion estimates. This lesson runs a complete two-group analysis, from a count matrix to a ranked table of differentially expressed genes.

Building the dataset

library(DESeq2)

counts  <- as.matrix(read.csv("counts.csv", row.names = 1))
coldata <- read.csv("samples.csv", row.names = 1)

# DESeq2 needs a factor; the reference level is the comparison baseline
coldata$condition <- relevel(factor(coldata$condition), ref = "control")

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

# drop genes with almost no signal before fitting
keep <- rowSums(counts(dds) >= 10) >= 3
dds  <- dds[keep, ]

The design formula ~ condition tells DESeq2 which column of colData holds the groups to compare. The reference level, set with relevel, becomes the denominator of every fold change, so treated versus control means log2(treated / control). The columns of the count matrix must be in the same order as the rows of colData, because DESeqDataSetFromMatrix pairs them by position and does not reorder by name.

Fitting the model

dds <- DESeq(dds)
res <- results(dds, alpha = 0.05)
summary(res)

out of 11416 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 1108, 9.7%
LFC < 0 (down)     : 987, 8.6%
outliers [1]       : 0, 0%
low counts [2]     : 3567, 31%
(mean count < 6)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

Each row of res describes one gene. baseMean is the average of normalized counts across all samples, log2FoldChange is the effect size on a log2 scale, and lfcSE is its standard error. The stat column is the Wald statistic, pvalue is the raw p-value, and padj is the Benjamini-Hochberg adjusted p-value that controls the false discovery rate. Setting alpha = 0.05 makes independent filtering optimize for that threshold; DESeq2 removes low-count genes with little chance of significance to boost power, and those filtered genes plus any Cook's-distance outliers receive padj of NA.

Contrasts and shrinkage

resultsNames(dds)

# explicit contrast: c(variable, numerator, denominator)
res <- results(dds, contrast = c("condition", "treated", "control"), alpha = 0.05)

# shrink noisy log2 fold changes for ranking and MA plots
# apeglm needs coef (a name from resultsNames), not contrast
library(apeglm)
resLFC <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")
resLFC <- resLFC[order(resLFC$padj), ]
head(resLFC)

log2 fold change (MAP): condition treated vs control
Wald test p-value: condition treated vs control
DataFrame with 6 rows and 5 columns
             baseMean log2FoldChange     lfcSE       pvalue         padj
            <numeric>      <numeric> <numeric>    <numeric>    <numeric>
FBgn0039155   730.596       -3.98423  0.171607 4.75119e-122 8.13421e-118
FBgn0025111  1501.411        2.88690  0.126582 5.13310e-119 4.39383e-115
FBgn0003360  4342.832       -3.17567  0.142754 8.75565e-112 3.74767e-108
FBgn0029167  3706.117       -2.19341  0.096696 8.06713e-108 4.60411e-104
FBgn0035085   638.234       -2.55651  0.137368  1.11172e-79  3.80587e-76
FBgn0029896   277.552       -2.51631  0.148814  5.24967e-66  1.49876e-62