Differential expression in R with DESeq2, start to finish

Every few months someone asks me for “the DESeq2 script”. There isn’t really one script, but there is a sequence of steps I go through every time, and a short list of mistakes I keep seeing. This post writes that sequence down once, so next time I can send a link.

I use the airway dataset from Bioconductor throughout. It is small, it is real, and it has a paired design, which happens to be exactly the kind of detail people forget to model. Everything below runs on a laptop in a few minutes.

Setup

DESeq2 and friends live on Bioconductor, and as of this year the recommended installer is BiocManager (it replaces the old biocLite() script).

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install(c("DESeq2", "airway", "apeglm",
                       "org.Hs.eg.db", "AnnotationDbi"))
install.packages(c("ggplot2", "pheatmap"))

library(DESeq2)
library(airway)
library(ggplot2)
library(pheatmap)
library(AnnotationDbi)
library(org.Hs.eg.db)

The data

The airway experiment (Himes et al. 2014) used four primary human airway smooth muscle cell lines. Each cell line was either treated with dexamethasone, a glucocorticoid, or left untreated. That gives eight samples: four cell lines times two conditions. The reads have already been aligned and summarized to genes, and the package ships the result as a RangedSummarizedExperiment.

data("airway")
airway
colData(airway)[, c("cell", "dex")]
head(assay(airway), 3)
colSums(assay(airway)) / 1e6   # library sizes in millions of reads

The assay is a matrix with one row per Ensembl gene and one column per sample. Each entry is an integer: the number of reads (or read pairs) assigned to that gene in that sample. The colData holds the sample table, and the two columns we care about are cell (which cell line) and dex (trt or untrt).

Why DESeq2 wants raw counts

This is the first thing to get right, so it gets its own heading. DESeq2 models each gene’s counts with a negative binomial distribution. The variance of a count depends on its mean, and the model uses that relationship directly. Differences in sequencing depth are handled inside the model through size factors, which enter as offsets.

If you hand DESeq2 counts that are already normalized, or TPM, or FPKM, two things go wrong. The mean-variance relationship no longer matches what the model assumes, and the depth correction gets applied on top of a depth correction that already happened. DESeqDataSet() will complain if the values aren’t integers, which catches some of these cases, but rounding TPMs to make the error go away doesn’t fix the underlying problem. If you quantified with Salmon or kallisto, import with tximport and build the object with DESeqDataSetFromTximport(), which is designed for exactly that.

Building the DESeqDataSet

Because airway is already a SummarizedExperiment, building the DESeq2 object is one line. If you have a count matrix from featureCounts or HTSeq instead, DESeqDataSetFromMatrix(countData, colData, design) does the same job.

dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds$dex <- relevel(dds$dex, ref = "untrt")

Two decisions are packed into those lines.

The design is paired. Each cell line contributes one treated and one untreated sample, and the cell lines differ from each other at baseline. Putting cell in the design lets the model absorb those baseline differences, so the dex effect is estimated within each cell line. This is the same idea as a paired t-test versus an unpaired one. By convention the variable you care about goes last in the formula, which also makes it the default for results().

The reference level is set explicitly. R orders factor levels alphabetically, so without relevel() the reference would be trt and every fold change would come out with the sign flipped relative to what you’d expect. I always set the reference by hand, even when the alphabet happens to agree with me, because the next dataset won’t.

Pre-filtering low counts

Genes with almost no reads carry no information about differential expression. Removing them before fitting makes everything faster and the plots cleaner.

smallest_group_size <- 4
keep <- rowSums(counts(dds) >= 10) >= smallest_group_size
dds <- dds[keep, ]
nrow(dds)

This keeps genes with at least 10 reads in at least four samples (four being the size of the smaller condition group). The exact thresholds aren’t sacred. Pre-filtering is mostly a convenience; the statistically meaningful filtering happens later, inside results(), through independent filtering.

Running DESeq()

dds <- DESeq(dds)

This single call runs three steps, and it prints a message for each one: estimating size factors, estimating dispersions (gene-wise, then the fitted trend, then the final values), and fitting the model and running Wald tests. You can run them individually with estimateSizeFactors(), estimateDispersions() and nbinomWaldTest(), but there is rarely a reason to.

Size factors

Size factors account for differences in sequencing depth and RNA composition between samples. DESeq2 uses the median-of-ratios method. For every gene it computes the geometric mean across all samples, which acts as a pseudo-reference sample. For each sample it then takes the ratio of each gene’s count to that reference, and the median of those ratios is the size factor. Using the median means a few hundred strongly changed genes don’t drag the estimate around.

sizeFactors(dds)
colSums(counts(dds)) / mean(colSums(counts(dds)))

The two vectors usually track each other but are not identical, and the difference is the point: size factors respond to composition, not just total reads. Normalized counts are the raw counts divided by the size factors, available with counts(dds, normalized = TRUE). Those are for plotting and for sharing with collaborators, never for feeding back into DESeq2.

Dispersion estimates

In the negative binomial model, the variance of a gene’s count is mu + alpha * mu^2, where mu is the mean and alpha is the dispersion. Dispersion captures biological variability beyond Poisson noise. With only a handful of samples per condition, estimating it separately for each gene is very noisy, so DESeq2 borrows information across genes. It estimates a dispersion per gene, fits a smooth trend of dispersion against mean expression, and then shrinks each gene’s estimate toward the trend. Genes with plenty of evidence for being more variable than the trend are left alone.

plotDispEsts(dds)

In this plot the black points are gene-wise estimates, the red line is the fitted trend, and the blue points are the final estimates used for testing. Genes whose gene-wise estimate is far above the trend are drawn with a blue circle and are not shrunk. A healthy plot shows dispersion decreasing as mean count increases, with the cloud of points scattered around the red curve. If the trend is flat or the dispersions are very high everywhere, I go back and look for a sample swap, a hidden batch, or a design that is missing something.

Look at the samples before the genes

Before I read a single p-value, I want to know whether the samples group the way the experiment says they should. For distances and PCA, raw counts are a bad choice because a few very highly expressed genes dominate. A plain log2(count + 1) overcorrects and inflates noise at low counts. DESeq2 offers two transformations that make the variance roughly constant across the range of means: the variance stabilizing transformation (vst) and the regularized log (rlog).

vsd <- vst(dds, blind = FALSE)
# rld <- rlog(dds, blind = FALSE)   # slower; I use it mainly for small n

blind = TRUE (the default) ignores the design entirely, which is what you want for fully unsupervised QC. blind = FALSE reuses the dispersion trend already fitted with the design. It does not remove the effect of the design variables from the data; it just avoids treating expected large biological differences as noise. Either way, transformed values are for visualization and clustering. They are not the input for the differential expression test.

PCA plot

pcaData <- plotPCA(vsd, intgroup = c("dex", "cell"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(PC1, PC2, color = dex, shape = cell)) +
  geom_point(size = 3) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed() +
  theme_bw()

plotPCA() uses the 500 most variable genes by default (ntop = 500). In airway you should see the treated and untreated samples separate, and the cell lines show up as a second, separate source of variation. That second pattern is the visual argument for the paired design. If you want a PCA with the cell-line effect removed for display only, limma::removeBatchEffect() on assay(vsd) does that, but never give that adjusted matrix to DESeq2.

Sample distance heatmap

sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsd$dex, vsd$cell, sep = "_")
colnames(sampleDistMatrix) <- NULL

pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         color = colorRampPalette(c("navy", "white"))(100))

Darker cells mean more similar samples. I look for two things: samples clustering by condition (or by pair), and any single sample that is far from everything. An isolated sample is worth investigating before you go any further.

Getting results

resultsNames(dds)
res <- results(dds, contrast = c("dex", "trt", "untrt"), alpha = 0.05)
summary(res)
head(res[order(res$padj), ])

resultsNames() lists the model coefficients. Here you’ll see an intercept, three cell-line coefficients, and dex_trt_vs_untrt. The contrast argument says: for the factor dex, compare trt (numerator) to untrt (denominator). Positive log2 fold changes mean higher expression in treated samples.

The result table has six columns: baseMean (the mean of normalized counts across all samples), log2FoldChange, lfcSE (its standard error), stat (the Wald statistic), pvalue and padj. mcols(res)$description tells you exactly which comparison the table describes, and I check it every time.

alpha and independent filtering

The alpha argument is the false discovery rate you intend to use as your cutoff. The default is 0.1. If you are going to call genes at padj < 0.05, set alpha = 0.05, because results() uses that value to tune independent filtering.

Independent filtering removes genes with low mean normalized counts from the multiple testing correction. Those genes have little power to be detected anyway, and the filter statistic (mean count across all samples) doesn’t depend on the condition labels, so removing them doesn’t bias the test under the null. Fewer tests means a smaller multiple testing penalty for the genes that remain. results() tries a range of thresholds and picks the one that maximizes the number of genes with padj below alpha.

metadata(res)$filterThreshold
plot(metadata(res)$filterNumRej, type = "b",
     xlab = "quantile of filter", ylab = "number of rejections")
lines(metadata(res)$lo.fit, col = "red")
abline(v = metadata(res)$filterTheta)

Genes removed by the filter get padj = NA but keep their p-value. A pvalue of NA means something else: either all counts were zero or the gene was flagged as having an outlier (more on that below).

Multiple testing

With around twenty thousand genes tested, a raw p-value cutoff of 0.05 would let through about a thousand false positives even if nothing at all were changing. The padj column applies the Benjamini-Hochberg procedure, which controls the false discovery rate: among the genes you call significant at padj < 0.05, the expected fraction of false discoveries is at most 5%. Use padj for calling genes. The raw pvalue is useful for diagnostics such as a p-value histogram, and for plotting.

hist(res$pvalue[res$baseMean > 1], breaks = 50, col = "grey",
     main = "", xlab = "p-value")

A healthy histogram is flat with a spike near zero. A hump in the middle or a pile-up near one usually points to a problem with the model.

Testing against a fold change threshold

A common request is “genes with padj < 0.05 and at least a two-fold change”. Filtering the table on abs(log2FoldChange) > 1 after the fact is the usual approach, but the p-values in that table still answer a different question (is the fold change different from zero?). DESeq2 can test the question you actually care about:

res_lfc1 <- results(dds, contrast = c("dex", "trt", "untrt"),
                    alpha = 0.05, lfcThreshold = 1)
summary(res_lfc1)

With lfcThreshold = 1 and the default altHypothesis = "greaterAbs", the null hypothesis becomes “the absolute log2 fold change is at most 1”. You will get fewer genes than with post hoc filtering, and those genes come with p-values that mean what you say they mean.

Shrinking log fold changes

The maximum likelihood fold changes in res are noisy for genes with low counts or high dispersion. A gene with counts of 0, 1, 0, 0 in one group and 3, 5, 2, 4 in the other can show an enormous fold change that tells you almost nothing. Shrinkage pulls those estimates toward zero in proportion to their uncertainty, while well-measured genes barely move. The shrunken values are what I use for ranking genes, for MA and volcano plots, and for reporting effect sizes.

resLFC <- lfcShrink(dds, coef = "dex_trt_vs_untrt",
                    type = "apeglm", res = res)

apeglm (Zhu, Ibrahim and Love) uses a heavy-tailed prior, so it shrinks noisy estimates strongly but leaves large, well-supported effects mostly intact. That is a better trade-off than the original normal prior (type = "normal"), which can over-shrink genuinely large changes. Two practical notes. First, apeglm needs a coefficient name from resultsNames(dds) rather than a contrast, so if the comparison you want isn’t a coefficient, relevel the factor and rerun DESeq() (or nbinomWaldTest()). Second, passing res = res carries over the p-values and padj from the table you already made with alpha = 0.05, so shrinkage changes the fold changes and nothing else.

The difference is easiest to see side by side:

par(mfrow = c(1, 2))
plotMA(res, ylim = c(-4, 4), main = "MLE")
plotMA(resLFC, ylim = c(-4, 4), main = "apeglm")
par(mfrow = c(1, 1))

On the left, low-count genes fan out into a wide spray of large fold changes. On the right, that spray collapses toward zero and the genes that stand out are the ones with real support.

MA and volcano plots with ggplot2

plotMA() is fine for a quick look, but for figures I build both plots myself. First, add gene symbols (covered properly in the next section) so the labels mean something.

resLFC$symbol <- mapIds(org.Hs.eg.db, keys = rownames(resLFC),
                        column = "SYMBOL", keytype = "ENSEMBL",
                        multiVals = "first")

plot_df <- as.data.frame(resLFC)
plot_df$sig <- !is.na(plot_df$padj) & plot_df$padj < 0.05

The MA plot shows fold change against mean expression:

ggplot(plot_df, aes(baseMean, log2FoldChange, color = sig)) +
  geom_point(size = 0.6, alpha = 0.5) +
  scale_x_log10() +
  geom_hline(yintercept = 0) +
  scale_color_manual(values = c("grey60", "firebrick"),
                     labels = c("not significant", "padj < 0.05")) +
  labs(x = "mean of normalized counts", y = "log2 fold change (apeglm)",
       color = NULL) +
  theme_bw()

The volcano plot shows fold change against significance. I use the shrunken fold change on the x axis and the raw p-value on the y axis (the ordering is the same as with padj, and the axis is easier to read), and color by the adjusted value.

vol_df <- plot_df[!is.na(plot_df$padj), ]
vol_df$label <- ifelse(rank(vol_df$padj) <= 15, vol_df$symbol, NA)

ggplot(vol_df, aes(log2FoldChange, -log10(pvalue), color = sig)) +
  geom_point(size = 0.8, alpha = 0.6) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_text(aes(label = label), size = 3, vjust = -0.6,
            color = "black", check_overlap = TRUE, na.rm = TRUE) +
  scale_color_manual(values = c("grey60", "firebrick"), guide = "none") +
  labs(x = "log2 fold change (apeglm)", y = "-log10 p-value") +
  theme_bw()

ggrepel::geom_text_repel() does a nicer job with overlapping labels if you have it installed.

Looking at individual genes

Any gene you plan to talk about deserves a look at its actual counts. With a paired design, connecting the two samples from each cell line makes the within-pair change obvious.

top_gene <- rownames(res)[which.min(res$padj)]
d <- plotCounts(dds, gene = top_gene, intgroup = c("dex", "cell"),
                returnData = TRUE)

ggplot(d, aes(dex, count, group = cell, color = cell)) +
  geom_point(size = 3) +
  geom_line() +
  scale_y_log10() +
  labs(title = top_gene, x = NULL, y = "normalized count") +
  theme_bw()

Annotating genes

The row names in airway are Ensembl gene IDs without version suffixes, which is what org.Hs.eg.db expects. mapIds() returns one value per key, so it slots straight into a results table.

res$symbol <- mapIds(org.Hs.eg.db, keys = rownames(res),
                     column = "SYMBOL", keytype = "ENSEMBL",
                     multiVals = "first")
res$entrez <- mapIds(org.Hs.eg.db, keys = rownames(res),
                     column = "ENTREZID", keytype = "ENSEMBL",
                     multiVals = "first")

Some IDs won’t map and will come back as NA. That is normal (many are non-coding or newer annotations), so keep the Ensembl ID as the primary key and treat the symbol as a label. multiVals = "first" quietly picks one symbol when an ID maps to several; if that matters for your analysis, use multiVals = "list" and inspect them. If your IDs carry version suffixes like ENSG00000000003.14, strip them first with sub("\\..*$", "", ids). To see what else you can map to, run columns(org.Hs.eg.db).

Exporting results

I export one table with everything a collaborator needs: the shrunken fold change for effect size, the p-values, the annotation, and the gene ID as a real column instead of row names.

out <- as.data.frame(resLFC)
out$gene_id <- rownames(out)
out$entrez <- res$entrez[match(out$gene_id, rownames(res))]
out <- out[order(out$padj), c("gene_id", "symbol", "entrez",
                              "baseMean", "log2FoldChange", "lfcSE",
                              "pvalue", "padj")]
write.csv(out, "airway_dex_vs_untrt_deseq2.csv", row.names = FALSE)

norm_counts <- counts(dds, normalized = TRUE)
write.csv(norm_counts, "airway_normalized_counts.csv")

saveRDS(dds, "airway_dds.rds")
writeLines(capture.output(sessionInfo()), "sessionInfo.txt")

Saving the dds object means you can make new plots or contrasts later without refitting, and the session info records the package versions that produced the numbers.

Common mistakes

Feeding in normalized counts

Covered above, but it is the most common problem I see, so here it is again. DESeq2 takes raw integer counts. Not CPM, not TPM, not FPKM, not counts from counts(dds, normalized = TRUE) of some earlier object. For transcript-level quantifiers, use tximport with DESeqDataSetFromTximport().

Ignoring the pairing

It is easy to write ~ dex and move on. You can see what that costs by fitting both designs:

dds_unpaired <- dds
design(dds_unpaired) <- ~ dex
dds_unpaired <- DESeq(dds_unpaired)
res_unpaired <- results(dds_unpaired, contrast = c("dex", "trt", "untrt"),
                        alpha = 0.05)

sum(res$padj < 0.05, na.rm = TRUE)
sum(res_unpaired$padj < 0.05, na.rm = TRUE)

Without cell in the model, the baseline differences between cell lines are treated as noise, dispersions go up and power goes down. The same logic applies to batches, donors, sequencing runs or anything else that groups your samples. If it’s part of the experimental structure, it belongs in the design.

Confusing p-value and padj

A gene with pvalue < 0.05 is not “significant” in an experiment that tested tens of thousands of genes. Call genes on padj, and set alpha in results() to the same cutoff you use. If you need a fold change cutoff, use lfcThreshold instead of filtering after the fact.

Outliers and Cook’s distance

For every gene and sample, DESeq2 computes Cook’s distance, which measures how much that single sample influences the fitted coefficients. In conditions with three or more replicates, results() sets the p-value to NA for genes where a sample’s Cook’s distance exceeds a cutoff based on the F distribution. With seven or more replicates per condition, DESeq() goes further and replaces the outlier counts with a trimmed mean, then refits those genes.

boxplot(log10(assays(dds)[["cooks"]]), range = 0, las = 2)

If one sample has systematically higher Cook’s distances than the rest, that sample is the problem, not individual genes. Look at it in the PCA and the distance heatmap. You can turn the flagging off with results(dds, cooksCutoff = FALSE), but do that after looking at the counts for the flagged genes, not as a way to make a favorite gene come back.

Reading the sign backwards

If the reference level wasn’t set, or the contrast was written in the wrong order, every fold change in your table has the wrong sign. Check mcols(res)$description, then check one well-known gene by hand with plotCounts(). In airway, classic glucocorticoid-induced genes such as DUSP1, PER1 and KLF15, along with CRISPLD2 from the title of the original paper, should come out with positive fold changes.

Testing on transformed data

The vst and rlog values are for QC, clustering and heatmaps. Running t-tests or limma on them as a substitute for the DESeq2 test throws away the count model that makes the method work with small samples.

A checklist to copy

Start from raw integer counts (or tximport for Salmon and kallisto).
Put every known grouping (pairs, donors, batches) in the design, with the variable of interest last.
Set the reference level explicitly with relevel().
Pre-filter very low counts for speed.
Run DESeq(), then look at sizeFactors() and plotDispEsts().
Check vst() PCA and sample distances before reading any results.
Call results() with alpha set to your FDR cutoff, and lfcThreshold if you need a minimum effect size.
Shrink fold changes with lfcShrink(type = "apeglm") for plots and rankings.
Look at the counts of any gene you plan to talk about.
Export with gene IDs as a column, and save the dds object and sessionInfo().

The DESeq2 vignette (browseVignettes("DESeq2")) and the Bioconductor RNA-seq workflow by the same authors cover everything here in more depth, including interactions, likelihood ratio tests (DESeq(dds, test = "LRT", reduced = ~ cell)) and time course designs. When you have a question that isn’t in either, the Bioconductor support site is where the package authors actually answer.

References

Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15, 550 (2014).
Himes BE, et al. RNA-Seq transcriptome profiling identifies CRISPLD2 as a glucocorticoid responsive gene that modulates cytokine function in airway smooth muscle cells. PLoS One 9, e99625 (2014).
Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics 35, 2084-2092 (2019).
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B 57, 289-300 (1995).

← All writing