[scRNA-seq]doublets检测——DoubletFinder scrublet （下）

在上一篇文章里我们聊到为什么要进行doublet检测以及DoubletFinder的基本原理。在这篇文章里，我们就来聊聊另一个doublet检测工具——scrublet的基本算法、它和DoubletFinder的使用比较，以及使用这两个工具的代码范例。

上图是上篇文章聊到的DoubletFinder的算法示例图，scrublet的核心思想和DoubletFinder非常相近。它们的主要的区别在于第五步，即对于每一个细胞邻近细胞中人工模拟的doublet细胞分布的描述上。

scrublet的作者认为当我们的数据中出现doublet时，有两种情况。一种称之为Embedded errors，在这种情况下，doublet和真正存在的某种细胞类型有相似的基因表达，doublet会和这些细胞被聚类到一起，同时在分群结果中占某一个群的一小部分，不会对最终的分析结果产生严重的影响。另一种情况称之为Neotypic errors，在这种情况下，doublet会构成一个和现有的细胞类型基因表达非常不同的群，而这个新的群会严重影响到后续的分析结果。但不管在什么情况下，作者都假定doublet只占样本数据中很小的一部分。

基于上面的假设，scrublet对于第五步获得的每个细胞邻近细胞中人为模拟的doublet的分布用一个作者自定义的公式去进行拟合，如下图所示。其中Pobs为观察到的分布，PD代表是doublet的可能性，PS代表是单细胞的可能性，PD和PS的和为1。在算出蓝框内的值之后，我们可以将其与1比较。如果远大于1，则说明doublet的影响要远大于单细胞，即属于Neotypic errors。如果远小于1，则说明单细胞的起主要影响，属于Embedded errors。如果约等于1，则不属于doublet的两种情况，为真正存在的生物学发现。

DoubletFinder和scrublet实操比较

1. scrublet是基于python的。尽管它有非常详尽的tutorial，对于没有python基础的同学，还是有些许挑战。而DoubletFinder则是基于R的，使用R对于很多做scRNA-seq生信分析的同学来说，还是比较自然。

2. 因为scrublt是基于python的，相较于DoubletFinder，它的速度也要快很多。一个10x的样本只需要几分钟到十几分钟就能跑完，而DoubletFinder则需要几十分钟到一个多小时。

3. 最重要的一点是scrublet比DoubletFinder要准确得多。DoubletFinder最终的结果受设定的pANN阈值的影响，结果稳定性不高。小L曾经有一次得到超过90%数据都是doublet的结果。而scrublet的结果则要稳定得多，也符合我们的预期。

DoubletFinder代码范例

总结完毕，下面就是具体的实操了。正如小L在上篇推送中分析的，DoubletFinder的使用应该是针对每一个sample进行单独分析，并在用seurat对数据进行任何处理之前进行。同时因为DoubletFinder非常占用内存，我们需要及时保留我们需要的信息，即每个细胞的预测结果以及pANN值，并把不需要的部分删除。其中大写字母的部分需要大家根据自己的需求进行替换。

library(tidyverse)
library(Seurat)
library(DoubletFinder)sampleV = c('SAMPLE1', 'SAMPLE2', 'SAMPLE3')
dataL = list()
## create seurat object for all samples
for (sample in sampleV) {file_name = str_c('DATA_FOLDER', sample, '_cellranger/outs/filtered_feature_bc_matrix')seuratObj_counts <- Read10X(data.dir = file_name)seuratObj = CreateSeuratObject(counts = seuratObj_counts, project = str_c(sample, "1"))dataL[[sample]] <- seuratObj
}for (idx in 1:length(dataL) {## preprocessingdata <- NormalizeData(dataL[[idx]])data <- ScaleData(data, verbose = FALSE)data <- FindVariableFeatures(data, verbose = FALSE)data <- RunPCA(data, npcs = 30, verbose = FALSE)data <- RunUMAP(data, reduction = "pca", dims = 1:30)data <- FindNeighbors(data, reduction = "pca", dims = 1:30)data <- FindClusters(data, resolution = 0.5)sweep.res.list <- paramSweep_v3(data, PCs = 1:30, sct = FALSE)sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)bcmvn <- find.pK(sweep.stats)## Homotypic Doublet Proportion Estimate -------------------------------------------------------------------------------------annotations <- data@meta.data$ClusteringResultshomotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- sample14@meta.data$ClusteringResultsnExp_poi <- round(0.075*nrow(data@meta.data))  ## Assuming 7.5% doublet formation rate - tailor for your datasetnExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))## Run DoubletFinder with varying classification stringencies ----------------------------------------------------------------data <- doubletFinder_v3(data, PCs = 1:30, pN = 0.25, pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)## save resultsdataL[[idx]]$doubFind_res = data@meta.data %>% select(contains('DF.classifications'))dataL[[idx]]$doubFind_score = data@meta.data %>% select(contains('pANN'))
}

经过上述处理后，dataL就是由所有sample得到的seurat对象构成的list。其中每一个seurat对象的metadata里包含了pANN值（doubFind_score）以及DoubletFinder的预测结果（doubFind_score）。用dataL，我们就可以根据seurat的教程进行后续的filtering、integration等操作了。需要注意的是，如果sample的数量太多，我们可以通过改for循环里idx的范围，每次只处理3-5个sample并及时保存，来避免内存不够等问题。

scrublet代码范例

scrublet的安装可以在这里找到（https://github.com/swolock/scrublet），使用指南则在这里（https://github.com/swolock/scrublet/blob/master/examples/scrublet_basics.ipynb）。需要注意的是scrublet使用的是解压后的matrix.mtx文件和genes.tsv文件，而cellranger产生的是压缩文件，我们需要先提前解压一下。

在使用指南中更改好matrix.mtx和genes.tsv文件对应的位置之后，我们就可以按照里面的步骤一步步跑下来，得到doublet_scores和predicted_doublets。因为后续小L是用seurat对数据进行处理，所以就将doublet_scores和predicted_doublets单独保存下来，再放进seurat对象的metadata里。其中scrublet_res是scrublet的判断结果而scrublet_score是scurblet得到判断所基于的分数。


## 保存doublet_scores和predicted_doublets (python)
d = {'score': doublet_scores, 'prediction': predicted_doublets}
df = pd.DataFrame(data=d)
df.to_csv(DATA_FOLDER2 + sample + '_doublet.csv', index=False)## 导入seurat对象的metadata里 (R)
library(tidyverse)
library(Seurat)
library(DoubletFinder)sampleV = c('SAMPLE1', 'SAMPLE2', 'SAMPLE3')
dataL = list()
## create seurat object for all samples
for (sample in sampleV) {file_name = str_c('DATA_FOLDER', sample, '_cellranger/outs/filtered_feature_bc_matrix')seuratObj_counts <- Read10X(data.dir = file_name)seuratObj = CreateSeuratObject(counts = seuratObj_counts, project = str_c(sample, "1"))doublet_info = read.csv(str_c('DATA_FOLDER2', sample, '_doublet.csv'))seuratObj$scrublet_res = doublet_info$predictionseuratObj$scrublet_score = doublet_info$scoredataL[[sample]] <- seuratObj
}

不管是使用DoubletFinder还是scrublet，我们最终都获得了一个分数和一个基于分数得到的判断结果。之后的步骤就是大家看一看doublet的分布，和在每一个细胞类型中的比例，来判断doublet的检测结果是否可靠，以及是要移除整一个细胞群还是移除被判定是doublet的细胞了。

祝大家吃好喝好睡好，科研快乐~

欢迎关注微信公众号 “小L的读博日常”，第一时间获得更多和生物信息学相关的小tips。

Ref：

[1] Wolock, Samuel L., Romain Lopez, and Allon M. Klein. "Scrublet: computational identification of cell doublets in single-cell transcriptomic data." Cell Systems 8.4 (2019): 281-291.