代码分析 | 单细胞转录组数据整合详解

两种整合方法详解

NGS系列文章包括NGS基础、转录组分析（Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析（ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述）)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集）等内容。

单细胞转录组数据整合的方法多种多样，当然软件也是层出不穷，有的感觉矫正效应过强，导致不是一种细胞类型的细胞也会整合在一起，让人难以评判。前段时间有同学问我她有不同人相同肿瘤的样本，问我应该使用Merge还是使用CCA（单细胞分析Seurat使用相关的10个问题答疑精选！），我只能说实话我是不知道的，如果是我我都会试一试。由于肿瘤细胞的异质性过强，并且具有极大的样本差异性，如果使用CCA等进行整合，不知道会不会影响分析结果。并且在文章中我看的往往肿瘤样本会因为异质性导致无法重合，但其微环境如T、B等细胞还是可以有效的整合在一起，所以，保险起见，try a lot，select one!

本次接着上两节进行的芬兰CSC-IT科学中心主讲的生物信息课程(https://www.csc.fi/web/training/-/scrnaseq)视频，官网上还提供了练习素材以及详细代码，今天就来练习一下单细胞数据整合的过程。

在本教程中将探讨不同的整合多个单细胞RNA-seq数据集方法。我们将探索两种不同方法的校正整个数据集批次效应的效果并定量评估整合数据的质量。

数据集

在本教程中，我们将使用来自四种技术的3种不同的人类胰岛细胞数据集：CelSeq（GSE81076）、CelSeq2（GSE85241）、Fluidigm C1（GSE86469）和SMART-Seq2（E-MTAB-5061）。

原始数据矩阵和metadata下载链接：

https://www.dropbox.com/s/1zxbn92y5du9pu0/pancreas_v3_files.tar.gz?dl=1

加载所需R包

suppressMessages(require(Seurat))
suppressMessages(require(ggplot2))
suppressMessages(require(cowplot))
suppressMessages(require(scater))
suppressMessages(require(scran))
suppressMessages(require(BiocParallel))suppressMessages(require(BiocNeighbors))

Seurat (anchors and CCA)

我们将使用在文章Comprehensive Integration of Single Cell Data[1]中所提到的数据整合方法。

数据处理

加载表达式矩阵和metadata。metadata文件包含四个数据集中每个细胞所用技术平台和细胞类型注释。

pancreas.data <- readRDS(file = "session-integration_files/pancreas_expression_matrix.rds")metadata <- readRDS(file = "session-integration_files/pancreas_metadata.rds")

创建具有所有数据集的Seurat对象。

pancreas <- CreateSeuratObject(pancreas.data, meta.data = metadata)

在应用任何批次校正之前先看一下数据集。我们先执行标准预处理（log-normalization）并基于方差稳定化转换（“vst”）识别变量特征，接下来对集成数据进行归一化、运行PCA并使用UMAP可视化结果。集成数据集是按细胞类型而不是测序平台进行聚类。

# 标准化并查找可变基因
pancreas <- NormalizeData(pancreas, verbose = FALSE)
pancreas <- FindVariableFeatures(pancreas, selection.method = "vst", nfeatures = 2000, verbose = FALSE)

# 运行标准流程并进行可视化
pancreas <- ScaleData(pancreas, verbose = FALSE)
pancreas <- RunPCA(pancreas, npcs = 30, verbose = FALSE)
pancreas <- RunUMAP(pancreas, reduction = "pca", dims = 1:30)
p1 <- DimPlot(pancreas, reduction = "umap", group.by = "tech")
p2 <- DimPlot(pancreas, reduction = "umap", group.by = "celltype", label = TRUE, repel = TRUE) +
    NoLegend()plot_grid(p1, p2)

将合并的对象分成一个列表，每个数据集都作为一个元素。通过执行标准预处理（log-normalization）并基于方差稳定化转换（“vst”）分别为每个数据集查找变化的基因。

pancreas.list <- SplitObject(pancreas, split.by = "tech")

for (i in 1:length(pancreas.list)) {
    pancreas.list[[i]] <- NormalizeData(pancreas.list[[i]], verbose = FALSE)
    pancreas.list[[i]] <- FindVariableFeatures(pancreas.list[[i]], selection.method = "vst", nfeatures = 2000,
        verbose = FALSE)}

4个胰岛细胞数据集的整合

使用FindIntegrationAnchors函数来识别锚点（anchors），该函数的输入数据是Seurat对象的列表。

reference.list <- pancreas.list[c("celseq", "celseq2", "smartseq2", "fluidigmc1")]pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

## Computing 2000 integration features

## Scaling features for provided objects

## Finding all pairwise anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 3499 anchors

## Filtering anchors

##  Retained 2821 anchors

## Extracting within-dataset neighbors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 3515 anchors

## Filtering anchors

##  Retained 2701 anchors

## Extracting within-dataset neighbors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 6173 anchors

## Filtering anchors

##  Retained 4634 anchors

## Extracting within-dataset neighbors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2176 anchors

## Filtering anchors

##  Retained 1841 anchors

## Extracting within-dataset neighbors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2774 anchors

## Filtering anchors

##  Retained 2478 anchors

## Extracting within-dataset neighbors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2723 anchors

## Filtering anchors

##  Retained 2410 anchors
## Extracting within-dataset neighbors

然后将这些锚（anchors）传递给IntegrateData函数，该函数返回Seurat对象。

pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)

## Merging dataset 4 into 2

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

## Merging dataset 1 into 2 4

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

## Merging dataset 3 into 2 4 1

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights
## Integrating data

运行IntegrateData之后，Seurat对象将包含一个具有整合（或“批次校正”）表达矩阵的新Assay。请注意，原始值（未校正的值）仍存储在“RNA”分析的对象中，因此可以来回切换。

然后可以使用这个新的整合矩阵进行下游分析和可视化，我们在这里做了整合数据标准化、运行PCA并使用UMAP可视化结果。

# switch to integrated assay. The variable features of this assay are automatically set during
# IntegrateData
DefaultAssay(pancreas.integrated) <- "integrated"

# 运行标准流程并进行可视化
pancreas.integrated <- ScaleData(pancreas.integrated, verbose = FALSE)
pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 30, verbose = FALSE)
pancreas.integrated <- RunUMAP(pancreas.integrated, reduction = "pca", dims = 1:30)
p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "tech")
p2 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype", label = TRUE, repel = TRUE) +
    NoLegend()plot_grid(p1, p2)

Mutual Nearest Neighbor (MNN)

整合单细胞RNA-seq数据的另一种方法是使用相互最近邻（MNN）批次校正方法，可参考Haghverdi et al[2]。

首先可以直接从计数矩阵创建SingleCellExperiment（SCE）对象，也可以直接从Seurat转换为SCE。

celseq.data <- as.SingleCellExperiment(pancreas.list$celseq)
celseq2.data <- as.SingleCellExperiment(pancreas.list$celseq2)
fluidigmc1.data <- as.SingleCellExperiment(pancreas.list$fluidigmc1)smartseq2.data <- as.SingleCellExperiment(pancreas.list$smartseq2)

数据处理

查找共同基因并将每个数据集简化为那些共同基因：

keep_genes <- Reduce(intersect, list(rownames(celseq.data),rownames(celseq2.data),
                                     rownames(fluidigmc1.data),rownames(smartseq2.data)))
celseq.data <- celseq.data[match(keep_genes, rownames(celseq.data)), ]
celseq2.data <- celseq2.data[match(keep_genes, rownames(celseq2.data)), ]
fluidigmc1.data <- fluidigmc1.data[match(keep_genes, rownames(fluidigmc1.data)), ]smartseq2.data <- smartseq2.data[match(keep_genes, rownames(smartseq2.data)), ]

通过查找具有异常低的总计数或特征（基因）总数的异常值，使用calculateQCMetrics()计算QC来确定低质量细胞（对一篇单细胞RNA综述的评述：细胞和基因质控参数的选择）。

## celseq.data
celseq.data <- calculateQCMetrics(celseq.data)
low_lib_celseq.data <- isOutlier(celseq.data$log10_total_counts, type="lower", nmad=3)
low_genes_celseq.data <- isOutlier(celseq.data$log10_total_features_by_counts, type="lower", nmad=3)
celseq.data <- celseq.data[, !(low_lib_celseq.data | low_genes_celseq.data)]
## celseq2.data
celseq2.data <- calculateQCMetrics(celseq2.data)
low_lib_celseq2.data <- isOutlier(celseq2.data$log10_total_counts, type="lower", nmad=3)
low_genes_celseq2.data <- isOutlier(celseq2.data$log10_total_features_by_counts, type="lower", nmad=3)
celseq2.data <- celseq2.data[, !(low_lib_celseq2.data | low_genes_celseq2.data)]
## fluidigmc1.data
fluidigmc1.data <- calculateQCMetrics(fluidigmc1.data)
low_lib_fluidigmc1.data <- isOutlier(fluidigmc1.data$log10_total_counts, type="lower", nmad=3)
low_genes_fluidigmc1.data <- isOutlier(fluidigmc1.data$log10_total_features_by_counts, type="lower", nmad=3)
fluidigmc1.data <- fluidigmc1.data[, !(low_lib_fluidigmc1.data | low_genes_fluidigmc1.data)]
## smartseq2.data
smartseq2.data <- calculateQCMetrics(smartseq2.data)
low_lib_smartseq2.data <- isOutlier(smartseq2.data$log10_total_counts, type="lower", nmad=3)
low_genes_smartseq2.data <- isOutlier(smartseq2.data$log10_total_features_by_counts, type="lower", nmad=3)smartseq2.data <- smartseq2.data[, !(low_lib_smartseq2.data | low_genes_smartseq2.data)]

使用scran包的computeSumFactors()和normalize()函数计算大小（size）并对数据进行标准化：

# 计算尺寸因子（sizefactors）
celseq.data <- computeSumFactors(celseq.data)
celseq2.data <- computeSumFactors(celseq2.data)
fluidigmc1.data <- computeSumFactors(fluidigmc1.data)
smartseq2.data <- computeSumFactors(smartseq2.data)

# 标准化
celseq.data <- normalize(celseq.data)
celseq2.data <- normalize(celseq2.data)
fluidigmc1.data <- normalize(fluidigmc1.data)smartseq2.data <- normalize(smartseq2.data)

特征选择：我们使用TrendVar()和decomposeVar()函数来计算每个基因的变异（variance），并将其分为技术平台和生物两个部分的差异。

## celseq.data
fit_celseq.data <- trendVar(celseq.data, use.spikes=FALSE)
dec_celseq.data <- decomposeVar(celseq.data, fit_celseq.data)
dec_celseq.data$Symbol_TENx <- rowData(celseq.data)$Symbol_TENx
dec_celseq.data <- dec_celseq.data[order(dec_celseq.data$bio, decreasing = TRUE), ]

## celseq2.data
fit_celseq2.data <- trendVar(celseq2.data, use.spikes=FALSE)
dec_celseq2.data <- decomposeVar(celseq2.data, fit_celseq2.data)
dec_celseq2.data$Symbol_TENx <- rowData(celseq2.data)$Symbol_TENx
dec_celseq2.data <- dec_celseq2.data[order(dec_celseq2.data$bio, decreasing = TRUE), ]

## fluidigmc1.data
fit_fluidigmc1.data <- trendVar(fluidigmc1.data, use.spikes=FALSE)
dec_fluidigmc1.data <- decomposeVar(fluidigmc1.data, fit_fluidigmc1.data)
dec_fluidigmc1.data$Symbol_TENx <- rowData(fluidigmc1.data)$Symbol_TENx
dec_fluidigmc1.data <- dec_fluidigmc1.data[order(dec_fluidigmc1.data$bio, decreasing = TRUE), ]

## smartseq2.data
fit_smartseq2.data <- trendVar(smartseq2.data, use.spikes=FALSE)
dec_smartseq2.data <- decomposeVar(smartseq2.data, fit_smartseq2.data)
dec_smartseq2.data$Symbol_TENx <- rowData(smartseq2.data)$Symbol_TENx
dec_smartseq2.data <- dec_smartseq2.data[order(dec_smartseq2.data$bio, decreasing = TRUE), ]

#选择在所有数据集中共有的且信息最丰富的基因：
universe <- Reduce(intersect, list(rownames(dec_celseq.data),rownames(dec_celseq2.data),
                                   rownames(dec_fluidigmc1.data),rownames(dec_smartseq2.data)))
mean.bio <- (dec_celseq.data[universe,"bio"] + dec_celseq2.data[universe,"bio"] +
               dec_fluidigmc1.data[universe,"bio"] + dec_smartseq2.data[universe,"bio"])/4hvg_genes <- universe[mean.bio > 0]

将数据集合并到SingleCellExperiment中：

# 总原始counts的整合
counts_pancreas <- cbind(counts(celseq.data), counts(celseq2.data),
                         counts(fluidigmc1.data), counts(smartseq2.data))

# 总的标准化后的counts整合 (with multibatch normalization)
logcounts_pancreas <- cbind(logcounts(celseq.data), logcounts(celseq2.data),
                            logcounts(fluidigmc1.data), logcounts(smartseq2.data))

# 构建整合数据的sce对象
sce <- SingleCellExperiment(
    assays = list(counts = counts_pancreas, logcounts = logcounts_pancreas),
    rowData = rowData(celseq.data), # same as rowData(pbmc4k)
    colData = rbind(colData(celseq.data), colData(celseq2.data),
                    colData(fluidigmc1.data), colData(smartseq2.data))
)

# 将前面的hvg_genes存储到sce对象的metadata slot中:metadata(sce)$hvg_genes <- hvg_genes

用MNN处理批次效应之前先看一下这些datasets：

sce <- runPCA(sce,
              ncomponents = 20,
              feature_set = hvg_genes,
              method = "irlba")

names(reducedDims(sce)) <- "PCA_naive"

p1 <- plotReducedDim(sce, use_dimred = "PCA_naive", colour_by = "tech") +
    ggtitle("PCA Without batch correction")
p2 <- plotReducedDim(sce, use_dimred = "PCA_naive", colour_by = "celltype") +
    ggtitle("PCA Without batch correction")plot_grid(p1, p2)

使用MNN进行数据整合

scran软件包中的MNN方法利用一种新方法来调整批次效应-fastMNN()。

fastMNN()函数返回的是降维数据表示形式，该表示形式的使用与其他较低维度的表示形式（例如PCA）类似。

跑fastMNN()之前，我们需要先rescale每一个批次，来调整不同批次之间的测序深度。用scran包里的multiBatchNorm()函数对size factor进行调整后，重新计算log标准化的表达值以适应不同SingleCellExperiment对象的系统差异。之前的size factors仅能移除单个批次里细胞之间的bias。现在我们通过消除批次之间技术差异来提高校正的质量。

rescaled <- multiBatchNorm(celseq.data, celseq2.data, fluidigmc1.data, smartseq2.data)
celseq.data_rescaled <- rescaled[[1]]
celseq2.data_rescaled <- rescaled[[2]]
fluidigmc1.data_rescaled <- rescaled[[3]]smartseq2.data_rescaled <- rescaled[[4]]

跑fastMNN，把降维的MNN representation存在sce对象的reducedDimsslot里：

mnn_out <- fastMNN(celseq.data_rescaled,
                   celseq2.data_rescaled,
                   fluidigmc1.data_rescaled,
                   smartseq2.data_rescaled,
                   subset.row = metadata(sce)$hvg_genes,
                   k = 20, d = 50, approximate = TRUE,
                   # BPPARAM = BiocParallel::MulticoreParam(8),
                   BNPARAM = BiocNeighbors::AnnoyParam())
reducedDim(sce, "MNN") <- mnn_out$correct

注意：fastMNN（）不会生成批次校正的表达矩阵。因此，fastMNN（）的结果应仅被视为降维表示，适合直接绘图，如TSNE/ UMAP、聚类和依赖于此类结果的轨迹分析（NBT|45种单细胞轨迹推断方法比较，110个实际数据集和229个合成数据集）。

p1 <- plotReducedDim(sce, use_dimred = "MNN", colour_by = "tech") + ggtitle("MNN Ouput Reduced Dimensions")
p2 <- plotReducedDim(sce, use_dimred = "MNN", colour_by = "celltype") + ggtitle("MNN Ouput Reduced Dimensions")plot_grid(p1, p2)

Session info

sessionInfo()

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets## [8] methods   base

[1]：https://www.biorxiv.org/content/10.1101/460147v1[2]：https://www.nature.com/articles/nbt.4091

撰文：Tiger校对：生信宝典

你可能还想看

“harmony”整合不同平台的单细胞数据之旅

Cell 深度一套普遍适用于各类单细胞测序数据集的锚定整合方案

重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述）

如何使用Bioconductor进行单细胞分析？

如何火眼金睛鉴定那些单细胞转录组中的混杂因素

什么？你做的差异基因方法不合适？

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 158,847评论 4赞 362
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 67,208评论 1赞 292
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 108,587评论 0赞 243
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 43,942评论 0赞 205
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 52,332评论 3赞 287
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 40,587评论 1赞 218
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 31,853评论 2赞 312
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 30,568评论 0赞 198
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 34,273评论 1赞 242
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 30,542评论 2赞 246
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 32,033评论 1赞 260
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 28,373评论 2赞 253
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 33,031评论 3赞 236
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 26,073评论 0赞 8
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 26,830评论 0赞 195
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 35,628评论 2赞 274
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 35,537评论 2赞 269

代码分析 | 单细胞转录组数据整合详解

数据集

加载所需R包

Seurat (anchors and CCA)

数据处理

4个胰岛细胞数据集的整合

Mutual Nearest Neighbor (MNN)

数据处理

使用MNN进行数据整合

Session info

你可能还想看

推荐阅读更多精彩内容