劉小澤寫于2020.7.15
為何取名叫“交響樂”挤茄？因為單細(xì)胞分析就像一個大樂團(tuán)德澈，需要各個流程的協(xié)同配合
單細(xì)胞交響樂1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細(xì)胞交響樂2-scRNAseq從實驗到下游簡介
單細(xì)胞交響樂3-細(xì)胞質(zhì)控
單細(xì)胞交響樂4-歸一化
單細(xì)胞交響樂5-挑選高變化基因
單細(xì)胞交響樂6-降維
單細(xì)胞交響樂7-聚類分群
單細(xì)胞交響樂8-marker基因檢測
單細(xì)胞交響樂9-細(xì)胞類型注釋
單細(xì)胞交響樂10-數(shù)據(jù)集整合后的批次矯正
單細(xì)胞交響樂11-多樣本間差異分析
1 前言
scRNA中，doublets指的就是一個文庫中存在兩個細(xì)胞的情況译隘。一般是由于技術(shù)誤差導(dǎo)致的（比如細(xì)胞分選亲桥、捕獲），尤其在基于液滴的技術(shù)中（例如10X）比較突出固耘。它干擾了對單個細(xì)胞表達(dá)量以及細(xì)胞形態(tài)的判斷题篷，比如一個液滴中有兩個細(xì)胞，會通過誤導(dǎo)我們這個細(xì)胞可能處于分化的“過渡態(tài)”厅目。因此檢測和去除這部分的影響至關(guān)重要番枚。
一個比較通用的檢測方法是：單純根據(jù)表達(dá)量(Dahlin et al. 2018)法严。下面將利用一個10X數(shù)據(jù)來展示兩種檢測方法，它們的主要區(qū)別是我們是否需要提前知道分群的信息
數(shù)據(jù)準(zhǔn)備
#--- 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):
或者直接加載做好的數(shù)據(jù)：https://share.weiyun.com/ReZZnAMw
2 兩種檢測方法
2.1 基于分群結(jié)果的檢測
使用doubletCluster() 將任一cluster與其他另外兩個clusters的表達(dá)量進(jìn)行比較葫笼，即3個cluser為一組深啤，其中一個是query，另外兩個是source路星∷萁郑基于的原假設(shè)是：query cluster的細(xì)胞中如果包含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群，所以下面的結(jié)果也是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
返回的結(jié)果包括：
- 基因數(shù)量N：query cluster與另外兩個source cluster相比特有的差異基因堤尾，它的數(shù)量多少可以為拒絕原假設(shè)提供依據(jù)。這個基因數(shù)量越少迁客，query越有可能是doublets
- 文庫大小比例 lib.size：每個source cluster的各個細(xì)胞文庫大小中位數(shù) / query cluster中的各個細(xì)胞文庫大小中位數(shù)郭宝，因此兩個source就對應(yīng)兩個lib.size。我們知道doublets是一個文庫包含兩個細(xì)胞掷漱，因此它會比單個細(xì)胞的文庫更大粘室。這個值越小，query越可能是doublets
- 占全部細(xì)胞的百分比 prop：query中的細(xì)胞數(shù)量占全部細(xì)胞的百分比切威，一般這個值小于5%，不過也取決于10X機(jī)器上樣的數(shù)量
最后主要還是看N的數(shù)量丙号，可以將這個結(jié)果按照N排個序
library(scater)
chosen.doublet <- rownames(dbl.out)[isOutlier(dbl.out$N,
type="lower", log=TRUE)]
chosen.doublet
## [1] "6"
挑出來懷疑對象先朦，可以對cluster6進(jìn)一步檢查
比如找到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）中也有類似的表達(dá)模式喳魏？
另外，基于背景知識怀薛，沒有細(xì)胞會同時表達(dá)basal cells (Acta2) and alveolar cells (Csn2) 這兩個基因刺彩，但是看到cluster6中這兩個基因表達(dá)量都較高，再一次證明了我們的假設(shè)：cluster6是一個doublet混合體枝恋，而不是純粹的一個細(xì)胞類型
plotExpression(sce.mam, features=c("Acta2", "Csn2"),
x="label", colour_by="label")
注意
從上面也能看到创倔，doubletCluster()的優(yōu)點就是方便操作，并且結(jié)果比較好理解焚碌。一旦有懷疑的cluster畦攘，就可以立即檢查。但缺點是高度依賴分群的質(zhì)量十电，如果分群不好知押，那么這個結(jié)果的可信度就會大打折扣叹螟。另外，如果真的有某個亞群中細(xì)胞數(shù)量很少台盯，帶來的結(jié)果就是：N小罢绽，讓這個亞群很有可能成為懷疑對象。
不過静盅，隨著scRNA的技術(shù)改進(jìn)良价，這個doublets情況會逐漸好轉(zhuǎn)
2.2 基于模擬推斷的檢測
將利用來自scran包的doubletCells() ，它基于的假設(shè)是：模擬的doublets和真實的doublets接近
它的算法是：
- 隨機(jī)選取兩個原始單細(xì)胞表達(dá)量温亲，加和棚壁，當(dāng)做一個模擬的doublet，這樣操作上千次
- 對每一個原始細(xì)胞栈虚，看看它附近有多少的模擬的doublet袖外，并計算密度
- 對每一個原始細(xì)胞，同時計算它附近的其他原始細(xì)胞的數(shù)量魂务，并計算密度
- 對每一個細(xì)胞曼验，計算兩個密度的比值，作為“doublet score”
為了加快密度的計算粘姜，這個函數(shù)會進(jìn)行一個PCA以及l(fā)og轉(zhuǎn)換
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
然后把計算結(jié)果的doublet score畫出來鬓照，數(shù)值較高的細(xì)胞聚成一團(tuán)
sce.mam$DoubletScore <- log10(dbl.dens+1)
plotTSNE(sce.mam, colour_by="DoubletScore")
同時，結(jié)合之前doubletCluster()的結(jié)果看一眼：多么明顯的cluster6孤紧！
注意
doubletCells() 的優(yōu)點就是不需要依賴分群結(jié)果豺裆，降低了分群的影響。缺點是需要對doublets的模擬更精確号显，真的要保證它和真實的情況接近臭猜。
另外簡單去掉那些doublet scores較高的細(xì)胞有時也是不夠的，例如一個潛在的doublet cluster中只有一小部分的分值高押蚤，去掉它們剩下的細(xì)胞依然會干擾判斷蔑歌。那么問題來了，怎么樣才叫高的doublet scores呢揽碘？是不是得設(shè)置一個閾值次屠？就像剛才這種情況，如果把閾值降低會不會就多包括了一些潛在的doublets呢雳刺？其實這也是生信中的一個比較頭疼的問題劫灶，沒有一個固定的閾值或者范圍，一切都是相對的掖桦。
比較推薦的方法是將doubletCells() 的結(jié)果再用分群展示出來浑此，就像上面的圖，可以更直觀看到那些細(xì)胞影響較大滞详。
小結(jié)
- 檢測doublet操作必須要求數(shù)據(jù)來自同一批次凛俱，因為doublet也不會來自兩次捕獲的細(xì)胞紊馏，它畢竟是一個文庫，只是包含兩個細(xì)胞蒲犬，因此也不需要擔(dān)心檢測doublet時怎么去除批次效應(yīng)之類的問題朱监。最好還是在分析前最好先清楚實驗設(shè)計
- 如果數(shù)據(jù)中包含了細(xì)胞分化軌跡的信息，doublet就不好判斷了原叮，因為處于變化狀態(tài)的細(xì)胞就很像doublet赫编，容易被誤認(rèn)
歡迎關(guān)注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球奋隶，想讓它成為一個不拽術(shù)語擂送、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
