Seurat4.0系列教程22：空间转录组的分析、可视化与整合
Seurat4.0系列教程告一段落，但这决不是终点。这个系列教程是给大家打开一扇窗，让大家知道Seurat4.0有这些功能可用，少走弯路。后续还要自己深入研究。学，然后知不足。加油吧！少年！不足之处，恳请大家批评指正！
概述本教程演示了如何使用 Seurat 4.0 分析单细胞空间转录组数据。虽然分析流程类似于单细胞RNA-seq分析的Seurat工作流程，但我们引入了更新的交互和可视化工具，特别强调了空间和分子信息的整合。本教程将涵盖以下任务：
标准化降维和聚类检测空间高变基因交互式可视化与单细胞RNA-seq数据整合处理多个切片对于第一个教程，我们分析使用Visium 技术从 10 x genomics生成的数据集。我们将在不远的将来扩展seurat，以处理其他数据类型，包括SLIDE-Seq、STARmap和MERFISH。
首先，我们加载Seurat和所需的其他包。
library(Seurat)library(SeuratData)library(ggplot2)library(patchwork)library(dplyr)这里，我们将使用最近发布的使用Visium v1 chemistry产生的小鼠大脑切片的数据集。
可以在此处下载数据，然后使用Load10X_Spatial()将其加载到 Seurat 中。这一步读取spaceranger流程的输出，并返回一个包含 spot-level表达数据以及组织切片相关图像的 Seurat 对象。可以使用SeuratData 包轻松访问数据，如下所示。
InstallData("stxBrain")brain <- LoadData("stxBrain", type="anterior1")通过基因表达数据执行的预处理步骤类似于传统的 scRNA-seq 。首先需要使数据标准化，以便较少数据之间测序深度的差异。分子计数/点的差异对于空间数据集来说可能很大，特别是如果整个组织的细胞密度存在差异。在这里看到巨大的异质性，这需要有效的标准化。
plot1 <- VlnPlot(brain, features="nCount_Spatial", pt.size=0.1) + NoLegend()plot2 <- SpatialFeaturePlot(brain, features="nCount_Spatial") + theme(legend.position="right")wrap_plots(plot1, plot2)这些图表明，跨点分子计数的差异不仅在技术上，而且依赖于组织解剖学。例如，神经元（如皮质白质）缺失的组织区域，可反复显示较低的分子计数。因此，标准化方法（如函数LogNormalize()）可能会有问题，这些方法使每个数据点在标准化后具有相同的基础"大小"。
作为替代方案，我们建议使用 sctransform，该模型构建了基因表达的常规负二元模型，以便在保持生物差异的同时考虑技术。sctransform 使数据标准化，检测高变异基因，并将数据存储在SCTassay.
brain <- SCTransform(brain, assay="Spatial", verbose=FALSE)在 Seurat 中，我们具有探索空间数据固有的视觉性质并与之交互的功能。SpatialFeaturePlot()是Seurat 的FeaturePlot()扩展，并且可以在组织学上覆盖分子数据。例如，在小鼠大脑的这一数据集中，Hpca基因是一个强大的海马标记，Ttr是choroid plexus的标记。
SpatialFeaturePlot(brain, features=c("Hpca", "Ttr"))Seurat 中的默认参数强调分子数据的可视化。但是，您也可以通过更改以下参数来调整点的大小（及其透明度），以改善组织学图像的可视化：
pt.size.factor-这将扩展点的大小。默认值为1.6alpha-最低和最大的透明度。默认值为 c（1，1）。尝试设置为 c（0.1， 1），以降低低表达点的透明度alphap1 <- SpatialFeaturePlot(brain, features="Ttr", pt.size.factor=1)p2 <- SpatialFeaturePlot(brain, features="Ttr", alpha=c(0.1, 1))p1 + p2然后，我们可以继续运行 RNA 表达数据上的降维和聚类，使用与用于 scRNA-seq 分析相同的工作流。
brain <- RunPCA(brain, assay="SCT", verbose=FALSE)brain <- FindNeighbors(brain, reduction="pca", dims=1:30)brain <- FindClusters(brain, verbose=FALSE)brain <- RunUMAP(brain, reduction="pca", dims=1:30)然后，我们可以使用DimPlot()在 UMAP 空间中或使用SpatialDimPlot()可视化聚类的结果。
p1 <- DimPlot(brain, reduction="umap", label=TRUE)p2 <- SpatialDimPlot(brain, label=TRUE, label.size=3)p1 + p2由于有许多颜色，因此可视化哪个 voxel 属于哪个聚类可能具有挑战性。我们有几个策略来帮助解决此问题。设置label参数时，每个集群的中位数放置一个彩色框（参见上图）。
您也可以使用cells.highlight参数在SpatialDimPlot() 。这对于区分单个群的局部空间非常有用，如下所示：
SpatialDimPlot(brain, cells.highlight=CellsByIdentities(object=brain, idents=c(2, 1, 4, 3, 5, 8)), facet.highlight=TRUE, ncol=3)我们还建立了一些交互式绘图功能。SpatialDimPlot()和SpatialDimPlot()两者现在都有一个参数，当设置interactive为TRUE，将打开Rstudio窗格与shiny互动。下面的示例演示了一个交互式，您可以在点上悬停并查看细胞名称和当前marker。
SpatialDimPlot(brain, interactive=TRUE)对于SpatialFeaturePlot() ，设置交互式为TRUE，显示一个交互式窗格，您可以在其中调整点的透明度、点大小以及正在绘制的函数。在浏览数据后，选择已完成的按钮将返回最后一个活动图作为 ggplot 对象。
SpatialFeaturePlot(brain, features="Ttr", interactive=TRUE)LinkedDimPlot()函数将 UMAP 与组织图像连接，并允许进行交互式选择。例如，您可以选择 UMAP 绘图中的区域，并将突出显示图像表示中的相应点。
LinkedDimPlot(brain)Seurat 提供两个工作流程来识别与组织内空间位置相关的分子特征。首先是根据组织内预先注释的解剖区域进行差异表达分析，这些区域可能由不受监督的聚类或先验的知识确定。在这种情况下，此策略将起作用，因为上述集群表现出明确的空间限制。
de_markers <- FindMarkers(brain, ident.1=5, ident.2=6)SpatialFeaturePlot(object=brain, features=rownames(de_markers)[1:3], alpha=c(0.1, 1), ncol=3)另一种方法，是FindSpatiallyVariables()，是寻找空间模式中没有预先注释的基因。默认方法是method='markvariogram`。
brain <- FindSpatiallyVariableFeatures(brain, assay="SCT", features=VariableFeatures(brain)[1:1000], selection.method="markvariogram")现在，我们可视化前 6 个基因的表达。
top.features <- head(SpatiallyVariableFeatures(brain, selection.method="markvariogram"), 6)SpatialFeaturePlot(brain, features=top.features, ncol=3, alpha=c(0.1, 1))与single-cell对象一样，可以对象取子集。在这里，我们取出前额皮层。此过程还有助于将这些数据与下一节中的外皮层 scRNA-seq 数据集整合。首先，我们取群的子集，然后根据确切的位置进一步细分。subset后，我们可以在完整图像或裁剪图像上可视化皮质细胞。
cortex <- subset(brain, idents=c(1, 2, 3, 4, 6, 7))# now remove additional cells, use SpatialDimPlots to visualize what to remove# SpatialDimPlot(cortex,cells.highlight=WhichCells(cortex, expression=image_imagerow > 400 |# image_imagecol < 150))cortex <- subset(cortex, anterior1_imagerow > 400 | anterior1_imagecol < 150, invert=TRUE)cortex <- subset(cortex, anterior1_imagerow > 275 & anterior1_imagecol > 370, invert=TRUE)cortex <- subset(cortex, anterior1_imagerow > 250 & anterior1_imagecol > 440, invert=TRUE)p1 <- SpatialDimPlot(cortex, crop=TRUE, label=TRUE)p2 <- SpatialDimPlot(cortex, crop=FALSE, label=TRUE, pt.size.factor=1, label.size=3)p1 + p2在~50um时， visium assay检测的点将包含多个细胞的表达特征。随着越来越多的可用 scRNA-seq 数据，用户可能有兴趣对每个空间 voxels "去卷积"来预测细胞类型的潜在组成。在准备此教程时，我们测试了各种去卷积和整合方法，使用了由 Allen 研究所生成的 14，000 个成年鼠皮质细胞的scRNA-seq 数据集。我们发现使用整合方法（而不是去卷积方法）更优越，可能是因为空间和单细胞数据集的噪声模型存在较大差异，整合方法针对这些差异进行了稳健的设计。因此，我们应用 Seurat v3中引入的基于"锚点"的sctransform 整合工作流，从而能够将细胞从参考集转到查询集进行注释。
首先加载数据（此处下载），预处理 scRNA-seq 参考集，然后执行标签转移。该过程为每个点输出每个scRNA-seq细胞的概率分类。我们将这些预测添加到新的Seurat 对象进行分析。
allen_reference <- readRDS("../data/allen_cortex.rds")# note that setting ncells=3000 normalizes the full dataset but learns noise models on 3k cells# this speeds up SCTransform dramatically with no loss in performancelibrary(dplyr)allen_reference <- SCTransform(allen_reference, ncells=3000, verbose=FALSE) %>% RunPCA(verbose=FALSE) %>% RunUMAP(dims=1:30)# After subsetting, we renormalize cortexcortex <- SCTransform(cortex, assay="Spatial", verbose=FALSE) %>% RunPCA(verbose=FALSE)# the annotation is stored in the 'subclass' column of object metadataDimPlot(allen_reference, group.by="subclass", label=TRUE)anchors <- FindTransferAnchors(reference=allen_reference, query=cortex, normalization.method="SCT")predictions.assay <- TransferData(anchorset=anchors, refdata=allen_reference$subclass, prediction.assay=TRUE, weight.reduction=cortex[["pca"]], dims=1:30)cortex[["predictions"]] <- predictions.assay现在，我们得到每个类每个点的预测分数。最感兴趣的是frontal cortex region的laminar excitatory neurons。在这里，我们可以区分这些神经元亚型的不同顺序层，例如：
DefaultAssay(cortex) <- "predictions"SpatialFeaturePlot(cortex, features=c("L2/3 IT", "L4"), pt.size.factor=1.6, ncol=2, crop=TRUE)基于这些预测分数，我们还可以预测其位置受空间限制的细胞类型。使用基于标记点过程的相同方法来定义空间变异基因，但使用细胞类型预测分数作为"标记"，而不是基因表达。
cortex <- FindSpatiallyVariableFeatures(cortex, assay="predictions", selection.method="markvariogram", features=rownames(cortex), r.metric=5, slot="data")top.clusters <- head(SpatiallyVariableFeatures(cortex), 4)SpatialPlot(object=cortex, features=top.clusters, ncol=2)最后，我们表明，我们的整合程序能够恢复神经元和非神经元子集已知的空间定位模式，包括laminar excitatory, layer-1 astrocytes, and the cortical grey matter。
SpatialFeaturePlot(cortex, features=c("Astro", "L2/3 IT", "L4", "L5 PT", "L5 IT", "L6 CT", "L6 IT", "L6b", "Oligo"), pt.size.factor=1, ncol=2, crop=FALSE, alpha=c(0.1, 1))鼠大脑的这个数据集包含另一个与大脑另一半相对应的切片。在这里，我们读入它，并执行相同的初始标准化。
brain2 <- LoadData("stxBrain", type="posterior1")brain2 <- SCTransform(brain2, assay="Spatial", verbose=FALSE)为了处理同一 Seurat 对象中的多个切片，我们提供了merge。
brain.merge <- merge(brain, brain2)这便能在下面的RNA表达数据上实现整合降维和聚类。
DefaultAssay(brain.merge) <- "SCT"VariableFeatures(brain.merge) <- c(VariableFeatures(brain), VariableFeatures(brain2))brain.merge <- RunPCA(brain.merge, verbose=FALSE)brain.merge <- FindNeighbors(brain.merge, dims=1:30)brain.merge <- FindClusters(brain.merge, verbose=FALSE)brain.merge <- RunUMAP(brain.merge, dims=1:30)最后，数据可以在单个 UMAP 绘图中共同可视化。 默认情况下，将所有切片绘制为列，分组/基因作为行。
DimPlot(brain.merge, reduction="umap", group.by=c("ident", "orig.ident"))SpatialDimPlot(brain.merge)SpatialFeaturePlot(brain.merge, features=c("Hpca", "Plp1"))这里，我们将分析使用小鼠海马的Slide-seq v2数据集。此教程将遵循与 10x Visium 数据的空间教程大致相同的结构。
您可以使用SeuratData 包轻松访问数据，如下所示。安装数据集后，您可以键入?ssHippo以查看用于创建 Seurat 对象的命令。
InstallData("ssHippo")slide.seq <- LoadData("ssHippo")通过基因表达数据对bead 的初始预处理步骤类似于其他空间seurat分析和传统的scRNA-seq实验。在这里，我们注意到许多beads 包含特别低的 UMI 计数，我们选择保留所有检测到的beads 进行下游分析。
plot1 <- VlnPlot(slide.seq, features="nCount_Spatial", pt.size=0, log=TRUE) + NoLegend()slide.seq$log_nCount_Spatial <- log(slide.seq$nCount_Spatial)plot2 <- SpatialFeaturePlot(slide.seq, features="log_nCount_Spatial") + theme(legend.position="right")wrap_plots(plot1, plot2)然后，我们使用sctransform进行数据标准化，并执行标准的 scRNA-seq 降维和聚类工作流。
slide.seq <- SCTransform(slide.seq, assay="Spatial", ncells=3000, verbose=FALSE)slide.seq <- RunPCA(slide.seq)slide.seq <- RunUMAP(slide.seq, dims=1:30)slide.seq <- FindNeighbors(slide.seq, dims=1:30)slide.seq <- FindClusters(slide.seq, resolution=0.3, verbose=FALSE)然后，我们可以在 UMAP space或bead coordinate space中可视化聚类的结果。
plot1 <- DimPlot(slide.seq, reduction="umap", label=TRUE)plot2 <- SpatialDimPlot(slide.seq, stroke=0)plot1 + plot2SpatialDimPlot(slide.seq, cells.highlight=CellsByIdentities(object=slide.seq, idents=c(1, 6, 13)), facet.highlight=TRUE)为了促进Slide-seq数据集的细胞类型注释，我们利用了由Saunders *，Macosko *等人制作的现有小鼠单细胞RNA-seq海马数据集。数据可在此处作为已处理的Seurat对象下载，而原始计数矩阵可在DropViz网站上获得。
ref <- readRDS("../data/mouse_hippocampus_reference.rds")本文的原始注释在Seurat对象的metadata 中提供。这些注释以几种“resolutions”提供，从大类（ref$class）到亚群ref$subcluster。出于本教程的目的，我们将对细胞类型注释（ref$celltype）进行修改，使之达到良好的平衡。
我们将从运行Seurat标签转移方法开始，以预测每个bead的主要细胞类型。
anchors <- FindTransferAnchors(reference=ref, query=slide.seq, normalization.method="SCT", npcs=50)predictions.assay <- TransferData(anchorset=anchors, refdata=ref$celltype, prediction.assay=TRUE, weight.reduction=slide.seq[["pca"]], dims=1:50)slide.seq[["predictions"]] <- predictions.assay然后，我们可以可视化一些主要类型的预测分数。
DefaultAssay(slide.seq) <- "predictions"SpatialFeaturePlot(slide.seq, features=c("Dentate Principal cells", "CA3 Principal cells", "Entorhinal cortex", "Endothelial tip", "Ependymal", "Oligodendrocyte"), alpha=c(0.1, 1))slide.seq$predicted.id <- GetTransferPredictions(slide.seq)Idents(slide.seq) <- "predicted.id"SpatialDimPlot(slide.seq, cells.highlight=CellsByIdentities(object=slide.seq, idents=c("CA3 Principal cells", "Dentate Principal cells", "Endothelial tip")), facet.highlight=TRUE)正如Visium教程中提到的那样，我们可以通过两种通用方法来识别空间变异基因：预先标注的解剖区域之间的差异表达检测或基因特征对空间位置的依赖性的统计信息。
在这里，我们演示后者，通过在FindSpatiallyVariableFeatures()中设置method='moransi'。 Moran's I计算总体空间自相关性，并给出一个统计量（类似于相关系数），该统计量用于衡量要素对空间位置的独立性。这使我们能够根据基因表达的变化对基因进行排名。为了便于快速估算此统计信息，我们实施了一种基本的分级策略，该策略将在Slide-seqpuck 上绘制一个矩形网格，并对每个分级中的基因和位置取平均值。x和y方向上的bins 分别由x.cuts和y.cuts参数控制。此外，虽然不是必需的，但安装可选Rfast2软件包（install.packages('Rfast2')），将通过更有效的方式来减少运行时间。
DefaultAssay(slide.seq) <- "SCT"slide.seq <- FindSpatiallyVariableFeatures(slide.seq, assay="SCT", slot="scale.data", features=VariableFeatures(slide.seq)[1:1000], selection.method="moransi", x.cuts=100, y.cuts=100)现在，我们可视化由Moran's I识别的前6个基因的表达。
SpatialFeaturePlot(slide.seq, features=head(SpatiallyVariableFeatures(slide.seq, selection.method="moransi"), 6), ncol=3, alpha=c(0.1, 1), max.cutoff="q95")jing 生信人
空间转录组（ST）能够在单细胞或亚细胞水平上全面刻画组织结构，帮助理解传统技术无法获得的生物学机制。过去的几年已经开发出了各种各样空间转录组分析技术，这些技术被Nature Methods评为2020年的最佳方法。因此，今天小编和大家分享一篇今年10月发表在Genome Research(IF:9.043)杂志上的综述。文章总结了空间转录组数据分析方法的现状，并讨论了它们如何在不同的技术平台上运行，可以说一文带你了解空间转录组，感兴趣的小伙伴get起来吧。
Advances in spatial transcriptomic data analysis
空间转录组数据分析进展
一．文章背景
多细胞生物由组织和器官组成，组织和器官都专门从事一个生物过程。虽然正常细胞都拥有相同基因组，但它们的基因表达模式和形态可能截然不同。这种差异不仅是由内部基因调控引起，也是由外部组织环境发出的信号引起的。尽管几十年来基因组研究积累了大量细胞类型特异性基因调控的信息，但目前对细胞外部环境相互的理解仍然有限。不过随着近年来技术的爆炸式增长，已经实现了细胞异质性和组织及器官空间组织的系统级刻画。其中最值得注意的是单细胞RNA-seq技术(scRNA-seq)的快速发展，使得在一个器官中描绘和比较细胞的基因表达模式成为可能。scRNA-seq在发现新细胞类型方面发挥了关键作用，并为创建不同物种的综合细胞图谱奠定了基础。然而，由于实验过程通过机械和酶解步骤创造单细胞悬浮液，这不可避免地破坏了原始的组织结构。因此，仅从细胞成分重建组织结构是极其困难的。
二．空间转录组主要技术及数据
目前在单细胞分辨率和全基因组规模上研究组织的方法主要有三种，分别为序列荧光原位杂交(FISH)，及连接已知坐标的所有转录本的空间标记和选择基因用于原位测序(ISS)。文章主要使用了四个数据集（图1）来说明不同空间转录组分析的结果或方法:1）乳腺肿瘤生物标本切片的全基因空间转录组数据集。2）全小鼠冠状脑切片的亚细胞空间数据集。3）人类心脏全基因组空间数据集。4）老鼠体感皮层内数百个细胞的10000个基因的亚细胞空间数据。从原始ST数据集获取基因表达矩阵和相应的空间坐标不是一个简单的过程，它包含许多预处理步骤。这些步骤通常依赖于技术或平台，但有一些重复的步骤是固有的，例如基于成像数据的图像配准、拼接和细胞分割（图2）。对于基于图像的ST数据，如FISH和ISS技术，最常见的图像处理步骤是图像校正、拼接、配准、分割，然后定位和解码。另一方面，有些数据并不一定需要成像，而是通过在一个定义的空间单元内捕获转录本，并在测序步骤之前将它们与一个已知的坐标系统连接。因此，这些方法通常不依赖于上述图像处理步骤。然而，测序后还需要额外的步骤来将转录本映射回它们的空间坐标。当有组织图像时，可以将其与空间坐标系统叠加。不过无论技术如何差异，ST分析的一个共同目标是连接和整合来自基因表达和细胞转录本位置的信息，这对于提取有用的生物信息至关重要(图3)。
图1 文章使用的数据集
图2 原始空间转录组数据的预处理
图3 空间转录组学分析方法综述
三．ST数据识别细胞类型
细胞类型的识别和定位是ST数据分析最基本的任务，因此在这一部分作者对ST数据识别细胞类型进行了介绍。如果数据有单细胞分辨率，无监督聚类结合手动或自动注释是识别细胞类型的常见方法(图4A)。由于细胞类型识别不需要空间信息，因此与scRNA-seq分析相似，如基于community的方法，如Louvain和Leiden聚类。为了说明，作者使用MERFISH冠状切片数据集并应用Leiden聚类，总共得到19个不同的类。然后对这些簇进行注释并映射回空间坐标(图4B)。目前有两种估算细胞类型组成的方法(图4C)。第一种方法是评估每个位点上表达基因中细胞类型特异性标记的富集情况。这种方法速度很快，每次可以执行一种细胞类型。然而，结果是定性的，表明存在或不存在一种细胞类型。第二种方法是去卷积，目的是定量估计每个位置不同类型细胞的比例。已经开发了许多去卷积方法用于RNA-seq数据分析，原则上也可以应用于ST分析，不过ST数据具有某些不同的特性，如与每个位置关联的细胞数量很少。因此，通常使用专为ST分析设计的方法。在这些方法中，RCTD使用了基因计数的线性回归模型，进一步纳入了针对平台特异性变异的随机效应项，基因表达水平采用泊松分布模型。Cell2location使用类似的方法，但基因表达使用负二项分布模型，它还可以模拟特定平台和特定位置的效果。SpatialDWLS使用两个步骤来减少噪声：第一步使用富集分析来识别细胞类型，第二步使用衰减加权最小二乘法量化每种细胞类型的相对比例。SPOTlight使用种子非负矩阵分解(NMF)回归，并使用细胞类型marker基因和非负最小二乘(NNLS)进行初始化，用于后续的去卷积。DSTG使用基于图的卷积网络。DestVI使用变化推理方法去卷积。作者使用Visium心脏数据集和匹配的scRNA-seq数据进行细胞类型富集(图4D)和空间去卷积(图4E)，对这些方法具体说明。研究细胞类型定位的一种补充方法是使用scRNA-seq数据作为起点，然后根据空间表达轮廓的相似性重构空间信息。在ST技术爆炸之前，使用传统方法只能获得少数标志性基因的空间信息。利用这些有限的信息开发了tomo-seq和Geo-seq技术，能从2D切片中获得的基因表达谱重建3D模式。不过这些空间信息不是直接从数据中测量的，因此仍然是推测性的。随着ST技术在过去几年的快速发展，现在可以直接测量空间信息，并进一步与scRNA-seq数据整合。因此，新的方法以更平衡的方式整合scRNA-seq和ST数据。例如，使用了与平台无关的相互最近邻(MNN)方法来对齐这些数据类型，从而生成细胞位置映射。其中DEEPsc使用人工神经网络来预测空间位置；GLUER结合NMF、MNN算法和深度神经网络对数据进行对齐；Tangram将scRNA-seq和ST数据集对齐，同时优化了scRNA-seq数据和空间数据中每个基因之间的空间相关性。NovaSparc和D-CE也类似，不过对齐可以是概率的，也可以是确定性的。在已知细胞总数的前提下，Tangram对齐的确定性模式也可以作为一种去卷积方法。
图4 利用空间转录组数据进行细胞类型识别的策略
四．刻画转录组谱的空间模式
ST分析的关键不仅在于刻画细胞类型，还在于它们的空间组织。这对于研究组织结构和细胞间相互作用至关重要(图5A,C,E)。成对富集分析可以用于识别可能相邻的细胞类型对，为了说明如何研究空间网络模式和细胞邻近，作者使用MERFISH冠状切片数据，基于每个细胞的物理坐标创建了一个细胞邻近网络。细胞-细胞邻近网络连接热图显示不同细胞类型之间的连接性及所研究组织的空间拓扑结构(图5B)。图5D显示了对单个生态位的详细探索，将特定的细胞识别为源，然后描述它们与其他邻近细胞类型的连接。spatialDE使用一个随机效应模型，该模型包含两个terms，分别对应于空间成分和非空间成分。空间变量可以指定为各种形式，如线性、周期或高斯过程。空间变异性的程度由这两项所解释的方差比来量化。SOMDE就使用类似的方法，但通过使用自组织映射变换压缩空间信息来提高计算效率；Trendsceek将空间格局模型作为一个标记点过程；SPARK通过广义线性空间模型对空间计数数据进行建模；有些方法主要考虑局部连续性，如binSpect检测空间相干基因，即倾向于在相邻细胞中共同表达的基因。另一种方法是量化空间结构的扩散步骤。作者以MERFISH冠状脑切片数据为例，利用binSpect识别具有空间相干模式的基因，排名靠前的基因如图5F所示。
图5 空间模式分析
五．亚细胞结构分析
在这一部分作者使用小鼠体感皮层数据集来说明亚细胞数据分析的一些关键概念(图6)。在亚细胞分辨率的数据集中，每个点通常代表一个单一的转录本(图6A)。分析亚细胞基因表达模式可以作为空间分析的一种替代方法，也可以用于提高细胞分割的准确性(图6B)。通过分析共定位模式(图6C)和每个细胞内的转录动力学(图6D)，可以发现基因之间或基因与亚细胞结构之间的空间关系。许多方法尝试利用亚细胞基因表达模式来规避细胞分割。如SSAM直接将细胞类型标签分配给像素，而不进行细胞分割。stLearn使用类似的方法，但将空间近端像素进一步集群。Spage2vec也使用了类似的方法，但采用了一种神经网络公式。另外，也开发了基于已知细胞类型特定签名的监督细胞类型映射策略。如用朴素贝叶斯模型为HDST数据分配细胞类型。亚细胞基因表达模式可以反过来用于改善细胞分割。例如，Baysor模型利用马尔可夫随机场模型对亚细胞基因表达模式进行建模，并进一步整合细胞形状标记信息来提高细胞分割的准确性。Sparcle使用Dirichlet过程混合模型以及相邻细胞与相邻转录本之间的转录本距离来增强细胞分割。JTSA使用EM算法迭代改进像素及基因表达谱分类和细胞边界标注。分析基因表达的亚细胞模式也可以提供新的生物学见解。如已经开发了一种原位RNA速度方法，利用亚细胞RNA定位信息来推断转录速率。由于新转录的RNA在细胞核中积累，而成熟的mRNA需要被运输到细胞质中进行翻译(图6D)，每个基因相关的核转录本与细胞质转录本的相对组成可以用来估计转录活性。此外，利用过氧化物酶APEX2对RNA进行直接邻近标记(APEX-seq)，可以高分辨率地识别细胞质中共定位的mRNA种类。对结果数据的分析发现，共定位的RNA与已知的蛋白质共定位模式之间存在显著的对应关系，这表明RNA共定位可能有助于局部蛋白质翻译和复杂的形成(图6c)。
图6 亚细胞分辨率空间转录组学分析示意图
六．细胞与组织环境沟通
ST分析的一个重要目标是研究细胞如何与组织环境沟通(图7)。细胞行为可以通过直接的物理相互作用、分泌分子或与细胞外基质相互作用受到组织环境的影响(图7A)。细胞间的通信通常是空间协调的，并且具有高度的细胞类型特异性。因此，即使在相同的细胞类型中，细胞类型组成的变化也可能导致基因表达的显著变化(图7B,C)。R包Giotto采用双向比较，比较同一细胞类型被不同相邻细胞包围的基因表达模式，来识别相互作用改变基因。与单独使用基因表达信息相比，使用空间信息可以显著减少假阳性配体受体活性预测的数量，CellPhoneDB中使用了类似的方法。在这一部分作者使用Cell2location来推断不同细胞类型的位置，然后比较不同细胞邻区相关的基因表达模式。也有其他方法用于量化相邻细胞类型的影响，包括卷积神经网络。及将基因表达谱分解为空间分量和非空间分量，然后利用邻域内的细胞类型组成来估计空间分量的方法。此外，也有算法根据细胞与细胞的相互作用模式重建空间位置。
图7 从配体-受体相互作用推断细胞通信
七．空间数据分析和可视化工具
开发空间数据结构工具如今越来越重要，理想情况下导入原始数据(图8A)及图像分析(图8B)，和分析结果及数据准备发表(图8C)都应该是在个人电脑上。这一部分作者介绍了下游数据分析工具。这些工具基本都是用R或Python编写的。如R包Giotto可以用于各种各样的ST技术，核心由一个专门为空间数据设计的对象组成。同时，Giotto提供了一个基于浏览器的可视化工具，该工具允许用户导出他们获得的结果，并交互式的探索空间数据集。Seurat作为一个流行scRNA-seq分析R包为人所熟知，它也开始提供一些针对空间数据可视化和识别空间表达模式的功能。此外其他工具如STUtility和SPATA也建立在Seurat数据结构之上创建了更全面的管道，这些管道目前仅适用于ST技术。STUtility是专为ST技术开发的，提供了多种成像和数据分析方法。SPATA也关注ST数据，其开发目的是促进Seurat和Monocle集成。除了可视化和常见的数据分析功能，SPATA还可以识别或描绘空间轨迹。Squidpy是空间版的SCANPY, SCANPY是用于scRNA-seq分析的Python库，它还提供图像层面的分析。Stlearn是另一个用于ST数据分析的Python库，它专注于将基因表达和图像信息集成。这些包或工具大多数都是独立的实验室开发的，导致这些数据结构不一定共享相同的数据格式。因此R/Bioconductor团队精心设计普遍适用的数据结构，并于最近发布了spatialExperiment class的第一个版本。这是一种新的S4类，它扩展了singleCellExperexperiment class，能用于操作几种类型的ST数据集。一些R包已经开始使用这种数据结构，如SpatialLIBD和Spaniel，它们都擅长创建交互式程序来可视化ST数据集。总之，这些努力可以在未来促进不同工具之间的交互操作性。
图8 交互式探索性分析管道概述
到这里这篇文章的主要内容就介绍完了，文章对空间转录组数据分析方法进行了总结，并使用空间转录数据集进行了举例说明。近几年来空间转录组领域快速发展，是一个热点方向，感兴趣的小伙伴不要错过呀。
参考文献
1.Advances in spatial transcriptomic data analysis；
2.High-throughput spatial mapping of single-cell RNA-seq data to tissue of origin；
上一篇：5个半月的宝宝学什么好
下一篇：上联-母爱深如海，如何对下联？
发表评论