#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 7: R for Next Gen Biologist
# Script: r07_genomicranges_s4.R
# Description: Bioconductor Core: S4 Classes & GenomicRanges:
#              Understanding S4 objects, GRanges genomic intervals arithmetic,
#              finding overlaps (ChIP-seq peaks vs promoters), and SummarizedExperiment.
# Usage: Rscript r07_genomicranges_s4.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 7\n")
cat(" Bioconductor S4 Architecture & GenomicRanges (GRanges)\n")
cat("=======================================================\n\n")

# Check if GenomicRanges is available
has_granges <- requireNamespace("GenomicRanges", quietly = TRUE) && requireNamespace("IRanges", quietly = TRUE)

if (has_granges) {
  library(GenomicRanges, quietly = TRUE)
  library(IRanges, quietly = TRUE)
  cat("[Bioconductor Mode] GenomicRanges package loaded.\n")
} else {
  cat("[Standalone Mode] Bioconductor GenomicRanges not installed.\n")
  cat("Simulating S4 GRanges interval arithmetic logic natively in R.\n")
}

# 1. Defining Gene Models & Promoters
cat("\n[Step 1] Constructing Gene Annotation Intervals...\n")

if (has_granges) {
  genes_gr <- GRanges(
    seqnames = c("chr1", "chr1", "chr2", "chr17"),
    ranges   = IRanges(start = c(10000, 50000, 25000, 7668402),
                       end   = c(15000, 58000, 32000, 7687550)),
    strand   = c("+", "-", "+", "-"),
    gene_id  = c("GENE_A", "GENE_B", "GENE_C", "TP53"),
    score    = c(100, 250, 80, 500)
  )
  cat("GenomicRanges (Genes):\n")
  print(genes_gr)
  
  # Calculate Promoters (2000bp upstream, 200bp downstream of TSS)
  promoters_gr <- promoters(genes_gr, upstream = 2000, downstream = 200)
  cat("\nPromoter Regions (TSS-relative):\n")
  print(promoters_gr)
  
} else {
  # Base R S4 representation simulation
  genes_df <- data.frame(
    seqnames = c("chr1", "chr1", "chr2", "chr17"),
    start    = c(10000, 50000, 25000, 7668402),
    end      = c(15000, 58000, 32000, 7687550),
    strand   = c("+", "-", "+", "-"),
    gene_id  = c("GENE_A", "GENE_B", "GENE_C", "TP53"),
    stringsAsFactors = FALSE
  )
  
  calc_promoters <- function(df, upstream = 2000, downstream = 200) {
    p_start <- ifelse(df$strand == "+", pmax(1, df$start - upstream), df$end - downstream)
    p_end   <- ifelse(df$strand == "+", df$start + downstream, df$end + upstream)
    data.frame(seqnames = df$seqnames, start = p_start, end = p_end, strand = df$strand, gene_id = df$gene_id)
  }
  
  promoters_df <- calc_promoters(genes_df)
  cat("Promoter Regions Table (Simulated):\n")
  print(promoters_df)
}

# 2. Defining ChIP-seq Peaks and Computing Overlaps
cat("\n[Step 2] Finding Overlaps Between Transcription Factor Peaks and Promoters...\n")

if (has_granges) {
  # Simulated ChIP-seq transcription factor binding peaks
  tf_peaks_gr <- GRanges(
    seqnames = c("chr1", "chr1", "chr17"),
    ranges   = IRanges(start = c(8500, 62000, 7686500),
                       end   = c(9200, 63000, 7688000)),
    peak_id  = c("Peak_01", "Peak_02", "Peak_03"),
    qvalue   = c(14.5, 8.2, 32.1)
  )
  cat("ChIP-seq Peaks:\n")
  print(tf_peaks_gr)
  
  # Find overlaps
  overlaps <- findOverlaps(tf_peaks_gr, promoters_gr)
  cat("\nOverlapping Pairs (Hits Object):\n")
  print(overlaps)
  
  # Annotate peaks
  query_idx <- queryHits(overlaps)
  subj_idx  <- subjectHits(overlaps)
  cat(sprintf("\nMatched Peak '%s' to Promoter of '%s'!\n",
              tf_peaks_gr$peak_id[query_idx],
              promoters_gr$gene_id[subj_idx]))

} else {
  # Native interval overlap algorithm
  peaks_df <- data.frame(
    seqnames = c("chr1", "chr1", "chr17"),
    start    = c(8500, 62000, 7686500),
    end      = c(9200, 63000, 7688000),
    peak_id  = c("Peak_01", "Peak_02", "Peak_03"),
    stringsAsFactors = FALSE
  )
  
  find_overlaps_base <- function(query, subject) {
    matches <- list()
    for (i in 1:nrow(query)) {
      for (j in 1:nrow(subject)) {
        if (query$seqnames[i] == subject$seqnames[j]) {
          # Overlap condition: max(start) <= min(end)
          if (max(query$start[i], subject$start[j]) <= min(query$end[i], subject$end[j])) {
            matches[[length(matches) + 1]] <- data.frame(
              Peak = query$peak_id[i],
              Gene = subject$gene_id[j],
              Overlap_Chr = query$seqnames[i]
            )
          }
        }
      }
    }
    do.call(rbind, matches)
  }
  
  overlap_results <- find_overlaps_base(peaks_df, promoters_df)
  cat("Interval Overlap Hits Table:\n")
  print(overlap_results)
}

# 3. SummarizedExperiment Architecture Overview
cat("\n[Step 3] The SummarizedExperiment Architecture:\n")
cat("  - assays()    : Count matrices (Rows = Genomic Features, Cols = Samples)\n")
cat("  - rowRanges() : GRanges objects containing coordinate intervals for each gene\n")
cat("  - colData()   : Clinical metadata data.frame with phenotype annotations\n")
cat("  - metadata()  : Arbitrary list with pipeline parameters, genome assembly, etc.\n")

cat("\n[SUCCESS] Chapter 7 GenomicRanges script executed successfully.\n")