前言
好久没写简书了，这两年多一直在做着单细胞的数据分析，有单细胞数据分析需求的朋友可以私信我。
目前单细胞转录组单个样本的成本不断降低，发表一篇高分文章的样本要求也越来越高，不可避免的我们会遇到多样本合并的问题。今天来介绍下目前单细胞转录组分析使用Seurat进行多样本合并的几种方法。
方法选择
当我们手里拿到多个样本数据时首先需要判断批次效应的大小，判断方法是我们先不考虑批次效应直接合并走单样本聚类分群流程。可以参考我之前的文章Seurat学习与使用（一）。
左图是三个样本直接合并按照单样本分析流程的聚类图，右图是三个样本合并后使用Harmony方法的聚类图。从两个图可以看到，左图样本间存在比较大的批次效应，相同细胞类型基本按照样本来分群。而右图可以看到相同细胞类型各样本很好的重叠在一起。目前scRNA数据校正批次效应的算法有很多种：主流的有CCA、MNN、CCA+MNN、RPCA、Harmony，当然还有其他的一些如Scanorama、LIGER、Conos等。下面我主要介绍CCA+MNN、MNN、RPCA、Harmony四种常用的方法。
CCA/CCA+MNN
CCA是早期Seurat（v2）整合数据的核心算法。文章发表在2018年的nature biotechnology，作者是Seurat的开发者Andrew Butler。同年在NBT杂志上Haghverdi等人开发了MNN的算法校正批次效应。之后在2019年Andrew等人在Cell发表文章，将CCA与MNN算法结合起来，并参考SNN算法的理念设计了“锚点”评分体系，使Seurat整合数据能力更强大更稳健。使用CCA方法时我们推荐CCA+MNN的方法，Seurat版本在v3及以上。CCA+MNN的使用方法可以参考我之前的文章Seurat学习与使用（一）中多个样本分析部分内容。
MNN
MNN，即相互近邻Mutual nearest neighbors，由Haghverdi提出。开发的思想是对一个批次中的细胞找到它们在另一个批次中对应的最相邻细胞集合。目前SeuratWrappers包已经包含了可以使用MNN方法的RunFastMNN函数。运行操作也很简单，如下：
library(SeuratWrappers)
merge_data<-merge(x = sample1, y=c(sample2,sample3) , add.cell.ids = c('S1','S1','S3'))#合并单个样本数据
merge_data<- RunFastMNN(object.list = SplitObject(merge_data, split.by = 'orig.ident'))#去批次
merge_data<- FindNeighbors(merge_data, reduction = 'mnn',dims = 1:15)#注意使用的是mnn数据
merge_data<- FindClusters(merge_data, resolution = 0.5)
merge_data<- RunUMAP(merge_data, reduction = 'mnn', dims = 1:15)#注意使用的是mnn数据
DimPlot(merge_data,label = T,reduction = 'umap')#可视化
RPCA
RPCA，即reciprocal PCA，是目前Seurat主推的一种去批次整合的方法。基于RPCA的整合运行得更快，同时与CCA方法相比，也是一种更保守的方法，不同生物状态的细胞在整合后不太可能“对齐”。在实际项目分析过程中我也发现RPCA的效果是要优于CCA，大家可以对比试试。详细使用方法可以点击reciprocal PCA链接。
Harmony
与其他整合方法相比，Harmony运行时比较快速，内存低，整合效果没CCA那么暴力，对稀有细胞的敏感性也比较好。同时在亚群细分时可以提取子集后再运行Harmony这样也比CCA的方法要方便。
Harmony算法概述如上图所示：用颜色表示数据集，用形状表示不同的细胞类型。首先，Harmony 应用PCA主成分分析将转录组范围的表达谱嵌入到降维空间中。然后我们应用一个迭代过程来移除数据集的特定效果。 (A) Harmony 将细胞概率性地分配给集群，最大限度地提高每个集群内数据集的多样性。 (B) Harmony 计算每个集群的所有数据集的全局质心，以及特定于数据集的质心。 (C) 在每个集群内，Harmony 根据质心计算每个数据集的校正因子。 (D) 最后，Harmony 使用基于 C 的细胞特定因子校正每个细胞。由于Harmony使用软聚类，每个单独的细胞可以通过与 A 中软聚类分配比例的多个因素的线性组合进行校正。 Harmony 重复步骤 A 到 D 直到收敛。集群分配和数据集之间的依赖性随着每一轮而减少。操作命令如下：
library(harmony)
merge_data<-merge(x = sample1, y=c(sample2,sample3) , add.cell.ids = c('S1','S1','S3'))#合并单个样本数据
merge_data<- merge_data %>% Seurat::NormalizeData(verbose = FALSE) %>% FindVariableFeatures(selection.method = 'vst', nfeatures = 2000) %>% ScaleData(verbose = FALSE) %>% RunPCA(verbose = FALSE)#标准化
merge_data<- RunHarmony(merge_data, 'orig.ident')#去批次
merge_data<- FindNeighbors(merge_data, reduction = 'harmony',dims = 1:15)#注意使用的是harmony数据
merge_data<- FindClusters(merge_data, resolution = 0.5)
merge_data<- RunUMAP(merge_data, reduction = 'harmony', dims = 1:15)#注意使用的是harmony数据
DimPlot(merge_data,label = T,reduction = 'umap')#可视化
总结
每种去批次的方法各有优劣，需要根据自己样本的实际情况选定。如果自己样本批次效应非常大（例如不同平台、不同组织），建议选用CCA的方法，如果批次效应较小（同一平台、相同组织），可以选择RPCA或者Harmony的方法。当然，样本整合的最终效果还得看分出来的群是否清晰无重叠、群与群之间的marker很特异，这都需要花费大量时间反复调整整合参数，最终才能得到理想的分群结果。
