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

  • 2020 年 3 月 30 日
  • 笔记

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

Seurat 学习笔记

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

我们一起来看一下。

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

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

步骤如下:

数据预处理

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

 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。

 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函数识别锚点。参数默认。

 reference.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")]   pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

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

 pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)

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

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

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

 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函数基于参考数据寻找细胞。

 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%的细胞被正确标记。

 pancreas.query$prediction.match <- pancreas.query$predicted.id == pancreas.query$celltype   table(pancreas.query$prediction.match)
 table(pancreas.query$predicted.id)
 VlnPlot(pancreas.query, c("REG1A", "PPY", "SST", "GHRL", "VWF", "SOX10"), group.by = "predicted.id")

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

未完待续…