一、前言

左边的图简单地把多个单细胞的数据合并在一起，不考虑去除批次效应，样本之间有明显的分离现象。右边的图是使用算法校正批次效应，不同的样本基本融和在一起了。scRNA数据校正批次效应的算法有很多：MNN, CCA+MNN, Harmony, Scanorama, scMerge等，本文推荐发表在Cell上的CCA+MNN方法，通过Seurat包就可以实现。

1. Seurat数据整合功能简介

Seurat早期版本整合数据的核心算法是CCA，文章发表在2018年的nature biotechnology，作者是Seurat的开发者Andrew Butler。同年Haghverdi等人开发了MNN算法校正批次效应，文章也发表在了nature biotechnology。2019年Andrew等人将CCA与MNN算法结合起来，并参考SNN算法的理念设计了“锚点”评分体系，使Seurat整合数据更强大更稳健。它不仅可以校正实验的批次效应，还能跨平台整合数据，例如将10x单细胞数据、BD单细胞数据和SMART单细胞数据整合在一起；也能整合单细胞多组学数据，例如将单细胞ATAC、空间转录组与单细胞转录组数据整合在一起。本文只讨论多样本数据的合并与校正批次效应，多组学数据的整合以后专门写篇文章介绍。

2. Seurat整合流程与原理

（1）使用CCA分析将两个数据集降维到同一个低维空间，因为CCA降维之后的空间距离不是相似性而是相关性，所以相同类型与状态的细胞可以克服技术偏倚重叠在一起。CCA分析效果见下图：

左图使用PCA降维，细胞之间的距离体现的是转录特征相似性，批次效应引入的系统误差会使样本分离。右图使用CCA降维，细胞之间的距离体现的是转录特征相关性，因此同类型且同状态的细胞可以跨越技术差异重叠在一起。

（2）CCA降维之后细胞在低维空间有了可以度量的“距离”，MNN(mutual nearest neighbor)算法以此找到两个数据集之间互相“距离”最近的细胞，Seurat将这些相互最近邻细胞称为“锚点细胞”。我们用两个数据集A和B来说明锚点，假设：

• 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两个数据集中的锚点之一。实际数据中，两个数据集之间的锚点可能有几百上千个，如下图所示：

图中每条线段连接的都是相互最近邻细胞

（3）理想情况下相同类型和状态的细胞才能构成配对锚点细胞，但是异常的情况也会出现，如上图中query数据集中黑色的细胞团。它在reference数据集没有相同类型的细胞，但是它也找到了锚点配对细胞（红色连线）。Seurat会通过两步过滤这些不正确的锚点：

1. 在CCA低维空间找到的锚点，返回到基因表达数据构建的高维空间中验证，如果它们的转录特征相似性高则保留，否则过滤此锚点。

2. 检查锚点细胞所在数据集最邻近的30个细胞，查看它们重叠的锚点配对细胞的数量，重叠越多分值越高，代表锚点可靠性更高。原理见下图：

左边query数据集的一个锚点细胞能在reference数据集邻近区域找到多个配对锚点细胞，可以得到更高的锚点可靠性评分；右边一个锚点细胞只能在reference数据集邻近区域找到一个配对锚点细胞，锚点可靠性评分则较低。

4、经过层层过滤剩下的锚点细胞对，可以认为它们是相同类型和状态的细胞，它们之间的基因表达差异是技术偏倚引起的。Seurat计算它们的差异向量，然后用此向量校正这个锚点锚定的细胞子集的基因表达值。校正后的基因表达值即消除了技术偏倚，实现了两个单细胞数据集的整合。

深究技术细节的朋友可以参阅原文：Tim S, Andrew Butler, Paul Hoffman , et al. Comprehensive integration of single cell data[J].Cell,2019.

二、整合

1. 读入数据

GSE162631，4个胶质瘤样本，总计5万多个细胞。

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



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

2. 读取样本，创建 Seurat 对象

在读取样本的同时直接设置样本的筛选条件，如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. 对数据先进行标准化，并识别 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. 整合数据

Seurat提供了一组Integration方法来去除批次效应，这些方法首先识别处于匹配生物学状态（anchors）的跨数据集细胞对，然后基于这些anchors校正数据集之间的批次效应

首先提取用来进行Integration的基因，然后找到anchors，基于anchors进行批次效应矫正

使用FindIntegrationAnchors函数识别锚点。参数默认。

然后我们将这些锚点传递给IntegrateData函数，该函数返回一个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 是两个样品经过批次效应矫正后合并的Seurat 对象，对这个对象进行分群分析

然后我们可以使用这个新的表达矩阵进行下游分析和可视化。

5. 下游分析

按照单样本的流程继续分析，包括进行标准化，运行PCA，并使用UMAP可视化结果。

# 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

可以看到经过批次矫正后，4个样品的大部分细胞在UMAP上都是重叠的。

6. 注释

测试样品经过批次矫正后，相同类型的细胞都聚到了一起，这样更有利于对细胞类型进行鉴定。实际分析时，我们可以对自己的样品同时做这两种分析，根据结果来判断是否需要进行批次矫正。

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

参考：
https://www.jianshu.com/p/9d798e95fb65
https://www.bilibili.com/read/cv15514134
https://www.jianshu.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