在我们上一个视频中进行了单细胞转录组测序数据分析中 最重要 的单样本分析后，本次课程将进入 多样本数据整合 模块的学习。在各位 实际 数据分析的过程中，实验设计往往是多样本，多组别的，因此，因此多样本数据的 科学整合与管理 是十分有必要。这里我们会针对不同的 数据特征，考虑到 批次效应 、 计算机资源利用 特征教大家如何进行多样本 Seurat 对象或矩阵的整合与拆分。
不同实验批次、不同样本、不同实验条件、不同技术平台的单细胞数据整合是单细胞数据分析长期以来的重难点。让不同数据间相同的细胞类型进行高质量、具有置信度的匹配能够让后续的基因差异分析更加的科学。因此，Seurat v5 致力于移除数据间的条件差异而非生物学差异，较V4而言更是优化了 内存 使用与 计算速度 等方面。
具体流程与原来的多样本整合中演示步骤具有一定的差别，详细差别可参考：快速上手Seurat V5(订阅用户请在龙年订阅资料下载连接中找，无需再次付费)。1. 视频及测试文件 1.1 视频教程 1.2 测试文件 本手册下载时包含测试数据，若无测试数据，可在推送中获取。
推送链接 2. 代码 # 加载R包：library (Seurat) library (multtest) library (dplyr) library (ggplot2) library (patchwork) if ( ! require (SeuratData))devtools :: install_github ( 'satijalab/seurat-data' ) 2.1. 单纯的merge 2.1.1 准备用于拆分的数据集 你如果对拆分的过程或 ident 的管理不理解，例如你不理解为什么下文中有 group 变量，可参考以下两讲内容：Seurat中分类变量处理技巧 答读者问
（六）：Seurat中如何让细胞听你指挥 单细胞的多样本整合十分耗费计算资源，当各位实战时遇到内存不足的问题，可以考虑入一个具有性价比的服务器：足够支持你完成硕博生涯的生信环境: # 如果你的内存不是很充裕，可以进行抽样，降低计算量 # pbmc <- subset(pbmc, downsample = 50) ifnb <- readRDS ( './pbmcrenamed.rds' ) # 这其实是一个已经整合好的Seurat对象 class (ifnb) # [1] "Seurat" # attr(,"package") # [1] "Seurat" # 这是一个很老的数据，以Seurat V3版本生成的，V4、V5可以向下兼容 SeuratObject :: Version (ifnb) # [1] '3.2.2' # 其中包含了三个分组：unique (ifnb $ group) # [1] "C57" "AS1" "P3" # 那我们就按照分组将这个Seurat对象拆开，split的方式在Seurat V5中有所变化：# ifnb.list <- SplitObject(ifnb, split.by = "group") # 原来的方式 ifnb[[ "RNA" ]] <- split (ifnb[[ "RNA" ]], f = ifnb $ group) class (ifnb) # [1] "Seurat" # attr(,"package") # [1] "Seurat" # 分别取出这三个对象 C57 <- ifnb[,ifnb $ group == 'C57' ] C57 # An object of class Seurat # 22254 features across 309 samples within 2 assays # Active assay: RNA (20254 features, 0 variable features) # 2 layers present: counts.C57, data.C57 # 1 other assay present: integrated # 3 dimensional reductions calculated: pca, umap, tsne AS1 <- ifnb[,ifnb $ group == 'AS1' ] AS1 # An object of class Seurat # 22254 features across 287 samples within 2 assays # Active assay: RNA (20254 features, 0 variable features) # 2 layers present: counts.AS1, data.AS1 # 1 other assay present: integrated # 3 dimensional reductions calculated: pca, umap, tsne P3 <- ifnb[,ifnb $ group == 'P3' ] P3 # An object of class Seurat # 22254 features across 304 samples within 2 assays # Active assay: RNA (20254 features, 0 variable features) # 2 layers present: counts.P3, data.P3 # 1 other assay present: integrated # 3 dimensional reductions calculated: pca, umap, tsne # 注意，这里就等同于我们之前读入单个文件。
每一个对象相当于读取了一次一个样本的单细胞数据后生成的单样本矩阵。单样本矩阵读取后，参考我们前面介绍的单样本读取流程，接下来创建每个单样本seurta对象，即可继续接下来的merge合并对象流程。
2.1.2 merge # 简单利用merge整合这两个数据 pbmc <- merge (C57, y = c (AS1), add.cell.ids = c ( "C57" , "AS1" ), project = "ALL" ) # 这时我们就获得了一个新的Seurat对象 pbmc # An object of class Seurat # 22254 features across 596 samples within 2 assays # Active assay: RNA (20254 features, 0 variable features) # 4 layers present: counts.C57.1, data.C57.1, counts.AS1.2, data.AS1.2 # 1 other assay present: integrated # 可以发现这是barcode的名称发生了变化：head ( colnames (pbmc)) # [1] "C57_AAAGGGCTCGCAGAGA-1_1" "C57_AACACACAGCATCCCG-1_1" # [3] "C57_AACAGGGTCAACACCA-1_1" "C57_AACCCAAAGGCCTTGC-1_1" # [5] "C57_AACCCAAAGTGGACGT-1_1" "C57_AACCCAAGTCCGAAGA-1_1" unique ( sapply ( X = strsplit ( colnames (pbmc), split = "_" ), FUN = "[" , 1 )) # [1] "C57" "AS1" table (pbmc $ orig.ident) # # AS1 C57 # 287 309 # 以此类推，合并多个对象：merged_obj <- merge ( x = C57, y = c (AS1,P3), add.cell.ids = c ( "C57" , "AS1" , "P3" ), project = "MergedProject" ) merged_obj # An object of class Seurat # 22254 features across 900 samples within 2 assays # Active assay: RNA (20254 features, 0 variable features) # 6 layers present: counts.C57.1, data.C57.1, counts.AS1.2, data.AS1.2, counts.P3.3, data.P3.3 # 1 other assay present: integrated rm (merged_obj) # 接下来的分析流程可参考上一讲单样本分析 Seurat V5 中与 Seurat V4 不同的是，现在 list 可以直接参与 Seurat 中的计算函数，无需先合并各样本，也无需写循环采用遍历的形式。
# 标准分析流程：ifnb <- NormalizeData (ifnb) # 标准化 ifnb <- FindVariableFeatures (ifnb) # 计算高变基因 ifnb <- ScaleData (ifnb) # 缩放高变基因 ifnb <- RunPCA (ifnb) # PCA分析 ifnb <- FindNeighbors (ifnb, dims = 1 : 30 , reduction = "pca" ) # 计算最近临近图 ifnb <- FindClusters (ifnb, resolution = 0.4 , cluster.name = "unintegrated_clusters" ) # 分群 # Modularity Optimizer version 1.3.0 by Ludo Waltmand Nees Jan van Eck # # Number of nodes: 900 # Number of edges: 23290 # # Running Louvain algorithm... # Maximum modularity in 10 random starts: 0.9462 # Number of communities: 10 # Elapsed time: 0 seconds ifnb <- RunUMAP (ifnb, dims = 1 : 30 , reduction = "pca" , reduction.name = "umap.unintegrated" ) # 利用不去批次的数据进行降维 DimPlot (ifnb, reduction = "umap.unintegrated" , group.by = c ( "group" , "seurat_clusters" )) # 降维图展示批次效应 table ( Idents (ifnb),ifnb @ meta.data $ group) # # AS1 C57 P3 # 0 13 71 44 # 1 56 48 5 # 2 47 40 15 # 3 46 33 23 # 4 28 0 69 # 5 49 17 28 # 6 18 52 22 # 7 7 46 26 # 8 19 2 42 # 9 4 0 30 2.1.3 接下来我们尝试使用一个新的数据做多样本merge 数据下载地址：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE139107 library (Matrix) # read.delim直接读取压缩包,而且相较于read.table更为简洁，默认分隔符就是"\t" # 两个文件细胞数较多，读取时间约3min test1 <- read.delim ( "./GSE139107/GSE139107_MouseIRI_4hours.dge.txt.gz" ) test2 <- read.delim ( "./GSE139107/GSE139107_MouseIRI_12hours.dge.txt.gz" ) # 创建对象 obj4 <- CreateSeuratObject ( counts = test1) obj5 <- CreateSeuratObject ( counts = test2) obj <- merge (obj4, y = obj5, project = "ALL" ) head (obj @ meta.data) table (obj @ meta.data $ orig.ident) head ( colnames (obj4)) head ( colnames (obj5)) # 标准化 obj <- NormalizeData (obj) # 计算高变基因 obj <- FindVariableFeatures (obj) # 缩放高变基因 obj <- ScaleData (obj) # PCA分析 obj <- RunPCA (obj) # 计算最近临近图 obj <- FindNeighbors (obj, dims = 1 : 30 , reduction = "pca" ) # 分群 obj <- FindClusters (obj, resolution = 0.4 , cluster.name = "unintegrated_clusters" ) # 利用不去批次的数据进行降维 obj <- RunUMAP (obj, dims = 1 : 30 , reduction = "pca" , reduction.name = "umap.unintegrated" ) # saveRDS(obj,"iri.intergrate.havebatch.rds") # saveRDS(obj4,"iri.obj4.rds") # saveRDS(obj5,"iri.obj5.rds") obj4 <- readRDS ( "iri.obj4.rds" ) obj5 <- readRDS ( "iri.obj5.rds" ) obj <- readRDS ( "iri.intergrate.havebatch.rds" ) cols <- c ( " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " , " ) # 保存色板，后续会使用到这个色板 # save(cols,file = "color.Rdata") p1 <- DimPlot (obj, reduction = "umap.unintegrated" , group.by = "orig.ident" , cols = cols) # 降维图展示批次效应 # 指定颜色 p2 <- DimPlot (obj, reduction = "umap.unintegrated" , group.by = "seurat_clusters" , label = T, cols = cols ) p1 | p2 table ( Idents (obj),obj @ meta.data $ orig.ident) # # IRI12h1b1 IRI12h1b2 IRI12h2 IRI12h3 IRI4h1 IRI4h2 IRI4h3 # 0 16 7 42 2 1775 29 1758 # 1 576 45 469 442 873 201 711 # 2 634 586 425 492 205 65 431 # 3 573 560 503 484 18 211 50 # 4 1203 82 52 996 11 18 15 # 5 1147 29 32 1049 6 7 6 # 6 7 35 962 0 8 5 982 # 7 339 327 197 282 207 50 211 # 8 2 1560 7 1 0 0 4 # 9 0 1405 13 0 0 1 0 # 10 7 1 2 5 1239 9 54 # 11 248 287 239 177 163 62 101 # 12 209 184 175 178 230 63 236 # 13 7 5 1 4 1186 36 18 # 141 122 126 88 527 36 192 # 15 4 3 1173 0 1 0 8 # 16 229 169 167 164 133 49 252 # 17 6 9 11 2 414 11 377 # 18 4 0 1 0 0 792 0 # 19 200 196 65 129 5 3 15 # 20 1 0 0 2 0 297 0 # 21 40 27 37 32 52 25 42 可以看到的是，数据存在一定的批次效应，这也就是为什么我们需要在多样本整合的过程中使用下面几节的去批次算法。
2.2. anchor(CCA) 此方法为最耗费内存以及计算时间的方法，但去批次效果最为强劲，且能返回 integrated 去批次整合后获得的表达矩阵。
Canonical Correlation Analysis (CCA)去批次比较耗时且占内存，我们可以用以下方式开启并行计算：# 先设置并行计算，加速计算过程：makecore <- function (workcore,memory){ if ( ! require (Seurat)) install.packages ( 'Seurat' ) if ( ! require (future)) install.packages ( 'future' ) plan ( "multisession" , workers = workcore) options ( future.globals.maxSize= memory * 1024 * 1024 ** 2 ) } makecore ( 4 , 4 ) # 以下是Seurat V4的运行方式，这里不做运行 library (Seurat) library (tidyverse) # 定义一个函数帮助预处理数据 myfunction1 <- function (testA.seu){ testA.seu <- NormalizeData (testA.seu, normalization.method = "LogNormalize" , scale.factor = 10000 ) testA.seu <- FindVariableFeatures (testA.seu, selection.method = "vst" , nfeatures = 2000 ) return (testA.seu) } # 执行对每个数据的NormalizeData（标准化表达矩阵）与FindVariableFeatures（高变基因计算），这两部操作的过程可参考第三讲：相关链接见专题页 obj4_test <- myfunction1 (obj4) obj5_test <- myfunction1 (obj5) # Integration步骤，真正耗时耗电的步骤：testAB.anchors <- FindIntegrationAnchors ( object.list = list (obj4_test,obj5_test), dims = 1 : 20 ) # 寻找数据锚定点,约25min testAB.integrated <- IntegrateData ( anchorset = testAB.anchors, dims = 1 : 20 ) # 根据锚定点整合数据 # 需要注意的是：上面的整合步骤相对于**harmony**整合方法更耗费计算资源，尤其是较大的数据集（几万个细胞）非常消耗内存和时间，大约**9G**的数据在**32G**的内存下就已经无法运行；
当存在某一个`Seurat`对象细胞数很少（印象中200以下这样子），会报错，这时建议用第二种整合方法。
# 下游计算 # # 将后续计算用默认矩阵由"RNA"改为"integrated" DefaultAssay (testAB.integrated) <- "integrated" # 整合后的数据从scale开始运行，后续基本与单样本分析部分无异 testAB.integrated <- ScaleData (testAB.integrated, features = rownames (testAB.integrated)) testAB.integrated <- RunPCA (testAB.integrated, npcs = 50 , verbose = FALSE ) testAB.integrated <- FindNeighbors (testAB.integrated, dims = 1 : 30 ) testAB.integrated <- FindClusters (testAB.integrated, resolution = 0.4 ) testAB.integrated <- RunUMAP (testAB.integrated, dims = 1 : 30 ) testAB.integrated <- RunTSNE (testAB.integrated, dims = 1 : 30 ) DimPlot (testAB.integrated, reduction= "umap" ) DimPlot (testAB.integrated, reduction= "tsne" ) # saveRDS(testAB.integrated,"iri.intergrate.havebatch.v4CCA.rds") Seurat V5 中的运行则会简单许多：# 多个layer之间直接去批次整合：obj_test1 <- IntegrateLayers ( object = obj, method = CCAIntegration, # 以CCA的方式整合 orig.reduction = "pca" , new.reduction = "integrated.cca" , # 整合后的降维命名，收录在'obj@reductions$integrated.cca'中 verbose = FALSE ) # 可以看到过程非常简单，省去了原来寻找anchor的步骤 # 把两个Seurat对象的layer合并：obj_test1[[ "RNA" ]] <- JoinLayers (obj_test1[[ "RNA" ]]) obj_test1 <- FindNeighbors (obj_test1, reduction = "integrated.cca" , dims = 1 : 30 ) # 利用去批次后的数据计算最近临近图 obj_test1 <- FindClusters (obj_test1, resolution = 0.4 ) # 利用去批次后的数据分群 obj_test1 <- RunUMAP (obj_test1, dims = 1 : 30 , reduction = "integrated.cca" ) # 利用去批次后的数据降维 # 保存Seurat对象 # saveRDS(obj_test1,"iri.intergrate.havebatch.v5CCA.rds") obj_test1 <- readRDS ( "iri.intergrate.havebatch.v5CCA.rds" ) # 可视化：p3 <- DimPlot (obj_test1, reduction = "umap" , group.by = "orig.ident" , cols = cols) # 查看降维结果 p4 <- DimPlot (obj_test1, reduction = "umap" , group.by = "seurat_clusters" , cols = cols) p3 | p4 table ( Idents (obj_test1),obj_test1 @ meta.data $ orig.ident) # # IRI12h1b1 IRI12h1b2 IRI12h2 IRI12h3 IRI4h1 IRI4h2 IRI4h3 # 0 1050 136 63 973 2867 640 128 # 1 12 1038 1081 5 53 68 1647 # 2 1207 73 27 998 1215 299 63 # 3 7 1126 966 2 7 27 972 # 4 580 560 503 485 459 214 240 # 5 297 565 490 204 364 74 696 # 6 560 487 393 439 199 76 420 # 7 334 504 207 222 414 165 322 # 8 327 305 185 265 219 51 283 # 9 212 188 175 180 230 63 236 # 10 248 288 238 177 165 66 102 # 11 282 7 6 239 487 116 8 # 12 220 168 173 161 123 48 193 # 13 142 121 129 91 173 36 84 # 14 75 46 26 56 26 2 28 # 15 40 27 37 32 52 25 41 # 现在也不会生成新的'integrated'矩阵来干扰大家，assay中仍是只有'RNA'矩阵 names (obj_test1 @ assays) # [1] "RNA" 2.3. harmony 速度快、内存少 2.3.1 预处理 if ( ! require (harmony))devtools :: install_github ( "immunogenomics/harmony" ) # harmony前需要完成标准化、高变基因计算、scale、PCA等分析 obj_test2 <- merge (obj4, y = obj5, project = "ALL" ) obj_test2 <- obj_test2 %>% NormalizeData %>% FindVariableFeatures ( selection.method = "vst" , nfeatures = 2000 ) %>% ScaleData obj_test2 <- RunPCA (obj_test2, npcs = 50 , verbose = FALSE ) 2.3.2 harmony run 到 PCA 再进行 harmony，harmony 实际上相当于一种降维方法 obj_test2 <- obj_test2 %>% RunHarmony ( "orig.ident" , plot_convergence = TRUE ) # 以Harmony的结果进行后续的降维与分群 # UMAP及分群：obj_test2 <- obj_test2 %>% RunUMAP ( reduction = "harmony" , dims = 1 : 30 ) %>% FindNeighbors ( reduction = "harmony" , dims = 1 : 30 ) %>% FindClusters ( resolution = 0.4 ) %>% identity # TSNE, 运行时间比umap长 obj_test2 <- obj_test2 %>% RunTSNE ( reduction = "harmony" , dims = 1 : 30 ) # 保存seurat对象 # saveRDS(obj_test2,"iri.intergrate.harmony.rds") # 创建图片对象 obj_test2 <- readRDS ( "iri.intergrate.harmony.rds" ) p5 <- DimPlot (obj_test2, reduction = "umap" , group.by = "orig.ident" , pt.size= 0.5 , cols = cols) + theme ( axis.line = element_blank , axis.ticks = element_blank , axis.text = element_blank ) p6 <- DimPlot (obj_test2, reduction = "umap" , group.by = "ident" , pt.size= 0.5 , label = TRUE , repel = TRUE , cols = cols) + theme ( axis.line = element_blank , axis.ticks = element_blank , axis.text = element_blank ) # 看一下降维结果展示 p5 | p6