跟着Seurat官网学空间转录组（10X数据集为例，包含细胞映射，无监督学习）

这是一个使用Seurat进行空间转录组分析的教程的摘要，包括每个部分的代码：

本教程演示如何使用 Seurat (>=3.2) 分析空间分辨的 RNA-seq 数据。虽然分析流程类似于用于单细胞 RNA-seq 分析的 Seurat 工作流程，但我们引入了更新的交互和可视化工具，特别强调空间和分子信息的整合。本教程将涵盖以下任务，我们相信这些任务对于许多空间分析来说都是常见的：

正常化
降维和聚类
检测空间变化的特征
交互式可视化
与单细胞 RNA-seq 数据整合
使用多个切片
对于我们的第一个小插图，我们分析了使用 10x Genomics 的Visium 技术生成的数据集。我们将在不久的将来扩展 Seurat 以使用其他数据类型，包括SLIDE-Seq、STARmap和MERFISH。

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

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

数据集

在这里，我们将使用最近发布的使用 Visium v​​1 化学生成的小鼠矢状脑切片数据集。有两个连续的前部分和两个（匹配的）连续的后部分。

您可以在此处下载数据，并使用Load10X_Spatial()函数将其加载到 Seurat 中。这将读取spaceranger管道的输出，并返回一个 Seurat 对象，该对象包含点级表达数据以及组织切片的关联图像。您还可以使用我们的SeuratData 包来轻松访问数据，如下所示。安装数据集后，您可以键入?stxBrain以了解更多信息。

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

数据预处理

我们通过基因表达数据在现场执行的初始预处理步骤类似于典型的 scRNA-seq 实验。我们首先需要对数据进行归一化，以解释数据点之间测序深度的差异。我们注意到，对于空间数据集，分子计数/斑点的差异可能很大，尤其是在整个组织的细胞密度存在差异的情况下。我们在这里看到了大量的异质性，这需要有效的标准化。

这些图表明，跨点分子计数的差异不仅在本质上是技术性的，而且还取决于组织解剖学。例如，神经元耗尽的组织区域（如皮质白质）可重复地表现出较低的分子计数。因此，LogNormalize()强制每个数据点在规范化后具有相同的基础“大小”的标准方法（例如函数）可能会出现问题。

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)

作为替代方案，我们建议使用 sctransform（Hafemeister 和 Satija，Genome Biology 2019），它构建正则化的基因表达负二项式模型，以便在保留生物学差异的同时考虑技术人工制品。有关 sctransform 的更多详细信息，请参阅此处的论文和此处的Seurat vignette 。sctransform 规范化数据，检测高方差特征，并将数据存储在分析中SCT。

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

基因表达可视化

在 Seurat 中，我们具有探索空间数据固有视觉特性并与之交互的功能。SpatialFeaturePlot()Seurat 中的功能扩展了，FeaturePlot()可以将分子数据叠加在组织组织学之上。例如，在小鼠大脑的这个数据集中，基因 Hpca 是海马体的强标记，Ttr 是脉络丛的标记。

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

library(ggplot2)
plot <- SpatialFeaturePlot(brain, features = c("Ttr")) + theme(legend.text = element_text(size = 0),
    legend.title = element_text(size = 20), legend.key.size = unit(1, "cm"))
jpeg(filename = "../output/images/spatial_vignette_ttr.jpg", height = 700, width = 1200, quality = 50)
print(plot)
dev.off()

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

pt.size.factor- 这将缩放斑点的大小。默认为 1.6
alpha- 最小和最大透明度。默认为 c(1, 1)。
尝试设置为alphac(0.1, 1)，以降低具有较低表达的点的透明度

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)
p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p1 + p2

DimPlot()然后，我们可以在 UMAP 空间（使用）或使用 覆盖在图像上可视化聚类的结果SpatialDimPlot()。

由于有多种颜色，因此很难可视化哪个体素属于哪个簇。我们有一些策略可以帮助解决这个问题。设置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 查看器窗格。下面的示例演示了一个交互SpatialDimPlot()，您可以在其中将鼠标悬停在点上并查看单元格名称和当前标识类（类似于之前的do.hover行为）。

SpatialDimPlot(brain, interactive = TRUE)

对于SpatialFeaturePlot()，将交互式设置为TRUE调出一个交互式窗格，您可以在其中调整点的透明度、点大小以及Assay绘制的特征。探索数据后，选择完成按钮将返回最后一个活动图作为 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) 受到Trendsceek的启发，它将空间转录组学数据建模为标记点过程并计算“变差函数”，它识别表达水平取决于其空间位置的基因。更具体地说，此过程计算 gamma® 值，测量相隔一定“r”距离的两个点之间的相关性。默认情况下，我们在这些分析中使用“5”的 r 值，并且只计算可变基因的这些值（其中变异是独立于空间位置计算的）以节省时间。

我们注意到文献中有多种方法可以完成此任务，包括SpatialDE和Splotch。我们鼓励感兴趣的用户探索这些方法，并希望在不久的将来增加对它们的支持。

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

现在我们可视化该度量确定的前 6 个特征的表达。

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

划分出解剖区域

与单细胞对象一样，您可以对对象进行子集化以专注于数据子集。在这里，我们大致划分了额叶皮层。此过程还有助于将这些数据与下一节中的皮质 scRNA-seq 数据集集成。首先，我们取一个集群子集，然后根据确切位置进一步分割。子集化后，我们可以在完整图像或裁剪图像上可视化皮质细胞。

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

与单细胞数据整合

在 ~50um，来自 visium 测定的斑点将包含多个细胞的表达谱。对于越来越多的可用 scRNA-seq 数据的系统列表，用户可能有兴趣“解卷积”每个空间体素以预测细胞类型的潜在组成。在准备这个小插图时，我们使用参考 scRNA-seq 数据集测试了多种去卷积和整合方法来自艾伦研究所的约 14,000 个成年小鼠皮层细胞分类学，使用 SMART-Seq2 协议生成。我们始终发现使用积分方法（与反卷积方法相反）具有优越的性能，这可能是因为表征空间和单细胞数据集的噪声模型有很大不同，并且积分方法专门设计为对这些差异具有鲁棒性。因此，我们应用了Seurat v3中引入的基于“锚点”的集成工作流，该工作流可以实现注释从引用到查询集的概率转移。因此，我们遵循此处介绍的标签传输工作流程，利用 sctransform 规范化，但预计将开发新方法来完成此任务。

我们首先加载数据（可在此处下载），预处理 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

现在我们得到每个班级每个地点的预测分数。对额叶皮层区域特别感兴趣的是层状兴奋性神经元。在这里，我们可以区分这些神经元亚型的不同顺序层，例如：

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 = "moransi",
    features = rownames(cortex), r.metric = 5, slot = "data")
top.clusters <- head(SpatiallyVariableFeatures(cortex, selection.method = "moransi"), 4)
SpatialPlot(object = cortex, features = top.clusters, ncol = 2)

最后，我们表明我们的综合程序能够恢复神经元和非神经元子集的已知空间定位模式，包括层状兴奋性、第 1 层星形胶质细胞和皮质灰质。

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))