1. 我为什么要把质控这件事单独写一篇如果你跟我一样接手过几套单细胞转录组数据那你大概率体会过这种感觉数据下载完暴力跑完Seurat全流程满心欢喜打开UMAP图发现细胞团全部糊成一团marker基因的热图也画不出层次感。这个时候你第一反应是这个样本生物学就是这么复杂其实大概率是刚开始的质控步骤没做对。三年前我第一次处理10X平台的单细胞数据时就是这副惨样。当时项目催得急我拿着默认参数就直接开跑过滤后的对象里有12000多个细胞看起来数量挺漂亮但整个数据怎么调都出不了清晰的分群。后来静下心来重新看数据发现光质控这一关就滤掉了接近三分之一——细胞数从12000掉到8000出头但重新跑出来的UMAP分群干净利落marker基因的表达模式终于跟文献对上了。那次返工让我彻底明白了一个道理单细胞数据分析的成败一半决定在质控。不是说后面的聚类、差异分析不重要而是如果进来的细胞本身就有大量的破损细胞、双重体、空滴后面所有结果都会受到污染。这篇文章我不打算讲整套单细胞分析流程而是把用Seurat过滤低质量细胞这一件事拆开揉碎从指标含义、阈值设定到完整代码、常见翻车场景一次性讲透。这篇文章适合刚开始接触单细胞转录组数据分析的研究生、想优化自己质控流程的生信从业者也适合那些跑完流程但结果总是不对劲、怀疑自己哪里做错了的同学。2. 数据到手先别急先理解Seurat对象里的三个质控指标很多新手拿到10X的filtered_feature_bc_matrix文件夹第一件事就是赶紧读进Seurat开始跑这其实是在跳过最关键的一步你得先知道你手里的数据是什么样子的、哪些细胞其实根本不能用。在Seurat里做质控翻来覆去就三个核心指标nFeature_RNA、nCount_RNA、percent.mt。下面一个一个说清楚。2.1 nFeature_RNA每个细胞检测到的基因数nFeature_RNA直译过来就是检测到的基因数量代表这个细胞里到底有多少种基因被捕捉到了。它是判断细胞状态最敏感的指标之一。一个健康完整的细胞在10X平台上通常能检测到几百到几千个基因。如果检测到的基因数特别低比如只有二三十个那这个所谓的细胞大概率是一个被抽破了、mRNA几乎漏光了的空泡或者是一个即将凋亡、转录本已经降解得差不多的垂死细胞。反过来如果nFeature_RNA高得离谱比如超过一万那它很可能是两个甚至多个细胞被同一个GEM微滴包裹了也就是所谓双重体doublet或多细胞聚合体。这里有一个很容易被忽略的细节nFeature_RNA的判断标准跟细胞类型是强相关的。中性粒细胞、某些静止的T细胞亚群它们本身转录本就很低nFeature只有几百是正常的而肿瘤细胞或高度活化的B细胞nFeature高也说得通。所以你不能一概而论地认为nFeature小于500就是垃圾细胞也不能认为nFeature大于6000就一定是双重体。质控的精髓恰恰在于结合你的样本类型和数据分布去判断。2.2 nCount_RNA每个细胞的UMI总数nCount_RNA代表这个细胞里所有基因的UMI计数总和可以粗略理解成这个细胞被捕获到的转录本总量。它和nFeature_RNA是高度相关的——转录本总量越高自然检出的基因种类也会更多。你在做质控的时候通常会看到nCount_RNA低得离谱的那一批往往跟nFeature低的细胞高度重合也就是同一条低质量线。但两边的信息不完全冗余有些细胞尽管nFeature看着还行但nCount极低说明这个细胞整体的转录活跃度已经很低了很可能是处于应激或者濒死状态。Seurat官方的质控教程里也特别强调nFeature和nCount要组合起来看这就是为什么要画FeatureScatter的原因——两个指标做散点图能很直观地看出哪些点偏离了正常的相关性区间。2.3 线粒体基因比例判断细胞状态的重要线索percent.mt指的是线粒体基因的UMI计数占总UMI的比例。这个指标算是质控里的黄金指标因为它在生物学机制上特别说得通一个健康细胞拥有完整的细胞膜和细胞质细胞质里的mRNA整体是保留完好的而一个细胞膜破损、胞质丢失的垂死细胞细胞质里的mRNA大量流失但线粒体因为存在于相对稳定的细胞器结构中往往留在细胞碎片里于是线粒体转录本占总转录本的比例就会异常升高。用人话翻译一下线粒体比例高说明这个细胞的浆漏了大概率处于凋亡或机械破损状态。这个指标在过滤时极其有效因为它是独立于nFeature和nCount的另一个维度能把前面两个指标看不出来的漏浆细胞揪出来。这也是为什么Seurat官方教程里会算出一个percent.mt来几乎没有一个单细胞项目会跳过它。Seurat里有一行固定操作pbmc[[percent.mt]] - PercentageFeatureSet(pbmc, pattern ^MT-)。注意这行的人源数据是^MT-小鼠则是^mt-大小写不要搞错否则算出来的线粒体比例永远是0。这个低级错误我在多个同学的脚本里见过血的教训。3. 完整代码从Cell Ranger输出到一个干净的单细胞对象理解了指标含义之后就该上手操作了。下面这套代码是我在实际项目中打磨过的版本可以直接复制修改使用但我不建议无脑跑因为每个步骤背后的参数含义你需要看懂。3.1 准备环境与加载数据质控这一步对R和Seurat的版本要求并不苛刻R 4.2以上、Seurat 4.x或5.x都行。我目前主要用Seurat 5但下面这套代码在4.x上同样能跑。需要的基础依赖包括dplyr和Matrix通常装Seurat的时候会自动带上来。library(Seurat) library(dplyr) library(Matrix) # 读取Cell Ranger的filtered_feature_bc_matrix目录 # 目录下应该有barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz三个文件 counts - Read10X(data.dir filtered_feature_bc_matrix/) # 如果拿到的是h5文件这样读 # counts - Read10X_h5(filtered_feature_bc_matrix.h5)这里有个小知识点Read10X()读进来的是一个稀疏矩阵行是基因列是细胞barcode。Cell Ranger在得到这个矩阵之前其实做了一步空滴过滤EmptyDrops把很大一部分没有细胞但背景RNA污染产生的barcode去掉了。但即便经过了这一步矩阵里剩下的barcode也仍然混杂着大量低质量细胞——Cell Ranger的空滴过滤只能筛掉绝对空筛不掉破了但还有残留转录本的细胞。这就是为什么我们还需要Seurat来再做一轮质控。3.2 创建Seurat对象并计算质控指标# 创建Seurat对象 sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)min.cells 3的意思是一个基因至少在3个细胞里检测到才被保留。这是为了过滤掉那些极低表达、基本属于背景噪声的基因一般设在3~10之间太小吃不干净噪声太大会误伤一些真实但低表达的基因。min.features 200的意思是一个细胞至少检测到200个基因才被保留这一步做的是一次非常粗的预过滤先把那些一看就是空泡的垃圾barcode剔除掉。注意这个200只是一个保底的初始值真正的精细过滤在下一步。然后是计算线粒体比例# 计算线粒体基因比例人源用 ^MT-小鼠用 ^mt- sce[[percent.mt]] - PercentageFeatureSet(sce, pattern ^MT-)PercentageFeatureSet()的本质就是把匹配了^MT-这些基因的总UMI计数求和再除以这个细胞所有基因的总UMI最后乘以100。所以这个指标是一个百分比数值范围从0到100。3.3 质控指标可视化拿到三个指标之后先别急着设阈值一定要可视化。不看图就设阈值约等于闭着眼开车。下面这几行代码就能帮你快速了解数据长什么样# 3个质控指标的分布小提琴图 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3, pt.size 0.01) # nCount与percent.mt的散点图 plot1 - FeatureScatter(sce, feature1 nCount_RNA, feature2 percent.mt) # nCount与nFeature的散点图 plot2 - FeatureScatter(sce, feature1 nCount_RNA, feature2 nFeature_RNA) plot1 plot2小提琴图能让你一眼看到三个指标的分布形态和异常尾巴。散点图则能帮你检查这些指标的联动关系是否合理。正常情况下nCount和nFeature应该呈明显的正相关如果你看到一大团细胞nCount很低但nFeature却很高或者反过来那都是信号异常需要格外注意。3.4 核心过滤操作看完了分布图才轮到核心的过滤代码。以下是我在多数外周血和肿瘤组织样本上验证过的一套初始参数# 过滤 # 阈值需要结合你的样本类型和数据分布灵活调整下面只是常见起始值 sce_filtered - subset(sce, subset nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15)注意subset()函数会把不符合条件的细胞直接从对象中移除。过滤完之后建议再看一次小提琴图和散点图确认过滤后的数据分布是否合理细胞数变化多少。这一步是所有质控中最核心的一环也是最容易被忽视的一环——很多人过滤完就直接进入NormalizeData了根本不做验证。3.5 过滤后的输出与常用检查# 查看过滤前后细胞数变化 ncol(sce) ncol(sce_filtered) # 查看过滤后对象的详细信息 sce_filtered # 查看每个细胞的质控指标统计 summary(sce_filtered$nFeature_RNA) summary(sce_filtered$nCount_RNA) summary(sce_filtered$percent.mt)如果过滤比例适中一般滤掉10%~30%说明数据质量还可以。如果一次性滤掉了一半以上就要回头想想是组织解离过程出了问题还是阈值设得太狠还是这个样本本身就比较脏这些都需要结合实际情况判断而不是简单调低阈值把细胞硬留下来。4. 阈值不是拍脑袋基于分布特征确定质控参数上面的代码里给的nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15是一个起步参考值离精准还有距离。为什么这么说因为这个参数组合是我在特定样本类型上调出来的换一个样本类型、换一个测序平台它就可能完全不适用。Seurat官方教程里的经典参数是线粒体比例小于5%那是基于PBMC数据集给出的示例很多初学者直接把它套到自己的肿瘤或组织样本上结果滤掉了大量真实细胞。4.1 默认参数为什么经常翻车单细胞转录组的本质是抽样测序抽样的效率跟建库批次、10X试剂盒版本v2、v3、Chromium X、测序深度都有关系。同样的细胞类型在不同批次里测出来的nFeature中位数可能相差上千。如果你拿一个中位nFeature只有800的冷冻肿瘤样本去套用PBMC的200~2500区间那你的前半段正常细胞就会被无情误杀。更让人头疼的是细胞类型差异。我处理过一个中性粒细胞比例很高的样本中性粒细胞nFeature中位数不到500按nFeature必须大于600的常规思路过滤整个中性粒细胞群几乎被团灭。可是人家中性粒细胞本来就是短命、转录本少的细胞你把它们全删了后续分析还怎么做这就是精准两个字的分量过滤不是要套公式而是要理解你的数据里有哪些细胞、它们的指标形态是怎样的。4.2 用MAD法则和分布形态来定阈值那怎么定阈值才算靠谱我自己的习惯是两步走先看分布再用分位数辅助判断。第一步用核密度图看nFeature和nCount的整体形态特别注意有没有一个矮子峰# 画核密度图观察nFeature_RNA的分布 plot(density(sce$nFeature_RNA), main nFeature density, col steelblue) # 观察percent.mt的分布 plot(density(sce$percent.mt), main percent.mt density, col firebrick)如果nFeature_RNA的密度图在左边低基因数区域出现一个独立的小峰或者峰型很扁那说明数据里确实混了一批低质量的细胞。而percent.mt的密度图如果出现一个拖得很长的尾巴或者主峰明显右偏那线粒体高比例的细胞一定不在少数。第二步就是更定量的做法。很多人会采用中位数±3倍绝对中位差MAD来帮助确定阈值。它的逻辑跟平均值±3倍标准差类似但MAD对异常值更不敏感更适合单细胞这种厚尾分布的数据# 用中位数±3倍MAD估计合理范围仅供初步参考 lower_bound - median(sce$nFeature_RNA) - 3 * mad(sce$nFeature_RNA) upper_bound - median(sce$nFeature_RNA) 3 * mad(sce$nFeature_RNA) cat(nFeature_RNA 参考范围:, lower_bound, -, upper_bound, \n)这套方法的好处是阈值完全由你自己的数据驱动而不是生搬硬套别人的经验。但它同样有局限性——如果你的数据里低质量细胞占比太高MAD本身就会被它们拉着走导致算出来的范围过宽。所以我的建议是MAD算出来的范围用来做参考最终阈值还是回到密度图上人工确认一下看看切的位置是否正好落在主峰和低质量小峰之间的洼地上。4.3 不同样本类型的初始阈值参考为了让你有个大概的手感我把自己处理过的几类样本的初始阈值整理成了表格。注意这只是第一轮粗过滤的参考值不代表所有项目通用样本类型nFeature_RNA 参考下限nFeature_RNA 参考上限percent.mt 参考上限PBMC 外周血200~4004000~500010~15肿瘤组织新鲜解离300~5006000~800015~20冷冻组织200~3004000~600020~25中性粒细胞富集样本100~2003000~400010~15snRNA-seq核测序500~8008000~100001~5表格里的数值跨度比较大因为每批数据的质量波动实在太大。我一直强调一个观点表格里的数字只是用来入门的不是你交论文时的最终答案。最终答案永远要靠你对自己样本的理解和数据分布的实际观察来回答。这也是精准质控和应付差事式质控的核心差别。5. 过滤完别急着开心怎么验证质控效果过滤代码跑通了细胞数也减少了并不代表质控做完了。很多人在这时候马上进入NormalizeData、PCA、UMAP的下一个阶段结果最后发现分群有问题才想起回头找质控的问题。我更建议你在过滤之后用几分钟做几项验证把异常在早期就暴露出来。5.1 过滤前后指标的对比最直接的验证就是对比过滤前后的指标分布# 过滤前后指标分布对比 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3) VlnPlot(sce_filtered, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)过滤后你应该看到nFeature_RNA和nCount_RNA的低值尾被切掉了整体分布更集中percent.mt的高值尾巴被削平大部分细胞集中在一个相对正常的区间小提琴图整体变得干净不再有明显的离群点。如果过滤后你的nFeature主峰反而被砍掉一大半或者percent.mt的分布几乎没有变化那说明阈值设得有问题要么误伤了大群正常细胞要么压根没把低质量细胞滤除干净。5.2 检查高表达基因列表是否合理这是一个非常容易被忽略但很聪明的检查方法。正常情况下一个健康单细胞数据集中高表达基因应当是各类细胞类型的功能基因比如PBMC里应该是MALAT1、B2M、各类核糖体蛋白基因等。但如果你的数据里混入了大量垂死细胞或破裂细胞高表达基因排行榜会被线粒体基因霸榜——因为线粒体基因在破损细胞里占比太高了。我有个习惯性操作质控后跑一下每个细胞的top基因看看排名靠前的是不是有大量MT-基因。具体可以这样# 检查每个细胞中最显著高表达的基因 # 取前500个高表达基因 top_genes - rownames(sce_filtered)[ order(Matrix::rowSums(sce_filteredassays$RNA$counts), decreasing TRUE) ] head(top_genes, 30)如果你看到head出来的30个基因里有大量MT-开头的基因那就要警惕线粒体基因过滤可能还不够干净尤其在有核红细胞、心肌组织等线粒体天然丰富的样本里。至于要不要进一步往上提高percent.mt阈值同样还是回到分布图和生物学实际情况来决定。5.3 双重体到底要不要处理过滤验证阶段最容易被忽略的其实是双重体doublets。Seurat的质控指标可以筛掉低质量细胞但对双重体的效果非常有限——因为一个双重体的nFeature通常并不是特别离谱它可能是两个中等基因数的细胞粘在一起生成的nFeature值恰好落在正常范围的高值区蒙混过关。要严格去双就得用专门的工具。目前最常用的是DoubletFinder# 这个步骤建议在质控完成、数据标准化之后运行 # 简单用法示意 # library(DoubletFinder) # sce_filtered - NormalizeData(sce_filtered) # sce_filtered - FindVariableFeatures(sce_filtered) # sce_filtered - ScaleData(sce_filtered) # sce_filtered - RunPCA(sce_filtered) # sce_filtered - FindNeighbors(sce_filtered, dims 1:20) # sce_filtered - FindClusters(sce_filtered, resolution 0.5) # sce_filtered - RunUMAP(sce_filtered, dims 1:20) # # 然后基于聚类结果预测doublet # sweep.res - paramSweep(sce_filtered, PCs 1:20, sct FALSE) # sweep.stats - summarizeSweep(sweep.res, GT FALSE) # bcmvn - find.pK(sweep.stats) # ...DoubletFinder的原理是通过人工模拟doublet看哪些细胞在特征空间上与模拟的doublet最相似从而打分预测。这个工具能筛掉相当一部分漏网的doublet但它也有误杀风险而且是基于聚类结果预测的如果前面聚类本身就不好预测的可靠性也会受影响。所以我通常建议在分群结果出来后、确定marker基因之前结合doublet注释结果做一次排查尽量不要在只有几十个细胞的支持下强行去双。5.4 与后续分析流程的正确衔接质控是起点不是终点。过滤干净之后后续步骤必须衔接好# 完整的下游分析衔接 sce_filtered - NormalizeData(sce_filtered, normalization.method LogNormalize, scale.factor 10000) sce_filtered - FindVariableFeatures(sce_filtered, selection.method vst, nfeatures 2000) sce_filtered - ScaleData(sce_filtered) sce_filtered - RunPCA(sce_filtered, npcs 30) sce_filtered - RunUMAP(sce_filtered, dims 1:20) sce_filtered - FindNeighbors(sce_filtered, dims 1:20) sce_filtered - FindClusters(sce_filtered, resolution 0.5) # 后续差异分析 markers - FindAllMarkers(sce_filtered, only.pos TRUE, min.pct 0.25, logfc.threshold 0.25)这里我不展开讲每个步骤的细节但要注意的是质控越靠前后面全部步骤都能受益如果你在整合多个样本最好先分别对每个样本做各自的质控再进入整合流程。把不同质量的样本混在一起再统一过滤往往会出现好的样本被拉到差的样本水平这种尴尬处境。6. 实战中的特殊样本那些容易翻车的质控场景写代码部分的质控大家都能跑通真正拉开差距的是处理那些非典型样本。我在多个项目里踩过坑下面几个场景几乎每个做单细胞的人都可能遇到专门拿出来说。6.1 免疫细胞样本的线粒体比例天然偏高如果你是做免疫细胞、尤其是组织浸润免疫细胞的一定不要迷信Seurat官方教程的percent.mt 5。5%这个阈值来自PBMC外周血的样本背景干净、细胞状态比较均一。但组织里的T细胞、巨噬细胞在解离过程中经历了机械剪切、酶消化应激状态下线粒体活性会升高percent.mt就很容易冲到10%~20%甚至更高。我处理过一批肿瘤浸润T细胞样本percent.mt的中位数就有14%如果按5%切直接砍掉一大半T细胞后面根本没法分析。这种情况下我的建议是把阈值放到20%甚至25%同时结合nFeature和nCount来判断——一个T细胞如果percent.mt是20%但nFeature有1200它大概率还是一条好汉但如果percent.mt 20%同时nFeature只有150那它大概率确实快死了。6.2 肿瘤样本低质量与恶性信号混杂的难题肿瘤组织是单细胞项目里最复杂的一类样本。肿瘤细胞本身存在拷贝数变异、倍性异常、增殖活跃它们的nFeature和nCount往往很高高到如果用nFeature_RNA 5000一刀切会把真正的肿瘤细胞当成doublet去掉。另一方面肿瘤微环境里又有大量浸润的免疫细胞、基质细胞各类细胞之间的nFeature跨度极大。在这种样本里做质控我的策略是在解离后先做一个偏保守的粗过滤确保把明显的垃圾先清掉然后在聚类之后用inferCNV等工具确认肿瘤细胞群单独看它们的质控指标分布再做一轮针对性的二次检视。这一步很关键——如果你在前期就下手太狠可能把一个重要的肿瘤亚群给做没了后面再怎么注释都找不回来。6.3 snRNA-seq核测序千万不要当成普通单细胞处理单细胞核测序snRNA-seq这几年特别火尤其在冷冻组织中优势明显因为核不容易受到解离损伤的影响。但它的质控指标跟普通单细胞完全不一样核里面不含完整的细胞质线粒体基因本来就极少percent.mt通常只有1%以下同时核转录本的捕获效率更高nFeature中位数往往比完整细胞更高。如果有人拿着普通单细胞的标准去过滤snRNA-seq数据大概率会把数据搞坏。snRNA-seq的质控重点是排除带大量细胞质残留的完整细胞碎片和空的细胞核而不是纠结线粒体比例。另外snRNA-seq中低nFeature并不一定代表低质量——有些静止细胞的核转录本本来就少。如果后面还要结合ATAC数据做多组学分析质控更要谨慎。6.4 污染与背景RNA的影响做过组织样本的朋友应该都有体会组织解离之后总有一些细胞碎片和游离RNA它们会被随机包裹进GEM微滴。这种背景RNA污染会造成一个现象所有细胞的nFeature和nCount都被稀释了——具体来说垃圾细胞带着高背景RNA会让整体数据出现一个低质量团块。Seurat的普通质控对这种污染有一定的清理效果但未必能完全解决。如果背景污染比较严重我更推荐用CellBender或者SoupX这类专门的去污染工具在质控之前先把背景RNA信号扣除掉。这个操作能让后续阈值判断更精准但要注意这类工具对计算资源和参数设置有一定要求不是所有项目都必需。一般建议是如果聚类后所有细胞类型都带有同样一组高表达基因、UMAP图分群模糊并且高表达基因列表里有明显的环境RNA特征比如各种血红蛋白基因才需要考虑去污染。7. 写在最后的实操心得文章写到这儿代码和理论都讲得差不多了。最后分享几点只有自己上手跑过几批数据才会总结出来的体会。第一质控阈值不要追求一次到位。我自己现在跑一个新数据集很少会第一次就定出完美的阈值。通常是先按照分布图设一个大概范围过滤后看看细胞数和UMAP效果再回来微调。这个跑一遍—看结果—回来调的循环在单细胞项目里极其正常不要觉得麻烦。第二质控日志一定要记。具体包括过滤前多少细胞、过滤后多少细胞、三个指标的阈值各是多少、当时看了哪些图、为什么这么定。不要高估自己的记忆力。一篇论文从投稿到返修可能隔了大半年如果到时候审稿人问你你的过滤标准是怎么确定的而你连跑时的参数都忘了那场面会很尴尬。我现在每个项目都会留一个quality_control_notes.md文件随手记录当时的判断逻辑。第三不同样本之间不要套用同一个阈值。哪怕两个样本来自同一个病人的同一类组织只要批次不同、建库时间不同就应该各看各的分布再分别定阈值。把两个样本合并成一个对象后统一过滤是新手最容易犯的错误之一它会直接导致高质量的样本被你误伤。第四质控的上限不只是不滤错还要考虑下一步需要什么。比如你后面想做拟时序分析那G0/G1期的细胞就要尽量保留如果后面想做转录因子活性分析转录本种类不够丰富的细胞可能支撑不起结论。每一种下游分析对数据质量的要求都不同。质控参数不是越严格越好而是越契合你的分析目标越好。单细胞转录组质控这条路说简单也简单就是看三个指标、设三个阈值说复杂也复杂因为每个数据集都有它自己的脾气。如果你把这篇内容看懂了再回到自己的数据上去实践几轮应该能明显感受到精准质控和套参数跑流程之间的差距。希望这些经验能帮你少走一点弯路。 单细胞转录组质控全解析：用Seurat精准过滤低质量细胞
1. 我为什么要把质控这件事单独写一篇如果你跟我一样接手过几套单细胞转录组数据那你大概率体会过这种感觉数据下载完暴力跑完Seurat全流程满心欢喜打开UMAP图发现细胞团全部糊成一团marker基因的热图也画不出层次感。这个时候你第一反应是这个样本生物学就是这么复杂其实大概率是刚开始的质控步骤没做对。三年前我第一次处理10X平台的单细胞数据时就是这副惨样。当时项目催得急我拿着默认参数就直接开跑过滤后的对象里有12000多个细胞看起来数量挺漂亮但整个数据怎么调都出不了清晰的分群。后来静下心来重新看数据发现光质控这一关就滤掉了接近三分之一——细胞数从12000掉到8000出头但重新跑出来的UMAP分群干净利落marker基因的表达模式终于跟文献对上了。那次返工让我彻底明白了一个道理单细胞数据分析的成败一半决定在质控。不是说后面的聚类、差异分析不重要而是如果进来的细胞本身就有大量的破损细胞、双重体、空滴后面所有结果都会受到污染。这篇文章我不打算讲整套单细胞分析流程而是把用Seurat过滤低质量细胞这一件事拆开揉碎从指标含义、阈值设定到完整代码、常见翻车场景一次性讲透。这篇文章适合刚开始接触单细胞转录组数据分析的研究生、想优化自己质控流程的生信从业者也适合那些跑完流程但结果总是不对劲、怀疑自己哪里做错了的同学。2. 数据到手先别急先理解Seurat对象里的三个质控指标很多新手拿到10X的filtered_feature_bc_matrix文件夹第一件事就是赶紧读进Seurat开始跑这其实是在跳过最关键的一步你得先知道你手里的数据是什么样子的、哪些细胞其实根本不能用。在Seurat里做质控翻来覆去就三个核心指标nFeature_RNA、nCount_RNA、percent.mt。下面一个一个说清楚。2.1 nFeature_RNA每个细胞检测到的基因数nFeature_RNA直译过来就是检测到的基因数量代表这个细胞里到底有多少种基因被捕捉到了。它是判断细胞状态最敏感的指标之一。一个健康完整的细胞在10X平台上通常能检测到几百到几千个基因。如果检测到的基因数特别低比如只有二三十个那这个所谓的细胞大概率是一个被抽破了、mRNA几乎漏光了的空泡或者是一个即将凋亡、转录本已经降解得差不多的垂死细胞。反过来如果nFeature_RNA高得离谱比如超过一万那它很可能是两个甚至多个细胞被同一个GEM微滴包裹了也就是所谓双重体doublet或多细胞聚合体。这里有一个很容易被忽略的细节nFeature_RNA的判断标准跟细胞类型是强相关的。中性粒细胞、某些静止的T细胞亚群它们本身转录本就很低nFeature只有几百是正常的而肿瘤细胞或高度活化的B细胞nFeature高也说得通。所以你不能一概而论地认为nFeature小于500就是垃圾细胞也不能认为nFeature大于6000就一定是双重体。质控的精髓恰恰在于结合你的样本类型和数据分布去判断。2.2 nCount_RNA每个细胞的UMI总数nCount_RNA代表这个细胞里所有基因的UMI计数总和可以粗略理解成这个细胞被捕获到的转录本总量。它和nFeature_RNA是高度相关的——转录本总量越高自然检出的基因种类也会更多。你在做质控的时候通常会看到nCount_RNA低得离谱的那一批往往跟nFeature低的细胞高度重合也就是同一条低质量线。但两边的信息不完全冗余有些细胞尽管nFeature看着还行但nCount极低说明这个细胞整体的转录活跃度已经很低了很可能是处于应激或者濒死状态。Seurat官方的质控教程里也特别强调nFeature和nCount要组合起来看这就是为什么要画FeatureScatter的原因——两个指标做散点图能很直观地看出哪些点偏离了正常的相关性区间。2.3 线粒体基因比例判断细胞状态的重要线索percent.mt指的是线粒体基因的UMI计数占总UMI的比例。这个指标算是质控里的黄金指标因为它在生物学机制上特别说得通一个健康细胞拥有完整的细胞膜和细胞质细胞质里的mRNA整体是保留完好的而一个细胞膜破损、胞质丢失的垂死细胞细胞质里的mRNA大量流失但线粒体因为存在于相对稳定的细胞器结构中往往留在细胞碎片里于是线粒体转录本占总转录本的比例就会异常升高。用人话翻译一下线粒体比例高说明这个细胞的浆漏了大概率处于凋亡或机械破损状态。这个指标在过滤时极其有效因为它是独立于nFeature和nCount的另一个维度能把前面两个指标看不出来的漏浆细胞揪出来。这也是为什么Seurat官方教程里会算出一个percent.mt来几乎没有一个单细胞项目会跳过它。Seurat里有一行固定操作pbmc[[percent.mt]] - PercentageFeatureSet(pbmc, pattern ^MT-)。注意这行的人源数据是^MT-小鼠则是^mt-大小写不要搞错否则算出来的线粒体比例永远是0。这个低级错误我在多个同学的脚本里见过血的教训。3. 完整代码从Cell Ranger输出到一个干净的单细胞对象理解了指标含义之后就该上手操作了。下面这套代码是我在实际项目中打磨过的版本可以直接复制修改使用但我不建议无脑跑因为每个步骤背后的参数含义你需要看懂。3.1 准备环境与加载数据质控这一步对R和Seurat的版本要求并不苛刻R 4.2以上、Seurat 4.x或5.x都行。我目前主要用Seurat 5但下面这套代码在4.x上同样能跑。需要的基础依赖包括dplyr和Matrix通常装Seurat的时候会自动带上来。library(Seurat) library(dplyr) library(Matrix) # 读取Cell Ranger的filtered_feature_bc_matrix目录 # 目录下应该有barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz三个文件 counts - Read10X(data.dir filtered_feature_bc_matrix/) # 如果拿到的是h5文件这样读 # counts - Read10X_h5(filtered_feature_bc_matrix.h5)这里有个小知识点Read10X()读进来的是一个稀疏矩阵行是基因列是细胞barcode。Cell Ranger在得到这个矩阵之前其实做了一步空滴过滤EmptyDrops把很大一部分没有细胞但背景RNA污染产生的barcode去掉了。但即便经过了这一步矩阵里剩下的barcode也仍然混杂着大量低质量细胞——Cell Ranger的空滴过滤只能筛掉绝对空筛不掉破了但还有残留转录本的细胞。这就是为什么我们还需要Seurat来再做一轮质控。3.2 创建Seurat对象并计算质控指标# 创建Seurat对象 sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)min.cells 3的意思是一个基因至少在3个细胞里检测到才被保留。这是为了过滤掉那些极低表达、基本属于背景噪声的基因一般设在3~10之间太小吃不干净噪声太大会误伤一些真实但低表达的基因。min.features 200的意思是一个细胞至少检测到200个基因才被保留这一步做的是一次非常粗的预过滤先把那些一看就是空泡的垃圾barcode剔除掉。注意这个200只是一个保底的初始值真正的精细过滤在下一步。然后是计算线粒体比例# 计算线粒体基因比例人源用 ^MT-小鼠用 ^mt- sce[[percent.mt]] - PercentageFeatureSet(sce, pattern ^MT-)PercentageFeatureSet()的本质就是把匹配了^MT-这些基因的总UMI计数求和再除以这个细胞所有基因的总UMI最后乘以100。所以这个指标是一个百分比数值范围从0到100。3.3 质控指标可视化拿到三个指标之后先别急着设阈值一定要可视化。不看图就设阈值约等于闭着眼开车。下面这几行代码就能帮你快速了解数据长什么样# 3个质控指标的分布小提琴图 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3, pt.size 0.01) # nCount与percent.mt的散点图 plot1 - FeatureScatter(sce, feature1 nCount_RNA, feature2 percent.mt) # nCount与nFeature的散点图 plot2 - FeatureScatter(sce, feature1 nCount_RNA, feature2 nFeature_RNA) plot1 plot2小提琴图能让你一眼看到三个指标的分布形态和异常尾巴。散点图则能帮你检查这些指标的联动关系是否合理。正常情况下nCount和nFeature应该呈明显的正相关如果你看到一大团细胞nCount很低但nFeature却很高或者反过来那都是信号异常需要格外注意。3.4 核心过滤操作看完了分布图才轮到核心的过滤代码。以下是我在多数外周血和肿瘤组织样本上验证过的一套初始参数# 过滤 # 阈值需要结合你的样本类型和数据分布灵活调整下面只是常见起始值 sce_filtered - subset(sce, subset nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15)注意subset()函数会把不符合条件的细胞直接从对象中移除。过滤完之后建议再看一次小提琴图和散点图确认过滤后的数据分布是否合理细胞数变化多少。这一步是所有质控中最核心的一环也是最容易被忽视的一环——很多人过滤完就直接进入NormalizeData了根本不做验证。3.5 过滤后的输出与常用检查# 查看过滤前后细胞数变化 ncol(sce) ncol(sce_filtered) # 查看过滤后对象的详细信息 sce_filtered # 查看每个细胞的质控指标统计 summary(sce_filtered$nFeature_RNA) summary(sce_filtered$nCount_RNA) summary(sce_filtered$percent.mt)如果过滤比例适中一般滤掉10%~30%说明数据质量还可以。如果一次性滤掉了一半以上就要回头想想是组织解离过程出了问题还是阈值设得太狠还是这个样本本身就比较脏这些都需要结合实际情况判断而不是简单调低阈值把细胞硬留下来。4. 阈值不是拍脑袋基于分布特征确定质控参数上面的代码里给的nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15是一个起步参考值离精准还有距离。为什么这么说因为这个参数组合是我在特定样本类型上调出来的换一个样本类型、换一个测序平台它就可能完全不适用。Seurat官方教程里的经典参数是线粒体比例小于5%那是基于PBMC数据集给出的示例很多初学者直接把它套到自己的肿瘤或组织样本上结果滤掉了大量真实细胞。4.1 默认参数为什么经常翻车单细胞转录组的本质是抽样测序抽样的效率跟建库批次、10X试剂盒版本v2、v3、Chromium X、测序深度都有关系。同样的细胞类型在不同批次里测出来的nFeature中位数可能相差上千。如果你拿一个中位nFeature只有800的冷冻肿瘤样本去套用PBMC的200~2500区间那你的前半段正常细胞就会被无情误杀。更让人头疼的是细胞类型差异。我处理过一个中性粒细胞比例很高的样本中性粒细胞nFeature中位数不到500按nFeature必须大于600的常规思路过滤整个中性粒细胞群几乎被团灭。可是人家中性粒细胞本来就是短命、转录本少的细胞你把它们全删了后续分析还怎么做这就是精准两个字的分量过滤不是要套公式而是要理解你的数据里有哪些细胞、它们的指标形态是怎样的。4.2 用MAD法则和分布形态来定阈值那怎么定阈值才算靠谱我自己的习惯是两步走先看分布再用分位数辅助判断。第一步用核密度图看nFeature和nCount的整体形态特别注意有没有一个矮子峰# 画核密度图观察nFeature_RNA的分布 plot(density(sce$nFeature_RNA), main nFeature density, col steelblue) # 观察percent.mt的分布 plot(density(sce$percent.mt), main percent.mt density, col firebrick)如果nFeature_RNA的密度图在左边低基因数区域出现一个独立的小峰或者峰型很扁那说明数据里确实混了一批低质量的细胞。而percent.mt的密度图如果出现一个拖得很长的尾巴或者主峰明显右偏那线粒体高比例的细胞一定不在少数。第二步就是更定量的做法。很多人会采用中位数±3倍绝对中位差MAD来帮助确定阈值。它的逻辑跟平均值±3倍标准差类似但MAD对异常值更不敏感更适合单细胞这种厚尾分布的数据# 用中位数±3倍MAD估计合理范围仅供初步参考 lower_bound - median(sce$nFeature_RNA) - 3 * mad(sce$nFeature_RNA) upper_bound - median(sce$nFeature_RNA) 3 * mad(sce$nFeature_RNA) cat(nFeature_RNA 参考范围:, lower_bound, -, upper_bound, \n)这套方法的好处是阈值完全由你自己的数据驱动而不是生搬硬套别人的经验。但它同样有局限性——如果你的数据里低质量细胞占比太高MAD本身就会被它们拉着走导致算出来的范围过宽。所以我的建议是MAD算出来的范围用来做参考最终阈值还是回到密度图上人工确认一下看看切的位置是否正好落在主峰和低质量小峰之间的洼地上。4.3 不同样本类型的初始阈值参考为了让你有个大概的手感我把自己处理过的几类样本的初始阈值整理成了表格。注意这只是第一轮粗过滤的参考值不代表所有项目通用样本类型nFeature_RNA 参考下限nFeature_RNA 参考上限percent.mt 参考上限PBMC 外周血200~4004000~500010~15肿瘤组织新鲜解离300~5006000~800015~20冷冻组织200~3004000~600020~25中性粒细胞富集样本100~2003000~400010~15snRNA-seq核测序500~8008000~100001~5表格里的数值跨度比较大因为每批数据的质量波动实在太大。我一直强调一个观点表格里的数字只是用来入门的不是你交论文时的最终答案。最终答案永远要靠你对自己样本的理解和数据分布的实际观察来回答。这也是精准质控和应付差事式质控的核心差别。5. 过滤完别急着开心怎么验证质控效果过滤代码跑通了细胞数也减少了并不代表质控做完了。很多人在这时候马上进入NormalizeData、PCA、UMAP的下一个阶段结果最后发现分群有问题才想起回头找质控的问题。我更建议你在过滤之后用几分钟做几项验证把异常在早期就暴露出来。5.1 过滤前后指标的对比最直接的验证就是对比过滤前后的指标分布# 过滤前后指标分布对比 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3) VlnPlot(sce_filtered, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)过滤后你应该看到nFeature_RNA和nCount_RNA的低值尾被切掉了整体分布更集中percent.mt的高值尾巴被削平大部分细胞集中在一个相对正常的区间小提琴图整体变得干净不再有明显的离群点。如果过滤后你的nFeature主峰反而被砍掉一大半或者percent.mt的分布几乎没有变化那说明阈值设得有问题要么误伤了大群正常细胞要么压根没把低质量细胞滤除干净。5.2 检查高表达基因列表是否合理这是一个非常容易被忽略但很聪明的检查方法。正常情况下一个健康单细胞数据集中高表达基因应当是各类细胞类型的功能基因比如PBMC里应该是MALAT1、B2M、各类核糖体蛋白基因等。但如果你的数据里混入了大量垂死细胞或破裂细胞高表达基因排行榜会被线粒体基因霸榜——因为线粒体基因在破损细胞里占比太高了。我有个习惯性操作质控后跑一下每个细胞的top基因看看排名靠前的是不是有大量MT-基因。具体可以这样# 检查每个细胞中最显著高表达的基因 # 取前500个高表达基因 top_genes - rownames(sce_filtered)[ order(Matrix::rowSums(sce_filteredassays$RNA$counts), decreasing TRUE) ] head(top_genes, 30)如果你看到head出来的30个基因里有大量MT-开头的基因那就要警惕线粒体基因过滤可能还不够干净尤其在有核红细胞、心肌组织等线粒体天然丰富的样本里。至于要不要进一步往上提高percent.mt阈值同样还是回到分布图和生物学实际情况来决定。5.3 双重体到底要不要处理过滤验证阶段最容易被忽略的其实是双重体doublets。Seurat的质控指标可以筛掉低质量细胞但对双重体的效果非常有限——因为一个双重体的nFeature通常并不是特别离谱它可能是两个中等基因数的细胞粘在一起生成的nFeature值恰好落在正常范围的高值区蒙混过关。要严格去双就得用专门的工具。目前最常用的是DoubletFinder# 这个步骤建议在质控完成、数据标准化之后运行 # 简单用法示意 # library(DoubletFinder) # sce_filtered - NormalizeData(sce_filtered) # sce_filtered - FindVariableFeatures(sce_filtered) # sce_filtered - ScaleData(sce_filtered) # sce_filtered - RunPCA(sce_filtered) # sce_filtered - FindNeighbors(sce_filtered, dims 1:20) # sce_filtered - FindClusters(sce_filtered, resolution 0.5) # sce_filtered - RunUMAP(sce_filtered, dims 1:20) # # 然后基于聚类结果预测doublet # sweep.res - paramSweep(sce_filtered, PCs 1:20, sct FALSE) # sweep.stats - summarizeSweep(sweep.res, GT FALSE) # bcmvn - find.pK(sweep.stats) # ...DoubletFinder的原理是通过人工模拟doublet看哪些细胞在特征空间上与模拟的doublet最相似从而打分预测。这个工具能筛掉相当一部分漏网的doublet但它也有误杀风险而且是基于聚类结果预测的如果前面聚类本身就不好预测的可靠性也会受影响。所以我通常建议在分群结果出来后、确定marker基因之前结合doublet注释结果做一次排查尽量不要在只有几十个细胞的支持下强行去双。5.4 与后续分析流程的正确衔接质控是起点不是终点。过滤干净之后后续步骤必须衔接好# 完整的下游分析衔接 sce_filtered - NormalizeData(sce_filtered, normalization.method LogNormalize, scale.factor 10000) sce_filtered - FindVariableFeatures(sce_filtered, selection.method vst, nfeatures 2000) sce_filtered - ScaleData(sce_filtered) sce_filtered - RunPCA(sce_filtered, npcs 30) sce_filtered - RunUMAP(sce_filtered, dims 1:20) sce_filtered - FindNeighbors(sce_filtered, dims 1:20) sce_filtered - FindClusters(sce_filtered, resolution 0.5) # 后续差异分析 markers - FindAllMarkers(sce_filtered, only.pos TRUE, min.pct 0.25, logfc.threshold 0.25)这里我不展开讲每个步骤的细节但要注意的是质控越靠前后面全部步骤都能受益如果你在整合多个样本最好先分别对每个样本做各自的质控再进入整合流程。把不同质量的样本混在一起再统一过滤往往会出现好的样本被拉到差的样本水平这种尴尬处境。6. 实战中的特殊样本那些容易翻车的质控场景写代码部分的质控大家都能跑通真正拉开差距的是处理那些非典型样本。我在多个项目里踩过坑下面几个场景几乎每个做单细胞的人都可能遇到专门拿出来说。6.1 免疫细胞样本的线粒体比例天然偏高如果你是做免疫细胞、尤其是组织浸润免疫细胞的一定不要迷信Seurat官方教程的percent.mt 5。5%这个阈值来自PBMC外周血的样本背景干净、细胞状态比较均一。但组织里的T细胞、巨噬细胞在解离过程中经历了机械剪切、酶消化应激状态下线粒体活性会升高percent.mt就很容易冲到10%~20%甚至更高。我处理过一批肿瘤浸润T细胞样本percent.mt的中位数就有14%如果按5%切直接砍掉一大半T细胞后面根本没法分析。这种情况下我的建议是把阈值放到20%甚至25%同时结合nFeature和nCount来判断——一个T细胞如果percent.mt是20%但nFeature有1200它大概率还是一条好汉但如果percent.mt 20%同时nFeature只有150那它大概率确实快死了。6.2 肿瘤样本低质量与恶性信号混杂的难题肿瘤组织是单细胞项目里最复杂的一类样本。肿瘤细胞本身存在拷贝数变异、倍性异常、增殖活跃它们的nFeature和nCount往往很高高到如果用nFeature_RNA 5000一刀切会把真正的肿瘤细胞当成doublet去掉。另一方面肿瘤微环境里又有大量浸润的免疫细胞、基质细胞各类细胞之间的nFeature跨度极大。在这种样本里做质控我的策略是在解离后先做一个偏保守的粗过滤确保把明显的垃圾先清掉然后在聚类之后用inferCNV等工具确认肿瘤细胞群单独看它们的质控指标分布再做一轮针对性的二次检视。这一步很关键——如果你在前期就下手太狠可能把一个重要的肿瘤亚群给做没了后面再怎么注释都找不回来。6.3 snRNA-seq核测序千万不要当成普通单细胞处理单细胞核测序snRNA-seq这几年特别火尤其在冷冻组织中优势明显因为核不容易受到解离损伤的影响。但它的质控指标跟普通单细胞完全不一样核里面不含完整的细胞质线粒体基因本来就极少percent.mt通常只有1%以下同时核转录本的捕获效率更高nFeature中位数往往比完整细胞更高。如果有人拿着普通单细胞的标准去过滤snRNA-seq数据大概率会把数据搞坏。snRNA-seq的质控重点是排除带大量细胞质残留的完整细胞碎片和空的细胞核而不是纠结线粒体比例。另外snRNA-seq中低nFeature并不一定代表低质量——有些静止细胞的核转录本本来就少。如果后面还要结合ATAC数据做多组学分析质控更要谨慎。6.4 污染与背景RNA的影响做过组织样本的朋友应该都有体会组织解离之后总有一些细胞碎片和游离RNA它们会被随机包裹进GEM微滴。这种背景RNA污染会造成一个现象所有细胞的nFeature和nCount都被稀释了——具体来说垃圾细胞带着高背景RNA会让整体数据出现一个低质量团块。Seurat的普通质控对这种污染有一定的清理效果但未必能完全解决。如果背景污染比较严重我更推荐用CellBender或者SoupX这类专门的去污染工具在质控之前先把背景RNA信号扣除掉。这个操作能让后续阈值判断更精准但要注意这类工具对计算资源和参数设置有一定要求不是所有项目都必需。一般建议是如果聚类后所有细胞类型都带有同样一组高表达基因、UMAP图分群模糊并且高表达基因列表里有明显的环境RNA特征比如各种血红蛋白基因才需要考虑去污染。7. 写在最后的实操心得文章写到这儿代码和理论都讲得差不多了。最后分享几点只有自己上手跑过几批数据才会总结出来的体会。第一质控阈值不要追求一次到位。我自己现在跑一个新数据集很少会第一次就定出完美的阈值。通常是先按照分布图设一个大概范围过滤后看看细胞数和UMAP效果再回来微调。这个跑一遍—看结果—回来调的循环在单细胞项目里极其正常不要觉得麻烦。第二质控日志一定要记。具体包括过滤前多少细胞、过滤后多少细胞、三个指标的阈值各是多少、当时看了哪些图、为什么这么定。不要高估自己的记忆力。一篇论文从投稿到返修可能隔了大半年如果到时候审稿人问你你的过滤标准是怎么确定的而你连跑时的参数都忘了那场面会很尴尬。我现在每个项目都会留一个quality_control_notes.md文件随手记录当时的判断逻辑。第三不同样本之间不要套用同一个阈值。哪怕两个样本来自同一个病人的同一类组织只要批次不同、建库时间不同就应该各看各的分布再分别定阈值。把两个样本合并成一个对象后统一过滤是新手最容易犯的错误之一它会直接导致高质量的样本被你误伤。第四质控的上限不只是不滤错还要考虑下一步需要什么。比如你后面想做拟时序分析那G0/G1期的细胞就要尽量保留如果后面想做转录因子活性分析转录本种类不够丰富的细胞可能支撑不起结论。每一种下游分析对数据质量的要求都不同。质控参数不是越严格越好而是越契合你的分析目标越好。单细胞转录组质控这条路说简单也简单就是看三个指标、设三个阈值说复杂也复杂因为每个数据集都有它自己的脾气。如果你把这篇内容看懂了再回到自己的数据上去实践几轮应该能明显感受到精准质控和套参数跑流程之间的差距。希望这些经验能帮你少走一点弯路。