各位小伙伴，好久不见！最近我忙完了博士毕业相关的事，离开武汉刚回到杭州，六年的硕博生活终于结束啦！还在坑里的朋友们要坚持住！
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数据和单细胞数据。
- 整合单细胞和空转数据多种方法之CellTrek
- 整合单细胞和空转数据多种方法之Cell2location
- 整合单细胞和空转数据多种方法之RCTD
- 整合单细胞和空转数据多种方法之R包semla
上述分析算法其实并不局限于HD这个技术。可以说单细胞数据分析、传统的Visium空转分析、HD空转分析甚至其他平台的空转，都有非常多交叉的内容，我们在学习的过程中需要融会贯通，举一反三。
一. 安装/更新Seurat v5
这里使用的是Seurat v5版本：
代码语言：javascript复制#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"指定默认使用的插槽表达矩阵。
检查一下加载好的数据：
代码语言：javascript复制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，其他指标和华大的过滤标准一致。”
在这里我贴一下我质控过滤的代码，仅供参考：
代码语言：javascript复制## 华大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的单细胞）：
代码语言：javascript复制# 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
当然也可以映射到空间水平：
代码语言：javascript复制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】：
学习一下如何解读和描述单细胞/空转亚群富集分析的结果：
本期Seurat官网的HD标准分析流程介绍到此就结束了！下期推文我介绍一下HD空转的空间组织域识别和反卷积算法，这里Seurat使用的是BANKSY算法以及RCTD算法。
