思考题：

+ 如何将空间数据与表达数据关联在一起？
+ 有了空间转录组数据，如何与单细胞转录组数据联用？
+ 做了多层切片如何展示真实的三维空间的转录本信息？

随着转录组技术的发展，空间转录组已经正式走向商业化时代，作为单细胞数据分析的工具箱的Seurat与时俱进，也相应地开发了空间转录组分析的一套函数，让我们跟随卑微小王看看Seurat官网教程吧。

本教程演示如何使用Seurat v3.2分析空间解析的RNA-seq数据。虽然分析流程类似于Seurat的单细胞RNA-seq分析流程，但我们引入了交互可视化工具，特别强调了空间和分子信息的集成。本教程将介绍以下任务，我们相信这些任务在许多空间分析中都很常见:

• 归一化
• 降维与聚类
• 检测spatially-variable特性
• 交互式可视化
• 与单细胞RNA-seq数据集成
• 处理多个片（multiple slices）

我们分析了使用来自10x Genomics 的Visium技术（Visium technology）生成的数据集。我们将在不久的将来扩展Seurat以处理其他数据类型，包括SLIDE-Seq、STARmap和MERFISH。

安装

devtools::install_github("satijalab/seurat", ref = "spatial")

这种方法，只能只能好运。直接下载安装包，本地安装。在遇到问题的时候，假装这是你写的R包参考《R包开发》这本书，尝试本地重构：

devtools::load_all('H:\\singlecell\\Seurat\\seurat-master\\seurat')

但是，一般用不到这么复杂，下载安装包，本地安装基本就可以了：

[https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial](https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial)  # 下载地址
install.packages("H:/singlecell/Seurat/satijalab-seurat-v3.1.1-302-g1cb8a3d.tar.gz", repos = NULL, type = "source")

注意，别和之前的包安装冲突啦，这个3.2还处在开发中（2019-12-19）

library(Seurat)
library(SeuratData)
library(ggplot2)
library(cowplot)
library(dplyr)

数据集

和以往一样，seurat为了使其教程的可用，会提供测试的已发表的数据集。在这里，我们将使用最新发布的使用Visium v1化学生成的sagital小鼠大脑切片数据集。有两个连续的前段和两个(匹配的)连续的后段。

您可以在这里here下载数据，并使用Load10X_Spatial函数将其加载到Seurat。这将读取spaceranger管道的输出，并返回Seurat对象，该对象包含spot级别的表达数据以及相关的组织切片图像。您还可以使用我们的SeuratData包方便地访问数据，如下所示。安装数据集之后，可以键入?stxBrain以了解更多信息。

InstallData("stxBrain")

每次我在国内直接用这种方法下载数据集都没有成功，如上，下载安装包，本地安装：

 install.packages("H:/singlecell/Seurat/stxBrain.SeuratData_0.1.1.tar.gz", repos = NULL, type = "source")
brain <- LoadData("stxBrain", type = "anterior1")

当然，我们不推荐这种读取数据的方法，毕竟没有人把我们的数据也打包成一个R包的样子，我们拿到的是10 X Space Ranger的输出结果：

├── analysis
│?? ├── clustering
│?? ├── diffexp
│?? ├── pca
│?? ├── tsne
│?? └── umap
├── cloupe.cloupe
├── filtered_feature_bc_matrix
│?? ├── barcodes.tsv.gz
│?? ├── features.tsv.gz
│?? └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5
├── metrics_summary.csv
├── molecule_info.h5
├── possorted_genome_bam.bam
├── possorted_genome_bam.bam.bai
├── raw_feature_bc_matrix
│?? ├── barcodes.tsv.gz
│?? ├── features.tsv.gz
│?? └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5
├── spatial  # 空间信息全在这 :这些文件是用户提供的原始全分辨率brightfield图像的下采样版本。
下采样是通过box滤波实现的，它对全分辨率图像中像素块的RGB值进行平均，得到下采样图像中一个像素点的RGB值。
│?? ├── aligned_fiducials.jpg   这个图像的尺寸是tissue_hires_image.png。由基准对齐算法发现的基准点用红色高亮显示。此文件对于验证基准对齐是否成功非常有用。
│?? ├── detected_tissue_image.jpg
│?? ├── scalefactors_json.json
│?? ├── tissue_hires_image.png 图像的最大尺寸为2,000像素 
│?? ├── tissue_lowres_image.png 图像的最大尺寸为600像素。
│?? └── tissue_positions_list.csv
└── web_summary.html

其实只要 Space Ranger的输出结果 filtered_feature_bc_matrix.h5文件和一个spatial文件夹就可以了，一如？stxBrain中给出的示例：

setwd("H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/")
  # Load the expression data
  expr.url <- 'H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5'
  expr.data <- Seurat::Read10X_h5(filename =  expr.url )
  anterior1 <- Seurat::CreateSeuratObject(counts = expr.data, project = 'anterior1', assay = 'Spatial')
  anterior1$slice <- 1
  anterior1$region <- 'anterior'
  # Load the image data
  img.url <- 'H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz'
  untar(tarfile =  img.url)
  img <- Seurat::Read10X_Image(image.dir = 'H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/spatial')
  Seurat::DefaultAssay(object = img) <- 'Spatial'
  img <- img[colnames(x = anterior1)]
  anterior1[['image']] <- img

  anterior1

An object of class Seurat 
31053 features across 2696 samples within 1 assay 
Active assay: Spatial (31053 features)

}

如果这两个文件在同一个文件下的话，也可以这样读入：

list.files("H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio") # 注意命名要正确

# [1] "filtered_feature_bc_matrix.h5"                                  "spatial"                                                       
# [3] "V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5" "V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz"  

brain<-Seurat::Load10X_Spatial("H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio") 
？Load10X_Spatial
brain


An object of class Seurat 
31053 features across 2696 samples within 1 assay 
Active assay: Spatial (31053 features)

空间数据如何存储在Seurat中?

来自10x的visium数据包括以下数据类型:

• 通过基因表达矩阵得到一个点(spot )
• 组织切片图像(采集数据时H&E染色)
• 用于显示的原始高分辨率图像与低分辨率图像之间的比例因子。

在Seurat对象中，spot by基因表达矩阵与典型的“RNA”分析类似，但包含spot水平，而不是单细胞水平的数据。图像本身存储在Seurat对象中的一个images 槽（slot ）中。图像槽还存储必要的信息，以将斑点与其在组织图像上的物理位置相关联。

数据预处理

在spot中基因表达数据进行初始的预处理步骤与典型的scRNA-seq相似。我们首先需要对数据进行规范化，以考虑数据点之间测序深度的差异。我们注意到，对于空间数据集来说，分子数/点的差异可能是巨大的，特别是在整个组织的细胞密度存在差异的情况下。我们在这里看到大量的异质性，这需要有效的标准化。

熟悉的函数名~

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

这些图表明，分子计数（molecular counts）在点间的差异不仅是技术性的，而且还依赖于组织解剖。例如，组织中神经元消耗的区域(如皮层白质)，可再生地显示出较低的分子计数。因此，标准方法(如LogNormalize函数)可能会有问题，因为它会强制每个数据点在标准化之后具有相同的底层“大小”。

作为一种替代方法，我们推荐使用sctransform (Hafemeister和Satija，已出版)，它构建了基因表达的正则化负二项模型，以便在保留生物差异的同时考虑技术因素。有关sctransform的更多信息，请参见 here的预印和here的Seurat教程。sctransform将数据归一化，检测高方差特征，并将数据存储在SCT分析中。

brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)

sct 与log-归一化相比，结果如何?

为了探究规范化方法的不同，我们研究了sctransform和log规范化结果如何与UMIs的数量相关。为了进行比较，我们首先重新运行sctransform来存储所有基因的值(这将会比较慢)，并通过NormalizeData运行一个log规范化过程。

# rerun normalization to store sctransform residuals for all genes
brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
# also run standard log normalization for comparison
brain <- NormalizeData(brain, verbose = FALSE, assay = "Spatial")

# Computes the correlation of the log normalized data and sctransform residuals with the number
# of UMIs
brain <- GroupCorrelation(brain, group.assay = "Spatial", assay = "Spatial", slot = "data", do.plot = FALSE)
brain <- GroupCorrelation(brain, group.assay = "Spatial", assay = "SCT", slot = "scale.data", do.plot = FALSE)
p1 <- GroupCorrelationPlot(brain, assay = "Spatial", cor = "nCount_Spatial_cor") + ggtitle("Log Normalization") + 
    theme(plot.title = element_text(hjust = 0.5))
p2 <- GroupCorrelationPlot(brain, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransform Normalization") + 
    theme(plot.title = element_text(hjust = 0.5))
plot_grid(p1, p2)

对于上面的箱形图，我们计算每个特征(基因)与UMIs数量(这里的nCount_Spatial变量)的相关性。然后，我们根据基因的平均表达将它们分组，并生成这些相关性的箱形图。您可以看到，log-normalization未能充分规范化前三组的基因，这表明技术因素影响高表达基因的规范化表达估计值。相反，sctransform规范化降低了这种效果。差别真的很大么？读者诸君自行判断。

基因表达可视化

在Seurat v3.2中，我们加入了新的功能来探索和与空间数据固有的可视化特性。Seurat的SpatialFeaturePlot功能扩展了FeaturePlot，可以将表达数据覆盖在组织组织上。例如，在这组小鼠大脑数据中，Hpca基因是一个强的海马marker ，Ttr是一个脉络丛marker 。

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

Seurat的默认参数强调分子数据的可视化。然而，你也可以调整斑点的大小(和它们的透明度)来改善组织学图像的可视化，通过改变以下参数:

• pt.size。因素-这将比例大小的斑点。默认为1.6
• alpha -最小和最大透明度。默认是c(1,1)
• 尝试设置为alpha c(0.1, 1)，以降低表达式较低的点的透明度

降维、聚类和可视化

然后，我们可以使用与scRNA-seq分析相同的工作流，对RNA表达数据进行降维和聚类。我们可以在UMAP空间(使用DimPlot)或使用SpatialDimPlot将分群的结果显示在图像上。

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)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
22:08:28 UMAP embedding parameters a = 0.9922 b = 1.112
22:08:28 Read 2696 rows and found 30 numeric columns
22:08:28 Using Annoy for neighbor search, n_neighbors = 30
22:08:28 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:08:30 Writing NN index file to temp file C:\Users\ADMINI~1\AppData\Local\Temp\Rtmp08OC2k\file3ddc1e802bb6
22:08:30 Searching Annoy index using 1 thread, search_k = 3000
22:08:31 Annoy recall = 100%
22:08:32 Commencing smooth kNN distance calibration using 1 thread
22:08:33 Initializing from normalized Laplacian + noise
22:08:33 Commencing optimization for 500 epochs, with 106630 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
22:08:50 Optimization finished

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

因为有许多颜色，所以可视化哪个体素属于哪个簇是很有挑战性的。我们有一些策略来帮助解决这个问题。通过设置label参数，可以在每个集群的中间位置放置一个彩色框(参见上面的图)以及do.hover。SpatialDimPlot的悬停参数允许交互式查看当前的spot标识。

# move your mouse
SpatialDimPlot(brain, do.hover = TRUE)

Warning messages:
1: In if (is.na(col)) { :
  the condition has length > 1 and only the first element will be used
2: In if (is.na(col)) { :
  the condition has length > 1 and only the first element will be used
3: `error_y.color` does not currently support multiple values. 
4: `error_x.color` does not currently support multiple values. 
5: `line.color` does not currently support multiple values. 
6: The titlefont attribute is deprecated. Use title = list(font = ...) instead. 

你也可以使用cells.highlight，用于在空间坐标图上划分感兴趣的特定单元格。这对于区分单个集群的空间定位非常有用，如下所示:

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

LinkedDimPlot和LinkedFeaturePlot函数支持交互式可视化。这些图将UMAP表示与组织图像表示联系起来，并允许交互选择。例如，您可以在UMAP图中选择一个区域，图像表示中相应的点将突出显示。

LinkedDimPlot(brain)

空间变量特征的识别

Seurat提供了两个工作流程来识别与组织空间位置相关的分子特征。第一种是根据组织内预先标注的解剖区域进行差异表达，这种差异表达可以通过非监督聚类或先验知识来确定。这种策略在这种情况下有效，因为上面的集群显示出明显的空间差异。

de_markers <- FindMarkers(brain, ident.1 = 4, ident.2 = 6)

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=37s  

SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)

在FindSpatiallyVariables中实现的另一种方法是在没有预注释的情况下搜索显示空间模式的特性。默认的方法(method = 'markvariogram ')受到 Trendsceek,的启发，后者将空间转录组学数据建模为标记点过程，并计算一个' variogram '，它识别其表达水平取决于其空间位置的基因。更具体地说，这个过程计算伽玛(r)值，测量两个点之间一定的“r”距离的相关性。默认情况下，我们在这些分析中使用的r值为‘5’，并且只计算可变基因的这些值(其中的变异是独立于空间位置计算的)，以节省时间。

现在，我们可视化的表达前6个特征确定了这一措施。

brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000], 
    selection.method = "markvariogram")
top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))

可视化解剖区域的子集

与单细胞对象一样，您可以对该对象进行子集设置，以将重点放在数据的子集上。在这里，我们大致划分了额叶皮质。这个过程也促进了这些数据与下一节的皮层scRNA-seq数据集的整合。首先，我们取集群的一个子集，然后根据精确的位置进一步细分。设置好亚组后，我们可以在完整图像或裁剪图像上看到皮质细胞。

cortex <- subset(brain, idents = c(1, 2, 3, 5, 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, image_imagerow > 400 | image_imagecol < 150, invert = TRUE)
cortex <- subset(cortex, image_imagerow > 275 & image_imagecol > 370, invert = TRUE)
cortex <- subset(cortex, image_imagerow > 250 & image_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)
plot_grid(p1, p2)

与单细胞数据关联分析（空间细胞类型定义）

在~50um时，visium检测的斑点将包含多个细胞的表达谱。对于越来越多可用scRNA-seq数据的系统，用户可能有兴趣“反卷积”每个空间体素，以预测单元类型的底层组成。在准备这篇文章的过程中，我们测试了各种各样的脱卵方法和整合方法（decovonlution and integration methods），使用的是来自Allen研究所的参考scRNA-seq数据集（reference scRNA-seq dataset），其中包含了大约14000只成年小鼠的皮层细胞分类，并使用SMART-Seq2 protocol 生成。

我们一致认为，使用集成方法(与反褶积方法相反)可以获得更好的性能，这可能是因为空间和单细胞数据集的噪声模型本质上是不同的，而集成方法的特殊设计是为了对这些差异具有鲁棒性。因此我们应用“锚”的集成工作流一如最近在 Seurat v3介绍的，使注释的概率从一个reference 数据query 数据集转移。因此,我们遵循这里介绍的标签转移流程,利用sctransform正常化,但预测新方法被开发来完成这项任务。

我们首先加载数据( here提供下载，我只能假装我下载成功了^_)，预处理scRNA-seq 参考数据集，然后执行标签传输。该过程为每个spot输出每个scRNA-seq派生类的概率分类。我们将这些预测添加到Seurat对象中作为新的试验。

allen_reference <- readRDS("~/Downloads/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"]])
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", features = rownames(cortex), 
    r.metric = 5, slot = "data")
top.clusters <- head(SpatiallyVariableFeatures(cortex), 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))

处理Seurat中的多个slices

这个老鼠大脑的数据集包含了另一个与大脑另一半相对应的切片。这里我们读取它并执行相同的初始规范化。

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图中共同可视化数据。SpatialDimPlot和SpatialFeaturePlot将默认将所有片绘制为列，将groupings/features 绘制为行。

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

SpatialDimPlot(brain.merge)   

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

我们要感谢Nigel Delaney和Stephen Williams对Seurat 分析空间数据代码的有益反馈和贡献。

我是教程官网，点我直达