引言

本系列讲解 空间转录组学 (Spatial Transcriptomics) 相关基础知识与数据分析教程,持续更新,欢迎关注,转发,文末有交流群!
简介
Visium HD 是由 10x Genomics 开发的下一代空间转录组技术,旨在为整个组织切片提供单细胞分辨率的空间基因表达数据。
自 2024 年起商业化,Visium HD 将该平台的分辨率从 55 µm 点推进到亚细胞分辨率(2 µm)的 bins。
在此工作流程中,我们将使用结直肠癌的 Visium HD 数据集演示常见的分析步骤。
依赖
library(Banksy)
library(BiocParallel)
library(dplyr)
library(DropletUtils)
library(ggplot2)
library(ggspavis)
library(igraph)
library(magick)
library(OSTA.data)
library(patchwork)
library(pheatmap)
library(scater)
library(scran)
library(scuttle)
library(sf)
library(spacexr)
library(SpatialExperiment)
library(SpotSweeper)
library(Statial)
library(tidyr)
library(VisiumIO)
数据集
在本演示中,我们对 16 µm 的 bins 进行去卷积,并将其与 10x Genomics 提供的 8 µm bins 的结果进行一致性比较。
# retrieve dataset from OSF repo
id <- "VisiumHD_HumanColon_Oliveira"
pa <- OSTA.data_load(id)
dir.create(td <- tempfile())
unzip(pa, exdir=td)
# read 8um bins into 'SpatialExperiment'
vhd8 <- TENxVisiumHD(
spacerangerOut=td,
processing="filtered",
format="h5",
images="lowres",
bin_size="008") |>
import()
# use gene symbols as feature names
gs <- rowData(vhd8)$Symbol
rownames(vhd8) <- make.unique(gs)
vhd8

这些数据附带由 spacexr 中实现的 RCTD去卷积估算得到的 bin 级注释。具体而言,提供两组标签,分别对应每个 bin 内 38 种细胞类型中出现频率最高和第二高的类型:
# retrieve annotations
gz <- "binned_outputs/square_008um/deconvolution.csv.gz"
df <- read.csv(file.path(td, gz), row.names=2)
head(df <- df[complete.cases(df), -1])

# keep only annotated bins
vhd8 <- vhd8[, rownames(df)]
names(df) <- c("DeconClass", "DeconLabel1", "DeconLabel2")
colData(vhd8)[names(df)] <- df
lab <- grep("Label", names(cd <- colData(vhd8)), value=TRUE)
pal <- hcl.colors(length(unique(unlist(cd[lab]))), "Spectral")
lapply(lab, \(.) {
plotCoords(vhd8, annotate=., point_size=0.2, point_shape = 15) + ggtitle(.)
}) |>
wrap_plots(nrow=1, guides="collect") &
scale_color_manual(NULL, values=pal) &
theme(legend.key.size=unit(0, "lines"))

具有 16 µm bin 的 VisiumHD 的另一种分辨率可以如下读取:
# read 16um bins into 'SpatialExperiment'
vhd16 <- TENxVisiumHD(
spacerangerOut=td,
processing="filtered",
format="h5",
images="lowres",
bin_size="016") |>
import()
# use symbols as feature names
gs <- rowData(vhd16)$Symbol
rownames(vhd16) <- make.unique(gs)
vhd16

Space Ranger 生成的聚类结果适用于两种分辨率。
csv <- list.files(td, "clustering.csv.gz", recursive=TRUE, full.names=TRUE)
dfs <- lapply(csv, read.csv, row.names="barcode")
colData(vhd8)$cluster <- factor(dfs[[1]][colnames(vhd8), "cluster"])
colData(vhd16)$cluster <- factor(dfs[[2]][colnames(vhd16), "cluster"])
我们可以在 8 µm 和 16 µm 分辨率下可视化给定的聚类结果。稍后,我们将比较 16 µm 处的聚类结果与 Banksy 的聚类结果。
(plotCoords(vhd8,
annotate="cluster", point_size=0.05, point_shape=15,
pal=unname(pals::kelly())) + ggtitle("8 µm")) |
(plotCoords(vhd16[, !is.na(vhd16$cluster)],
annotate="cluster", point_size=0.1, point_shape=15,
pal=unname(pals::kelly())) + ggtitle("16 µm"))

为保持运行时间较短,我们的下游分析将仅在 16 µm bins 的一个子区域(占数据覆盖面积的 1/64)上进行,我们在此读取该子区域;
.rng <- \(spe) {
xy <- spatialCoords(spe)*scaleFactors(spe)
xs <- range(xy[, 1]); ys <- nrow(imgRaster(spe))-range(xy[, 2])
x1 <- xs[1] + 4*(xs[2]-xs[1])/8; x2 <- xs[2] - 3*(xs[2]-xs[1])/8
y1 <- ys[1] + 3*(ys[2]-ys[1])/8; y2 <- ys[2] - 4*(ys[2]-ys[1])/8
list(box=c(x1, x2, y1, y2), cov=c(xs, ys))
}
vhd8r <- .rng(vhd8)
# vhd16r <- .rng(vhd16)
# use 8um bounding boxes to subset spes at both resolutions;
# it is similar to 16um's box range
.box <- \(roi, lty, lwd, col) {
geom_rect(
xmin=vhd8r[[roi]][1], xmax=vhd8r[[roi]][2],
ymin=vhd8r[[roi]][3], ymax=vhd8r[[roi]][4],
col=col, fill=NA, linetype=lty, linewidth=lwd)
}
cov <- .box(roi="cov", lty=2, lwd=1/2, col="grey")
box <- .box(roi="box", lty=4, lwd=2/3, col="black")
# plotting
.lim <- \(spe) list(
xlim(spe[["box"]][c(1, 2)]),
ylim(spe[["box"]][c(4, 3)]))
plotVisium(vhd8, spots=FALSE, point_shape=22) + cov + box +
ggtitle("grey: data coverage\n black: zoomed region") +
plotVisium(vhd8, point_size=0.8, zoom=TRUE, point_shape=22) + .lim(vhd8r) +
ggtitle("8 µm bins in black box") +
plotVisium(vhd16, point_size=1.6, zoom=TRUE, point_shape=22) + .lim(vhd8r) +
ggtitle("16 µm bins in black box") +
plot_layout(nrow=1, widths=c(1.5, 1, 1)) & facet_null() &
theme(plot.title=element_text(hjust=0.5, vjust=0.5))
