在單細胞分析當(dāng)中瞭吃，經(jīng)常會遇到整合分析的問題，即去除多樣本數(shù)據(jù)之間的批次效應(yīng)（batch effect）杜跷，那么什么是批次效應(yīng)呢傍念？簡而言之，批次效應(yīng)就是由于不同時間葛闷、不同實驗人員憋槐、不同儀器等因素造成的實驗性誤差区赵，而非本身的生物學(xué)差異泻肯。如果我們不去除批次效應(yīng)业踢，那么這些差異就會和本身的生物學(xué)差異相混淆陕悬。但是隨著測序成本的降低，單細胞測序已經(jīng)“深入尋常百姓家”搞糕，所以在追求大數(shù)據(jù)量的同時稿湿，肯定會伴隨著batch effect的產(chǎn)生集峦，自然batch effect的去除就成為單細胞數(shù)據(jù)分析的重要技能旷赖。2020年發(fā)表在Genome Biology上的一篇文章系統(tǒng)性總結(jié)了目前的batch effect去除方法顺又。
今天給大家分享幾種目前使用比較廣泛的單細胞數(shù)據(jù)整合分析的方法。本次演示所使用的示例數(shù)據(jù)如有需要等孵，可在留言區(qū)留言獲取。
merge()
首先是直接使用merge()函數(shù)對兩個單細胞數(shù)據(jù)進行直接整合蹂空，這時我們需要準備的輸入文件為一個由需要去除batch effect的Seurat對象組成的列表俯萌，那么如何實現(xiàn)呢？
library(Seurat)
library(ggplot2)
library(data.table)
setwd('GSE129139-mouse-T_cells')
samples=list.files('GSE129139_RAW/')
注意上枕，我們這里的數(shù)據(jù)是怎么存放的咐熙，我們在GSE129139_RAW/這個文件夾下面存放著我們需要去除batch effect的樣品數(shù)據(jù)，一個樣品辨萍，一個文件夾棋恼，每個文件夾里面是什么就不用說了吧！
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()輸入文件的文件夾锈玉，并依次創(chuàng)建了Seurat對象爪飘，存放在一個名為sceList的列表中。
然后我們利用merge()函數(shù)進行數(shù)據(jù)的整合：
sce.all = merge(sceList[[1]], y = sceList[[2]],
add.cell.ids = samples)
需要注意的是：（1）我們想把sample信息添加到cell barcode上拉背，只需要添加add.cell.ids參數(shù)即可师崎，這個參數(shù)賦給它一個向量；（2）上述的merge()默認只會簡單整合源數(shù)據(jù)（raw data）椅棺，如果你的Seurat對象是已經(jīng)經(jīng)過NormalizeData()的犁罩，可以直接添加merge.data = TRUE齐蔽，來merge標(biāo)準化后的數(shù)據(jù)。
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為了適應(yīng)大需求添加的新功能床估，錨點整合是從Seurat3開始上線的含滴，其原理在這里不贅述，放出原始論文鏈接Stuart, Butler, et al., Cell 2019 [Seurat V3]
同樣是需要由幾個Seurat對象組成的列表作為輸入丐巫，不同的是谈况，我們需要提前對數(shù)據(jù)進行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)
需要注意的是，從這里開始鞋吉，后面的數(shù)據(jù)分析請指定assay為integrated鸦做，否則你還在用原始的RNA assay進行分析，等于沒整合谓着。你可以通過以下命令更改默認assay泼诱，這樣就不用每次都進行聲明！
DefaultAssay(immune.combined) <- "integrated"
harmony
harmony單細胞數(shù)據(jù)整合方法于2019年發(fā)表在Nature Methods上赊锚，題為Fast, sensitive and accurate integration of single-cell data with Harmony治筒。harmony整合方法算得上是一種比較好的方法，目前應(yīng)用也是比較多的舷蒲，原理見文章耸袜，這里繼續(xù)展示具體流程：
#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還有參數(shù)lambda牲平，這個參數(shù)越小堤框，整合力度越大，一般按默認
#……
需要注意的是纵柿，如果你用harmony整合蜈抓，后續(xù)的下游分析，請指定 reduction = 'harmony' 昂儒，否則你的整合沒有意義沟使。
