OSCA单细胞数据分析笔记14—Empty/Doublet droplet

对应原版教程第15、16章：http://bioconductor.org/books/release/OSCA/overview.html

现行主流的Droplet-based单细胞测序技术主要思路是一个磁珠捕获一个细胞置于油包水的腔室里完成添加标签、建库操作。但在磁珠捕获的过程会出现未捕获到细胞或者两个细胞的异常情况。这就需要我们在分析单细胞数据中识别、过滤掉这些bad barcode(cell)。

1、Empty droplet

1.1 关于Empty

• 空液滴，即磁珠未捕获到细胞，是Droplet-based单细胞测序技术的常见现象。一般Cellranger分析得到的单细胞表达矩阵往往是过滤了Empty droplet之后的，一般只有数千至上万个细胞。如果下载的表达矩阵里有几十万的细胞，那么很有可能是未过滤empty droplet的原始表达矩阵。例如GSE138665数据集。
• 空液滴并非没有任何基因表达数据。由于溶液环境中的ambient RNA造成空液滴实际上会有少数基因的结果(lower library)。所以最简单的方式就是直接设定一个阈值过滤lower library的droplet/barcode/cell。但这就可能误删本身基因表达量就少的细胞类型。
• 如下实例数据集是一个未过滤emoty droplet的原始测序数据

#--- loading ---#
library(DropletTestFiles)
raw.path <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/raw.tar.gz")
out.path <- file.path(tempdir(), "pbmc4k")
untar(raw.path, exdir=out.path)
library(DropletUtils)
fname <- file.path(out.path, "raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
sce.pbmc
#class: SingleCellExperiment 
#dim: 33694 737280

我们可以根据每个细胞的counts数大小从高到底排序，观察是否有明显的断崖趋势。

library(DropletUtils)
bcrank <- barcodeRanks(counts(sce.pbmc))
# Only showing unique points for plotting speed.
uniq <- !duplicated(bcrank$rank)
plot(bcrank$rank[uniq], bcrank$total[uniq], log="xy",
     xlab="Rank", ylab="Total UMI count", cex.lab=1.2)
abline(h=metadata(bcrank)$inflection, col="darkgreen", lty=2)
abline(h=metadata(bcrank)$knee, col="dodgerblue", lty=2)
legend("bottomleft", legend=c("Inflection", "Knee"), 
       col=c("darkgreen", "dodgerblue"), lty=2, cex=1.2)

如下图呈现结果，在count数1000左右，细胞数量突然减少可能就意味着empty droplet。

1.2 去除empty droplet

• 如上所说，直接设定一个low count阈值去除empty droplet有些太过武断。针对上述例子：主要在于判断counts数处于500~1000的droplet是受ambient RNA影响的empty droplet还是自身表达量就很低的细胞类型 。
• 这里我们使用DropletUtils包的emptyDrops()函数尽可能区分这两种混淆的droplet。
• 其主要思路：首先设定counts量明显过低(<100)一定是empty droplet，然后(空假设)根据计算其余droplet的表达模式是否明显偏离empty droplet的表达模式。最后我们保留p值显著(<0.05)的结果，认为是捕获到cell的droplet，用于后续的数据分析。

set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
summary(e.out$FDR <= 0.001)
##    Mode   FALSE    TRUE    NA's 
## logical     989    4300  731991
#最后保留得到4300个细胞
#NA值表示count < 100的empty droplet(default parameter，可自定义修改)

• 观察去除empty droplet 之后的counts分布情况

retained <- sce.pbmc[,which(e.out$FDR <= 0.001)]
bcrank <- barcodeRanks(counts(retained))
uniq <- !duplicated(bcrank$rank)
plot(bcrank$rank[uniq], bcrank$total[uniq], log="xy",
     xlab="Rank", ylab="Total UMI count", cex.lab=1.2)

• 可以看到依旧存在明显的拐点，但转点阈值已经提高到2000左右，而且只有不到500个细胞，可能是测序不合格的细胞，待下一步的QC解决。

• 最后值得一提的是emptyDrops()进行假设检验的方法是蒙特卡洛抽样法(Monte Carlo simulations)。一般来说抽样次数越多，p值越低(前提是该droplet的确有可能不符合空假设)。所有对于emptyDrops()返回结果的limited为logical，表示是否可以随着增加抽样次数，进一步降低p值，从而挽救可能是cell captured droplet，但却被判定为empty droplet。但对于上述的分析结果，如下所示无该风险，可放心进行过滤。

table(Sig=e.out$FDR <= 0.001, Limited=e.out$Limited)
##        Limited
## Sig     FALSE TRUE
##   FALSE   989    0
##   TRUE   1728 2572

2、Doublet droplet

2.1 关于Doublet droplet

• Doublet droplet，即一个磁珠同时捕获到了两个细胞，很有可能同时含有两种细胞类型的表达特征，会对后续的分析造成很大的混淆。
• 目前cellranger软件尚不能识别、过滤Doublet droplet。下面介绍了基于表达矩阵的，两种判断是否为Doublet方法。
• 如下为示例数据

sce.mam
#class: SingleCellExperiment 
#dim: 27998 2772

2.2 方法一：Detection with clusters

• scDblFinder包的findDoubletClusters()函数。这种方法根据聚类结果，判断出是否由Doublet聚成的1个 cluster。
• 该方法的空假设为：一个cluster(potential Doublet)的cell 是由两个 origin cluster里的cell混合而成的。所以该方法对于一个cluster考虑其余所有cluster的两两组合视为potential origin cluster。
• 该结果为每个cluster找到了最有可能的2个 origin cluster，供分析者判断。

library(scDblFinder)
dbl.out <- findDoubletClusters(sce.mam)
dbl.out

如下结果含义可通过?findDoubletClusters查看。关键的几个指标有 num.de: query cluster与 source cluster的差异基因数(越少越有可能为doublet)，lib.size: source cluster的median cell library/query cluster的median cell library(越少越有可能为doublet，因为doublet cell的library size会大一些)。此外prop表示query cluster的细胞数占总细胞数的比例，这可为我们对预期的doublet 数量提供判断(一般小于5%)。

• 根据上述结果，进一步判断最可能的doublet cluster

library(scater)
chosen.doublet <- rownames(dbl.out)[isOutlier(dbl.out$num.de, type="lower", log=TRUE)]
chosen.doublet
dbl.out[chosen.doublet,]
#DataFrame with 1 row and 9 columns
#      source1     source2    num.de median.de        best    p.value lib.size1 lib.size2      prop
#  <character> <character> <integer> <numeric> <character>  <numeric> <numeric> <numeric> <numeric>
#6           2           1        13     507.5       Pcbp2 0.00128336  0.811531  0.516399  0.030303

• 根据来自两种细胞的marker gene进一步验证该cluster是否由doublet聚成的

#basal cells (Acta2) and alveolar cells (Csn2)
plotExpression(sce.mam, features=c("Acta2", "Csn2"), 
    x="label", colour_by="label")

如下图，的确cluster6高表达两种细胞类型的marker gene，很有可能就是doublet cluster。在后续分析中予以去除。

最后关于这种基于cluster的判断doublet方法比较依赖合适的聚类结果。如教程所说：Clusters that are too coarse will fail to separate doublets from other cells, while clusters that are too fine will complicate interpretation.

2.2 方法二：Detection by simulation

• scran包的computeDoubletDensity()函数。这种方法会为每一个细胞计算doublet score，越高表明越可能是doublet。
• 其主要思路是：首先从数据集中通过合并随机两个细胞的表达谱，模拟上千个doublet添加到表达矩阵中。然后根据PCA降维结果，分别以每一个细胞为中心，计算该细胞周围(高维距离)模拟doublet的密度与真实droplet的密度之间的比例。score比值越高表明该细胞越接近simulated doublet。

library(BiocSingular)
set.seed(100)
# Setting up the parameters for consistency with denoisePCA();
# this can be changed depending on your feature selection scheme.
dbl.dens <- computeDoubletDensity(sce.mam, subset.row=top.mam, 
    d=ncol(reducedDim(sce.mam)))
summary(dbl.dens)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.249   0.527   1.041   1.188  14.570

• 可视化 DoubletScore

sce.mam$DoubletScore <- dbl.dens
plotTSNE(sce.mam, colour_by="DoubletScore")

plotColData(sce.mam, x="label", y="DoubletScore", colour_by="label")

• 如上看到还是主要集中第六群最值得怀疑。直接删去high score的suspected doublet缺乏可解释性，因为很难定义一个合理的阈值(predicted doublet proportion？)。这里还是建议结合聚类结果进行判断出suspected cluster。

最后强调一点：无论是empty droplet还是doublet，都应当基于单批次的样本进行操作处理。因为每个批次的测序环境都是不一致的，不能够符合统一的假设。即对于多批次的样本，应该进行分别的鉴别、过滤，再进行后续的分析流程。