上次介绍了scRNA-seq的常规流程，但是往往一个数据集的数据丰富度并不能满足我们日常分析的需要，那么我们肯定会选择多数据集联合分析来论证我们的结论，这样势必会面临着数据集融合的烦恼，这篇晨曦scRNA-seq笔记，将为你扫清多单细胞数据集联合分析的一切障碍~
首次揭秘！不做实验也能发10+SCI，CNS级别空间转录组套路全解析（附超详细代码！）
#准备工作，这里我们的示例数据需要下载并安装SeuratData包library(Seurat)library(SeuratData)library(patchwork)# install dataset（对网速要求比较严格）InstallData("ifnb")# load datasetLoadData("ifnb")#官网的示例代码给的是下载这个数据集，但是查看大小为300MB左右，很大，所以建议网速不是很给力的可以尝试使用panc8这个数据集，只有100MB左右，后面的分析也将使用panc8这个示例数据
晨曦解读
这一步其实特别考验网速，我这里也是捣鼓了半天，换了很多时间点才下载成功，所以如果下载示例数据失败也不要气馁，先理清思路即可~
#上面是官网提供的示例代码，下面是根据电脑情况设置的代码#安装并加载SeuratData包devtools::install_github('satijalab/seurat-data')library(SeuratData)library(Seurat)# 下载安装SeuratData包收集的特定数据集InstallData("panc8")# 加载数据集library(panc8.SeuratData)#查看一下示例数据data("panc8")#panc8#An object of class Seurat #34363 features across 14890 samples within 1 assay #Active assay: RNA (34363 features)
晨曦解读
既然这里我们是多数据集，所以我们就要把数据集进行分割，分割成多个，好在这个数据集本身也是存在着分组观测可以允许我们进行分割（这里我们采用meta信息中不用测序技术对Seurat对象进行分割）
#然后我们来具体看一下我们的数据view(panc8)
晨曦解读
这里面我们重点关注tech变量，因为它是我们分类的依据，根据查阅相关数据集信息，我们知道这个数据集是通过四种不同测序技术，包括CelSeq、CelSeq2、Fluidigm以及Smart-seq生成的人类胰岛细胞数据集，当然这里的技术我们也不是很清楚，但是我们最起码有一个最基本的认识，就是如果技术、平台、时间、操作手段不同的情况下，数据与数据之间是会存在批次效应的，所以我们学习如何合并数据集其实就是学习如何消除批次效应~
#进行分割数据集pancreas.list <- SplitObject(panc8, split.by = "tech")#然后我们看一下数据集view(pancreas.list)
晨曦解读
这样一个看似完整的数据集，就被我们按照技术的不同，也可以理解为平台的不同划分出来了5个数据集
#选出四个数据集来作为我们后续合并数据集的数据pancreas.list <- pancreas.list[c("celseq", "celseq2", "smartseq2", "fluidigmc1")]#这里因为对变量进行了重复赋值，所以以新的赋值变量为准#然后我们再来看一下我们现在的数据集列表view(pancreas.list)
晨曦解读
有可能有的同学看到这里会有疑惑，我们为什么要把数据集切分呢？其实这里是为了模拟我们工作中的场景，在我们平时测试去除批次效应方法的时候，这中切分数据集的方法也是一种很常用的用来准备输入数据的技巧~
#接下来对每个数据集进行标准的预处理，因为已经是Seurat对象，所以就是走常规的标准化以及寻找高可变基因pancreas.list <- lapply(X = pancreas.list, FUN = function(x){ x <- NormalizeData(x) x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)})
晨曦解读
下一步我们就要把这四个数据集进行合并，也是我们这里介绍的第一种方法CCA法
先介绍一下背景知识
CCA——这种合并数据集的方法的原理是，假设你要合并两个数据集，你找到两个数据集内相同的一个基因，这个基因需要满足稳定表达，然后以这个基因作为基准，来整合这两个数据集，是不是觉得很抽象？
那么我们现在来举个例子，有两个班级，我们把它称为A班和B班，其中A班是运动班，B班是学习班，有一天这两个班级需要整合成一个班级，这时候我们把这个班级的所有同学一起放在一个教室里，虽然表面上是一个班级，但是彼此之间的隔阂仍旧存在，也就是批次效应，这时候该如何消除隔阂呢？
这时候因为已经是一个班级了，我们找一个学习和体育都还差不多的同学来设定班级活动，这个同学肯定就会设计出学习和运动的活动数目都差不多，这样久而久之班级就融合了起来，这就是CCA
然后再来介绍一下后面会用到的概念
锚点——由一组共同的分子特征定义的两个细胞（每个数据集一个），将对应关系表示锚定，共同分子称为锚点
参考——具有批次效应的数据集整合起来后的大数据集
查询——还没有整合前单个的数据集
转化学习——生成在参考数据集基础上的模型，可以把这个模型在投射到查询上，以此来进行查询数据集的注释
#注意CCA的基础是寻找锚点，这个方法随着数据集的增加，时间是成倍的增加的，所以这里为了方便演示，我们选择三个数据集来进行操作reference.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")]pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)#这里数据集的维数可以选择默认，也可以自己调节，一般来说#这里的中间的数据列表不需要多再意，这是一个中间过程，然后把上一步计算出来的锚点信息带入到下面的函数，并得到整合数据集pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)
晨曦解读
如果大家确实认真阅读并实践过官方教程，这里大家可能会有疑问，请仔细思考一下下面代码的区别
#代码一pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)#代码二pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, anchor.features = features)
晨曦解读
相信大家已经看到区别了，这时候肯定会想知道代码二中的参数是什么意思，当然最快的学习方法自然是看帮助文档，但是这里我就直接公布答案了，代码一是我们自己设定维度，代码二则是通过调取前面默认的维度通过一系列计算给你确定维度，这里建议就自己设置就好，能更直观点~
最后我们得到了整合后的数据集，让我们来看一下整合后数据集的样子
View(pancreas.integrated)
晨曦解读
请大家仔细看我上面用红箭头标注的地方，再回想一下我上个笔记讲解过Seurat对象的特点，是不是有一种恍然大悟的感觉呢？
#数据集整合完后，我们自然是要进行可视化，来看一看结果#首先我们要先设定一下，我们接下来要使用整合后的数据，为什么要指定，请看上面的数据格式DefaultAssay(pancreas.integrated) <- "integrated"#然后进行标准的后续分析pancreas.integrated <- ScaleData(pancreas.integrated, verbose = FALSE)pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 30, verbose = FALSE)pancreas.integrated <- RunUMAP(pancreas.integrated, reduction = "pca", dims = 1:30)pancreas.integrated <- FindNeighbors(pancreas.integrated, reduction = "pca", dims = 1:30)pancreas.integrated <- FindClusters(pancreas.integrated, resolution = 0.5)
晨曦解读
基础要牢固哦~
#下面可视化，看看结果p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "tech")p2 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend()p1 + p2#NoLegend()参数的意思是去掉图例

晨曦解读
可以很清楚的看到不同tech的数据被很好的融合在了一起并没有独立的分开，说明我们确实去除了不同处理条件和技术下数据集的批次效应
#我们把不同处理条件放在并成一张图来看DimPlot(pancreas.integrated, reduction = "umap", split.by = "tech")[![2suyjI.png](https://z3.ax1x.com/2021/06/08/2suyjI.png
晨曦解读
到这一步，其实我们数据集整合就完成了，但是对于多数据集联合可不仅仅是光整合就够了哦，难道你不想探索一下，即使平台、技术、刺激方式等外在因素不同，仍然有些基因的表达是稳定的，也就是所谓的保守细胞markers，这些保守的marker可不可以为我们不管是后续的分析还是注释提供思路呢？
答案显然是可以的，这类基因不论外界的条件如何变化，它们的表达都是恒定的，我们可以从这些基因中选择表达强烈的基因来作为我们后续细胞亚群注释的marker基因
#首先我们这里使用原始数据集DefaultAssay(pancreas.integrated) <- "RNA"one.markers <- FindConservedMarkers(pancreas.integrated, ident.1 = 1, grouping.var = "tech", verbose = FALSE)
晨曦解读
这里解释一下参数，ident.1参数想必大家应该很熟悉了，就是指代我们这个函数作用哪个细胞亚群，grouping.var参数则是指代外部变化的条件，这里因为我们前面是按照不同技术划分的分组所以这里要保持和上面一样
这时候可能会有同学问，老师如果是不同外界刺激条件呢，那么可以把代码更改一下，如下所示
one.markers <-FindConservedMarkers(pancreas.integrated, ident.1 = 1,grouping.var = "stim", verbose = FALSE)
晨曦解读
然后我们来看一下运算后的结果~
晨曦解读
大家可以从这里的变量中看到规律，基本上这些变量大家都是知道什么意思的（不知道？快去看晨曦scRNA-seq笔记（1））然后我这里说一下，如果有几个外界变量这个变量大重复就会多重复一次~
#然后我们来探索一下不同Gene在细胞亚群中的表达FeaturePlot(pancreas.integrated, features = c("MAFA", "ADCYAP1", "WSCD2", "CYYR1"), min.cutoff = "q9")
晨曦解读
显然上面的基因都是我随便选的，所以看表达量不是特别的高，一般来说选择表达量大于3的两个或三个标记基因就可以作为我们后续细胞注释的marker基因了（这里并不是细胞注释的新方法，还是基于差异表达找寻marker基因的思路）
#上面我们已经探讨了刺激下表达恒定的基因，接下来我们就要来探讨不同刺激下，表达出现差异的基因#首先加载一些R包并进行一些准备工作library(ggplot2)library(cowplot)theme_set(theme_cowplot())#然后一下我们开始进行数据的整理#这里我们往往选择的是对照组和实验组的细胞亚群来进行分析#我们通过计算实验组AND对照组中细胞亚群中Gene的平均表达量并绘制散点图来寻找散点图中视觉异常值的基因#因为这里我们并没有进行细胞的注释，所以我们用编号来指代我们想要分析的细胞亚群one.cells <- subset(pancreas.integrated, idents = "1")Idents(one.cells) <- "tech"#这里填写的是具有分组效力的变量avg.one.cells <- as.data.frame(log1p(AverageExpression(one.cells, verbose = FALSE)$RNA))#这一步其实就是计算出我们这个细胞亚群的平均表达量，可以作为一个比较的基准avg.one.cells$gene <- rownames(avg.one.cells)
晨曦解读
这里就解释一下为什么使用log1p函数？
在进行数据预处理的时候，首先可以对偏度较大的数据用log1p函数进行转化，使其更加服从高斯分布，此步处理可能会使我们后续的分类结果得到一个更好的结果。同时使用log1p能避免复值（复值指一个自变量对应多个因变量）的问题
#在选择另一个细胞亚群，进行上述同样的操作two.cells <- subset(pancreas.integrated, idents = "2")Idents(two.cells) <- "tech"avg.two.cells <- as.data.frame(log1p(AverageExpression(two.cells, verbose = FALSE)$RNA))avg.two.cells$gene <- rownames(avg.two.cells)
晨曦解读
首先我们继续来观察一下数据，毕竟数据是我们可视化的基础，这里我们以avg.one.cells数据为例，大家可以很清楚的看到，针对不同处理方式，我们计算了其Gene的平均表达，这里如果替换成实验组和对照住可能会更加直观~
#然后我们进行简单的可视化#genes.to.label = c()#如果有想要标注在图上的基因可以把Gene symbol填写在这里p1 <- ggplot(avg.one.cells, aes(celseq, celseq2)) + geom_point() + ggtitle("one.cells")#p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)p1
晨曦解读
大家可以很清楚的看到尽管技术不同，但是基本上Gene的表达都是呈现出正相关也就是即技术之间存在差异但是并不是十分的相距甚远，这里如果换成实验组和对照组，则是说明在刺激条件下，细胞内大部分基因其实是不改变的，只有少部分基因发生改变以适应环境，这可能代表细胞保守途径中的一环~
#在这之前我们必须确保我们可以熟练的识别不同刺激条件下，细胞与细胞差异基因，因为这个毕竟是我们细胞注释的基础，然后我们才可以来探究不同条件下，相同细胞类型的差异基因#这里注意区分，我们在前面讲解的是同一个样本内不同细胞亚群之间的差异基因来作为我们的marker并进行细胞注释，这里则是两个样本（不同条件），然后相同细胞亚群之间的差异基因#首先在meta.data中创建一个变量，用来保存细胞类型和刺激信息，并将当前数据切换到该列pancreas.integrated$cell.tech <- paste(Idents(pancreas.integrated), pancreas.integrated$tech, sep = "_")pancreas.integrated$cell.type <- Idents(pancreas.integrated)Idents(pancreas.integrated) <- "cell.tech"#这句代码的意思其实就是把细胞亚型用我们第一句代码创建的形式来表示
晨曦解读
让我们来看一下我们在Seurat对象中新创建的变量
#随后我们通过FindMarkers()函数就可以找到不同外在条件下某亚群细胞之间不同基因one.interferon.response <- FindMarkers(pancreas.integrated, ident.1 = "0_celseq" , ident.2 = "1_celseq", verbose = FALSE)
晨曦解读
这样我们就得到了不用技术下相同细胞类型之间的差异表达，当然也有单纯只是可视化的方式也可以表明不同条件下相同细胞亚群中相同Gene的表达情况，留给同学们自己探究（示例代码在下面，大家可以尝试自己修改并且可以在评论区积极讨论哦~）
#这里把示例代码给大家，大家可以参照着来理解代码的含义FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,    cols = c("grey", "red"))plots <- VlnPlot(immune.combined, features = c("LYZ", "ISG15", "CXCL10"), split.by = "stim", group.by = "celltype", pt.size = 0, combine = FALSE)wrap_plots(plots = plots, ncol = 1)
晨曦解读
前面介绍了CCA的方法来联合多数据集进行分析，那么下面再介绍一种联合多数据集的方法——SCT，基本上这两种方法就是我们目前常见的两种联合多数据集的方法
#加载示例数据library(panc8.SeuratData)#拆分示例数据pancreas.list <- SplitObject(panc8, split.by = "tech")#对数据集进行标准化pancreas.list <- lapply(X = pancreas.list, FUN = SCTransform)#找寻数据集的特征基因（这一步其实很像scRNA-seq寻找高可变基因只不过原理是不同的，后面会解释）features <- SelectIntegrationFeatures(object.list = pancreas.list, nfeatures = 3000)#这一步是后面识别锚点之前需要运行的准备工作pancreas.list <- PrepSCTIntegration(object.list = pancreas.list, anchor.features = features)
晨曦解读
有意思的是SCTransform也可以取代我们常规的scRNA-seq三步骤即标化数据→高可变基因的筛选→标化数据（线性），这次将尽可能保证官网教程的完整性，后期将会列举四大方法联合多数据集分析流程~
#常规运行下面三句代码即可完成数据整合，因为SCT中已经包含ScaleData，所以不需要再次标准化数据，这里需要区别于CCA，CCA整合完数据后是需要走寻找到高可变基因后（不包括寻找高可变基因，因为CCA之前在整合数据集之前就已经进行过这一步骤了）的一系列常规操作的pancreas.anchors <- FindIntegrationAnchors(object.list = pancreas.list, normalization.method = "SCT", anchor.features = features)pancreas.combined.sct <- IntegrateData(anchorset = pancreas.anchors, normalization.method = "SCT")
晨曦解读
注意和CCA的步骤进行对比哦，SCT是不需要再次进行ScaleData的
#下面就开始对整合数据进行常规的操作pancreas.combined.sct <- RunPCA(pancreas.combined.sct, verbose = FALSE)pancreas.combined.sct <- RunUMAP(pancreas.combined.sct, reduction = "pca", dims = 1:30)p1 <- DimPlot(pancreas.combined.sct, reduction = "umap", group.by = "tech")p2 <- DimPlot(pancreas.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE, repel = TRUE)p1 + p2
晨曦解读
总体来说还不错，尽管因为不同技术造成了批次效应，但是通过SCT或者CCA的方法确实可以消除这样的批次效应，至于算法本身的差异大家感兴趣可以自行了解，下次将带大家再多解锁一些多数据集整合的方法，并且会出一次方法与方法之间的比较笔记，感兴趣的小伙伴可以关注下期的晨曦单细胞笔记
我是晨曦，我们下次再见~
回复“晨曦02”，即可获得今天的范例数据和代码哦~
—
END
—

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