单细胞交响乐12-检测Doublet

刘小泽写于2020.7.15
为何取名叫“交响乐”？因为单细胞分析就像一个大乐团，需要各个流程的协同配合
单细胞交响乐1-常用的数据结构SingleCellExperiment
单细胞交响乐2-scRNAseq从实验到下游简介
 单细胞交响乐3-细胞质控
 单细胞交响乐4-归一化
 单细胞交响乐5-挑选高变化基因
 单细胞交响乐6-降维
 单细胞交响乐7-聚类分群
 单细胞交响乐8-marker基因检测
 单细胞交响乐9-细胞类型注释
 单细胞交响乐10-数据集整合后的批次矫正
 单细胞交响乐11-多样本间差异分析

1 前言

scRNA中，doublets指的就是一个文库中存在两个细胞的情况。一般是由于技术误差导致的（比如细胞分选、捕获），尤其在基于液滴的技术中（例如10X）比较突出。它干扰了对单个细胞表达量以及细胞形态的判断，比如一个液滴中有两个细胞，会通过误导我们这个细胞可能处于分化的“过渡态”。因此检测和去除这部分的影响至关重要。

一个比较通用的检测方法是：单纯根据表达量(Dahlin et al. 2018)。下面将利用一个10X数据来展示两种检测方法，它们的主要区别是我们是否需要提前知道分群的信息

数据准备

#--- 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)

参考帮助文档，它的具体做法是：
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 intermediate between the two sources after library size normalization.

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

返回的结果包括：

基因数量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排个序

library(scater)
chosen.doublet <- rownames(dbl.out)[isOutlier(dbl.out$N, 
    type="lower", log=TRUE)]
chosen.doublet
## [1] "6"

挑出来怀疑对象，可以对cluster6进一步检查

比如找到cluster6的marker基因

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）中也有类似的表达模式？

另外，基于背景知识，没有细胞会同时表达basal cells (Acta2) and alveolar cells (Csn2) 这两个基因，但是看到cluster6中这两个基因表达量都较高，再一次证明了我们的假设：cluster6是一个doublet混合体，而不是纯粹的一个细胞类型

plotExpression(sce.mam, features=c("Acta2", "Csn2"), 
    x="label", colour_by="label")

注意

从上面也能看到，doubletCluster()的优点就是方便操作，并且结果比较好理解。一旦有怀疑的cluster，就可以立即检查。但缺点是高度依赖分群的质量，如果分群不好，那么这个结果的可信度就会大打折扣。另外，如果真的有某个亚群中细胞数量很少，带来的结果就是：N小，让这个亚群很有可能成为怀疑对象。

不过，随着scRNA的技术改进，这个doublets情况会逐渐好转

2.2 基于模拟推断的检测

将利用来自scran包的doubletCells() ，它基于的假设是：模拟的doublets和真实的doublets接近

它的算法是：

随机选取两个原始单细胞表达量，加和，当做一个模拟的doublet，这样操作上千次
对每一个原始细胞，看看它附近有多少的模拟的doublet，并计算密度
对每一个原始细胞，同时计算它附近的其他原始细胞的数量，并计算密度
对每一个细胞，计算两个密度的比值，作为“doublet score”

为了加快密度的计算，这个函数会进行一个PCA以及log转换

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画出来，数值较高的细胞聚成一团

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

同时，结合之前doubletCluster()的结果看一眼：多么明显的cluster6！

注意

doubletCells() 的优点就是不需要依赖分群结果，降低了分群的影响。缺点是需要对doublets的模拟更精确，真的要保证它和真实的情况接近。

另外简单去掉那些doublet scores较高的细胞有时也是不够的，例如一个潜在的doublet cluster中只有一小部分的分值高，去掉它们剩下的细胞依然会干扰判断。那么问题来了，怎么样才叫高的doublet scores呢？是不是得设置一个阈值？就像刚才这种情况，如果把阈值降低会不会就多包括了一些潜在的doublets呢？其实这也是生信中的一个比较头疼的问题，没有一个固定的阈值或者范围，一切都是相对的。

比较推荐的方法是将doubletCells() 的结果再用分群展示出来，就像上面的图，可以更直观看到那些细胞影响较大。

小结

检测doublet操作必须要求数据来自同一批次，因为doublet也不会来自两次捕获的细胞，它毕竟是一个文库，只是包含两个细胞，因此也不需要担心检测doublet时怎么去除批次效应之类的问题。最好还是在分析前最好先清楚实验设计
如果数据中包含了细胞分化轨迹的信息，doublet就不好判断了，因为处于变化状态的细胞就很像doublet，容易被误认

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!