Doublets检测1----DoubletFinder

目前比较常用的单细胞测序技术，比如油包滴技术，是通过形成一个个液滴，让一个液滴内包含一个细胞，对细胞进行裂解，释放mRNA，添加barcode再逆转录成cDNA最后进行批量测序的一个过程。因为我们假定一个液滴内仅包含一个细胞，因而拥有相同barcode、来自同一个液滴的mRNA，我们在分析中也会当成同一个细胞的数据。但在实际情况中，有一些液滴会包含超过两个或两个以上细胞，被称之为doublet或multiplet。因为同一个液滴内的mRNA会拥有一样的barcode，当液滴内包含超过一个细胞时，它们依旧会被当做一个细胞处理，在后续分析中，引入一些错误的信息。

我们在进行分析前通常会使用细胞表达的总reads数、细胞表达的基因数以及reads中属于线粒体基因的比例这类信息去对细胞做一个基本的筛选，去除掉一些低质量的数据。在这一过程中，我们能去除掉部分拥有高于正常值的总reads数和基因数的doublet。

但通过这样简单的筛选，有的时候并没有办法完全去除doublet和multiplet。一些分化程度高、只拥有有限、固定功能的细胞，可能会有比较低的总reads数及基因数。它们只拥有比较固定的功能，因此不需要表达大量不同的基因。当样本中包含这类细胞时，他们与其它细胞在一个液滴中形成的doublet就不会有太高的总reads数及基因数。在这种情况下，我们需要运用一些其它的方法检测出可能的doublet。其中DoubletFinder和scrublet是各个帖子比较推荐的两种。

今天我们主要来学习DoubletFinder。

官网地址：https://github.com/chris-mcginnis-ucsf/DoubletFinder

DoubletFinder的原理是从现有的矩阵的细胞中根据我们预先定义好的细胞类型模拟一些双细胞出来（比如单核和T细胞的双细胞、B细胞和中性粒细胞的双细胞等等），将模拟出的双细胞和原有矩阵的细胞混合在一起，进行降维聚类，原则上人工模拟的doublets会与真实的doublets距离较近，因此计算每个细胞K最近邻细胞中人工模拟doublets的比例 (pANN)，就可以根据pANN值对每个barcode的doublets概率进行排序。另外依据泊松分布的统计原理可以计算每个样本中doublets的数量，结合之前的细胞pANN值排序，就可以过滤doublets了。

所以，整体来说，DoubletFinder主要分为以上这几个步骤：

1. 首先，我们有一个未分析的scRNA数据。

2. 因为doublet是两个细胞被当成了一个细胞处理，DoubletFinder通过随机将两个细胞的基因表达数据相加来模拟可能出现的doublet的基因表达情况，并把这些模拟出来的doublet也当作正常的细胞进行后续的处理。

3. 通过常规的PCA对reads数矩阵（cell×gene matrix）进行降维，更好地对信息进行压缩和提取。

4. 找出每一个细胞的邻近细胞（neighbors）

5. 在这些细胞中，包含了我们数据中原有的细胞，也包括模拟出来的doublet。DoubletFinder会计算出每个细胞的邻近细胞中包含多大比例模拟出来的doublet。比例高即说明这个细胞和模拟出来的doublet相似，即更有可能是doublet。DoubletFinder用一个用户自己设定的值pANN作为分界线，高于设定的pANN即为doublet。

6. 去除DoubletFinder认定为doublet的细胞。

注：双细胞去除与Seurat的交互：在seurat标准流程进行到TSNE和UMAP降维之后，FindAllMarkers之前进行DoubletFinder操作。

========安装============

devtools::install_github('chris-mcginnis-ucsf/DoubletFinder')

======使用测试======

library(tidyverse)

library(Seurat)

library(DoubletFinder)

###我们还用经常用的pbmc数据来测试一下

pbmc <- readRDS("pbmc.rds")

###寻找最优pK值

sweep.res.list <- paramSweep_v3(pbmc, PCs = 1:10, sct = FALSE)

sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)

bcmvn_pbmc <- find.pK(sweep.stats)

###检测双细胞

通常情况下，双细胞有两种，同源双细胞和异源双细胞。DoubletFinder只能检测异源双细胞。所以需要把同源双细胞可能的比率去除掉，以优化期望的doublets数量。

DoubletRate= 0.075        #直接查表，10000细胞对应的doublets rate是～7.6%

DoubletRate= ncol(pbmc)*8*1e-6      #更通用的一种参数设置

#估计同源双细胞比例

homotypic.prop <- modelHomotypic(pbmc$seurat_clusters)

nExp_poi <- round(DoubletRate*ncol(pbmc))      #计算双细胞的比例

nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))    #校正

#开始鉴定双细胞

pbmc <- doubletFinder_v3(pbmc, PCs= 1:10, pN = 0.25, pK = pk_bcmvn, nExp = nExp_poi.adj, reuse.pANN = F, sct = F)

主要参数如下所述：

http://events.jianshu.io/p/6770c6a05287

从下面的结果，可以看出预测完之后，亏有pANN和DF.classifications两个量。

####下面，我们看下预测出来的双细胞在图中分分布。因为，我们这个是一个测试数据，不具备实际意义。

DimPlot(pbmc, reduction = "umap", group.by = "DF.classifications_0.25_0.02_46")

确定是否剔除预测双细胞的一般步骤：

1. 查看降维图上预测出来的双细胞位置。在降维图上，真实的双细胞更倾向于处于主群边缘或游离出来。比如，我们用的这个测试数据。

2. 如果预测的双细胞单独聚成群，查看其marker基因是否具有明显的多种细胞marker。

3. 不是我们主要研究对象的细胞群中的双细胞可以不去除。