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

单细胞数据预处理：原始数据格式与质量控制

好，咱们正式进入单细胞数据分析的第一步——数据预处理。说实话，这一步做得怎么样，直接决定了后面所有分析的成败。我见过太多人拿到数据就急着跑聚类、找Marker基因，结果QC没做好，后面全是垃圾。

数据预处理的核心就两件事：读懂原始数据格式，然后做好质量控制（QC）。这两件事搞明白了，后面的分析才有意义。

一、原始数据格式：你要读懂的"三件套"

单细胞测序最常用的平台是10x Genomics。你从测序公司拿到的数据，通常不是直接能用的表达矩阵。我刚开始接触时也纳闷，怎么给了一堆文件？

其实10x的标准输出是三个文件，缺一不可：

文件	内容	说明
barcodes.tsv	所有细胞条形码列表	每个条形码代表一个微滴/细胞
features.tsv	基因ID和基因名	Ensembl ID + 基因符号
matrix.mtx	稀疏表达矩阵	Market Exchange Format，只存非零值

这三个文件必须放在同一个目录下。读取的时候，大多数工具（比如Seurat、Scanpy）会自动识别。我个人习惯用Seurat的Read10X()函数，它能一次性把三个文件都读进来。

# R语言 - 读取10x数据
library(Seurat)

# 假设你的数据放在"raw_data/filtered_feature_bc_matrix/"目录下
data_dir <- "raw_data/filtered_feature_bc_matrix/"

# Read10X会自动找三个文件
counts <- Read10X(data.dir = data_dir)

# 创建Seurat对象
seurat_obj <- CreateSeuratObject(counts = counts,
                                 project = "my_project",
                                 min.cells = 3,
                                 min.features = 200)

为什么用稀疏矩阵？ 单细胞数据里，大部分基因在单个细胞中表达量为0。假设你有20000个基因、5000个细胞，全量矩阵有1亿个值，但其中90%以上都是0。用稀疏格式存储，内存占用能减少一个数量级。

小技巧： 如果你的数据是原始的raw_feature_bc_matrix，那里面包含了所有微滴的数据，包括空液滴。我用filtered_feature_bc_matrix，因为10x已经帮你过滤掉了明显是空的液滴。但说实话，10x的过滤标准很宽松，我们自己还得再做一轮QC。

二、质量控制（QC）：你为什么不能跳过这一步？

QC的本质就是筛选出高质量的真实细胞，剔除那些不合格的。不合格的包括哪些？

死细胞或濒死细胞：线粒体基因比例高
空液滴或双细胞：基因数异常低或异常高
技术噪音：UMI总数过低或过高

我自己总结了一套QC指标，分为必查项和参考项：

三个核心QC指标

指标	意义	常见阈值	过滤原因
nFeature_RNA	每个细胞检测到的基因数	200 < x < 2500/5000	过低→空液滴/死细胞；过高→双细胞
nCount_RNA	每个细胞的UMI总数	一般与nFeature正相关	过低→捕获效率差；过高→双细胞
percent.mt	线粒体基因UMI占比	< 5% ~ 20%	过高→细胞膜破损，胞质RNA泄漏

这里要注意，阈值不是死的。我记得有一次做人胰腺数据，线粒体比例天然就高，用到20%的阈值才合适。如果是脑组织，5%都不一定能过。所以我的建议是：先看分布，再定阈值。

计算QC指标的代码（R + Seurat）

# 计算线粒体基因比例
# 基因名以"MT-"开头的就是线粒体基因
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj,
                                                   pattern = "^MT-")

# 查看QC统计
head(seurat_obj@meta.data)

# 可视化QC
VlnPlot(seurat_obj,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        ncol = 3)

我踩过的坑： 曾经有一批数据，我直接用默认的200~2500过滤基因数，结果把一类真正的低表达细胞（比如中性粒细胞）全过滤掉了。后来我才意识到，不同细胞类型的基因数差异巨大。所以呢？先画图看看分布，再决定阈值。

三、实际操作：QC过滤怎么做？

好，理论说完了，咱们直接上代码。这就是我每次跑单细胞数据的标准流程：

# 设定阈值
min_features <- 200
max_features <- 2500
max_mt <- 5  # 百分比

# 执行过滤
seurat_obj_filtered <- subset(seurat_obj,
                              subset = nFeature_RNA > min_features &
                                       nFeature_RNA < max_features &
                                       percent.mt < max_mt)

# 检查过滤后的细胞数
dim(seurat_obj_filtered)

过滤后，我通常会再跑一次QC图，确保过滤效果合理。如果一下子砍掉了超过30%的细胞，我会回头检查阈值是不是太严格了。

四、双细胞检测：一个容易被忽略的环节

我认为单细胞QC中最容易被跳过、但也后果最严重的一步就是双细胞检测。两个细胞挤在一个微滴里，测出来的表达谱就是一半一半，这种数据混进去，后面聚类时会凭空多出一个"杂交"群来。

常用的工具是DoubletFinder（R包）或scrublet（Python）。我的做法是在过滤完线粒体和低质量细胞后，先对数据进行初步归一化，再跑双细胞检测。

# 使用DoubletFinder检测双细胞
library(DoubletFinder)

# 先做标准化和PCA
seurat_obj_filtered <- NormalizeData(seurat_obj_filtered)
seurat_obj_filtered <- FindVariableFeatures(seurat_obj_filtered)
seurat_obj_filtered <- ScaleData(seurat_obj_filtered)
seurat_obj_filtered <- RunPCA(seurat_obj_filtered)

# 参数优化
sweep_res <- paramSweep(seurat_obj_filtered, PCs = 1:10)
sweep_stats <- summarizeSweep(sweep_res, GT = FALSE)
bcmvn <- find.pK(sweep_stats)

# 预测双细胞
optimal_pk <- as.numeric(as.character(bcmvn$pK[which.max(bcmvn$BCmetric)]))
annotations <- seurat_obj_filtered$seurat_clusters
homotypic_prop <- modelHomotypic(annotations)
nExp_poi <- round(0.05 * nrow(seurat_obj_filtered@meta.data))
nExp_poi_adj <- round(nExp_poi * (1 - homotypic_prop))

seurat_obj_filtered <- doubletFinder(seurat_obj_filtered,
                                     PCs = 1:10,
                                     pN = 0.25,
                                     pK = optimal_pk,
                                     nExp = nExp_poi_adj,
                                     reuse.pANN = FALSE)

# 查看双细胞标签
table(seurat_obj_filtered@meta.data$DF.classifications_0.25_0.004_49)

⚠️ 注意： 双细胞检测依赖于聚类结果，所以它是个迭代过程。我建议第一次先用宽松的QC做初步聚类，检测完双细胞去除后，再用更严格的QC跑第二轮。这样能最大程度减少假阳性。

五、完整的QC流程图

下面这张图展示了从原始数据到高质量数据的完整QC流程。重点在于每一步的递进关系：先读数据、再看指标、再过滤、再检测双细胞、最后检查。

六、常见问题与避坑

做QC这么多年，我整理了一些常见问题，你也可以对照着看看自己是不是也遇到过：

过滤后细胞太少？ 大概率是阈值设得太严格了。我建议先拿小提琴图看看分布，尤其注意"自然断点"的位置——阈值就设在断点处。
线粒体比例阈值怎么定？ 没有统一的答案。脑组织通常温和，5-10%算高；肝脏组织代谢活跃，15-20%可能都正常。一个笨办法：看分布图，找到大部分细胞的集中区域，阈值设在尾部开始的地方。
要不要过滤UMI总数？ 我个人习惯是过滤掉nFeature_RNA极低的细胞后，nCount_RNA一般也会跟着正常。除非特别异常（比如突然有个细胞的UMI数比中位数高出10倍），否则我很少单独设nCount的阈值。

⚠️ 必须记住： QC不是越严格越好。过分严格的QC会把稀有细胞类型过滤掉，导致后面的分析丢失生物学信息。我通常宁可"漏杀"也不要"错杀"——尤其是那些基因数低但生物学意义明确的细胞类型，比如红细胞、血小板、粒细胞等。

七、QC后的数据检查

过滤完别急着跑聚类。我习惯用下面这几步做个最终确认：

看过滤前后的细胞数变化——记录下过滤掉了多少细胞，心里有数。
重新跑一次QC小提琴图——确保剩下的细胞都在合理范围内。
检查基因-UMI关系——正常情况下，基因数和UMI数应该正相关，而且没有明显的离群点。
看看双细胞比例——一般控制在5%以内比较合理，太高说明数据质量本身就有问题。

好了，这就是我处理单细胞数据预处理的全部流程。你想想看，如果花半小时把这些步骤做扎实了，后面几天的分析都会顺畅得多。相反，要是图省事跳过了QC，等你聚类出来发现一堆奇怪的"新细胞类型"时，多半就是双细胞或者死细胞在捣乱了。