Seurat 基础：创建对象 · 质量控制

Seurat 对象是这个工具链里的核心容器——说白了，所有单细胞数据最后都塞进这个结构里。我这些年接手的项目，十有八九第一步就是检查数据能不能正确塞进 Seurat 对象。这一步搞砸了，后面所有分析都白搭。

下面拆开讲：怎么创建对象，怎么看质控指标，怎么避开我踩过的坑。

1. 从表达矩阵到 Seurat 对象

原始数据通常是个计数矩阵：行是基因，列是细胞。我们要把它封装成 Seurat 对象。

最常用的函数是 CreateSeuratObject。下面是个标准流程：

library(Seurat)

# 假设 data 是行基因×列细胞的稀疏矩阵
seurat_obj <- CreateSeuratObject(
  counts = data,
  project = "my_project",
  min.cells = 3,    # 至少在3个细胞中表达的基因才保留
  min.features = 200 # 至少表达200个基因的细胞才保留
)

参数 min.cells 和 min.features 是预过滤。我一般设得很宽松，后面会做更精细的过滤。曾经有个项目我贪快把 min.features 设成 500，结果发现丢弃了大量真实细胞——因为那个样本本身是低转录活性的中性粒细胞。所以预过滤门槛宁可低一些。

2. 质量控制的第一步：理解元数据

创建对象后，Seurat 会自动计算两个字段：nCount_RNA（每个细胞的总 UMI 数）和 nFeature_RNA（每个细胞检测到的基因数）。它们存在对象的元数据（meta.data）里。

查看元数据：

head(seurat_obj@meta.data)

细胞ID	nCount_RNA	nFeature_RNA
AAACCCAAGCATACT-1	2974	1128
AAACCCACACCGTAG-1	1876	862
AAACCCAGTCAGATG-1	4512	1653

这两个字段就是质控的原始素材。你想想看，如果一个细胞只测到几十个基因，大概率是空液滴或者死细胞。

3. 计算线粒体比例——这步最常出错

线粒体基因比例高（通常 >20%）意味着细胞膜破损，胞质 RNA 流失，只剩线粒体 RNA 被包在双层膜里。计算方式是用正则匹配：

seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")

这里我提醒一句：人源数据线粒体基因前缀是 ^MT-，小鼠是 ^mt-。我曾经在一个项目里直接复制人性代码跑小鼠数据，结果线粒体比例全是 0，排查了半天才发现是大小写问题。务必检查物种的前缀规则。

4. 可视化质控指标——快速定位异常

用 VlnPlot 看分布最直观：

VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

一般效果像这样（示意图，非真实数据）：

5. 设定过滤阈值——经验与统计结合

过滤没有绝对标准，但我总结了一套经验值：

指标	常见下限	常见上限	说明
nFeature_RNA	200 – 500	5000 – 8000	太低→空液滴；太高→可能是双细胞
nCount_RNA	500	20000 – 50000	双细胞往往总 UMI 数异常高
percent.mt	–	15% – 25%	超过 25% 大概率是破损细胞

实际操作中我喜欢先用 VlnPlot 看分布，再结合 FeatureScatter 看相关性：

FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "percent.mt")
FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

我曾经在一个肿瘤样本里发现 nCount_RNA 和 percent.mt 有明显正相关——那些 UMI 极高的细胞全是线粒体基因，最后确认是红细胞污染。如果不看散点图，单纯用阈值过滤会漏掉这种偏置。

6. 执行过滤——选择 subset 还是直接索引？

两种方式都行，我推荐用 subset 函数，可读性强：

seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 20)

注意：过滤条件里不能有 | 或复杂的逻辑嵌套？其实可以，但建议拆成多步防止误删。

另一种直接操作元数据：

cells_keep <- which(seurat_obj$nFeature_RNA > 200 & seurat_obj$nFeature_RNA < 5000 & seurat_obj$percent.mt < 20)
seurat_obj <- seurat_obj[, cells_keep]

7. 过滤后检查——别跳步骤

过滤完再用 VlnPlot 和 FeatureScatter 看一次，确认异常群体被移除。

VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

这一步看似多余，但我吃过亏——有一次过滤阈值设得太宽，实际只删了 2% 的细胞，可视化才发现分布根本没变。那意味着阈值本身无效，需要重新查原始数据的分布。

8. 保存预处理的 Seurat 对象

到这一步我就用 saveRDS 把干净的对象存起来，作为后续分析的起点。

saveRDS(seurat_obj, "data/seurat_raw_qc.rds")

attr(seurat_obj, "qc_params") <- list(
  min_features = 200,
  max_features = 5000,
  max_mt = 20
)
saveRDS(seurat_obj, "data/seurat_qc_params.rds")

小结

创建 Seurat 对象和质控是单细胞分析里最琐碎但最关键的一步。我见过太多项目因为线粒体比例算错、阈值设错导致后续聚类全歪。把流程固定下来，用 Snakemake 自动化——这就是这个课程后面要做的。不过先别急，先把这一章的手工操作练熟，后面集成到管道时才不容易翻车。