欢迎关注同名公主号：BBio

流程包括单切片的基本分析、scRNA整合、多切片整合。

参考：https://satijalab.org/seurat/articles/spatial_vignette.html

//DataSet: 10x Visium

https://www.10xgenomics.com/resources/datasets/mouse-brain-serial-section-1-sagittal-anterior-1-standard-1-0-0

#image
wget https://cf.10xgenomics.com/samples/spatial-exp/1.0.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz
#filtered_feature_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/spatial-exp/1.0.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5

//Data preprocessing

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

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

pic <- function(pic, png, width=480, height=480){
    png(png, width=width, height=height)
    print(pic)
    dev.off()
}

brain <- Load10X_Spatial("data", 
    filename = "filtered_feature_bc_matrix.h5", 
    assay = "Spatial",
    slice = "slice1",
    filter.matrix = TRUE)

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

image-20220515144640284.png

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

作为一种替代方法，我们建议使用sctransform (Hafemeister和Satija, Genome Biology 2019)，该方法构建正则化负二项基因表达模型，以便在保留生物方差的同时解释技术伪像。有关sctransform的更多细节，请参阅这里的论文和Seurat插图。sctransform将数据归一化，检测高方差特征，并将数据存储在SCT分析中。

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

//How do results compare to log-normalization?

为了探讨标准化方法的差异，我们研究sctransform和log-normalization的结果如何与UMIs的数量相关。为了进行比较，我们首先重新运行sctransform来存储所有基因的值，并通过NormalizeData()运行标准化。

# 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))
pic(p1+p2, "CorrelationPlot.png", width=960)

image-20220515144725646.png

对于上面的箱线图，我们计算每个特征(基因)与UMIs数量(这里的nCount_Spatial变量)的相关性。然后，我们根据基因的平均表达将它们分组，并生成这些相关性的箱线图。可以看到， log-normalization未能充分地将前三组中的基因归一化，这表明技术因素继续影响高表达基因的归一化表达估计。相反，sctransform的规范化大大减轻了这种影响。

//Gene expression visualization

例如，在这组小鼠大脑数据中，Hpca基因是一个强大的海马体标记，而Ttr是脉络膜丛的标记。

p <- SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))
pic(p, "FeaturePlot_Hpca_Ttr.png", width=960)
p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
pic(p1+p2, "FeaturePlot_Ttr.png", width=960)

image-20220515144757399.png

image-20220515144823034.png

//Dimensionality reduction, clustering, and visualization

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)
pic(p1+p2, "cluster.png", width=960)

p <- SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3,5, 8)), facet.highlight = TRUE, ncol = 3)
pic(p, "cluster_highlight.png", width=960)

image-20220515144901723.png

image-20220515144919739.png

//Identification of Spatially Variable Features

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

de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 6)
p <- SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)
pic(p, "FindMarkers.png", width=480*3)

image-20220515144954018.png

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

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

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

top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 6)
p <- SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))
pic(p, "FindSpatiallyVariableFeatures.png", width=480*1.5)

image-20220515145027974.png

//Subset out anatomical regions

与单格对象一样，您可以对对象进行子集化，以关注数据的子集。这里，我们大致划分了额叶皮层。这一过程也有助于下一节中这些数据与皮质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)
pic(p1+p2, "cortex.png", width=480*2)

image-20220515145051589.png

//Integration with single-cell data

在~50um时，来自visium实验的斑点将包含多个细胞的表达谱。对于可获得scRNA-seq数据的越来越多的系统，用户可能有兴趣对每个空间体素进行“解卷积”，以预测细胞类型的底层组成。在准备这篇插图时，我们使用艾伦研究所使用SMART-Seq2协议生成的14000个成年小鼠皮质细胞分类的参考scRNA-seq数据集，测试了多种解卷积和整合方法。我们始终发现，使用集成方法(与反卷积方法相反)可以获得更好的性能，这可能是因为空间和单细胞数据集的噪声模型存在本质上的差异，而集成方法是专门设计来应对这些差异的。因此，我们应用了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
p <- DimPlot(allen_reference, group.by = "subclass", label = TRUE)
pic(p, "ref.png")

image-20220515145126636.png

现在我们得到每个班级每个位置的预测分数。在额叶皮层区域特别有趣的是层流兴奋神经元。在这里，我们可以区分这些神经元亚型的不同顺序层，例如:

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"
p <- SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
pic(p, "order.png")

image-20220515145220362.png

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

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

image-20220515145252635.png

最后，我们展示了我们的整合过程能够恢复已知的神经元和非神经元亚群的空间定位模式，包括层兴奋性、第1层星形胶质细胞和皮质灰质。

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

image-20220515145500145.png

//Working with multiple slices in Seurat

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

#https://www.10xgenomics.com/resources/datasets/mouse-brain-serial-section-1-sagittal-posterior-1-standard-1-1-0
brain2 <- Load10X_Spatial("posterior1", 
    filename = "filtered_feature_bc_matrix.h5", 
    assay = "Spatial",
    slice = "slice2",
    filter.matrix = TRUE)

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

brain.merge <- merge(brain, brain2)

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)

p <- DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))
pic(p, "merge.png", width=480*2)

image-20220515145618942.png

p <- SpatialDimPlot(brain.merge)
pic(p, "merge1.png", width=480*2)
p <- SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))
pic(p, "merge2.png", width=480*2)

image-20220515145706020.png

image-20220515145725774.png