📚 课程目录第2章单细胞数据加载与预处理

第二章：单细胞数据加载与预处理 — 从原始数据到干净矩阵

做单细胞分析，第一件事就是拿到原始数据。我见过不少新手直接拿 Count 矩阵就跑下游分析，结果质控没做好，后面聚类怎么调都不对。这一章咱们就把数据加载和预处理这块地基打扎实。

说白了，预处理的四个核心步骤就是：读取 → 质控 → 标准化 → 归一化。每一步都有坑，我来一个一个讲。

2.1 读取数据：建立 Seurat 对象

Seurat 支持多种数据格式。最常用的是 10X Genomics 的 Cell Ranger 输出，也就是三个文件：barcodes.tsv、features.tsv（旧版叫 genes.tsv）、matrix.mtx。我建议直接从这仨文件创建对象，因为最不容易出错。

标准做法：用 Read10X() 读取，再传给 CreateSeuratObject()。

# 读取 10X 数据
data_dir <- "path/to/filtered_feature_bc_matrix/"
counts <- Read10X(data.dir = data_dir)

# 创建 Seurat 对象，只保留至少在 3 个细胞中表达的基因
seurat_obj <- CreateSeuratObject(
  counts = counts,
  project = "pbmc_tutorial",
  min.cells = 3,
  min.features = 200
)

min.cells 和 min.features 是两个常用过滤参数。前者过滤掉只在极少数细胞中表达的基因（噪音），后者过滤掉基因数过少的细胞（可能是空液滴或死细胞）。

我的习惯：先设 min.cells = 3，等质控完成后再根据情况决定是否收紧。不要一开始就过滤得太狠，你永远不知道那些低表达的基因是不是关键标记。

# 查看对象基本信息
seurat_obj

如果手头是 H5 格式文件（比如从 GEO 下载的），就用 Read10X_h5()。

counts <- Read10X_h5("filtered_feature_bc_matrix.h5")
seurat_obj <- CreateSeuratObject(counts = counts, project = "h5_example")

嗯，这里要注意一个容易忽略的点：CreateSeuratObject 默认把基因名放在行名、细胞 barcode 放在列名。后续所有操作都基于这个矩阵结构，所以一开始就要确保行列名没有乱码或重复。

2.2 质量控制（QC）：筛掉坏细胞

质控的目标是去掉三样东西：空液滴（基因数太少）、双细胞/多细胞（基因数异常多）、死细胞（线粒体基因比例高）。

Seurat 提供了 PercentageFeatureSet() 来计算线粒体基因占比。我建议再顺手算一个红细胞相关基因的占比，用于排除红细胞污染。

核心指标：nFeature_RNA（基因数）、nCount_RNA（UMI 总数）、percent.mt（线粒体占比）。

# 计算线粒体基因比例
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")

# 可选：计算核糖体基因比例
seurat_obj[["percent.ribo"]] <- PercentageFeatureSet(seurat_obj, pattern = "^RP[SL]")

# 查看质控指标
head(seurat_obj@meta.data)

拿到指标后，怎么确定阈值？我个人的做法是先画小提琴图看看分布。

# 可视化三个关键指标
VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        ncol = 3, pt.size = 0.01)

你会看到大部分细胞集中在某个范围，但也有一些离群值。比如 nFeature_RNA 超过 5000 的很可能是双细胞，percent.mt 超过 20% 的基本是死细胞。

我曾经在一个肝组织样本上吃过亏：默认用了 percent.mt < 25，结果下游发现一群 hepatocyte 死活聚类不对。后来才发现肝细胞的线粒体比例天生就高，当时把正常细胞也过滤掉了。所以阈值一定要根据组织类型和实验条件灵活调整。

避坑指南：不要复制粘贴别人的阈值。每种组织、每个建库批次都可能不同。建议先画图看一下分布，再定阈值。

# 手动过滤：以常见 PBMC 数据为例
seurat_obj <- subset(seurat_obj,
                     subset = nFeature_RNA > 200 & nFeature_RNA < 5000 &
                              nCount_RNA > 500 & nCount_RNA < 20000 &
                              percent.mt < 15)

2.3 标准化：让细胞之间可比

标准化是为了消除细胞间测序深度的差异。说白了，一个细胞测到 10 万 UMI，另一个只测到 1 万，直接用原始 Count 比较肯定有问题。

Seurat 默认用 NormalizeData()，方法选 "LogNormalize"：先对每个细胞的 UMI 总数做缩放，再取自然对数。

# Seurat 默认标准化
seurat_obj <- NormalizeData(seurat_obj,
                            normalization.method = "LogNormalize",
                            scale.factor = 10000)

公式很简单：log1p( (gene_count / total_umi_per_cell) * 10000 )。log1p 加 1 是为了避免 log(0) 导致负无穷。

你问我为什么默认是 10000？ 这个值相当于假设每个细胞有 1 万条 UMI，是个经验值。如果你的数据平均 UMI 深度是 5000，可以改成 5000，结果差别不大。

2.4 鉴定高变基因：抓住最值得关注的基因

全基因组几万个基因，不是每个都对聚类有用。高变基因（HVG）就是那些在不同细胞间表达差异最大的基因。我一般取前 2000 个 HVG 做下游分析，这个数量基本够用。

Seurat 的 FindVariableFeatures() 提供了三种选择方法："vst"（推荐）、"mean.var.plot"、"dispersion"。

# 选择前 2000 个高变基因，采用 vst 方法
seurat_obj <- FindVariableFeatures(seurat_obj,
                                   selection.method = "vst",
                                   nfeatures = 2000)

# 查看选中的基因
head(VariableFeatures(seurat_obj))

做这一步之前别忘了先做 NormalizeData，否则原始 Count 的方差会被测序深度噪声严重干扰。

重要提醒：NormalizeData 之后，FindVariableFeatures 会认为你已经完成了标准化。如果你先做了 NormalizeData 又做了其他操作，顺序要捋清楚。

2.5 归一化（Scale Data）：给 PCA 做准备

归一化和标准化不是一回事。标准化解决的是细胞间深度差异，归一化解决的是基因表达量量纲不统一的问题。PCA 对量纲敏感，所以这一步给每个基因做 Z-score 变换，让所有基因的均值 = 0，方差 = 1。

# 对所有基因做归一化
seurat_obj <- ScaleData(seurat_obj, features = rownames(seurat_obj))

# 也可以只对高变基因做归一化（省内存）
seurat_obj <- ScaleData(seurat_obj,
                        features = VariableFeatures(seurat_obj))

注意：ScaleData 默认对全部基因做。如果数据集有 3 万个基因，这一步会非常慢。我通常只传入 VariableFeatures，毕竟聚类主要靠高变基因。

ScaleData 里还有一个 vars.to.regress 参数，可以用来去除批次效应或细胞周期影响的信号。不过这一般放到后面批次校正部分再做，这里先记住有这么一个功能就行。

核心逻辑总结：

读取数据 → 创建 Seurat 对象
QC 过滤 → 根据 nFeature、percent.mt 等指标去掉坏细胞
标准化（Normalize）→ 按细胞 UMI 总数缩放 + 对数变换
鉴定 HVG → 挑出最值得关注的前 2000 个基因
归一化（Scale）→ Z-score 变换，让基因间量纲统一

2.6 完整流程串联：一步到位

实际项目中我不会一句一句敲。写一个管道函数更高效。下面是常见 PBMC 样本的完整流程：

# ==== 完整数据预处理管道 ====
run_preprocessing <- function(data_dir, project_name = "project") {
  # 1. 读取
  counts <- Read10X(data.dir = data_dir)
  
  # 2. 创建对象 + 基础过滤
  obj <- CreateSeuratObject(
    counts = counts,
    project = project_name,
    min.cells = 3,
    min.features = 200
  )
  
  # 3. QC 指标
  obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
  
  # 4. 画出 QC 小提琴图（可选）
  # VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"))
  
  # 5. 过滤
  obj <- subset(obj,
                subset = nFeature_RNA > 200 & nFeature_RNA < 5000 &
                         nCount_RNA > 500 & nCount_RNA < 20000 &
                         percent.mt < 15)
  
  # 6. 标准化
  obj <- NormalizeData(obj, normalization.method = "LogNormalize",
                       scale.factor = 10000)
  
  # 7. 高变基因
  obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 2000)
  
  # 8. 归一化
  obj <- ScaleData(obj, features = VariableFeatures(obj))
  
  return(obj)
}

# 使用示例
seurat_clean <- run_preprocessing("path/to/filtered_feature_bc_matrix/",
                                  project_name = "PBMC_clean")

2.7 本章知识结构图

我画了一张图来梳理这五个主要步骤的依赖关系。从原始数据到可以做 PCA 的标准化矩阵，每一步都不能跳。