#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 2: R for Next Gen Biologist
# Script: r02_data_structures_factors.R
# Description: R Data Structures for Biology:
#              Factors in clinical designs, expression count matrices,
#              heterogeneous biological lists, and clinical metadata data.frames.
# Usage: Rscript r02_data_structures_factors.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 2\n")
cat(" R Data Structures & Experimental Design Factors\n")
cat("=======================================================\n\n")

# 1. Atomic Vectors in Genomics
cat("[Step 1] Atomic Vectors in Genomics...\n")
gene_symbols <- c("TP53", "EGFR", "BRCA1", "MYC", "KRAS", "PTEN")
log2_fc <- c(2.45, -1.82, 3.12, -0.45, 1.95, -2.10)
p_values <- c(1.2e-6, 4.5e-4, 8.9e-9, 0.42, 3.1e-5, 2.8e-5)
is_significant <- (p_values < 0.001) & (abs(log2_fc) >= 1.5)

cat(sprintf("Total Genes Examined: %d\n", length(gene_symbols)))
cat(sprintf("Significant Genes   : %d\n", sum(is_significant)))
for (i in seq_along(gene_symbols)) {
  status <- ifelse(is_significant[i], "[SIGNIFICANT]", "[NS]")
  cat(sprintf("  %-6s: log2FC=%+5.2f | p-val=%.2e %s\n", 
              gene_symbols[i], log2_fc[i], p_values[i], status))
}

# 2. Factors: Crucial for Experimental Design in RNA-seq
cat("\n[Step 2] Biological Factors in Experimental Design...\n")
# Raw character vector of sample conditions
raw_conditions <- c("Control", "Treated", "Control", "Treated", "Control", "Treated")
# In DESeq2, the baseline condition MUST be the first level of the factor
condition_factor <- factor(raw_conditions, levels = c("Control", "Treated"))

tissue_type <- factor(c("Liver", "Liver", "Kidney", "Kidney", "Lung", "Lung"))
tumor_stage <- factor(c("Stage_I", "Stage_IV", "Stage_II", "Stage_III", "Stage_I", "Stage_IV"),
                      levels = c("Stage_I", "Stage_II", "Stage_III", "Stage_IV"),
                      ordered = TRUE)

cat("Sample Condition Factor Levels (Baseline is first):\n")
print(levels(condition_factor))
cat("Tumor Stage (Ordered Factor):\n")
print(table(tumor_stage))

# 3. Gene Expression Count Matrix
cat("\n[Step 3] Constructing Gene Expression Count Matrix...\n")
set.seed(42)
num_genes <- 5
num_samples <- 6
sample_names <- paste0("Sample_", 1:num_samples)
genes <- c("TP53", "EGFR", "BRCA1", "MYC", "KRAS")

# Simulate raw RNA-seq read counts
counts_matrix <- matrix(
  rpois(num_genes * num_samples, lambda = c(120, 850, 45, 3200, 410)),
  nrow = num_genes,
  ncol = num_samples,
  dimnames = list(genes, sample_names)
)

cat("Raw Count Matrix (Genes x Samples):\n")
print(counts_matrix)

# Computing basic row-wise and column-wise statistics
cat("\nLibrary Sizes (Total Reads per Sample):\n")
print(colSums(counts_matrix))

cat("\nMean Read Count per Gene Across Cohort:\n")
print(round(rowMeans(counts_matrix), 1))

# 4. Data Frames: Integrating Clinical Metadata & Experimental Covariates
cat("\n[Step 4] Clinical & Phenotypic Metadata DataFrame...\n")
col_data <- data.frame(
  SampleID    = sample_names,
  Condition   = condition_factor,
  Tissue      = tissue_type,
  Stage       = tumor_stage,
  Batch       = factor(c("Batch_1", "Batch_1", "Batch_2", "Batch_2", "Batch_1", "Batch_2")),
  Age         = c(54, 62, 47, 71, 58, 65),
  stringsAsFactors = FALSE
)
rownames(col_data) <- col_data$SampleID

print(col_data)

# 5. Lists: Heterogeneous Biological Data Container
cat("\n[Step 5] Lists for Storing Multi-Omics Pipeline Results...\n")
experiment_container <- list(
  experiment_name = "Lung Adenocarcinoma Transcriptomic Screen",
  date_conducted  = Sys.Date(),
  platform        = "Illumina NovaSeq 6000 Paired-End 150bp",
  raw_counts      = counts_matrix,
  sample_metadata = col_data,
  qc_metrics      = list(
    mean_q30_pct  = 94.6,
    adapter_dimers = 0.2,
    mapped_pct    = 91.8
  )
)

cat(sprintf("Experiment Container Created: '%s'\n", experiment_container$experiment_name))
cat(sprintf("Platform Used               : %s\n", experiment_container$platform))
cat(sprintf("Mean Q30 Quality            : %.1f%%\n", experiment_container$qc_metrics$mean_q30_pct))

cat("\n[SUCCESS] Chapter 2 data structures script executed successfully.\n")