欢迎关注同名公主号：BBio
10X也是基于此文献的，考古一下，学习区分空载和真实细胞大体思路。
当时已经存在的一些方法会假设含有细胞的GEMs会有更高的UMI总数，并以UMI数目指标筛选细胞，但是这种方法难以区分本就存在的小细胞和空载。
文章开发了一种新方法，首先评估ambiant RNA的表达特征，然后检验每个barcode和ambiant RNA的差异，有显著差异的barcode认为是一个真实的细胞，并结合barcode rank曲线的拐点，保证总UMI数较多的barcode始终保留。
//文献
EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data
//原理
- 构建ambient RNA池的表达谱
首先确定一个UMI阈值，默认为100，UMI数目低于阈值的定义为ambient RNA。每个droplet中相同基因的UMI数目总和为ambient RNA表达谱中该基因的UMI数目，得到所有基因的UMI数目。使用Good-Turing算法处理，生成每个基因的UMI数目比例的期望值。
假设溶液中的转录本随机的封装到空载中，对于每个droplet来说，每个基因的转录本被抽到的概率和期望值相同。使用Monte Carlo计算每个barcode的p-value。
- 检测barcode rank曲线的拐点
使用p-value可以筛选和ambient RNA有显著差异的barcode，但是有些情况下还能存在问题。ambient RNA是有很多破裂的细胞组成的，很难代表任何单一的细胞，但是当细胞群高度均匀，或者一个更易裂解的细胞亚群不成比例地贡献ambient RNA时，就可能存在barcode和ambient RNA表达相似。barcode序列测序错误也可能会对ambient RNA的估计产生偏差，原因是将包含细胞的droplet的UMI数错误地分配给具有低UMI总数的barcode。
通过绘制barcode rank plot，并计算曲线的拐点，第一个拐点标志UMI总数从高到低的快速转变。UMI总数较高的barcode都应该认定为是一个真实细胞。以拐点的UMI数为阈值，凡是大于阈值的barcode都认为是一个真实细胞。拐点下方的细胞也能因为和ambient RNA的显著差异认定为细胞，这是其它方法做不到的。
//10X官网的cell calling方法描述
- It uses a cutoff based on total UMI counts of each barcode to identify cells. This step identifies the primary mode of high RNA content cells.
- Then the algorithm uses the RNA profile of each remaining barcode to determine if it is an “empty" or a cell containing partition. This second step captures low RNA content cells whose total UMI counts may be similar to empty GEMs.
- Therefore, starting from Cell Ranger 6.1, it is recommended to run all analyses with the --expect-cells option with a reasonable estimate of recovered cells, especially for higher cell load experiments.
首先cellranger软件expect-cells参数(默认3000)作为期望细胞数，对这些细胞的UMI总数进行排序，并以99%分位数除以10作为UMI阈值。第二步的描述也就和EmptyDrops方法相同了。
补充说明里对于expect-cells参数的选择应该有合理的评估，但评估出准确的阈值似乎很难。当把两个步骤鉴定的细胞数的并集作为最终的细胞时，依赖expect-cells的第一步也显得很重要。当细胞上样量过大或者过小时，需要谨慎选择expect-cells。
//R包DropletUtils测试
模拟10X数据，并绘制barcode rank plot。
set.seed(0)
my.counts <- DropletUtils:::simCounts()
br.out <- barcodeRanks(my.counts)
names(br.out)
# Making a plot.
plot(br.out$rank, br.out$total, log="xy", xlab="Rank", ylab="Total")
o <- order(br.out$rank)
lines(br.out$rank[o], br.out$fitted[o], col="red")
abline(h=metadata(br.out)$knee, col="dodgerblue", lty=2)
abline(h=metadata(br.out)$inflection, col="forestgreen", lty=2)
legend("bottomleft", lty=2, col=c("dodgerblue", "forestgreen"),
legend=c("knee", "inflection"))
测试emptyDrops函数。
out <- emptyDrops(my.counts)
out
#DataFrame with 11100 rows and 5 columns
# Total LogProb PValue Limited FDR
# <integer> <numeric> <numeric> <logical> <numeric>
#1 2 NA NA NA NA
#2 9 NA NA NA NA
#3 20 NA NA NA NA
#4 20 NA NA NA NA
#5 1 NA NA NA NA
#... ... ... ... ... ...
#11096 215 -246.428 9.999e-05 TRUE 0.00013799
#11097 201 -250.234 9.999e-05 TRUE 0.00013799
#11098 247 -275.905 9.999e-05 TRUE 0.00013799
#11099 191 -228.763 9.999e-05 TRUE 0.00013799
#11100 198 -233.043 9.999e-05 TRUE 0.00013799
is.cell <- out$FDR <= 0.001
sum(is.cell, na.rm=TRUE)
#[1] 942
