#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 6: R for Next Gen Biologist
# Script: r06_ggplot2_publication_figures.R
# Description: Publication-Quality Biological Visualizations with ggplot2:
#              Volcano plots, jittered boxplots for gene expression,
#              custom themes, and high-DPI publication vector export (PDF/PNG).
# Usage: Rscript r06_ggplot2_publication_figures.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 6\n")
cat(" Publication Figures: Volcano Plots & Expression Boxplots\n")
cat("=======================================================\n\n")

temp_dir <- tempdir()
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)

# 1. Generating Mock Differential Expression Dataset for Volcano Plot
set.seed(42)
n_genes <- 1000
log2FC <- rnorm(n_genes, mean = 0, sd = 1.2)
p_vals <- rbeta(n_genes, shape1 = 0.5, shape2 = 5) # Heavy tail towards 0
p_adj <- p.adjust(p_vals, method = "BH")

gene_names <- paste0("GENE_", sprintf("%04d", 1:n_genes))
# Give recognizable names to top outliers
gene_names[1:4] <- c("TP53", "EGFR", "CDKN2A", "MYC")
log2FC[1:4] <- c(3.5, -2.8, -3.2, 4.1)
p_adj[1:4] <- c(1e-12, 1e-9, 1e-11, 1e-15)

volcano_df <- data.frame(
  Gene = gene_names,
  log2FC = log2FC,
  padj = p_adj,
  neg_log10_padj = -log10(p_adj),
  stringsAsFactors = FALSE
)

# Categorize significance
volcano_df$Significance <- "Not Significant"
volcano_df$Significance[volcano_df$padj < 0.05 & volcano_df$log2FC >= 1.5] <- "Upregulated"
volcano_df$Significance[volcano_df$padj < 0.05 & volcano_df$log2FC <= -1.5] <- "Downregulated"
volcano_df$Significance <- factor(volcano_df$Significance, levels = c("Downregulated", "Not Significant", "Upregulated"))

cat("[Step 1] Volcano Plot Data Summary:\n")
print(table(volcano_df$Significance))

# 2. Rendering Volcano Plot
cat("\n[Step 2] Generating Volcano Plot...\n")
volcano_pdf <- file.path(temp_dir, "volcano_plot_publication.pdf")
volcano_png <- file.path(temp_dir, "volcano_plot_publication.png")

if (has_ggplot2) {
  library(ggplot2, quietly = TRUE)
  
  color_map <- c("Downregulated" = "#0284C7", "Not Significant" = "#94A3B8", "Upregulated" = "#EF4444")
  
  p_volcano <- ggplot(volcano_df, aes(x = log2FC, y = neg_log10_padj, color = Significance)) +
    geom_point(alpha = 0.65, size = 1.8) +
    scale_color_manual(values = color_map) +
    geom_vline(xintercept = c(-1.5, 1.5), linetype = "dashed", color = "#475569", linewidth = 0.5) +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "#475569", linewidth = 0.5) +
    labs(
      title = "Transcriptomic Differential Expression Volcano Plot",
      subtitle = "Treatment vs Control (FDR < 0.05, |log2FC| >= 1.5)",
      x = expression(Log[2] ~ "Fold Change"),
      y = expression(-Log[10] ~ "Adjusted P-Value (BH)"),
      caption = "Source: Omics Portal Bioinformatics Suite"
    ) +
    theme_minimal(base_size = 12) +
    theme(
      legend.position = "top",
      panel.grid.minor = element_blank(),
      plot.title = element_text(face = "bold", size = 14),
      plot.subtitle = element_text(color = "#64748B", size = 10)
    )
  
  ggsave(volcano_png, plot = p_volcano, width = 7, height = 6, dpi = 300)
  cat(sprintf("High-resolution PNG saved: %s\n", volcano_png))
  
} else {
  # Base R high-quality fallback
  png(volcano_png, width = 2100, height = 1800, res = 300)
  cols <- ifelse(volcano_df$Significance == "Upregulated", "#EF4444",
          ifelse(volcano_df$Significance == "Downregulated", "#0284C7", "#94A3B8"))
  plot(volcano_df$log2FC, volcano_df$neg_log10_padj, col = cols, pch = 16, cex = 0.8,
       xlab = "Log2 Fold Change", ylab = "-Log10 Adjusted P-Value",
       main = "Volcano Plot (Base R High-Res)")
  abline(v = c(-1.5, 1.5), lty = 2, col = "#475569")
  abline(h = -log10(0.05), lty = 2, col = "#475569")
  dev.off()
  cat(sprintf("Base R PNG saved: %s\n", volcano_png))
}

# 3. Generating Expression Boxplot with Jitter
cat("\n[Step 3] Generating Gene Expression Boxplot across Clinical Stages...\n")
boxplot_png <- file.path(temp_dir, "expression_boxplot_publication.png")

# Mock expression values across 60 clinical samples
expr_df <- data.frame(
  Sample = paste0("S", 1:60),
  Stage = factor(rep(c("Stage I", "Stage II", "Stage III", "Stage IV"), each = 15),
                 levels = c("Stage I", "Stage II", "Stage III", "Stage IV")),
  TP53_TPM = c(rnorm(15, 45, 6), rnorm(15, 38, 5), rnorm(15, 26, 7), rnorm(15, 14, 4))
)

if (has_ggplot2) {
  p_box <- ggplot(expr_df, aes(x = Stage, y = TP53_TPM, fill = Stage)) +
    geom_boxplot(alpha = 0.7, outlier.shape = NA, width = 0.5) +
    geom_jitter(width = 0.15, alpha = 0.8, size = 2, color = "#1E293B") +
    scale_fill_brewer(palette = "Set2") +
    labs(
      title = "TP53 Expression Progression across Tumor Stages",
      y = "Normalized Expression (TPM)",
      x = "Pathological Tumor Stage"
    ) +
    theme_light(base_size = 12) +
    theme(legend.position = "none", plot.title = element_text(face = "bold"))
  
  ggsave(boxplot_png, plot = p_box, width = 6.5, height = 5, dpi = 300)
  cat(sprintf("Expression Boxplot saved: %s\n", boxplot_png))
} else {
  png(boxplot_png, width = 1950, height = 1500, res = 300)
  boxplot(TP53_TPM ~ Stage, data = expr_df, col = "#38BDF8",
          main = "TP53 Expression across Stages", ylab = "TPM", xlab = "Stage")
  dev.off()
  cat(sprintf("Base R Boxplot saved: %s\n", boxplot_png))
}

cat("\n[SUCCESS] Chapter 6 visualization script executed successfully.\n")