scRNA多数据集联合的进阶探索
大家好，我是晨曦，不知道上周的推文大家学习的如何呢？
有学习上的困惑也欢迎大家在评论区或者酸谈社群里提问哦，晨曦只要看到都会回复的。那么今天我们依旧是打基础的环节，我们继续来把我们多数据集分析的最后一块缺口填补上，进行一下联合分析的总结，这也可能是我们基础阶段的倒数几篇推文了，也希望大家在打基础的过程中，可以更加意识到基础的重要性。
那么话不多说，我们这就开始吧~
首次揭秘！不做实验也能发10+SCI，CNS级别空间转录组套路全解析（附超详细代码！）
过关神助！99%审稿人必问，多数据集联合分析，你注意到这点了吗？
太猛了！万字长文单细胞分析全流程讲解，看完就能发文章！建议收藏！（附代码）
秀儿！10+生信分析最大的难点在这里！30多种方法怎么选？今天帮你解决！
图好看易上手！没有比它更适合小白入手的单细胞分析了！老实讲，这操作很sao！
毕业救星！这个R包在高分文章常见，实用！好学！
#准备工作library(Seurat)library(SeuratData)library(patchwork)#加载数据scRNA <- LoadData("ifnb")#查看数据head(scRNA)# orig.ident nCount_RNA nFeature_RNA stim seurat_annotations#AAACATACATTTCC.1 IMMUNE_CTRL 3017 877 CTRL CD14 Mono#AAACATACCAGAAA.1 IMMUNE_CTRL 2481 713 CTRL CD14 Mono#AAACATACCTCGCT.1 IMMUNE_CTRL 3420 850 CTRL CD14 Mono#AAACATACCTGGTA.1 IMMUNE_CTRL 3156 1109 CTRL pDC#AAACATACGATGAA.1 IMMUNE_CTRL 1868 634 CTRL CD4 Memory T#AAACATACGGCATT.1 IMMUNE_CTRL 1581 557 CTRL CD14 Mono#AAACATACTGCGTA.1 IMMUNE_CTRL 2747 980 CTRL T activated#AAACATACTGCTGA.1 IMMUNE_CTRL 1341 581 CTRL CD4 Naive T#AAACATTGAGTGTC.1 IMMUNE_CTRL 2155 880 CTRL CD8 T#AAACATTGCTTCGC.1 IMMUNE_CTRL 2536 669 CTRL CD14 Mono#切割数据ifnb.list <- SplitObject(ifnb, split.by = "stim")#对每个数据集进行标准的操作，包括标准化和寻找高可变Geneifnb.list <- lapply(X = ifnb.list, FUN = function(x) { x <- NormalizeData(x) x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)})#至于上面分割数据集后为什么要进行这类操作#首先官网在分割数据集后执行了这类操作这是第一#第二认为如果对于数据集进行改动后，为了平衡后续操作，使用标准化可以平衡数据与数据之间的改动#这里我们首先先对数据集进行合并，使用merge函数scRNA2 <- merge(ifnb.list[["CTRL"]],ifnb.list[["STIM"]])#对这个数据集走标准流程，这里可以知道，其实我们进行完数据集整合后，也是走标准流程scRNA1 <- NormalizeData(scRNA1)scRNA1 <- FindVariableFeatures(scRNA1, selection.method = "vst")scRNA1 <- ScaleData(scRNA1, features = VariableFeatures(scRNA1))scRNA1 <- RunPCA(scRNA1, features = VariableFeatures(scRNA1))plot1 <- DimPlot(scRNA1, reduction = "pca", group.by="orig.ident")plot2 <- ElbowPlot(scRNA1, ndims=30, reduction="pca") plotc <- plot1+plot2
晨曦解读
其实从PC来看就已经产生了明显的批次效应，但是我们还是应该以最后的聚类图为准~
#实际情况画tsne或者UMAP一个就可以，但是这里绘制全面，然后集中展示#设置降维数pc.num=1:30#聚类标准操作scRNA2 <- FindNeighbors(scRNA2, dims = pc.num) scRNA2 <- FindClusters(scRNA2, resolution = 0.5)#tSNEscRNA2 = RunTSNE(scRNA2, dims = pc.num)plot1 = DimPlot(scRNA1, reduction = "tsne", label=T) plot2 = DimPlot(scRNA1, reduction = "tsne", group.by='orig.ident')plot3 <- plot1+plot2plot3
晨曦解读
很明显的批次效应
#UMAP#UMAPscRNA2 <- RunUMAP(scRNA2, dims = pc.num)#group_by_clusterplot4 = DimPlot(scRNA2, reduction = "umap", label=T) #group_by_sampleplot5 = DimPlot(scRNA2, reduction = "umap", group.by='orig.ident')#combinateplot6 <- plot4+plot5plot6ggsave("UMAP_cluster.png", plot = plotc, width = 10, height = 5)
晨曦解读
很明显的批次效应~
#总览可视化plota <- plot2+plot5+ plot_layout(guides = 'collect')plota
晨曦解读
可以看到如果只是单纯的合并，凭借着Seurat对于scRNA-seq的处理，不足以消除数据与数据之间的批次效应，所以这里我们需要额外的算法来消除数据与数据之间的批次效应
接下来我们将介绍几种批次效应的方法，并且着重在运行速度和结果上来区分这些算法的差异，之所以不从算法本身来区分，首先是因为我本人并不是计算机领域出身，对于算法的研究并不是十分的擅长，而且我认为绝大多数同学其实考虑上面我说的这两方面即可~
方法一：CCA
#准备工作library(Seurat)library(SeuratData)library(patchwork)ifnb.list <- SplitObject(ifnb, split.by = "stim")#分别对两个数据集进行标准的处理流程ifnb.list <- lapply(X = ifnb.list, FUN = function(x) { x <- NormalizeData(x) x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)})
晨曦解读
上面的流程都是标准化的流程，没有什么需要讲解的地方
#跑一下标准流程（这个就是标准的流程，你只需要套入自己的数据即可）Sys.time()ifnb.anchors <- FindIntegrationAnchors(object.list = ifnb.list, dims = 1:30)ifnb.integrated <- IntegrateData(anchorset = ifnb.anchors, dims = 1:30)Sys.time()
晨曦解读
请大家仔细看我用箭头标注的位置，用时5min，时间上来说其实没有比较也不能说是快还是慢，但是这个方法大家直观看就知道是特别消耗内存的，因为这才是两个数据集联合，而且这两个数据集本身也不是很大，就找到了16393个锚点，所以这种方法我们可以简单得出一个结论：内存消耗较大，用时肯定也是相对来说比较长~
接下来来看一下合并数据集后的结果
#标准后续流程DefaultAssay(ifnb.integrated) <- "integrated"ifnb.integrated <- ScaleData(ifnb.integrated, verbose = FALSE)ifnb.integrated <- RunPCA(ifnb.integrated, npcs = 30, verbose = FALSE)ifnb.integrated <- RunUMAP(ifnb.integrated, reduction = "pca", dims = 1:30)ifnb.integrated <- FindNeighbors(ifnb.integrated, reduction = "pca", dims = 1:30)ifnb.integrated <- FindClusters(ifnb.integrated, resolution = 0.5)#可视化p1 <- DimPlot(ifnb.integrated, reduction = "umap", group.by = "stim")p2 <- DimPlot(ifnb.integrated, reduction = "umap", group.by = "seurat_annotations", label = TRUE, repel = TRUE) + NoLegend()p1 + p2
晨曦解读
批次效应被很好的去除，证明我们这种算法是有效的
同时，我们后续的分析就可以基于ifnb.integrated数据来进行展开，上面可以看到，我们是选择这个数据集来进行展开的
方法二：SCT
#这里我们依旧使用我们上面所用到的数据集，来保证我们的公平性#对数据集进行标准化Sys.time()ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)#找寻数据集的特征基因（这一步其实很像scRNA-seq寻找高可变基因只不过原理是不同的，后面会解释）features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)#这一步是后面识别锚点之前需要运行的准备工作ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)#ifnb.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT", anchor.features = features)ifnb.combined <- IntegrateData(anchorset = ifnb.anchors, normalization.method = "SCT")Sys.time()
晨曦解读
SCT内置的计算较多，所以图片都截不全，我这里只是对比一下时间
> Sys.time()——Start[1] "2021-07-01 20:09:48 CST"> Sys.time()——End[1] "2021-07-01 20:21:41 CST"
晨曦解读
运算时间上来说，用时12min，是明显比CCA的用时要长的，在不考虑结果的情况下，显然这种方法我们是不太建议选择的，但是作为结果，我这里有一点我自己的看法，本身我们本不是生信出身，对于各种算法，我们只需要会用就可以，既然CCA和SCT，或者是后面的算法能够被我们知道，能够被我们使用，证明这些算法肯定有可取之处，所以如果一再纠结于用哪个，倒不如全用一遍，毕竟代码都是标准化的流程，然后我们选择可视化结果好的即可~
#标准的后续分析流程ifnb.combined <- RunPCA(ifnb.combined, verbose = FALSE)ifnb.combined <- RunUMAP(ifnb.combined, reduction = "pca", dims = 1:30)p1 <- DimPlot(ifnb.combined, reduction = "umap", group.by = "stim")p2 <- DimPlot(ifnb.combined, reduction = "umap", group.by = "seurat_annotations", label = TRUE, repel = TRUE)p1 + p2
晨曦解读
很明显SCT的批次效应的去除要比我们CCA的感觉上要好上一点，具体为什么？
对比CCA和SCT的批次效应去除情况，SCT会上不同数据集之间融合的更加彻底，都是运行时间也相对较长，但是都是对的方法可以自行选择~
方法三：harmony
晨曦解读
这里参考网上教程得知，大多数人建议如果要进行跨物种整合或者不同数据类型整合的时候，我们可以使用CCA或者SCT，但是如果单纯的单细胞数据整合，我们可以考虑harmony
#Harmony接受的数据是PCA后的数据，而且代码也比较简单，所以我们这里直接列出来#而且Harmony也不用多个数据集分别操作，只需要融合在一起即可library(devtools)install_github("immunogenomics/harmony")scRNA_harmony <- merge(ifnb.list[["CTRL"]],ifnb.list[["STIM"]])#跑一下标准流程，这里涉及到管道符号可以加载tidyverse包，这里往下涉及到两种方法来获得我们harmony的输入数据，这里都进行一下演示#方法一library(tidyverse)ifnb.list <- SplitObject(ifnb, split.by = "stim")ifnb_harmony <- merge(ifnb.list[["CTRL"]],ifnb.list[["STIM"]])scRNA_harmony <- NormalizeData(ifnb_harmony) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)#因为我们知道SCT可以取代标准化-寻找高可变基因-线性标准化，这三个步骤，所以我们完全可以分成两步，一步SCT，再加上一步PCA即可获得Harmony的输入数据#方法二ifnb_harmony <- SCTransform(ifnb_harmony)ifnb_harmony <- RunPCA(ifnb_harmony，verbose=FALSE)#上面两种方法都可以达到我们的目的
晨曦解读
上面的流程说的已经很详细了，就是告诉我们如何准备Harmony的输入数据，但是这里我遇见了一个报错，也是在这里和大家分享一下，报错解决的流程
#扩展ifnb_harmony <- NormalizeData(ifnb_harmony) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)#Error: SCT assay is comprised of multiple SCT models. To change the variable features, please set manually with VariableFeatures<-
晨曦解读
首先解决报错最直接方法就是把这个报错信息复制，然后粘贴到任何一个浏览器的搜索框中，然后进行搜索，但是我搜索得到的结果如下
很遗憾上面的方法并没有帮助我解决报错，那么接下来，我就从报错信息入手来解决报错
报错信息提示我们在进行这一系列操作的时候，有重复的地方，那么就可以推断我们的输入数据有问题，我们就开始逐一运行函数来进行矫正
最后我发现，报错的原因来自FindVariableFeatures()函数，然后开始查看函数的帮助文档，因为这里我并没有使用任何参数，所以参数问题可以忽略，然后函数的使用环境也是正确的，那么就只有一种情况，就是这个函数的输入数据有问题
但是都是标准的Seurat对象怎么会出现问题呢？
最后经过仔细比对每一步信息，我发现了答案
#因为在前一步骤我们运行了这句函数ifnb.list <- SplitObject(ifnb, split.by = "stim")ifnb.list <- lapply(X = ifnb.list, FUN = function(x) { x <- NormalizeData(x) x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)})#这个本来是我们为了创建两个独立数据集而做的准备，都是这样显然是不可以的，因为这样我们显然是对每一个独立数据集进行了高可变基因的筛选，这个和我们后续的步骤是有重复的，所以我这里把这句函数注释掉以后，报错信息消失#这个看似不经意的报错就花费了我1个小时的时间
晨曦解读
这里想分享的是，生信分析有的地方需要严肃有的地方需要宽松，要有灵活的思维~
#那么我们继续正题library(harmony)ifnb_harmony <- RunHarmony(ifnb_harmony, group.by.vars = "stim")#group.by.vars参数为分组，我们要把哪个分组变量进行整合#这里面还有一个参数：max.iter.harmony是设置迭代次数，这个不需要深究，默认即可#lambde参数为整合力度，一般取值在0.5到2之间，越小，整合力度越大#然后整合完后我们使用整合后的数据进行后续的操作ifnb_harmony <- RunUMAP(ifnb_harmony, reduction = "harmony", dims = 1:30)ifnb_harmony <- FindNeighbors(ifnb_harmony, reduction = "harmony", dims = 1:30) %>% FindClusters()plot1 = DimPlot(ifnb_harmony, reduction = "umap", label=T，group.by = "seurat_annotations") plot2 = DimPlot(ifnb_harmony, reduction = "umap", group.by='stim')plot1+plot2#保存当前环境（通用）saveRDS(ifnb_harmony, 'ifnb_harmony.rds')
晨曦解读
这里大家需要捋清楚一个细节，我们整合数据，影响的是聚类的结果，但是我们差异分析是基于原始的表达矩阵，所以不同数据集整合的方法不会影响我们后续差异分析的结果
当然对这里的理解我也不是十分的到位，后期会专门解析一篇综述来帮助大家更好的理解，文献名字我放在下面，方便大家课余时间进行阅读
晨曦解读
我们可以从图上很清除的看到，批次效应已经被去除了，这里说一下harmony的优点，也方便各位同学进行选择
（1）整合数据的同时对稀有细胞的敏感性依然很好；

（2）省内存；

（3）适合于更复杂的单细胞分析实验设计，可以比较来自不同供体，组织和技术平台的细胞。
经过简单的评测可以知道，harmony的整合效果可以和CCA媲美，且运行速度有明显的优势
基本上我们常规掌握这三种整合数据的方式就可以了，毕竟算法日新月异，掌握几个主流的即可
下面介绍一个技巧：如果对于多个数据集，如何进行指定数据集的整合？
#因为这里需要多个数据集，我就换了一个数据集，并且进行了切分pancreas.list <- SplitObject(panc8, split.by = "tech")
#然后我们只需要进行下面一步提取子集的操作，即可通过改变变量的方式获取到我们想要的数据集列表#然后接下来就是执行标准的数据集整合流程pancreas.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")] #当然我们也可以通过添加参数来让这个流程更快一点#因为我们知道如果单纯运行CCA相当于做了一个数数比较，就是假设有5个数据集，A、B、C、D、E，那么运行CCA就会出现A和B，A和C，A和D等等这样一系列比较，这样的运行如果内存不够，往往都会失败#所以我们采用reference参数来人为设置参考数据集，即如果我设置1，3，则这个函数会先把1和3数据集合并，并且设置为ref，然后后续就变成了ref和其他数据集的一一合并，即ref同数据集A合并后，在和数据集B合并这样#这样的时间会大幅度减少#方法一：pancreas.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")] pancreas.list <- lapply(X = pancreas.list, FUN = function(x) { x <- NormalizeData(x) x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)})scRNA.anchors <- FindIntegrationAnchors(object.list = pancreas.list2, dims = 1:30,reference = c(1,3))pancreas.integrated <- IntegrateData(anchorset = scRNA.anchors, dims = 1:30)pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)#方法二：pancreas.list2 <- pancreas.list[c("celseq", "celseq2", "smartseq2")]pancreas.list2 <- lapply(X=pancreas.list, FUN=function(x) SCTransform(x))scRNA.anchors <- FindIntegrationAnchors(object.list = pancreas.list2, dims = 1:30,reference = c(1,3))pancreas.integrated <- IntegrateData(anchorset = scRNA.anchors, dims = 1:30)pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)#上面两种方法的区别在于是否用SCT来代替我们常规的标准化流程，然后针对使用SCT标准化方法后，合并数据集选择CCA的时候就不要进行scaleData的步骤，但是如果是走方法一，则在进行后续分析的时候需要加上scaleData的步骤#接下来把数据集设置为integrated，然后就可以走标准流程DefaultAssay(pancreas.integrated) <- "integrated"
晨曦解读
这里介绍了三种方法来整合数据集，亚群聚类使用integrated，差异分析依旧使用RNA，如果感到不解可以DefaultAssay()看一下究竟使用了什么数据集，具体的细节感兴趣可以参考上面那篇综述，但是我觉得其实注意一下这个点就可以，因为代码基本都是标准流程，因为在进行差异比较的时候，我们需要的输入数据是一个Seurat对象，所以具体细节的不明确，不太影响我们整个的分析流程~
晨曦解读
至于三种方法的选择就因人而异了，推荐harmony，CCA运行较慢但是比较经典，一点点拙见，也希望大家可以把自己的想法发在评论区，我们可以一起探讨，互相学习~
那么今天的推文到这里就结束啦，下一期估计就会迎来我们单细胞第一个高阶分析——monocle分析，也就是所谓的细胞轨迹分析，敬请期待吧~
我是晨曦，我们下次再见
主页回复“晨曦07”，即可获得今天的范例数据和代码哦~
—
END
—

撰文丨晨   曦
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文