单细胞转录组整合分析——seurat包

2020-03-30 11:52:34 浏览数 (1)

Seurat是一个分析转录组数据的R包,我们之前的推文对其进行过描述:

Seurat 学习笔记

该包于去年新推出了整合功能。文章19年6月份发表于cell杂志,原文题目为:Comprehensive Integration of Single-Cell Data 被引量超过300次

我们一起来看一下。

该方法的目的是识别不同数据集中存在的共享细胞状态,即使它们是从不同的个体、实验条件、技术甚至物种中收集来的。

重点是找到不同数据集中的锚点anchors,这些“锚点”然后用于协调数据集,或将信息从一个数据集传输到另一个数据集。

步骤如下:

数据预处理

作者把单细胞数据放在了SeuratData等一系列包中,如果你的网速不行,可以直接到网页下载数据。

代码语言:javascript复制
 library(Seurat)
 #devtools::install_github('satijalab/seurat-data')
 library(SeuratData)
 #InstallData("panc8")
 #data("panc8")
 load('panc8.SeuratData/data/panc8.rda')
 
 #To construct a reference, we will identify ‘anchors’ between the individual datasets.
 #首先,将组合的数据分成列表,每个数据集是单独的元素
 pancreas.list <- SplitObject(panc8, split.by = "tech")
 pancreas.list <- pancreas.list[c("celseq", "celseq2", "fluidigmc1", "smartseq2")]

对数据先进行标准化,并识别variable feature。

代码语言:javascript复制
 for (i in 1:length(pancreas.list)) {
   pancreas.list[[i]] <- NormalizeData(pancreas.list[[i]], verbose = FALSE)
   pancreas.list[[i]] <- FindVariableFeatures(pancreas.list[[i]], selection.method = "vst",
                                              nfeatures = 2000, verbose = FALSE)
 }

整合3个胰岛细胞数据集

整合三个数据集作为参考,并使用FindIntegrationAnchors函数识别锚点。参数默认。

代码语言:javascript复制
 reference.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")]
 pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

然后我们将这些锚点传递给IntegrateData函数,该函数返回一个Seurat对象。

代码语言:javascript复制
 pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)

现在我们得到了seurat对象——一个整合后的表达矩阵pancreas.integrated。

然后我们可以使用这个新的表达矩阵进行下游分析和可视化。

包括进行标准化,运行PCA,并使用UMAP可视化结果。

代码语言:javascript复制
 library(ggplot2)
 library(cowplot)
 # switch to integrated assay. The variable features of this assay are automatically
 # set during IntegrateData
 DefaultAssay(pancreas.integrated) <- "integrated"
 
 # Run the standard workflow for visualization and clustering
 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)
 
 p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "tech")
 p2 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype", label = TRUE,
               repel = TRUE)   NoLegend()
 plot_grid(p1, p2)

左图按照技术聚类,右图按照细胞类型聚类。

使用参考数据集进行细胞类型分类

找到锚点之后,我们使用TransferData函数基于参考数据寻找细胞。

代码语言:javascript复制
 pancreas.query <- pancreas.list[["fluidigmc1"]]
 pancreas.anchors <- FindTransferAnchors(reference = pancreas.integrated, query = pancreas.query,
     dims = 1:30)
 predictions <- TransferData(anchorset = pancreas.anchors, refdata = pancreas.integrated$celltype,
     dims = 1:30)
 pancreas.query <- AddMetaData(pancreas.query, metadata = predictions)

因为我们有来自完整整合分析的原始标签注释,所以我们可以评估我们预测的细胞类型注释与完整参考的匹配程度。在这个例子中,我们发现在细胞类型分类上有很高的一致性,超过97%的细胞被正确标记。

代码语言:javascript复制
 pancreas.query$prediction.match <- pancreas.query$predicted.id == pancreas.query$celltype
 table(pancreas.query$prediction.match)
代码语言:javascript复制
 table(pancreas.query$predicted.id)
代码语言:javascript复制
 VlnPlot(pancreas.query, c("REG1A", "PPY", "SST", "GHRL", "VWF", "SOX10"), group.by = "predicted.id")

可以看到这几个基因在水平表达量的高低。

未完待续...

0 人点赞