| description |
|---|
刘小泽写于2020.7.21 |
2016年很多研究都对人类胰腺细胞的scRNA很感兴趣，也因此发表了很多文章：(Muraro et al. 2016; Grun et al. 2016; Lawlor et al. 2017; Segerstolpe et al. 2016)。这些不同作者不同技术手段得到的数据，也给数据整合带来了不小的挑战。相比于之前的PBMC数据整合，这里更为复杂，因为包含的建库、测序方法多样，供体的类型、数量更是不一致。
首先拿技术相似的两套数据来做：分别是 Muraro et al. 2016; Grun et al. 2016，采用了CEL-seq 和 CEL-seq2
load('final.sce.grun.Rdata')
load('final.sce.muraro.RData')
final.sce.grun
# class: SingleCellExperiment
# dim: 17548 1063
# metadata(0):
# assays(2): counts logcounts
# rownames(17548): ENSG00000268895 ENSG00000121410 ...
# ENSG00000074755 ENSG00000036549
# rowData names(2): symbol chr
# colnames(1063): D2ex_1 D2ex_2 ... D17TGFB_94
# D17TGFB_95
# colData names(3): donor sample sizeFactor
# reducedDimNames(0):
# altExpNames(1): ERCC
final.sce.muraro
# class: SingleCellExperiment
# dim: 16940 2299
# metadata(0):
# assays(2): counts logcounts
# rownames(16940): ENSG00000268895 ENSG00000121410 ...
# ENSG00000159840 ENSG00000074755
# rowData names(2): symbol chr
# colnames(2299): D28-1_1 D28-1_2 ... D30-8_93
# D30-8_94
# colData names(4): label donor plate sizeFactor
# reducedDimNames(0):
# altExpNames(1): ERCC首先获得交集基因
universe <- intersect(rownames(final.sce.grun), rownames(final.sce.muraro))
> nrow(final.sce.grun);nrow(final.sce.muraro);length(universe)
[1] 17548
[1] 16940
[1] 15974对数据集取子集
sce.grun2 <- final.sce.grun[universe,]
sce.muraro2 <- final.sce.muraro[universe,]既然是经过处理后的数据，那么就略过了之前介绍的质控步骤
首先测序深度导致的文库大小差异是批次效应的一个重要来源，因此可以先对不同的批次进行文库矫正。会以文库最小的批次为基准，对其他批次进行文库归一化。最后返回一个列表
使用一个归一化函数：
multiBatchNorm，它应用的就是最简单的library size normalization归一化方法： Perform scaling normalization within each batch to provide comparable results to the lowest-coverage batch.
既然是要处理文库差异，那就先看看各自原本的文库大小
summary(colSums(logcounts(sce.muraro2)))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 4433 8041 8680 8558 9193 10594
summary(colSums(logcounts(sce.grun2)))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 2347 4141 4565 4538 4969 5992然后进行处理，可以看看前后的变化，就明白了这个函数做了什么事情
library(batchelor)
normed.pancreas <- multiBatchNorm(sce.grun2, sce.muraro2)
sce.grun3 <- normed.pancreas[[1]]
sce.muraro3 <- normed.pancreas[[2]]
summary(colSums(logcounts(sce.muraro3)))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 3072 4735 4994 4954 5204 5864
summary(colSums(logcounts(sce.grun3)))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 2347 4141 4565 4538 4969 5992看到混合以后的sce.muraro3相对于之前独立的sce.muraro2的变化了吧
首先对表达量变化模型取子集
## 原来的sce.grun模型
library(scran)
block1 <- paste0(final.sce.grun$sample, "_", final.sce.grun$donor)
dec.grun <- modelGeneVarWithSpikes(final.sce.grun, spikes="ERCC", block=block1)
# 取子集
dec.grun2 <- dec.grun[universe,]
## 原来的sce.muraro模型
block2 <- paste0(final.sce.muraro$plate, "_", final.sce.muraro$donor)
dec.muraro <- modelGeneVarWithSpikes(final.sce.muraro, "ERCC", block=block2)
# 取子集
dec.muraro2 <- dec.muraro[universe,]之后组合两组的结果
使用
combineVar，它的作用是：Combine the results of multiple variance decompositions, usually generated for the same genes across separate batches of cells.
library(scran)
combined.pan <- combineVar(dec.grun2, dec.muraro2)
# 把更有可能代表生物差异的基因选出来，用于下游的PCA和聚类
chosen.genes <- rownames(combined.pan)[combined.pan$bio > 0]之前在：3.8 批次效应处理 中介绍过：
bulk mRNA转录组中常用的矫正批次效应方法就是线性回归，对每个基因表达量拟合一个线性模型。例如limma的
removeBatchEffect()(Ritchie et al. 2015) 、sva的comBat()( Leek et al. 2012)。如果要使用这类方法，就需要假设：批次间的细胞组成相同。另外的一个假设是：批次效应的累积的，对于任何给定的基因，在不同亚群中经过任何因素诱导的表达变化倍数是相同的。（其实，从这两个假设就看出来，这个方法不适合我们的单细胞数据，但还是要继续了解下去）
先来看看基于线性回归的rescaleBatches()
它也是对每个基因的log表达量进行了线性回归，并提高了一些运行性能。另外与removeBatchEffect()不同的是，rescaleBatches()保持了数据的稀疏性，而removeBatchEffect()会破坏稀疏性
library(scater)
rescaled.pancreas <- rescaleBatches(sce.grun2, sce.muraro2)
set.seed(100101)
rescaled.pancreas <- runPCA(rescaled.pancreas, subset_row=chosen.genes,
exprs_values="corrected")
rescaled.pancreas <- runTSNE(rescaled.pancreas, dimred="PCA")
plotTSNE(rescaled.pancreas, colour_by="batch")不过结果并不尽如人意，两个批次还是分得很开，表明处理有效果但不彻底 。影响效果的原因是：我们的数据违背了这个方法的假设
再来使用fastMNN()
与线性回归方法相比，MNN方法不会假设细胞群组成相同或者事先已知。MNN会自己学习细胞群的结构并进行估计.
可能之前听过：mnnCorrect()这个方法，它是Haghverdi et al. (2018) 提出来的，之前也介绍过：单细胞转录组数据校正批次效应实战
它和fastMNN()原理类似，但速度会慢很多，总之它们的不同可以概括为：
For scRNA-seq data,
fastMNN()tends to be both faster and better at achieving a satisfactory merge.mnnCorrect()is mainly provided here for posterity’s sake, though it is more robust thanfastMNN()to certain violations of the orthogonality assumptions.
set.seed(1011011)
mnn.pancreas <- fastMNN(sce.grun2, sce.muraro2, subset.row=chosen.genes)
snn.gr <- buildSNNGraph(mnn.pancreas, use.dimred="corrected")
clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab <- table(Cluster=clusters, Batch=mnn.pancreas$batch)
tab
## Batch
## Cluster 1 2
## 1 239 281
## 2 312 257
## 3 200 837
## 4 56 193
## 5 37 1
## 6 24 108
## 7 109 391
## 8 63 80
## 9 18 115
## 10 0 17
## 11 5 19再做个图
mnn.pancreas <- runTSNE(mnn.pancreas, dimred="corrected")
plotTSNE(mnn.pancreas, colour_by="batch")效果有明显改进
前面是将一个CEL-seq和一个CEL-seq2数据整合，总的说来还不是很复杂
但这里，会再加上两个数据，分别来自Lawlor et al. 2017 和 Segerstolpe et al. 2016 的数据进行整合，这样四个数据就会包括不同的技术、不同的UMI、不同的表达量、更加不同的供体
rm(list = ls())
load('final.sce.grun.RData')
load('final.sce.muraro.RData')
load('final.sce.lawlor.RData')
load('final.sce.seger.RData')
sce.grun=final.sce.grun
sce.muraro=final.sce.muraro
sce.grun
# class: SingleCellExperiment
# dim: 17548 1063
# metadata(0):
# assays(2): counts logcounts
# rownames(17548): ENSG00000268895 ENSG00000121410 ...
# ENSG00000074755 ENSG00000036549
# rowData names(2): symbol chr
# colnames(1063): D2ex_1 D2ex_2 ... D17TGFB_94
# D17TGFB_95
# colData names(3): donor sample sizeFactor
# reducedDimNames(0):
# altExpNames(1): ERCC
sce.muraro
# class: SingleCellExperiment
# dim: 16940 2299
# metadata(0):
# assays(2): counts logcounts
# rownames(16940): ENSG00000268895 ENSG00000121410 ...
# ENSG00000159840 ENSG00000074755
# rowData names(2): symbol chr
# colnames(2299): D28-1_1 D28-1_2 ... D30-8_93
# D30-8_94
# colData names(4): label donor plate sizeFactor
# reducedDimNames(0):
# altExpNames(1): ERCC
sce.lawlor
# class: SingleCellExperiment
# dim: 26616 604
# metadata(0):
# assays(2): counts logcounts
# rownames(26616): ENSG00000229483 ENSG00000232849 ...
# ENSG00000251576 ENSG00000082898
# rowData names(2): SYMBOL SEQNAME
# colnames(604): 10th_C11_S96 10th_C13_S61 ...
# 9th-C96_S81 9th-C9_S13
# colData names(9): title age ... Sex sizeFactor
# reducedDimNames(0):
# altExpNames(0):
sce.seger
# class: SingleCellExperiment
# dim: 25454 2090
# metadata(0):
# assays(2): counts logcounts
# rownames(25454): ENSG00000118473 ENSG00000142920 ...
# ENSG00000278306 eGFP
# rowData names(2): symbol refseq
# colnames(2090): HP1502401_H13 HP1502401_J14 ...
# HP1526901T2D_N8 HP1526901T2D_A8
# colData names(4): CellType Donor Quality sizeFactor
# reducedDimNames(0):
# altExpNames(1): ERCC首先获得交集基因
all.sce <- list(Grun=sce.grun, Muraro=sce.muraro,
Lawlor=sce.lawlor, Seger=sce.seger)
universe <- Reduce(intersect, lapply(all.sce, rownames))
nrow(sce.grun);nrow(sce.muraro);nrow(sce.lawlor);nrow(sce.seger);length(universe)
# [1] 17548
# [1] 16940
# [1] 26616
# [1] 25454
# [1] 15231 #共有基因再分别取子集
all.sce <- lapply(all.sce, "[", i=universe,)normed.pancreas <- do.call(multiBatchNorm, all.sce)首先获得各个数据集的表达量变化模型
library(scran)
# sce.grun
block1 <- paste0(sce.grun$sample, "_", sce.grun$donor)
dec.grun <- modelGeneVarWithSpikes(sce.grun, spikes="ERCC", block=block1)
# sce.muraro
block2 <- paste0(sce.muraro$sample, "_", sce.muraro$donor)
dec.muraro <- modelGeneVarWithSpikes(sce.muraro, spikes="ERCC", block=block2)
# sce.lawlor：没有ERCC信息，就用modelGeneVar()
dec.lawlor <- modelGeneVar(sce.lawlor, block=sce.lawlor$`islet unos id`)
# sce.seger
for.hvg <- sce.seger[,librarySizeFactors(altExp(sce.seger)) > 0
& sce.seger$Donor!="AZ"]
dec.seger <- modelGeneVarWithSpikes(for.hvg, "ERCC", block=for.hvg$Donor)
# 整合起来
all.dec <- list(Grun=dec.grun, Muraro=dec.muraro,
Lawlor=dec.lawlor, Seger=dec.seger)再取子集
all.dec <- lapply(all.dec, "[", i=universe,)再组合四组的结果
combined.pan <- do.call(combineVar, all.dec)
# 把更有可能代表生物差异的基因选出来，用于下游的PCA和聚类
chosen.genes <- rownames(combined.pan)[combined.pan$bio > 0]fastMNN也包含了PCA的操作
set.seed(1011110)
mnn.pancreas <- fastMNN(normed.pancreas)snn.gr <- buildSNNGraph(mnn.pancreas, use.dimred="corrected", k=20)
clusters <- igraph::cluster_walktrap(snn.gr)$membership
clusters <- factor(clusters)
tab <- table(Cluster=clusters, Batch=mnn.pancreas$batch)
tab
## Batch
## Cluster Grun Lawlor Muraro Seger
## 1 77 33 677 211
## 2 311 26 254 383
## 3 104 244 390 180
## 4 225 16 246 138
## 5 125 203 158 108
## 6 56 17 196 109
## 7 0 0 0 43
## 8 0 0 0 42
## 9 0 2 17 4
## 10 69 12 82 16
## 11 24 19 108 55
## 12 0 0 0 50
## 13 27 0 1 0
## 14 22 6 34 49
## 15 18 18 117 157
## 16 0 0 0 208
## 17 0 0 0 26
## 18 0 0 0 108
## 19 0 0 0 186
## 20 5 8 19 17看到一个批次中都包含了很多clusters，说明数据整合的效果还不错，批次效应没有很强；当然有些clusters只显示在了seger批次中（比如7、8、12、16、17、18、19），那究竟这些clusters到底是不是seger数据特有的细胞类型呢？这个还有待考证
注意其中使用到了一个很有趣的函数
I()，简单的一个字母，它是base包里的函数
因为我们是根据mnn.pancreas进行作图的，但clusters这个向量是根据mnn.pancreas创建的，但又不直接存在于mnn.pancreas（不像batch一样存在于mnn.pancreas中）。
因此要从外部把它导入到作图函数中，就可以用这个I()
mnn.pancreas <- runTSNE(mnn.pancreas, dimred="corrected")
gridExtra::grid.arrange(
plotTSNE(mnn.pancreas, colour_by="batch", text_by=I(clusters)),
plotTSNE(mnn.pancreas, colour_by=I(clusters), text_by=I(clusters)),
ncol=2
)
上面我们粗略根据四个数据集看了下批次效应，发现Seger这个数据还有点特殊，因为很多clusters只存在于Seger中
批次效应的来源除了表面上的4个数据集整合，还有一个重点考虑对象是：供体的种类
看一下来自各个数据中供体的批次效应
首先检查一下各个数据的供体信息，看到Seger最多，因此它的风险也最大。但这个怀疑到底对不对，还要做个图看看
seger.donors <- donors
seger.donors[mnn.pancreas$batch!="Seger"] <- NA
grun.donors <- donors
grun.donors[mnn.pancreas$batch!="Grun"] <- NA
lawlor.donors <- donors
lawlor.donors[mnn.pancreas$batch!="Lawlor"] <- NA
muraro.donors <- donors
muraro.donors[mnn.pancreas$batch!="Muraro"] <- NA
gridExtra::grid.arrange(
plotTSNE(mnn.pancreas, colour_by=I(muraro.donors))+ ggtitle('muraro.donors'),
plotTSNE(mnn.pancreas, colour_by=I(lawlor.donors))+ ggtitle('lawlor.donors'),
plotTSNE(mnn.pancreas, colour_by=I(grun.donors))+ ggtitle('grun.donors'),
plotTSNE(mnn.pancreas, colour_by=I(seger.donors))+ ggtitle('seger.donors'),
ncol=2
)看到图中Seger的供体各个细胞最为分散，因此它的供体批次效应是最强的
虽说供体也是生物信息，但它对于后续的细胞类型注释没有直接的帮助，相反还会产生混淆的作用（比如前面看到很多clusters只存在于Seger中，说不定就是由于Seger的供体批次效应导致的）
因此，一个更为谨慎的操作是：除了去除数据集之间的批次效应以外，还要将每个数据内部的供体信息作为另一个批次效应，分别处理掉
将原来的4个分开的数据聚合在一起，使用noCorrect 进行简单的聚合
它的含义是：This function is effectively equivalent to cbinding the matrices together without any correction.
combined <- noCorrect(normed.pancreas)
assayNames(combined) <- "logcounts"
combined$donor <- donors
对这个数据进行批次矫正，但首先要把数据批次和供体批次分开
donors.per.batch <- split(combined$donor, combined$batch)
# 获得每个数据批次下的供体批次
donors.per.batch <- lapply(donors.per.batch, unique)
donors.per.batch
## $Grun
## [1] "D2" "D3" "D7" "D10" "D17"
##
## $Lawlor
## [1] "ACIW009" "ACJV399" "ACCG268" "ACCR015A" "ACEK420A" "ACEL337" "ACHY057"
## [8] "ACIB065"
##
## $Muraro
## [1] "D28" "D29" "D31" "D30"
##
## $Seger
## [1] "HP1502401" "HP1504101T2D" "AZ" "HP1508501T2D" "HP1506401"
## [6] "HP1507101" "HP1509101" "HP1504901" "HP1525301T2D" "HP1526901T2D"依然是使用fastMNN
set.seed(1010100)
# batch信息使用全部的供体
# 增加一步指定weights：可以理解为哪些供体属于哪个数据集
multiout <- fastMNN(combined, batch=combined$donor,
subset.row=chosen.genes, weights=donors.per.batch)
# 将两大批次信息记录在新的矫正结果multiout中
multiout$dataset <- combined$batch
multiout$donor <- multiout$batch检查一下聚类结果
从下面的结果中可以看到单独属于Seger的cluster没有了
library(scater)
g <- buildSNNGraph(multiout, use.dimred=1, k=20)
clusters <- igraph::cluster_walktrap(g)$membership
tab <- table(clusters, multiout$dataset)
tab
##
## clusters Grun Lawlor Muraro Seger
## 1 246 20 278 187
## 2 200 239 835 862
## 3 171 254 473 294
## 4 315 27 260 387
## 5 57 17 193 108
## 6 24 18 107 55
## 7 26 0 0 0
## 8 5 9 19 17
## 9 19 19 118 176
## 10 0 1 16 4作图结果也发现数据混合更理想了
multiout <- runTSNE(multiout, dimred="corrected")
gridExtra::grid.arrange(
plotTSNE(multiout, colour_by="dataset", text_by=I(clusters)),
plotTSNE(multiout, colour_by=I(seger.donors)),
ncol=2
)
最后，和已发表的细胞类型做对比
由于这些数据都已发表，数据集中也包含了最后作者注释的细胞类型（sce.grun除外）
因此，可以将我们自己整合后又矫正的分群，与发表的细胞分群对比，来说明批次处理质量
# 获得已发表的细胞类型信息
proposed <- c(rep(NA, ncol(sce.grun)),
sce.muraro$label,
sce.lawlor$`cell type`,
sce.seger$CellType)
看到其中大小写参差不齐，可以全变成小写
proposed <- tolower(proposed)
# 并根据原文章修改一下细胞类型名称
proposed[proposed=="gamma/pp"] <- "gamma"
proposed[proposed=="pp"] <- "gamma"
proposed[proposed=="duct"] <- "ductal"
proposed[proposed=="psc"] <- "stellate"最后检查一下
table(proposed, clusters)
看到，我们处理完批次效应后的分群结果可以比较好的匹配到真实的细胞类型，因此说明了这里使用的批次矫正的方法的力度刚刚好