Reading Volcano and MA Plots

Two plots, two questions

Both plots are drawn from the same results table, but they answer different questions:

• The volcano plot asks: "which genes are both big and trustworthy?" It puts effect size on one axis and confidence on the other, so the genes that matter shoot up into the corners.
• The MA plot asks: "is my data healthy, and are changes spread fairly across expression levels?" It plots fold change against how strongly each gene is expressed, which is a great sanity check.

You'll use the volcano to pick out hits and the MA plot to trust your experiment. Both belong in any RNA-seq report.

1Get a results table to plot

We'll rebuild the results from the previous lesson so this notebook stands on its own. We keep both the raw res and the shrunken res_shrunk, because each plot prefers a different one.

library(DESeq2)

base <- "https://raw.githubusercontent.com/owkin/PyDESeq2/main/datasets/synthetic/"
counts <- as.matrix(read.csv(paste0(base, "test_counts.csv"), row.names = 1))
metadata <- read.csv(paste0(base, "test_metadata.csv"), row.names = 1)
metadata$condition <- factor(metadata$condition)
metadata <- metadata[colnames(counts), , drop = FALSE]

dds <- DESeqDataSetFromMatrix(counts, metadata, design = ~ condition)
dds <- DESeq(dds)

res        <- results(dds, contrast = c("condition", "B", "A"))
res_shrunk <- lfcShrink(dds, coef = "condition_B_vs_A", type = "apeglm")

2The MA plot (the health check)

DESeq2 has this one built in. Use the shrunken results so low-count noise doesn't spray dots all over the edges.

plotMA(res_shrunk, ylim = c(-3, 3), main = "MA plot (shrunken)")

How to read it: each dot is a gene. The horizontal axis is the average expression level (low on the left, high on the right). The vertical axis is the log2 fold change. Dots far from the centre line changed a lot; DESeq2 colours the significant ones. A healthy MA plot is a fat cloud hugging the zero line, fanning out a little on the low-expression left side (where there's naturally more noise). If it looks lopsided or tilted, something may be off with normalization.

3The volcano plot (the hit finder)

There's no single built-in function, so we draw it ourselves. The recipe: fold change on the x-axis, and -log10(padj) on the y-axis so that the most significant genes sit highest.

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

# a gene is a "hit" if it's significant AND changed at least two-fold
df$sig <- df$padj < 0.05 & abs(df$log2FoldChange) > 1

plot(df$log2FoldChange, -log10(df$padj),
     pch = 20,
     col = ifelse(df$sig, "#a855f7", "#cbd5e1"),
     xlab = "log2 fold change (B vs A)",
     ylab = "-log10 adjusted p-value",
     main = "Volcano plot")

# threshold guides: significance line and two-fold lines
abline(h = -log10(0.05), lty = 2, col = "#64748b")
abline(v = c(-1, 1), lty = 2, col = "#64748b")

# label the five most significant genes
top <- head(df[order(df$padj), ], 5)
text(top$log2FoldChange, -log10(top$padj), labels = rownames(top), pos = 3, cex = 0.7)

4How to read a volcano plot

Picture the plot as four regions:

• Top-right: strongly up in condition B, and highly significant. Your most interesting "up" genes.
• Top-left: strongly down in B, highly significant. Your most interesting "down" genes.
• Bottom middle: little change and low confidence. The quiet majority of genes.
• Far sides but low down: big fold change but weak significance. Treat these with caution, they're often low-count genes (this is exactly what shrinkage helps with).

The dashed lines are your thresholds: above the horizontal line means padj < 0.05, and outside the vertical lines means at least a two-fold change. The coloured dots in the upper corners are the genes you'd carry forward into the next step, functional enrichment.

Check your understanding

Next in Track 3

Heatmaps and clustering that mean something →