#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 5: R for Next Gen Biologist
# Script: r05_genomic_file_io.R
# Description: Biological File I/O & Bioconductor Genomic Formats:
#              Fast tabular reading, parsing FASTA records in R,
#              working with BED & GTF coordinates, and saving/loading RDS objects.
# Usage: Rscript r05_genomic_file_io.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 5\n")
cat(" Biological File I/O & Genomic Format Handlers\n")
cat("=======================================================\n\n")

# Temporary working directory for file demo
temp_dir <- tempdir()
cat(sprintf("Working with temporary sandbox: %s\n", temp_dir))

# 1. Generating & Writing a Mock FASTA File
cat("\n[Step 1] Creating and Writing a Genomic FASTA File...\n")
mock_fasta_path <- file.path(temp_dir, "reference_contigs.fasta")
fasta_lines <- c(
  ">chr17:7668402-7687550 TP53 Tumor Suppressor Exon 1-11",
  "ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCAGACCTATGG",
  "AAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATGGATGATTTGATGCTGTCC",
  ">chr7:55019017-55211628 EGFR Receptor Tyrosine Kinase",
  "ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTGCTGGCTGCGCTCTGCCCGGCGAGTCGG",
  "GCTCTGGAGGAAAAGAAAGTTTGCCAAGGCACGAGTAACAAGCTCACGCAGTTGGGCACTTTTGAAGAT",
  ">chr13:32315474-32400266 BRCA2 DNA Repair Associated",
  "ATGCCTATTGGATCCAAAGAGAGGCCAACATTTTTTGAAATTTTTAAGACACGCTGCAACAAAGCAGAT"
)
writeLines(fasta_lines, mock_fasta_path)
cat(sprintf("Wrote 3 FASTA contigs to: %s\n", mock_fasta_path))

# 2. Parsing FASTA Records in Base R
cat("\n[Step 2] Custom Base R FASTA Parser (Biostrings-independent)...\n")
read_fasta_base <- function(filepath) {
  lines <- readLines(filepath)
  headers <- c()
  sequences <- c()
  current_seq <- ""
  
  for (line in lines) {
    line <- trimws(line)
    if (nchar(line) == 0) next
    if (startsWith(line, ">")) {
      if (current_seq != "") {
        sequences <- c(sequences, current_seq)
        current_seq <- ""
      }
      headers <- c(headers, substring(line, 2))
    } else {
      current_seq <- paste0(current_seq, line)
    }
  }
  if (current_seq != "") {
    sequences <- c(sequences, current_seq)
  }
  
  data.frame(
    Header = headers,
    Sequence = sequences,
    Length = nchar(sequences),
    stringsAsFactors = FALSE
  )
}

parsed_fasta <- read_fasta_base(mock_fasta_path)
print(parsed_fasta[, c("Header", "Length")])

# 3. Generating and Parsing a BED (Browser Extensible Data) File
cat("\n[Step 3] BED Format Specification & Parsing...\n")
mock_bed_path <- file.path(temp_dir, "exons.bed")
bed_data <- data.frame(
  chrom      = c("chr17", "chr17", "chr7", "chr7", "chr13"),
  chromStart = c(7668402, 7673770, 55019017, 55086713, 32315474), # 0-based
  chromEnd   = c(7668550, 7673890, 55019180, 55086850, 32315600), # 1-based
  name       = c("TP53_exon1", "TP53_exon2", "EGFR_exon1", "EGFR_exon2", "BRCA2_exon1"),
  score      = c(1000, 950, 1000, 880, 1000),
  strand     = c("-", "-", "+", "+", "+"),
  stringsAsFactors = FALSE
)

# Standard BED is tab-separated with no header
write.table(bed_data, mock_bed_path, sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE)
cat(sprintf("BED file written: %s\n", mock_bed_path))

# Read BED file into R
read_bed <- function(filepath) {
  read.table(filepath, sep = "\t", header = FALSE,
             col.names = c("chrom", "chromStart", "chromEnd", "name", "score", "strand"),
             stringsAsFactors = FALSE)
}
loaded_bed <- read_bed(mock_bed_path)
loaded_bed$width <- loaded_bed$chromEnd - loaded_bed$chromStart
cat("Loaded BED Intervals with Calculated Widths (bp):\n")
print(loaded_bed[, c("chrom", "chromStart", "chromEnd", "name", "width", "strand")])

# 4. Saving & Loading R Binary Objects (.rds)
cat("\n[Step 4] Serialization via saveRDS() and readRDS()...\n")
# Saving preprocessed genomic datasets preserves factors, S4 classes, and exact metadata types
rds_path <- file.path(temp_dir, "genomic_bundle.rds")
bundle <- list(
  reference_fasta = parsed_fasta,
  annotation_bed  = loaded_bed,
  build_timestamp = Sys.time(),
  genome_assembly = "GRCh38.p14"
)
saveRDS(bundle, rds_path, compress = "xz")
file_size <- file.info(rds_path)$size
cat(sprintf("Saved compressed bundle to RDS (%d bytes): %s\n", file_size, rds_path))

# Reload and verify
restored_bundle <- readRDS(rds_path)
cat(sprintf("Successfully restored bundle! Assembly: %s, Features: %d\n",
            restored_bundle$genome_assembly, nrow(restored_bundle$annotation_bed)))

cat("\n[SUCCESS] Chapter 5 biological file I/O script executed successfully.\n")