摘要：不同的测序平台或同一平台不同样本得到的测序数据需要去除批次效应后进行整合分析，拟采用Seurat2程序包对CelSeq，CelSeq2，Fluidigm C1，SMART-Seq2四种测序技术产生的胰岛单细胞测序数据进行整合分析并去除批次效应。
1 数据获取
从https://satijalab.org/seurat/pancreas_integration_label_transfer.html处下载得到四种测序技术的胰岛单细胞测序数据进行后续分析。
2 分析过程
2.1 数据预处理
library(Seurat)
### 构建Seurat对象
pancreas.data <- readRDS(file = "pancreas_expression_matrix.rds")
metadata <- readRDS(file = "pancreas_metadata.rds")
pancreas <- CreateSeuratObject(raw.data = pancreas.data, meta.data = metadata)
### 分离4种技术测序数据，构建独立对象
pancreas.list <- SplitObject(object = pancreas,attribute.1 = "tech")
celseq <- pancreas.list$celseq
celseq2 <- pancreas.list$celseq2
fluidigmc1 <- pancreas.list$fluidigmc1
smartseq2 <- pancreas.list$smartseq2
### QC and normalize, preprocessing
# for celseq data
celseq <- FilterCells(celseq, subset.names = "nGene", low.thresholds = 1750)
celseq <- NormalizeData(celseq)
celseq <- FindVariableGenes(celseq, do.plot = F, display.progress = F)
celseq <- ScaleData(celseq)
# for celseq2 data
celseq2 <- FilterCells(celseq2, subset.names = "nGene", low.thresholds = 2500)
celseq2 <- NormalizeData(celseq2)
celseq2 <- FindVariableGenes(celseq2, do.plot = F, display.progress = F)
celseq2 <- ScaleData(celseq2)
# for fluidigmc1 data, already filtered and no need to FilterCells again
fluidigmc1 <- NormalizeData(fluidigmc1)
fluidigmc1 <- FindVariableGenes(fluidigmc1, do.plot = F, display.progress = F)
fluidigmc1 <- ScaleData(fluidigmc1)
# for smartseq2 data
smartseq2 <- FilterCells(smartseq2, subset.names = "nGene", low.thresholds = 2500)
smartseq2 <- NormalizeData(smartseq2)
smartseq2 <- FindVariableGenes(smartseq2, do.plot = F, display.progress = F)
smartseq2 <- ScaleData(smartseq2)
### 对数据进行采样，每个数据集随机选取500个细胞进行下游分析，以减少内存消耗，缩短分析时间
celseq <- SetAllIdent(celseq, id = "tech")
celseq <- SubsetData(celseq, max.cells.per.ident = 500, random.seed = 1)
celseq2 <- SetAllIdent(celseq2, id = "tech")
celseq2 <- SubsetData(celseq2, max.cells.per.ident = 500, random.seed = 1)
fluidigmc1 <- SetAllIdent(fluidigmc1, id = "tech")
fluidigmc1 <- SubsetData(fluidigmc1, max.cells.per.ident = 500, random.seed = 1)
smartseq2 <- SetAllIdent(smartseq2, id = "tech")
smartseq2 <- SubsetData(smartseq2, max.cells.per.ident = 500, random.seed = 1)
# 将预处理好的数据存储备用
save(celseq, celseq2, fluidigmc1, smartseq2, file='splited.4.seurat.object.RData')
2.2 不进行批次效应的矫正，直接将数据汇总之后进行分析
# Merge Seurat objects
gcdata <- MergeSeurat(celseq, celseq2, do.normalize=F)
gcdata <- MergeSeurat(gcdata, fluidigmc1, do.normalize=F)
gcdata <- MergeSeurat(gcdata, smartseq2, do.normalize=F)
# 查看不同4种技术之间在测得基因数量上的差异，如图1所示
VlnPlot(gcdata, "nGene", group.by = "tech")
图1. 不同数据集单个细胞测得基因数量.png
# 对汇总后的数据进行矫正和归一化
gcdata <- NormalizeData(object = gcdata)
gcdata <- FindVariableGenes(gcdata, do.plot = F, display.progress = F)
gcdata <- ScaleData(gcdata, genes.use = gcdata@var.genes)
# PCA降维
gcdata <- RunPCA(gcdata, pc.genes = gcdata@var.genes,
pcs.compute = 30, do.print = TRUE, pcs.print = 5, genes.print = 5)
# 查看细胞在PCA前两个维度空间中的分布，观察有无批次效应，如图2所示
# 可以用看出细胞分布受批次效应影响较大，尤其是Fluidigm C1数据集与其他数据集存在明显偏离
DimPlot(gcdata, reduction.use = "pca", dim.1 = 1, dim.2 = 2, group.by = "tech")
图2. 未去除批次效应下细胞PCA降维分布.png
# 聚类并进行tSNE降维
gcdata <- FindClusters(gcdata, reduction.type = "pca", dims.use = 1:20,
print.output = F, save.SNN = T, force.recalc = T, random.seed = 100)
gcdata <- RunTSNE(gcdata, dims.use = 1:20, do.fast = T, reduction.use = "pca", perplexity = 30)
# 在tSNE降维空间查看细胞聚类结果与批次分布，如图3所示。
# 可看出细胞聚类结果受到批次效应影响较大，同一细胞群基本来自同一批次
p1 <- DimPlot(gcdata, reduction.use = "tsne", dim.1 = 1, dim.2 = 2, group.by = "res.0.8")
p2 <- DimPlot(gcdata, reduction.use = "tsne", dim.1 = 1, dim.2 = 2, group.by = "tech")
plot_grid(p1, p2)
图3. 未去除批次效应下细胞聚类结果.png
# 根据Adjusted rand index（ARI）系数对批次效应进行评价
# ARI对当前聚类结果与4种技术标签进行比较，ARI值在-1到1之间
# ARI接近于1说明两种聚类标签具有很高的相似性
# ARI接近于0说明一种聚类标签相对于另一种为随机产生
# ARI小于0说明两种聚类标签相似性甚至低于随机水平(此处存疑)
# 具体计算过程见https://davetang.org/muse/2017/09/21/adjusted-rand-index/
# 自然，我们希望当前聚类结果与4种技术标签互相独立，ARI值应接近于0
# 结果结果显示ARI为0.3430069，说明当前聚类结果受到批次效应影响
ari <- dplyr::select(gcdata@meta.data, tech, res.0.8)
ari$tech <- plyr::mapvalues(ari$tech, from = c("celseq", "celseq2", "fluidigmc1", "smartseq2"), to = c(0, 1, 2, 3))
mclust::adjustedRandIndex(as.numeric(ari$tech), as.numeric(ari$res.0.8))
# 保存当前结果
save(gcdata, file = “gcdata.Rdata”)
2.3 通过CCA算法进行批次效应的校正后进行下游分析
# 鉴定出至少在两个数据集中同时出现的高变异基因，使用这些基因进行CCA分析
ob.list <- list(celseq, celseq2, fluidigmc1, smartseq2)
genes.use <- c()
for (i in 1:length(ob.list)) {
genes.use <- c(genes.use, head(rownames(ob.list[[i]]@hvg.info), 1000))
}
genes.use <- names(which(table(genes.use) > 1))
for (i in 1:length(ob.list)) {
genes.use <- genes.use[genes.use %in% rownames(ob.list[[i]]@scale.data)]
}
# 使用上述程序鉴定到的高变异基因，对4个数据集进行CCA分析，并计算前15个CC维度
gcdata.CCA <- RunMultiCCA(ob.list, genes.use = genes.use, num.ccs = 15)
# 查看细胞在CCA前两个维度空间中的分布，观察有无批次效应，如图4所示
# 可以看出与图2相比，批次效应得到一定程度的缓解
DimPlot(gcdata.CCA, reduction.use = "cca", group.by = "tech", pt.size = 0.5)
图4. CCA降维后细胞分布.png
# 查看在前几个CC维度上重要性较大的部分基因，这些基因的表达可能在细胞类别的区分上权重较大，如图5所示
DimHeatmap(object = gcdata.CCA, reduction.type = "cca", cells.use = 500, dim.use = 1:9, do.balanced = TRUE)
图5. 前9个CC维度上基因得分.png
# 使用MetageneBicorPlot功能，选择下游分析使用的CC维度，如图6所示
# 可选择前10-12个CC维度进行下游分析
MetageneBicorPlot(gcdata.CCA, grouping.var = "tech", dims.eval = 1:15)
图6. MetageneBicorPlot.png
# 将4个数据集的CCA维度空间进行比对，产生新的降维空间cca.aligned
# 在新的空间对细胞进行聚类分析
# 使用AlignSubspace，同时声明分组变量与比对维度
gcdata.CCA <- AlignSubspace(gcdata.CCA, grouping.var = "tech", dims.align = 1:12)
# 查看比对之后的维度分数分布，如图7所示
VlnPlot(gcdata.CCA, features.plot = c("ACC1","ACC2"),, group.by = "tech")
图7. CCA比对后维度分数分布.png
# 与未比对之前的CC维度分数分布进行比较，如图8所示
# 图7、图8比较，可见经过CCA比对之后，细胞分布收到批次效应影响的程度更小
VlnPlot(gcdata.CCA, features.plot = c("CC1","CC2"),, group.by = "tech")
图8. CCA比对前维度分数分布.png
# 在比对后的ACC空间进行聚类分析与tSNE降维分析
gcdata.CCA <- FindClusters(gcdata.CCA, reduction.type = "cca.aligned", dims.use = 1:12,
save.SNN = T, random.seed = 100)
gcdata.CCA <- RunTSNE(gcdata.CCA, reduction.use = "cca.aligned", dims.use = 1:12,
do.fast = TRUE, seed.use = 1)
# 在tSNE降维空间查看细胞聚类结果与批次分布，如图9所示
# 与图3比较，可见批次效应得到很大程度的缓解，各个批次的细胞呈均匀分布
p1 <- DimPlot(gcdata.CCA, reduction.use = "tsne", dim.1 = 1, dim.2 = 2,
group.by = "res.0.8", do.return = T)
p2 <- DimPlot(gcdata.CCA, reduction.use = "tsne", dim.1 = 1, dim.2 = 2,
group.by = "tech", do.return = T)
plot_grid(p1, p2)
图9. 去除批次效应后细胞聚类结果.png
# 计算矫正批次效应后的ARI，结果显示为0.04356477，基本接近与0
# 且与未矫正批次效应的ARI（0.3430069）相比，有明显降低，说明矫正后聚类结果基本不受批次效应影响
ari <- dplyr::select(gcdata.CCA@meta.data, tech, res.0.8)
ari$tech <- plyr::mapvalues(ari$tech,
from = c("celseq", "celseq2", "fluidigmc1", "smartseq2"),
to = c(0, 1, 2, 3))
mclust::adjustedRandIndex(as.numeric(ari$tech), as.numeric(ari$res.0.8))
sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.2
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] Seurat_2.3.4 Matrix_1.2-17 cowplot_0.9.4 ggplot2_3.1.1
loaded via a namespace (and not attached):
[1] tsne_0.1-3 segmented_0.5-4.0 nlme_3.1-139
[4] bitops_1.0-6 RColorBrewer_1.1-2 httr_1.4.0
[7] prabclus_2.3-1 tools_3.6.0 backports_1.1.4
[10] irlba_2.3.3 R6_2.4.0 rpart_4.1-15
[13] KernSmooth_2.23-15 Hmisc_4.2-0 lazyeval_0.2.2
[16] colorspace_1.4-1 nnet_7.3-12 npsurv_0.4-0
[19] withr_2.1.2 tidyselect_0.2.5 gridExtra_2.3
[22] compiler_3.6.0 htmlTable_1.13.1 hdf5r_0.0.0-9000
[25] labeling_0.3 diptest_0.75-7 caTools_1.17.1.2
[28] scales_1.1.0 checkmate_1.9.3 lmtest_0.9-37
[31] DEoptimR_1.0-8 mvtnorm_1.0-10 robustbase_0.93-5
[34] ggridges_0.5.1 pbapply_1.4-0 dtw_1.20-1
[37] proxy_0.4-23 stringr_1.4.0 digest_0.6.19
[40] mixtools_1.1.0 foreign_0.8-71 R.utils_2.8.0
[43] base64enc_0.1-3 pkgconfig_2.0.2 htmltools_0.3.6
[46] bibtex_0.4.2 htmlwidgets_1.3 rlang_0.4.5
[49] rstudioapi_0.10 farver_2.0.3 jsonlite_1.6
[52] zoo_1.8-6 ica_1.0-2 mclust_5.4.3
[55] gtools_3.8.1 acepack_1.4.1 dplyr_0.8.1
[58] R.oo_1.22.0 magrittr_1.5 modeltools_0.2-22
[61] Formula_1.2-3 lars_1.2 Rcpp_1.0.1
[64] munsell_0.5.0 reticulate_1.12 ape_5.3
[67] lifecycle_0.1.0 R.methodsS3_1.7.1 stringi_1.4.3
[70] gbRd_0.4-11 MASS_7.3-51.4 flexmix_2.3-15
[73] gplots_3.0.1.1 Rtsne_0.15 plyr_1.8.4
[76] grid_3.6.0 parallel_3.6.0 gdata_2.18.0
[79] crayon_1.3.4 doSNOW_1.0.16 lattice_0.20-38
[82] splines_3.6.0 SDMTools_1.1-221.1 knitr_1.23
[85] pillar_1.4.1 igraph_1.2.4.1 fpc_2.2-1
[88] reshape2_1.4.3 codetools_0.2-16 stats4_3.6.0
[91] glue_1.3.1 lsei_1.2-0 metap_1.1
[94] latticeExtra_0.6-28 data.table_1.12.2 png_0.1-7
[97] Rdpack_0.11-0 foreach_1.4.4 tidyr_0.8.3
[100] gtable_0.3.0 RANN_2.6.1 purrr_0.3.2
[103] kernlab_0.9-27 assertthat_0.2.1 xfun_0.7
[106] class_7.3-15 survival_2.44-1.1 tibble_2.1.2
[109] snow_0.4-3 iterators_1.0.10 cluster_2.0.8
[112] fitdistrplus_1.0-14 ROCR_1.0-7
主要参考资料：
