#!/usr/bin/env Rscript
# =============================================================================
# BioF3 Module 10 (教程 09): 空间转录组真实示例 — Visium 小鼠脑 (Anterior)
# =============================================================================
#
# 这个脚本在单细胞实践教程 09 章（空间转录组学）里用到。它：
#   1. 通过 SeuratData::stxBrain 加载 10x Visium 小鼠脑 Anterior 切片
#   2. QC（nCount/nFeature 在切片上的分布）
#   3. SCTransform + PCA + 聚类 + UMAP
#   4. 聚类在组织上的空间分布
#   5. 空间变异基因（Moran's I）
#   6. 对几个已知解剖学 marker 画空间分布（Hpca 海马 / Ttr 脉络丛 / Mbp 胼胝体 / Cux2 皮层）
#   7. 输出 6 张真实 PNG 到 static/img/tutorial/single-cell/module10/
#
# 数据来源：
#   SeuratData::stxBrain anterior1 (10x Visium 小鼠脑前部，~136 MB 的数据包)
#   包是通过 install.packages('stxBrain.SeuratData',repos='http://seurat.nygenome.org/') 装
#   脚本里用 requireSpatial() 自动检查/提示
#
# 用法：
#   Rscript scripts/module10_spatial_sci.R
#
# 可选环境变量：
#   BIOF3_OUTPUT_DIR=/path/to/static/img/tutorial/single-cell/module10
#
# 依赖：
#   Seurat, SeuratData, stxBrain.SeuratData, ggplot2, dplyr, patchwork, RColorBrewer
# =============================================================================

options(stringsAsFactors = FALSE)
set.seed(42)

# ---- 路径 ----------------------------------------------------------------
args <- commandArgs(trailingOnly = FALSE)
file_arg <- grep("^--file=", args, value = TRUE)
script_path <- if (length(file_arg) > 0) {
  normalizePath(sub("^--file=", "", file_arg[[1]]), mustWork = TRUE)
} else {
  normalizePath("scripts/single-cell/sc10_spatial_sci.R", mustWork = FALSE)
}
script_dir <- dirname(script_path)
project_root <- if (basename(script_dir) == "scripts" && basename(dirname(script_dir)) == "static") {
  dirname(dirname(script_dir))
} else {
  dirname(script_dir)
}
output_dir <- Sys.getenv(
  "BIOF3_OUTPUT_DIR",
  file.path(project_root, "static", "img", "tutorial", "modules", "module10")
)
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)

# ---- 依赖 ----------------------------------------------------------------
suppressPackageStartupMessages({
  library(Seurat)
  library(SeuratData)
  library(ggplot2)
  library(dplyr)
  library(patchwork)
  library(RColorBrewer)
})

# 检查 stxBrain 是否已安装
if (!"stxBrain" %in% InstalledData()$Dataset) {
  message("Installing stxBrain via SeuratData (~136 MB)...")
  InstallData("stxBrain")
}

biof3_colors <- c(
  green = "#0f766e", mint  = "#43d1ae", blue  = "#2563eb",
  amber = "#f59e0b", red   = "#dc2626", slate = "#475569",
  violet = "#7c3aed", pink = "#ec4899"
)

theme_biof3 <- function(base_size = 12) {
  theme_classic(base_size = base_size) +
    theme(
      axis.text  = element_text(color = "#111827"),
      axis.title = element_text(color = "#111827", face = "bold"),
      plot.title = element_text(color = "#12342f", face = "bold", hjust = 0),
      plot.subtitle = element_text(color = "#526760", hjust = 0),
      legend.title = element_text(face = "bold"),
      legend.position = "right",
      strip.background = element_blank(),
      strip.text = element_text(color = "#12342f", face = "bold")
    )
}

save_plot <- function(p, name, width = 8, height = 5) {
  ggsave(file.path(output_dir, name), p, width = width, height = height, dpi = 300, bg = "white")
}

# ---- 加载 -----------------------------------------------------------------
message("加载 stxBrain anterior1 ...")
brain <- LoadData("stxBrain", type = "anterior1")
message("spots: ", ncol(brain), "  features: ", nrow(brain))

# ---- 图 1：QC 在切片上的空间分布 ----------------------------------------
message("图 1：nCount / nFeature 在组织上的分布")
p1a <- SpatialFeaturePlot(brain, features = "nCount_Spatial") +
  labs(title = "Total UMI per spot") +
  theme_biof3() + theme(axis.line = element_blank(), axis.ticks = element_blank(),
                        axis.text = element_blank(), axis.title = element_blank())
p1b <- SpatialFeaturePlot(brain, features = "nFeature_Spatial") +
  labs(title = "Genes detected per spot") +
  theme_biof3() + theme(axis.line = element_blank(), axis.ticks = element_blank(),
                        axis.text = element_blank(), axis.title = element_blank())
p1 <- (p1a | p1b) + plot_annotation(
  title = "Spatial QC: sequencing depth across tissue",
  subtitle = "Check for edge dropout or local artifacts (bubbles, contamination)",
  theme = theme(plot.title = element_text(face = "bold", color = "#12342f"),
                plot.subtitle = element_text(color = "#526760"))
)
save_plot(p1, "01-qc-spatial.png", width = 12, height = 5.5)

# ---- SCTransform + PCA + 聚类 + UMAP -------------------------------------
message("SCTransform ...")
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
message("PCA / Neighbors / Clusters / UMAP ...")
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30, verbose = FALSE)
brain <- FindClusters(brain, verbose = FALSE, resolution = 0.5)
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30, verbose = FALSE)

# ---- 图 2：UMAP 聚类 ------------------------------------------------------
message("图 2：UMAP 聚类")
p2 <- DimPlot(brain, reduction = "umap", label = TRUE, repel = TRUE) +
  labs(title = "UMAP: SCT-normalized clustering",
       subtitle = paste0(length(levels(brain$seurat_clusters)), " clusters")) +
  theme_biof3() + theme(legend.position = "right")
save_plot(p2, "02-umap-clusters.png", width = 8.5, height = 6)

# ---- 图 3：聚类在组织切片上的空间分布 -----------------------------------
message("图 3：聚类在组织上的分布")
p3 <- SpatialDimPlot(brain, label = TRUE, label.size = 3) +
  labs(title = "Cluster distribution on tissue section",
       subtitle = "Cortex, hippocampus, corpus callosum, and choroid plexus map to distinct clusters") +
  theme(plot.title = element_text(face = "bold", color = "#12342f"),
        plot.subtitle = element_text(color = "#526760"))
save_plot(p3, "03-clusters-spatial.png", width = 9, height = 7)

# ---- 图 4：解剖学 marker 的空间表达 -------------------------------------
message("图 4：已知解剖 marker 空间表达")
# 常见 marker：
#   Hpca  = 海马
#   Ttr   = 脉络丛
#   Mbp   = 胼胝体 / 白质
#   Cux2  = 皮层 II/III 层
marker_genes <- intersect(c("Hpca", "Ttr", "Mbp", "Cux2"), rownames(brain))
message("使用的 marker：", paste(marker_genes, collapse = ", "))
p4 <- SpatialFeaturePlot(brain, features = marker_genes, ncol = 2, alpha = c(0.1, 1)) &
  theme(plot.title = element_text(face = "bold", color = "#12342f"))
save_plot(p4, "04-anatomy-markers.png", width = 10, height = 9)

# ---- 图 5：空间变异基因 -------------------------------------------------
message("图 5：空间变异基因（Moran's I）")
# Seurat 5.x 的 FindSpatiallyVariableFeatures 对 Assay5 默认用 "moransi"
brain_svg <- tryCatch(
  FindSpatiallyVariableFeatures(
    brain,
    assay            = "SCT",
    features         = VariableFeatures(brain)[1:1000],
    selection.method = "moransi"
  ),
  error = function(e) {
    message("FindSpatiallyVariableFeatures 失败：", conditionMessage(e))
    message("退回到用 VariableFeatures 的前 6 个画")
    NULL
  }
)

if (!is.null(brain_svg)) {
  top_svg <- tryCatch(
    head(SpatiallyVariableFeatures(brain_svg, selection.method = "moransi"), 6),
    error = function(e) head(VariableFeatures(brain_svg), 6)
  )
} else {
  brain_svg <- brain
  top_svg <- head(VariableFeatures(brain), 6)
}
message("Top 6 空间变异基因：", paste(top_svg, collapse = ", "))
p5 <- SpatialFeaturePlot(brain_svg, features = top_svg, ncol = 3, alpha = c(0.1, 1))
save_plot(p5, "05-spatially-variable.png", width = 13, height = 8)

# ---- 图 6：单一基因 UMAP + 空间双图 -------------------------------------
message("图 6：单基因 UMAP + 空间联合视图 (Mbp)")
g_focus <- if ("Mbp" %in% rownames(brain)) "Mbp" else top_svg[1]
p6a <- FeaturePlot(brain, features = g_focus, reduction = "umap") +
  scale_color_gradient(low = "lightgrey", high = biof3_colors["violet"]) +
  labs(title = paste0(g_focus, " (UMAP)")) +
  theme_biof3()
p6b <- SpatialFeaturePlot(brain, features = g_focus, alpha = c(0.1, 1)) +
  labs(title = paste0(g_focus, " (tissue section)")) +
  theme(plot.title = element_text(face = "bold", color = "#12342f"))
p6 <- (p6a | p6b) + plot_annotation(
  title = paste0(g_focus, ": UMAP vs tissue section"),
  subtitle = "Linking UMAP clusters to anatomical structures",
  theme = theme(plot.title = element_text(face = "bold", color = "#12342f"),
                plot.subtitle = element_text(color = "#526760"))
)
save_plot(p6, "06-gene-umap-vs-spatial.png", width = 12, height = 5.5)

# ---- 输出清单 -----------------------------------------------------------
message("=== 已生成 module10 空间图 ===")
print(list.files(output_dir, pattern = "\\.png$", full.names = FALSE))
message("输出目录：", normalizePath(output_dir))
message("完成。")

sessionInfo()