跟着Seurat官网学Visium HD空转分析（一）标准分析流程

生信菜鸟团

发布于 2024-07-31 19:23:53

6410

发布于 2024-07-31 19:23:53

文章被收录于专栏：生信菜鸟团

各位小伙伴，好久不见！最近我忙完了博士毕业相关的事，离开武汉刚回到杭州，六年的硕博生活终于结束啦！还在坑里的朋友们要坚持住！

Visium HD在国内已发售近半年了，全球首篇HD文章也已预印，题为《Characterization of immune cell populations in the tumor microenvironment of colorectal cancer using high definition spatial profiling》，小伙伴们开始上手HD实战了吗？三月底的时候我测试了10X官网发布的示例数据【全网首发 | Visium HD空转数据开箱测试】，那会Seurat并不兼容HD分析。而在五月份的时候，Seurat针对HD进行了更新，并在官网发布了一个分析教程，详见https://satijalab.org/seurat/articles/visiumhd_analysis_vignette：

这里Seurat官网再次强调8 um分辨率的HD数据是比较适合进行下游分析的，但Seurat支持同时加载不同bin的数据，并将它们存储在一个Seurat对象中作为多个assays。

此外，Seurat提供了一些HD的空间分析流程，特别是:

Unsupervised clustering
Identification of spatial tissue domains (BANKSY)
Subsetting spatial regions
Integration with scRNA-seq data (RCTD)
Comparing the spatial localization of different cell types

熟悉我推文的小伙伴会发现，这些内容其实我基本都介绍过：

无监督聚类，即标准流程，包括降维聚类分析、差异表达分析以及亚群注释等等；
空间组织域的识别，这里使用的是BANKSY分析流程，可参考【空转数据分析之生态位聚类算法Banksy】；
整合scRNA-seq数据，即反卷积算法，这个我写了一系列的教程，例如：这里Seurat官方使用了RCTD来整合HD数据和单细胞数据。

上述分析算法其实并不局限于HD这个技术。可以说单细胞数据分析、传统的Visium空转分析、HD空转分析甚至其他平台的空转，都有非常多交叉的内容，我们在学习的过程中需要融会贯通，举一反三。

一. 安装/更新Seurat v5

这里使用的是Seurat v5版本：

#Seurat v5
install.packages('Seurat')
library(Seurat)

#We recommend users install these along with Seurat
setRepositories(ind = 1:3, addURLs = c('https://satijalab.r-universe.dev', 'https://bnprks.r-universe.dev/'))
install.packages(c("BPCells", "presto", "glmGamPoi"))

#We also recommend installing these additional packages, which are used in our vignettes, and enhance the functionality of Seurat:
# Install the remotes package
if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}
install.packages('Signac')
remotes::install_github("satijalab/seurat-data", quiet = TRUE)
remotes::install_github("satijalab/azimuth", quiet = TRUE)
remotes::install_github("satijalab/seurat-wrappers", quiet = TRUE)

# packages required for Visium HD
install.packages("hdf5r")
install.packages("arrow")

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

二. 数据加载

示例数据是小鼠脑的HD空转数据，下载自10X官网：https://support.10xgenomics.com/spatial-gene-expression/datasets

这里下载了8 um和16 um的数据，可以通过bin.size进行指定读入同一个样本不同分辨率的数据：

localdir <- "/brahms/lis/visium_hd/mouse/new_mousebrain/"
object <- Load10X_Spatial(data.dir = localdir, bin.size = c(8, 16))

# Setting default assay changes between 8um and 16um binning
Assays(object)
DefaultAssay(object) <- "Spatial.008um"

这里和常规的单细胞数据分析一样，使用DefaultAssay(object) <- "Spatial.008um"指定默认使用的插槽表达矩阵。

检查一下加载好的数据：

vln.plot <- VlnPlot(object, features = "nCount_Spatial.008um", pt.size = 0) + theme(axis.text = element_text(size = 4)) + NoLegend()
count.plot <- SpatialFeaturePlot(object, features = "nCount_Spatial.008um") + theme(legend.position = "right")

# note that many spots have very few counts, in-part
# due to low cellular density in certain tissue regions
vln.plot | count.plot

三. 标准分析流程

质控？

这里Seurat官网没有给出一个质控的标准，但我在【全网首发 | Visium HD空转数据开箱测试】，提到：

“Visium HD空转数据的测序深度似乎比10X单细胞低一些。另外，在过滤低质量细胞过程中，我发现按照10X单细胞的过滤标准对于Visium HD空转数据而言不太合适。然后我查阅了华大的过滤标准，华大对percent.mt的设置为5，在这里我设置为15，其他指标和华大的过滤标准一致。”

在这里我贴一下我质控过滤的代码，仅供参考：

## 华大QC流程：
CRC.qc <- subset(CRC_data, subset = nFeature_Spatial >= 20 & 
                    nCount_Spatial>= 3 & percent.mt < 15)
CRC.qc #过滤后 18085 features across 497201 samples

sketch算法是Seurat v5的一个特色，特别适合用于大图谱数据的降维聚类分析，可参考：https://satijalab.org/seurat/articles/parsebio_sketch_integration。简要来说，sketch算法通过抽样取小子集的方式，可以对个小子集进行基础分析，以获得一个聚类和降维结果，然后通过ProjectData函数从抽样细胞中学习到的聚类标签和降维结果投影到整个数据集，这样可以大大节省内存和时间。

因此，这里可以重点学习一下Seurat v5是如何进行HD数据的标准分析流程（毕竟一个HD样本就拥有50W的单细胞）：

# 3.1 normalize 8um bins
DefaultAssay(object) <- "Spatial.008um"
object <- object %>% NormalizeData(verbose = F) %>%
  FindVariableFeatures(verbose = F) %>% 
  ScaleData(verbose = F)

# 3.2 Unsupervised clustering
# we select 50,0000 cells and create a new 'sketch' assay
object <- SketchData(
  object = object,
  ncells = 50000,
  method = "LeverageScore",
  sketched.assay = "sketch"
)

# switch analysis to sketched cells
DefaultAssay(object) <- "sketch"

# perform clustering workflow
object <- object %>% FindVariableFeatures() %>% 
  ScaleData() %>% 
  RunPCA(assay = "sketch", reduction.name = "pca.sketch") %>% 
  FindNeighbors(assay = "sketch", reduction = "pca.sketch", dims = 1:50) %>% 
  FindClusters(cluster.name = "seurat_cluster.sketched", resolution = 3) %>% 
  RunUMAP(reduction = "pca.sketch", reduction.name = "umap.sketch", return.model = T, dims = 1:50)

可以看到，在进行sketch进行抽样运算之后得到一个新的插槽sketch，然后对这个新的插槽走标准分析流程即可。

然后，我们可以使用ProjectData函数从抽样细胞（示例代码是50000个）中学习到的聚类标签和降维(PCA和UMAP)投影到整个数据集：

object <- ProjectData(
  object = object,
  assay = "Spatial.008um",
  full.reduction = "full.pca.sketch",
  sketched.assay = "sketch",
  sketched.reduction = "pca.sketch",
  umap.model = "umap.sketch",
  dims = 1:50,
  refdata = list(seurat_cluster.projected = "seurat_cluster.sketched")
)

然后对子集结果和完整数据的结果进行可视化比较，这里要充分理解DefaultAssay(object)函数的目的和意义：

DefaultAssay(object) <- "sketch"
Idents(object) <- "seurat_cluster.sketched"
p1 <- DimPlot(object, reduction = "umap.sketch", label = F) + ggtitle("Sketched clustering (50,000 cells)") + theme(legend.position = "bottom")

# switch to full dataset
DefaultAssay(object) <- "Spatial.008um"
Idents(object) <- "seurat_cluster.projected"
p2 <- DimPlot(object, reduction = "full.umap.sketch", label = F) + ggtitle("Projected clustering (full dataset)") + theme(legend.position = "bottom")

p1 | p2

当然也可以映射到空间水平：

SpatialDimPlot(object, label = T, repel = T, label.size = 4)

由于亚群数量太多了，在空间水平很难观察到具体某个亚群的分布情况，我们可以借助CellsByIdentities函数进行可视化：

Idents(object) <- "seurat_cluster.projected"
cells <- CellsByIdentities(object, idents = c(0, 4, 32, 34, 35))
p <- SpatialDimPlot(object,
  cells.highlight = cells[setdiff(names(cells), "NA")],
  cols.highlight = c("#FFFF00", "grey50"), facet.highlight = T, combine = T
) + NoLegend()
p

我们还可以使用FindAllMarkers函数进行差异表达分析，用于识别并可视化每个亚群的top基因表达情况，但是有单细胞分析经验的小伙伴应该都知道，单细胞数量多的情况下，FindAllMarkers函数运行非常慢，因此Seurat流程使用了抽样的方式运行FindAllMarkers函数。另外，这里可使用BuildClusterTree函数根据亚群的相似性进行排序。以下是HD差异表达分析的完整代码：

# Crete downsampled object to make visualization either
DefaultAssay(object) <- "Spatial.008um"
Idents(object) <- "seurat_cluster.projected"
object_subset <- subset(object, cells = Cells(object[["Spatial.008um"]]), downsample = 1000)

# Order clusters by similarity
DefaultAssay(object_subset) <- "Spatial.008um"
Idents(object_subset) <- "seurat_cluster.projected"
object_subset <- BuildClusterTree(object_subset, assay = "Spatial.008um", reduction = "full.pca.sketch", reorder = T)

markers <- FindAllMarkers(object_subset, assay = "Spatial.008um", only.pos = TRUE)
markers %>%
  group_by(cluster) %>%
  dplyr::filter(avg_log2FC > 1) %>%
  slice_head(n = 5) %>%
  ungroup() -> top5

object_subset <- ScaleData(object_subset, assay = "Spatial.008um", features = top5$gene)
p <- DoHeatmap(object_subset, assay = "Spatial.008um", features = top5$gene, size = 2.5) + theme(axis.text = element_text(size = 5.5)) + NoLegend()
p

四. 补充分析

这里除了Seurat官网使用的FindAllMarkers函数，我也非常推荐大家使用R包COSG识别亚群的top基因【速度上吊打FindAllMarkers的单细胞亚群特异性高表达基因查询算法】，哪怕是近40W的细胞，COSG的运行速度也是非常快的，非常适合单细胞以及HD分析：

#remotes::install_github(repo = 'genecell/COSGR')
library(COSG)
marker_cosg <- cosg(
  object,
  groups='all',
  assay='RNA',
  slot='data',
  mu=1,
  n_genes_user=200)

然后就可以利用识别到的marker进行手动注释了！虽然手动注释是比较麻烦的，好在各种组织的单细胞图谱当前都已被绘制，所以有大量的单细胞/空转文献可供参考，我们可以从这些文献中获取大量相关的细胞类型标志物。此外，也非常推荐结合自动注释的方案进行自动注释进行交叉验证，例如singleR或者Celltypist【Celltypist：超越singleR的单细胞注释工具】。

识别到的marker除了用于注释以外，还可以用作富集分析【R包BioEnricher：一键式富集分析和可视化】，用于了解感兴趣亚群的功能特征，例如这篇文章的Fig. 4F【https://www.journal-of-hepatology.eu/article/S0168-8278(20)30360-3/fulltext】：