劉小澤寫于2020.7.7
為何取名叫“交響樂”？因?yàn)閱渭?xì)胞分析就像一個大樂團(tuán)，需要各個流程的協(xié)同配合
單細(xì)胞交響樂1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細(xì)胞交響樂2-scRNAseq從實(shí)驗(yàn)到下游簡介
單細(xì)胞交響樂3-細(xì)胞質(zhì)控
單細(xì)胞交響樂4-歸一化
單細(xì)胞交響樂5-挑選高變化基因
單細(xì)胞交響樂6-降維
單細(xì)胞交響樂7-聚類分群
單細(xì)胞交響樂8-marker基因檢測
單細(xì)胞交響樂9-細(xì)胞類型注釋
1 前言
大型的scRNA分析一般都需要整合多個批次的數(shù)據(jù)集，但每個批次的質(zhì)量可能存在不可調(diào)和的差異艳丛，例如儀器的調(diào)整政供、試劑質(zhì)量的差異。結(jié)果就導(dǎo)致了不同批次細(xì)胞中基因表達(dá)量的系統(tǒng)差異椒振，稱之為“批次效應(yīng)”昭伸。之前在單細(xì)胞交響樂4-歸一化 中也介紹了一些批次效應(yīng)的事情，它最大的隱患就是：如果它這個差異成為整個數(shù)據(jù)的主力軍澎迎，那么真正的生物學(xué)差異就難以檢測庐杨，對結(jié)果的解讀變得更難。
因此現(xiàn)在出了一些多數(shù)據(jù)整合時批次效應(yīng)處理工具夹供，早期的方法 (Ritchie et al. 2015; Leek et al. 2012)是基于線性模型灵份，它假設(shè)細(xì)胞群體的組成要么是已知的，要么是在批次間也是同類型的罩引。但很明顯這個假設(shè)有局限各吨，我們不能保證這個假設(shè)成立枝笨，于是定制的算法產(chǎn)生(Haghverdi et al. 2018; Butler et al. 2018; Lin et al. 2019)袁铐，它不再需要對細(xì)胞組成有一個先驗(yàn)知識，保證了對未知的探索過程横浑。
數(shù)據(jù)準(zhǔn)備
將使用兩個不同批次的10X PBMC數(shù)據(jù)進(jìn)行整合剔桨，它們都是從 TENxPBMCData這個數(shù)據(jù)包獲得的，分別進(jìn)行了前期的數(shù)據(jù)處理徙融。各自前期處理的一個好處洒缀，就是可以根據(jù)各自的特點(diǎn)先過濾掉一些低質(zhì)量數(shù)據(jù)（比如各自根據(jù)QC結(jié)果對低質(zhì)量細(xì)胞進(jìn)行過濾）
常規(guī)的前期處理步驟還是：質(zhì)控=》歸一化=》找HVGs=》降維=》聚類
只是注意下面的lapply用法，多個數(shù)據(jù)全部基于list處理
# 數(shù)據(jù)下載（數(shù)據(jù)我已做好欺冀，鏈接在：https://share.weiyun.com/iE0VjVD0）
library(TENxPBMCData)
all.sce <- list(
pbmc3k=TENxPBMCData('pbmc3k'),
pbmc4k=TENxPBMCData('pbmc4k'),
pbmc8k=TENxPBMCData('pbmc8k')
)
# 批量質(zhì)控
library(scater)
stats <- high.mito <- list()
for (n in names(all.sce)) {
current <- all.sce[[n]] #接下來的每步還是和之前一樣
is.mito <- grep("MT", rowData(current)$Symbol_TENx)
stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito))
high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher")
all.sce[[n]] <- current[,!high.mito[[n]]]
}
# 批量歸一化
all.sce <- lapply(all.sce, logNormCounts)
# 批量找HVGs
library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)
> length(all.hvgs[[1]])
[1] 1001
> length(all.hvgs[[2]])
[1] 1352
# 批量降維(使用了一個更高級的mapply：Apply a Function to Multiple List or Vector Arguments)
library(BiocSingular)
set.seed(10000)
all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs,
MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()),
SIMPLIFY=FALSE)
set.seed(100000)
all.sce <- lapply(all.sce, runTSNE, dimred="PCA")
set.seed(1000000)
all.sce <- lapply(all.sce, runUMAP, dimred="PCA")
# 批量聚類
for (n in names(all.sce)) {
g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA')
clust <- igraph::cluster_walktrap(g)$membership
colLabels(all.sce[[n]]) <- factor(clust)
}
# 把數(shù)據(jù)提取出來
pbmc3k <- all.sce$pbmc3k
dec3k <- all.dec$pbmc3k
pbmc4k <- all.sce$pbmc4k
dec4k <- all.dec$pbmc4k
為了方便下面的批次矯正树绩，先做幾件事
1 選出兩個數(shù)據(jù)集的基因交集，并各自取子集
這里兩個數(shù)據(jù)都是Ensembl ID隐轩，因此這個過程很簡單
universe <- intersect(rownames(pbmc3k), rownames(pbmc4k))
> length(rownames(pbmc3k));length(rownames(pbmc4k));length(universe)
[1] 32738
[1] 33694
[1] 31232
然后取子集
pbmc3k <- pbmc3k[universe,]
pbmc4k <- pbmc4k[universe,]
dec3k <- dec3k[universe,]
dec4k <- dec4k[universe,]
2 矯正批次間測序深度的影響
利用multiBatchNorm()函數(shù)饺饭，將兩個批次放在一起，重新計(jì)算log-count职车，可以提高下游批次矯正的質(zhì)量
還記得之前的size factor這個名詞嗎瘫俊？
它是為了去除一個批次數(shù)據(jù)中的各個細(xì)胞之間的文庫差異
而這里，是為了去除兩個批次之間的差異
library(batchelor)
rescaled <- multiBatchNorm(pbmc3k, pbmc4k)
> class(rescaled)
[1] "list"
pbmc3k <- rescaled[[1]]
pbmc4k <- rescaled[[2]]
3 將之前單獨(dú)挑出來的模型構(gòu)建結(jié)果再次整合
當(dāng)然還也可以把每個批次的HVGs取交集或并集悴灵，選多少基因這個沒有成文規(guī)定扛芽。
這里只是因?yàn)閟cran提供了這么一個函數(shù)，并可以提供更多一些的基因积瞒，畢竟現(xiàn)在基因多一點(diǎn)沒壞處川尖，省的后面捉襟見肘
library(scran)
combined.dec <- combineVar(dec3k, dec4k)
chosen.hvgs <- combined.dec$bio > 0
sum(chosen.hvgs)
## [1] 13431
2 檢查批次效應(yīng)
在進(jìn)行批次矯正之前茫孔，肯定要先看看有沒有批次效應(yīng)以及有多嚴(yán)重吧
最簡單的方法是把兩個數(shù)據(jù)合并筐钟，然后PCA+t-SNE
合并數(shù)據(jù)
如果簡單的cbind(pbmc3k, pbmc4k)壹将，會報(bào)錯
> cbind(pbmc3k, pbmc4k)
Error in FUN(X[[i]], ...) :
column(s) 'Symbol_TENx' in ‘mcols’ are duplicated and the data do not match
# 看似問題出現(xiàn)在mcols中的Symbol_TENx暴区，看一眼
看樣子的確存在不匹配的問題刃唤，第一行就暴露了辐真，同一個Ensembl ID的Symbol就不同楔脯。這個問題先不深究偎箫。先看看兩個數(shù)據(jù)的Ensembl 是否一致吧
> identical(mcols(pbmc3k)[,1],mcols(pbmc4k)[,1])
[1] TRUE
# 一致，那么就把它們的變成一樣就好了皆串，symbol的問題我們不考慮
rowData(pbmc3k) <- rowData(pbmc4k)
pbmc3k$batch <- "3k"
pbmc4k$batch <- "4k"
uncorrected <- cbind(pbmc3k, pbmc4k)
進(jìn)行PCA
使用的HVGs也是合并兩個數(shù)據(jù)集后得到的chosen.hvgs
library(scater)
set.seed(0010101010)
uncorrected <- runPCA(uncorrected, subset_row=chosen.hvgs,
BSPARAM=BiocSingular::RandomParam())
再看聚類分群結(jié)果
如果說淹办，這兩個批次是真的重復(fù)，差異不大的話恶复，那么分群的cluster應(yīng)該在兩個批次中包含差不多數(shù)量的細(xì)胞怜森。但是下面發(fā)現(xiàn)基本clusters主要是僅僅包含了一個批次
library(scran)
snn.gr <- buildSNNGraph(uncorrected, use.dimred="PCA")
clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab <- table(Cluster=clusters, Batch=uncorrected$batch)
tab
## Batch
## Cluster 3k 4k
## 1 1 781
## 2 0 1309
## 3 0 535
## 4 14 51
## 5 0 605
## 6 489 0
## 7 0 184
## 8 1272 0
## 9 0 414
## 10 151 0
## 11 0 50
## 12 155 0
## 13 0 65
## 14 0 61
## 15 0 88
## 16 30 0
## 17 339 0
## 18 145 0
## 19 11 3
## 20 2 36
看t-SNE結(jié)果
set.seed(1111001)
uncorrected <- runTSNE(uncorrected, dimred="PCA")
plotTSNE(uncorrected, colour_by="batch")
雖然說t-SNE圖上的點(diǎn)大小、遠(yuǎn)近說明不了問題谤牡，但是兩個批次分的這么開副硅，也提供了兩個批次的證據(jù)
當(dāng)然，還可以繼續(xù)找證據(jù)翅萤。利用上一次的細(xì)胞類型注釋恐疲，看看兩個批次中的細(xì)胞類型，是不是差異很大套么。
3 矯正批次效應(yīng)之線性回歸
bulk mRNA轉(zhuǎn)錄組中常用的矯正批次效應(yīng)方法就是線性回歸培己，對每個基因表達(dá)量擬合一個線性模型。例如limma的removeBatchEffect() (Ritchie et al. 2015) 违诗、sva的comBat() (Leek et al. 2012)漱凝。如果要使用這類方法疮蹦，就需要假設(shè)：批次間的細(xì)胞組成相同诸迟。另外的一個假設(shè)是：批次效應(yīng)的累積的，對于任何給定的基因愕乎，在不同亞群中經(jīng)過任何因素誘導(dǎo)的表達(dá)變化倍數(shù)是相同的阵苇。（其實(shí)，從這兩個假設(shè)就看出來感论，這個方法不適合我們的單細(xì)胞數(shù)據(jù)绅项，但還是要繼續(xù)了解下去）
下面將使用batchelor 包中的rescaleBatches()函數(shù)進(jìn)行處理
它也是對每個基因的log表達(dá)量進(jìn)行了線性回歸，并提高了一些運(yùn)行性能比肄。另外與removeBatchEffect()不同的是快耿，rescaleBatches()保持了數(shù)據(jù)的稀疏性，而removeBatchEffect()會破壞稀疏性
library(batchelor)
rescaled <- rescaleBatches(pbmc3k, pbmc4k)
# 依然進(jìn)行PCA芳绩、聚類
set.seed(1010101010)
rescaled <- runPCA(rescaled, subset_row=chosen.hvgs, exprs_values="corrected")
snn.gr <- buildSNNGraph(rescaled, use.dimred="PCA")
clusters.resc <- igraph::cluster_walktrap(snn.gr)$membership
tab.resc <- table(Cluster=clusters.resc, Batch=rescaled$batch)
tab.resc
## Batch
## Cluster 1 2
## 1 278 525
## 2 16 23
## 3 337 606
## 4 43 748
## 5 604 529
## 6 22 71
## 7 188 48
## 8 25 49
## 9 263 0
## 10 123 135
## 11 16 85
## 12 11 57
## 13 116 6
## 14 455 1035
## 15 6 31
## 16 89 187
## 17 3 36
## 18 3 8
## 19 11 3
效果還是有的掀亥，現(xiàn)在大部分clusters都包含了兩個批次的細(xì)胞，但是依然有存在批次差異的cluster（看cluster9妥色，在batch2中細(xì)胞數(shù)為0）搪花，表明處理有效果但不徹底 。影響效果的原因是：我們的數(shù)據(jù)違背了這個方法的假設(shè)
再做個圖看看
rescaled <- runTSNE(rescaled, dimred="PCA")
rescaled$batch <- factor(rescaled$batch)
plotTSNE(rescaled, colour_by="batch")
注意
除了這個函數(shù)，還可以嘗試更常規(guī)的regressBatches() 撮竿，只不過它的表現(xiàn)可能會更差吮便，因?yàn)闀G失數(shù)據(jù)稀疏性
4 MNN矯正
4.1 先了解一下這個方法
試想batchA中有一個細(xì)胞a，然后想根據(jù)挑選的特征基因（如HVGs）去在batchB中鑒定與a細(xì)胞空間相近的細(xì)胞們幢踏。同樣的髓需，對batchB中的細(xì)胞b 重復(fù)這個過程，也鑒定它在batchA中的近鄰房蝉。于是這個MNN就是：Mutual nearest neighbors 授账，是對一個批次的細(xì)胞找到它們在另一個批次中對應(yīng)的最相鄰細(xì)胞集合。
Mutual nearest neighbors are pairs of cells from different batches that belong in each other s set of nearest neighbors. Haghverdi et al. (2018)
在進(jìn)行批次矯正前惨驶，MNN找到的一對對細(xì)胞都是生物學(xué)狀態(tài)相似的白热，因此如果再有不同，那么就姑且認(rèn)為是外部的批次效應(yīng)導(dǎo)致的粗卜，也就方便了去除屋确。
與線性回歸方法相比，MNN方法不會假設(shè)細(xì)胞群組成相同或者事先已知续扔。MNN會自己學(xué)習(xí)細(xì)胞群的結(jié)構(gòu)并進(jìn)行估計(jì)攻臀。
4.2 還是應(yīng)用到PBMC數(shù)據(jù)
batchelor 這個包除了線性模型，還提供了MNN 的函數(shù)fastMNN()纱昧，但這個函數(shù)與單純的MNN算法還有不同刨啸。fastMNN()先進(jìn)行PCA降維，以加速下面的檢測細(xì)胞近鄰识脆。
set.seed(1000101001)
# 這里的d指的是前多少個主成分设联，k指近鄰的數(shù)量（k增大會導(dǎo)致數(shù)據(jù)集合并更激進(jìn)；一般不要設(shè)置太大灼捂，如果看到相同的細(xì)胞類型沒有在批次之間充分合并离例，可以適當(dāng)增加k）
mnn.out <- fastMNN(pbmc3k, pbmc4k, d=50, k=20, subset.row=chosen.hvgs,
BSPARAM=BiocSingular::RandomParam(deferred=TRUE))
mnn.out
## class: SingleCellExperiment
## dim: 13431 6791
## metadata(2): merge.info pca.info
## assays(1): reconstructed
## rownames(13431): ENSG00000239945 ENSG00000228463 ... ENSG00000198695
## ENSG00000198727
## rowData names(1): rotation
## colnames: NULL
## colData names(1): batch
## reducedDimNames(1): corrected
## altExpNames(0):
其中mnn.out的每列表示某個批次中的一個細(xì)胞（具體存儲在batch），行表示chosen.hvgs中的基因
head(mnn.out$batch)
## [1] 1 1 1 1 1 1
看到降維相關(guān)模塊reducedDimNames中多了一個corrected悉稠，它是根據(jù)50個主成分得到的所有6791在低維空間的坐標(biāo)
dim(reducedDim(mnn.out, "corrected"))
## [1] 6791 50
> reducedDim(mnn.out, "corrected")[1:3,1:3]
[,1] [,2] [,3]
[1,] -0.12783777 0.0409469 -0.000116194
[2,] -0.03364614 -0.1466529 0.161690348
[3,] -0.09895631 0.1219669 0.010321992
看到表達(dá)量模塊assays() 中多了一個reconstructed矩陣宫蛆，是每個細(xì)胞中每個基因矯正后的表達(dá)量
> dim(assay(mnn.out, "reconstructed"))
[1] 13431 6791
> assay(mnn.out, "reconstructed")[1:4,1:4]
<4 x 4> matrix of class LowRankMatrix and type "double":
[,1] [,2]
ENSG00000239945 -2.522191e-06 -1.851424e-06
ENSG00000228463 -6.626821e-04 -6.724341e-04
ENSG00000237094 -8.077231e-05 -8.038006e-05
ENSG00000229905 3.838135e-06 6.179994e-06
[,3] [,4]
ENSG00000239945 -1.198984e-05 -3.192269e-06
ENSG00000228463 -4.820230e-04 -6.330560e-05
ENSG00000237094 -9.630608e-05 -6.855333e-05
ENSG00000229905 5.432122e-06 -2.118592e-05
4.3 矯正后的檢查
這次函數(shù)已經(jīng)做過PCA了，那就繼續(xù)分群
library(scran)
snn.gr <- buildSNNGraph(mnn.out, use.dimred="corrected")
clusters.mnn <- igraph::cluster_walktrap(snn.gr)$membership
tab.mnn <- table(Cluster=clusters.mnn, Batch=mnn.out$batch)
tab.mnn
## Batch
## Cluster 1 2
## 1 337 606
## 2 289 542
## 3 152 181
## 4 12 4
## 5 517 467
## 6 17 19
## 7 313 661
## 8 162 118
## 9 11 56
## 10 547 1083
## 11 17 59
## 12 16 58
## 13 144 93
## 14 67 191
## 15 4 36
## 16 4 8
可以看到這里的效果不錯的猛，沒有哪個cluster在某個batch中細(xì)胞數(shù)為0
再作圖看看
library(scater)
set.seed(0010101010)
mnn.out <- runTSNE(mnn.out, dimred="corrected")
mnn.out$batch <- factor(mnn.out$batch)
plotTSNE(mnn.out, colour_by="batch")
注意
fastMNN()函數(shù)其實(shí)也自帶了檢測數(shù)據(jù)
metadata(mnn.out)$merge.info$lost.var
## [,1] [,2]
## [1,] 0.006617087 0.003315395
它的含義是：the proportion of variance within each batch that is lost during MNN correction耀盗。如果這個值很大（比如>10%），說明有些”矯枉過正“卦尊，把一些內(nèi)在的生物學(xué)差異也給矯正了叛拷。
5 更為細(xì)致的檢查
5.1 同一批次內(nèi)clusters的比較
兩個批次數(shù)據(jù)混合后進(jìn)行矯正，還是應(yīng)該保留每個批次的特性猫牡。一個很不想看到的情況是：看上去細(xì)胞混合效果很好胡诗，以為矯正的效果不錯邓线，但沒想到矯正過了，細(xì)胞中生物學(xué)異質(zhì)性也被抹除了煌恢。
我們想在整合+矯正批次后骇陈，看到矯正后的一個cluster中既包含了原來的batch1信息，也包含batch2瑰抵。當(dāng)然表現(xiàn)的差異越大越好你雌，體現(xiàn)在下圖中就是：after1（即矯正后的cluster1）對應(yīng)的PBMC3k 的clusters與PBMC 4k的clusters最好別重復(fù)，并且差別越大越好
library(pheatmap)
# 畫出二者混合后的cluster與各自之前的cluster對比
# batch 1
tab <- table(paste("after", clusters.mnn[rescaled$batch==1]),
paste("before", colLabels(pbmc3k)))
heat3k <- pheatmap(log10(tab+10), cluster_row=FALSE, cluster_col=FALSE,
main="PBMC 3K comparison", silent=TRUE)
# batch 2
tab <- table(paste("after", clusters.mnn[rescaled$batch==2]),
paste("before", colLabels(pbmc4k)))
heat4k <- pheatmap(log10(tab+10), cluster_row=FALSE, cluster_col=FALSE,
main="PBMC 4K comparison", silent=TRUE)
gridExtra::grid.arrange(heat3k[[4]], heat4k[[4]])
就以after5來說二汛，分別對應(yīng)3k混合之前的cluster3婿崭、4、7肴颊；以及4k混合前的cluster1氓栈、10、3婿着，差別蠻大授瘦，說明混合后依然保持了各自的個性
除了看圖，還能通過數(shù)據(jù)體現(xiàn)——計(jì)算Rand index
這個指標(biāo)是矯正后各個batch原來差異的保留度竟宋，這個數(shù)值越大越好（但如果某個批次矯正方法不給力提完，這個值也是大的，所以還是要慎重對待）
library(fossil)
# batch1
ri3k <- rand.index(as.integer(clusters.mnn[rescaled$batch==1]),
as.integer(colLabels(pbmc3k)))
ri3k
## [1] 0.9335226
# batch2
ri4k <- rand.index(as.integer(clusters.mnn[rescaled$batch==2]),
as.integer(colLabels(pbmc4k)))
ri4k
## [1] 0.9575746
5.2 利用marker基因檢查
目的就是看數(shù)據(jù)整合+矯正后丘侠，原始batch的內(nèi)部結(jié)構(gòu)是否還完整
首先分別從原來的兩個批次數(shù)據(jù)中尋找marker基因徒欣，因?yàn)閙arker基因的存在就保證了這個批次數(shù)據(jù)集中”生物學(xué)的有趣性“存在
stats3 <- pairwiseWilcox(pbmc3k, direction="up")
markers3 <- getTopMarkers(stats3[[1]], stats3[[2]], n=10)
stats4 <- pairwiseWilcox(pbmc4k, direction="up")
markers4 <- getTopMarkers(stats4[[1]], stats4[[2]], n=10)
然后用雙方的marker基因合集，作為矯正的HVGs蜗字，看看算法會不會把這部分”有趣性“去除
marker.set <- unique(unlist(c(unlist(markers3), unlist(markers4))))
length(marker.set)
## [1] 314
set.seed(1000110)
mnn.out2 <- fastMNN(pbmc3k, pbmc4k, subset.row=marker.set,
BSPARAM=BiocSingular::RandomParam(deferred=TRUE))
最后作圖看看
mnn.out2 <- runTSNE(mnn.out2, dimred="corrected")
gridExtra::grid.arrange(
plotTSNE(mnn.out2[,mnn.out2$batch==1], colour_by=I(colLabels(pbmc3k))),
plotTSNE(mnn.out2[,mnn.out2$batch==2], colour_by=I(colLabels(pbmc4k))),
ncol=2
)
看樣子根據(jù)這些marker基因進(jìn)行的批次矯正打肝，并沒有去掉marker基因背后的生物學(xué)差異，因?yàn)閮蓚€批次數(shù)據(jù)還是能分出來群的秽澳，并且分群結(jié)果看似還不錯
6 矯正后數(shù)據(jù)的應(yīng)用
數(shù)據(jù)整合并矯正后闯睹，對合并后數(shù)據(jù)的cluster分析即可，不用再對每個batch的各個cluster單獨(dú)分析担神。另一個好處是：batch整合后作為一整個數(shù)據(jù)集，細(xì)胞數(shù)量增加始花，相當(dāng)于增加了群體結(jié)構(gòu)的分辨率妄讯，方便下游的細(xì)胞層面分析（比如后面的軌跡推斷）。
可能你還想利用合并后的數(shù)據(jù)進(jìn)行基因?qū)用娴姆治隹嵯绮町惙治鰉arker基因鑒定亥贸。但一般不推薦，因?yàn)樗惴▽Χ鄠€批次進(jìn)行合并時浇垦，不會保留每個基因表達(dá)的差異炕置。如果要探索基因表達(dá)相關(guān)，最好還是要在未矯正數(shù)據(jù)基礎(chǔ)上進(jìn)行，并且還要把批次信息封鎖掉（還記得之前findMarkers中的block參數(shù)嗎朴摊？）
# 例如
m.out <- findMarkers(uncorrected, clusters.mnn, block=uncorrected$batch,
direction="up", lfc=1, row.data=rowData(uncorrected)[,3,drop=FALSE])
歡迎關(guān)注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩默垄，打造生信星球，想讓它成為一個不拽術(shù)語甚纲、通俗易懂的生信知識平臺口锭。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
