# ============================================================
# 03_preprocess.R —— 预处理与质控（对应教程第 08 章）
# 用法：Rscript code/03_preprocess.R（或 RStudio 逐段运行）
# 前置：已运行 02_download.R，数据在 data/raw/visium/
# ============================================================

library(Seurat)
library(patchwork)

## 1. 读入 Visium 数据
# 注意：Load10X_Spatial 要求目录下有 filtered_feature_bc_matrix.h5 和 spatial/ 子目录
obj <- Load10X_Spatial(
  data.dir = "data/raw/visium",
  filename = "filtered_feature_bc_matrix.h5"
)
obj

## 2. QC 指标计算
# nCount_Spatial：每个 spot 的 UMI 总数（文库深度）
# nFeature_Spatial：每个 spot 检测到的基因数
# percent.mt：线粒体基因比例（组织坏死/细胞状态指标）
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")

## 3. 先看分布，再定阈值（不要拍脑袋！）
p1 <- VlnPlot(obj, features = c("nCount_Spatial", "nFeature_Spatial", "percent.mt"),
              ncol = 3, pt.size = 0) & theme(legend.position = "none")
p2 <- SpatialFeaturePlot(obj, features = c("nCount_Spatial", "nFeature_Spatial", "percent.mt"))
print(p1 / p2)

## 4. 过滤低质量 spot（阈值以分布为准，这里是常见示例值）
obj <- subset(obj,
              subset = nFeature_Spatial > 200 &
                       nCount_Spatial > 500 &
                       percent.mt < 25)
cat("过滤后 spot 数：", ncol(obj), "\n")

## 5. 归一化 + 高变基因（SCTransform 一步完成）
obj <- SCTransform(obj, assay = "Spatial", verbose = FALSE)

## 6. PCA → 聚类 → UMAP
obj <- RunPCA(obj, npcs = 30)
obj <- FindNeighbors(obj, dims = 1:30)
obj <- FindClusters(obj, resolution = 0.6)
obj <- RunUMAP(obj, dims = 1:30)

## 7. 可视化
p_umap  <- DimPlot(obj, reduction = "umap", label = TRUE) +
           ggtitle("聚类（UMAP）")
p_spat  <- SpatialDimPlot(obj, label = TRUE, label.size = 3) +
           ggtitle("聚类（空间）")
print(p_umap | p_spat)

## 8. 保存中间结果（后续 04/05/06 章都从这里继续）
dir.create("data/processed", showWarnings = FALSE)
saveRDS(obj, "data/processed/visium_preprocessed.rds")
cat("预处理完成，已保存 data/processed/visium_preprocessed.rds\n")
