0. 原理简介

空间转录组的概念：将基因表达与组织切片的免疫组化图像相结合，从而将组织内不同细胞的基因表达信息定位到组织的原始空间上，达到直观检测组织中不同部位基因表达差异的目的。可以避免单细胞转录组样本消化对细胞的影响，结合组织定位也可以更好的判断细胞类型。
核心：单细胞转录组测序+组织定位

空间转录组实验流程

空转芯片：
• 一个载玻片可容纳4个切片(捕获区域)，制备4个文库；
• 每个6.5mm x 6.5mm捕获区域, 有5000个barcode标记的ST位点；
• 每个ST位点有百万个UMI探针；
• 每个位点直径55μm，点与点之间相距100μm;
• 每个位点覆盖1-10个细胞，检测数千个基因;

10x公司的Visium空间转录组方法与10x的单细胞测序在原理上很相似。单细胞测序是用胶珠和油包水方法，把细胞分隔开，同时又用DNA条形码保留单细胞信息。Visium空间转录组则是把切片在芯片上展开，在空间上用条形码来保留切片上每个小点的空间位置信息。

1. 载入数据

这里使用的是Visium v1得到的小鼠脑组织切片演示数据集。包含两个连续的anterior切片和2个连续的posterior切片。如果是自己的数据，可以使用 Load10X_Spatial() 函数读入（filtered_feature_bc_matrix.h5）。

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

InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
#只读取了anterior1这个sample的数据

我们先来看一下空间转录组数据的Seurat对象
有2696个spots（一共是5000个）

View(brain)

和单细胞转录组数据是差不多的

2. 数据预处理

没有做质控哦（过滤方法跟单细胞转录组一样的，过滤掉的细胞在切片上会看到spots变灰了）

brain[['percent.mt']] <- PercentageFeatureSet(brain,pattern = '^mt-')
brain[['percent.rb']] <- PercentageFeatureSet(brain,pattern = '^Rp[ls]')
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)

小提琴图展示了每个spot测到的reads数，一个点是一个spot。右图是看每个spot测到的reads在空间上的表达，和FeaturePlot实际上是一样的。可以看到不同解剖部位的细胞的基因表达量还是存在明显差异的

上图表明，spot测得reads数的变化不仅是技术上的问题，而且还取决于组织的解剖结构。As a result, standard approaches (such as the LogNormalize function), which force each data point to have the same underlying ‘size’ after normalization, can be problematic.

作为替代方案，我们建议使用sctransform，它构建了基因表达的正则化负二项式模型，以便在保留生物学差异的同时考虑技术伪像。sctransform可以对数据进行归一化，检测高变异特征并将数据存储在SCTassay中。

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

拓展：
How do results compare to log-normalization
为了探究标准化方式的差异，作者检测了sctransform 和log normalization和UMIs的相关性

# 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))
p1 + p2

可以看到SCTransform表现的更好

3. 基因表达可视化

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

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

4. 降维，聚类和可视化

和scRNA-seq一模一样（Louvain聚类）

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)

这种聚类方法最容易实现，但是不够精确。（BayesSpace聚类等效果会更好）

p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p3 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)&ggsci::scale_fill_igv()
p1|p2|p3 #对接ggsci自定义颜色

SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3,
    5, 8)), facet.highlight = TRUE, ncol = 3)
# 相当于单细胞的split.by

5. 鉴定空间高变基因

FindMarkers找的是上面不同分群间的差异基因

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

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

FindSpatiallyVariableFeatures找的是空间高变基因

空间高变基因指的是表达与空间位置相关的基因，也可以说是具有空间偏好性的基因。空间转录组学允许研究人员调查基因表达趋势如何在空间上变化，从而确定基因表达的空间模式。寻找空间高变基因的算法有很多，比如SPARK，SpatialDE等。

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

6. 细分解剖区域

这个其实用Loupe Browser手动圈会更方便一点

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

7. 和单细胞数据整合

导入单细胞数据

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)

对空间数据里面选出来的亚群（brain里面挑出来的cortex）数据重新进行SCT标准化

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

整合单细胞和空间组学数据（基于锚点转换的去卷积）

注意这里是按组织结构去做注释的（从空转数据中提出了cortex）。一般不建议使用大的单细胞数据集去做参考。原因一是速度会很慢，二是容易引入批次效应让结果不准确。所以尽量还是找一对一匹配的数据集做锚点转换的去卷积，也就是测了空转的数据，最好有相应的同一个样本的单细胞数据做映射。整合了多个样品的单细胞数据集虽然可能细胞类型更全，但会引入批次效应。

##寻找锚点
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

得到的prediction矩阵是行为输入的单细胞的细胞类型(24种)，列是每一个spot，矩阵的内容是每个spot的该细胞类型评分

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)

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

可以看到不同的单细胞细胞群在空间数据上的映射都存在一个评分。这个评分其实不是可信度，可以将其理解为细胞占比。我们知道每个spot其实是多个细胞的混合，一个spot中有些细胞占比多，有些细胞占比少。比如说上图第一行第三张图的L4范围只有0.1-0.3，就说明最红的这些spot中，每个spot中可能最多也只有30%的细胞是L4，其他70%可能是其他细胞，需要结合其他细胞群的图来一起看。

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