空间单细胞｜基于图像的空间数据分析（2）
## 引言在这篇[指南](https://satijalab.org/seurat/articles/seurat5_spatial_vignette_2 "Source")中，我们介绍了Seurat的一个新扩展功能，用以分析新型的空间解析数据，将重点介绍由不同成像技术生成的三个公开数据集。
- Vizgen MERSCOPE（用于小鼠大脑研究）
- Nanostring CosMx空间分子成像仪（用于FFPE人类肺组织）
- Akoya CODEX（用于人类淋巴结研究）
## 小鼠大脑：10x Genomics Xenium In Situ
在本节中，我们将对Xenium平台生成的数据进行分析。本指南展示了如何导入Xenium平台输出的每转录本位置信息、细胞与基因的表达矩阵、细胞边界分割以及细胞中心点数据。得到的Seurat对象将包含每个细胞的基因表达概况、细胞的中心点和边界位置，以及每个被检测到的转录本的具体位置。这些细胞层面的基因表达概况与标准的单细胞RNA测序数据相似，可以运用相同的分析工具进行处理。
这里使用的是10x Genomics公司为Xenium Explorer Demo提供的“Tiny”子集数据集，该数据集源自“Fresh Frozen Mouse Brain”，下载方式如下。这些分析步骤同样适用于更大的“完整冠状切面”数据集，但处理时间会更长。
```sh
wget https://cf.10xgenomics.com/samples/xenium/1.0.2/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP_outs.zip
unzip Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP_outs.zip
```
首先，我们读取数据集并创建一个 Seurat 对象。提供 Xenium 运行的数据文件夹的路径作为输入路径。 RNA 数据存储在 Seurat 对象的 Xenium 分析中。此步骤大约需要一分钟。
```R
path <- "/brahms/hartmana/vignette_data/xenium_tiny_subset"
# Load the Xenium data
xenium.obj <- LoadXenium(path, fov = "fov")
# remove cells with 0 counts
xenium.obj <- subset(xenium.obj, subset = nCount_Xenium > 0)
```
空间数据被嵌入到Seurat对象的特定字段中，这些字段以加载的“视野”（Field of View，简称FOV）名称来标识。一开始，所有的数据都被加载到默认的FOV中，名为fov。接下来，我们将创建一个新的视野，这个视野将聚焦于我们特别关注的区域。
Xenium分析工具提供了Seurat的标准质量控制图表。例如，我们可以通过小提琴图来查看每个细胞中的基因数量（nFeature_Xenium）和每个细胞中的转录本计数（nCount_Xenium）。这些图表帮助我们对数据集的质量有一个直观的了解。
```R
VlnPlot(xenium.obj, features = c("nFeature_Xenium", "nCount_Xenium"), ncol = 2, pt.size = 0)
```
接下来，我们使用 ImageDimPlot() 在组织上绘制全抑制神经元标记 Gad1、抑制神经元亚型标记 Pvalb 和 Sst 以及星形胶质细胞标记 Gfap 的位置。
```R
ImageDimPlot(xenium.obj, fov = "fov", molecules = c("Gad1", "Sst", "Pvalb", "Gfap"), nmols = 20000)
```
在这里，我们使用 ImageFeaturePlot() 在每个细胞水平上可视化一些关键层标记基因的表达水平，这类似于用于可视化 2D 嵌入表达的 FeaturePlot() 函数。我们手动将每个基因的 max.cutoff 调整到其计数分布的大约第 90 个百分位数（可以使用 max.cutoff='q90' 指定）以提高对比度。
```R
ImageFeaturePlot(xenium.obj, features = c("Cux2", "Rorb", "Bcl11b", "Foxp2"), max.cutoff = c(25,
35, 12, 10), size = 0.75, cols = c("white", "red"))
```
我们可以使用 Crop() 函数放大所选区域。放大后，我们可以可视化细胞分割边界以及单个分子。
```R
cropped.coords <- Crop(xenium.obj[["fov"]], x = c(1200, 2900), y = c(3750, 4550), coords = "plot")
xenium.obj[["zoom"]] <- cropped.coords
# visualize cropped area with cell segmentations & selected molecules
DefaultBoundary(xenium.obj[["zoom"]]) <- "segmentation"
ImageDimPlot(xenium.obj, fov = "zoom", axes = TRUE, border.color = "white", border.size = 0.1, cols = "polychrome",
coord.fixed = FALSE, molecules = c("Gad1", "Sst", "Npy2r", "Pvalb", "Nrn1"), nmols = 10000)
```
接下来，我们使用 SCTransform 进行归一化，然后进行标准降维和聚类。此步骤从开始到结束大约需要 5 分钟。
```R
xenium.obj <- SCTransform(xenium.obj, assay = "Xenium")
xenium.obj <- RunPCA(xenium.obj, npcs = 30, features = rownames(xenium.obj))
xenium.obj <- RunUMAP(xenium.obj, dims = 1:30)
xenium.obj <- FindNeighbors(xenium.obj, reduction = "pca", dims = 1:30)
xenium.obj <- FindClusters(xenium.obj, resolution = 0.3)
```
然后，我们可以使用 DimPlot() 在 UMAP 空间中根据每个细胞的聚类对每个单元格进行着色，或者使用 ImageDimPlot() 覆盖在图像上，从而可视化聚类结果。
```R
DimPlot(xenium.obj)
```
我们可以在 UMAP 坐标上可视化之前查看的标记的表达水平。
```R
FeaturePlot(xenium.obj, features = c("Cux2", "Bcl11b", "Foxp2", "Gad1", "Sst", "Gfap"))
```
现在，我们可以使用 ImageDimPlot() 为上一步中确定的簇标签着色的细胞位置着色。
```R
ImageDimPlot(xenium.obj, cols = "polychrome", size = 0.75)
```
利用每个细胞的定位信息，我们能够计算出它们各自的空间生态位。为了对细胞进行注释，我们采用了Allen Brain Institute提供的大脑皮层参考数据，因此首先需要将整个数据集的范围限定在大脑皮层。可以通过提供的链接安装Allen Brain的参考数据。
在下面的分析中，我们利用Slc17a7基因的表达情况来辅助识别大脑皮层的具体区域。
```R
xenium.obj <- LoadXenium("/brahms/hartmana/vignette_data/xenium_tiny_subset")
p1 <- ImageFeaturePlot(xenium.obj, features = "Slc17a7", axes = TRUE, max.cutoff = "q90")
p1
```
```R
crop <- Crop(xenium.obj[["fov"]], x = c(600, 2100), y = c(900, 4700))
xenium.obj[["crop"]] <- crop
p2 <- ImageFeaturePlot(xenium.obj, fov = "crop", features = "Slc17a7", size = 1, axes = TRUE, max.cutoff = "q90")
p2
```
`FindTransferAnchors`函数能够用来整合空间转录组数据集中的单个斑点数据。而Seurat v5版本进一步支持了一种名为Robust Cell Type Decomposition（稳健细胞类型分解，简称RCTD）的计算方法，该方法能够在提供单细胞RNA测序（scRNA-seq）参考数据的基础上，对空间数据集中的斑点数据进行解卷积分析。RCTD已被证实能够有效地对来自SLIDE-seq、Visium以及10x Genomics公司的Xenium原位空间平台等多种技术的空间数据进行注释。
为了执行RCTD分析，我们首先需要从GitHub上安装名为spacexr的R包，该包提供了RCTD功能的实现。
```R
devtools::install_github("dmcable/spacexr", build_vignettes = FALSE)
```
从 Seurat 查询和参考对象中提取计数、聚类和点信息，以构建 RCTD 用于注释的参考和 SpatialRNA 对象。然后将注释的输出添加到 Seurat 对象。
```R
library(spacexr)
query.counts <- GetAssayData(xenium.obj, assay = "Xenium", slot = "counts")[, Cells(xenium.obj[["crop"]])]
coords <- GetTissueCoordinates(xenium.obj[["crop"]], which = "centroids")
rownames(coords) <- coords$cell
coords$cell <- NULL
query <- SpatialRNA(coords, query.counts, colSums(query.counts))
# allen.corted.ref can be downloaded here:
# https://www.dropbox.com/s/cuowvm4vrf65pvq/allen_cortex.rds?dl=1
allen.cortex.ref <- readRDS("/brahms/shared/vignette-data/allen_cortex.rds")
allen.cortex.ref <- UpdateSeuratObject(allen.cortex.ref)
Idents(allen.cortex.ref) <- "subclass"
# remove CR cells because there aren't enough of them for annotation
allen.cortex.ref <- subset(allen.cortex.ref, subset = subclass != "CR")
counts <- GetAssayData(allen.cortex.ref, assay = "RNA", slot = "counts")
cluster <- as.factor(allen.cortex.ref$subclass)
names(cluster) <- colnames(allen.cortex.ref)
nUMI <- allen.cortex.ref$nCount_RNA
names(nUMI) <- colnames(allen.cortex.ref)
nUMI <- colSums(counts)
levels(cluster) <- gsub("/", "-", levels(cluster))
reference <- Reference(counts, cluster, nUMI)
# run RCTD with many cores
RCTD <- create.RCTD(query, reference, max_cores = 8)
RCTD <- run.RCTD(RCTD, doublet_mode = "doublet")
annotations.df <- RCTD@results$results_df
annotations <- annotations.df$first_type
names(annotations) <- rownames(annotations.df)
xenium.obj$predicted.celltype <- annotations
keep.cells <- Cells(xenium.obj)[!is.na(xenium.obj$predicted.celltype)]
xenium.obj <- subset(xenium.obj, cells = keep.cells)
```
与以往的独立细胞分析不同，空间数据让我们能够更全面地理解细胞，不仅考虑它们周围的局部环境，还包括它们在整个空间中的背景。在Seurat v5版本中，我们新增了对空间数据进行“生态位”分析的功能，这种分析可以识别出组织中的特定区域（即“生态位”），每个区域都有其独特的空间邻近细胞类型组合。这一方法受到了Goltsev等人在2018年发表于《细胞》杂志和He等人在2022年发表于《自然生物技术》杂志的研究方法的启发，我们为每个细胞定义了一个“局部邻域”，这包括了与其在空间上距离最近的k个邻居，并统计这个邻域内每种细胞类型的频率。接着，我们应用k均值聚类算法，将具有相似邻域特征的细胞归类到相同的空间生态位中。
在Seurat中，我们通过调用BuildNicheAssay函数来创建一个新的分析模块，称为“生态位”，它包含了每个细胞周围空间的细胞类型组成信息。此外，该函数还返回了一个元数据列“niches”，该列展示了基于生态位分析的细胞聚类结果。
```R
xenium.obj <- BuildNicheAssay(object = xenium.obj, fov = "crop", group.by = "predicted.celltype",
niches.k = 5, neighbors.k = 30)
```
然后，我们可以根据细胞类型身份或利基身份对细胞进行分组。确定的生态位清楚地划分了皮质中的神经元层。
```R
celltype.plot <- ImageDimPlot(xenium.obj, group.by = "predicted.celltype", size = 1.5, cols = "polychrome",
dark.background = F) + ggtitle("Cell type")
niche.plot <- ImageDimPlot(xenium.obj, group.by = "niches", size = 1.5, dark.background = F) + ggtitle("Niches") +
scale_fill_manual(values = c("#442288", "#6CA2EA", "#B5D33D", "#FED23F", "#EB7D5B"))
celltype.plot | niche.plot
```
此外，我们观察到每个生态位的组成对于不同的细胞类型都是丰富的。
```R
table(xenium.obj$predicted.celltype, xenium.obj$niches)
##
## 1 2 3 4 5
## Astro 115484241146 89
## Endo 27250 66 62 45
## L2-3 IT 0 1749 14 5 6
## L4 2 2163 3 94 14
## L5 IT 2 66 2627 84
## L5 PT 2 92 4711 21
## L6 CT 82 2 0 34857
## L6 IT 4 22 1 56275
## L6b 95 0 0 2 2
## Lamp5 4 77 61 8 7
## Macrophage 7 48 15 16 10
## Meis2 0 0 0 0 0
## NP 0 1 0 78 9
## Oligo 305130 42 76 70
## Peri 6 40 4 10 12
## Pvalb 5155 0 75 54
## Serpinf1 0 5 0 1 1
## SMC 2 34 2 12 2
## Sncg 3 15 1 0 5
## Sst 0 53 0 55 26
## Vip 3 85 5 14 4
## VLMC 2 31257 7 2
```
页:
[1]