#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 12: R for Next Gen Biologist
# Script: r12_single_cell_seurat_pipeline.R
# Description: Introduction to Single-Cell RNA-Seq Analysis with Seurat:
#              Cell filtering, mitochondrial QC, log-normalization,
#              PCA dimensionality reduction, and cluster biomarker identification.
# Usage: Rscript r12_single_cell_seurat_pipeline.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 12\n")
cat(" Single-Cell RNA-Seq Analysis: Quality Control & Clustering\n")
cat("=======================================================\n\n")

has_seurat <- requireNamespace("Seurat", quietly = TRUE)

# 1. Simulating 10x Genomics Single-Cell UMI Count Matrix
cat("[Step 1] Simulating 10x Genomics Cell-by-Gene Sparse Matrix...\n")
set.seed(42)
n_cells <- 300
n_genes <- 60

# 3 distinct immunological cell populations
# Cluster 1: T Cells (CD3D, CD3E, CD8A, IL7R)
# Cluster 2: B Cells (CD19, MS4A1, CD79A, HLA-DRA)
# Cluster 3: Monocytes (CD14, LYZ, FCGR3A, CST3)
cell_types <- c(rep("T_Cell", 120), rep("B_Cell", 100), rep("Monocyte", 80))
cell_barcodes <- paste0("AAACCCAAG", sprintf("%04d", 1:n_cells), "-1")

genes <- c(
  "CD3D", "CD3E", "CD8A", "IL7R",
  "CD19", "MS4A1", "CD79A", "HLA-DRA",
  "CD14", "LYZ", "FCGR3A", "CST3",
  paste0("MT-", c("CO1", "CO2", "ATP6", "ND1")), # Mitochondrial genes for QC
  paste0("GENE_", 1:(n_genes - 16))
)

umi_matrix <- matrix(rpois(n_genes * n_cells, lambda = 0.5), nrow = n_genes, ncol = n_cells,
                     dimnames = list(genes, cell_barcodes))

# Inject cell-type specific marker expression
umi_matrix[1:4, cell_types == "T_Cell"]   <- umi_matrix[1:4, cell_types == "T_Cell"] + rpois(4 * 120, lambda = 12)
umi_matrix[5:8, cell_types == "B_Cell"]   <- umi_matrix[5:8, cell_types == "B_Cell"] + rpois(4 * 100, lambda = 14)
umi_matrix[9:12, cell_types == "Monocyte"] <- umi_matrix[9:12, cell_types == "Monocyte"] + rpois(4 * 80, lambda = 10)

# Inject elevated mitochondrial reads in a few apoptotic cells
dead_cells <- c(10, 11, 140, 250)
umi_matrix[13:16, dead_cells] <- umi_matrix[13:16, dead_cells] + rpois(4 * length(dead_cells), lambda = 45)

cat(sprintf("Simulated Single-Cell Matrix: %d Genes across %d Single Cells\n", n_genes, n_cells))

# 2. Single-Cell Quality Control Metrics
cat("\n[Step 2] Computing Single-Cell QC Metrics (nUMI, nGene, %% Mitochondrial)...\n")
nUMI <- colSums(umi_matrix)
nGene <- colSums(umi_matrix > 0)

mt_genes <- grep("^MT-", rownames(umi_matrix), value = TRUE)
mt_counts <- colSums(umi_matrix[mt_genes, , drop = FALSE])
percent_mt <- (mt_counts / nUMI) * 100

cell_qc_df <- data.frame(
  Barcode = cell_barcodes,
  True_Type = cell_types,
  nUMI = nUMI,
  nGene = nGene,
  percent_mt = round(percent_mt, 2)
)

# Filtering thresholds: min 50 UMIs, max 15% mitochondrial reads
is_pass_qc <- (nUMI >= 50) & (percent_mt <= 15.0)
cat(sprintf("Cells Passing Clinical QC: %d / %d (%.1f%% Pass Rate)\n",
            sum(is_pass_qc), n_cells, (sum(is_pass_qc) / n_cells) * 100))
cat(sprintf("Flagged Low-Quality / Apoptotic Cells: %d cells removed\n", sum(!is_pass_qc)))

filtered_matrix <- umi_matrix[, is_pass_qc]

# 3. Log-Normalization (Total Count Scaling by 10,000)
cat("\n[Step 3] Log-Normalization: log1p(UMI / Total * 10,000)...\n")
total_counts <- colSums(filtered_matrix)
norm_matrix <- log1p(sweep(filtered_matrix, 2, total_counts, "/") * 10000)

# 4. Dimensionality Reduction & Biomarker Detection
cat("\n[Step 4] PCA Dimensionality Reduction & Top Cluster Biomarkers...\n")
pca_fit <- prcomp(t(norm_matrix), scale. = TRUE)
cat(sprintf("Top 2 PCs Account for %.1f%% of Single-Cell Transcriptomic Variation\n",
            sum(pca_fit$sdev[1:2]^2) / sum(pca_fit$sdev^2) * 100))

# Identify top differentially expressed biomarkers for T Cells vs B Cells
t_cell_cols <- cell_types[is_pass_qc] == "T_Cell"
b_cell_cols <- cell_types[is_pass_qc] == "B_Cell"

t_means <- rowMeans(norm_matrix[, t_cell_cols])
b_means <- rowMeans(norm_matrix[, b_cell_cols])
log2_fc_markers <- log2((t_means + 1e-3) / (b_means + 1e-3))

top_t_markers <- sort(log2_fc_markers, decreasing = TRUE)[1:4]
top_b_markers <- sort(log2_fc_markers, decreasing = FALSE)[1:4]

cat("\nTop Validated Biomarkers for T-Cells (vs B-Cells):\n")
for (gene in names(top_t_markers)) {
  cat(sprintf("  Marker: %-8s | Log2FC: %+5.2f\n", gene, top_t_markers[gene]))
}

cat("\nTop Validated Biomarkers for B-Cells (vs T-Cells):\n")
for (gene in names(top_b_markers)) {
  cat(sprintf("  Marker: %-8s | Log2FC: %+5.2f\n", gene, top_b_markers[gene]))
}

cat("\n[SUCCESS] Chapter 12 single-cell RNA-seq script executed successfully.\n")