scRNA-seq数据处理—demultiplexing

2020-03-27 13:22:27 浏览数 (2)

正文

处理原始scRNA-seq数据

3.4 demultiplexing

demultiplexing，指的就是将multiplexed的reads根据index从不同或者同一个lane中分出，生成sample对应的fastq文件

根据使用的protocol和完整pipeline中的的特定pipeline，demultiplexing的方式不同。我们所知道的最灵活的demultiplexing pipeline是zUMI，它可用于demultiplexing和大多数基于UMI的protocol的比对。对于Smartseq2或其他pair-end全转录的protocol，数据通常已经被分解。诸如GEO或ArrayExpress之类的公共存储库需要在上传之前对基于小规模/基于plate的scRNASeq数据进行分解，并且许多测序设备将在将数据返回给您之前自动demultiplexing。如果您没有使用已发布的pipeline，并且数据未被测序工具demultiplexing，则您必须自己做。这通常需要编写自定义脚本，因为barcode可能具有不同的长度和在read中有不同位置，具体取决于所使用的protocol。

对于所有数据类型“demultiplexing”涉及从一个或两个reads中识别和移除细胞barcode序列。如果提前知道预期的细胞barcode，即数据来自基于PCR-plate，那么所需要做的是将每个细胞barcode与预期的carcode进行比较，并让相关的reads和最相近的细胞-barcode（最大不匹配为1或2，取决于细胞-barcode的设计）进行比对。这些数据通常在比对之前被demultiplexing，作为平行比对步骤的简单方法。

我们有公开可用的 perl脚本，能够使用单个细胞barcode demultiplexing任何scRNASeq数据，无论它有或没有UMI用于plate-based的protocol。这些可以这样使用：

代码语言：javascript复制
perl 1_Flexible_UMI_Demultiplexing.pl 10cells_read1.fq 10cells_read2.fq "C12U8" 10cells_barcodes.txt 2 Ex
代码语言：javascript复制
## 
##  Doesn't match any cell: 0
##  Ambiguous: 0
##  Exact Matches: 400
##  Contain mismatches: 0
##  Input Reads: 400
##  Output Reads: 400
## Barcode Structure: 12 bp CellID followed by 8 bp UMI
代码语言：javascript复制
perl 1_Flexible_FullTranscript_Demultiplexing.pl 10cells_read1.fq 10cells_read2.fq "start" 12 10cells_barcodes.txt 2 Ex
代码语言：javascript复制
## 
## Doesn't match any cell: 0
## Ambiguous: 0
## Exact Matches: 400
## Contain Mismatches: 0
## Input Reads: 400
## Output Reads: 400

对于包含UMI的数据，demultiplexing包括将UMI代码附加到包含read的gene-body的read名称。如果数据来自基于droplet的protocol或SeqWell，其中预期barcode的数量远远高于预期的数量，那么通常细胞-barcode也将附加到read名称以避免产生大量的数据文件。在这些情况下，在量化步骤期间将发生demultiplexing，以便于识别对应于完整细胞而不是背景噪声的细胞-barcode。

3.4.1 识别含有细胞的droplet/microwells

对于基于液滴（droplet）的方法，仅一小部分液滴包含珠子（beads）和完整细胞。然而，生物学实验是混乱的，一些RNA会从死亡/受损细胞中泄漏出来。因此，没有完整细胞的液滴可能捕获少量环境RNA，这些RNA将最终进入测序文库，并为最终的测序输出提供read。液滴的大小，扩增效率和测序的变化将导致“背景”和真实细胞具有多种文库大小。已经使用各种方法来试图区分那些对应于真实细胞的那些细胞-barcode。

大多数方法使用每个barcode的所有分子（可以应用于总reads）并尝试通常是细胞加上一些背景的大文库和被认为是纯背景的小文库之间找到“断点”。加载一些包含大小细胞的示例模拟数据：

代码语言：javascript复制
umi_per_barcode <- read.table("droplet_id_example_per_barcode.txt.gz")
truth <- read.delim("droplet_id_example_truth.gz", sep=",")
练习

检测到多少个独特的barcode？数据中有多少真细胞？为简化此部分的计算，排除所有分子总数少于10的barcode。

答案

一种方法是寻找每个条形码的总分子突然下降的拐点：

代码语言：javascript复制
barcode_rank <- rank(-umi_per_barcode[,2])
plot(barcode_rank, umi_per_barcode[,2], xlim=c(1,8000))

在这里，我们可以看到库大小的粗略指数曲线，因此为了简单起见，我们可以对它们进行对数转换。

代码语言：javascript复制
log_lib_size <- log10(umi_per_barcode[,2])
plot(barcode_rank, log_lib_size, xlim=c(1,8000))

分布中的“膝盖”更加明显。我们可以手动估计“膝盖”的位置，但通过算法识别这个点。

代码语言：javascript复制
# inflection point
o <- order(barcode_rank)
log_lib_size <- log_lib_size[o]
barcode_rank <- barcode_rank[o]

rawdiff <- diff(log_lib_size)/diff(barcode_rank)
inflection <- which(rawdiff == min(rawdiff[100:length(rawdiff)], na.rm=TRUE))

plot(barcode_rank, log_lib_size, xlim=c(1,8000))
abline(v=inflection, col="red", lwd=2)
代码语言：javascript复制
threshold <- 10^log_lib_size[inflection]

cells <- umi_per_barcode[umi_per_barcode[,2] > threshold,1]
TPR <- sum(cells %in% truth[,1])/length(cells)
Recall <- sum(cells %in% truth[,1])/length(truth[,1])
c(TPR, Recall)
代码语言：javascript复制
## [1] 1.0000000 0.7831707

另一种方法是建造一个混合模型，找出较高和较低分布相交的位置。但是，数据可能不是很适合假定的分布:

代码语言：javascript复制
set.seed(-92497)
# mixture model
require("mixtools")
代码语言：javascript复制
## Loading required package: mixtools
代码语言：javascript复制
## mixtools package, version 1.1.0, Released 2017-03-10
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.
代码语言：javascript复制
mix <- normalmixEM(log_lib_size)
代码语言：javascript复制
## number of iterations= 43
代码语言：javascript复制
plot(mix, which=2, xlab2="log(mol per cell)")
代码语言：javascript复制
p1 <- dnorm(log_lib_size, mean=mix$mu[1], sd=mix$sigma[1])
p2 <- dnorm(log_lib_size, mean=mix$mu[2], sd=mix$sigma[2])
if (mix$mu[1] < mix$mu[2]) {
    split <- min(log_lib_size[p2 > p1])
} else {
    split <- min(log_lib_size[p1 > p2])
}
练习

使用此分割点识别cell并计算TPR和Recall

答案

CellRanger使用的第三种方法是假定真实细胞的文库大小范围为~10-fold，并使用预期的细胞数估计该范围。

代码语言：javascript复制
n_cells <- length(truth[,1])
# CellRanger
totals <- umi_per_barcode[,2]
totals <- sort(totals, decreasing = TRUE)
# 99th percentile of top n_cells divided by 10
thresh = totals[round(0.01*n_cells)]/10
plot(totals, xlim=c(1,8000))
abline(h=thresh, col="red", lwd=2)
练习

使用此阈值识别cell并计算TPR和Recall。

答案

最后（EmptyDrops）[ https://github.com/MarioniLab/DropletUtils ]，目前正在进行beta测试，它使用所有液滴的完整基因的细胞分子计数矩阵，并估计来自那些分布计数极低的液滴的“背景”RNA的谱，然后寻找具有与背景不同的基因表达谱的细胞。这与拐点方法相结合，因为背景RNA通常看起来非常类似于群体中较大细胞的表达谱。因此，EmptyDrops是唯一能够识别高度多样化样本中非常小的细胞barcode的方法。

下面我们提供了当前如何运行此方法的代码:(我们将在方法正式发布时更新此页面）

代码语言：javascript复制
require("Matrix")
raw.counts <- readRDS("droplet_id_example.rds")

require("DropletUtils")
# emptyDrops
set.seed(100)
e.out <- emptyDrops(my.counts)
is.cell <- e.out$FDR <= 0.01
sum(is.cell, na.rm=TRUE)
plot(e.out$Total, -e.out$LogProb, col=ifelse(is.cell, "red", "black"),
    xlab="Total UMI count", ylab="-Log Probability")

cells <- colnames(raw.counts)[is.cell]

TPR <- sum(cells %in% truth[,1])/length(cells)
Recall <- sum(cells %in% truth[,1])/length(truth[,1])
c(TPR, Recall)
barcode cell na pipeline umi

0 人点赞