Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（一）
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（二）
Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（三）
收藏|北大生信平臺(tái)"單細(xì)胞分析蓖谢、染色質(zhì)分析"視頻和PPT分享
生信老司機(jī)以中心法則為主線講解組學(xué)技術(shù)的應(yīng)用和生信分析心得 - 限時(shí)免費(fèi)
測(cè)序文庫(kù)拆分 (Demultiplexing)
文庫(kù)拆分因使用的前期Protocol不同或構(gòu)建的流程不同需要有對(duì)應(yīng)的處理方式捂蕴。我們認(rèn)為最靈活可用的文庫(kù)拆分工具是zUMIs (https://github.com/sdparekh/zUMIs/wiki/Usage)，可以用來(lái)拆分和比對(duì)大部分基于UMI的建庫(kù)方式闪幽。對(duì)于Smartseq2或其他雙端全長(zhǎng)轉(zhuǎn)錄本方案啥辨，數(shù)據(jù)通常已經(jīng)拆分好了。例如GEO或ArrayExpress之類的公共數(shù)據(jù)存儲(chǔ)庫(kù)會(huì)要求小規(guī)亩㈦纾或plate-based scRNASeq數(shù)據(jù)拆分好再上傳溉知，并且很多測(cè)序服務(wù)商提供的數(shù)據(jù)都是自動(dòng)拆分好的。如果使用的分析流程依賴于拆分好的數(shù)據(jù)但測(cè)序服務(wù)商提供的數(shù)據(jù)沒有拆分時(shí)就需要自己拆分腕够。因?yàn)椴煌慕◣?kù)方案引入的barcode序列的長(zhǎng)度和位置不同级乍，通常都需要自己寫腳本解決。
對(duì)于所有數(shù)據(jù)類型帚湘，”文庫(kù)拆分”涉及從一端或雙端短序列中識(shí)別和移除細(xì)胞條形碼序列 (cell barcode)玫荣。如果提前知道加入的細(xì)胞條形碼，比如數(shù)據(jù)來(lái)自基于PCR板的方案大诸，只需要找到條形碼并與條形碼庫(kù)作比對(duì)捅厂，歸類于與之最相似的那個(gè)就可以 (根據(jù)條形碼的設(shè)計(jì)，一般允許最多1-2個(gè)錯(cuò)配)底挫。這些數(shù)據(jù)通常在比對(duì)之前先做拆分恒傻，從而可以并行比對(duì)，提高效率建邓。
我們有公開可用 (<>)的 perl腳本盈厘，可以拆分任何plate-based的建庫(kù)方案生成的數(shù)據(jù)，不管有沒有UMI官边。操作如下：
# https://github.com/tallulandrews/scRNASeqPipeline
# perl 1_Flexible_UMI_Demultiplexing.pl read1.fq read2.fq b_structure index mismatch prefix\n";
# print STDERR "
# read1.fq : barcode/umi containing read
# read2.fq : non-barcode containing read
# b_structure : a single string of the format C##U# or U#C##
# where C## is the cell-barcode and U# is the UMI.
# e.g. C10U4 = a 10bp cell barcode followed by a 4bp UMI
# index : file containg a single column of expected cell-barcodes.
# if equal to \"UNKNOWN\" script will output read counts for each unique barcode.
# mismatch : maximum number of permitted mismatches (recommend 2)
# prefix : prefix for output fastq files.
perl 1_Flexible_UMI_Demultiplexing.pl 10cells_read1.fq 10cells_read2.fq "C12U8" 10cells_barcodes.txt 2 Ex
運(yùn)行過程輸出如下：
## 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
# https://github.com/tallulandrews/scRNASeqPipeline
# perl 1_Flexible_FullTranscript_Demultiplexing.pl read1.fq read2.fq b_pos b_len index mismatch prefix
# read1.fq : barcode containing read
# read2.fq : non-barcode containg read
# b_pos : position of cell-barcode in the read. [\"start\" or \"end\"]
# b_len : length of cell-barcode (bp)
# index : file contain a single column of expected barcodes
# mismatch : maximum number of permitted mismatches (recommend 2)
# prefix : prefix for output fq files.
perl 1_Flexible_FullTranscript_Demultiplexing.pl 10cells_read1.fq 10cells_read2.fq "start" 12 10cells_barcodes.txt 2 Ex
##
## Doesn't match any cell: 0
## Ambiguous: 0
## Exact Matches: 400
## Contain Mismatches: 0
## Input Reads: 400
## Output Reads: 400
對(duì)于包含UMI的數(shù)據(jù)沸手，文庫(kù)拆分包含把UMI code加到來(lái)源于基因區(qū)的序列的名字上。如果數(shù)據(jù)來(lái)源于droplet-based protocol或者SeqWell注簿，實(shí)際加入的barcode序列的種類多于捕獲到的細(xì)胞數(shù)時(shí)契吉，為了避免生成特別多的文件，一般也把cell-barcode加到測(cè)序reads的名字后面诡渴。在這種情況下捐晶，文庫(kù)拆分是在定量過程中進(jìn)行的，有利于識(shí)別來(lái)源于完整細(xì)胞而不是背景噪聲中的cell-barcode序列妄辩。
鑒定含有細(xì)胞的液滴/微孔
液滴型scRNA-seq方法中只有一小部分的液滴包含珠狀物和一個(gè)完整細(xì)胞惑灵。然而生物實(shí)驗(yàn)不會(huì)那么理想，有些RNA會(huì)從死細(xì)胞或破損細(xì)胞中漏出來(lái)眼耀。所以沒有完整細(xì)胞的液滴有可能捕獲周圍環(huán)境游離出的少了RNA并且走完測(cè)序環(huán)節(jié)出現(xiàn)在最終測(cè)序結(jié)果中英支。液滴大小、擴(kuò)增效率和測(cè)序環(huán)節(jié)中的波動(dòng)會(huì)導(dǎo)致”背景”和真實(shí)細(xì)胞最終獲得的文庫(kù)大小變化很大哮伟，使得區(qū)分哪些文庫(kù)來(lái)源于背景哪些來(lái)源于真實(shí)細(xì)胞變得復(fù)雜干花。
大多數(shù)方法使用每個(gè)barcode對(duì)應(yīng)的總分子數(shù)(如果是UMI)或總reads數(shù)的分布來(lái)尋找一個(gè)”break point”區(qū)分來(lái)自于真實(shí)細(xì)胞的較大的文庫(kù)和來(lái)自于背景的較小的文庫(kù)妄帘。下面加載一些包含大細(xì)胞和細(xì)胞的模擬數(shù)據(jù)實(shí)際操作下：
umi_per_barcode <- read.table("droplet_id_example_per_barcode.txt.gz")
truth <- read.delim("droplet_id_example_truth.gz", sep=",")
練習(xí): 有多少唯一的條形碼被檢測(cè)到？數(shù)據(jù)里多少來(lái)自真細(xì)胞池凄？為了簡(jiǎn)化計(jì)算抡驼，寫代碼排除掉少于10個(gè)分子的條形碼。
答案:
# 每一行代碼對(duì)應(yīng)上面一個(gè)問題
dim(umi_per_barcode)[1]
dim(truth)[1]
umi_per_barcode <- umi_per_barcode[umi_per_barcode[,2] > 10,]
一個(gè)找尋每個(gè)條形碼對(duì)應(yīng)的分子數(shù)突然下降的拐點(diǎn)的方法：
barcode_rank <- rank(-umi_per_barcode[,2])
plot(barcode_rank, umi_per_barcode[,2], xlim=c(1,8000))
可以看到文庫(kù)大小分布類似一條指數(shù)曲線肿仑，為了看的更清楚婶恼，對(duì)數(shù)據(jù)進(jìn)行log轉(zhuǎn)換。
log_lib_size <- log10(umi_per_barcode[,2])
plot(barcode_rank, log_lib_size, xlim=c(1,8000))
轉(zhuǎn)換后柏副，拐點(diǎn)更明顯了勾邦，我們可以手動(dòng)估計(jì)拐點(diǎn)的位置來(lái)區(qū)分背景噪音和真實(shí)細(xì)胞，但是用算法來(lái)識(shí)別要更精確割择，可重用性更強(qiáng)眷篇。
# 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)
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)
## [1] 1.0000000 0.7831707
另一種方法是構(gòu)建混合模型，找出較高和較低分布相交的位置荔泳。但是數(shù)據(jù)可能不會(huì)像假定的分布那么好：
set.seed(-92497)
# mixture model
require("mixtools")
mix <- normalmixEM(log_lib_size)
## number of iterations= 43
plot(mix, which=2, xlab2="log(mol per cell)")
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])
}
練習(xí) 使用分割點(diǎn)識(shí)別細(xì)胞并計(jì)算TPR和Recall蕉饼。
答案
cells <- umi_per_barcode[umi_per_barcode[,2] > 10^split,1]
TPR <- sum(cells %in% truth[,1])/length(cells)
Recall <- sum(cells %in% truth[,1])/length(truth[,1])
c(TPR, Recall)
第三種方法，CellRanger假設(shè)真實(shí)細(xì)胞文庫(kù)大小變化在10倍以內(nèi)玛歌，用期望的細(xì)胞數(shù)目估計(jì)區(qū)間的分布昧港。
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)
練習(xí) 用這個(gè)閾值識(shí)別細(xì)胞并計(jì)算TPR和Recall
答案
cells <- umi_per_barcode[umi_per_barcode[,2] > thresh,1]
TPR <- sum(cells %in% truth[,1])/length(cells)
Recall <- sum(cells %in% truth[,1])/length(truth[,1])
c(TPR, Recall)
最后一個(gè)工具，EmptyDrops (https://github.com/MarioniLab/DropletUtils, 目前還是測(cè)試版支子。它的輸入是所有液滴的 基因 x 細(xì)胞分子計(jì)數(shù)矩陣创肥，從具有極低reads計(jì)數(shù)的液滴中估計(jì)background RNA的分布，然后尋找與background不同的基因表達(dá)譜的液滴視為真實(shí)細(xì)胞值朋。background RNA通程局叮看起來(lái)非常類似于群體中最大的細(xì)胞的表達(dá)譜，這個(gè)方法通常也會(huì)和拐點(diǎn)結(jié)合昨登。因此趾代，EmptyDrops是唯一能夠在高度多樣化樣本中鑒定非常小的細(xì)胞的方法。
下面提供了這個(gè)方法的運(yùn)行代碼：（我們會(huì)根據(jù)官方及時(shí)更新代碼）
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)
