预测单细胞表达谱中双细胞的软件有不少，上一篇用Scrublet这个工具预测了一下(https://www.jianshu.com/p/dc44038e9928)，这回试一下DoubletFinder，并对两者的预测结果做个对比。
DoubletFinder原理和Scrublet非常相似（下图），也是通过模拟出双细胞表达谱，对原始细胞和模拟细胞一起聚类，然后通过判断和原始细胞聚集在一起的模拟双细胞的多少对原始细胞打分。通常来说得分高的细胞就会被判定为双细胞。
Screen Shot 2021-01-20 at 5.22.29 PM.png
首先读入同一个PBMC数据并简单跑一下Seurat流程。
data<-Read10X(data.dir="./test_data/pbmc8k_10X/filtered_gene_bc_matrices/GRCh38/")
Seu_obj <- CreateSeuratObject(counts = data)
Seu_obj <- NormalizeData(Seu_obj)
Seu_obj <- FindVariableFeatures(Seu_obj, selection.method = "vst", nfeatures = 2000)
Seu_obj <- ScaleData(Seu_obj)
Seu_obj <- RunPCA(Seu_obj)
Seu_obj <- RunTSNE(Seu_obj, dims = 1:30)
Seu_obj <- FindNeighbors(Seu_obj, dims = 1:30)
Seu_obj <- FindClusters(Seu_obj, resolution = 0.6)
然后跑DoubletFinder流程。主要引入三个参数pN, pK 和 pANN。pN代表要模拟的双细胞数。 pK是对每个细胞选取多少与之邻近的细胞进行下一步计算。pANN则是在pK范围内有多少模拟的双细胞。
sweep.res.list_pbmc <- paramSweep_v3(Seu_obj, PCs = 1:30, sct = FALSE)
sweep.stats_pbmc <- summarizeSweep(sweep.res.list_pbmc, GT = FALSE)
bcmvn_pbmc <- find.pK(sweep.stats_pbmc)
opt_pK <- as.numeric(as.vector(bcmvn_pbmc$pK[which.max(bcmvn_pbmc$BCmetric)]))
print(opt_pK)
# [1] 0.09
annotations <- Seu_obj@meta.data$seurat_clusters
homotypic.prop <- modelHomotypic(annotations)
nExp_poi <- round(0.06*nrow(Seu_obj@meta.data))
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
Seu_obj <- doubletFinder_v3(Seu_obj,
PCs = 1:30,
pN = 0.25,
pK = opt_pK,
nExp = nExp_poi,
reuse.pANN = FALSE,
sct = FALSE)
比较重要的一步是找到合适的pK值。作者引入了一个参数BCmetric，是他们通过分析人和老鼠混合单细胞样本（知道真实双细胞率的样本）得出的经验参数，其最大值对应合适的pK值。如下图所示，该样本的pK值为0.09
image.png
最后我们读入上次用Scrublet对同一个数据的预测结果并进行比较。
scrublet_pred <- read.csv("./scrublet/doublet_prediction.txt", header=T, sep=",")
scrublet_pred[scrublet_pred$prediction=="False",]$prediction <- "Singlet"
scrublet_pred[scrublet_pred$prediction=="True",]$prediction <- "Doublet"
Seu_obj@meta.data$scrublet <- scrublet_pred$prediction
p0 <- DimPlot(Seu_obj, reduction = "tsne", label=T)
p1 <- DimPlot(Seu_obj, reduction = "tsne", group.by = "DF.classifications_0.25_0.22_503",order = T) + ggtitle("DoubletFinder")
p2 <- DimPlot(Seu_obj, reduction = "tsne", group.by = "scrublet", order = T) + ggtitle("Scrublet")
p0+p1+p2
可以看到两者结果还挺一致的。主要在cluster 0和6有双细胞。
image.png
最后附一张pANN在pN为0.25时选取不同pK值时的分布。有意思的是当pK很小的时候是一个均值较小的单峰，随着pK增大开始出现双峰甚至三个峰。其中pN_0.25_pK_0.09就是最终挑选出的参数值，可以看到此时还没有出现双峰， 这点让人比较费解。
所以这两个软件相比较Scrublet更为直观一些～
pANN.partial.png
Reference:
https://github.com/chris-mcginnis-ucsf/DoubletFinder
网友评论