单细胞数据分析必备：Scrublet多细胞检测与Python模块实现实战

去年处理一个10x单细胞样本的时候聚类注释做到一半卡住了有一群细胞同时高表达T细胞的CD3D和髓系的LYZ怎么看都不像任何一种已知细胞类型。当时第一反应是数据有污染排查到最后才发现罪魁祸首是多细胞doublet——同一个液滴里包进了两个细胞测出来就变成了“四不像”。从那以后我把多细胞去除正式做成了单细胞数据分析流程里固定的一步并用Python封装成了一个小模块。这篇文章就把整个模块的设计思路、实现代码和踩坑经验一次讲清楚适合用Scanpy做单细胞数据分析、但多细胞清理还停留在“听说过”层面的同学参考。1. 单细胞数据里的“多细胞”问题为什么必须专门处理1.1 多细胞是什么又是怎么混进去的单细胞转录组测序scRNA-seq现在的标配是10x Genomics这类微液滴平台。它的基本原理是把单个细胞、凝胶珠和反转录试剂一起包进一个微小液滴里每个液滴理论上只应该捕获一个细胞这样测出来的转录本就能对应到单个细胞。但现实中没有那么理想——液滴包裹是一个随机过程偶尔会有一个液滴同时装进两个甚至多个细胞的情况这就是“多细胞”doublet。多细胞通常有三个来源一是上样细胞浓度太高细胞被液滴随机捕获时“撞车”概率增加二是细胞悬液里有聚团团块被当成一个颗粒包进去了三是液滴生成过程中的偶然误差。很多做数据分析的同学一听到“多细胞”就以为是小概率事件不值得专门处理但实际上这个比例并不低。以10x官方技术文档给出的参考来看多细胞率通常按回收细胞数估算大致是每回收1000个细胞约有0.4%~0.8%的多细胞不同试剂盒版本略有差异。也就是说一个回收了5000个细胞的样本保守估计也有2%~4%的细胞是多细胞。单独看比例不大但在10万个细胞的项目里就是两三千个污染细胞。1.2 不清理会怎样它不是“噪音”而是会制造假信号多细胞最麻烦的地方在于它的转录组是两种甚至多种细胞的表达谱直接叠加。假设一个T细胞和一个巨噬细胞被包在同一个液滴里测出来的结果会同时出现T细胞的markerCD3D和巨噬细胞的markerLYZ。这种混合表达谱会带来一系列连锁问题。第一是聚类结果出现无法解释的“中间态”细胞群。正常的细胞类型聚类应该是界限清楚的而多细胞由于表达谱介于两种细胞之间往往会在UMAP上聚成一团或散落在多种类型交界处。很多新手分析到这里就开始怀疑自己的数据有问题实际上只是没有做多细胞去除。第二是marker基因鉴定被污染。如果你把一个多细胞簇当成新的细胞亚型去注释找出来的“特异marker”实际上只是两个细胞类型marker的混合后面所有基于这个注释的差异表达分析、富集分析都会跟着出错。第三是拟时序分析被严重干扰。拟时序算法会把表达谱连续的细胞串联起来多细胞那种“既像A又像B”的谱型很容易被当成发育过渡态制造出根本不存在的分化轨迹。在有明显祖细胞/终末细胞结构的数据里多细胞造成的伪轨迹尤其隐蔽因为看起来非常“合理”。1.3 为什么常规QC过滤不掉多细胞有人可能会问多细胞的总转录本量不是更高吗用常规质控比如过滤总counts过高、基因数过高的细胞不就能去掉这个问题我一开始也这么想过实际试过之后发现行不通。常规单细胞QC主要过滤三类细胞总counts过低的空滴、基因数过少的破损细胞、线粒体比例过高的濒死细胞。多细胞在总counts上确实偏高但它往往高得并不离谱。尤其当两个被包在一起的细胞转录本本来就少叠加之后可能只相当于一个中等表达水平的正常细胞完全落在QC过滤阈值之内。反过来巨噬细胞、浆细胞这类本身转录本就很丰富的真实单细胞总counts天然就比其他细胞高如果你为了滤除多细胞把counts上限卡得很死会误杀一大批高表达量的正常细胞。所以想靠“表达量高就删掉”的思路处理多细胞基本是捡了芝麻丢西瓜。真正可行的做法必须从表达谱的“混合特征”入手。这也引出了多细胞去除模块的核心思路不是看单个细胞的表达量高低而是看它的表达谱是否像是两个已知细胞类型的“拼接体”。2. 主流多细胞检测思路与Python生态选型2.1 三条技术路线模拟混合邻居分类是最容易上手的目前主流的多细胞检测工具算法上大致能分成三条路线。第一条是“模拟混合邻居分类”路线代表就是Scrublet。它的思路非常直观既然真实多细胞是两个细胞的counts求和那我干脆随机挑选两两细胞把它们的counts矩阵相加人为制造出一批“已知的多细胞”当作训练参考。然后基于真实细胞和模拟多细胞的表达特征构建细胞间的k近邻图计算每个细胞到其邻居的距离。真实单细胞通常能在图中找到大量紧密邻居而多细胞因为表达谱是混合的它到任何单一类型细胞的邻居距离都会偏大不容易形成紧密聚集。这个距离特征经打分后就是每个细胞的doublet score。这个思路的好处是无需真实多细胞的标签纯靠数据内部结构就能学计算开销也可控。第二条是“人工邻域”路线代表是DoubletFinder。它先做PCA降维和聚类构造人工多细胞然后统计每个细胞附近人工多细胞的比例pANN比例高的判为多细胞。这个思路最早是R生态的标配效果不错但调参比较繁琐尤其是pK参数在不同数据集上波动很大需要反复grid search不太适合纯Python流水线。第三条是“深度生成模型”路线代表是Solo、scDblFinder。它们本质上是训练一个神经网络来区分真实细胞和模拟多细胞对复杂数据的表达模式捕捉能力更强但训练时间更长可解释性也差一些而且对电脑配置有一定要求。2.2 Python生态里到底选哪个直接给一张对比表把我实际用过的几个工具列一下工具语言算法核心优势主要注意点ScrubletPython模拟多细胞kNN打分轻量、与Scanpy配合好、速度快对复杂连续发育数据会有一定误判DoubletFinderR模拟多细胞pANN邻域比例在R生态成熟、处理批次设计有额外考量需要调pK参数敏感进入Python流程要引入rpy2scDblFinderR模拟多细胞集成分类阈值自适应较强同样是R包Python调用麻烦SoloPython深度生成模型复杂数据表现稳健训练慢需要GPU更舒服解释性弱如果项目整体是用Scanpy在Python里做的我的建议很明确用Scrublet作为主力工具。原因很简单它可以直接吃AnnData里的counts矩阵输出doublet score和预测标签整个流程不需要跨语言调用出报告也方便。Solo可以留着做备份验证在Scrublet结果存疑的时候跑一遍交叉确认。至于DoubletFinder除非团队里已经有一整套R流程否则不建议专门为了它去折腾跨语言环境。2.3 一个需要理解的关键点为什么“模拟多细胞”能代表真实多细胞Scrublet这类工具能成立有一个隐含假设模拟多细胞的表达谱分布和真实多细胞足够接近。这个假设在大部分样本里是成立的因为真实多细胞本质上就是两个随机单细胞转录本的相加而Scrublet也是随机抽两个细胞做counts相加两者的统计规律几乎一致。但也有例外。如果样本里某个细胞类型占比特别高比如肿瘤组织里T细胞占了60%那么随机抽两个细胞相加时大概率模拟出的是“TT”多细胞但真实多细胞里同样会有大量“T肿瘤细胞”“T巨噬细胞”这类跨类型组合。为了覆盖这些情况Scrublet专门提供了一个参数sim_doublet_ratio用于控制模拟多细胞的数量。在细胞组成比较复杂的样本里我通常会把默认的2.0调高到3.0甚至4.0确保模拟出来的多细胞类型多样性足够。这也是后面模块设计里一个比较重要的参数。3. 把多细胞去除做成一个可复用Python模块3.1 模块的定位输入什么、输出什么、放在流水线哪一步既然要做一个“多细胞去除的模块”就不能只是一段跑完就扔的脚本而是要设计成可以反复调用、能对接不同项目的小工具。先说清楚它的边界。我的设计思路是输入一个已经做完基础QC过滤的AnnData对象里面需要有一个obs列用于区分不同样本/文库比如sample模块会对每个样本分别做多细胞检测避免批次效应干扰输出则在原来的AnnData对象上新增两列一列是doublet_score连续打分一列是predicted_doublet布尔值同时返回过滤后的AnnData。整个模块只负责“检测并去除多细胞”不掺和归一化、聚类、注释这些步骤。放在流水线的位置很关键。我推荐放在基础QC之后、归一化之前。原因有两个第一Scrublet需要在counts空间做模拟混合归一化之后的log数据破坏了counts的数值含义直接拿来模拟会失真第二基础QC先把严重破损、线粒体占比过高的细胞去掉可以避免这些异常细胞干扰后续的邻居距离计算。但要注意这一阶段只能做轻度QC不要用总counts上限去卡细胞否则会把真正的高表达量细胞误杀这个问题后面第5章还会展开讲。3.2 完整代码DoubletCleaner类的实现下面给出我目前在实际项目中使用的模块实现。为了照顾不同数据规模我做了两个关键设计一是按样本分批次检测二是跳过已经跑过Scrublet的样本通过检查obs中的列是否存在方便重复运行。import numpy as np import pandas as pd import scanpy as sc import scrublet as scr import matplotlib.pyplot as plt from scipy import sparse class DoubletCleaner: 多细胞去除模块 参数 ---- adata : AnnData 已做完基础QC的counts表达矩阵需确保 .X 为原始counts batch_key : str 用于按样本/文库分组的obs列名建议填样本编号 expected_doublet_rate : float 预期的多细胞比例按上样细胞数和平台文档估算 def __init__(self, adata, batch_keysample, expected_doublet_rate0.06): self.adata adata self.batch_key batch_key self.expected_doublet_rate expected_doublet_rate self.doublet_score pd.Series(np.nan, indexadata.obs_names) self.predicted_doublet pd.Series(False, indexadata.obs_names) def _prepare_input_matrix(self, sub_adata): 把AnnData的X转成Scrublet需要的矩阵格式 X sub_adata.X if sparse.issparse(X): # Scrublet接受稀疏矩阵但内部一些操作对稠密矩阵更友好 # 细胞数不大时直接转稠密细胞数很大时保留稀疏 if sub_adata.n_obs 50000: X X.toarray() return X def detect(self, sim_doublet_ratio2.0, min_counts3, min_cells3, min_genes200, verboseFalse): 对每个样本分别执行多细胞检测 参数 ---- sim_doublet_ratio : float 模拟多细胞数量与真实细胞数量之比复杂样本建议调高 min_counts / min_cells / min_genes : int 传递给Scrublet内部的基因/细胞过滤参数 if self.batch_key not in self.adata.obs.columns: # 如果没有分样本信息就整体跑一遍 self.adata.obs[self.batch_key] all for batch in self.adata.obs[self.batch_key].astype(str).unique(): idx self.adata.obs[self.batch_key].astype(str) batch sub self.adata[idx].copy() if sub.n_obs 100: print(f[warning] {batch} 细胞数过少跳过多细胞检测) continue X self._prepare_input_matrix(sub) scrub scr.Scrublet( X, expected_doublet_rateself.expected_doublet_rate ) doublet_scores, predicted scrub.scrub_doublets( min_countsmin_counts, min_cellsmin_cells, min_genesmin_genes, sim_doublet_ratiosim_doublet_ratio, verboseverbose, ) self.doublet_score.loc[sub.obs_names] doublet_scores self.predicted_doublet.loc[sub.obs_names] predicted self.adata.obs[doublet_score] self.doublet_score.astype(float) self.adata.obs[predicted_doublet] self.predicted_doublet.astype(bool) return self.adata def filter_doublets(self): 依据predicted_doublet列过滤多细胞返回新的AnnData return self.adata[~self.adata.obs[predicted_doublet]].copy() def plot_score_distribution(self, bins50, saveNone): 画出每个样本的doublet score分布辅助检查阈值 n_batch self.adata.obs[self.batch_key].nunique() fig, axes plt.subplots(1, n_batch, figsize(5 * n_batch, 4)) if n_batch 1: axes [axes] for ax, (batch, group) in enumerate( self.adata.obs.groupby(self.batch_key, observedTrue) ): scores group[doublet_score].dropna() ax.hist(scores, binsbins, color#4C72B0, alpha0.7) ax.set_title(f{batch}\nn{len(scores)}) ax.set_xlabel(doublet score) ax.set_ylabel(cell count) fig.tight_layout() if save: fig.savefig(save, dpi150, bbox_inchestight) return fig def summarize(self): 输出多细胞去除的摘要统计 total self.adata.n_obs detected int(self.adata.obs[predicted_doublet].sum()) ratio detected / total * 100 print(f总细胞数: {total}) print(f预测多细胞数: {detected} ({ratio:.2f}%)) return {total: total, detected: detected, ratio: ratio}3.3 集成进标准Scanpy流程的调用方式模块写完之后在实际Scanpy流程里调用非常直接。下面是一个典型的10x数据入口示例import scanpy as sc from doublet_cleaner import DoubletCleaner # 读取10x的h5文件假设里面有多个样本已经给obs加了sample列 adata sc.read_10x_h5(sample_filtered_feature_bc_matrix.h5) adata.obs[sample] sample_01 # 第一步基础QC只做轻量过滤 adata adata[adata.obs[total_counts] 500].copy() adata adata[adata.obs[n_genes_by_counts] 200].copy() # 第二步多细胞检测与去除 cleaner DoubletCleaner( adata, batch_keysample, expected_doublet_rate0.06, ) adata_with_label cleaner.detect(sim_doublet_ratio2.5) cleaner.plot_score_distribution(savedoublet_score_dist.png) cleaner.summarize() # 过滤掉多细胞后继续常规流程 adata_clean cleaner.filter_doublets() sc.pp.normalize_total(adata_clean, target_sum1e4) sc.pp.log1p(adata_clean) sc.pp.highly_variable_genes(adata_clean, n_top_genes2000) sc.pp.pca(adata_clean, n_comps50) sc.pp.neighbors(adata_clean) sc.tl.umap(adata_clean) sc.tl.leiden(adata_clean)注意一个细节detect()方法会在原来的adata上新增doublet_score和predicted_doublet两列你不需要再去合并结果。跑完之后如果想核对具体的多细胞barcode列表直接取adata_clean.obs_names的补集就行。4. 阈值怎么设才靠谱双峰分布、auto_threshold和人工确认4.1 不要盲信auto_thresholdScrublet在scrub_doublets()运行时会自动算出一个阈值用来把continuous score切成“多细胞/非多细胞”。很多教程到这里就结束了仿佛这个阈值自带权威性。但实际项目中auto_threshold只能作为起点不能作为终点。之所以这么说是因为auto_threshold的算法本质上是在最小化假阳性和假阴性之间的一个平衡。当数据分布比较理想、真实细胞主峰和模拟多细胞峰明显分离时auto_threshold表现很好但当细胞类型高度连续、测序深度不均匀或者存在较强批次效应时score分布可能变成一个大土包这时候auto_threshold硬切出来的阈值要么漏掉大批多细胞要么误杀一片正常细胞。所以我强烈建议每次跑完多细胞检测第一件事是画score分布直方图。理想情况下应该看到两个峰左侧是真实细胞右侧是模拟多细胞也就是被识别出的多细胞阈值落在两个峰之间的低谷。如果没看到这个双峰结构说明数据可能存在问题需要进一步排查而不是直接拿auto_threshold的结果去过滤。4.2 两种有效的人工确认策略看完分布之后还需要从生物学角度做二次确认。我常用的有两个策略。第一个是marker基因共表达检查。多细胞最典型的特征是跨谱系marker同时高表达。比如你之前注释出了T细胞、B细胞、髓系细胞那么预测为多细胞的细胞里应该能观察到“CD3D CD79A”“CD3D LYZ”这类跨系别共表达。可以在Scanpy里直接按预测标签分组看marker表达markers [CD3D, CD79A, LYZ, CD68, NKG7] sc.pl.dotplot(adata, var_namesmarkers, groupbypredicted_doublet)第二个是聚类联动验证。先不做多细胞去除正常跑完聚类和UMAP然后看predicted_doublet细胞在UMAP上的分布。如果多细胞细胞集中出现在两类细胞交界处或者聚成一个独立的混乱簇那基本可以确认检测结果是对的。如果预测多细胞分散在几乎所有簇里均匀地占据5%~10%的空间那就要警惕auto_threshold是不是定得太低了。4.3 参数调整的基本逻辑和常见选择根据我处理几十个样本的经验参数调整主要围绕三个旋钮。第一个是expected_doublet_rate。这个值尽量根据上样回收细胞数来定。10x官方文档对不同版本试剂盒给出的比例略有差异常见说法是每回收1000个细胞对应约0.4%~0.8%的多细胞你按自己的回收细胞数简单乘一下就能得到一个先验值。比如回收5000个细胞0.06左右是一个比较稳的起点如果上样浓度偏高或者拿到了很久的旧数据可以放宽到0.08甚至0.10。第二个是sim_doublet_ratio。复杂样本建议不低于2.0。什么叫复杂样本肿瘤组织、发育组织、脑组织这种细胞类型多且比例悬殊的样本2.0的模拟量不够覆盖跨类型多细胞的多样性我会直接调到3.0或4.0。不过要注意这个值调得越大计算时间越长内存占用也越高。第三个是过滤强度。min_counts3、min_cells3这类参数会先对基因做一轮过滤相当于轻量清洗。如果你前面QC阶段已经做了同样的过滤这里保持一致即可如果你希望Scrublet在更完整的数据上运行可以把min_counts设为1。5. 实战中反复踩过的坑与排查思路5.1 全样本合并跑 vs 按样本分跑这是最容易踩的坑没有之一。Scrublet的doublet score是基于全局邻居图计算的如果你的AnnData里有多个样本直接合并跑不同样本之间的测序深度差异、批次效应都会混进邻居关系里导致score基线完全错位。我见过一个项目合并跑之后某个低测序深度样本的所有细胞都被判成多细胞另一个人高表达量的样本反而一个都没检测出来。正确做法就是模块里写的按样本分跑。每个样本独立模拟多细胞、独立设定阈值最后再拼回AnnData。如果同一个样本还分了多个捕获泳道比如CTRL_A、CTRL_B也建议把它们分别作为batch不要合并。这一步虽然听起来多此一举但在多批次大项目里能省下大量返工时间。5.2 高表达量细胞和低RNA细胞容易被误判Scrublet的判定逻辑依赖邻居距离而高表达量细胞巨噬细胞、浆细胞、肝细胞由于转录本总量大在邻居关系上天然“离谁都远”容易被模拟多细胞带偏被判成多细胞。反过来转录本很少的静止T细胞、静息B细胞就算真是多细胞叠加后的表达量也不一定够高容易被漏掉。这个问题没有完美的自动解法。我的处理方式是先不急着删把predicted_doublet按细胞类型分组统计一次。如果发现某类细胞被预测为多细胞的比例明显高于其他类型而且这类细胞本身表达量又很高就先保留它们只删掉那些出现在“跨类型marker混合簇”里的多细胞。耐心一点宁可少删几个也不要误杀一群真实细胞后续聚类注释比少几个多细胞重要得多。5.3 score分布没有双峰、auto_threshold失灵怎么办有些项目跑完Scrubletscore直方图是一个平缓的单峰根本找不到双峰间隙。这种情况我遇到过大概三分之一。原因往往是数据里的细胞类型呈连续分化状态比如发育中的造血细胞、神经元分化轨迹它们的真实单细胞表达谱本身就介于两种状态之间和模拟多细胞的谱型差异变小双峰自然就消失了。遇到这种情况我一般做三件事第一把sim_doublet_ratio调大到3.0以上增加模拟多细胞的多样性第二把min_genes这类过滤参数放宽避免过滤掉过多低RNA细胞导致邻居关系失真第三用Solo或者DoubletFinder做一次交叉验证看两个工具预测结果的重合度。如果重合度很低说明数据里的doublet信号本身较弱与其硬切阈值不如回到marker共表达检查手工识别明确的多细胞簇去除。5.4 先过滤低质量细胞还是先检测多细胞这个问题我纠结过很久现在固定下来的顺序是先做基础QC再做多细胞检测。基础QC指的是去掉总counts过低、基因数过少、线粒体比例过高的明显破损细胞。这些细胞的表达谱残缺不全如果让它们进入多细胞模拟环节会被随机采来和其他细胞混合制造出一堆虚假的“模拟多细胞”显著抬高整个数据的噪声。但千万注意这个阶段的过滤只做“保底”不要做“上限”。我见过有些流程会在多细胞检测前顺手把总counts排名前5%的细胞删掉理由是“多细胞表达量高”。这个操作等于把最可能含有真实多细胞的高表达量细胞直接抹掉了剩下的数据自然检测不出几个多细胞看似“干净”实际上是把问题藏起来了。5.5 对结果做最终审查UMAP可视化与比例报表最后不管前面参数调得多漂亮多细胞去除这一步一定要有可视化输出和数字报表方便你自己审也方便后续复查。我在模块里已经写了plot_score_distribution()和summarize()实际项目中还会额外加一张UMAP图把predicted_doublet细胞标出来看空间分布。adata.obs[predicted_doublet] adata.obs[predicted_doublet].astype(str) sc.pl.umap(adata, color[predicted_doublet, doublet_score], palette[#4C72B0, #C44E52], size12, showFalse)总结一下我实际项目中的操作习惯多细胞去除和线粒体QC一样已经成了每次单细胞分析的默认定式步骤。它不能直接让你的聚类变得“完美”但它能避免你在注释阶段面对一堆“四不像”细胞群时抓耳挠腮。把这个模块固化下来之后我每次拿到新数据只需要改样本名、改预期多细胞率剩下的交给代码跑然后看一眼score分布和UMAP图就行。这套流程我用了大半年几十个样本下来基本没有再因为多细胞问题返工过。如果你也在Python单细胞流程里为doublet头疼建议直接把这个模块抄进项目里跑一遍你就知道它有多香。