Lesson 11 of 13 · 12 min

Visualizing Results: Volcano and MA Plots

DESeq2 gives you a results table with thousands of genes, but a table hides the shape of the experiment. The MA plot and volcano plot summarize per-gene effects, while PCA and a sample-distance heatmap check whether the samples themselves behave. This lesson assumes a DESeqDataSet named dds that you already ran through DESeq().

library(DESeq2)

# dds was already built and run: dds <- DESeq(dds)
resultsNames(dds)

res <- results(dds, contrast = c("condition", "treated", "control"), alpha = 0.05)
summary(res)

[1] "Intercept"                    "condition_treated_vs_control"

out of 18324 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 1284, 7%
LFC < 0 (down)     : 1121, 6.1%
outliers [1]       : 42, 0.23%
low counts [2]     : 3556, 19%
(mean count < 5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

MA plots

An MA plot shows log2 fold change on the y-axis against the mean of normalized counts on the x-axis, one point per gene. Low-count genes on the left have noisy, inflated fold changes, so DESeq2 recommends shrinking the log2 fold changes with lfcShrink before plotting. Shrinkage pulls unreliable estimates toward zero while leaving well-measured genes essentially unchanged.

library(apeglm)

# apeglm shrinkage uses the coefficient name shown by resultsNames(dds)
resLFC <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")

plotMA(resLFC, ylim = c(-4, 4), alpha = 0.05)
abline(h = c(-1, 1), col = "dodgerblue", lty = 2)

Volcano plots

library(ggplot2)

df <- as.data.frame(res)
df <- df[!is.na(df$padj), ]

# classify each gene by threshold: |log2FC| > 1 and padj < 0.05
df$status <- "NS"
df$status[df$padj < 0.05 & df$log2FoldChange >  1] <- "Up"
df$status[df$padj < 0.05 & df$log2FoldChange < -1] <- "Down"

ggplot(df, aes(log2FoldChange, -log10(padj), color = status)) +
  geom_point(alpha = 0.6, size = 1) +
  scale_color_manual(values = c(Up = "#d73027", Down = "#4575b4", NS = "grey70")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  labs(x = "log2 fold change", y = "-log10 adjusted p-value", color = NULL) +
  theme_bw()

Threshold on padj, the Benjamini-Hochberg adjusted p-value, not the raw pvalue column. DESeq2 sets padj to NA for genes dropped by independent filtering or flagged as Cook's-distance outliers, so remove those NA rows before you plot or count hits, exactly as the volcano code does above.

Sample-level QC

library(pheatmap)

# variance-stabilizing transform; blind = FALSE keeps the design in dispersion
# estimation so real group effects don't inflate the mean-variance trend
vsd <- vst(dds, blind = FALSE)

# PCA on the transformed counts, colored by condition
plotPCA(vsd, intgroup = "condition")

# sample-to-sample Euclidean distances on transposed matrix (rows = samples)
sampleDists <- dist(t(assay(vsd)))
pheatmap(as.matrix(sampleDists),
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         main = "Sample-to-sample distances")

Try it yourself: extend the volcano plot to label the strongest hits. Add a gene column with df$gene <- rownames(df), then split by status and take the five smallest-padj genes from each of the Up and Down groups into a data frame called top_genes. Layer ggrepel::geom_text_repel(data = top_genes, aes(label = gene)) onto the plot and confirm the labeled points land in the expected upper-left and upper-right corners.

Differential Expression with DESeq2 Functional and GO Enrichment