兩種整合方法詳解
NGS系列文章包括NGS基礎、轉(zhuǎn)錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學設計實驗GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。
單細胞轉(zhuǎn)錄組數(shù)據(jù)整合的方法多種多樣，當然軟件也是層出不窮，有的感覺矯正效應過強，導致不是一種細胞類型的細胞也會整合在一起，讓人難以評判。前段時間有同學問我她有不同人相同腫瘤的樣本，問我應該使用Merge還是使用CCA（單細胞分析Seurat使用相關的10個問題答疑精選！），我只能說實話我是不知道的，如果是我我都會試一試。由于腫瘤細胞的異質(zhì)性過強，并且具有極大的樣本差異性，如果使用CCA等進行整合，不知道會不會影響分析結果。并且在文章中我看的往往腫瘤樣本會因為異質(zhì)性導致無法重合，但其微環(huán)境如T、B等細胞還是可以有效的整合在一起，所以，保險起見，try a lot，select one!
本次接著上兩節(jié)進行的芬蘭CSC-IT科學中心主講的生物信息課程(https://www.csc.fi/web/training/-/scrnaseq)視頻，官網(wǎng)上還提供了練習素材以及詳細代碼，今天就來練習一下單細胞數(shù)據(jù)整合的過程。
在本教程中將探討不同的整合多個單細胞RNA-seq數(shù)據(jù)集方法。我們將探索兩種不同方法的校正整個數(shù)據(jù)集批次效應的效果并定量評估整合數(shù)據(jù)的質(zhì)量。
數(shù)據(jù)集
在本教程中，我們將使用來自四種技術的3種不同的人類胰島細胞數(shù)據(jù)集：CelSeq（GSE81076）、CelSeq2（GSE85241）、Fluidigm C1（GSE86469）和SMART-Seq2（E-MTAB-5061）。
原始數(shù)據(jù)矩陣和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]中所提到的數(shù)據(jù)整合方法。
數(shù)據(jù)處理
加載表達式矩陣和metadata。metadata文件包含四個數(shù)據(jù)集中每個細胞所用技術平臺和細胞類型注釋。
pancreas.data <- readRDS(file = "session-integration_files/pancreas_expression_matrix.rds")metadata <- readRDS(file = "session-integration_files/pancreas_metadata.rds")創(chuàng)建具有所有數(shù)據(jù)集的Seurat對象。
pancreas <- CreateSeuratObject(pancreas.data, meta.data = metadata)在應用任何批次校正之前先看一下數(shù)據(jù)集。我們先執(zhí)行標準預處理（log-normalization）并基于方差穩(wěn)定化轉(zhuǎn)換（“vst”）識別變量特征，接下來對集成數(shù)據(jù)進行歸一化、運行PCA并使用UMAP可視化結果。集成數(shù)據(jù)集是按細胞類型而不是測序平臺進行聚類。
# 標準化并查找可變基因
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)將合并的對象分成一個列表，每個數(shù)據(jù)集都作為一個元素。通過執(zhí)行標準預處理（log-normalization）并基于方差穩(wěn)定化轉(zhuǎn)換（“vst”）分別為每個數(shù)據(jù)集查找變化的基因。
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個胰島細胞數(shù)據(jù)集的整合
使用FindIntegrationAnchors函數(shù)來識別錨點（anchors），該函數(shù)的輸入數(shù)據(jù)是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函數(shù)，該函數(shù)返回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”分析的對象中，因此可以來回切換。
然后可以使用這個新的整合矩陣進行下游分析和可視化，我們在這里做了整合數(shù)據(jù)標準化、運行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數(shù)據(jù)的另一種方法是使用相互最近鄰（MNN）批次校正方法，可參考Haghverdi et al[2]。
首先可以直接從計數(shù)矩陣創(chuàng)建SingleCellExperiment（SCE）對象，也可以直接從Seurat轉(zhuǎn)換為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)數(shù)據(jù)處理
查找共同基因并將每個數(shù)據(jù)集簡化為那些共同基因：
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)), ]通過查找具有異常低的總計數(shù)或特征（基因）總數(shù)的異常值，使用calculateQCMetrics()計算QC來確定低質(zhì)量細胞（對一篇單細胞RNA綜述的評述：細胞和基因質(zhì)控參數(shù)的選擇）。
## 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()函數(shù)計算大?。?code>size）并對數(shù)據(jù)進行標準化：
# 計算尺寸因子（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()函數(shù)來計算每個基因的變異（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), ]
#選擇在所有數(shù)據(jù)集中共有的且信息最豐富的基因：
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]將數(shù)據(jù)集合并到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))
# 構建整合數(shù)據(jù)的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進行數(shù)據(jù)整合
scran軟件包中的MNN方法利用一種新方法來調(diào)整批次效應-fastMNN()。
fastMNN()函數(shù)返回的是降維數(shù)據(jù)表示形式，該表示形式的使用與其他較低維度的表示形式（例如PCA）類似。
跑fastMNN()之前，我們需要先rescale每一個批次，來調(diào)整不同批次之間的測序深度。用scran包里的multiBatchNorm()函數(shù)對size factor進行調(diào)整后，重新計算log標準化的表達值以適應不同SingleCellExperiment對象的系統(tǒng)差異。之前的size factors僅能移除單個批次里細胞之間的bias?，F(xiàn)在我們通過消除批次之間技術差異來提高校正的質(zhì)量。
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個實際數(shù)據(jù)集和229個合成數(shù)據(jù)集）。
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校對：生信寶典
