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

一、前言

左邊的圖簡單地把多個單細胞的數(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單細胞數(shù)據(jù)、BD單細胞數(shù)據(jù)和SMART單細胞數(shù)據(jù)整合在一起；也能整合單細胞多組學(xué)數(shù)據(jù)，例如將單細胞ATAC、空間轉(zhuǎn)錄組與單細胞轉(zhuǎn)錄組數(shù)據(jù)整合在一起。本文只討論多樣本數(shù)據(jù)的合并與校正批次效應(yīng)，多組學(xué)數(shù)據(jù)的整合以后專門寫篇文章介紹。

2. Seurat整合流程與原理

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

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

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

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

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

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

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

  1. 在CCA低維空間找到的錨點，返回到基因表達數(shù)據(jù)構(gòu)建的高維空間中驗證，如果它們的轉(zhuǎn)錄特征相似性高則保留，否則過濾此錨點。

  2. 檢查錨點細胞所在數(shù)據(jù)集最鄰近的30個細胞，查看它們重疊的錨點配對細胞的數(shù)量，重疊越多分值越高，代表錨點可靠性更高。原理見下圖：

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

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

深究技術(shù)細節(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萬多個細胞。

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ù)先進行標(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ù)集細胞對，然后基于這些anchors校正數(shù)據(jù)集之間的批次效應(yīng)


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

使用FindIntegrationAnchors函數(shù)識別錨點。參數(shù)默認。

然后我們將這些錨點傳遞給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 對象，對這個對象進行分群分析

然后我們可以使用這個新的表達矩陣進行下游分析和可視化。

5. 下游分析

按照單樣本的流程繼續(xù)分析，包括進行標(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個樣品的大部分細胞在UMAP上都是重疊的。

6. 注釋

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

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.itdecent.cn/p/9d798e95fb65
https://www.bilibili.com/read/cv15514134
http://www.itdecent.cn/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

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容