# 3.10 检测Doublet ## 1 前言 scRNA中，**doublets指的就是一个文库中存在两个细胞的情况**。一般是由于技术误差导致的（比如细胞分选、捕获），尤其在包含几千个细胞的基于液滴的技术中比较突出 ([Zheng et al. 2017](https://pubmed.ncbi.nlm.nih.gov/28091601/))。它干扰了对单个细胞表达量以及细胞形态的判断，比如一个液滴中有两个细胞，会通过误导我们这个细胞可能处于分化的“过渡态”。因此检测和去除这部分的影响至关重要。 一个比较通用的检测方法是：单纯根据表达量([Dahlin et al. 2018](https://pubmed.ncbi.nlm.nih.gov/29588278/))。下面将利用一个10X数据来展示两种检测方法，它们的主要区别是我们是否需要提前知道分群的信息 ### **数据准备** ```r #--- loading ---# library(scRNAseq) sce.mam <- BachMammaryData(samples="G_1") #--- gene-annotation ---# library(scater) rownames(sce.mam) <- uniquifyFeatureNames( rowData(sce.mam)$Ensembl, rowData(sce.mam)$Symbol) library(AnnotationHub) ens.mm.v97 <- AnnotationHub()[["AH73905"]] rowData(sce.mam)$SEQNAME <- mapIds(ens.mm.v97, keys=rowData(sce.mam)$Ensembl, keytype="GENEID", column="SEQNAME") #--- quality-control ---# is.mito <- rowData(sce.mam)$SEQNAME == "MT" stats <- perCellQCMetrics(sce.mam, subsets=list(Mito=which(is.mito))) qc <- quickPerCellQC(stats, percent_subsets="subsets_Mito_percent") sce.mam <- sce.mam[,!qc$discard] #--- normalization ---# library(scran) set.seed(101000110) clusters <- quickCluster(sce.mam) sce.mam <- computeSumFactors(sce.mam, clusters=clusters) sce.mam <- logNormCounts(sce.mam) #--- variance-modelling ---# set.seed(00010101) dec.mam <- modelGeneVarByPoisson(sce.mam) top.mam <- getTopHVGs(dec.mam, prop=0.1) #--- dimensionality-reduction ---# library(BiocSingular) set.seed(101010011) sce.mam <- denoisePCA(sce.mam, technical=dec.mam, subset.row=top.mam) sce.mam <- runTSNE(sce.mam, dimred="PCA") #--- clustering ---# snn.gr <- buildSNNGraph(sce.mam, use.dimred="PCA", k=25) colLabels(sce.mam) <- factor(igraph::cluster_walktrap(snn.gr)$membership) sce.mam ## class: SingleCellExperiment ## dim: 27998 2772 ## metadata(0): ## assays(2): counts logcounts ## rownames(27998): Xkr4 Gm1992 ... Vmn2r122 CAAA01147332.1 ## rowData names(3): Ensembl Symbol SEQNAME ## colnames: NULL ## colData names(5): Barcode Sample Condition sizeFactor label ## reducedDimNames(2): PCA TSNE ## altExpNames(0): ``` 或者[直接加载分享的数据](https://share.weiyun.com/ReZZnAMw) ## 2 两种检测方法 ### **2.1 基于分群结果的检测** 使用`doubletCluster()` 将任一cluster与其他另外两个clusters的表达量进行比较，即3个cluser为一组，其中一个是query，另外两个是source。基于的原假设是：query cluster的细胞中如果包含doublet，那么它是来自两个source clusters。 ([Bach et al. 2017](https://pubmed.ncbi.nlm.nih.gov/29225342/)) > 参考帮助文档，它的具体做法是： For each “query” cluster, we examine **all possible pairs of “source” clusters,** hypothesizing that the **query consists of doublets formed from the two sources**. > > If so, gene expression in the query cluster should be strictly i**ntermediate between the two sources** after library size normalization. ```r library(scran) # sce.mam一共分了10群，所以下面的结果也是10行，每一群都做了一次检测 > table(sce.mam$label) 1 2 3 4 5 6 7 8 9 10 550 799 716 452 24 84 52 39 32 24 dbl.out <- doubletCluster(sce.mam) > dim(dbl.out) [1] 10 9 ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-15-080546.png) 返回的结果包括： * 基因数量N：query cluster与另外两个source cluster相比特有的差异基因，它的数量多少可以为拒绝原假设提供依据。这个基因数量越少，query越有可能是doublets * 文库大小比例 lib.size：每个source cluster的各个细胞文库大小中位数 / query cluster中的各个细胞文库大小中位数，因此两个source就对应两个lib.size。我们知道doublets是一个文库包含两个细胞，因此它会比单个细胞的文库更大。这个值越小，query越可能是doublets * 占全部细胞的百分比 prop：query中的细胞数量占全部细胞的百分比，一般这个值小于5%，不过也取决于10X机器上样的数量 最后主要还是看N的数量，可以将这个结果按照N排个序 ```r library(scater) chosen.doublet <- rownames(dbl.out)[isOutlier(dbl.out$N, type="lower", log=TRUE)] chosen.doublet ## [1] "6" ``` 挑出来怀疑对象，可以对cluster6进一步检查 比如找到cluster6的marker基因 ```r markers <- findMarkers(sce.mam, direction="up") dbl.markers <- markers[[chosen.doublet]] > dim(dbl.markers) [1] 27998 13 # 然后找Top10基因 library(scater) chosen <- rownames(dbl.markers)[dbl.markers$Top <= 10] > length(chosen) [1] 43 # 最后热图 plotHeatmap(sce.mam, order_columns_by="label", features=chosen, center=TRUE, symmetric=TRUE, zlim=c(-5, 5)) ``` 看到这个cluster的marker基因，是不是在它的source cluster（即cluster1、2）中也有类似的表达模式？ ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-15-083507.png) 另外，基于背景知识，没有细胞会同时表达basal cells (**Acta2**) and alveolar cells (**Csn2**) 这两个基因，但是看到cluster6中这两个基因表达量都较高，再一次证明了我们的假设：cluster6是一个doublet混合体，而不是纯粹的一个细胞类型 ```r plotExpression(sce.mam, features=c("Acta2", "Csn2"), x="label", colour_by="label") ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-15-083828.png) **注意** 从上面也能看到，`doubletCluster()`的优点就是方便操作，并且结果比较好理解。一旦有怀疑的cluster，就可以立即检查。但缺点是高度依赖分群的质量，如果分群不好，那么这个结果的可信度就会大打折扣。另外，如果真的有某个亚群中细胞数量很少，带来的结果就是：N小，让这个亚群很有可能成为怀疑对象。 不过，随着scRNA的技术改进，这个doublets情况会逐渐好转 ### **2.2 基于模拟推断的检测** 将利用来自scran包的`doubletCells()` ，它基于的假设是：模拟的doublets和真实的doublets接近 它的算法是： * 随机选取两个原始单细胞表达量，加和，当做一个模拟的doublet，这样操作上千次 * 对每一个原始细胞，看看它附近有多少的模拟的doublet，并计算密度 * 对每一个原始细胞，同时计算它附近的其他原始细胞的数量，并计算密度 * 对每一个细胞，计算两个密度的比值，作为“doublet score” 为了加快密度的计算，这个函数会进行一个PCA以及log转换 ```r library(BiocSingular) set.seed(100) dbl.dens <- doubletCells(sce.mam, subset.row=top.mam, d=ncol(reducedDim(sce.mam))) summary(dbl.dens) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.00 7.63 21.04 395.42 49.30 39572.52 ``` 然后把计算结果的doublet score画出来，数值较高的细胞聚成一团 ```r sce.mam$DoubletScore <- log10(dbl.dens+1) plotTSNE(sce.mam, colour_by="DoubletScore") ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-15-085115.png) 同时，结合之前`doubletCluster()`的结果看一眼：**多么明显的cluster6！** ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-15-085508.png) **注意** `doubletCells()` 的优点就是不需要依赖分群结果，降低了分群的影响。缺点是需要对doublets的模拟更精确，真的要保证它和真实的情况接近。 另外简单去掉那些doublet scores较高的细胞有时也是不够的，例如一个潜在的doublet cluster中只有一小部分的分值高，去掉它们剩下的细胞依然会干扰判断。那么问题来了，**怎么样才叫高的doublet scores呢？是不是得设置一个阈值？**就像刚才这种情况，如果把阈值降低会不会就多包括了一些潜在的doublets呢？其实这也是生信中的一个比较头疼的问题，没有一个固定的阈值或者范围，一切都是相对的。 比较推荐的方法是将`doubletCells()` 的结果再用分群展示出来，就像上面的图，可以更直观看到那些细胞影响较大。 #### 小结 * 检测doublet操作必须要求数据来自同一批次，因为doublet也不会来自两次捕获的细胞，它毕竟是一个文库，只是包含两个细胞，因此也不需要担心检测doublet时怎么去除批次效应之类的问题。最好还是在分析前最好先清楚实验设计 * 如果数据中包含了细胞分化轨迹的信息，doublet就不好判断了，因为处于变化状态的细胞就很像doublet，容易被误认