# 第 08 章 预处理与质控

> 本章配套脚本：`code/03_preprocess.R`；本章图：`figures/fig07_qc.png`、`figures/fig08_clustering.png`

拿到 Visium 数据后，你手里是一张“spot × 基因”的计数矩阵，外加每个 spot 的组织坐标。直接分析它是个坏主意：捕获过程中混进了背景 RNA、低质量 spot，文库深度也参差不齐。预处理与质控（quality control, QC）就是把脏数据清洗成可靠数据的第一步，也是后面所有分析的地基。地基不稳，后面的聚类、注释、通讯推断全都白搭。这一章我们手把手把这一步做完。

## 本章目标

1. 学会用 `Load10X_Spatial` 读入 Visium 数据并理解它的目录结构。
2. 理解三个核心 QC 指标：`nCount_Spatial`（UMI 总数）、`nFeature_Spatial`（基因数）、`percent.mt`（线粒体比例）。
3. 掌握“先看分布、再定阈值”的过滤方法，而不是拍脑袋。
4. 学会归一化（SCTransform）、高变基因、PCA、聚类、UMAP 的完整流程。
5. 会用 `SpatialDimPlot` / `SpatialFeaturePlot` 做空间可视化，并正确解读 QC 图与聚类图。

## 8.1 读入数据：Load10X_Spatial

第 07 章我们从 10x Genomics 官网下载了人乳腺癌 demo 数据。下载完成后，`data/visium/` 目录里应该包含这些内容：

- `filtered_feature_bc_matrix.h5`：spot × 基因的计数矩阵（HDF5 格式，包含基因名、spot 条形码）；
- `spatial/` 目录：`tissue_positions_list.csv`（每个 spot 的组织坐标）、`scalefactors_json.json`（图像缩放因子）、`tissue_hires_image.png` 与 `tissue_lowres_image.png`（组织切片染色图）。

`Load10X_Spatial` 会把这些信息整合成一个 Seurat 对象。注意两个参数：`data.dir` 指向包含 `spatial/` 子目录的文件夹（而不是 h5 文件本身）；`filename` 是 h5 文件的文件名。第一个代码块先安装并加载 Seurat v5：

```r
# 若未安装：install.packages("Seurat") 或 BiocManager::install("Seurat")
library(Seurat)
library(dplyr)      # 管道操作 %>%

# 读入 Visium 数据：data.dir 必须指向含 spatial/ 子目录的文件夹
obj <- Load10X_Spatial(
  data.dir  = "data/visium/",
  filename  = "filtered_feature_bc_matrix.h5"
)

# 看一眼对象结构：应该显示 1 个 assay（Spatial）和约 5000 个 spot
obj
```

**运行结果解读**：终端会打印类似 `15000 features across 4780 samples within 1 assay` 的信息。“features”是基因数，“samples”在这里就是 spot 数（10x Visium 捕获区约 5000 个 spot，过滤后通常剩 4000–5000 个）。Seurat v5 会自动把 assay 命名为 `Spatial`，图像信息存在 `obj@images` 里，坐标存在 `obj@images$slice1@coordinates`。此时 `obj` 里还没有任何质控指标，下一步来算。

## 8.2 QC 指标的含义：先看懂三个数字

每个 spot 是一个直径约 55 μm 的小圆点，覆盖约 1–10 个细胞（见第 02 章）。由于捕获和建库过程的随机性，每个 spot 的质量参差不齐。我们关心三个指标：

**① nCount_Spatial：每个 spot 捕获到的 UMI 总数（唯一分子标识符，unique molecular identifier）**。UMI 是建库时给每条 mRNA 分子贴的“身份证号”，所以 UMI 总数反映该 spot 的文库深度（测到了多少转录本）。太低说明这个 spot 可能没捕获到组织（背景 spot 或空 spot），太高则要警惕是不是技术异常。

**② nFeature_Spatial：每个 spot 检测到的基因种类数**。它与 UMI 数正相关（测得多自然基因多），但两者并不完全等价——一个 spot 如果只测到寥寥几个基因，说明它信息量太少，无法支撑后续分析。

**③ percent.mt：线粒体基因（mitochondrial genes）占总 UMI 的比例**。在单细胞分析里，线粒体比例高常提示细胞破裂、胞浆 RNA 流失，只剩线粒体 RNA 残留，代表坏死或濒死细胞。在 Visium 里同样适用：组织边缘或坏死区域的 spot 往往 percent.mt 偏高。这个指标需要我们自己算，Seurat 提供了现成函数：

```r
# 计算每个 spot 的线粒体基因比例（人样本线粒体基因以 MT- 开头）
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")

# 查看前 6 个 spot 的质控指标（行是 spot，列是指标）
head(obj[[]])
```

**运行结果解读**：`obj[[]]` 是 Seurat 对象的元数据（metadata）表，现在包含 `orig.ident`、`nCount_Spatial`、`nFeature_Spatial`、`percent.mt` 四列。`head()` 会打印前 6 行，你可以看到不同 spot 之间数字差异很大：有的 spot UMI 数上万，有的只有几百。这就是为什么必须做归一化和过滤。

> **为什么模式是 `"^MT-"`？** 人线粒体基因名以 `MT-` 开头（如 MT-ND1、MT-CO1），正则表达式 `^MT-` 表示“以 MT- 开头”。小鼠数据要用 `"^mt-"`，分析前先确认物种，别搞混。

## 8.3 阈值设定：先看分布，不拍脑袋

新手最常见的错误是直接抄别人论文的阈值（比如“基因数 > 200 就留下”）。不同组织、不同建库批次的数据分布天差地别，抄来的阈值可能把好 spot 全过滤掉，也可能留下大片垃圾。正确做法分三步：**看分布 → 按低分位数粗定 → 结合生物学常识微调**。

```r
# 第一步：画分布，直观感受数据长什么样
# VlnPlot 是 Seurat 的默认"体检"图，横轴是 spot 分组（未分组时全是 sample1）
VlnPlot(obj, features = c("nCount_Spatial", "nFeature_Spatial", "percent.mt"),
        ncol = 3, pt.size = 0)
```

```r
# 直方图补充：更精确地看数量分布（用 ggplot2 画 nCount 的直方图）
library(ggplot2)
ggplot(obj[[]], aes(x = nCount_Spatial)) +
  geom_histogram(bins = 80, fill = "#5B9BD5", color = "white") +
  labs(x = "UMI 总数 (nCount_Spatial)", y = "spot 数量", title = "文库深度分布")
```

**运行结果解读**：VlnPlot 会画出三个小提琴图，小提琴的“肚子”显示大多数 spot 落在哪个范围。直方图则更清楚：如果绝大多数 spot 的 UMI 数集中在 1000–8000，只有极少数 spot 接近 0，那接近 0 的那一小撮就是可疑对象。图 8-1 上排就是“先看分布”这一步的成品（见 8.7 节）。

**第二步：定阈值。** 一个实用经验是取分布的低分位数（比如 5%–10%）作为下限，再结合组织类型常识微调。对 Visium 人乳腺癌 demo，教程统一采用下面的阈值（与第 00-TECH 速查第 3 节一致）：

```r
# 过滤：同时满足 基因数>200 且 UMI 数>500 且 线粒体比例<25%
obj <- subset(obj, subset = nFeature_Spatial > 200 &
                     nCount_Spatial > 500 & percent.mt < 25)

# 过滤后还剩多少 spot？拿过滤前的数字对比一下
ncol(obj)
```

**运行结果解读**：`ncol(obj)` 返回剩余的 spot 数量。以 10x 乳腺癌 demo 为例，通常从约 4700 个 spot 过滤到约 4400–4600 个，丢掉 2%–6% 的异常 spot。注意三个条件用 `&`（且）连接，表示“任何一个指标不合格就淘汰”——这是更严格的策略；如果你只想剔除“极端异常”的 spot，可以放宽阈值，让更多 spot 留在下游分析里，之后再验证它们是否带来噪声。

**第三步：验证过滤效果。** 过滤后再画一次 VlnPlot，确认留下的 spot 分布干净、没有拖尾。这一步很多人跳过，但它其实是最快的自检。

> **阈值怎么微调？** 判据：过滤后 spot 数量不应骤降（骤降说明阈值太狠，好数据被误杀）；低质量 spot 应主要位于组织边缘或组织外（对照 tissue 图像看，见 8.6 节）。这两个判据比“抄论文数字”靠谱得多。8.8 节的练习让你亲手体验阈值变化的影响。

## 8.4 归一化：LogNormalize 还是 SCTransform？

即使过滤干净，不同 spot 的文库深度仍然不同：一个 spot 测到 8000 个 UMI，另一个只测到 2000 个，前者每个基因的计数天然更大。如果不做归一化，聚类时“文库深度差异”会伪装成“表达差异”，把 spot 按深度而不是按细胞类型分开。归一化的目标就是消除这种技术差异。

Seurat 提供两条主流路线：

| 方法 | 做法 | 优点 | 缺点 |
| --- | --- | --- | --- |
| LogNormalize | 每个基因的计数除以该 spot 总 UMI，乘 10000，再取 log1p（`log(x+1)`） | 简单、快、兼容性好 | 两步走（归一化后再单独找高变基因）；对文库深度校正不够精细 |
| SCTransform | 用正则化负二项回归同时建模 UMI 深度等技术因素与生物学异质性，一步完成归一化 + 高变基因选择 | 校正更彻底，聚类更干净，是当前推荐默认 | 明显更慢、更吃内存 |

**选型建议**：本书统一用 SCTransform（更优）；只有当你的机器跑不动（见第 06 章内存提示：Visium 单切片约 2–6 GB），或你明确需要与旧脚本的 LogNormalize 结果对齐时，才退回 `NormalizeData` + `FindVariableFeatures`。SCTransform 的输出直接替代“归一化 + 高变基因”两步，所以第 8.5 节里不需要再单独找高变基因：

```r
# SCTransform：一步完成归一化与高变基因选择（自动把结果存进 SCT assay）
obj <- SCTransform(obj, assay = "Spatial", verbose = FALSE)
```

**运行结果解读**：运行后 `obj` 里会出现一个名为 `SCT` 的新 assay。之后 PCA、聚类、UMAP 默认都用 `SCT` 的数据。`verbose = FALSE` 只是关掉进度刷屏，不影响结果。

> **对比记忆**：`LogNormalize` 相当于“除以总量再取对数”，`SCTransform` 相当于“对每个基因拟合一个带文库深度协变量的回归，取残差”。后者更精细，代价是计算量。

## 8.5 高变基因、PCA 与聚类

归一化之后，标准流程是：**降维 → 建图 → 聚类 → 再降维可视化**。

- **PCA（主成分分析，principal component analysis）**：把上万维的基因空间压缩成几十个“主成分”，每个主成分是基因的线性组合，捕捉最大的表达变异。我们取前 30 个主成分，因为前几个主成分包含绝大部分信息，后面基本是噪声。
- **FindNeighbors**：基于 PCA 空间为每个 spot 找邻居，构建 KNN 图（K 近邻图）。
- **FindClusters**：在 KNN 图上做 Louvain 社区检测，把图分成若干“表达相似”的群落——这就是聚类（cluster）。
- **RunUMAP**：把高维数据压到 2 维，方便肉眼观察聚类结构（UMAP 只用于可视化，聚类本身基于 PCA 空间）。

```r
# 完整降维聚类管道（与 00-TECH 第 3 节一致）
obj <- RunPCA(obj) %>%
  FindNeighbors(dims = 1:30) %>%
  FindClusters(resolution = 0.6) %>%
  RunUMAP(dims = 1:30)

# 看 UMAP 散点：每个点一个 spot，按聚类着色
DimPlot(obj, reduction = "umap", label = TRUE)
```

**运行结果解读**：`DimPlot` 会画出 UMAP 散点图，通常能看到几个抱团的群落，每个团一个颜色、一个 cluster 编号（如 0–7）。`resolution = 0.6` 控制聚类“颗粒度”：数值越大聚类越多越碎，越小越少越粗。0.6 是 Visium 默认推荐值，先跑一遍看结果再调。聚类编号本身没有生物学含义——0 号不等于“最好”，它只是社区检测的序号，**注释**是第 09 章的事。

## 8.6 空间可视化入门

单细胞分析里我们看 UMAP，空间转录组里更关键的是**把结果映射回组织切片**，看它长在哪个位置。Seurat 提供两个专属函数：

- `SpatialDimPlot`：按聚类/分组着色，画在组织图像上；
- `SpatialFeaturePlot`：按某个基因的表达量着色，画在组织图像上。

```r
# 聚类结果映射回组织切片：颜色相同的 spot 属于同一个 cluster
SpatialDimPlot(obj, label = TRUE)

# 单个基因的空间表达：看 EPCAM（上皮/肿瘤 marker）长在哪
SpatialFeaturePlot(obj, features = "EPCAM", pt.size.factor = 1.2)
```

**运行结果解读**：`SpatialDimPlot` 输出的是一张组织切片图，每个 spot 一个色块，能直观看到哪些 cluster 聚在组织中央、哪些贴着边缘——这是后面识别“空间生态位”（第 10 章）的雏形。`SpatialFeaturePlot` 用颜色深浅表示表达量：EPCAM 高表达的区域（深色）大概率是肿瘤细胞聚集区。看到这里你应该能体会到空间转录组的魅力：UMAP 丢掉的“位置信息”在这里回来了。

## 8.7 结果解读：QC 图与聚类图

现在把本章的产出图放在一起看。图 8-1 是 QC 三联图（模拟数据演示）：

![图 8-1 质控指标分布（模拟数据演示）。上排：按空间生态位的 UMI 数、基因数、线粒体比例小提琴图；下排：三个指标映射回组织切片的空间散点图，低质量 spot 用灰色标出。](figures/fig07_qc.png){width=85%}

**怎么读这张图**：上排三个小提琴图告诉你指标的整体分布——小提琴越“胖”的地方 spot 越多。如果线粒体比例的小提琴在 25% 以上还有一大截，说明坏死组织占比高，过滤阈值要收紧；如果小提琴整体贴着 0，说明数据干净，阈值可以放宽。下排空间散点图是“位置版”的指标：正常组织里，低 UMI、低基因、高线粒体的 spot（灰色）应该零星散布在组织边缘或组织外——如果它们成片出现在组织中央，那更可能是切片或捕获环节出了问题，而不是单纯的质量差。**注意：本图为模拟数据演示，仅用于展示“真实分析会长什么样”，不代表任何真实样本的结论。**

图 8-2 是降维聚类结果（模拟数据演示）：

![图 8-2 聚类结果（模拟数据演示）。左：UMAP 降维散点图，共 8 个 cluster；右：同一聚类结果映射回组织切片的空间分布。](figures/fig08_clustering.png){width=85%}

**怎么读这张图**：左图是“无位置视角”——8 个 cluster 在 UMAP 上各自成团，说明它们的表达谱确实不同；右图是“有位置视角”——可以看到某些 cluster 的空间位置高度集中（比如一大片红色占据切片中央），这通常对应组织结构里的某一类区域。**注意**：聚类是纯数据驱动（只看表达，不看位置），它和后面 BayesSpace 的“空间聚类”（第 10 章，把位置信息也纳入）是两回事，别混淆。此时我们还不知道每个 cluster 是什么细胞——注释是第 09 章的任务。

## 8.8 动手练习：调整阈值，观察 spot 数量变化

过滤阈值没有标准答案，最好的老师是自己动手。做下面这个实验：

```r
# 练习：先记住过滤前的 spot 数量
n_before <- ncol(obj)

# 尝试三组阈值，分别记录剩余的 spot 数
thr1 <- subset(obj, subset = nFeature_Spatial > 200 &
                      nCount_Spatial > 500 & percent.mt < 25)
thr2 <- subset(obj, subset = nFeature_Spatial > 500 &
                      nCount_Spatial > 1000 & percent.mt < 10)
thr3 <- subset(obj, subset = nFeature_Spatial > 100 &
                      nCount_Spatial > 200 & percent.mt < 40)

# 打印对比结果（用英文列名，避免中文列名在部分环境被改写）
data.frame(
  plan    = c("过滤前", "宽松(thr3)", "教程默认", "严格(thr2)"),
  n_spots = c(n_before, ncol(thr3), ncol(thr1), ncol(thr2))
)
```

**运行结果解读**：你会看到“严格”方案可能丢掉 30% 以上的 spot，“宽松”方案几乎不丢。接着分别对 `thr2` 和 `thr3` 跑一次 8.5 节的聚类管道，对比 UMAP 和 `SpatialDimPlot`：如果“宽松”方案出现一个主要由低质量 spot 组成的小团，说明阈值太松；如果“严格”方案把组织边缘的有用信息也删了，说明阈值太狠。**练习的结论不是“哪个阈值对”，而是让你建立“阈值与下游结果联动”的感觉**——这正是做科研和跑脚本的区别。

## 本章小结

本章完成了空间转录组分析的第一块地基：用 `Load10X_Spatial` 读入数据，计算并理解三个 QC 指标（UMI 总数、基因数、线粒体比例），按“先看分布、再定阈值”的原则过滤低质量 spot，用 SCTransform 归一化，跑通 PCA → FindNeighbors → FindClusters → RunUMAP 的完整管道，最后学会把结果映射回组织切片。现在你手里有一个干净、归一化、带聚类的 Seurat 对象，下一步（第 09 章）就是回答最关键的问题：**每个 cluster、每个 spot 里到底是什么细胞？**

## 练习与思考

1. **概念题**：为什么 `nCount_Spatial` 低的 spot 要过滤掉？如果不过滤，它会对后续聚类造成什么影响？
2. **概念题**：线粒体比例高的 spot 在生物学上通常意味着什么？这个指标在 Visium 和单细胞（scRNA-seq）里的解读有什么异同？
3. **动手题**：按 8.8 节的三组阈值分别跑一遍过滤，记录 spot 数量变化，并解释“严格”与“宽松”两种策略各自的利弊。
4. **动手题**：对过滤后的对象画 `SpatialFeaturePlot` 展示基因 `EPCAM` 和 `COL1A1`，观察两者空间分布是否不同，想想这提示了什么（提示：与乳腺癌的组织结构有关）。
5. **思考题**：SCTransform 相比 LogNormalize 多做了什么？在什么情况下你会选择更慢的 SCTransform？（提示：从聚类稳定性角度思考。）
