第3章：质量控制——QC指标计算、小提琴图与散点图过滤、细胞/基因筛选

3.1 为什么要做质量控制？

单细胞数据刚下机时，里边混着不少“次品”——死细胞、双胞、甚至空液滴。如果你直接拿这些去聚类，出来的结果可能是一团浆糊，后续差异分析也全被噪声淹没。所以，QC这一步不是可选项，是必须要过的坎。

我记得刚入行时，接过一个肺癌样本，跑完聚类发现有一大坨“肝细胞”标记基因高表达。排查半天才意识到是死细胞释放的线粒体RNA在捣乱。从那以后，我给自己定了条规矩：数据再漂亮，也得先过QC这关。

3.2 QC指标计算——三把标尺

Seurat里常用的QC指标就三个：

• nCount_RNA：每个细胞的总UMI数。太低可能是空液滴或细胞破损，太高可能是双胞。
• nFeature_RNA：每个细胞检测到的基因数。原理跟UMI数类似，但更能反映基因捕获效率。
• percent.mt：线粒体基因占比。坏死或凋亡细胞会大量释放线粒体RNA，这比例通常很高。

我习惯用PercentageFeatureSet这个函数来计算线粒体比例。代码很简洁：

# 计算线粒体基因比例（标准前缀为MT-）
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

如果你用的人源数据，记得大小写——基因名里MT-是大写。鼠源的话是mt-，别搞混。

指标	含义	常用过滤范围	异常时会发生什么
nCount_RNA	总UMI数	500 – 8000	过低 → 空液滴；过高 → 双胞
nFeature_RNA	检测到的基因数	200 – 5000	过低 → 细胞碎片；过高 → 双胞
percent.mt	线粒体基因占比	0 – 20%	过高 → 坏死/凋亡细胞

3.3 小提琴图与散点图——用眼睛做第一步过滤

数值是一回事，但分布长什么样更重要。我通常先画小提琴图：

# 小提琴图展示三个指标的分布
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

小提琴图的“肚子”越鼓越好，说明大部分细胞集中在正常范围。如果尾巴出奇地长，或者出现多个峰，那里面有鬼。比如percent.mt的小提琴图如果右侧拖出老长一根须，那些细胞大概率是快死的。

接着画散点图，把两个指标关联起来看：

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

这张图我特别喜欢。正常情况下，nCount_RNA和percent.mt呈负相关——UMI数越少的细胞线粒体比例越高。如果你看到一大坨点集中在左上角（低UMI高线粒体），那准是空液滴或者碎片的信号。

还有一张经典组合图：

FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

正常细胞会沿着一条斜线分布。如果出现大量离群点——比如nCount_RNA很高但nFeature_RNA很低——那可能是双胞（两个细胞被当成一个了）。

3.4 过滤——动刀要稳准狠

看完分布，下一步就是过滤。Seurat里用的是subset函数：

# 根据常用阈值过滤细胞
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 &
               nCount_RNA > 500 & nCount_RNA < 8000 &
               percent.mt < 20)

过滤细胞的逻辑写在一个条件里，一目了然。我习惯把阈值设得比理论值宽松一点，因为过严的过滤会丢掉稀有细胞类型。比如某种细胞本身就很小，nFeature_RNA可能就只有300，你一棍子打死就再也找不回来了。

除了过滤细胞，基因也要过滤：

# 在尽可能多的细胞中表达的基因才保留
# Seurat默认在创建Seurat对象时已做过初步过滤（min.cells = 3）
# 但你可以手动再收紧
counts <- GetAssayData(pbmc, assay = "RNA", slot = "counts")
nonzero <- rowSums(counts > 0)
genes_to_keep <- names(which(nonzero > 10))  # 至少在10个细胞中表达
pbmc <- pbmc[genes_to_keep, ]

为什么过滤基因？那些只在1-2个细胞中表达的基因，几乎全是技术噪声。留着它们只会稀释真正的信号，增加多重假设检验的负担。

3.5 实战：一套完整的QC流程

下面是我个人习惯用的QC脚本，基本每个项目都会拿它来开刀：

# 加载数据后立即执行
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 可视化检查分布
p1 <- VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
p2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
p3 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
print(p1 + p2 + p3 + plot_layout(ncol = 3))

# 观察分布后手动调整阈值
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 &
               nCount_RNA > 500 & nCount_RNA < 8000 &
               percent.mt < 20)

# 基因过滤：保留至少在10个细胞中表达的基因
counts <- GetAssayData(pbmc, assay = "RNA", slot = "counts")
keep_genes <- rowSums(counts > 0) >= 10
pbmc <- pbmc[keep_genes, ]

# 再次可视化确认过滤效果
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

做完这一步，你的数据应该变得清爽很多。小提琴图的右尾被切掉了，散点图上的离群点也不见了。接下来就可以放心地去做标准化和降维了。