gabrielmoreira Version Compatibility
Reference examples tested with: DESeq2 1.42+, edgeR 4.0+, ggplot2 3.5+, limma 3.58+, matplotlib 3.8+
Before using code patterns, verify installed versions match. If versions differ:
- R:
packageVersion('<pkg>')then?function_nameto 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.
DE Visualization
Create visualizations for differential expression analysis using DESeq2 and edgeR built-in plotting functions.
Scope
This skill covers DE-specific built-in functions:
- DESeq2:
plotMA(),plotPCA(),plotDispEsts(),plotCounts() - edgeR:
plotMD(),plotBCV(),plotMDS() - Sample distance heatmaps and p-value distributions
For custom ggplot2/matplotlib implementations of volcano and MA plots, see data-visualization/volcano-and-ma-plots; for PCA / UMAP / t-SNE / PHATE, see data-visualization/dimensionality-reduction-plots.
Required Libraries
library(DESeq2)
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
library(ggrepel) # For labeled points
Installation
install.packages(c('ggplot2', 'pheatmap', 'RColorBrewer', 'ggrepel'))
# Optional: Enhanced volcano plots
BiocManager::install('EnhancedVolcano')
MA Plot
Goal: Visualize the relationship between mean expression and log fold change to assess DE results.
Approach: Plot log fold change against mean normalized counts, highlighting significant genes.
"Make an MA plot of my DE results" → Plot mean expression vs. fold change with significant genes colored, using plotMA or ggplot2.
DESeq2 MA Plot
# Built-in MA plot
plotMA(res, ylim = c(-5, 5), main = 'MA Plot')
# With custom alpha
plotMA(res, alpha = 0.05, ylim = c(-5, 5))
# Highlight specific genes
plotMA(res, ylim = c(-5, 5))
with(subset(res, padj < 0.01 & abs(log2FoldChange) > 2),
points(baseMean, log2FoldChange, col = 'red', pch = 20))
Custom ggplot2 MA Plot
res_df <- as.data.frame(res)
res_df$significant <- res_df$padj < 0.05 & !is.na(res_df$padj)
ggplot(res_df, aes(x = log10(baseMean), y = log2FoldChange, color = significant)) +
geom_point(alpha = 0.5, size = 1) +
scale_color_manual(values = c('grey60', 'red')) +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(x = 'log10(Mean Expression)', y = 'log2 Fold Change', title = 'MA Plot') +
theme_bw() +
theme(legend.position = 'bottom')
edgeR MA Plot
# Using plotMD (mean-difference plot)
plotMD(qlf, main = 'MD Plot')
abline(h = c(-1, 1), col = 'blue', lty = 2)
Volcano Plot
Goal: Display statistical significance against fold change magnitude to identify the most important DE genes.
Approach: Plot -log10(p-value) vs. log2 fold change with threshold lines and optional gene labels.
"Create a volcano plot of differentially expressed genes" → Scatter plot of fold change vs. significance with colored significance regions and labeled top hits.
Basic Volcano Plot
res_df <- as.data.frame(res)
res_df$significant <- res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1
ggplot(res_df, aes(x = log2FoldChange, y = -log10(pvalue), color = significant)) +
geom_point(alpha = 0.5, size = 1) +
scale_color_manual(values = c('grey60', 'red')) +
geom_vline(xintercept = c(-1, 1), linetype = 'dashed', color = 'blue') +
geom_hline(yintercept = -log10(0.05), linetype = 'dashed', color = 'blue') +
labs(x = 'log2 Fold Change', y = '-log10(p-value)', title = 'Volcano Plot') +
theme_bw()
Volcano with Gene Labels
res_df <- as.data.frame(res)
res_df$gene <- rownames(res_df)
res_df$significant <- res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1
# Label top genes
top_genes <- head(res_df[order(res_df$padj), ], 10)
ggplot(res_df, aes(x = log2FoldChange, y = -log10(pvalue))) +
geom_point(aes(color = significant), alpha = 0.5, size = 1) +
scale_color_manual(values = c('grey60', 'red')) +
geom_text_repel(data = top_genes, aes(label = gene),
size = 3, max.overlaps = 20) +
geom_vline(xintercept = c(-1, 1), linetype = 'dashed') +
geom_hline(yintercept = -log10(0.05), linetype = 'dashed') +
labs(x = 'log2 Fold Change', y = '-log10(p-value)') +
theme_bw()
EnhancedVolcano
library(EnhancedVolcano)
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
pCutoff = 0.05,
FCcutoff = 1,
title = 'Differential Expression',
subtitle = 'Treatment vs Control')
PCA Plot
Goal: Assess sample clustering and identify batch effects or outliers via dimensionality reduction.
Approach: Apply variance-stabilizing transformation then project samples onto principal components, coloring by experimental variables.
"Show me a PCA plot of my samples" → Perform PCA on transformed expression data and visualize sample separation by condition and batch.
DESeq2 PCA
# Variance stabilizing transformation first
vsd <- vst(dds, blind = FALSE)
# Basic PCA
plotPCA(vsd, intgroup = 'condition')
# With more options
plotPCA(vsd, intgroup = c('condition', 'batch'), ntop = 500)
Custom PCA with ggplot2
vsd <- vst(dds, blind = FALSE)
pca_data <- plotPCA(vsd, intgroup = c('condition', 'batch'), returnData = TRUE)
percentVar <- round(100 * attr(pca_data, 'percentVar'))
ggplot(pca_data, aes(x = PC1, y = PC2, color = condition, shape = batch)) +
geom_point(size = 4) +
xlab(paste0('PC1: ', percentVar[1], '% variance')) +
ylab(paste0('PC2: ', percentVar[2], '% variance')) +
ggtitle('PCA Plot') +
theme_bw() +
theme(legend.position = 'right')
edgeR PCA (via limma)
library(limma)
log_cpm <- cpm(y, log = TRUE)
plotMDS(log_cpm, col = as.numeric(group), pch = 16)
legend('topright', legend = levels(group), col = 1:nlevels(group), pch = 16)
Heatmaps
Goal: Visualize expression patterns of significant genes across samples to reveal clusters and condition effects.
Approach: Z-score normalize VST-transformed counts for significant genes and cluster with pheatmap, annotating by condition.
"Make a heatmap of the top differentially expressed genes" → Extract significant genes, z-score normalize, and create a clustered heatmap with sample annotations.
Top DE Genes Heatmap
library(pheatmap)
# Get top significant genes
sig_genes <- rownames(subset(res, padj < 0.01))
# Get normalized counts
vsd <- vst(dds, blind = FALSE)
mat <- assay(vsd)[sig_genes, ]
# Scale by row (z-score)
mat_scaled <- t(scale(t(mat)))
# Create annotation
annotation_col <- data.frame(
condition = colData(dds)$condition,
row.names = colnames(mat)
)
pheatmap(mat_scaled,
annotation_col = annotation_col,
show_rownames = FALSE,
clustering_distance_rows = 'correlation',
clustering_distance_cols = 'correlation',
color = colorRampPalette(c('blue', 'white', 'red'))(100),
main = 'Top DE Genes')
Sample Distance Heatmap
vsd <- vst(dds, blind = FALSE)
# Calculate sample distances
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
# Annotation
annotation <- data.frame(
condition = colData(dds)$condition,
row.names = colnames(dds)
)
pheatmap(sampleDistMatrix,
annotation_col = annotation,
annotation_row = annotation,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
color = colorRampPalette(c('white', 'steelblue'))(100),
main = 'Sample Distance Matrix')
Gene Expression Heatmap
# Select genes of interest
genes_of_interest <- c('gene1', 'gene2', 'gene3', 'gene4', 'gene5')
mat <- assay(vsd)[genes_of_interest, ]
pheatmap(mat,
scale = 'row',
annotation_col = annotation_col,
show_rownames = TRUE,
cluster_cols = TRUE,
cluster_rows = TRUE,
main = 'Genes of Interest')
Dispersion Plot
Goal: Assess the fit of the dispersion model to verify DE analysis assumptions.
Approach: Plot gene-wise, fitted, and final dispersion estimates against mean expression.
DESeq2
plotDispEsts(dds, main = 'Dispersion Estimates')
edgeR
plotBCV(y, main = 'Biological Coefficient of Variation')
Counts Plot for Individual Genes
Goal: Visualize expression of a specific gene across samples and conditions.
Approach: Extract per-sample counts for a gene and plot by condition using plotCounts or ggplot2.
DESeq2
# Plot counts for a specific gene
plotCounts(dds, gene = 'GENE_NAME', intgroup = 'condition')
# With ggplot2
d <- plotCounts(dds, gene = 'GENE_NAME', intgroup = 'condition', returnData = TRUE)
ggplot(d, aes(x = condition, y = count, color = condition)) +
geom_point(position = position_jitter(width = 0.1), size = 3) +
scale_y_log10() +
ggtitle('GENE_NAME Expression') +
theme_bw()
edgeR
# Get CPM for a gene
gene_idx <- which(rownames(y) == 'GENE_NAME')
cpm_gene <- cpm(y)[gene_idx, ]
# Plot
df <- data.frame(cpm = cpm_gene, group = group)
ggplot(df, aes(x = group, y = cpm, color = group)) +
geom_point(position = position_jitter(width = 0.1), size = 3) +
scale_y_log10() +
labs(y = 'CPM', title = 'GENE_NAME Expression') +
theme_bw()
P-value Histogram
Goal: Diagnose the quality of the DE analysis by examining the raw p-value distribution.
Approach: Histogram of raw p-values; a uniform distribution with a peak near zero indicates a well-calibrated test.
# Check p-value distribution (should be uniform under null with peak near 0)
res_df <- as.data.frame(res)
ggplot(res_df, aes(x = pvalue)) +
geom_histogram(bins = 50, fill = 'steelblue', color = 'white') +
labs(x = 'P-value', y = 'Frequency', title = 'P-value Distribution') +
theme_bw()
Saving Plots
Goal: Export publication-quality plots in vector or raster formats.
Approach: Use pdf/png devices or ggsave with appropriate resolution and dimensions.
# Save as PDF (vector)
pdf('volcano_plot.pdf', width = 8, height = 6)
# ... plot code ...
dev.off()
# Save as PNG (raster)
png('volcano_plot.png', width = 800, height = 600, res = 150)
# ... plot code ...
dev.off()
# Using ggsave for ggplot objects
p <- ggplot(...) + ...
ggsave('plot.pdf', p, width = 8, height = 6)
ggsave('plot.png', p, width = 8, height = 6, dpi = 300)
Color Palettes
Goal: Select appropriate color schemes for heatmaps and categorical data.
Approach: Use RColorBrewer palettes -- diverging for expression, sequential for distances, qualitative for groups.
# For heatmaps
library(RColorBrewer)
# Diverging (for expression: blue-white-red)
colorRampPalette(rev(brewer.pal(n = 7, name = 'RdBu')))(100)
# Sequential (for distances)
colorRampPalette(brewer.pal(n = 9, name = 'Blues'))(100)
# For categorical groups
brewer.pal(n = 8, name = 'Set1')
Quick Reference: Common Plots
| Plot | Purpose | Function |
|---|---|---|
| MA plot | LFC vs mean expression | plotMA(), plotMD() |
| Volcano | LFC vs significance | ggplot2, EnhancedVolcano |
| PCA | Sample clustering | plotPCA(), plotMDS() |
| Heatmap | Gene patterns | pheatmap() |
| Dispersion | Model fit | plotDispEsts(), plotBCV() |
| Counts | Individual genes | plotCounts() |
Diagnostic Interpretation
P-value Histogram
Check raw p-value distribution before trusting DE results:
| Shape | Meaning | Action |
|---|---|---|
| Uniform + spike near 0 | Correct: null genes uniform, true DE near 0 | Proceed normally |
| Anti-conservative (U-shape, spikes at 0 and 1) | Inflated significance; unmodeled batch or violated assumptions | Check for batch effects, verify model specification |
| Conservative (depleted near 0, spike near 1) | Over-correction; too many covariates or wrong dispersion | Simplify model, check dispersion plot |
| Spike at p = 1 only | Discrete artifact from low-count genes | Pre-filter more aggressively |
MA Plot Diagnostics
| Pattern | Meaning |
|---|---|
| Symmetric cloud centered at LFC = 0 | Correct normalization |
| Cloud shifted up or down | Normalization failure; majority-DE experiment may violate assumptions |
| Funnel shape widening at low expression | Expected — low-count genes have noisier fold changes |
| Discrete horizontal bands | Low-count artifacts; consider stronger pre-filtering |
Volcano Plot: Shrunken LFCs
Use shrunken LFCs (apeglm/ashr) on the x-axis and un-shrunken p-values on the y-axis. This combination gives stable fold change estimates while preserving the original significance assessment. Without shrinkage, low-count genes with extreme but unreliable fold changes dominate the plot edges.
Dispersion Plot Diagnostics
| Pattern | Meaning |
|---|---|
| Gene-wise points scattered around fitted line | Good model fit |
| Gene-wise points far above fitted line | Possible outlier genes or unmodeled batch effects |
| Fitted line flat (no trend) | Unusual — check if data is over-filtered or has unusual structure |
| Final estimates much lower than gene-wise | Expected — shrinkage toward the fitted trend |
PCA Diagnostics
| Pattern | Meaning | Action |
|---|---|---|
| Clear separation by condition on PC1/PC2 | Strong biological signal | Proceed |
| Separation by batch, not condition | Batch effect dominates | Include batch in model; do NOT use corrected counts for DE |
| One sample far from its group | Potential outlier or sample swap | Check library QC metrics; consider removing |
| No separation on PC1/PC2 but present on PC3+ | Subtle effects | May still find DE genes; check dispersion estimates |
Related Skills
- deseq2-basics - Generate DESeq2 results for visualization
- edger-basics - Generate edgeR results for visualization
- de-results - Filter genes before visualization
- data-visualization/volcano-and-ma-plots - Custom ggplot2/matplotlib volcano + MA with LFC shrinkage
- data-visualization/dimensionality-reduction-plots - PCA / UMAP / t-SNE / PHATE
- data-visualization/heatmaps-clustering - Advanced heatmap customization