# ============================================================
# 05_spatial_domain.R —— 空间生态位识别（对应教程第 10 章）
# 用法：Rscript code/05_spatial_domain.R（或 RStudio 逐段运行）
# 前置：已运行 04_annotation.R，读取 data/processed/visium_annotated.rds
# 说明：生态位识别三步：①BayesSpace 空间聚类分割 ②Squidpy 邻域富集验证
#       （Python，见 code/04b_cell2location_squidpy.py）③命名与特征化。
# ============================================================

library(Seurat)
library(BayesSpace)
library(SingleCellExperiment)
library(patchwork)

obj <- readRDS("data/processed/visium_annotated.rds")

## 1. 转成 SingleCellExperiment（BayesSpace 的输入）
sce <- as.SingleCellExperiment(obj)

## 2. 空间预处理（PCA 降到 15 维，供聚类使用）
sce <- spatialPreprocess(sce, platform = "Visium", n.PCs = 15)

## 3. 选择聚类数 q：先跑 qTune 看哪个 q 的伪似然最高（较慢，可选）
# sce <- qTune(sce, qs = seq(4, 8), platform = "Visium", d = 15)
# qPlot(sce)

## 4. 空间聚类（MRF 空间先验，让相邻 spot 更可能同类）
q <- 6   # 教程建议 4-7 之间多试几个，结合生物学可解释性选
sce <- spatialCluster(sce, q = q, platform = "Visium", d = 15, init = "kmeans")

## 5. 结果映射回 Seurat，并按生物学含义重命名生态位
obj$spatial_domain <- as.character(sce$spatial.cluster)

# 查看每个空间域的组成，据此命名（示例映射，实际请结合 marker 判断）
# 例：1=肿瘤核心，2=侵袭前沿，3=免疫浸润区，4=基质区……
# obj$niche <- dplyr::recode(obj$spatial_domain,
#                            "1" = "肿瘤核心", "2" = "侵袭前沿",
#                            "3" = "免疫浸润区", "4" = "基质区")

p1 <- SpatialDimPlot(obj, group.by = "spatial_domain", label = FALSE) +
      ggtitle("空间域（BayesSpace）")
print(p1)

## 6. 各空间域的特征基因（为命名提供证据）
Idents(obj) <- obj$spatial_domain
domain_markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25,
                                 logfc.threshold = 0.5, assay = "SCT")
print(domain_markers %>% dplyr::group_by(cluster) %>%
        dplyr::slice_max(n = 5, order_by = avg_log2FC) %>% as.data.frame())

## 7. 保存
obj$niche <- obj$spatial_domain   # 占位：实际分析中替换为命名后的生态位
saveRDS(obj, "data/processed/visium_niche.rds")
cat("生态位识别完成，已保存 data/processed/visium_niche.rds\n")
cat("下一步：运行 code/04b_cell2location_squidpy.py 做邻域富集验证，\n")
cat("然后依次运行 06-08 号脚本做分区细胞通讯分析。\n")
