Chapter 9: Publication-Quality DEG Visualizations: Enhanced Volcano, MA Plots & ComplexHeatmap

1. Visualizing Differential Expression for High-Impact Journals

High-profile biomedical journals (Cell, Nature, Science) demand rigorous, multi-faceted graphical summaries of transcriptomic experiments. Three visual representations form the cornerstone of publication figures:

1. The MA Plot: Plots \(\log_2 \text{Fold Change}\) (\(M\)) against mean expression intensity (\(A = \text{baseMean}\)). Demonstrates that fold changes are symmetric and checks for intensity-dependent bias.
2. The Enhanced Volcano Plot: Combines magnitude of biological change (\(x = \log_2\text{FC}\)) with statistical significance (\(y = -\log_{10} p\)), labeling key driver oncogenes and tumor suppressors.
3. The Complex Heatmap: Displays unsupervised hierarchical clustering of top DEGs across patient cohorts with annotated clinical sidebars (stage, sex, survival).

2. Complete Publication Graphics Script

# Install packages if missing
if (!requireNamespace("EnhancedVolcano", quietly = TRUE)) BiocManager::install("EnhancedVolcano")
if (!requireNamespace("ComplexHeatmap", quietly = TRUE)) BiocManager::install("ComplexHeatmap")

library(EnhancedVolcano)
library(ComplexHeatmap)
library(circlize)

# 1. Publication-Quality Enhanced Volcano Plot
# Labels specific candidate genes of interest
key_genes <- c("ENSG00000000001", "ENSG00000000002", "ENSG00000000003")

EnhancedVolcano(
    res_shrunk,
    lab = rownames(res_shrunk),
    x = 'log2FoldChange',
    y = 'padj',
    selectLab = key_genes,
    pCutoff = 0.05,
    FCcutoff = 1.0,
    pointSize = 2.5,
    labSize = 4.5,
    col = c("grey70", "#06b6d4", "#f59e0b", "#f43f5e"),
    colAlpha = 0.8,
    legendPosition = 'right',
    drawConnectors = TRUE,
    widthConnectors = 0.75,
    title = 'Volcano Plot: Treated vs. Control Cohort',
    subtitle = 'Significance threshold: FDR < 0.05 and |log2FC| >= 1.0'
)

# 2. MA Plot Diagnostic
plotMA(res_shrunk, ylim = c(-3, 3), main = "MA Plot with apeglm LFC Shrinkage")

# 3. Multi-Track ComplexHeatmap of Top 50 DEGs
top50_genes <- res_df %>% filter(Call != "NS") %>% arrange(padj) %>% head(50) %>% pull(Gene_Symbol)
vst_mat_top50 <- assay(vsd)[top50_genes, ]

# Row-wise Z-score transformation
z_score_mat <- t(scale(t(vst_mat_top50)))

# Clinical Column Annotations
col_annot <- HeatmapAnnotation(
    Condition = clinical_clean$Condition,
    Stage = clinical_clean$Tumor_Stage,
    Batch = clinical_clean$Extraction_Batch,
    col = list(
        Condition = c("Control" = "#06b6d4", "Treated" = "#f43f5e"),
        Stage = c("Normal" = "#10b981", "Stage_II" = "#fcd34d", "Stage_III" = "#f97316", "Stage_IV" = "#b91c1c"),
        Batch = c("Batch_1" = "#818cf8", "Batch_2" = "#c084fc", "Batch_3" = "#f472b6")
    )
)

# Color gradient for Z-scores (Blue = Down, White = Mean, Red = Up)
col_fun <- colorRamp2(c(-2, 0, 2), c("#2563eb", "#ffffff", "#dc2626"))

Heatmap(
    z_score_mat,
    name = "Z-Score",
    col = col_fun,
    top_annotation = col_annot,
    show_row_names = FALSE,
    show_column_names = TRUE,
    clustering_distance_rows = "pearson",
    clustering_method_rows = "complete",
    clustering_distance_columns = "euclidean",
    clustering_method_columns = "complete",
    column_title = "Unsupervised Hierarchical Clustering of Top 50 DEGs"
)