10. Publication-Grade Visualizations & STRING Functional Networks

Transforming statistical tables into publication-quality figures and biological insights:

PCA, Volcano Plot & Hierarchical Heatmap in R

library(EnhancedVolcano)
library(pheatmap)
library(RColorBrewer)

# 1. Variance Stabilizing Transformation (VST) for Visualization
vsd <- vst(dds, blind = FALSE)

# 2. Publication-Quality Principal Component Analysis (PCA)
pcaData <- plotPCA(vsd, intgroup = c("Condition", "Batch"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

p_pca <- ggplot(pcaData, aes(x = PC1, y = PC2, color = Condition, shape = Batch)) +
  geom_point(size = 4, alpha = 0.85) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  theme_bw(base_size = 14) +
  scale_color_manual(values = c("Control" = "#06b6d4", "Treated" = "#f43f5e"))
ggsave("06_deseq2/PCA_plot_publication.pdf", p_pca, width = 7, height = 5)

# 3. High-Impact Volcano Plot
pdf("06_deseq2/Volcano_plot.pdf", width = 8, height = 7)
EnhancedVolcano(res_shrink,
  lab = rownames(res_shrink),
  x = 'log2FoldChange',
  y = 'padj',
  pCutoff = 0.05,
  FCcutoff = 1.0,
  pointSize = 2.0,
  labSize = 4.0,
  title = 'Transcriptomic Response: Treated vs. Control',
  subtitle = 'Differential Expression (apeglm shrinkage)',
  col = c('grey60', '#10b981', '#06b6d4', '#f43f5e')
)
dev.off()

# 4. Top 50 Significant DEGs Hierarchical Clustering Heatmap
top50_genes <- head(rownames(sig_degs[order(sig_degs$padj), ]), 50)
mat <- assay(vsd)[top50_genes, ]
mat <- mat - rowMeans(mat) # Mean-center z-scores

pheatmap(mat,
  annotation_col = as.data.frame(colData(dds)[, c("Condition")]),
  color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100),
  show_colnames = TRUE,
  fontsize_row = 7,
  filename = "06_deseq2/Heatmap_top50_DEGs.pdf"
)

Downstream Functional Network Analysis: STRING & clusterProfiler

A gene list is not a biological conclusion. We feed the list of significantly upregulated and downregulated genes into STRING-DB and clusterProfiler to extract protein-protein interaction (PPI) subnetworks, molecular hub bottlenecks, Gene Ontology (Biological Process, Molecular Function, Cellular Component) terms, and KEGG/Reactome signaling pathways:

library(clusterProfiler)
library(org.Hs.eg.db)

# Extract significant upregulated Entrez IDs
up_genes <- rownames(subset(sig_degs, log2FoldChange > 1))
up_entrez <- bitr(up_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

# GO Biological Process Over-Representation Enrichment
ego <- enrichGO(
  gene          = up_entrez,
  OrgDb         = org.Hs.eg.db,
  ont           = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.01,
  qvalueCutoff  = 0.05,
  readable      = TRUE
)

# Dotplot of Top Enriched Pathways
dotplot(ego, showCategory = 15, title = "Enriched GO Biological Processes")
ggsave("06_deseq2/GO_Enrichment_Dotplot.pdf", width = 8, height = 6)