10X空间转录组聚类分析之BayesSpace算法聚类

hello，大家好，要过礼拜了，但是今天还是要好好工作，今天我们来分享一个针对空间转录组的聚类算法，目前大多数多空间转录组的聚类还是基于找高变基因降维聚类，但是文献中采用这个方法的很少，我们来看看新方法有什么改进。文献在Spatial transcriptomics at subspot resolution with BayesSpace,发表于Nature Biotechnology（IF54，相当高了）。

单细胞转录组测序（scRNA-seq）实现了高通量和高分辨率的基因表达谱分析，但由于样品制备过程中组织被分离，空间信息并没有被保留。近年来，新兴起的空间转录组学技术能够在保留空间背景的同时全面检测转录组谱，使人们对转录组表达的空间位置有了更深的认知，并且为生物学功能和病理提供重要的见解。

目前，利用空间转录组学，人们已经对肿瘤异质性、脑功能等不同领域产生了更加深刻的认识。但这些空间基因表达平台仍然存在分辨率的技术限制。虽然后续开发的结合荧光原位杂交（FISH）技术的空间转录组测序方法以及其他方法，如Slide-seq和ZipSeq等提供了更高的分辨率，但大多数方法仍然存在通量较低、敏感性较低、依赖某种特定条件或不能够被广泛使用等限制问题。

此外，在测序技术发展的同时，还需要新的统计方法来分析空间基因表达数据，从而有效地利用现有的空间信息进行科学研究。其中，聚类分群是此类数据分析的重要一步，也是下游深入分析的基础，后续研究诸如细胞类型或组织注释、差异表达等均依赖有效的分群。有学者提出了一种隐马尔科夫随机场模型（HMRF），通过对基因表达和空间邻域结构联合建模，将低分辨率原位杂交数据聚类到不同的空间域，选择空间差异表达基因进行聚类。但是，现有的大多数空间基因表达数据分析方法往往依赖于非空间scRNA-seq数据的聚类方法，难以精确反映空间转录组所提供的信息。

近日，美国西雅图弗雷德·哈钦森癌症研究中心、华盛顿大学等单位的联合研究团队在Nature Biotechnology发表了题为“Spatial transcriptomics at subspot resolution with BayesSpace”的文章。该研究提出了BayesSpace算法，这是一种利用空间转录组数据中的邻域结构来增加子点级别分辨率的计算方法，通过使用贝叶斯统计来实现超分辨率图像分析。研究人员利用BayesSpace进行空间聚类改进了对空间分布的组织域的识别，提高了基因表达图谱的分辨率，并可以重现接近单细胞分辨率的真实空间结构。

主要研究内容

BayesSpace基本原理

文章介绍，BayesSpace通过对低维的基因表达矩阵进行建模，并通过空间先验知识诱导真实的邻近点聚集，以此进行推广，从而实现空间聚类。这种方法来自于此前开发的用于图像分析和微阵列数据的空间统计方法。与已有的方法相比，BayesSpace允许对聚类结构和错误项进行更灵活的调整和规范。

图1. BayesSpace基本原理，

BayesSpace空间聚类提高了对脑组织中已知皮层的认识

为了检测BayesSpace的性能，研究人员使用Maynard等人公开发表的12个背外侧前额叶皮层 （DLPFC） 样本的Visium空间表达谱数据，以及每个样本的6个皮质层和白质的手工注释，这些是作为R包spatialLIBD的一部分。借助该数据集，研究人员评估了BayesSpace识别不同皮层特定表达轮廓的能力，并将其性能与其他空间和非空间聚类方法进行比较。

结果显示，BayesSpace大大优于原始的spatialLIBD聚类分区，以及为空间转录组数据开发的其他所有非空间聚类算法和空间聚类方法。除BayesSpace外，大多数聚类分区表现出大量的噪声，且聚类之间缺乏明确的空间分隔。相比之下，BayesSpace利用了空间信息平滑数据，并提供不同的集群分层。此外，BayesSpace的运行时和内存占用与其他空间聚类方法相当。

图2. BayesSpace性能评估及比较

BayesSpace能够鉴定出易被其他方法遗漏的组织结构

随后，研究人员还使用BayesSpace分析了由Thrane等人首次注释和描述的黑色素瘤空间转录组样本，这些数据包括了人工标注识别的黑色素瘤、间质和淋巴组织的区域。结果显示，利用BayesSpace得到的4个空间聚类与手工标注的组织类型相吻合。此外，BayesSpace增强的空间聚类提供了更高分辨率的组织类型图，例如增强识别的肿瘤边缘的淋巴样区域和可能的免疫浸润肿瘤区域。这些区域在最初的分辨率下是无法识别的，并且这些区域在很大程度上也没有被其他聚类方法识别。

差异表达分析结果提示，淋巴区域有一个独特的表达谱，其淋巴细胞标志物CD52和MS4A1的表达升高，黑色素瘤标志物MCAM和SPP1的表达相对于周围肿瘤边界的表达降低，四个聚类间的增强分辨率差异表达分析强调了基因表达的额外空间变异。

图3. BayesSpace鉴定出黑色素瘤样本中的肿瘤近端淋巴组织结构

BayesSpace能够区分浸润性导管癌的瘤内异质性

接下来，研究人员进一步分析了浸润性导管癌（乳腺癌的一种病理分型）的组织切片以确定聚类分群的生物学相关性。每一张组织切片，病理学家注释了主要的浸润性癌、原位癌和良性增生的区域，从中可以得到每个点的真实标签。结果显示，BayesSpace聚类与组织病理学注释基本一致。

此外，如果没有苏木精和伊红（H&E）染色或免疫荧光染色作为肿瘤标志物，肿瘤-基质界面在组织学上不能被完全描述。BayesSpace增强的聚类可以识别出组织内的异质性，并得到关键肿瘤标记基因的明确数据支持。也就是说，已知肿瘤标记基因的空间表达模式和这些聚类之间的差异表达分析与临床和组织病理学注释基本一致。例如，在整个肿瘤群体中观察到ERBB2和ESR1基因的高表达水平；非肿瘤细胞群1、7和10的特征是免疫基因的表达，如PTPRC（白细胞共同抗原CD45） 的高表达。这些空间表达模式表明，侵袭性肿瘤存在明显瘤内转录异质性，目前的组织病理学分析方法无法识别这些差异，表明空间转录组数据相对于单纯免疫荧光的优势。

图4. BayesSpace区分浸润性导管癌的瘤内异质性

研究总结

综上所述，该研究首次报道了BayesSpace这一基于空间转录组模型的聚类方法。BayesSpace使用t分布错误模型来识别空间聚类，这些聚类对于技术噪声引起的离群值的存在更加鲁棒。随后的应用及分析结果证明，BayesSpace在识别具有相似表达谱的空间群体和提高空间转录组的分辨率方面有较高的实用性。

我们来看看代码

准备数据

BayesSpace supports three ways of loading a SingleCellExperiment for analysis.（这个是R语言常见的数据结构）

Visium datasets processed with Space Ranger can be loaded directly via the readVisium() function. This function takes only the path to the Space Ranger output directory (containing the spatial/ and filtered_feature_bc_matrix/ subdirectories) and returns a SingleCellExperiment.（直接读取10XSpaceranger的文件）.

sce <- readVisium("path/to/spaceranger/outs/")

其次，为 BayesSpace manuscript分析的所有数据集都可以通过 getRDS() 函数轻松访问。 这个函数有两个参数——数据集的名称和数据集中样本的名称。

melanoma <- getRDS(dataset="2018_thrane_melanoma", sample="ST_mel1_rep2")

也可以从计数矩阵和行列数据表手动构建 SingleCellExperiment 对象。 BayesSpace 仅要求将点数组坐标作为 colData 中名为 row 和 col 的列提供。 （注意，Visium 数据集的增强还需要组织图像中每个点的像素坐标，但在这种情况下，数据集应使用 readVisium() 加载，它会自动加载这些数据。）

library(Matrix)

rowData <- read.csv("path/to/rowData.csv", stringsAsFactors=FALSE)
colData <- read.csv("path/to/colData.csv", stringsAsFactors=FALSE, row.names=1)
counts <- read.csv("path/to/counts.csv.gz",
                   row.names=1, check.names=F, stringsAsFactors=FALSE))

sce <- SingleCellExperiment(assays=list(counts=as(counts, "dgCMatrix")),
                            rowData=rowData,
                            colData=colData)

前处理

BayesSpace 需要最少的数据预处理，提供了一个辅助函数来自动化它。

spatialPreprocess() 对计数矩阵进行对数归一化，并对顶部 n.HVGs 高度可变的基因执行 PCA，保留顶部 n.PCs 主成分。 此外，空间排序平台作为meta数据添加到SingleCellExperiment 中，用于下游分析。 如果不想重新运行 PCA，运行带有skip.PCA=TRUE 标志的 spatialPreprocess() 只会添加 BayesSpace 需要的meta数据。

在这里，省略了对数归一化，因为所有通过 getRDS() 可用的数据集都已经包含对数归一化计数。

set.seed(102)
melanoma <- spatialPreprocess(melanoma, platform="ST", 
                              n.PCs=7, n.HVGs=2000, log.normalize=FALSE)

聚类

Selecting the number of clusters

我们可以使用 qTune() 和 qPlot() 函数来帮助选择 q，即我们分析中要使用的聚类数。

• qTune() 为 q 的多个指定值（默认为 3 到 7）运行 BayesSpace 聚类算法，并计算它们的平均伪对数似然。 它接受spatialCluster() 的任何参数。
• qPlot() 将伪对数似然绘制为 q 的函数； 我们建议在该图的肘部周围选择一个 q。

melanoma <- qTune(melanoma, qs=seq(2, 10), platform="ST", d=7)
qPlot(melanoma)

Clustering with BayesSpace

spatialCluster() 函数对spot进行聚类，并将预测的聚类标签添加到SingleCellExperiment。 通常，建议至少运行 10,000 次迭代 (nrep=10000)，但为了运行示例中使用了 1,000 次迭代。 （注意必须设置随机种子才能使结果可重现。）

set.seed(149)
melanoma <- spatialCluster(melanoma, q=4, platform="ST", d=7,
                           init.method="mclust", model="t", gamma=2,
                           nrep=1000, burn.in=100,
                           save.chain=TRUE)

mclust 初始化 (cluster.init) 和 BayesSpace 集群分配 (spatial.cluster) 现在都可以在 SingleCellExperiment 的 colData 中使用。

head(colData(melanoma))
#> DataFrame with 6 rows and 5 columns
#>            row       col sizeFactor cluster.init spatial.cluster
#>      <integer> <integer>  <numeric>    <numeric>       <numeric>
#> 7x15         7        15   0.795588            1               1
#> 7x16         7        16   0.307304            1               1
#> 7x17         7        17   0.331247            2               2
#> 7x18         7        18   0.420747            3               2
#> 8x13         8        13   0.255453            1               1
#> 8x14         8        14   1.473439            1               1

可视化

clusterPlot(melanoma)

由于 clusterPlot() 返回一个 ggplot 对象，因此可以通过组合熟悉的 ggplot2 函数对其进行自定义。 此外，参数调色板设置用于每个簇的颜色，clusterPlot() 将附加参数传递给 geom_polygon()，例如大小或颜色，以控制斑点边界。

clusterPlot(melanoma, palette=c("purple", "red", "blue", "yellow"), color="black") +
  theme_bw() +
  xlab("Column") +
  ylab("Row") +
  labs(fill="BayesSpace\ncluster", title="Spatial clustering of ST_mel1_rep2")

Enhanced resolution（加强精度）

Clustering at enhanced resolution

spatialEnhance() 函数将提高主成分的分辨率，并将这些 PC 以及子点分辨率的预测聚类标签添加到新的 SingleCellExperiment。 与我们上面对 spatialCluster() 的演示一样，我们在本示例中使用的迭代次数 (nrep=1000) 比我们在实践中推荐的次数 (nrep=100000 或更大) 少。

melanoma.enhanced <- spatialEnhance(melanoma, q=4, platform="ST", d=7,
                                    model="t", gamma=2,
                                    jitter_prior=0.3, jitter_scale=3.5,
                                    nrep=1000, burn.in=100,
                                    save.chain=TRUE)

The enhanced SingleCellExperiment includes an index to the parent spot in the original sce (spot.idx), along with an index to the subspot. It adds the offsets to the original spot coordinates, and provides the enhanced cluster label (spatial.cluster).

clusterPlot(melanoma.enhanced)

Enhancing the resolution of gene expression

ayesSpace 对基因表达矩阵的主要成分进行运算，因此 spatialEnhance() 计算增强分辨率的 PC 向量。增强的基因表达不是直接计算的，而是使用回归算法估算的。对于每个基因，训练使用每个点的 PC 向量的模型来预测点级别的基因表达，并使用拟合模型来预测来自subspot PC 的subspot表达。

基因表达增强在enhanceFeatures() 函数中实现。 BayesSpace 默认使用 xgboost 预测表达式，但也可以通过模型参数使用线性和狄利克雷回归。使用 xgboost 时，建议通过将 nrounds 参数设置为 0 来自动调整它，尽管这是以增加运行时间为代价的（实际上比预先指定的 nrounds 慢约 4 倍）。

EnhanceFeatures() 可用于估算所有基因或感兴趣基因子集的subspot水平表达。在这里，我们将通过增强四种标记基因的表达来证明：PMEL（黑色素瘤）、CD2（T 细胞）、CD19（B 细胞）和 COL1A1（成纤维细胞）。

markers <- c("PMEL", "CD2", "CD19", "COL1A1")
melanoma.enhanced <- enhanceFeatures(melanoma.enhanced, melanoma,
                                     feature_names=markers,
                                     nrounds=0)

By default, log-normalized expression (logcounts(sce)) is imputed, although other assays or arbitrary feature matrices can be specified.

logcounts(melanoma.enhanced)[markers, 1:5]
#>        subspot_1.1 subspot_2.1 subspot_3.1 subspot_4.1 subspot_5.1
#> PMEL     2.6428437   1.8550344   2.4704804   3.1827958   2.3572879
#> CD2      0.3489273   0.6066852   0.2315192   0.2210583   0.3489273
#> CD19     0.6170074   0.6528370   0.4164957   0.2558892   0.6259792
#> COL1A1   0.0000000   2.9053805   1.1940085   0.1023711   0.8157212

Diagnostic measures from each predictive model, such as rmse when using xgboost, are added to the rowData of the enhanced dataset.

rowData(melanoma.enhanced)[markers, ]
#> DataFrame with 4 rows and 4 columns
#>                gene_id   gene_name    is.HVG enhanceFeatures.rmse
#>            <character> <character> <logical>            <numeric>
#> PMEL   ENSG00000185664        PMEL      TRUE             0.804628
#> CD2    ENSG00000116824         CD2      TRUE             0.614575
#> CD19   ENSG00000177455        CD19      TRUE             0.697328
#> COL1A1 ENSG00000108821      COL1A1      TRUE             0.704845

可视化enhanced gene expression

featurePlot(melanoma.enhanced, "PMEL")

enhanced.plots <- purrr::map(markers, function(x) featurePlot(melanoma.enhanced, x))
patchwork::wrap_plots(enhanced.plots, ncol=2)

spot.plots <- purrr::map(markers, function(x) featurePlot(melanoma, x))
patchwork::wrap_plots(c(enhanced.plots, spot.plots), ncol=4)

方法真的不错，最后，大家周末愉快

生活很好，有你更好