空间转录组: 降维聚类+差异分析

引言

本系列讲解 空间转录组学 (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")))