📚 课程目录第3章数据预处理

数据预处理：读入数据、过滤低质量Spot、标准化、Log转换

好，咱们进入正题。数据预处理这一步，说白了就是给原始计数矩阵“洗个澡”，把脏东西去掉，再整理成适合后续分析的样子。很多新手一上来就急着跑差异分析，结果跑出来的都是噪声——我自己早期就吃过这个亏，后来才明白，预处理做不好，后面全是白费功夫。

读入数据：从文件夹到内存对象

空间转录组的数据格式一堆，但最常用的还是10X Visium那一套。你拿到一个样本目录，里面通常有三个文件：raw_feature_bc_matrix.h5（或 mtx）、tissue_positions_list.csv、scalefactors_json.json。我习惯用 Seurat 的 Load10X_Spatial() 或者 SpatialExperiment 的 read10xVisium() 来读。

下面是一个读入示例，用 SpatialExperiment 包：

library(SpatialExperiment)
library(Matrix)

spe <- read10xVisium(
  samples = "path/to/your/sample",
  sample_id = "sample1",
  type = "sparse",       # 原始稀疏矩阵
  data = "raw",          # 读取原始计数
  images = "lowres",     # 低分辨率组织图像
  load = TRUE
)

读进来以后，spe 对象里就包含了计数矩阵、空间坐标、图像数据。这时候最好看一眼总UMI数、基因数、Spot数，心里有个底。我曾经见过一个样本，读进来总UMI只有几千——后来发现是压缩包没解压完全，白白浪费了半天排查。

提示：如果使用 Seurat，记得加载 SeuratData 并指定 assay = "Spatial"。两种工具底层数据结构不同，但核心逻辑一致。

过滤低质量Spot：宁缺毋滥

空间转录组里，一个Spot对应一个直径约55微米的圆点，里面可能包含几个到十几个细胞。如果这个Spot没捕获到多少mRNA，或者大部分来自线粒体，那它要么是组织边缘的空泡，要么是死细胞区。留着它只会增加噪声。

我一般看三个指标：

总UMI计数（nCount_Spatial）：通常取几千到几万，过低说明捕获失败，过高可能来自组织折叠。
检测到的基因数（nFeature_Spatial）：和UMI数正相关，但更敏感，过低直接丢掉。
线粒体基因比例（percent.mt）：高比例（比如 >20%）表示细胞膜破损，RNA泄漏。

举个例子，我手头一个乳腺癌样本，初始Spot有5000个，但过滤完只剩4200个——那800个都是组织边缘的脏点。下面是过滤代码：

# 计算线粒体基因比例（基因名以MT-开头）
spe <- addPerCellQCMetrics(spe, subsets = list(Mt = grep("^MT-", rownames(spe))))

# 设定阈值，这里按经验取
qc_lib_size <- spe$sum > 500
qc_detected <- spe$detected > 300
qc_mt <- spe$subsets_Mt_percent < 20

# 合并筛选
discard <- !(qc_lib_size & qc_detected & qc_mt)
spe_filtered <- spe[, !discard]

注意：阈值不能死搬硬套。不同组织、不同平台（10X v1 vs v2）的分布差异很大。我建议先画小提琴图看分布，再手动定边界。有一次我直接用了默认值，结果把真正的腺体区域给过滤掉了——气得我重新跑了一整天。

标准化：让数据可比

过滤完，每个Spot的测序深度不一样——深的几十万UMI，浅的可能只有几百。不标准化的话，差异分析全是深度差异，而不是真正的生物学差异。常用方法有两种：

方法	核心思想	适用场景	缺点
LogNormalize (Seurat / SCANPY)	每个Spot的UMI除以总UMI，乘以缩放因子，再log1p	快速、通用	对高表达基因压制强，可能丢失弱信号
SCTransform	基于负二项回归，同时标准化并方差稳定	数据异方差严重时更稳健	计算慢，内存大

我个人偏好 SCTransform。为啥？因为它把标准化和方差稳定合并了，还自动处理了技术噪声。下面是用Seurat的做法：

# 如果用的是Seurat
seurat_obj <- SCTransform(seurat_obj, assay = "Spatial", verbose = FALSE)

如果数据量很大（比如超过10万Spot），SCTransform会跑几分钟甚至更久。这时我会退回到 NormalizeData() + FindVariableFeatures() 的经典流程，虽然粗糙一点，但够用。

Log转换：让数据更接近正态

标准化之后，数据仍然偏右。计数数据天然是泊松/负二项分布，做log转换能拉平方差，让高表达和低表达基因的差异不那么极端。常用的就是 log1p：

# 手算log转换
log_counts <- log1p(normalized_counts)

在Seurat里，如果你用了 NormalizeData()，默认就已经做了log1p。如果你用SCTransform，它会在内部做类似的处理，但更精细。不过要注意——log1p之前的标准化因子选择会影响结果。我见过有人把scale factor设成10000，结果高表达基因被压得太狠，后续聚类全混在一起。

  推荐参数：对于10X Visium数据，scale factor取10000比较稳妥，这是Seurat的默认值。如果你用SCTransform，它的内部参数调校更好，几乎不用动。

整体流程一览

我把这四步串起来，画了个流程图，拿走不谢：

从读取原始数据开始，经过过滤、标准化、Log转换，最终得到干净、可比较的表达矩阵。后面的批次效应消除、聚类、差异分析全都要基于这一步的输出。

小技巧：每次预处理完，我都习惯用 pairs() 画一下几个QC指标的散点图矩阵，看看过滤阈值是否合理。如果看到一刀切把大片连续区域切掉了，那就得放宽标准，或者改用更软性的迭代过滤。

嗯，到这里咱们的数据预处理就完成了。数据干净了，后面才有底气去处理批次效应。记住：投喂给算法的数据质量，决定了你能挖出多少真东西。