# Cell Omics Guide / Shengxin Manmanxue: corrected teaching adaptation.
# Not validated on the full dataset; not an exact Figure 1 reproduction.
# Run each stage from a new working directory; review outputs before continuing.
# Packages: Seurat 5, Matrix, data.table, R.utils, sctransform, future,
# ggplot2, patchwork, SingleCellExperiment, scDblFinder; optional harmony.
args <- commandArgs(trailingOnly = TRUE)
stage <- if (length(args)) args[[1]] else "help"
allowed <- c("import", "qc-audit", "clean", "cluster", "harmony", "compose")
if (!stage %in% allowed) stop(paste("Choose a stage:", paste(allowed, collapse = ", ")))
suppressPackageStartupMessages({ library(Seurat); library(Matrix); library(ggplot2); library(patchwork) })
stopifnot(packageVersion("Seurat") >= "5.0.0")
future::plan(future::sequential)
set.seed(1224)
out <- "lung-output"
dir.create(out, showWarnings = FALSE)
p <- function(name) file.path(out, name)
put <- function(x, name) write.csv(x, p(name), row.names = FALSE)
load_object <- function(name) {
  stopifnot(file.exists(p(name)))
  readRDS(p(name))
}
save_plot <- function(x, name, width = 14, height = 8) ggsave(p(name), x, width = width, height = height, limitsize = FALSE)
normal <- paste0("LUNG_N", c("01", "06", "08", "09", "18", "19", "20", "28", "30", "31", "34"))
tumor <- paste0("LUNG_T", c("06", "08", "09", "18", "19", "20", "25", "28", "30", "31", "34"))
target <- c(normal, tumor)
qc_features <- c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rbc")
add_qc <- function(x) {
  mt <- grep("^MT-", rownames(x), value = TRUE)
  hb <- intersect(c("HBA1", "HBA2", "HBB", "HBD", "HBE1", "HBG1", "HBG2", "HBM", "HBQ1", "HBZ"), rownames(x))
  if (!length(mt) || !length(hb)) stop("Missing MT-/hemoglobin symbols: audit feature naming before applying percentage thresholds.")
  x[["percent.mt"]] <- PercentageFeatureSet(x, features = mt)
  x[["percent.rbc"]] <- PercentageFeatureSet(x, features = hb)
  writeLines(c("Mitochondrial genes:", mt, "Hemoglobin genes:", hb), p("qc-gene-sets.txt"))
  x
}
read_design <- function() {
  d <- read.csv(p("sample-design.csv"), colClasses = "character", check.names = FALSE)
  stopifnot(all(c("Sample", "group", "capture_id", "capture_evidence", "donor_id") %in% names(d)),
    !anyDuplicated(d$Sample), setequal(d$Sample, target))
  d <- d[match(target, d$Sample), , drop = FALSE]
  stopifnot(identical(d$group, ifelse(d$Sample %in% normal, "Normal", "Tumor")))
  if (anyNA(d$capture_id) || anyNA(d$capture_evidence) || any(!nzchar(trimws(d$capture_id))) || any(!nzchar(trimws(d$capture_evidence)))) {
    stop("Fill verified capture_id and capture_evidence in sample-design.csv. If a sample spans captures, first add a cell-level capture map; this teaching script assumes one capture per Sample.")
  }
  d
}
check_samples <- function(x) {
  missing <- setdiff(target, unique(x$Sample))
  if (length(missing)) stop(paste("Samples with no retained cells; review QC/design before proceeding:", paste(missing, collapse = ", ")))
}

if (stage == "import") {
  paths <- c("GSE131907_Lung_Cancer_raw_UMI_matrix.rds.gz", "GSE131907_Lung_Cancer_raw_UMI_matrix.rds")
  counts_file <- paths[file.exists(paths)][1]
  if (is.na(counts_file)) stop("Download the raw UMI RDS from GEO first.")
  anno_file <- "GSE131907_Lung_Cancer_cell_annotation.txt.gz"
  stopifnot(file.exists(anno_file))
  x <- readRDS(counts_file) # readRDS supports gzip-compressed RDS files.
  if (!inherits(x, "sparseMatrix")) stop("Expected a sparse count matrix. Inspect the downloaded RDS instead of densifying it.")
  stopifnot(!is.null(rownames(x)), !is.null(colnames(x)), !anyDuplicated(rownames(x)), !anyDuplicated(colnames(x)))
  stopifnot(all(is.finite(x@x)), all(x@x >= 0), all(x@x == floor(x@x)))
  anno <- data.table::fread(anno_file, data.table = FALSE)
  stopifnot(all(c("Index", "Sample") %in% names(anno)), !anyNA(anno$Index), !anyDuplicated(anno$Index))
  selected <- anno[anno$Sample %in% target, , drop = FALSE]
  stopifnot(setequal(unique(selected$Sample), target), all(selected$Index %in% colnames(x)))
  cells <- colnames(x)[colnames(x) %in% selected$Index]
  selected <- selected[match(cells, selected$Index), , drop = FALSE]
  rownames(selected) <- selected$Index
  stopifnot(identical(rownames(selected), cells))
  selected$group <- ifelse(selected$Sample %in% normal, "Normal", "Tumor")
  lung <- CreateSeuratObject(x[, cells, drop = FALSE], meta.data = selected, project = "GSE131907")
  stopifnot(identical(rownames(lung[[]]), colnames(lung)))
  saveRDS(lung, p("lung-raw.rds"))
  put(as.data.frame(table(lung$Sample, lung$group)), "cohort-counts.csv")
  if (!file.exists(p("sample-design.csv"))) put(data.frame(Sample = target,
    group = ifelse(target %in% normal, "Normal", "Tumor"), capture_id = "", capture_evidence = "", donor_id = ""), "sample-design.csv")
  capture.output(list(source_dimensions = dim(x), subset_dimensions = dim(lung),
    checksums = tools::md5sum(c(counts_file, anno_file)), metadata_aligned = identical(rownames(lung[[]]), colnames(lung))), file = p("input-audit.txt"))
}
if (stage == "qc-audit") {
  lung <- add_qc(load_object("lung-raw.rds"))
  put(data.frame(barcode = colnames(lung), lung[[]][, c("Sample", qc_features)]), "qc-metrics.csv")
  save_plot(VlnPlot(lung, features = qc_features, group.by = "Sample", pt.size = 0, ncol = 2), "QC-before.pdf", height = 10)
  if (!file.exists(p("qc-settings.csv"))) put(data.frame(Sample = target, min_genes = 200, max_genes = 6000,
    min_umi = 200, max_umi = 40000, max_mt = 20, max_hb = 5, reason = ""), "qc-settings.csv")
}
if (stage == "clean") {
  d <- read_design()
  settings <- read.csv(p("qc-settings.csv"), check.names = FALSE)
  numeric_fields <- c("min_genes", "max_genes", "min_umi", "max_umi", "max_mt", "max_hb")
  stopifnot(all(c("Sample", "reason", numeric_fields) %in% names(settings)), !anyDuplicated(settings$Sample), setequal(settings$Sample, target))
  if (anyNA(settings) || any(!nzchar(trimws(settings$reason)))) stop("Review and document every sample QC threshold first.")
  stopifnot(all(vapply(settings[numeric_fields], is.numeric, logical(1))), all(is.finite(as.matrix(settings[numeric_fields]))),
    all(settings$min_genes >= 0), all(settings$min_umi >= 0), all(settings$min_genes < settings$max_genes),
    all(settings$min_umi < settings$max_umi), all(settings$max_mt > 0 & settings$max_mt <= 100), all(settings$max_hb > 0 & settings$max_hb <= 100))
  lung <- add_qc(load_object("lung-raw.rds"))
  lung$capture_id <- d$capture_id[match(lung$Sample, d$Sample)]
  thresholds <- settings[match(lung$Sample, settings$Sample), , drop = FALSE]
  keep <- with(lung[[]], nFeature_RNA > thresholds$min_genes & nFeature_RNA < thresholds$max_genes &
    nCount_RNA > thresholds$min_umi & nCount_RNA < thresholds$max_umi & percent.mt < thresholds$max_mt & percent.rbc < thresholds$max_hb)
  decisions <- data.frame(barcode = colnames(lung), Sample = lung$Sample, keep = keep)
  put(decisions, "qc-decisions.csv")
  put(settings, "qc-applied-settings.csv")
  put(d, "sample-design-applied.csv")
  raw_n <- table(factor(lung$Sample, levels = target))
  lung <- subset(lung, cells = colnames(lung)[keep])
  check_samples(lung)
  qc_n <- table(factor(lung$Sample, levels = target))
  sce <- SingleCellExperiment::SingleCellExperiment(assays = list(counts = GetAssayData(lung, assay = "RNA", layer = "counts")))
  sce <- scDblFinder::scDblFinder(sce, samples = lung$capture_id)
  meta <- as.data.frame(SummarizedExperiment::colData(sce))[, c("scDblFinder.class", "scDblFinder.score"), drop = FALSE]
  stopifnot(setequal(rownames(meta), colnames(lung)))
  lung <- AddMetaData(lung, metadata = meta)
  put(data.frame(barcode = colnames(lung), Sample = lung$Sample, capture_id = lung$capture_id,
    class = lung$scDblFinder.class, score = lung$scDblFinder.score), "doublet-decisions.csv")
  lung <- subset(lung, subset = scDblFinder.class == "singlet")
  check_samples(lung)
  put(data.frame(Sample = target, raw = as.integer(raw_n), after_qc = as.integer(qc_n),
    singlets = as.integer(table(factor(lung$Sample, levels = target)))), "retention.csv")
  saveRDS(lung, p("lung-clean.rds"))
  save_plot(VlnPlot(lung, features = qc_features, group.by = "Sample", pt.size = 0, ncol = 2), "QC-after.pdf", height = 10)
}
if (stage == "cluster") {
  if (file.exists(p("annotation-review.csv"))) stop("Archive annotation-review.csv before rerunning cluster. New clustering requires a newly generated and reviewed annotation table, even when numeric cluster IDs look unchanged.")
  lung <- load_object("lung-clean.rds")
  if (ncol(lung) < 60) stop("Too few cells for this teaching configuration.")
  lung <- SCTransform(lung, assay = "RNA", conserve.memory = TRUE, verbose = FALSE)
  lung <- RunPCA(lung, npcs = 30, verbose = FALSE)
  lung <- FindNeighbors(lung, reduction = "pca", dims = 1:20)
  lung <- FindClusters(lung, resolution = 0.3, random.seed = 1224)
  lung <- RunUMAP(lung, reduction = "pca", dims = 1:20, reduction.name = "umap.pca", reduction.key = "uPCA_", seed.use = 1224)
  save_plot(ElbowPlot(lung, ndims = 30) | VizDimLoadings(lung, dims = 1:2, reduction = "pca"), "pca-diagnostics.pdf")
  save_plot(DimPlot(lung, reduction = "umap.pca", group.by = c("seurat_clusters", "Sample", "group"), label = FALSE), "embedding.pdf")
  markers <- FindAllMarkers(lung, assay = "SCT", only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
  put(markers, "markers.csv")
  # These are evidence prompts, not a fixed dictionary of cluster labels.
  panel <- list(T = c("CD3D", "CD3E"), NK = c("NKG7", "GNLY"), Myeloid = c("LYZ", "CD68", "S100A9"),
    Epithelial = c("EPCAM", "KRT18", "KRT19"), B = c("MS4A1", "CD79A"), Mast = c("TPSAB1", "KIT"),
    Fibroblast = c("COL1A1", "DCN"), Endothelial = c("PECAM1", "VWF"), Plasma = c("JCHAIN", "MZB1"))
  panel <- lapply(panel, intersect, y = rownames(lung))
  panel <- panel[lengths(panel) > 0]
  stopifnot(length(panel) > 0)
  lung <- NormalizeData(lung, assay = "RNA", verbose = FALSE)
  save_plot(DotPlot(lung, assay = "RNA", features = panel) + RotatedAxis(), "marker-evidence.pdf")
  saveRDS(lung, p("lung-clusters.rds"))
  if (!file.exists(p("annotation-review.csv"))) put(data.frame(cluster = levels(Idents(lung)), celltype = "uncertain", evidence = "", reviewed = FALSE), "annotation-review.csv")
}
if (stage == "harmony") {
  stopifnot(requireNamespace("harmony", quietly = TRUE))
  lung <- load_object("lung-clusters.rds")
  lung <- harmony::RunHarmony(lung, group.by.vars = "Sample", reduction.use = "pca", dims.use = 1:20)
  lung <- RunUMAP(lung, reduction = "harmony", dims = 1:20, reduction.name = "umap.harmony", reduction.key = "uHarmony_", seed.use = 1224)
  saveRDS(lung, p("lung-harmony-comparison.rds"))
  save_plot((DimPlot(lung, reduction = "umap.pca", group.by = "Sample") | DimPlot(lung, reduction = "umap.harmony", group.by = "Sample")) /
    (DimPlot(lung, reduction = "umap.pca", group.by = "group") | DimPlot(lung, reduction = "umap.harmony", group.by = "group")), "harmony-comparison.pdf", height = 12)
}
if (stage == "compose") {
  lung <- load_object("lung-clusters.rds")
  a <- read.csv(p("annotation-review.csv"), colClasses = "character", check.names = FALSE)
  stopifnot(all(c("cluster", "celltype", "evidence", "reviewed") %in% names(a)), !anyDuplicated(a$cluster), setequal(a$cluster, levels(Idents(lung))))
  if (anyNA(a) || any(!nzchar(trimws(a$celltype))) || any(!nzchar(trimws(a$evidence))) || any(toupper(trimws(a$reviewed)) != "TRUE")) stop("Review every cluster; uncertain is an acceptable label with documented evidence.")
  lung$celltype <- a$celltype[match(as.character(Idents(lung)), a$cluster)]
  check_samples(lung)
  counts <- table(Sample = factor(lung$Sample, levels = target), celltype = factor(lung$celltype, levels = unique(a$celltype)))
  proportions <- prop.table(counts, margin = 1)
  stopifnot(all(abs(rowSums(proportions) - 1) < 1e-8))
  result <- as.data.frame(counts, responseName = "n")
  result$proportion <- as.vector(proportions)
  result$group <- ifelse(result$Sample %in% normal, "Normal", "Tumor")
  put(result, "sample-proportions.csv")
  put(a, "annotation-decisions.csv")
  saveRDS(lung, p("lung-annotated.rds"))
  save_plot(DimPlot(lung, reduction = "umap.pca", group.by = "celltype", label = TRUE, repel = TRUE), "annotated-UMAP.pdf")
  save_plot(ggplot(result, aes(x = Sample, y = proportion, fill = celltype)) + geom_col() +
    facet_grid(. ~ group, scales = "free_x", space = "free_x") + theme_bw() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)), "sample-composition.pdf")
  # Deliberately no condition test: donor IDs and the inference design are unverified.
}
capture.output(sessionInfo(), file = p(paste0("session-", stage, ".txt")))
message("Stage completed. Review artifacts in ", out, " before marking the stage verified.")