空间转录组: 反卷积

引言

本系列讲解 空间转录组学 (Spatial Transcriptomics) 相关基础知识与数据分析教程，持续更新，欢迎关注，转发，文末有交流群！

简介

基于测序的空间转录组（ST）数据在每个 spot 中可能包含 0 到多个细胞，这些细胞可能完全被 spot 覆盖，也可能只是部分被覆盖，具体取决于平台的空间分辨率以及组织细胞的密度。数据的这一特点意味着一个 spot 内可能存在细胞类型的混合，因此也会出现转录程序的混合。

上图展示了 10x Genomics Visium spot（直径为 55 µm；紫色线）叠加在苏木精-伊红（H&E）染色图像上的情况。spot 之间的中心距为 100 µm（黄线），典型免疫细胞的直径约为 10 µm（青线）。

为了帮助理解这些混合物，目前已有至少 20 种反卷积技术被提出，用于 spot 级别的 ST 数据。某些方法需要从单细胞 RNA 测序（scRNA-seq）参考数据集中借用信息，而另一些方法则无需参考。根据底层算法将这些方法分为五大类：

• probabilistic-based：使用贝叶斯推断、似然估计或概率建模来估算细胞类型组成，同时纳入不确定性。可用工具包括 R 中的 spacexr（即 RCTD）、CARD、SpatialDecon 和 STdeconvolve，以及 Python 中的 std-poisson、stereoscope、scvi-tools（即 DestVI）、STRIDE 和 cell2location。
• non-negative matrix factorization (NMF)-based：将基因表达数据分解为代表不同细胞类型的潜在成分。可用工具包括 R 中的 Giotto（即 SpatialDWLS）和 SPOTlight，以及 Python 中的 NMFReg。
• graph-based：使用图神经网络或基于图的优化来建模空间关系。可用工具包括 R 中的 SD2，以及 Python 中的 DSTG 和 SpiceMix。
• optimal transport (OT)-based：通过将 scRNA-seq 与 ST 数据进行映射，推断空间基因表达分布。可用工具包括 Python 中的 SpaOTsc 和 novosparc。
• deep learning-based：利用神经网络对齐并整合单细胞与空间转录组数据。例如 Python 中的 Tangram。

其中，std-poisson、STdeconvolve 和 SpiceMix 为无需参考的方法。纳入空间位置信息的方法有 CARD、DSTG、SD2、Tangram、cell2location、DestVI、std-poisson 和 SpiceMix。

本文将展示在每个 spot 上进行细胞类型反卷积，使用 RCTD 在 Visium 和 VisiumHD 数据集上的应用。

依赖

library(OSTA.data)
library(VisiumIO)
library(SpatialExperiment)
library(ggspavis)
library(spacexr)
library(CARDspa)
library(ggplot2)
library(patchwork)
library(DropletUtils)
library(BiocParallel)
library(scran)
library(scater)
library(pheatmap)

在 Visium 乳腺癌数据中，我们在没有单细胞（Chromium）参考的情况下进行细胞类型去卷积，并将结果与 10x Genomics 提供的 Visium 注释进行一致性比较。

# retrieve dataset from OSF repository
id <- "Visium_HumanBreast_Janesick"
pa <- OSTA.data_load(id)
dir.create(td <- tempfile())
unzip(pa, exdir = td)

# read into 'SpatialExperiment'
vis <- TENxVisium(
    spacerangerOut = td, 
    processing = "filtered", 
    format = "h5", 
    images = "lowres") |> 
    import()

# retrieve spot annotations & add as metadata
df <- read.csv(file.path(td, "annotation.csv"))
cs <- match(colnames(vis), df$Barcode)
vis$anno <- factor(df$Annotation[cs])

# set gene symbols as feature names
rownames(vis) <- make.unique(rowData(vis)$Symbol)
vis

xy <- spatialCoords(vis) * scaleFactors(vis)
xs <- range(xy[, 1])
ys <- nrow(imgRaster(vis)) - range(xy[, 2])
box <- geom_rect(
    xmin = xs[1], xmax = xs[2], ymin = ys[1], ymax = ys[2], 
    col = "black", fill = NA, linetype = 2, linewidth = 2/3)
plotVisium(vis, spots = FALSE, point_size = 1) + box + 
    plotVisium(vis, point_size = 1, zoom = TRUE) + 
    plot_layout(nrow = 1) & facet_null()

反卷积在质量控制之后执行，通常使用未归一化且未转换的（即原始的）计数矩阵。

此处，我们快速检查一些典型的 spot 级别指标。

sub <- list(mt = grep("^MT-", rownames(vis)))
vis <- addPerCellQCMetrics(vis, subsets = sub)

vis$log_sum <- log1p(vis$sum)
plotCoords(vis, annotate = "log_sum") + ggtitle("log library size") + 
    plotCoords(vis, annotate = "subsets_mt_percent") + ggtitle("% mitochondrial") + 
    plot_layout(nrow = 1) & theme(
      legend.key.width = unit(0.5, "lines"), 
      legend.key.height = unit(1, "lines")) & 
    scale_color_gradientn(colors = pals::jet())

ggplot(data.frame(colData(vis))) + 
  geom_point(aes(x=sum, y=subsets_mt_percent)) + 
  scale_x_log10() + 
  scale_y_sqrt()

一些spots的library sizes较低，可以删除。

vis <- vis[,vis$sum > 1000]

我们首先可视化10x Genomics提供的spot-level 细胞类型注释。

plotCoords(vis, 
           annotate = "anno", point_size = 1, 
           pal = unname(pals::trubetskoy())) + 
    theme(legend.key.size = unit(0, "lines"))

现在，我们为 Visium 数据集加载单细胞（Chromium）参考数据。为了简化演示，我们将 10x Genomics 提供的若干细胞类型注释（即 Annotation）合并为更宽泛的类别（即 Annogrp）。

# retrieve dataset from OSF repo
id <- "Chromium_HumanBreast_Janesick"
pa <- OSTA.data_load(id)
dir.create(td <- tempfile())
unzip(pa, exdir = td)

# read into 'SingleCellExperiment'
fs <- list.files(td, full.names = TRUE)
h5 <- grep("h5$", fs, value = TRUE)
sce <- read10xCounts(h5, col.names = TRUE)

# use gene symbols as feature names
rownames(sce) <- make.unique(rowData(sce)$Symbol)

# retrieve cell type labels
csv <- grep("csv$", fs, value = TRUE)
cd <- read.csv(csv, row.names = 1)

# ignore mixtures
lab <- cd$Annotation
lab[grepl("Hyb", lab)] <- NA 

# simplify annotations
pat <- c(
    "B Cell" = "B", "T Cell" = "T", "Mac" = "macro", 
    "Mast" = "mast", "DCs" = "dendritic", 
    "Peri" = "perivas", "End" = "endo", "Str" = "stromal", 
    "Inv" = "tumor", "Myo" = "myoepi")
for (. in names(pat)) 
    lab[grep(., lab)] <- pat[.]
lab <- gsub("\\s", "", lab)

# add as cell metadata
table(cd$Annogrp <- lab)

colData(sce)[names(cd)] <- cd[colnames(sce), ]

我们只保留带有注释且未被标记为 “Hybrid” 的 Chromium 数据，因为这些对应的是混合亚群。

sce <- sce[, !is.na(sce$Annogrp)]
dim(sce)

##  [1] 18082 26031