【原】Seurat4.0系列教程22：空间转录组的分析、可视化与整合

健明 2021-08-12

展开全文

Seurat4.0系列教程告一段落，但这决不是终点。这个系列教程只是给大家打开一扇窗，知道Seurat4.0有这些功能可用，少走弯路，起到一个抛砖引玉的作用，后续还要自己深入研究。不要像我当初入门单细胞之时，为了找到整合方法耗费大量时间及不必要的金钱。君子生非异也，善假于物也。学，然后知不足。加油吧！少年！用好单细胞测序及分析这个技术，为人类健康研究做出自己的贡献！不足之处，恳请批评指正！

概述

本教程演示了如何使用 Seurat 4.0 分析单细胞空间转录组数据。虽然分析流程类似于单细胞RNA-seq分析^[1]的Seurat工作流程，但我们引入了更新的交互和可视化工具，特别强调了空间和分子信息的整合。本教程将涵盖以下任务：

标准化
降维和聚类
检测空间高变基因
交互式可视化
与单细胞RNA-seq数据整合
处理多个切片

对于第一个教程，我们分析使用Visium 技术^[2]从 10 x genomics生成的数据集。我们将在不远的将来扩展seurat，以处理其他数据类型，包括SLIDE-Seq^[3]、STARmap^[4]和MERFISH。^[5]

首先，我们加载Seurat和所需的其他包。

library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)

10x Visium

数据

这里，我们将使用已发布的使用Visium v1 chemistry产生的小鼠大脑切片的数据集。

可以在此处^[6]下载数据，然后使用Load10X_Spatial()将其加载到 Seurat 中。这一步读取spaceranger^[7]流程的输出，并返回一个包含 spot-level表达数据以及组织切片相关图像的 Seurat 对象。可以使用SeuratData 包^[8]轻松访问数据，如下所示。

InstallData("stxBrain")

brain <- LoadData("stxBrain", type = "anterior1")

数据预处理

通过基因表达数据执行的预处理步骤类似于传统的 scRNA-seq 。首先需要使数据标准化，以便较少数据之间测序深度的差异。分子计数/点的差异对于空间数据集来说可能很大，特别是如果整个组织的细胞密度存在差异。在这里看到巨大的异质性，这需要有效的标准化。

plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")
wrap_plots(plot1, plot2)

这些图表明，跨点分子计数的差异不仅在技术上，而且依赖于组织解剖学。例如，神经元（如皮质白质）缺失的组织区域，可反复显示较低的分子计数。因此，标准化方法（如函数LogNormalize()）可能会有问题，这些方法使每个数据点在标准化后具有相同的基础"大小"。

作为替代方案，我们建议使用 sctransform，该模型构建了基因表达的常规负二元模型，以便在保持生物差异的同时考虑技术。sctransform 使数据标准化，检测高变异基因，并将数据存储在SCTassay.

brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)

基因表达可视化

在 Seurat 中，我们具有探索空间数据固有的视觉性质并与之交互的功能。SpatialFeaturePlot()是Seurat 的FeaturePlot()扩展，并且可以在组织学上覆盖分子数据。例如，在小鼠大脑的这一数据集中，Hpca基因是一个强大的海马标记，Ttr是choroid plexus的标记。

SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))

Seurat 中的默认参数强调分子数据的可视化。但是，您也可以通过更改以下参数来调整点的大小（及其透明度），以改善组织学图像的可视化：

pt.size.factor-这将扩展点的大小。默认值为1.6
alpha-最低和最大的透明度。默认值为 c（1，1）。
尝试设置为 c（0.1， 1），以降低低表达点的透明度alpha

p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
p1 + p2

降维、聚类和可视化

然后，我们可以继续运行 RNA 表达数据上的降维和聚类，使用与用于 scRNA-seq 分析相同的工作流。

brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30)
brain <- FindClusters(brain, verbose = FALSE)
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)

然后，我们可以使用DimPlot()在 UMAP 空间中或使用SpatialDimPlot()可视化聚类的结果。

p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p1 + p2

由于有许多颜色，因此可视化哪个 voxel 属于哪个聚类可能具有挑战性。我们有几个策略来帮助解决此问题。设置label参数时，每个集群的中位数放置一个彩色框（参见上图）。

您也可以使用cells.highlight参数在SpatialDimPlot()。这对于区分单个群的局部空间非常有用，如下所示：

SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3, 
    5, 8)), facet.highlight = TRUE, ncol = 3)

交互式绘图

我们还建立了一些交互式绘图功能。SpatialDimPlot()和SpatialFeaturePlot()两者现在都有一个参数，当设置interactive为TRUE，将打开Rstudio窗格与shiny互动。下面的示例演示了一个交互式，您可以在点上悬停并查看细胞名称和当前marker。

SpatialDimPlot(brain, interactive = TRUE)

对于SpatialFeaturePlot()，设置交互式为TRUE，显示一个交互式窗格，您可以在其中调整点的透明度、点大小以及正在绘制的函数。在浏览数据后，选择已完成的按钮将返回最后一个活动图作为 ggplot 对象。

SpatialFeaturePlot(brain, features = "Ttr", interactive = TRUE)

LinkedDimPlot()函数将 UMAP 与组织图像连接，并允许进行交互式选择。例如，您可以选择 UMAP 绘图中的区域，并将突出显示图像表示中的相应点。

LinkedDimPlot(brain)

空间高变基因的识别

Seurat 提供两个工作流程来识别与组织内空间位置相关的分子特征。首先是根据组织内预先注释的解剖区域进行差异表达分析，这些区域可能由不受监督的聚类或先验的知识确定。在这种情况下，此策略将起作用，因为上述集群表现出明确的空间限制。

de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 6)
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)

另一种方法，是FindSpatiallyVariables()，是寻找空间模式中没有预先注释的基因。默认方法是method = 'markvariogram`。

brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000], 
    selection.method = "markvariogram")

现在，我们可视化前 6 个基因的表达。

top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))

解剖区域取子集

与single-cell对象一样，可以对象取子集。在这里，我们取出前额皮层。此过程还有助于将这些数据与下一节中的外皮层 scRNA-seq 数据集整合。首先，我们取群的子集，然后根据确切的位置进一步细分。subset后，我们可以在完整图像或裁剪图像上可视化皮质细胞。

cortex <- subset(brain, idents = c(1, 2, 3, 4, 6, 7))
# now remove additional cells, use SpatialDimPlots to visualize what to remove
# SpatialDimPlot(cortex,cells.highlight = WhichCells(cortex, expression = image_imagerow > 400 |
# image_imagecol < 150))
cortex <- subset(cortex, anterior1_imagerow > 400 | anterior1_imagecol < 150, invert = TRUE)
cortex <- subset(cortex, anterior1_imagerow > 275 & anterior1_imagecol > 370, invert = TRUE)
cortex <- subset(cortex, anterior1_imagerow > 250 & anterior1_imagecol > 440, invert = TRUE)

p1 <- SpatialDimPlot(cortex, crop = TRUE, label = TRUE)
p2 <- SpatialDimPlot(cortex, crop = FALSE, label = TRUE, pt.size.factor = 1, label.size = 3)
p1 + p2

与single-cell数据的整合

在~50um时， visium assay检测的点将包含多个细胞的表达特征。随着越来越多的可用 scRNA-seq 数据，用户可能有兴趣对每个空间 voxels "去卷积"来预测细胞类型的潜在组成。在准备此教程时，我们测试了各种去卷积和整合方法，使用了由 Allen 研究所生成的 14，000 个成年鼠皮质细胞的scRNA-seq 数据集^[9]。我们发现使用整合方法（而不是去卷积方法）更优越，可能是因为空间和单细胞数据集的噪声模型存在较大差异，整合方法针对这些差异进行了稳健的设计。因此，我们应用[ Seurat v3](https://www./cell/fulltext/S0092-8674(19 " Seurat v3")30559-8)中引入的基于"锚点"的sctransform 整合工作流，从而能够将细胞从参考集转到查询集进行注释。

首先加载数据（此处下载^[10]），预处理 scRNA-seq 参考集，然后执行标签转移。该过程为每个点输出每个scRNA-seq细胞的概率分类。我们将这些预测添加到新的Seurat 对象进行分析。

allen_reference <- readRDS("../data/allen_cortex.rds")

# note that setting ncells=3000 normalizes the full dataset but learns noise models on 3k cells
# this speeds up SCTransform dramatically with no loss in performance
library(dplyr)
allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>% 
    RunUMAP(dims = 1:30)

# After subsetting, we renormalize cortex
cortex <- SCTransform(cortex, assay = "Spatial", verbose = FALSE) %>% RunPCA(verbose = FALSE)
# the annotation is stored in the 'subclass' column of object metadata
DimPlot(allen_reference, group.by = "subclass", label = TRUE)

anchors <- FindTransferAnchors(reference = allen_reference, query = cortex, normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors, refdata = allen_reference$subclass, prediction.assay = TRUE, 
    weight.reduction = cortex[["pca"]], dims = 1:30)
cortex[["predictions"]] <- predictions.assay

现在，我们得到每个类每个点的预测分数。最感兴趣的是frontal cortex region的laminar excitatory neurons。在这里，我们可以区分这些神经元亚型的不同顺序层，例如：

DefaultAssay(cortex) <- "predictions"
SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)

基于这些预测分数，我们还可以预测其位置受空间限制的细胞类型。使用基于标记点过程的相同方法来定义空间变异基因，但使用细胞类型预测分数作为"标记"，而不是基因表达。

cortex <- FindSpatiallyVariableFeatures(cortex, assay = "predictions", selection.method = "markvariogram", 
    features = rownames(cortex), r.metric = 5, slot = "data")
top.clusters <- head(SpatiallyVariableFeatures(cortex), 4)
SpatialPlot(object = cortex, features = top.clusters, ncol = 2)

最后，我们表明，我们的整合程序能够恢复神经元和非神经元子集已知的空间定位模式，包括laminar excitatory, layer-1 astrocytes, and the cortical grey matter。

SpatialFeaturePlot(cortex, features = c("Astro", "L2/3 IT", "L4", "L5 PT", "L5 IT", "L6 CT", "L6 IT", 
    "L6b", "Oligo"), pt.size.factor = 1, ncol = 2, crop = FALSE, alpha = c(0.1, 1))

在seurat处理多个切片

鼠大脑的这个数据集包含另一个与大脑另一半相对应的切片。在这里，我们读入它，并执行相同的初始标准化。

brain2 <- LoadData("stxBrain", type = "posterior1")
brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)

为了处理同一 Seurat 对象中的多个切片，我们提供了merge。

brain.merge <- merge(brain, brain2)

这便能在下面的RNA表达数据上实现整合降维和聚类。

DefaultAssay(brain.merge) <- "SCT"
VariableFeatures(brain.merge) <- c(VariableFeatures(brain), VariableFeatures(brain2))
brain.merge <- RunPCA(brain.merge, verbose = FALSE)
brain.merge <- FindNeighbors(brain.merge, dims = 1:30)
brain.merge <- FindClusters(brain.merge, verbose = FALSE)
brain.merge <- RunUMAP(brain.merge, dims = 1:30)

最后，数据可以在单个 UMAP 绘图中共同可视化。默认情况下，将所有切片绘制为列，分组/基因作为行。

DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))

SpatialDimPlot(brain.merge)

SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))

Slide-seq

数据

这里，我们将分析使用小鼠海马的Slide-seq v2^[11]数据集。此教程将遵循与 10x Visium 数据的空间教程大致相同的结构。

您可以使用SeuratData 包轻松访问数据，如下所示。安装数据集后，您可以键入?ssHippo以查看用于创建 Seurat 对象的命令。

InstallData("ssHippo")

slide.seq <- LoadData("ssHippo")

数据预处理

通过基因表达数据对bead 的初始预处理步骤类似于其他空间seurat分析和传统的scRNA-seq实验。在这里，我们注意到许多beads 包含特别低的 UMI 计数，我们选择保留所有检测到的beads 进行下游分析。

plot1 <- VlnPlot(slide.seq, features = "nCount_Spatial", pt.size = 0, log = TRUE) + NoLegend()
slide.seq$log_nCount_Spatial <- log(slide.seq$nCount_Spatial)
plot2 <- SpatialFeaturePlot(slide.seq, features = "log_nCount_Spatial") + theme(legend.position = "right")
wrap_plots(plot1, plot2)

然后，我们使用sctransform^[12]进行数据标准化，并执行标准的 scRNA-seq 降维和聚类工作流。

slide.seq <- SCTransform(slide.seq, assay = "Spatial", ncells = 3000, verbose = FALSE)
slide.seq <- RunPCA(slide.seq)
slide.seq <- RunUMAP(slide.seq, dims = 1:30)
slide.seq <- FindNeighbors(slide.seq, dims = 1:30)
slide.seq <- FindClusters(slide.seq, resolution = 0.3, verbose = FALSE)

然后，我们可以在 UMAP space或bead coordinate space中可视化聚类的结果。

plot1 <- DimPlot(slide.seq, reduction = "umap", label = TRUE)
plot2 <- SpatialDimPlot(slide.seq, stroke = 0)
plot1 + plot2

SpatialDimPlot(slide.seq, cells.highlight = CellsByIdentities(object = slide.seq, idents = c(1, 
    6, 13)), facet.highlight = TRUE)

与scRNA-seq参考整合

为了促进Slide-seq数据集的细胞类型注释，我们利用了由Saunders *，Macosko *等人^[13]制作的现有小鼠单细胞RNA-seq海马数据集^[14]。数据可在此处^[15]作为已处理的Seurat对象下载，而原始计数矩阵可在DropViz网站上获得^[16]。

ref <- readRDS("../data/mouse_hippocampus_reference.rds")

本文的原始注释在Seurat对象的metadata 中提供。这些注释以几种“resolutions”提供，从大类（ref$class）到亚群ref$subcluster。出于本教程的目的，我们将对细胞类型注释（ref$celltype）进行修改，使之达到良好的平衡。

我们将从运行Seurat标签转移方法开始，以预测每个bead的主要细胞类型。

anchors <- FindTransferAnchors(reference = ref, query = slide.seq, normalization.method = "SCT", 
    npcs = 50)
predictions.assay <- TransferData(anchorset = anchors, refdata = ref$celltype, prediction.assay = TRUE, 
    weight.reduction = slide.seq[["pca"]], dims = 1:50)
slide.seq[["predictions"]] <- predictions.assay

然后，我们可以可视化一些主要类型的预测分数。

DefaultAssay(slide.seq) <- "predictions"
SpatialFeaturePlot(slide.seq, features = c("Dentate Principal cells", "CA3 Principal cells", "Entorhinal cortex", 
    "Endothelial tip", "Ependymal", "Oligodendrocyte"), alpha = c(0.1, 1))

slide.seq$predicted.id <- GetTransferPredictions(slide.seq)
Idents(slide.seq) <- "predicted.id"
SpatialDimPlot(slide.seq, cells.highlight = CellsByIdentities(object = slide.seq, idents = c("CA3 Principal cells", 
    "Dentate Principal cells", "Endothelial tip")), facet.highlight = TRUE)

识别空间变异基因

正如Visium教程中提到的那样，我们可以通过两种通用方法来识别空间变异基因：预先标注的解剖区域之间的差异表达检测或基因特征对空间位置的依赖性的统计信息。

在这里，我们演示后者，通过在FindSpatiallyVariableFeatures()中设置method = 'moransi'。

Moran's I计算总体空间自相关性，并给出一个统计量（类似于相关系数），该统计量用于衡量要素对空间位置的独立性。这使我们能够根据基因表达的变化对基因进行排名。为了便于快速估算此统计信息，我们实施了一种基本的分级策略，该策略将在Slide-seqpuck 上绘制一个矩形网格，并对每个分级中的基因和位置取平均值。x和y方向上的bins 分别由x.cuts和y.cuts参数控制。此外，虽然不是必需的，但安装可选Rfast2软件包（install.packages('Rfast2')），将通过更有效的方式来减少运行时间。

DefaultAssay(slide.seq) <- "SCT"
slide.seq <- FindSpatiallyVariableFeatures(slide.seq, assay = "SCT", slot = "scale.data", features = VariableFeatures(slide.seq)[1:1000], 
    selection.method = "moransi", x.cuts = 100, y.cuts = 100)

现在，我们可视化由Moran's I识别的前6个基因的表达。

SpatialFeaturePlot(slide.seq, features = head(SpatiallyVariableFeatures(slide.seq, selection.method = "moransi"), 
    6), ncol = 3, alpha = c(0.1, 1), max.cutoff = "q95")