起因：最近有个问题样本，跑完cellranger，样本的cellranger结果如下，细胞数目极高(3W+)。在后续数据质控分析中，线粒体基因占比和双细胞率均很高，用scDblFinder进行双细胞预测，双细胞占率竟然高达34%。我很好奇，双细胞率为什么会这么高，如何审视这个结果？决定看看scDblFinder的细节。
image.png
问题1：如何理解Doublets？
在scRNA-seq的细胞捕获步骤中，两个或多个细胞聚集成单个液滴（双联体/多联体）会导致混合的转录组，也就是两个或多个细胞共用一个barcode，称为doublets或multiplets（后面统称为doublets）。它是基于液滴的单细胞测序的技术副产品。双细胞会造成每个“细胞”的高UMI计数，改变cluster的细胞类型鉴定干扰到下游分析。这会导致对稀有细胞类型、中间细胞状态和疾病相关转录组学特征的人为错误发现。双细胞率已被证明与捕获的细胞数量成正比（Bloom 2018; Kang et al. 2018）。
双细胞可以分为同型（相同细胞类型）或异型（不同细胞类型）。
问题2： Cell Ranger可以自动剔除doublets和multiplets吗?
答：目前没有方法可以识别与双细胞中单个细胞相关的转录本信息。
10X官网对双细胞率的相关答复，我们目前没有一种方法通过算法识别单个物种的单细胞基因表达数据，barcode是否包含多个细胞。
目前，Cell Ranger 软件仅在barnyard实验或多物种实验中估计双细胞率。
对此，10X也给出三条参考意见：
- Infer doublets in a single species case if there are known cell-type specific markers. For example, the presence of T cell and B cell specific markers coming from a single barcode may indicate a GEM with both T and B cells.
- Evaluate the suggested workflow in the R package, Seurat, which has some functionality for flagging cells with a clear outlier number of UMIs or genes detected http://satijalab.org/seurat/pbmc3k_tutorial.html ("QC and selecting cells for further analysis.)
- Evaluate the doublet identification method in the R package DoubletFinder.
他的意思是：1) 通过已知细胞类型的marker基因来鉴别双细胞，比如T/B细胞的marker基因，在同一个barcode细胞中同时高表达就可判定为双细胞；2) seurat标准分析流程，质控环节通过UMI和gene指标过滤；3）运用scDblFinder双细胞预测软件。
问题3：10X Genomics 单细胞实验中估计的双细胞率是多少？
假设不存在细胞结团，可使用下表（取自10X基因组学用户指南）来估计单细胞实验中估计的双细胞率。
Recommended starting point is to load ~1,600 cells per reaction, resulting in recovery of ~1000 cells, and a multiplet rate of ~0.8%. The optimal input cell concentration is 700-1,200 cells/µl
双细胞率（0.8%/1000cells），如果细胞数为1W，双细胞率为7.6%，约8%。
image.png
但是，在scDblFinder的文章中，提到在单细胞实验中，双细胞率为10-20%。
image.png
It is therefore at present common in singlecell experiments to have 10-20% doublets, making accurate doublet detection critical.
单细胞实验双细胞率为10-20%，这个数值明显高于上面10X给出的双细胞率(~8%)。这个怎么理解呢？
我想到的几个原因：
1）10X给出的是理想情况的双细胞率(细胞不结团)，用标准样本做基准比较；
2）双细胞率跟实验环节中的样本处理和细胞上样量都有关系；
3）还取决如何计算双细胞，双细胞率=双细胞数目/总细胞数；因计算的细胞总体不同而不同；
我们一般会进行QC细胞质控，用UMI/genes指标过滤掉低质量细胞和异常值细胞；
如果QC细胞质控后计算双细胞率，和QC细胞质控前计算双细胞率，预测的双细胞率会不一致；
最近，在网上找到一个10x Genomics 提供的估计双细胞率是：
image.png
比如1W个细胞，双细胞率为：0.008*(10000/1000)=0.08=8%。
下面我们看看scDblFinder软件是如何具体执行的。
双细胞在单细胞测序数据中很普遍，可能会导致人为错误的发现。
目前，实验层面还是无法检测同一样本的细胞形成的双细胞，包括异型双细胞。
算法层，已经开发了许多计算方法来根据转录谱识别双细胞。大多数这些方法依赖于通过对真实细胞求和或求平均来生成人造双细胞，并对它们与真实细胞之间的相似性进行评分。 例如，DoubletFinder在真实细胞和人工双细胞的合集上生成k最近邻近图(kNN) ，并估计每个细胞附近人造双细胞的密度（McGinnis, Murrow, and Gartner 2019）。 以类似的方式，Bais和Kostka (2020) 提出的bcds算法和共表达评分cxds算法。
Xi 和 Li (2021a) 最近发表的文章中对双细胞检测方法进行基准测试，使用模拟数据和包含双细胞实验标记的真实数据数据集，发现DoubletFinder的算法最为准确。 但是，基准测试也发现，没有一种方法在所有数据集上都是系统性最优，强调在各种数据集上测试和基准测试方法的必要性，并表明某些算法在不同情况下可能具有优势和劣势。
没有一种算法是完美不缺的，特别是预测模型算法。
scDblFinder achieves the highest mean AUPRC and AUROC values, and it is also the top method in terms of the precision, recall, and TNR under the 10% identification rate.
下图比较了scDblFinder 这个包中一些方法（以粗体显示）与其他方法：
image.png
因此，我们在单细胞转录组数据质控过滤时，会考虑到双细胞的因素，通过相关软件进行预测双细胞。常用的软件有scDblFinder(R语言)和Scrublet(python)。这里仅讨论scDblFinder。
image.png
1.scDblFinder安装
安装scDblFinder需要满足R >= 4.0 和 Bioconductor >= 3.12
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("scDblFinder")
# or, to get that latest developments:
BiocManager::install("plger/scDblFinder")
2.scDblFinder使用
scDblFinder的输入数据是SingleCellExperiment对象（空的drops已经移出），至少要包含counts矩阵(assay ‘counts’)。即sce对象都不应该包含空滴，但不应该经过非常严格的过滤（这会影响双细胞率的估计）。
如果还包含归一化矩阵 (assay ‘logcounts’) 和PCA (reducedDim ‘PCA’)，可以使用scDblFinder的cluster模式（不常见）。
理解SingleCellExperiment对象？
SingleCellExperiment对象是Bioconductor针对单细胞scRNA-seq，定义的通用的标准数据对象。scater和seurat等均可操作此对象。基因表达数据通常储存样本-基因定量表达矩阵。在scRNA-seq分析中，我们通常从计数矩阵开始，即每个细胞的特定特征的reads/UMIs数目。这个特征可以是基因、亚型或外显子等，通常我们在基因水平进行分析。
除了基因量化矩阵，我们还可能有关于基因的信息（例如基因组位置、它们是哪种基因类型、它们的长度等）和细胞的信息（例如它们的起源组织、患者供体、处理批次、疾病状态、治疗暴露等）。我们还可以从原始计数数据中生成其他矩阵，例如归一化计数矩阵。由于单细胞实验数据的维数非常高（有数千个细胞和数千个基因），我们经常采用降维技术来捕捉较低维数据中的主要变化。
SingleCellExperiment（简称 sce）是R语言中特殊的S4对象，以同步方式存储所有这些信息。对象的不同部分被称为“slot”：
- 一个或多个表达式矩阵 - 存储在‘assay’ slot中；
- 关于基因（对象的行）的信息 - 存储在‘rowData’ slot中；
有关细胞（对象的列）的信息 - 存储在‘colData’ slot中；
2.1 单样本
set.seed(123)
library(scDblFinder)
# we create a dummy dataset; since it's small we set a higher doublet rate
sce <- mockDoubletSCE(dbl.rate=0.1)
# we run scDblFinder (providing the unsually high doublet rate)
sce <- scDblFinder(sce, dbr=0.1)
对于 10x 数据，通常将dbr留空是安全的，它会自动估计。 scDblFinder的输出会在sce的colData中添加一些以‘scDblFinder’为前缀的列，其中最重要的是：
- sce$scDblFinder.score : the final doublet score（最终的doublet分数）
- sce$scDblFinder.class : the classification (doublet or singlet)判定的分类(双细胞或单细胞)
table(truth=sce$type, call=sce$scDblFinder.class)
## call
## truth singlet doublet
## singlet 499 1
## doublet 0 34
2.2 多样本
如果你有多个样本（理解为不同的细胞捕获），那么最好为每个样本分别进行双细胞识别（对于cell hashes实验中的多重样本，那意味着每个批次）。 可以通过简单地向scDblFinder的samples参数提供样本 id来完成，或者，将样本信息存储在colData列中，提供列名即可。另外，还可以考虑使用BPPARAM参数对其进行多线程处理（假设有足够的RAM）。 例如：
library(BiocParallel)
sce <- scDblFinder(sce, samples="sample_id", BPPARAM=MulticoreParam(3))
table(sce$scDblFinder.class)
2.3 案例实战
我们用之前案例中的数据测试下scDblFinder函数。
单个样本
rat <- CreateSeuratObject(adj.matrix,project = "pbmc10k")
srat
## Detect Doublets
set.seed(123)
sce <- as.SingleCellExperiment(srat)
sce <- scDblFinder(sce, dbr=0.1)
sce$doublet_logic <- ifelse(sce$scDblFinder.class == "doublet", TRUE, FALSE)
plotDoubletMap(sce)
table(sce$scDblFinder.class)
# singlet doublet
# 9209 985
该样本共10194个细胞，其中968个细胞被预测为双细胞，双细胞率为9.5%
image.png
image.png
问题样本
拿最开始提到的问题样本，不设置dbr，双细胞率为34.1%；设置dbr=0.1，双细胞率为11.4%。很显然，对于这个样本，不同参数下，预测双细胞数值差异很大。
scrna.matrix <- Read10X("data/sample1/filtered_feature_bc_matrix")
seurat_obj <- CreateSeuratObject(scrna.matrix,project = "sample", min.cells = 3, min.features = 200)
dim(seurat_obj)
# [1] 18975 31340
## Detect Doublets
set.seed(123)
sce <- as.SingleCellExperiment(seurat_obj)
sce1 <- scDblFinder(sce, dbr=0.1)
# 3574 (11.4%) doublets called
table(sce1$scDblFinder.class)
# singlet doublet
# 27766 3574
sce2 <- scDblFinder(sce)
# 10698 (34.1%) doublets called
table(sce2$scDblFinder.class)
# singlet doublet
# 20642 10698
那么该如何审视这个结果，选择怎样的参数？
我想到的是，对于一个预测模型来说，调参的意义不大，我们对结果没有预判，修改参数，都会出现不同的预测值。另外，对于一个常规的10X单细胞转录组数据，我们对双细胞率是有一定预判的，10X的实验步骤大致固定，cellranger的细胞数大致1W，双细胞率大致10%，我们知道预测的边界。我们有粗略的“标尺”。
但是现在，cellranger给出的细胞数是3W+，我们其实是不清楚双细胞率的边界，细胞数“超纲”了，只知道细胞数越多，双细胞也会越多。这类样本太少，我们没有横向可参考的实例。如果出现这种结果，最应该审视的是实验端出了什么问题，线粒体基因占比也非常高。
2.4 基于cluster的双细胞识别
scDblFinder有两种生成人造双细胞的主要模式：随机模式（scDblFinder.random，clusters=FALSE, 默认方式）和基于cluster的模式（scDblFinder.clusters，clusters=TRUE 或提供你自定义的cluster - 以前版本的方法）。在实际中，我们观察到两种方法都表现良好（比其他方法要好）。当数据集被分成清晰的cluster时，教程建议使用基于cluster的方法（例如发展轨迹），否则使用随机模式。
suppressPackageStartupMessages({
library(scDblFinder)
library(SingleCellExperiment)
})
sce <- readRDS("GSE162690.SCE.rds")
sce <- scDblFinder(sce, clusters=TRUE)
sce$scDblFinder.raw <- sce$scDblFinder.score
sce <- scDblFinder(sce, clusters=TRUE, aggregateFeatures=TRUE, nfeatures=25,
processing="normFeatures")
CD <- colData(sce)
saveRDS(CD, file="atac.colData.rds")
3.scDblFinder算法的实现
双细胞分为同型双细胞（'Homotypic' doublets）和异性双细胞（ 'Heterotypic' doublets）。同型双细胞由相同类型的细胞（即相似的转录状态）组成，仅根据它们的转录组信息很难辨识。而且，它们对于大多数分析来说也相对无害，因为它们看起来与单细胞高度相似。 相反，异型双细胞（由具有不同转录状态的细胞形成）表现为一种人为的新型细胞类型，会影响下游分析。
scDblFinder只关注异性双细胞。
step1：将数据集缩减到仅高表达的基因（默认为 1000）； 如果使用基于cluster的方法，则会选择每个cluster的表达靠前的基因。另外使用基于cluster的方法（而不是人为指定cluster），将会执行快速聚类（请参阅fastcluster）。
step2：通过组合不同cluster的细胞来创建人工双细胞，创建的人工双细胞数与cluster的数目成比例。 我们主要关注不同cluster间的双细胞，我们不会试图识别同型双细胞，无论如何，它们实际上无法识别且对下游分析相对无害。因此，我们减少了人工双细胞的必要数量， 也防止分类器被训练以识别与单细胞无法区分的细胞（因此将单细胞称为双细胞）。 scDblFinder另一种策略是生成完全随机的人工双细胞，并使用迭代程序从训练中排除无法识别的人工双细胞。 在实践中，这两种方法具有相当的性能，它们也可以结合使用。
step3：然后对真实细胞和人工双细胞的合集进行降维，并生成最近邻网络。 接着使用网络来估计每个细胞的许多特征，特别是最近邻居中人工双细胞的比例。 该比率不是选择特定的邻域大小，该比率是在不同的 k 值下计算的，通过使用多个预测变量创建分类器。预测变量还包括距离加权比，进一步添加了的细胞层面上的预测变量：对主成分的预测；文库大小； 和共表达分数（基于Bais 和 Kostka2020 的变体）。 然后 scDblFinder训练梯度提升树（ gradient boosted trees ），以根据这些特征区分来自真实细胞的人工双细胞。 最后，阈值程序通过同时最小化错误分类率和预期的双细胞率来决定调用细胞的分数（参见Thresholding）。
step4：基于分类器方法的一个关键问题是一些真实细胞被错误标记，从某种意义上说，它们实际上是双细胞，但被标记为单细胞。这些会误导分类器。出于这个原因，分类和阈值处理以迭代方式执行：在每一轮中，从下一轮的训练集中删除识别为双细胞的真实细胞。
Overview of the scDblFinder method.
方法说明
3.1 拆分捕获
双细胞只能出现在给定的样本或某次捕获中，因此需要为每个样本单独进行双细胞判别，这也加快了分析速度。如果给定samples参数，scDblFinder将利用该参数将细胞拆分为单个样本/捕获，并在给出BPPARAM参数的情况下并行分析。分类器将在全局范围内进行训练，但阈值将在每个样本的基础上进行优化。如果你的样品是多标签，即不同的样品混合在不同的批次中，那么需要提供批次信息。
3.2 数据的降维和聚类
通过将数据集减少到仅高表达的基因（由nfeatures参数控制），可以大大加快分析速度，即使会稍微影响到准确度。然后，根据cluster参数，将执行最终的PCA和聚类（使用内部 fastcluster函数）。基于cluster方法的基本原理是同型双细胞几乎不可能根据它们的转录组进行区分，因此创建这种双细胞是一种计算资源的浪费，而且还会误导分类器标记为单细胞。
然而，另一种方法是随机生成双细胞（将clusters设置为 FALSE 或 NULL），并使用迭代方法从训练中排除无法识别的人工双细胞。
3.3 生成人工双细胞
根据cluster和propRandom参数，将通过合并随机细胞和/或不相同的cluster对的细胞合并，形成人工双细胞（这可以使用getArtificialDoublets函数手动执行）。 一部分双细胞将简单地使用组成细胞的counts总和，而其余的将进行调整文库大小和进行泊松重采样，数据校正。
3.4 检查每个细胞的k最近邻 (kNN)
对真实细胞和人工细胞的组合执行新的PCA，从中生成 kNN网络。 使用这个 kNN，为每个细胞收集了许多参数，例如 KNN中双细胞的比例、到最近双细胞和最近非双细胞的距离之比等。在输出中报告了一些带有“scDblFinder.”前缀的功能，例如：
- distanceToNearest: 到最近细胞的距离（真实或人工细胞）
- ratio: KNN 的双细胞比例
- weighted: 为双细胞的KNN的比例，按其距离加权（对孤立细胞有用）
4.其他重要参数
scDblFinder有相当多的参数来控制预处理、双细胞的生成、分类等（参见?scDblFinder）。 我们仅对重要参数进行说明。
4.1 双细胞期望检出率
双细胞的期望检出率对邻域中人造双细胞的密度没有影响，但会影响分类器的分数，特别是分类临界值。是通过dbr和dbr.sd参数指定（dbr.sd指定dbr周围的 +/- 范围，在该范围内与dbr的偏差将被视为空）。
对于10x数据，捕获的细胞越多，产生双细胞的概率越大，Chromium文档表明每1000个细胞捕获的双细胞率大约为 1%（因此对于 5000 个细胞，(0.015)5000 = 250 个双细胞） ，scDblFinder默认的预期双细胞率将设置为0.1（默认标准偏差为 0.015）。但是请注意，不同的实验方案可能会产生更多的双细胞率，因此需要相应地更新。如果不确定双细胞率，您可能会考虑增加 dbr.sd，以便大多数/纯粹从错误分类错误中估计它。
5.高频问题
问题1：我捕获到太多双细胞-这是怎么回事？
那么你很可能有错误的双细胞率。 如果你没有提供dbr参数，双细胞率将使用10X Genomics预期双细胞率自动计算，这意味着捕获的细胞越多，双细胞率就越高。 如果你认为不适用于你的数据，可手动设置dbr。
如果出现意外高的双细胞率最常见原因是，1）你有一个多样本数据集并且没有按样本进行拆分。 scDblFinder会认为数据是具有大量细胞的单次捕获，因此具有非常高的双细胞率。 按样本拆分应该可以解决问题。
阈值根据预期双细胞数量和错误分类（即人造双细胞）试图最小化方差，这意味着有效（即最终）双细胞率将与给定的不同。 scDblFinder还认为假阳性比假阴性对后续的分析问题要小些。你可以通过设置 [dbr.sd=0]在一定程度上减少与输入双细胞率的偏差。
问题2：我应该使用基于cluster的双细胞判别吗？
虽然这两种方法在基准测试中的表现非常相似，但随机生成方法在复杂数据集中通常略胜一筹。 如果你的数据被非常清晰地划分为cluster，或者你对双细胞的起源感兴趣，则基于cluster的方法更可取。 这也将能够更准确地计算同型双细胞率，因此略好于阈值法(thresholding)。 否则，特别是如果你的数据没有非常清楚地划分为cluster，则随机方法（例如clusters= FALSE）更可取。
问题3：cluster没有任何意义!
如果你在多样本数据集上运行scDblFinder但是未提供cluster标签，而是基于特定样本的标签（意味着一个样本中的标签“1”可能与另一个样本中的标签“1”无关），并且 在 tSNE上绘制它们看起来没有意义。 出于这个原因，当运行多个样本时，建议首先将所有样本聚集在一起（例如使用 sce$cluster <- fastcluster(sce)），然后将cluster信息提供给 scDblFinder。
问题4：‘Size factors should be positive’报错
如果某些细胞的读数为零（或非常接近于零），则会出现‘Size factors should be positive’此错误。 过滤掉这些细胞后，错误应该消失了。 但是请注意，我们建议在运行 scDblFinder之前不要进行过于严格的过滤。
问题5：我怎样才能重现我的结果？
由于它依赖于人工双细胞的部分随机生成，因此对同一数据多次运行scDblFinder会产生略有不同的结果。你可以使用 set.seed() 确保可重复性，但是在多线程时使用 set.seed() 是不行的。请使用以下程序：
bp <- MulticoreParam(3, RNGseed=1234)
bpstart(bp)
sce <- scDblFinder(sce, clusters="cluster", samples="sample", BPPARAM=bp)
bpstop(bp)
问题6：我可以将它与Seurat或其他工具结合使用吗？
如果输入的sce对象已经包含归一化矩阵( logcounts)或名为“PCA”的reducedDim数据，scDblFinder将使用它们进行聚类分析。 此外，可以使用 scDblFinder() 函数的 cluster参数手动指定。 通过这种方式，seurat聚类可以例如用于创建人造双细胞（参见 ?Seurat::as.SingleCellExperiment.Seurat for conversion to SCE）。
在人造双细胞生成之后，真实和人造双细胞的计数必须一起重新处理（即归一化和 PCA），在内部使用scater执行的。 如果您希望以不同的方式执行此步骤，您可以提供自定义函数来执行此操作（请参阅 ?scDblFinder 中的处理参数）。 然而，我们注意到，这一步的变化对双细胞检测的影响不大。 事实上，例如，根本不执行任何归一化会降低双峰识别的准确性，但也是一点点。
问题7：scDblFinder可以应用到scATACseq数据吗？
可以，专门处理峰值数据。由于单细胞ATAC-seq数据的稀疏性比转录组要大得多，而且scDblFinder需要处理一系列基因，因此使用默认的标准参数，运行的性能较差(执行慢)。因此，我们推荐使用aggregateFeatures=TRUE，这将在正常的scDblFinder 过程之前聚合相似的基因（而不是选择基因），会产生不错的结果。如果基因足够少，我们推荐直接基于距离计算而不是通过SVD步骤获得，如下所示：
sce <- scDblFinder(sce, aggregateFeatures=TRUE, nfeatures=25, processing="normFeatures")
问题8：我是在双细胞鉴定前还是鉴定后运行QC细胞质控？
cDblFinder的输入数据不应包含空液滴，并且可能需要移除覆盖率非常低的细胞以避免错误（例如 <200 reads）。
进一步的细胞质控应该在双细胞识别的下游进行，有两个理由：
1.默认的预期双细胞率是根据给定的细胞计算的，如果你排除了很多质量低的细胞，scDblFinder可能会认为双细胞率应该低于实际值。
2.剔除所有低质量细胞可能会妨碍我们检测由高质量细胞和低质量细胞组合形成的双细胞的能力。
话虽如此，这些主要是理论依据，除非你的QC过滤非常严格（而且不应该如此！），否则结果不太可能有很大的不同。
scDblFinder运行之前要做一些初次过滤，但不要太严格（例如 <200 UMI）
质控粗略过滤->运行scDblFinder->较严格过滤
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# 粗略过滤
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
后记：
之前在群里看到有人问，双细胞质控后，拿质控后的数据重新跑scDblFinder，还是会有大量双细胞被检出？
这个是必然的，由scDblFinder的算法决定，它是由输入数据建立起预测分类模型，不像singleR有另外一套参考数据集，拿输入数据去map到参考数据集。
只要你喂给scDblFinder数据，它都会给输入的细胞一个score值，然后设置阈值进行分类，值高的为双细胞。
scDblFinder强烈依赖输入数据，有它使用的范畴，所以反复进行双细胞scDblFinder质控，用法是不对的。
拿到单细胞转录组数据，先质控粗略过滤->运行scDblFinder进行双细胞质控->较严格质控过滤。
其实最好有一个相互验证的过程。
网友评论