单细胞数据拿到手第一件事往往不是跑聚类而是先想清楚一个问题手里这一管细胞到底有多少是“一个液滴里掉进了两个细胞”的产物这类事件在单细胞建库过程中几乎无法避免如果不去处理后续的差异分析、细胞类型注释都可能被这些“缝合怪”带偏。这篇就围绕python生态里做多细胞去除的模块展开从原理讲到实操再把踩过的坑和判断依据一并说清楚。1. 多细胞到底是什么为什么必须处理1.1 一个液滴里掉进两个细胞是怎么发生的先解释一下背景。以10X Chromium这类基于微流控的平台为例建库时是把单个细胞、凝胶微珠和反应试剂一起包裹在油滴里理论上一个液滴对应一个细胞但实际上是个泊松分布低浓度加载时大部分液滴是空的小部分液滴会同时包进两个甚至更多个细胞。这个概率是固定的细胞悬液浓度越高多细胞比例越高。业界常用的估算方式是大约每多加载1000个细胞doublet率增加0.8%左右所以8000个细胞上机大概会有6%到8%的液滴含有两个细胞。这些双细胞液滴在测序后会被当成一个“细胞”进入表达矩阵它们的转录组是两个细胞表达谱的叠加。问题在于这个叠加不是简单的平均值而是两个细胞mRNA的混合后续标准化、聚类、找marker基因时这种混合信号会把结果搅乱。尤其是异型doublet比如一个巨噬细胞和一个T细胞混在一起聚类时可能形成一群“什么都表达”的过渡状态细胞注释时特别难处理经常会误判成一个新的细胞亚群。同型doublet更隐蔽两个同类细胞混在一起表达谱和正常的单细胞放在一起几乎看不出来通常只能通过表达量倍数异常来隐约察觉。好在实际影响最大的还是异型doublet它对下游细胞类型注释的干扰是毁灭性的。1.2 为什么不能只靠表达量阈值过滤很多人第一反应是doublet的UMI数和基因数肯定比单细胞高一倍直接设个阈值过滤不就行了我在早期做分析时也这么干过实际效果并不理想。原因有两个。第一doublet不一定都是高表达如果两个细胞本身的转录本都很少混合后可能刚好落在正常范围内特别是低质量细胞和高质量细胞形成的doublet总表达量可能跟一个中等质量的单细胞完全一样。第二巨噬细胞这类大转录组细胞在doublet中往往占主导总UMI可能非常高但另一个细胞的信号被完全掩盖这种非对称的doublet用阈值识别也不可靠。所以成熟的方案都是用专门的算法来检测而不是靠简单粗暴的过滤。检测的核心逻辑不是看表达量高低而是看一个细胞的表达谱是否是“两个表达谱的混合体”这就引入了python生态里几个专用模块。2. python生态里可选的多细胞去除模块2.1 主流工具横向对比目前python环境下常用的多细胞检测模块有三类分别是Scrublet、DoubletDetection和SOLO。它们在原理和使用场景上差别比较大先放一个对比表工具语言环境核心原理输出结果适用特点ScrubletPython模拟doublet 近邻分类打分doublet_score、predicted_doublets速度快、参数直观、生态最成熟DoubletDetectionPython迭代聚类 分类器投票每个细胞的标签纯python实现但维护较少SOLOPython (scvi-tools)变分自编码器 半监督分类doublet概率数据噪声大时更稳但训练耗时还有个DoubletFinder是R语言的虽然很多人用但既然标题强调python这里重点说python工具链。从我的实际使用体验来说Scrublet是首选它在运算速度和结果可解释性上表现最好而且和scanpy配合非常顺滑。2.2 Scrublet的核心原理拆解Scrublet的做法从原理上看其实不复杂但设计得很巧妙。它先根据真实数据构造一批人工doublet做法是随机抽取两个细胞的表达向量求和得到一个合成的表达谱。这里有个关键点它并不是全局随机配对而是在每个细胞的一定邻域内配对模拟这样可以保证模拟出的doublet在统计特征上更接近真实情况。有了真实单细胞和模拟doublet两组样本后它通过k近邻图构建分类器对每个细胞计算一个doublet score这个分数越高代表该细胞的表达谱越像doublet。这个方案有两点我很欣赏。第一它不需要任何先验的marker基因完全数据驱动适用性广第二它输出的是连续分数而不是一刀切的标签方便用户根据数据分布自定义阈值。实际使用时阈值的选择非常影响最终结果这个后面在实操部分详细讲。2.3 什么时候可以试试SOLOScrublet虽然不是万能的但绝大多数场景下够用了。如果跑出来的结果不理想比如doublet score的分布很混乱真实doublet和单细胞区分不明显这时候我会换SOLO再跑一次。SOLO是构建在scvi-tools生态里的深度学习方法它用变分自编码器先把转录组压缩到低维潜空间再做doublet分类判断。由于有深度模型的非线性拟合能力SOLO对表达谱噪声高、dropout严重的样本适应性更好代价是训练时间长很多而且对显卡要求更高。我的经验是常规10X数据Scrublet足够如果是FLAIR-seq这类高噪声数据才值得上SOLO。3. Scrublet实战从原始矩阵到干净数据3.1 数据准备和基本质控在跑Scrublet之前最好先做一轮基础质控把明显不合格的细胞清掉。这不只是为了减少计算量更重要的是如果输入数据里存在大量低质量细胞Scrublet模拟doublet的背景也会被污染阈值判断会失准。我一般用scanpy完成这一步骤import scanpy as sc import scrublet as scr import numpy as np # 读取10X的h5文件或者自己构建adata adata sc.read_10x_h5(filtered_feature_bc_matrix.h5) adata.var_names_make_unique() # 计算质控指标 adata.var[mt] adata.var_names.str.startswith(mt-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, inplaceTrue) # 基础过滤基因数太少或太多线粒体比例过高 adata adata[adata.obs.n_genes_by_counts 200, :] adata adata[adata.obs.n_genes_by_counts 6000, :] adata adata[adata.obs.pct_counts_mt 20, :].copy()这里的阈值要结合数据灵活调。比如某些组织类型线粒体比例普遍偏高硬套20%可能会误杀一堆细胞。看过分布再定阈值是最稳妥的。3.2 转置这个细节新手最容易掉坑接下来就是Scrublet的核心环节。这里有一个最常见的坑Scrublet的输入矩阵格式是基因乘以细胞也就是每一行是一个基因每一列是一个细胞。而scanpy中adata.X的默认方向是细胞乘以基因每行是一个细胞。如果不做转置直接丢进去程序会报shape不匹配的错误或者跑出来的结果完全不可信。我自己早期就因为这个浪费了半小时排查。正确的写法是# counts_matrix要求是 genes x cells counts_matrix adata.X.T scrub scr.Scrublet(counts_matrix, expected_doublet_rate0.06) doublet_scores, predicted_doublets scrub.scrub_doublets( min_counts2, min_cells3, min_gene_variability_pctl85, n_prin_comps30 ) # 把结果存回adata对象 adata.obs[doublet_scores] doublet_scores adata.obs[predicted_doublets] predicted_doubletsexpected_doublet_rate这个参数我一般会先算一下。按照经验公式加载1万个细胞时doublet率大约在8%左右如果上机目标是回收1万个细胞就填0.08。如果样本情况比较复杂可以填得稍微高一点让算法更敏感宁可后续多过滤一些也不要放过异型doublet。3.3 阈值怎么选才能不过度杀戮Scrublet跑完后会输出每个细胞的doublet score同时算法会自动给出一个最佳阈值。程序里会生成两个可视化图表一个是doublet score的直方图一个是基于UMAP嵌入的双色标记图。自动阈值并不是永远合理的。我在实际操作中见过两种情况第一种是直方图呈现典型的双峰分布左边是单细胞峰右边是doublet峰这时候自动阈值基本可靠第二种是只有一个平台期或者山脊状分布没有明显分界这时候自动阈值往往会定得偏高导致detection rate很低大量真实doublet被漏掉。遇到这种情况我的做法是手动调低阈值直到在UMAP上看到两边有合理的分离。# 查看阈值 print(scrub.threshold_) # 手动调整阈值例如比自动阈值稍低 doublet_scores adata.obs[doublet_scores].values manual_threshold 0.25 predicted doublet_scores manual_threshold adata.obs[predicted_doublets] predicted这里我特别想强调的是阈值一定不是越小越好。太激进的过滤会把边界区细胞尤其是那些转录组丰度高的巨噬细胞、正在增殖的细胞全部误杀。一个比较实用的技巧是把doublet score的分布按样本或者按批次分别画直方图因为不同批次的多细胞率不一样统一阈值对某些批次太严对另一些又太松。3.4 过滤后的数据要重新做PCA和聚类过滤完doublet后一定要重新跑一遍标准化、PCA、UMAP和聚类不要直接在旧的降维结果上做后续分析。因为去除了一批细胞后特征基因的选择和主成分的计算都会变化特别是如果doublet占了相当比例它们的信号会在计算过程中拉偏主成分方向影响所有细胞的低维表示。# 过滤 adata_clean adata[~adata.obs[predicted_doublets], :].copy() # 重新标准化 sc.pp.normalize_total(adata_clean, target_sum1e4) sc.pp.log1p(adata_clean) sc.pp.highly_variable_genes(adata_clean, min_mean0.0125, max_mean3, min_disp0.5) adata_clean adata_clean[:, adata_clean.var.highly_variable] sc.pp.scale(adata_clean, max_value10) sc.tl.pca(adata_clean, n_comps30, svd_solverarpack) sc.pp.neighbors(adata_clean, n_neighbors10, n_pcs30) sc.tl.umap(adata_clean) sc.tl.leiden(adata_clean, resolution0.5)注意标准化要在过滤后的高变基因筛选之前做顺序不能乱。如果不做标准化直接筛选高变基因高表达基因会垄断高变基因列表影响后续聚类分辨率。4. 结果验证和问题排查实录4.1 怎么判断去除效果到底好不好处理完之后最重要的一件事是验证。不要只看doublet比例降下来了就觉得自己做完了要做交叉验证。我常用的验证方式有三个。第一看UMAP上被标记为doublet的细胞的分布位置。如果doublet细胞比较均匀地散布在各个细胞群边缘说明它们和单细胞的区分度还不够高可能存在漏检如果doublet比较集中地聚成一团或两团且这团细胞表达的marker基因比较杂乱那基本可以确定它们是高置信度的doublet。第二用marker基因做交叉验证。对怀疑是doublet的细胞检查它们是否同时表达两个本该互斥的marker比如同一个细胞同时高表达T细胞的CD3D和B细胞的MS4A1这种特征组合在真实单细胞中几乎不可能出现基本就是doublet。第三如果条件允许和DoubletFinder的结果对照。虽然是R包但跑一次比对也很快。两个方法都标记为doublet的细胞是强阳性基本上可以放心删掉只被一个方法标记的则处于灰色地带可以结合marker判断去留。4.2 常见报错和低质量样本的处理实际跑数据时总有各种意外我挑常见的几个说。Shape mismatch报错。上面提过就是忘了转置矩阵检查counts_matrix的维度即可。细胞数量太少导致报错。Scrublet对细胞数量有最低要求如果样本只有几百个细胞运行时会直接报错或者结果很不稳定。这种情况我一般会降低min_gene_variability_pctl或者索性用更保守的QC过滤不做算法检测了。提示expected_doublet_rate设置过高。如果你填了0.5这种离谱的值算法会警告。一般不要超过0.2正常10X数据在0.05到0.1之间。稀疏矩阵内存爆掉。counts_matrix如果直接转成稠密矩阵会非常占内存尤其在大样本上。Scrublet内部其实会做处理但我建议输入用稀疏格式如果自定义读取要注意类型转换。一条很重要的经验是一批样本最好分开跑Scrublet不要全部合并到一起跑。原因在于不同样本的doublet率不同加载浓度也不同mixed在一起会导致模拟分布错乱结果方差很大。正确做法是对每个样本逐个跑然后把双重打分合并回统一的adata里再做下游分析。4.3 什么时候不推荐用模块检测Scrublet确实好用但也不是万能的。有一种情况我会直接放弃算法检测那就是数据本身非常稀疏且细胞量极少时比如只有两三千细胞、还是高dropout的10X数据。此时Scrublet模拟出的doublet和真实低质量细胞在转录特征上高度相似分类器基本分不出来跑了也只会给你一个悬浮在0.5附近的score分布毫无参考价值。这种情况下我宁可多做一轮严格QC把线粒体比例高、基因数极低的细胞清掉靠表达量差距来兜底。还有一类情况也容易误报处于连续分化轨迹上的细胞比如造血系统里从造血干细胞逐步向各谱系分化的中间态细胞它们的表达谱本身就是过渡态和doublet模拟出的“两个谱系混合”信号很像很容易被识别成doublet。如果关注的是分化轨迹建议先保留这些细胞跑完轨迹分析再根据分化轨迹和Marker验证哪些是真实doublet哪些是过渡态。5. 一些进阶场景和扩展思路5.1 多组学数据和CITE-seq里的doublet处理现在越来越多数据是CITE-seq即同一反应体系里同时检测转录组和表面蛋白。CITE-seq里doublet的判断有一个天然优势蛋白标签的background水平可以辅助验证。实际处理时转录组层面的Scrublet结果可以和蛋白层面的异常值互为印证如果某个细胞转录组score偏高同时蛋白标签出现多标签叠加那基本就是doublet。虽然处理流程上还是先跑Scrublet但验证环节可以多一重保障。5.2 基因表达层面的doublet检测并不可靠有些教程会说DNA拷贝数变异也能用来检测doublet这在肿瘤数据里确实有应用比如inferCNV可以看拷贝数异常。但用python生态的通用流程跑非肿瘤数据时这个方法意义不大因为正常组织细胞没有拷贝数变异检测不出差异。不要硬套术业有专攻。5.3 把Scrublet整合进Snakemake流程如果样本量很大比如几十个样本要批量处理建议把Scrublet封装成函数写进Snakemake或Nextflow流程里。每个样本独立执行输出doublet特征和标记结果方便统一追溯。我自己的一个简化版本是这样组织的# pseudo code for batch processing def run_scrublet(adata_sample, expected_rate): counts adata_sample.X.T scrub scr.Scrublet(counts, expected_doublet_rateexpected_rate) scores, pred scrub.scrub_doublets() adata_sample.obs[doublet_scores] scores adata_sample.obs[predicted_doublets] pred return adata_sample批量跑的时候建议把每个样本的score直方图保存下来之后统一核对看看有没有批次效应导致score分布偏移。这一步对大规模项目尤其关键只报一个细胞数不报分布形态复盘时很难定位问题。6. SOLO作为备选方案的具体操作如果Scrublet在当前数据集上表现确实不好SOLO是一个值得尝试的替代方案。它的代码流程也很清晰基于scvi-tools生态import scvi import torch # 准备数据 scvi.model.SOLO.setup_anndata(adata, batch_keysample) model scvi.model.SOLO(adata) model.train(max_epochs200) # 预测 df model.predict() df[prediction] df.idxmax(axis1) adata.obs[solo_doublet] df[prediction].valuesSOLO输出的是“doublet”和“singlet”两类概率预测结果比较直观。但从我的使用感受来说它在中小数据集上并没有比Scrublet显著占优运算时间却长得多。如果你只是想要一个快速可靠的结果先跑ScrubletSOLO作为交叉验证工具比作为首选更合理。从整体分析流程来看多细胞去除属于单细胞数据处理的中上游环节它做得好不好直接影响后面所有分析的可靠性。所以在这个环节多花点时间把原理搞透、把阈值调对后面做注释做差异分析时能省很多事。我个人在反复跑不同数据集之后最大的体会是doublet detection不是一个“跑完代码就算完成”的步骤它的核心价值在验证环节。只要筛选条件、阈值设定稍有偏差就会向聚类结果里埋下一颗雷往往直到画marker基因图的时候才爆炸。所以每次跑完Scrublet我都会强制自己在UMAP上多花十分钟把doublet的分布、marker的共表达都检查一遍确认这批数据是干净的再放心往下走。