#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 11: R for Next Gen Biologist
# Script: r11_pathway_enrichment_clusterprofiler.R
# Description: Functional Enrichment Analysis & Pathway Mining:
#              Over-Representation Analysis (ORA), Hypergeometric Fisher's test,
#              Gene Ontology (GO) terms, KEGG pathways, and GSEA ranking in R.
# Usage: Rscript r11_pathway_enrichment_clusterprofiler.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 11\n")
cat(" Pathway Enrichment Analysis & Gene Ontology Mining\n")
cat("=======================================================\n\n")

has_cp <- requireNamespace("clusterProfiler", quietly = TRUE)

# 1. Candidate Differentially Expressed Gene Set
cat("[Step 1] Query Gene Set (Significantly Upregulated in Tumor Cohort)...\n")
candidate_genes <- c("TP53", "CDKN1A", "MDM2", "BAX", "PUMA", "ATM", "CHEK2", "FAS", "APAF1", "GADD45A")
cat(sprintf("Input Query Genes (%d total): %s\n", length(candidate_genes), paste(candidate_genes, collapse = ", ")))

# 2. Defining Curated Biological Pathways (Gene Sets)
pathway_database <- list(
  "p53 Signaling Pathway" = c("TP53", "CDKN1A", "MDM2", "BAX", "PUMA", "ATM", "CHEK2", "FAS", "APAF1", "GADD45A", "CCND1", "CDK4", "RPRM", "SIAH1"),
  "Apoptosis & Cell Death" = c("BAX", "FAS", "APAF1", "CASP3", "CASP8", "CASP9", "BCL2", "CYCS", "BID", "PUMA", "NOXA"),
  "Cell Cycle G1/S Transition" = c("CDKN1A", "CCND1", "CDK4", "CDK6", "RB1", "E2F1", "CCNE1", "CDKN2A"),
  "DNA Double-Strand Break Repair" = c("ATM", "CHEK2", "BRCA1", "BRCA2", "RAD51", "TP53BP1", "MRE11", "NBN", "RAD50"),
  "PI3K-Akt Signaling Pathway" = c("PIK3CA", "AKT1", "PTEN", "MTOR", "FOXO3", "BAD", "GSK3B", "TSC1", "TSC2"),
  "Wnt Signaling Pathway" = c("CTNNB1", "APC", "AXIN1", "GSK3B", "WNT1", "WNT3A", "FZD1", "LRP6")
)

# Genome background size (all measured human protein-coding genes)
genome_universe_size <- 20000

# 3. Custom Hypergeometric / Fisher's Exact Test Engine (ORA)
cat("\n[Step 2] Executing Over-Representation Analysis (Hypergeometric / Fisher)...\n")
ora_results <- data.frame(
  Pathway = character(),
  Overlap = integer(),
  Pathway_Size = integer(),
  Fold_Enrichment = numeric(),
  p_value = numeric(),
  stringsAsFactors = FALSE
)

n_query <- length(candidate_genes)

for (pathway_name in names(pathway_database)) {
  pathway_genes <- pathway_database[[pathway_name]]
  overlap_genes <- intersect(candidate_genes, pathway_genes)
  k <- length(overlap_genes) # White balls drawn
  m <- length(pathway_genes) # White balls in urn
  n <- genome_universe_size - m # Black balls in urn
  
  if (k == 0) next
  
  # Hypergeometric p-value P(X >= k)
  p_val <- phyper(k - 1, m, n, n_query, lower.tail = FALSE)
  
  # Expected overlap by random chance = (m / universe) * n_query
  expected_overlap <- (m / genome_universe_size) * n_query
  fold_enrichment  <- round(k / expected_overlap, 1)
  
  ora_results <- rbind(ora_results, data.frame(
    Pathway = pathway_name,
    Overlap = k,
    Pathway_Size = m,
    Fold_Enrichment = fold_enrichment,
    p_value = p_val,
    stringsAsFactors = FALSE
  ))
}

# Adjust p-values for multiple hypothesis testing (Benjamini-Hochberg)
ora_results$FDR <- p.adjust(ora_results$p_value, method = "BH")
ora_results <- ora_results[order(ora_results$FDR), ]

cat("\nOver-Representation Analysis Results Table:\n")
for (i in 1:nrow(ora_results)) {
  cat(sprintf("  %-32s | Overlap: %2d/%2d | Fold: %5.1fx | p-val: %.2e | FDR: %.2e\n",
              ora_results$Pathway[i],
              ora_results$Overlap[i],
              ora_results$Pathway_Size[i],
              ora_results$Fold_Enrichment[i],
              ora_results$p_value[i],
              ora_results$FDR[i]))
}

# 4. GSEA Concept Overview
cat("\n[Step 3] Gene Set Enrichment Analysis (GSEA) Methodology:\n")
cat("  1. Continuous Ranking Metric: Compute signal-to-noise ratio or signed -log10(p) * sign(log2FC)\n")
cat("  2. Walking Down the Ranked List: Calculate running sum score (Enrichment Score, ES)\n")
cat("  3. Permutation Testing: Permute sample condition labels to estimate null distribution\n")
cat("  4. Normalized Enrichment Score (NES): Normalize ES across different gene set sizes\n")

cat("\n[SUCCESS] Chapter 11 pathway enrichment script executed successfully.\n")