单细胞数据整合分析——批次效应（batch effect）去除

在单细胞分析当中，经常会遇到整合分析的问题，即去除多样本数据之间的批次效应（batch effect），那么什么是批次效应呢？简而言之，批次效应就是由于不同时间、不同实验人员、不同仪器等因素造成的实验性误差，而非本身的生物学差异。如果我们不去除批次效应，那么这些差异就会和本身的生物学差异相混淆。但是随着测序成本的降低，单细胞测序已经“深入寻常百姓家”，所以在追求大数据量的同时，肯定会伴随着batch effect的产生，自然batch effect的去除就成为单细胞数据分析的重要技能。2020年发表在Genome Biology上的一篇文章系统性总结了目前的batch effect去除方法。

今天给大家分享几种目前使用比较广泛的单细胞数据整合分析的方法。本次演示所使用的示例数据如有需要，可在留言区留言获取。

merge()

首先是直接使用merge()函数对两个单细胞数据进行直接整合，这时我们需要准备的输入文件为一个由需要去除batch effect的Seurat对象组成的列表，那么如何实现呢？

library(Seurat)
library(ggplot2)
library(data.table)
setwd('GSE129139-mouse-T_cells')
samples=list.files('GSE129139_RAW/')

注意，我们这里的数据是怎么存放的，我们在GSE129139_RAW/这个文件夹下面存放着我们需要去除batch effect的样品数据，一个样品，一个文件夹，每个文件夹里面是什么就不用说了吧！

sceList = lapply(samples, function(pro){ 
  folder=file.path('GSE129139_RAW', pro) #file.path(a,b) means the path: a/b
  sce=CreateSeuratObject(counts = Read10X(folder), 
                         project = pro)
  return(sce)
})

上面的code实际上做了这样的一件事：按顺序读取了存放着三个Read10X()输入文件的文件夹，并依次创建了Seurat对象，存放在一个名为sceList的列表中。
然后我们利用merge()函数进行数据的整合：

sce.all = merge(sceList[[1]], y = sceList[[2]], 
                add.cell.ids = samples)

需要注意的是：（1）我们想把sample信息添加到cell barcode上，只需要添加add.cell.ids参数即可，这个参数赋给它一个向量；（2）上述的merge()默认只会简单整合源数据（raw data），如果你的Seurat对象是已经经过NormalizeData()的，可以直接添加merge.data = TRUE，来merge标准化后的数据。

By default, merge() will combine the Seurat objects based on the raw count matrices, erasing any previously normalized and scaled data matrices. If you want to merge the normalized data matrices as well as the raw count matrices, simply pass merge.data = TRUE. This should be done if the same normalization approach was applied to all objects.

Seurat锚点整合

这是Seurat为了适应大需求添加的新功能，锚点整合是从Seurat3开始上线的，其原理在这里不赘述，放出原始论文链接Stuart, Butler, et al., Cell 2019 [Seurat V3]
同样是需要由几个Seurat对象组成的列表作为输入，不同的是，我们需要提前对数据进行NormalizeData()和FindVariableFeatures()处理：

sceList <- lapply(sceList, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = sceList)
sce <- FindIntegrationAnchors(object.list = sceList, anchor.features = features)
# this command creates an 'integrated' data assay
sce <- IntegrateData(anchorset = sce)

需要注意的是，从这里开始，后面的数据分析请指定assay为integrated，否则你还在用原始的RNA assay进行分析，等于没整合。你可以通过以下命令更改默认assay，这样就不用每次都进行声明！

DefaultAssay(immune.combined) <- "integrated"

harmony

harmony单细胞数据整合方法于2019年发表在Nature Methods上，题为Fast, sensitive and accurate integration of single-cell data with Harmony。harmony整合方法算得上是一种比较好的方法，目前应用也是比较多的，原理见文章，这里继续展示具体流程：

#use the sce.all seurat object from the first part
sce.all <- NormalizeData(sce.all, normalization.method = "LogNormalize", scale.factor = 1e4) 
sce.all <- FindVariableFeatures(sce.all)
sce.all <- ScaleData(sce.all)
sce.all <- RunPCA(sce.all, features = VariableFeatures(object = sce.all))

library(harmony)
sce.all <- RunHarmony(sce.all, group.by.vars = "orig.ident") 
#group.by.vars表示按照什么来进行整合，亦即去掉什么方面的差异
#harmony还有参数lambda，这个参数越小，整合力度越大，一般按默认
#……

需要注意的是，如果你用harmony整合，后续的下游分析，请指定 reduction = 'harmony' ，否则你的整合没有意义。