#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 4: R for Next Gen Biologist
# Script: r04_tidyverse_wrangling.R
# Description: The Tidyverse for Biological Data Wrangling:
#              Pipes, dplyr filtering, mutating, grouping,
#              reshaping count matrices (pivot_longer/wider), and joining omics tables.
# Usage: Rscript r04_tidyverse_wrangling.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 4\n")
cat(" Data Wrangling with the Tidyverse & Base R Equivalents\n")
cat("=======================================================\n\n")

# Check if dplyr/tidyr are available, fallback gracefully to clean base R logic
has_tidyverse <- requireNamespace("dplyr", quietly = TRUE) && requireNamespace("tidyr", quietly = TRUE)

# 1. Generating Mock Differential Expression & Clinical Metadata Tables
set.seed(42)
genes <- c("TP53", "BRCA1", "EGFR", "PTEN", "MYC", "KRAS", "CDK4", "RB1", "PIK3CA", "ATM")
chromosomes <- c("chr17", "chr17", "chr7", "chr10", "chr8", "chr12", "chr12", "chr13", "chr3", "chr11")
log2FC <- c(2.8, 3.4, -2.1, -1.9, 0.4, 1.8, 2.2, -0.2, 1.5, -0.3)
baseMean <- c(1450, 890, 3420, 210, 4800, 1920, 560, 3100, 780, 2400)
padj <- c(1.2e-8, 4.5e-12, 2.1e-5, 3.8e-4, 0.45, 1.1e-3, 4.2e-5, 0.78, 2.9e-3, 0.62)

de_results <- data.frame(
  gene_symbol = genes,
  chr = chromosomes,
  baseMean = baseMean,
  log2FoldChange = log2FC,
  padj = padj,
  stringsAsFactors = FALSE
)

cat("[Step 1] Initial Differential Expression Table:\n")
print(head(de_results, 6))

# 2. Filtering & Mutating (Selecting Significant Hits)
cat("\n[Step 2] Filtering Significant Drivers (padj < 0.01 & |log2FC| >= 1.5)...\n")

if (has_tidyverse) {
  library(dplyr, quietly = TRUE)
  sig_genes <- de_results %>%
    filter(padj < 0.01, abs(log2FoldChange) >= 1.5) %>%
    mutate(
      direction = ifelse(log2FoldChange > 0, "UPREGULATED", "DOWNREGULATED"),
      neg_log10_padj = -log10(padj)
    ) %>%
    arrange(desc(abs(log2FoldChange)))
} else {
  # Base R equivalent
  sig_genes <- de_results[de_results$padj < 0.01 & abs(de_results$log2FoldChange) >= 1.5, ]
  sig_genes$direction <- ifelse(sig_genes$log2FoldChange > 0, "UPREGULATED", "DOWNREGULATED")
  sig_genes$neg_log10_padj <- -log10(sig_genes$padj)
  sig_genes <- sig_genes[order(-abs(sig_genes$log2FoldChange)), ]
}

print(sig_genes)

# 3. Grouping & Summarizing by Chromosome
cat("\n[Step 3] Grouping and Aggregating by Chromosome...\n")

if (has_tidyverse) {
  chr_summary <- sig_genes %>%
    group_by(chr) %>%
    summarise(
      gene_count = n(),
      mean_log2FC = mean(log2FoldChange),
      max_significance = max(neg_log10_padj)
    )
} else {
  chr_summary <- aggregate(
    cbind(log2FoldChange, neg_log10_padj) ~ chr, 
    data = sig_genes, 
    FUN = function(x) c(mean = mean(x), count = length(x), max = max(x))
  )
}
print(chr_summary)

# 4. Reshaping: Wide to Long Expression Matrix
cat("\n[Step 4] Reshaping Expression Data from Wide to Long Format...\n")
wide_counts <- data.frame(
  gene_id = c("TP53", "BRCA1", "EGFR"),
  Ctrl_Rep1 = c(120, 45, 890),
  Ctrl_Rep2 = c(135, 52, 910),
  Treat_Rep1 = c(450, 210, 310),
  Treat_Rep2 = c(480, 225, 290)
)
cat("Wide Matrix (Standard Count Table):\n")
print(wide_counts)

if (has_tidyverse) {
  library(tidyr, quietly = TRUE)
  long_counts <- wide_counts %>%
    pivot_longer(
      cols = -gene_id,
      names_to = "Sample_ID",
      values_to = "Read_Count"
    ) %>%
    mutate(
      Condition = ifelse(grepl("^Ctrl", Sample_ID), "Control", "Treated")
    )
} else {
  long_counts <- reshape(
    wide_counts,
    direction = "long",
    varying = list(2:5),
    v.names = "Read_Count",
    timevar = "Sample_ID",
    times = colnames(wide_counts)[2:5],
    idvar = "gene_id"
  )
  rownames(long_counts) <- NULL
  long_counts$Condition <- ifelse(grepl("^Ctrl", long_counts$Sample_ID), "Control", "Treated")
}

cat("\nLong Matrix (Ideal for ggplot2 & Linear Models):\n")
print(head(long_counts, 8))

# 5. Joining Genomic Tables (Annotating Genes with Pathways)
cat("\n[Step 5] Joining Differential Expression with Pathway Annotations...\n")
pathways <- data.frame(
  gene_symbol = c("TP53", "BRCA1", "EGFR", "PTEN", "MYC", "NOTCH1"),
  pathway = c("p53 Signaling", "Homologous Recombination", "RTK-RAS Pathway", "PI3K-Akt Signaling", "Cell Cycle", "Notch Signaling"),
  stringsAsFactors = FALSE
)

if (has_tidyverse) {
  annotated_hits <- sig_genes %>%
    left_join(pathways, by = "gene_symbol")
} else {
  annotated_hits <- merge(sig_genes, pathways, by = "gene_symbol", all.x = TRUE)
}

print(annotated_hits[, c("gene_symbol", "direction", "log2FoldChange", "pathway")])

cat("\n[SUCCESS] Chapter 4 data wrangling script executed successfully.\n")