#!/usr/bin/env Rscript

# Registered single-cell extension for FSHD GSE303359 and SMA GSE290980.
# Inference is restricted to donor/cell-line units. Cells and technical
# libraries are never treated as biological replicates.

options(stringsAsFactors = FALSE, future.globals.maxSize = 30 * 1024^3)
set.seed(20260805)

suppressPackageStartupMessages({
  library(Matrix)
  library(Seurat)
  library(edgeR)
  library(limma)
  library(ggplot2)
  library(patchwork)
  library(jsonlite)
  library(SingleCellExperiment)
  library(scuttle)
  library(scran)
})

ARGS <- commandArgs(trailingOnly = TRUE)
SELF_TEST <- "--self-test" %in% ARGS
FINALIZE_CANDIDATE_TABLES <- "--finalize-candidate-tables-only" %in% ARGS
ROOT_ARG <- grep("^--root=", ARGS, value = TRUE)
ROOT <- normalizePath(if (length(ROOT_ARG)) sub("^--root=", "", ROOT_ARG[[1]]) else getwd(), mustWork = TRUE)
DATA <- file.path(ROOT, "data", "external", "multidisease_v02")
OUT <- file.path(ROOT, "results", "multidisease_release_v02")
TABLES <- file.path(OUT, "tables")
FIGURES <- file.path(OUT, "figures")
LOGS <- file.path(OUT, "logs")

COLORS <- c(
  navy = "#18324A", blue = "#2A6F97", teal = "#278479", amber = "#D39C2E",
  purple = "#886173", red = "#C5534C", grey = "#7B8993", lightgrey = "#D5DFE3"
)

FSHD_MODULES <- list(
  DUX4_target = c("CCNA1", "DUX4", "DUXA", "KHDC1L", "LEUTX", "MBD3L2", "MBD3L3", "MBD3L5", "PRAMEF1", "PRAMEF2", "PRAMEF4", "RFPL2", "SLC34A2", "TRIM43", "TRIM49", "TRIM51", "ZSCAN4"),
  Membrane_repair = c("DYSF", "ANXA1", "ANXA2", "ANXA5", "ANXA6", "CAPN3", "TRIM72", "CAV3", "LAMP1", "CHMP3", "CHMP4B"),
  Ferroptosis = c("GPX4", "SLC7A11", "ACSL4", "TFRC", "FTH1", "FTL", "NCOA4", "HMOX1", "SAT1", "LPCAT3", "ALOX15"),
  Myogenesis = c("PAX7", "MYF5", "MYOD1", "MYOG", "DES", "ACTA1", "MYH3", "MYH8"),
  Stress_response = c("HSPA1A", "HSPA1B", "JUN", "FOS", "ATF3", "DDIT3", "HMOX1", "DNAJB1")
)

FSHD_STATE_MARKERS <- list(
  Cycling_myoblast = c("MKI67", "TOP2A", "CENPF", "TYMS", "UBE2C"),
  Basal_myoblast = c("PAX7", "MYF5", "VCAM1", "CDH15", "ITGA7"),
  Myogenic_differentiation = c("MYOD1", "MYOG", "DES", "ACTA1", "MYH3", "MYH8"),
  DUX4_target_high = FSHD_MODULES$DUX4_target,
  Stress_injury = c("HSPA1A", "HSPA1B", "JUN", "FOS", "ATF3", "DDIT3")
)

SMA_MODULES <- list(
  SMN_axis = c("SMN1", "SMN2", "GEMIN2", "GEMIN3", "GEMIN4", "GEMIN5", "GEMIN6", "GEMIN7", "GEMIN8", "STRAP"),
  Motor_neuron = c("MNX1", "ISL1", "ISL2", "CHAT", "SLC18A3", "PRPH", "NEFH", "NEFM"),
  Neuronal_maturation = c("RBFOX3", "MAP2", "TUBB3", "SNAP25", "SYT1", "NEFL", "DCX"),
  Synaptic = c("SYN1", "SYP", "SNAP25", "STX1A", "VAMP2", "DLG4", "SYT1"),
  Axonal_transport = c("KIF5A", "KIF5B", "KIF5C", "DYNC1H1", "DCTN1", "BICD2", "KLC1", "KLC2", "STMN2"),
  OXPHOS = c("NDUFA1", "NDUFA2", "NDUFA4", "NDUFB3", "NDUFS1", "SDHA", "SDHB", "UQCRC1", "UQCRC2", "COX4I1", "COX5A", "ATP5F1A", "ATP5F1B"),
  Cell_cycle = c("MKI67", "TOP2A", "TYMS", "UBE2C", "CENPF"),
  Stress_response = c("HSPA1A", "HSPA1B", "JUN", "FOS", "ATF3", "DDIT3")
)

SMA_STATE_MARKERS <- list(
  Cycling_progenitor = c("MKI67", "TOP2A", "TYMS", "UBE2C", "CENPF"),
  Neural_progenitor = c("SOX2", "NES", "VIM", "HES1", "HES5", "SOX1"),
  Motor_neuron_like = c("MNX1", "ISL1", "ISL2", "CHAT", "SLC18A3", "PRPH"),
  Neuron = c("RBFOX3", "MAP2", "TUBB3", "SNAP25", "SYT1"),
  Inhibitory_neuron = c("GAD1", "GAD2", "SLC32A1"),
  Excitatory_neuron = c("SLC17A6", "SLC17A7", "TBR1", "NEUROD2"),
  Glial_neural_crest = c("SOX10", "S100B", "PLP1", "SLC1A3", "FABP7"),
  Stress_state = c("HSPA1A", "HSPA1B", "JUN", "FOS", "ATF3", "DDIT3")
)

theme_nmd <- function(base_size = 8) {
  theme_minimal(base_size = base_size) +
    theme(
      plot.title = element_text(face = "bold", colour = COLORS[["navy"]], size = base_size + 1.2),
      plot.subtitle = element_text(colour = COLORS[["grey"]], size = base_size - 0.4),
      plot.caption = element_text(colour = COLORS[["grey"]], size = base_size - 1),
      axis.title = element_text(colour = COLORS[["navy"]]),
      axis.text = element_text(colour = "#31404B"),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(colour = "#E1E8EB", linewidth = 0.3),
      legend.title = element_text(face = "bold"),
      strip.text = element_text(face = "bold", colour = COLORS[["navy"]]),
      strip.background = element_rect(fill = "#F2F6F7", colour = NA)
    )
}

sanitize_tsv_text <- function(values) {
  gsub("[\\t\\r\\n]+", " ", as.character(values), perl = TRUE)
}

write_tsv <- function(x, path) {
  dir.create(dirname(path), recursive = TRUE, showWarnings = FALSE)
  connection <- if (grepl("\\.gz$", path)) gzfile(path, "wt") else file(path, "wt")
  on.exit(close(connection), add = TRUE)
  output <- x
  text_columns <- vapply(output, function(column) is.character(column) || is.factor(column), logical(1))
  output[text_columns] <- lapply(output[text_columns], sanitize_tsv_text)
  write.table(output, connection, sep = "\t", quote = FALSE, row.names = FALSE, na = "")
}

save_plot_family <- function(plot, stem, width, height) {
  dir.create(dirname(stem), recursive = TRUE, showWarnings = FALSE)
  ggsave(paste0(stem, ".pdf"), plot, width = width, height = height, units = "in", device = cairo_pdf)
  ggsave(paste0(stem, ".png"), plot, width = width, height = height, units = "in", dpi = 300, bg = "white")
  svg(paste0(stem, ".svg"), width = width, height = height, family = "sans")
  print(plot)
  dev.off()
}

pairwise_consistency <- function(case_values, control_values) {
  differences <- as.vector(outer(case_values, control_values, "-"))
  if (all(differences > 0)) return("UP_ALL_PAIRWISE")
  if (all(differences < 0)) return("DOWN_ALL_PAIRWISE")
  "MIXED"
}

pairwise_effect_table <- function(log_values, metadata, case_label, control_label, state = NULL) {
  case_columns <- metadata$group[metadata$genotype == case_label]
  control_columns <- metadata$group[metadata$genotype == control_label]
  pairs <- expand.grid(case = case_columns, control = control_columns, stringsAsFactors = FALSE)
  differences <- do.call(cbind, lapply(seq_len(nrow(pairs)), function(index) {
    as.numeric(log_values[, pairs$case[index]] - log_values[, pairs$control[index]])
  }))
  differences <- matrix(differences, nrow = nrow(log_values),
                         dimnames = list(rownames(log_values), paste(pairs$case, pairs$control, sep = "_minus_")))
  result <- data.frame(
    feature = rownames(log_values),
    log2_effect = rowMeans(log_values[, case_columns, drop = FALSE]) - rowMeans(log_values[, control_columns, drop = FALSE]),
    pairwise_min = apply(differences, 1, min),
    pairwise_max = apply(differences, 1, max),
    directional_consistency = apply(differences, 1, function(x) if (all(x > 0)) "UP_ALL_PAIRWISE" else if (all(x < 0)) "DOWN_ALL_PAIRWISE" else "MIXED"),
    p_value = NA_real_, fdr = NA_real_, inference = "directional_n2_lines_per_genotype_no_calibrated_p_value",
    stringsAsFactors = FALSE
  )
  if (!is.null(state)) result$marker_state <- state
  result
}

get_counts <- function(object) {
  tryCatch(
    GetAssayData(object, assay = "RNA", layer = "counts"),
    error = function(e) GetAssayData(object, assay = "RNA", slot = "counts")
  )
}

aggregate_counts <- function(object, groups) {
  counts <- get_counts(object)
  groups <- factor(as.character(groups), levels = unique(as.character(groups)))
  design <- sparse.model.matrix(~ 0 + groups)
  colnames(design) <- levels(groups)
  result <- counts %*% design
  colnames(result) <- levels(groups)
  result
}

complete_candidate_inventory <- function(effect_table, candidates, contrasts = NULL, effect_column) {
  grid <- if (is.null(contrasts)) {
    data.frame(gene = candidates, stringsAsFactors = FALSE)
  } else {
    expand.grid(gene = candidates, contrast = contrasts, stringsAsFactors = FALSE)
  }
  grid$.registered_order <- seq_len(nrow(grid))
  keys <- intersect(c("gene", "contrast"), names(grid))
  result <- merge(grid, effect_table, by = keys, all.x = TRUE, sort = FALSE)
  result <- result[order(result$.registered_order), , drop = FALSE]
  result$.registered_order <- NULL
  result$candidate_effect_status <- ifelse(
    is.na(result[[effect_column]]), "BELOW_EXPRESSION_FILTER", "ESTIMATED"
  )
  result$registered_before_analysis <- TRUE
  rownames(result) <- NULL
  result
}

center_nested_batch_pcs <- function(embedding, biological_line, technical_library) {
  if (nrow(embedding) != length(biological_line) || nrow(embedding) != length(technical_library)) {
    stop("Nested-batch PCA centering metadata length mismatch")
  }
  corrected <- embedding
  for (line in unique(biological_line)) {
    line_indices <- which(biological_line == line)
    line_centroid <- colMeans(embedding[line_indices, , drop = FALSE])
    for (library in unique(technical_library[line_indices])) {
      library_indices <- which(biological_line == line & technical_library == library)
      library_centroid <- colMeans(embedding[library_indices, , drop = FALSE])
      corrected[library_indices, ] <- sweep(
        embedding[library_indices, , drop = FALSE], 2, library_centroid - line_centroid, "-"
      )
    }
  }
  corrected
}

canonical_reduction_frame <- function(coordinates, prefix) {
  frame <- as.data.frame(coordinates, check.names = FALSE)
  names(frame) <- paste0(prefix, "_", seq_len(ncol(frame)))
  frame
}

canonicalize_cluster_columns <- function(average_expression, object_cluster_ids) {
  observed <- unique(as.character(object_cluster_ids))
  encoded <- colnames(average_expression)
  if (setequal(encoded, observed)) return(average_expression)
  if (setequal(encoded, paste0("g", observed))) {
    colnames(average_expression) <- sub("^g", "", encoded)
    return(average_expression)
  }
  stop(sprintf(
    "Average-expression cluster columns do not map to object clusters: average=%s; object=%s",
    paste(encoded, collapse = ","), paste(observed, collapse = ",")
  ))
}

run_self_test <- function() {
  stopifnot(pairwise_consistency(c(2, 3), c(-1, 0)) == "UP_ALL_PAIRWISE")
  stopifnot(pairwise_consistency(c(-2, -3), c(0, 1)) == "DOWN_ALL_PAIRWISE")
  stopifnot(pairwise_consistency(c(-1, 2), c(0, 1)) == "MIXED")
  toy <- Matrix(matrix(c(1, 0, 2, 0, 3, 0, 0, 4, 0, 5, 0, 6), nrow = 3), sparse = TRUE)
  rownames(toy) <- paste0("g", 1:3)
  colnames(toy) <- paste0("c", 1:4)
  groups <- factor(c("A", "A", "B", "B"), levels = c("A", "B"))
  design <- sparse.model.matrix(~ 0 + groups)
  colnames(design) <- levels(groups)
  summed <- toy %*% design
  stopifnot(all(dim(summed) == c(3, 2)), summed[1, "A"] == 1, summed[1, "B"] == 5)
  toy_log <- matrix(c(1, 2, 3, 4, 4, 3, 2, 1), nrow = 2, byrow = TRUE,
                    dimnames = list(c("up", "down"), c("CTRL1", "CTRL2", "SMA1", "SMA2")))
  toy_meta <- data.frame(group = colnames(toy_log), genotype = c("CTRL", "CTRL", "SMA", "SMA"))
  toy_effect <- pairwise_effect_table(toy_log, toy_meta, "SMA", "CTRL")
  stopifnot(nrow(toy_effect) == 2,
            toy_effect$directional_consistency[1] == "UP_ALL_PAIRWISE",
            toy_effect$directional_consistency[2] == "DOWN_ALL_PAIRWISE",
            toy_effect$pairwise_min[1] == 1, toy_effect$pairwise_max[1] == 3,
            toy_effect$pairwise_min[2] == -3, toy_effect$pairwise_max[2] == -1)
  toy_dge <- edgeR::DGEList(matrix(c(10, 20, 30, 40), nrow = 2))
  toy_cpm <- edgeR::cpm(toy_dge, log = TRUE, prior.count = 2)
  stopifnot(is.matrix(toy_cpm), all(dim(toy_cpm) == c(2, 2)))
  toy_embedding <- matrix(seq_len(6), nrow = 3, dimnames = list(NULL, c("umap_1", "umap_2")))
  canonical_embedding <- canonical_reduction_frame(toy_embedding, "UMAP")
  stopifnot(identical(names(canonical_embedding), c("UMAP_1", "UMAP_2")))
  toy_average <- matrix(seq_len(6), nrow = 3, dimnames = list(paste0("gene", 1:3), c("g0", "g1")))
  canonical_average <- canonicalize_cluster_columns(toy_average, c("0", "0", "1", "1"))
  stopifnot(identical(colnames(canonical_average), c("0", "1")))
  stopifnot(inherits(
    try(canonicalize_cluster_columns(toy_average, c("A", "B")), silent = TRUE),
    "try-error"
  ))
  toy_candidates <- complete_candidate_inventory(
    data.frame(gene = "A", contrast = "case", logFC = 1),
    c("A", "B"), c("case", "interaction"), "logFC"
  )
  stopifnot(
    nrow(toy_candidates) == 4,
    sum(toy_candidates$candidate_effect_status == "ESTIMATED") == 1,
    sum(toy_candidates$candidate_effect_status == "BELOW_EXPRESSION_FILTER") == 3,
    identical(unique(toy_candidates$gene), c("A", "B"))
  )
  toy_pcs <- matrix(c(0, 0, 2, 2, 10, 10, 14, 14), ncol = 2, byrow = TRUE)
  toy_centered <- center_nested_batch_pcs(toy_pcs, c("L1", "L1", "L2", "L2"), c("A", "B", "C", "D"))
  stopifnot(
    all(colMeans(toy_centered[1:2, , drop = FALSE]) == colMeans(toy_pcs[1:2, , drop = FALSE])),
    all(colMeans(toy_centered[3:4, , drop = FALSE]) == colMeans(toy_pcs[3:4, , drop = FALSE]))
  )
  stopifnot(identical(sanitize_tsv_text("line 1\nline 2\tfield"), "line 1 line 2 field"))
  cat(jsonlite::toJSON(list(status = "PASS", tests = 25, seed = 20260805), auto_unbox = TRUE, pretty = TRUE), "\n")
}

if (SELF_TEST) {
  run_self_test()
  quit(save = "no", status = 0)
}

dir.create(TABLES, recursive = TRUE, showWarnings = FALSE)
dir.create(FIGURES, recursive = TRUE, showWarnings = FALSE)
dir.create(LOGS, recursive = TRUE, showWarnings = FALSE)

if (FINALIZE_CANDIDATE_TABLES) {
  fshd_de <- read.delim(gzfile(file.path(TABLES, "fshd_gse303359_line_blocked_pseudobulk_de.tsv.gz")), check.names = FALSE)
  fshd_candidates <- complete_candidate_inventory(
    fshd_de,
    c("DUX4", "DUXA", "ZSCAN4", "LEUTX", "MAPK8", "MAPK14"),
    c("FSHD_vs_Healthy_baseline", "Injury_in_Healthy", "Injury_in_FSHD", "Genotype_by_injury_interaction"),
    "logFC"
  )
  write_tsv(fshd_candidates, file.path(TABLES, "fshd_gse303359_registered_candidate_effects.tsv"))
  sma_effects <- read.delim(gzfile(file.path(TABLES, "sma_gse290980_line_level_directional_effects.tsv.gz")), check.names = FALSE)
  sma_candidates <- complete_candidate_inventory(sma_effects, c("SMN1", "SMN2", "KIF5A"), effect_column = "log2_effect")
  write_tsv(sma_candidates, file.path(TABLES, "sma_gse290980_registered_candidate_effects.tsv"))
  finalization <- list(
    schema = "nmd-vcell-registered-candidate-table-finalization/1.0",
    status = "PASS",
    operation = "explicit placeholder rows for prespecified candidates below the unchanged expression filters",
    fshd_expected_rows = 24,
    fshd_estimated_rows = sum(fshd_candidates$candidate_effect_status == "ESTIMATED"),
    fshd_filtered_rows = sum(fshd_candidates$candidate_effect_status == "BELOW_EXPRESSION_FILTER"),
    sma_expected_rows = 3,
    sma_estimated_rows = sum(sma_candidates$candidate_effect_status == "ESTIMATED"),
    sma_filtered_rows = sum(sma_candidates$candidate_effect_status == "BELOW_EXPRESSION_FILTER"),
    analysis_design_changes = list(), threshold_changes = list(), candidate_rule_changes = list()
  )
  analysis_plan_path <- file.path(OUT, "analysis_plan.json")
  if (file.exists(analysis_plan_path)) {
    locked_plan <- jsonlite::read_json(analysis_plan_path, simplifyVector = FALSE)
    locked_plan$candidate_output_contract <- list(
      GSE303359 = c("DUX4", "DUXA", "ZSCAN4", "LEUTX", "MAPK8", "MAPK14"),
      GSE290980 = c("SMN1", "SMN2", "KIF5A"),
      filtered_candidate_policy = "retain explicit BELOW_EXPRESSION_FILTER rows with no imputed effect"
    )
    jsonlite::write_json(locked_plan, analysis_plan_path, pretty = TRUE, auto_unbox = TRUE)
  }
  jsonlite::write_json(finalization, file.path(OUT, "candidate_table_finalization.json"), pretty = TRUE, auto_unbox = TRUE)
  cat(jsonlite::toJSON(finalization, pretty = TRUE, auto_unbox = TRUE), "\n")
  quit(save = "no", status = 0)
}

manifest_path <- file.path(DATA, "download_manifest.json")
metadata_path <- file.path(DATA, "sample_metadata.tsv")
if (!file.exists(manifest_path) || !file.exists(metadata_path)) stop("Download manifest or sample metadata is missing")
download_manifest <- jsonlite::read_json(manifest_path, simplifyVector = TRUE)
for (accession in c("GSE303359", "GSE290980")) {
  ready <- isTRUE(download_manifest$accessions[[accession]]$analysis_ready)
  if (!ready) stop(sprintf("%s is not analysis-ready according to the file manifest", accession))
}

sample_metadata <- read.delim(metadata_path, sep = "\t", check.names = FALSE, na.strings = c("", "NA"))
sample_metadata$biological_line <- ifelse(
  sample_metadata$accession == "GSE303359",
  sub("_.*$", "", sample_metadata$sample_title),
  sample_metadata$cell_line
)
sample_metadata$treatment[is.na(sample_metadata$treatment)] <- "not_applicable"
sample_metadata$batch[is.na(sample_metadata$batch)] <- "not_reported"

validate_registered_design <- function(metadata, accession) {
  subset <- metadata[metadata$accession == accession, , drop = FALSE]
  if (nrow(subset) != 8 || length(unique(subset$sample_accession)) != 8) {
    stop(sprintf("%s must contain exactly eight distinct registered libraries", accession))
  }
  line_sizes <- table(subset$biological_line)
  if (length(line_sizes) != 4 || any(line_sizes != 2)) {
    stop(sprintf("%s must resolve to four biological lines with two libraries/conditions per line", accession))
  }
  genotype_per_line <- tapply(subset$genotype, subset$biological_line, function(value) length(unique(value)))
  lines_per_genotype <- table(unique(subset[, c("biological_line", "genotype")])$genotype)
  expected <- if (accession == "GSE303359") c(FSHD = 2L, Healthy = 2L) else c(SMA = 2L, CTRL = 2L)
  observed <- as.integer(lines_per_genotype[names(expected)])
  if (any(genotype_per_line != 1) || anyNA(observed) || !identical(observed, unname(expected))) {
    stop(sprintf("%s genotype-to-line mapping disagrees with the prespecified two-lines-per-genotype design", accession))
  }
  if (accession == "GSE303359") {
    treatments_per_line <- tapply(subset$treatment, subset$biological_line, function(value) setequal(value, c("un-injured", "injured")))
    if (!all(treatments_per_line)) stop("Each GSE303359 line must have baseline and injured libraries")
  }
  invisible(TRUE)
}

validate_registered_design(sample_metadata, "GSE303359")
validate_registered_design(sample_metadata, "GSE290980")
write_tsv(sample_metadata, file.path(LOGS, "registered_sample_metadata_locked.tsv"))

analysis_plan <- list(
  schema = "nmd-vcell-registered-scrna-analysis-plan/2.0",
  release = "multidisease_release_v02",
  seed = 20260805,
  datasets = list(
    GSE303359 = list(biological_units = "two FSHD and two healthy cell lines; baseline/injury repeated within line", inference = "voom-limma pseudobulk with line blocking; n=2 lines/genotype is directional"),
    GSE290980 = list(biological_units = "two SMA and two control donor-derived lines; two technical libraries per line", inference = "line-aggregated directional effects; no calibrated p-values")
  ),
  qc = list(min_features = 300, upper_complexity = "sample median + 5 MAD capped by the 99.5th percentile with conservative floor", mitochondrial = "sample-adaptive, maximum 25%", doublets = "sample-wise scran::doubletCells plus high-complexity outlier filter"),
  visualization = "PCA batch correction is visualization/clustering only; all effects use raw-count pseudobulk",
  candidate_output_contract = list(
    GSE303359 = c("DUX4", "DUXA", "ZSCAN4", "LEUTX", "MAPK8", "MAPK14"),
    GSE290980 = c("SMN1", "SMN2", "KIF5A"),
    filtered_candidate_policy = "retain explicit BELOW_EXPRESSION_FILTER rows with no imputed effect"
  ),
  claim_boundary = "No cells as replicates; no patient efficacy, target ranking or calibrated clinical prediction"
)
jsonlite::write_json(analysis_plan, file.path(OUT, "analysis_plan.json"), pretty = TRUE, auto_unbox = TRUE)

discover_mex <- function(accession) {
  directory <- file.path(DATA, accession, "mex")
  matrices <- sort(list.files(directory, pattern = "matrix\\.mtx\\.gz$", full.names = TRUE))
  if (length(matrices) != 8) stop(sprintf("%s expected 8 matrices; found %d", accession, length(matrices)))
  records <- lapply(matrices, function(matrix_path) {
    feature_path <- sub("matrix\\.mtx\\.gz$", "features.tsv.gz", matrix_path)
    barcode_path <- sub("matrix\\.mtx\\.gz$", "barcodes.tsv.gz", matrix_path)
    if (!file.exists(feature_path) || !file.exists(barcode_path)) stop(sprintf("Incomplete MEX trio for %s", matrix_path))
    sample_accession <- regmatches(basename(matrix_path), regexpr("GSM[0-9]+", basename(matrix_path)))
    data.frame(sample_accession = sample_accession, matrix = matrix_path, features = feature_path, barcodes = barcode_path)
  })
  do.call(rbind, records)
}

read_mex <- function(record) {
  features <- read.delim(gzfile(record$features), header = FALSE, sep = "\t", quote = "", comment.char = "")
  barcodes <- readLines(gzfile(record$barcodes))
  counts <- Matrix::readMM(gzfile(record$matrix))
  symbols <- as.character(features[[if (ncol(features) >= 2) 2 else 1]])
  if (nrow(counts) == length(barcodes) && ncol(counts) == length(symbols)) counts <- t(counts)
  if (nrow(counts) != length(symbols) || ncol(counts) != length(barcodes)) stop(sprintf("MEX dimensions disagree for %s", record$sample_accession))
  counts <- as(counts, "dgCMatrix")
  rownames(counts) <- make.unique(symbols)
  colnames(counts) <- paste(record$sample_accession, barcodes, sep = "_")
  counts
}

safe_upper <- function(values, floor_value) {
  med <- median(values)
  spread <- mad(values, constant = 1.4826)
  candidate <- if (is.finite(spread) && spread > 0) med + 5 * spread else as.numeric(quantile(values, 0.995))
  max(floor_value, min(candidate, as.numeric(quantile(values, 0.995))))
}

qc_sample <- function(accession, record, meta_row) {
  counts <- read_mex(record)
  object <- CreateSeuratObject(counts = counts, project = accession, min.cells = 3, min.features = 200)
  object[["percent.mt"]] <- PercentageFeatureSet(object, pattern = "^MT-")
  raw_cells <- ncol(object)
  upper_features <- safe_upper(object$nFeature_RNA, 2500)
  upper_counts <- safe_upper(object$nCount_RNA, 10000)
  mt_limit <- min(25, max(12, median(object$percent.mt) + 5 * mad(object$percent.mt, constant = 1.4826)))

  doublet_score <- rep(NA_real_, ncol(object))
  doublet_call <- rep(FALSE, ncol(object))
  doublet_method <- "scran_doubletCells"
  try_doublet <- tryCatch({
    sce <- SingleCellExperiment(list(counts = get_counts(object)))
    sce <- scuttle::logNormCounts(sce)
    score <- scran::doubletCells(sce)
    call <- scuttle::isOutlier(score, type = "higher", nmads = 3)
    list(score = as.numeric(score), call = as.logical(call))
  }, error = function(e) {
    doublet_method <<- paste0("scran_failed_predeclared_complexity_fallback:", sanitize_tsv_text(conditionMessage(e)))
    NULL
  })
  if (!is.null(try_doublet)) {
    doublet_score <- try_doublet$score
    doublet_call <- try_doublet$call
  }
  object$doublet_score <- doublet_score
  object$doublet_call <- doublet_call
  low_feature_fail <- object$nFeature_RNA < 300
  high_feature_fail <- object$nFeature_RNA > upper_features
  high_count_fail <- object$nCount_RNA > upper_counts
  mitochondrial_fail <- object$percent.mt > mt_limit
  keep <- !low_feature_fail & !high_feature_fail & !high_count_fail & !mitochondrial_fail & !doublet_call
  retained <- sum(keep)
  if (retained < 100) stop(sprintf("QC retained fewer than 100 cells for %s", record$sample_accession))
  object <- subset(object, cells = colnames(object)[keep])
  object$accession <- accession
  object$sample_accession <- meta_row$sample_accession
  object$sample_title <- meta_row$sample_title
  object$biological_line <- meta_row$biological_line
  object$genotype <- meta_row$genotype
  object$treatment <- meta_row$treatment
  object$batch <- meta_row$batch
  object$condition <- if (accession == "GSE303359") paste(meta_row$genotype, meta_row$treatment, sep = " · ") else meta_row$genotype
  qc <- data.frame(
    accession = accession, sample_accession = meta_row$sample_accession, sample_title = meta_row$sample_title,
    biological_line = meta_row$biological_line, genotype = meta_row$genotype, treatment = meta_row$treatment,
    batch = meta_row$batch, raw_cells = raw_cells, retained_cells = retained, retained_fraction = retained / raw_cells,
    min_features = 300, max_features = upper_features, max_counts = upper_counts, max_percent_mt = mt_limit,
    removed_low_features = sum(low_feature_fail), removed_high_features = sum(high_feature_fail),
    removed_high_counts = sum(high_count_fail), removed_high_mitochondrial = sum(mitochondrial_fail),
    model_doublets_before_filter = sum(doublet_call), doublet_method = doublet_method,
    stringsAsFactors = FALSE
  )
  list(object = object, qc = qc)
}

correct_pca_for_visualization <- function(object, accession) {
  embedding <- Embeddings(object, "pca")
  md <- object@meta.data
  if (accession == "GSE303359") {
    design <- model.matrix(~ genotype * treatment, md)
    batch <- factor(md$batch)
    method <- "limma_removeBatchEffect_on_PCs_preserving_genotype_by_treatment"
    corrected <- tryCatch(
      t(limma::removeBatchEffect(t(embedding), batch = batch, design = design)),
      error = function(e) {
        method <<- paste0("uncorrected_PCA_predeclared_rank_fallback:", conditionMessage(e))
        embedding
      }
    )
  } else {
    corrected <- center_nested_batch_pcs(embedding, md$biological_line, md$sample_accession)
    method <- "within_line_technical_library_PC_centering_preserving_line_centroids"
  }
  object[["pca.corrected"]] <- CreateDimReducObject(embeddings = corrected, key = "PCC_", assay = "RNA")
  object@misc$registered_visualization_correction <- method
  object
}

zscore <- function(values) {
  spread <- stats::sd(values, na.rm = TRUE)
  if (length(values) < 2 || !is.finite(spread) || spread == 0) return(rep(0, length(values)))
  as.numeric(scale(values))
}

annotate_clusters <- function(object, marker_sets) {
  average <- AverageExpression(object, assays = "RNA", group.by = "seurat_clusters", slot = "data", verbose = FALSE)$RNA
  average <- canonicalize_cluster_columns(average, object$seurat_clusters)
  raw_scores <- sapply(marker_sets, function(genes) {
    present <- intersect(genes, rownames(average))
    if (!length(present)) return(rep(NA_real_, ncol(average)))
    colMeans(as.matrix(average[present, , drop = FALSE]))
  })
  if (is.null(dim(raw_scores))) raw_scores <- matrix(raw_scores, ncol = 1)
  rownames(raw_scores) <- colnames(average)
  standardized <- apply(raw_scores, 2, zscore)
  if (is.null(dim(standardized))) standardized <- matrix(standardized, ncol = 1)
  rownames(standardized) <- colnames(average)
  colnames(standardized) <- names(marker_sets)
  labels <- colnames(standardized)[max.col(standardized, ties.method = "first")]
  score_table <- data.frame(cluster = rownames(standardized), marker_state = labels, standardized, check.names = FALSE)
  lookup <- setNames(score_table$marker_state, score_table$cluster)
  assigned <- unname(lookup[as.character(object$seurat_clusters)])
  if (anyNA(assigned)) stop("Marker-state assignment produced missing labels after cluster-ID normalization")
  object$marker_state <- assigned

  marker_rows <- lapply(seq_len(ncol(average)), function(index) {
    other <- if (ncol(average) > 1) rowMeans(as.matrix(average[, -index, drop = FALSE])) else rep(0, nrow(average))
    delta <- as.numeric(average[, index]) - other
    selected <- order(delta, decreasing = TRUE)[seq_len(min(25, length(delta)))]
    data.frame(cluster = colnames(average)[index], gene = rownames(average)[selected], average_log_expression = as.numeric(average[selected, index]), descriptive_delta = delta[selected])
  })
  list(object = object, scores = score_table, markers = do.call(rbind, marker_rows))
}

process_core <- function(accession, marker_sets) {
  records <- discover_mex(accession)
  objects <- vector("list", nrow(records))
  qc_rows <- vector("list", nrow(records))
  for (index in seq_len(nrow(records))) {
    sample_id <- records$sample_accession[index]
    meta_row <- sample_metadata[sample_metadata$accession == accession & sample_metadata$sample_accession == sample_id, , drop = FALSE]
    if (nrow(meta_row) != 1) stop(sprintf("Metadata mismatch for %s", sample_id))
    result <- qc_sample(accession, records[index, ], meta_row)
    objects[[index]] <- result$object
    qc_rows[[index]] <- result$qc
    message(sprintf("%s %s: %d cells retained", accession, sample_id, ncol(result$object)))
  }
  object <- merge(objects[[1]], y = objects[-1], merge.data = FALSE)
  if ("JoinLayers" %in% getNamespaceExports("SeuratObject")) object <- SeuratObject::JoinLayers(object)
  object <- NormalizeData(object, verbose = FALSE)
  object <- FindVariableFeatures(object, selection.method = "vst", nfeatures = 3000, verbose = FALSE)
  object <- ScaleData(object, features = VariableFeatures(object), vars.to.regress = "percent.mt", verbose = FALSE)
  object <- RunPCA(object, npcs = 30, features = VariableFeatures(object), verbose = FALSE)
  object <- correct_pca_for_visualization(object, accession)
  dimensions <- seq_len(min(30, ncol(Embeddings(object, "pca.corrected"))))
  object <- FindNeighbors(object, reduction = "pca.corrected", dims = dimensions, verbose = FALSE)
  object <- FindClusters(object, resolution = 0.35, random.seed = 20260805, verbose = FALSE)
  object <- RunUMAP(object, reduction = "pca.corrected", dims = dimensions, seed.use = 20260805, verbose = FALSE)
  annotation <- annotate_clusters(object, marker_sets)
  list(object = annotation$object, qc = do.call(rbind, qc_rows), annotation = annotation$scores, markers = annotation$markers)
}

module_matrix <- function(logcpm, gene_sets) {
  values <- lapply(gene_sets, function(genes) {
    present <- intersect(genes, rownames(logcpm))
    if (!length(present)) return(rep(NA_real_, ncol(logcpm)))
    colMeans(logcpm[present, , drop = FALSE])
  })
  result <- do.call(rbind, values)
  rownames(result) <- names(gene_sets)
  colnames(result) <- colnames(logcpm)
  result
}

embedding_frame <- function(object, max_cells = 50000) {
  frame <- cbind(object@meta.data, canonical_reduction_frame(Embeddings(object, "umap"), "UMAP"))
  frame$cell_id <- rownames(frame)
  if (nrow(frame) > max_cells) frame <- frame[sample(seq_len(nrow(frame)), max_cells), , drop = FALSE]
  frame
}

cell_metadata_frame <- function(object) {
  frame <- cbind(
    object@meta.data,
    canonical_reduction_frame(Embeddings(object, "umap"), "UMAP"),
    canonical_reduction_frame(Embeddings(object, "pca.corrected")[, 1:10, drop = FALSE], "PCC")
  )
  frame$cell_id <- rownames(frame)
  frame
}

pseudobulk_pca_plot <- function(logcpm, metadata, title, color_values) {
  variable <- order(apply(logcpm, 1, var), decreasing = TRUE)[seq_len(min(2000, nrow(logcpm)))]
  pcs <- prcomp(t(logcpm[variable, , drop = FALSE]), center = TRUE, scale. = FALSE)$x[, 1:2, drop = FALSE]
  frame <- cbind(metadata[match(rownames(pcs), metadata$group), , drop = FALSE], as.data.frame(pcs))
  ggplot(frame, aes(PC1, PC2, colour = genotype, label = biological_line)) + geom_hline(yintercept = 0, colour = "#D7E0E4", linewidth = 0.3) + geom_vline(xintercept = 0, colour = "#D7E0E4", linewidth = 0.3) + geom_point(size = 2.5) + geom_text(nudge_y = 0.4, size = 2.3, show.legend = FALSE) + scale_colour_manual(values = color_values) + labs(title = title, subtitle = "Raw-count pseudobulk; labels are biological lines") + theme_nmd()
}

analyze_fshd <- function() {
  core <- process_core("GSE303359", FSHD_STATE_MARKERS)
  object <- core$object
  write_tsv(core$qc, file.path(TABLES, "fshd_gse303359_sample_qc.tsv"))
  write_tsv(core$annotation, file.path(TABLES, "fshd_gse303359_cluster_annotation_scores.tsv"))
  write_tsv(core$markers, file.path(TABLES, "fshd_gse303359_descriptive_cluster_markers.tsv.gz"))

  pb <- aggregate_counts(object, object$sample_accession)
  smeta <- unique(object@meta.data[, c("sample_accession", "sample_title", "biological_line", "genotype", "treatment", "batch")])
  smeta <- smeta[match(colnames(pb), smeta$sample_accession), , drop = FALSE]
  smeta$group <- smeta$sample_accession
  smeta$genotype <- factor(smeta$genotype, levels = c("Healthy", "FSHD"))
  smeta$treatment <- factor(smeta$treatment, levels = c("un-injured", "injured"))
  smeta$batch <- factor(smeta$batch)
  design <- model.matrix(~ batch + genotype * treatment, smeta)
  required <- c("genotypeFSHD", "treatmentinjured", "genotypeFSHD:treatmentinjured")
  if (!all(required %in% colnames(design))) stop("FSHD design columns differ from the prespecified contrast matrix")
  dge <- edgeR::DGEList(pb)
  keep <- edgeR::filterByExpr(dge, design = design)
  dge <- edgeR::calcNormFactors(dge[keep, , keep.lib.sizes = FALSE])
  voom_object <- limma::voom(dge, design, plot = FALSE)
  correlation <- limma::duplicateCorrelation(voom_object, design, block = smeta$biological_line)$consensus
  fit <- limma::lmFit(voom_object, design, block = smeta$biological_line, correlation = correlation)
  contrast <- matrix(0, nrow = ncol(design), ncol = 4, dimnames = list(colnames(design), c("FSHD_vs_Healthy_baseline", "Injury_in_Healthy", "Injury_in_FSHD", "Genotype_by_injury_interaction")))
  contrast["genotypeFSHD", "FSHD_vs_Healthy_baseline"] <- 1
  contrast["treatmentinjured", "Injury_in_Healthy"] <- 1
  contrast["treatmentinjured", "Injury_in_FSHD"] <- 1
  contrast["genotypeFSHD:treatmentinjured", "Injury_in_FSHD"] <- 1
  contrast["genotypeFSHD:treatmentinjured", "Genotype_by_injury_interaction"] <- 1
  fit_contrast <- limma::eBayes(limma::contrasts.fit(fit, contrast), robust = TRUE)
  de <- do.call(rbind, lapply(colnames(contrast), function(name) {
    table <- limma::topTable(fit_contrast, coef = name, number = Inf, sort.by = "none")
    data.frame(gene = rownames(table), contrast = name, table, biological_unit = "cell_line_with_repeated_injury_state", inference = "directional_n2_lines_per_genotype", row.names = NULL)
  }))
  write_tsv(de, file.path(TABLES, "fshd_gse303359_line_blocked_pseudobulk_de.tsv.gz"))
  fshd_candidates <- complete_candidate_inventory(
    de,
    c("DUX4", "DUXA", "ZSCAN4", "LEUTX", "MAPK8", "MAPK14"),
    colnames(contrast),
    "logFC"
  )
  write_tsv(fshd_candidates, file.path(TABLES, "fshd_gse303359_registered_candidate_effects.tsv"))

  logcpm <- edgeR::cpm(dge, log = TRUE, prior.count = 2)
  modules <- module_matrix(logcpm, FSHD_MODULES)
  module_fit <- limma::lmFit(modules, design, block = smeta$biological_line, correlation = correlation)
  module_fit <- limma::eBayes(limma::contrasts.fit(module_fit, contrast), robust = TRUE)
  module_effects <- do.call(rbind, lapply(seq_len(ncol(contrast)), function(index) {
    estimate <- module_fit$coefficients[, index]
    se <- module_fit$stdev.unscaled[, index] * sqrt(module_fit$s2.post)
    critical <- qt(0.975, df = module_fit$df.total)
    data.frame(module = rownames(module_fit$coefficients), contrast = colnames(contrast)[index], estimate = estimate, standard_error = se, ci_low = estimate - critical * se, ci_high = estimate + critical * se, p_value = module_fit$p.value[, index], fdr = p.adjust(module_fit$p.value[, index], method = "BH"), biological_lines_per_genotype = 2, inference = "directional_small_n")
  }))
  write_tsv(module_effects, file.path(TABLES, "fshd_gse303359_module_effects.tsv"))
  module_long <- do.call(rbind, lapply(rownames(modules), function(module) data.frame(module = module, sample_accession = colnames(modules), score = as.numeric(modules[module, ]))))
  module_long <- merge(module_long, smeta, by = "sample_accession", all.x = TRUE)
  write_tsv(module_long, file.path(TABLES, "fshd_gse303359_module_by_sample.tsv"))

  write_tsv(cell_metadata_frame(object), file.path(TABLES, "fshd_gse303359_cell_metadata_and_embedding.tsv.gz"))
  plot_data <- embedding_frame(object)
  condition_colors <- c("Healthy · un-injured" = "#AAB5BB", "Healthy · injured" = COLORS[["blue"]], "FSHD · un-injured" = "#E8BD62", "FSHD · injured" = COLORS[["amber"]])
  p1 <- ggplot(plot_data, aes(UMAP_1, UMAP_2, colour = condition)) + geom_point(size = 0.08, alpha = 0.45) + scale_colour_manual(values = condition_colors) + labs(title = "FSHD membrane-injury atlas", subtitle = "Visualization only; effects use line-blocked pseudobulk", colour = "Condition") + theme_nmd() + theme(axis.text = element_blank(), axis.ticks = element_blank())
  p2 <- ggplot(module_long[module_long$module %in% c("DUX4_target", "Membrane_repair", "Ferroptosis"), ], aes(treatment, score, group = biological_line, colour = genotype)) + geom_line(linewidth = 0.55, alpha = 0.8) + geom_point(size = 1.6) + facet_wrap(~ module, scales = "free_y", nrow = 1) + scale_colour_manual(values = c(Healthy = COLORS[["grey"]], FSHD = COLORS[["amber"]])) + labs(title = "Biological-line module trajectories", subtitle = "Two lines per genotype; points are sample pseudobulks", x = NULL, y = "Mean logCPM") + theme_nmd() + theme(axis.text.x = element_text(angle = 25, hjust = 1), legend.position = "bottom")
  interaction_effects <- module_effects[module_effects$contrast == "Genotype_by_injury_interaction", ]
  p3 <- ggplot(interaction_effects, aes(estimate, reorder(module, estimate), xmin = ci_low, xmax = ci_high)) + geom_vline(xintercept = 0, linetype = 2, colour = COLORS[["grey"]]) + geom_errorbarh(height = 0.15, colour = COLORS[["amber"]]) + geom_point(size = 2, colour = COLORS[["amber"]]) + labs(title = "Genotype × injury module effect", subtitle = "95% CI shown; n=2 lines/genotype, directional", x = "Interaction effect (logCPM)", y = NULL) + theme_nmd()
  main_plot <- (p1 | p2 | p3) + plot_annotation(title = "FSHD GSE303359 · registered single-cell extension", subtitle = "Sample-wise QC, visualization-only batch correction and line-blocked raw-count pseudobulk", theme = theme(plot.title = element_text(face = "bold", colour = COLORS[["navy"]], size = 12), plot.subtitle = element_text(colour = COLORS[["grey"]], size = 8)))

  qc_long <- rbind(data.frame(sample = core$qc$sample_accession, status = "Raw", cells = core$qc$raw_cells), data.frame(sample = core$qc$sample_accession, status = "Retained", cells = core$qc$retained_cells))
  q1 <- ggplot(qc_long, aes(sample, cells, fill = status)) + geom_col(position = "dodge") + scale_fill_manual(values = c(Raw = "#CBD5DA", Retained = COLORS[["amber"]])) + labs(title = "Cells before/after QC", x = NULL, y = "Cells") + theme_nmd() + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  q2 <- ggplot(plot_data, aes(sample_accession, percent.mt, fill = genotype)) + geom_violin(scale = "width", linewidth = 0.2) + scale_fill_manual(values = c(Healthy = COLORS[["grey"]], FSHD = COLORS[["amber"]])) + labs(title = "Mitochondrial fraction after QC", x = NULL, y = "% mitochondrial") + theme_nmd() + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
  q3 <- ggplot(plot_data, aes(UMAP_1, UMAP_2, colour = sample_accession)) + geom_point(size = 0.07, alpha = 0.4) + labs(title = "Library mixing audit", colour = "Library") + theme_nmd() + theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "bottom")
  q4 <- pseudobulk_pca_plot(logcpm, smeta, "Pseudobulk PCA", c(Healthy = COLORS[["grey"]], FSHD = COLORS[["amber"]]))
  qc_plot <- (q1 | q2) / (q3 | q4) + plot_annotation(title = "FSHD GSE303359 · QC and biological-unit audit")
  save_plot_family(main_plot, file.path(FIGURES, "figure_fshd_registered_singlecell_extension_v02"), 10.8, 4.2)
  save_plot_family(qc_plot, file.path(FIGURES, "figure_s25_fshd_singlecell_qc_v02"), 7.2, 7.2)

  result <- list(main_plot = main_plot, qc_plot = qc_plot, raw_cells = sum(core$qc$raw_cells), retained_cells = sum(core$qc$retained_cells), module_effects = module_effects, top_interaction = de[de$contrast == "Genotype_by_injury_interaction", ][order(de$adj.P.Val[de$contrast == "Genotype_by_injury_interaction"]), ][1:10, c("gene", "logFC", "P.Value", "adj.P.Val")], visualization_correction = object@misc$registered_visualization_correction)
  rm(object, core)
  gc()
  result
}

analyze_sma <- function() {
  core <- process_core("GSE290980", SMA_STATE_MARKERS)
  object <- core$object
  write_tsv(core$qc, file.path(TABLES, "sma_gse290980_sample_qc.tsv"))
  write_tsv(core$annotation, file.path(TABLES, "sma_gse290980_cluster_annotation_scores.tsv"))
  write_tsv(core$markers, file.path(TABLES, "sma_gse290980_descriptive_cluster_markers.tsv.gz"))

  pb <- aggregate_counts(object, object$biological_line)
  lmeta <- unique(object@meta.data[, c("biological_line", "genotype")])
  lmeta <- lmeta[match(colnames(pb), lmeta$biological_line), , drop = FALSE]
  lmeta$group <- lmeta$biological_line
  dge <- edgeR::DGEList(pb)
  keep <- rowSums(edgeR::cpm(dge) > 1) >= 2
  dge <- edgeR::calcNormFactors(dge[keep, , keep.lib.sizes = FALSE])
  logcpm <- edgeR::cpm(dge, log = TRUE, prior.count = 2)
  gene_effects <- pairwise_effect_table(logcpm, lmeta, "SMA", "CTRL")
  names(gene_effects)[names(gene_effects) == "feature"] <- "gene"
  write_tsv(gene_effects, file.path(TABLES, "sma_gse290980_line_level_directional_effects.tsv.gz"))
  sma_candidates <- complete_candidate_inventory(
    gene_effects, c("SMN1", "SMN2", "KIF5A"), effect_column = "log2_effect"
  )
  write_tsv(sma_candidates, file.path(TABLES, "sma_gse290980_registered_candidate_effects.tsv"))

  modules <- module_matrix(logcpm, SMA_MODULES)
  module_effects <- pairwise_effect_table(modules, lmeta, "SMA", "CTRL")
  names(module_effects)[names(module_effects) == "feature"] <- "module"
  module_long <- do.call(rbind, lapply(rownames(modules), function(module) data.frame(module = module, biological_line = colnames(modules), score = as.numeric(modules[module, ]))))
  line_genotype <- setNames(as.character(lmeta$genotype), as.character(lmeta$biological_line))
  if (anyDuplicated(names(line_genotype))) stop("SMA biological-line genotype map is not unique")
  module_long$genotype <- unname(line_genotype[as.character(module_long$biological_line)])
  module_long$group <- module_long$biological_line
  if (anyNA(module_long$genotype)) stop("SMA module table contains an unmapped biological line")
  write_tsv(module_effects, file.path(TABLES, "sma_gse290980_module_directional_effects.tsv"))
  write_tsv(module_long, file.path(TABLES, "sma_gse290980_module_by_line.tsv"))

  composition_counts <- as.data.frame.matrix(table(object$biological_line, object$marker_state))
  if (!nrow(composition_counts) || !ncol(composition_counts)) stop("SMA marker-state composition table is empty")
  composition <- prop.table(as.matrix(composition_counts), 1)
  composition_long <- do.call(rbind, lapply(colnames(composition), function(state) data.frame(biological_line = rownames(composition), marker_state = state, proportion = composition[, state])))
  composition_long$genotype <- unname(line_genotype[as.character(composition_long$biological_line)])
  composition_long$group <- composition_long$biological_line
  if (anyNA(composition_long$genotype)) stop("SMA composition table contains an unmapped biological line")
  composition_effects <- do.call(rbind, lapply(unique(composition_long$marker_state), function(state) {
    subset <- composition_long[composition_long$marker_state == state, ]
    case <- subset$proportion[subset$genotype == "SMA"]
    control <- subset$proportion[subset$genotype == "CTRL"]
    data.frame(marker_state = state, sma_mean = mean(case), control_mean = mean(control), proportion_difference = mean(case) - mean(control), directional_consistency = pairwise_consistency(case, control), p_value = NA_real_, fdr = NA_real_, inference = "directional_n2_lines_per_genotype")
  }))
  write_tsv(composition_long, file.path(TABLES, "sma_gse290980_cell_state_composition_by_line.tsv"))
  write_tsv(composition_effects, file.path(TABLES, "sma_gse290980_cell_state_composition_effects.tsv"))

  state_line_counts <- as.data.frame(table(
    biological_line = factor(object$biological_line, levels = unique(lmeta$biological_line)),
    marker_state = object$marker_state
  ), stringsAsFactors = FALSE)
  state_minimums <- aggregate(Freq ~ marker_state, state_line_counts, min)
  names(state_minimums)[2] <- "minimum_cells_across_four_lines"
  state_minimums$eligible_for_state_effect <- state_minimums$minimum_cells_across_four_lines >= 30
  write_tsv(state_minimums, file.path(TABLES, "sma_gse290980_state_effect_eligibility.tsv"))
  eligible_states <- as.character(state_minimums$marker_state[state_minimums$eligible_for_state_effect])
  if (!length(eligible_states)) stop("No SMA marker-defined state has at least 30 cells in every biological line")
  state_cells <- colnames(object)[object$marker_state %in% eligible_states]
  state_groups <- paste(object$biological_line[state_cells], object$marker_state[state_cells], sep = "__")
  state_pb <- aggregate_counts(subset(object, cells = state_cells), state_groups)
  state_dge <- edgeR::calcNormFactors(edgeR::DGEList(state_pb))
  state_keep <- rowSums(edgeR::cpm(state_dge) > 1) >= 2
  state_dge <- state_dge[state_keep, , keep.lib.sizes = FALSE]
  state_logcpm <- edgeR::cpm(state_dge, log = TRUE, prior.count = 2)
  parsed <- do.call(rbind, strsplit(colnames(state_logcpm), "__", fixed = TRUE))
  state_meta <- data.frame(group = colnames(state_logcpm), biological_line = parsed[, 1], marker_state = parsed[, 2])
  state_meta$genotype <- unname(line_genotype[as.character(state_meta$biological_line)])
  if (anyNA(state_meta$genotype)) stop("SMA state-effect table contains an unmapped biological line")
  state_effects <- do.call(rbind, lapply(unique(state_meta$marker_state), function(state) {
    subset <- state_meta[state_meta$marker_state == state, ]
    if (sum(subset$genotype == "SMA") < 2 || sum(subset$genotype == "CTRL") < 2) return(NULL)
    values <- state_logcpm[, subset$group, drop = FALSE]
    metadata <- subset
    metadata$group <- subset$group
    pairwise_effect_table(values, metadata, "SMA", "CTRL", state)
  }))
  names(state_effects)[names(state_effects) == "feature"] <- "gene"
  write_tsv(state_effects, file.path(TABLES, "sma_gse290980_state_specific_directional_effects.tsv.gz"))
  write_tsv(cell_metadata_frame(object), file.path(TABLES, "sma_gse290980_cell_metadata_and_embedding.tsv.gz"))

  plot_data <- embedding_frame(object)
  state_levels <- sort(unique(plot_data$marker_state))
  state_palette <- setNames(grDevices::hcl.colors(length(state_levels), "Dark 3"), state_levels)
  p1 <- ggplot(plot_data, aes(UMAP_1, UMAP_2, colour = marker_state)) + geom_point(size = 0.08, alpha = 0.5) + scale_colour_manual(values = state_palette) + labs(title = "SMA spinal-organoid atlas", subtitle = "Marker-defined states; labels require orthogonal confirmation", colour = "State") + theme_nmd() + theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "bottom")
  p2 <- ggplot(composition_long, aes(biological_line, proportion, fill = marker_state)) + geom_col(width = 0.75) + scale_fill_manual(values = state_palette) + facet_grid(~ genotype, scales = "free_x", space = "free_x") + labs(title = "Cell-state composition by line", subtitle = "Technical libraries aggregated before comparison", x = "Biological line", y = "Proportion", fill = "State") + theme_nmd() + theme(legend.position = "none")
  candidate_state <- state_effects[state_effects$gene %in% c("SMN1", "SMN2", "KIF5A", "STMN2", "CHAT"), ]
  if (!nrow(candidate_state)) candidate_state <- state_effects[order(abs(state_effects$log2_effect), decreasing = TRUE)[1:min(20, nrow(state_effects))], ]
  p3 <- ggplot(candidate_state, aes(marker_state, gene, fill = log2_effect)) + geom_tile(colour = "white", linewidth = 0.2) + scale_fill_gradient2(low = COLORS[["purple"]], mid = "white", high = COLORS[["teal"]], midpoint = 0) + labs(title = "State-localized directional effects", subtitle = "SMA − control; no cell-level p-values", x = NULL, y = NULL, fill = "log2 effect") + theme_nmd() + theme(axis.text.x = element_text(angle = 35, hjust = 1))
  main_plot <- (p1 | p2 | p3) + plot_annotation(title = "SMA GSE290980 · registered single-cell extension", subtitle = "Technical libraries collapse to four donor-derived lines; all disease effects are directional", theme = theme(plot.title = element_text(face = "bold", colour = COLORS[["navy"]], size = 12), plot.subtitle = element_text(colour = COLORS[["grey"]], size = 8)))

  qc_long <- rbind(data.frame(sample = core$qc$sample_accession, status = "Raw", cells = core$qc$raw_cells), data.frame(sample = core$qc$sample_accession, status = "Retained", cells = core$qc$retained_cells))
  q1 <- ggplot(qc_long, aes(sample, cells, fill = status)) + geom_col(position = "dodge") + scale_fill_manual(values = c(Raw = "#CBD5DA", Retained = COLORS[["teal"]])) + labs(title = "Cells before/after QC", x = NULL, y = "Cells") + theme_nmd() + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  q2 <- ggplot(plot_data, aes(sample_accession, percent.mt, fill = genotype)) + geom_violin(scale = "width", linewidth = 0.2) + scale_fill_manual(values = c(CTRL = COLORS[["grey"]], SMA = COLORS[["teal"]])) + labs(title = "Mitochondrial fraction after QC", x = NULL, y = "% mitochondrial") + theme_nmd() + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
  q3 <- ggplot(plot_data, aes(UMAP_1, UMAP_2, colour = sample_accession)) + geom_point(size = 0.07, alpha = 0.4) + labs(title = "Technical-library mixing audit", colour = "Library") + theme_nmd() + theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "bottom")
  q4 <- pseudobulk_pca_plot(logcpm, lmeta, "Line-level pseudobulk PCA", c(CTRL = COLORS[["grey"]], SMA = COLORS[["teal"]]))
  qc_plot <- (q1 | q2) / (q3 | q4) + plot_annotation(title = "SMA GSE290980 · QC and biological-unit audit")
  save_plot_family(main_plot, file.path(FIGURES, "figure_sma_registered_singlecell_extension_v02"), 10.8, 4.2)
  save_plot_family(qc_plot, file.path(FIGURES, "figure_s26_sma_singlecell_qc_v02"), 7.2, 7.2)

  kif5a <- gene_effects[gene_effects$gene == "KIF5A", , drop = FALSE]
  result <- list(main_plot = main_plot, qc_plot = qc_plot, raw_cells = sum(core$qc$raw_cells), retained_cells = sum(core$qc$retained_cells), module_effects = module_effects, composition_effects = composition_effects, kif5a = kif5a, visualization_correction = object@misc$registered_visualization_correction)
  rm(object, core)
  gc()
  result
}

fshd <- analyze_fshd()
sma <- analyze_sma()

unified <- (fshd$main_plot / sma$main_plot) + plot_annotation(
  title = "Figure 9 | Registered single-cell extensions refine disease-state localization without pseudo-replication",
  subtitle = "FSHD effects use line-blocked repeated-injury pseudobulk; SMA effects aggregate technical libraries to two lines per genotype.",
  caption = "Cells support state annotation and visualization only. Biological lines—not cells—define disease contrasts; all n=2/group results remain directional.",
  theme = theme(plot.title = element_text(face = "bold", colour = COLORS[["navy"]], size = 13), plot.subtitle = element_text(colour = COLORS[["grey"]], size = 8), plot.caption = element_text(colour = COLORS[["red"]], size = 7))
)
save_plot_family(unified, file.path(FIGURES, "figure_9_registered_singlecell_extensions_v02"), 10.8, 8.5)

best_fshd <- fshd$module_effects[fshd$module_effects$contrast == "Genotype_by_injury_interaction", ]
best_fshd <- best_fshd[order(best_fshd$fdr), , drop = FALSE]
largest_sma_state <- sma$composition_effects[order(abs(sma$composition_effects$proportion_difference), decreasing = TRUE), , drop = FALSE]
kif5a_text <- if (nrow(sma$kif5a)) sprintf("KIF5A line-level log2 effect %.3f (%s)", sma$kif5a$log2_effect[[1]], sma$kif5a$directional_consistency[[1]]) else "KIF5A was below the line-level expression filter"

summary <- list(
  schema = "nmd-vcell-registered-singlecell-summary/2.0",
  release = "multidisease_release_v02",
  generated_at = format(Sys.time(), "%Y-%m-%dT%H:%M:%S%z"),
  status = "ANALYZED_WITH_DIRECTIONAL_SMALL_N_BOUNDARY",
  datasets = list(
    GSE303359 = list(status = "ANALYZED", libraries = 8, biological_lines = 4, lines_per_genotype = 2, raw_cells = fshd$raw_cells, retained_cells = fshd$retained_cells, statistical_strategy = "line-blocked voom-limma pseudobulk with repeated baseline/injury measures", visualization_correction = fshd$visualization_correction, leading_module_interaction = if (nrow(best_fshd)) as.list(best_fshd[1, c("module", "estimate", "ci_low", "ci_high", "p_value", "fdr")]) else NULL),
    GSE290980 = list(status = "ANALYZED", libraries = 8, biological_lines = 4, lines_per_genotype = 2, raw_cells = sma$raw_cells, retained_cells = sma$retained_cells, statistical_strategy = "technical-library aggregation followed by line-level directional effects without calibrated p-values", visualization_correction = sma$visualization_correction, largest_composition_shift = if (nrow(largest_sma_state)) as.list(largest_sma_state[1, c("marker_state", "proportion_difference", "directional_consistency")]) else NULL, kif5a = kif5a_text)
  ),
  candidate_update = list(promoted_to_validated_target = 0, interpretation = "The single-cell extensions localize states and refine experiment design; they do not promote any therapeutic target or create a cross-disease score."),
  claim_boundary = "n=2 biological lines per genotype; directional evidence only; no efficacy, clinical or patient-calibrated claims"
)
jsonlite::write_json(summary, file.path(OUT, "singlecell_summary.json"), pretty = TRUE, auto_unbox = TRUE, na = "null")

report <- c(
  "# Registered single-cell extension report v02",
  "",
  "## Design lock",
  "",
  "- GSE303359: eight libraries, four biological myoblast lines (two FSHD, two healthy), repeated baseline/6 h injury within line.",
  "- GSE290980: eight libraries, four donor-derived organoid lines (two SMA, two control), two technical libraries per line.",
  "- Cells are used for QC, visualization and marker-defined state localization. They are never inferential replicates.",
  "- FSHD uses line-blocked raw-count pseudobulk; SMA technical libraries are aggregated to line and reported directionally without calibrated p-values.",
  "",
  "## FSHD result",
  "",
  if (nrow(best_fshd)) sprintf("The leading prespecified module interaction was %s (effect %.3f, 95%% CI %.3f to %.3f, FDR %.3g). With two lines per genotype this remains directional and requires independent-line replication.", best_fshd$module[[1]], best_fshd$estimate[[1]], best_fshd$ci_low[[1]], best_fshd$ci_high[[1]], best_fshd$fdr[[1]]) else "No module interaction was estimable.",
  "",
  "## SMA result",
  "",
  if (nrow(largest_sma_state)) sprintf("The largest marker-defined composition difference was %s (SMA minus control %.3f; %s). Labels are transcript-marker definitions, not orthogonally validated cell identities.", largest_sma_state$marker_state[[1]], largest_sma_state$proportion_difference[[1]], largest_sma_state$directional_consistency[[1]]) else "No composition effect was estimable.",
  paste0("- ", kif5a_text, "."),
  "",
  "## Candidate consequence",
  "",
  "No candidate is promoted to VALIDATED. The new data refine state localization and the next experiment only.",
  "",
  "## Limitations",
  "",
  "Two biological lines per genotype cannot support calibrated population inference. Batch correction is restricted to PCA visualization/clustering. Raw-count pseudobulk retains the biological unit. Marker-defined states require orthogonal confirmation, and no result is a treatment recommendation."
)
writeLines(report, file.path(OUT, "REGISTERED_SINGLECELL_ANALYSIS_REPORT_v02.md"))
writeLines(capture.output(sessionInfo()), file.path(LOGS, "sessionInfo.txt"))
writeLines(capture.output(analysis_plan), file.path(LOGS, "analysis_plan_print.txt"))

cat(jsonlite::toJSON(list(status = "PASS", release = "multidisease_release_v02", figures = 5, datasets = summary$datasets, candidate_promotions = 0), auto_unbox = TRUE, pretty = TRUE), "\n")