引言
本系列讲解 空间转录组学 (Spatial Transcriptomics) 相关基础知识与数据分析教程，持续更新，欢迎关注，转发，文末有交流群！
降维
下一步的标准做法是使用主成分分析（PCA）进行降维，应用于 top HVGs（或 SVGs）集合。然后，top 主成分（PCs）的集合可以作为后续步骤的输入。
spe <- runPCA(spe, subset_row = sel)
此外，我们可以使用 UMAP 算法进行非线性降维，应用于 top PCs 的集合（默认 50 个）。前两个 UMAP 维度可以通过将它们绘制在 x 轴和 y 轴上用于可视化目的。
spe <- runUMAP(spe, dimred = "PCA")
聚类
接下来，我们在上述处理步骤之后，展示一些下游分析的简短示例。
在此，我们运行一种具有空间意识的聚类算法 BayesSpace，以识别空间域。
BayesSpace 专门为基于测序的 ST 数据开发，并考虑了测量的空间坐标。
# run BayesSpace clustering
.spe <- spatialPreprocess(spe, skip.PCA = TRUE)
.spe <- spatialCluster(.spe, nrep = 1000, burn.in = 100, q = 10, d = 20)
# cluster labels
table(spe$BayesSpace <- factor(.spe$spatial.cluster))
可视化
我们可以在 x-y 空间中把聚类标签可视化为空间域，并与该数据集附带的手工注释参考标签（ground_truth）进行比较。
# using plotting functions from ggspavis package
# and formatting using patchwork package
plotCoords(spe, annotate = "ground_truth") +
plotCoords(spe, annotate = "BayesSpace") +
plot_layout() &
theme(legend.key.size = unit(0, "lines")) &
scale_color_manual(values = unname(pals::trubetskoy()))
table(spe$BayesSpace, spe$ground_truth, useNA = "ifany")
通过查看 top PCs 的空间分布，我们可以观察到表达变异的主要驱动因素是白质（WM）与非 WM 皮质层之间的区别。
pcs <- reducedDim(spe, "PCA")
pcs <- pcs[, seq_len(4)]
lapply(colnames(pcs), \(.) {
spe[[.]] <- pcs[, .]
plotCoords(spe, annotate = .)
}) |>
wrap_plots(nrow = 1) & coord_equal() &
geom_point(shape = 16, stroke = 0, size = 0.2) &
scale_color_gradientn(colors = pals::jet(), n.breaks = 3) &
theme_void() & theme(
plot.title = element_text(hjust = 0.5),
legend.key.width = unit(0.2, "lines"),
legend.key.height = unit(0.8, "lines"))
差异表达
在将 spot 聚类以识别空间域之后，我们可以测试在空间域之间差异表达（DE）的基因；这些基因可视为空间域的标记基因。
在此，我们将使用成对 t 检验，并专门检测上调（而非下调），即所报告的标记基因应在该基因所属的聚类中表达更高。
# using scran package
mgs <- findMarkers(spe, groups = spe$BayesSpace, direction = "up")
top <- lapply(mgs, \(df) rownames(df)[df$Top <= 2])
length(top <- unique(unlist(top)))
我们可以将选定的标记基因可视化为热图（显示聚类水平的平均表达）：
# compute cluster-wise averages
pbs <- aggregateAcrossCells(spe,
ids = spe$BayesSpace, subset.row = top,
use.assay.type = "logcounts", statistics = "mean")
# use gene symbols as feature names
mtx <- t(assay(pbs))
colnames(mtx) <- rowData(pbs)$gene_name
# using pheatmap package
pheatmap(mat = mtx, scale = "column")
或空间图（显示 spot 水平的表达值）：
# gene-wise spatial plots
gs <- c("MBP", "PLP1", "NRGN", "SNAP25", "NEFL", "HPCAL1")
ps <- lapply(gs, \(.) {
plotCoords(spe,
annotate = .,
feature_names = "gene_name",
assay_name = "logcounts") })
# figure arrangement
wrap_plots(ps, nrow = 2) &
theme(legend.key.width = unit(0.4, "lines"),
legend.key.height = unit(0.8, "lines")) &
scale_color_gradientn(colors = rev(hcl.colors(9, "Rocket")))
