Seurat处理Visium数据

0. 原理简介

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

1. 载入数据

这里使用的是演示数据集，如果是自己的数据，可以使用 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对象

View(brain)

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

2. 数据预处理

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)

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一模一样

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

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. 细分解剖区域

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)

对空间组提出的亚群数据重新进行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)

整合单细胞和空间组学数据

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

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