一、前言

左邊的圖簡單地把多個單細(xì)胞的數(shù)據(jù)合并在一起趋惨，不考慮去除批次效應(yīng)，樣本之間有明顯的分離現(xiàn)象惦蚊。右邊的圖是使用算法校正批次效應(yīng)器虾，不同的樣本基本融和在一起了讯嫂。scRNA數(shù)據(jù)校正批次效應(yīng)的算法有很多：MNN, CCA+MNN, Harmony, Scanorama, scMerge等，本文推薦發(fā)表在Cell上的CCA+MNN方法兆沙，通過Seurat包就可以實現(xiàn)欧芽。

1. Seurat數(shù)據(jù)整合功能簡介

Seurat早期版本整合數(shù)據(jù)的核心算法是CCA，文章發(fā)表在2018年的nature biotechnology葛圃，作者是Seurat的開發(fā)者Andrew Butler千扔。同年Haghverdi等人開發(fā)了MNN算法校正批次效應(yīng)，文章也發(fā)表在了nature biotechnology库正。2019年Andrew等人將CCA與MNN算法結(jié)合起來昏鹃，并參考SNN算法的理念設(shè)計了“錨點”評分體系，使Seurat整合數(shù)據(jù)更強大更穩(wěn)健诀诊。它不僅可以校正實驗的批次效應(yīng)洞渤，還能跨平臺整合數(shù)據(jù)，例如將10x單細(xì)胞數(shù)據(jù)属瓣、BD單細(xì)胞數(shù)據(jù)和SMART單細(xì)胞數(shù)據(jù)整合在一起载迄；也能整合單細(xì)胞多組學(xué)數(shù)據(jù)，例如將單細(xì)胞ATAC抡蛙、空間轉(zhuǎn)錄組與單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)整合在一起护昧。本文只討論多樣本數(shù)據(jù)的合并與校正批次效應(yīng)，多組學(xué)數(shù)據(jù)的整合以后專門寫篇文章介紹粗截。

2. Seurat整合流程與原理

（1）使用CCA分析將兩個數(shù)據(jù)集降維到同一個低維空間惋耙，因為CCA降維之后的空間距離不是相似性而是相關(guān)性，所以相同類型與狀態(tài)的細(xì)胞可以克服技術(shù)偏倚重疊在一起熊昌。CCA分析效果見下圖：

左圖使用PCA降維绽榛，細(xì)胞之間的距離體現(xiàn)的是轉(zhuǎn)錄特征相似性，批次效應(yīng)引入的系統(tǒng)誤差會使樣本分離婿屹。右圖使用CCA降維灭美，細(xì)胞之間的距離體現(xiàn)的是轉(zhuǎn)錄特征相關(guān)性，因此同類型且同狀態(tài)的細(xì)胞可以跨越技術(shù)差異重疊在一起昂利。

（2）CCA降維之后細(xì)胞在低維空間有了可以度量的“距離”届腐，MNN(mutual nearest neighbor)算法以此找到兩個數(shù)據(jù)集之間互相“距離”最近的細(xì)胞，Seurat將這些相互最近鄰細(xì)胞稱為“錨點細(xì)胞”蜂奸。我們用兩個數(shù)據(jù)集A和B來說明錨點犁苏，假設(shè)：

A樣本中的細(xì)胞A3與B樣本中距離最近的細(xì)胞有3個（B1,B2,B3）
B樣本中的細(xì)胞B1與A樣本中距離最近的細(xì)胞有4個（A1,A2,A3,A4）
B樣本中的細(xì)胞B2與A樣本中距離最近的細(xì)胞有2個（A5,A6）
B樣本中的細(xì)胞B3與A樣本中距離最近的細(xì)胞有3個（A1,A2,A7）

那么A3與B1是相互最近鄰細(xì)胞，A3與B2扩所、B3不是相互最近鄰細(xì)胞围详，A3+B1就是A、B兩個數(shù)據(jù)集中的錨點之一碌奉。實際數(shù)據(jù)中短曾，兩個數(shù)據(jù)集之間的錨點可能有幾百上千個寒砖，如下圖所示：

圖中每條線段連接的都是相互最近鄰細(xì)胞

（3）理想情況下相同類型和狀態(tài)的細(xì)胞才能構(gòu)成配對錨點細(xì)胞，但是異常的情況也會出現(xiàn)嫉拐，如上圖中query數(shù)據(jù)集中黑色的細(xì)胞團哩都。它在reference數(shù)據(jù)集沒有相同類型的細(xì)胞，但是它也找到了錨點配對細(xì)胞（紅色連線）婉徘。Seurat會通過兩步過濾這些不正確的錨點：

在CCA低維空間找到的錨點漠嵌，返回到基因表達(dá)數(shù)據(jù)構(gòu)建的高維空間中驗證，如果它們的轉(zhuǎn)錄特征相似性高則保留盖呼，否則過濾此錨點儒鹿。
檢查錨點細(xì)胞所在數(shù)據(jù)集最鄰近的30個細(xì)胞，查看它們重疊的錨點配對細(xì)胞的數(shù)量几晤，重疊越多分值越高约炎，代表錨點可靠性更高。原理見下圖：

左邊query數(shù)據(jù)集的一個錨點細(xì)胞能在reference數(shù)據(jù)集鄰近區(qū)域找到多個配對錨點細(xì)胞蟹瘾，可以得到更高的錨點可靠性評分圾浅；右邊一個錨點細(xì)胞只能在reference數(shù)據(jù)集鄰近區(qū)域找到一個配對錨點細(xì)胞，錨點可靠性評分則較低憾朴。

4狸捕、經(jīng)過層層過濾剩下的錨點細(xì)胞對，可以認(rèn)為它們是相同類型和狀態(tài)的細(xì)胞众雷，它們之間的基因表達(dá)差異是技術(shù)偏倚引起的灸拍。Seurat計算它們的差異向量，然后用此向量校正這個錨點錨定的細(xì)胞子集的基因表達(dá)值砾省。校正后的基因表達(dá)值即消除了技術(shù)偏倚鸡岗，實現(xiàn)了兩個單細(xì)胞數(shù)據(jù)集的整合。

深究技術(shù)細(xì)節(jié)的朋友可以參閱原文：Tim S, Andrew Butler, Paul Hoffman , et al. Comprehensive integration of single cell data[J].Cell,2019.

二纯蛾、整合

1. 讀入數(shù)據(jù)

GSE162631纤房，4個膠質(zhì)瘤樣本，總計5萬多個細(xì)胞翻诉。

library(dplyr)
library(Seurat)
library(patchwork)



dirs = dir(pattern = "^R")
f = "dat.Rdata"
if(!file.exists(f)){
  scelist = list()
}

2. 讀取樣本，創(chuàng)建 Seurat 對象

在讀取樣本的同時直接設(shè)置樣本的篩選條件捌刮，如RNA碰煌，UMI，線粒體等绅作。

for(i in 1:length(dirs)){
  x = Read10X(data.dir = dirs[[i]])
  scelist[[i]] <- CreateSeuratObject(counts = x, project = paste0("R",i))
  scelist[[i]][["percent.mt"]] <- PercentageFeatureSet(scelist[[i]], pattern = "^MT-")
  scelist[[i]] <- subset(scelist[[i]], subset = percent.mt < 10)
}
names(scelist)  = paste0("R",1:4)
sum(sapply(scelist, function(x)ncol(x@assays$RNA@counts)))

3. 對數(shù)據(jù)先進(jìn)行標(biāo)準(zhǔn)化芦圾，并識別 variable feature

# normalize and identify variable features for each dataset independently
scelist <- lapply(X = scelist, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
  })

4. 整合數(shù)據(jù)

Seurat提供了一組Integration方法來去除批次效應(yīng)，這些方法首先識別處于匹配生物學(xué)狀態(tài)（anchors）的跨數(shù)據(jù)集細(xì)胞對俄认，然后基于這些anchors校正數(shù)據(jù)集之間的批次效應(yīng)

首先提取用來進(jìn)行Integration的基因个少，然后找到anchors洪乍，基于anchors進(jìn)行批次效應(yīng)矯正

使用FindIntegrationAnchors函數(shù)識別錨點。參數(shù)默認(rèn)夜焦。

然后我們將這些錨點傳遞給IntegrateData函數(shù)壳澳，該函數(shù)返回一個Seurat對象。

features <- SelectIntegrationFeatures(object.list = scelist)
immune.anchors <- FindIntegrationAnchors(object.list = scelist, anchor.features = features)
immune.combined <- IntegrateData(anchorset = immune.anchors)
DefaultAssay(immune.combined) <- "integrated"

immune.combined 是兩個樣品經(jīng)過批次效應(yīng)矯正后合并的Seurat 對象茫经，對這個對象進(jìn)行分群分析

然后我們可以使用這個新的表達(dá)矩陣進(jìn)行下游分析和可視化巷波。

5. 下游分析

按照單樣本的流程繼續(xù)分析，包括進(jìn)行標(biāo)準(zhǔn)化卸伞，運行PCA抹镊，并使用UMAP可視化結(jié)果。

# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)
save(immune.combined,file = f)


load(f)
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

可以看到經(jīng)過批次矯正后荤傲，4個樣品的大部分細(xì)胞在UMAP上都是重疊的垮耳。

6. 注釋

測試樣品經(jīng)過批次矯正后，相同類型的細(xì)胞都聚到了一起遂黍，這樣更有利于對細(xì)胞類型進(jìn)行鑒定终佛。實際分析時，我們可以對自己的樣品同時做這兩種分析妓湘，根據(jù)結(jié)果來判斷是否需要進(jìn)行批次矯正查蓉。

library(celldex)
library(SingleR)
#ref <- celldex::HumanPrimaryCellAtlasData()
ref <- get(load("../single_ref/ref_Human_all.RData"))
library(BiocParallel)
pred.scRNA <- SingleR(test = immune.combined@assays$integrated@data, 
                      ref = ref,
                      labels = ref$label.main, 
                      clusters = immune.combined@active.ident)
pred.scRNA$pruned.labels
##  [1] "Macrophage"        "Macrophage"        "Monocyte"         
##  [4] "Macrophage"        "Macrophage"        "Macrophage"       
##  [7] "Macrophage"        "Monocyte"          "Neutrophils"      
## [10] "Neutrophils"       "Endothelial_cells" "Monocyte"         
## [13] "Macrophage"        "Macrophage"        "Tissue_stem_cells"
## [16] "NK_cell"           "Monocyte"          "B_cell"
plotScoreHeatmap(pred.scRNA, clusters=pred.scRNA@rownames, fontsize.row = 9,show_colnames = T)

image

new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(immune.combined)
levels(immune.combined)
##  [1] "0"  "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
## [16] "15" "16" "17"
immune.combined <- RenameIdents(immune.combined,new.cluster.ids)
levels(immune.combined)
## [1] "Macrophage"        "Monocyte"          "Neutrophils"      
## [4] "Endothelial_cells" "Tissue_stem_cells" "NK_cell"          
## [7] "B_cell"
UMAPPlot(object = immune.combined, pt.size = 0.5, label = TRUE)

代碼主要來自：https://satijalab.org/seurat/articles/integration_introduction.html

參考：
http://www.reibang.com/p/9d798e95fb65
https://www.bilibili.com/read/cv15514134
http://www.reibang.com/p/d323d0291dd4
https://baijiahao.baidu.com/s?id=1707789923750679171&wfr=spider&for=pc
https://mp.weixin.qq.com/s/PA_xTrqVYiZCOscP43CgUw
https://mp.weixin.qq.com/s/FWvOEJSNb9epNZmxwFkhkg
https://www.genenergy.cn/news/156/2379.html

scRNA-Seq | 多樣本Seurat+CCA整合