2. Seurat 流程基础：数据加载与预处理
单细胞分析的第一步，说白了就是把测序仪吐出来的原始数据，变成你能在 R 里摆弄的 Seurat 对象。这一步看似简单，但坑不少。我个人习惯把整个预处理流程拆成 加载 → 创建 → 质检 → 过滤 → 标准化 → 找可变基因 → 缩放 七个步骤，每一步都有它存在的道理。
核心思路：预处理的目标是保留真实生物学信号，同时剔除技术噪声和死细胞/双细胞。你后面做的聚类、差异分析、降维，全部建立在这个“干净”的数据之上。
2.1 数据加载：从 10X 文件夹到矩阵
大部分单细胞数据来自 10X Genomics 平台。你拿到手的通常是一个文件夹，里面有 barcodes.tsv、features.tsv（或 genes.tsv）、matrix.mtx。Seurat 提供了 Read10X() 函数，一句代码就能读进来。
library(Seurat)
# 假设 data_dir 是包含三个文件的路径
counts <- Read10X(data.dir = "filtered_feature_bc_matrix/")注意：10X 现在默认输出 filtered_feature_bc_matrix（已过滤空液滴）和 raw_feature_bc_matrix（含空液滴）。我建议除非你要做空液滴分析，否则直接用 filtered 版本，能省很多过滤功夫。不过，有些公共数据集的文件夹名不规范，我曾经碰到过用 outs/ 目录下的 filtered_gene_bc_matrices 这种情况，你自己留个心眼。
如果是多个样本合并分析（比如你做完 Harmony 调参），可以这样先分别读取每个样本的 counts list，后续再用 MergeSeurat 或 CreateSeuratObject 的 merge 参数合并。
2.2 创建 Seurat 对象：把矩阵包起来
拿到矩阵后，用 CreateSeuratObject() 把它包装成一个 Seurat 对象。这时你可以指定 project 名称和 最小细胞数/最小基因数 的初筛参数——不过我个人更倾向于后面再统一用 QC 指标过滤，这里可以设得很宽松。
seurat_obj <- CreateSeuratObject(
counts = counts,
project = "my_project",
min.cells = 3, # 至少在3个细胞中表达
min.features = 200 # 至少检测到200个基因
)为什么设这俩值？min.cells = 3 是为了去除只在极个别细胞里出现的基因，那些通常是扩增噪声；min.features = 200 是为了剔除空液滴或破碎细胞。不过你看，这只是粗筛。
经验之谈：我刚开始做单细胞时，把 min.features 设成了 2000，结果 80% 的细胞都被剔除了……后来才发现数据来自冷冻切片，人家活性本来就低。所以，这个值要结合你的项目和物种来定，不要照搬教程。
2.3 质量控制（QC）：揪出坏细胞
QC 是预处理里最“玄学”也最关键的环节。Seurat 对象会自动计算三个核心指标：
- nCount_RNA：每个细胞的总 UMI 数（文库大小）
- nFeature_RNA：每个细胞检测到的基因数（基因丰富度）
- percent.mt：线粒体基因 UMI 占比（细胞健康状况指标）
你要手动算一下线粒体比例，因为 Seurat 不自动算这个。标准做法是找到所有以 MT- 开头的基因（人源），然后求和占比。
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(
seurat_obj, pattern = "^MT-"
)| 指标 | 正常范围（scRNA-seq 常用） | 可能的问题 |
|---|---|---|
| nCount_RNA | 500 – 50,000（因平台而异） | 过低：空液滴/死细胞；过高：双细胞 |
| nFeature_RNA | 200 – 5,000 | 过低：无活性细胞；过高：双细胞或异常细胞 |
| percent.mt | < 20%（有些组织 < 10%） | 过高：细胞膜破损，线粒体内容泄露 |
你可以用 VlnPlot 和 FeatureScatter 看一眼分布。但记住，图表是辅助，真实过滤阈值要靠你的知识和实验背景。我曾经在一个肺癌 T 细胞数据上，发现 CD8+ 效应 T 细胞的 nCount_RNA 天然就低于 B 细胞——如果我用统一阈值，就把它们全杀了。
避坑指南：不要同时用 nFeature 和 nCount 的上下限做“暴力过滤”。这两个指标高度正相关（UMI 越多，检测到的基因越多），你用一个就够了，另一个用来辅助判断双细胞。我一般用 nFeature 做主过滤，配合 percent.mt 做二次筛选。
2.4 过滤细胞：动刀前的最后确认
基于上面的 QC 指标，我们写几个过滤条件。注意你应该先画图、肉眼确认异常值，再写固定数值。常见写法：
# 示例过滤：保留 nFeature 在 200~5000 且 percent.mt < 25% 的细胞
seurat_obj <- subset(seurat_obj,
subset = nFeature_RNA > 200 & nFeature_RNA < 5000 &
percent.mt < 25)你是不是想问：为什么不设 nCount 的下限？我个人的方案是：nFeature 已经能区分空液滴，而 nCount 可能受批次影响大，所以我不爱用它做主过滤条件。不过如果你的数据中有明显的“高 UMI 低基因数”的异常细胞（比如血小板在测序中破裂），那就需要加一条 nCount_RNA < 30000 之类的条件。
2.5 标准化：消除文库大小差异
过滤完细胞，接下来要处理不同细胞间测序深度不一致的问题。Seurat 默认的标准化方法是 LogNormalize：用总 UMI 做缩放因子，再取 log1p。
seurat_obj <- NormalizeData(seurat_obj,
normalization.method = "LogNormalize",
scale.factor = 10000)这里 scale.factor = 10000 的意思是：每个细胞的所有基因先除以该细胞的总 UMI，再乘以 10000，最后做 log(x + 1)。得到的数值就是 logCPM（每百万计数对数）的一个变种。为什么取 log？因为原始 UMI 是偏态分布，取 log 后更接近正态，方便后续做线性降维。
要点：标准化的目的是让不同细胞之间可比，而不是修正批次效应（那是 Harmony 的工作）。别在这里就急着用 SCTransform，那是另一条路线。对于 Harmony 调参教程，我们先走经典的 LogNormalize 路线。
2.6 识别高变基因（HVG）：只挑最“有故事”的基因
全基因组有 2 万多个，很多基因在所有细胞里表达量都差不多，对区分细胞类型没贡献。所以我们只挑那些在细胞间表达量变异最大的基因——叫 高变基因（Highly Variable Genes, HVG）。
seurat_obj <- FindVariableFeatures(seurat_obj,
selection.method = "vst",
nfeatures = 2000)我一般选 2000 个，这是 Seurat 的默认值，也是社区经验值。不过你想想看：如果你的样本包含多种稀有细胞类型（比如神经元的亚型），可能需要上调到 3000。反之如果只是两种细胞混合，1000 就够。我在一个包含 10 种免疫细胞的数据集上试过，2000 和 3000 的结果几乎没区别，但 500 会丢失稀有亚群。
selection.method = "vst" 是目前最推荐的方法，它使用方差稳定变换来挑选基因。备选还有 "mean.var.plot" 和 "dispersion"，但我建议你就用 VST，又快又好。
2.7 缩放数据：为 PCA 做准备
最后一步是 ScaleData。它对每个基因做 z-score 标准化（中心化 + 缩放），让所有基因的均值为 0，方差为 1。这一步对 PCA 和后续的 Harmony 调整至关重要——否则高表达基因会主导主成分。
seurat_obj <- ScaleData(seurat_obj,
features = rownames(seurat_obj)等等！注意这里，默认 ScaleData 会对 全部基因 做缩放，但当你后续只打算用 HVG 跑 PCA 时，完全可以只对 变量基因 缩放，省内存又省时间。
seurat_obj <- ScaleData(seurat_obj,
features = VariableFeatures(seurat_obj))有没有副作用？其实没有。PCA 只输入变量基因，非变量基因缩放后也不会被用到。所以我个人一直用这种写法。不过要小心：如果你的下游想用 Heatmap 绘制全部基因的表达模式，就必须对所有基因缩放，这视情况而定。
一个小窍门：如果数据量很大（几万细胞以上），可以考虑用 ScaleData(..., do.center = TRUE, do.scale = TRUE, vars.to.regress = "percent.mt") 来回归掉线粒体比例这个技术变量。不过 在 Harmony 流水线里，我不推荐在预处理阶段回归太多变量，因为 Harmony 本身就处理批次效应，你提前回归可能干扰它的效果。把回归留着做“后处理”更安全。
好，到这一步，你的 Seurat 对象已经完成了预处理的所有标准步骤。接下来就是主成分分析（PCA）和聚类了。不过咱们得先把这一步流程理清楚，不然后面调参数容易迷失方向。
下面这张图把整个预处理流程串起来，你可以边看边核对代码里的每一步。