#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 1: R for Next Gen Biologist
# Script: r01_r_basics_genomics.R
# Description: Getting Started with R for Life Scientists:
#              Basic syntax, vectorization, nucleotide frequencies,
#              GC content calculation, and sequence reverse complementation.
# Usage: Rscript r01_r_basics_genomics.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 1\n")
cat(" Foundational R Syntax & Genomic Calculations\n")
cat("=======================================================\n\n")

# 1. R Environment & BiocManager Check
cat("[Step 1] Checking R Version and Bioconductor Manager...\n")
cat(sprintf("R Version: %s\n", getRversion()))
cat(sprintf("Platform : %s\n", R.version$platform))

# Helper to verify BiocManager presence
check_bioc <- function() {
  if (!requireNamespace("BiocManager", quietly = TRUE)) {
    cat("BiocManager is not installed. To install run: install.packages('BiocManager')\n")
  } else {
    cat(sprintf("BiocManager is ready (Version: %s)\n", as.character(packageVersion("BiocManager"))))
  }
}
check_bioc()

# 2. Defining Genomic Sequences in R
cat("\n[Step 2] Defining Sequences and Basic Variables...\n")
sample_dna <- "ATGCGATCGATCGTAGCTAGCTAGCTAGCTAATCGATCGATCGATCGTAAACCCGGGTTT"
gene_name <- "TP53_Exon_Fragment"
organism <- "Homo sapiens"

cat(sprintf("Gene ID    : %s\n", gene_name))
cat(sprintf("Organism   : %s\n", organism))
cat(sprintf("Sequence   : %s\n", sample_dna))
cat(sprintf("Length     : %d bp\n", nchar(sample_dna)))

# 3. Calculating Nucleotide Composition
cat("\n[Step 3] Nucleotide Composition Analysis...\n")
# Split sequence into single nucleotide character vector
bases <- strsplit(sample_dna, "")[[1]]
base_counts <- table(factor(bases, levels = c("A", "C", "G", "T")))
total_bases <- length(bases)

cat("Nucleotide Frequencies:\n")
for (base in names(base_counts)) {
  pct <- (base_counts[base] / total_bases) * 100
  cat(sprintf("  Base %s: %2d occurrences (%5.2f%%)\n", base, base_counts[base], pct))
}

# 4. GC Content Calculation Function
calc_gc_content <- function(dna_seq) {
  seq_clean <- toupper(gsub("[^ACGTN]", "", dna_seq))
  if (nchar(seq_clean) == 0) return(0.0)
  bases_vec <- strsplit(seq_clean, "")[[1]]
  g_count <- sum(bases_vec == "G")
  c_count <- sum(bases_vec == "C")
  gc_pct <- ((g_count + c_count) / length(bases_vec)) * 100
  return(round(gc_pct, 2))
}

gc_pct <- calc_gc_content(sample_dna)
cat(sprintf("\nCalculated GC Content: %5.2f%%\n", gc_pct))

# 5. Reverse Complement Generation
cat("\n[Step 4] Generating Reverse Complement (Antiparallel 5'->3')...\n")
dna_reverse_complement <- function(dna_seq) {
  # Standard IUPAC complement map
  comp_map <- c(A="T", T="A", C="G", G="C", N="N", a="t", t="a", c="g", g="c", n="n")
  bases_vec <- strsplit(dna_seq, "")[[1]]
  comp_bases <- comp_map[bases_vec]
  rc_seq <- paste(rev(comp_bases), collapse = "")
  return(rc_seq)
}

rev_comp <- dna_reverse_complement(sample_dna)
cat(sprintf("Original (5'->3'): %s\n", sample_dna))
cat(sprintf("RevComp  (5'->3'): %s\n", rev_comp))

# 6. Molecular Weight Approximation
calc_mw <- function(dna_seq) {
  # Average molecular weights of dNTPs (g/mol):
  # A: 313.21, T: 304.2, C: 289.18, G: 329.21
  mw_map <- c(A = 313.21, T = 304.20, C = 289.18, G = 329.21)
  bases_vec <- strsplit(toupper(dna_seq), "")[[1]]
  total_mw <- sum(mw_map[bases_vec], na.rm = TRUE) - 61.96 # Minus water & phosphate terminal loss
  return(round(total_mw, 2))
}

cat(sprintf("Approx. Molecular Weight: %s g/mol (ssDNA)\n", format(calc_mw(sample_dna), big.mark = ",")))

cat("\n[SUCCESS] Chapter 1 foundational script executed successfully.\n")