空间聚类算法BayesSpace

空间基因表达技术能够在保留空间背景信息的同时，全面测量转录组图谱。然而，现有的分析方法并没有解决技术分辨率有限或有效利用空间信息的问题。
BayesSpace是一种完全贝叶斯统计方法，用于聚类分析来增强空间转录组数据的分辨率，并且可以通过SingleCellExperiment对象无缝集成到当前转录组的分析工作流程中。BayesSpace利用多元t分布对基因表达的低维表示进行建模，然后通过Potts模型合并空间信息，这使得相邻点属于同一簇。这种方法借鉴了用于图像分析的成熟的空间统计方法。有效地利用空间信息可以极大地提高空间数据聚类的性能。

安装

BiocManager::install("BayesSpace")
library(SingleCellExperiment)
library(ggplot2)
library(BayesSpace)

读入数据

readVisium()读入的是Spatial Ranger的输出数据，格式参考：https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/output/overview

sce <- readVisium("path/to/spaceranger/outs/")
#得到的sce是一个SingleCellExperiment对象

数据预处理

BayesSpace需要很少的数据预处理，在这里这一步封装在了spatialPreprocess()这个函数。
spatialPreprocess()对count矩阵进行了log标准化( log.normalize=T)，并对top n.HVGs高变基因进行了PCA处理，保留了top n.PCs个主成分。此外，spatial sequencing platform被作为metadata的信息添加到SingleCellExperiment对象用于下游分析。如果不想进行PCA这一步，在进行spatialPreprocess()的时候可以设置skip.PCA=TRUE。

set.seed(102)
sce <- spatialPreprocess(sce, platform="Visium", n.PCs=15, n.HVGs=2000, log.normalize=T)
saveRDS(sce, "sce.bc1a.rds")

没有enhanced

聚类

选择cluster数目

我们使用qTune()和qPlot()功能来帮助选择q值，也就是用于后续分析的cluster数。（q值直接决定了后面聚出多少个cluster）

• qTune()：对q的多个指定值（默认：3-7）进行贝叶斯聚类算法并计算他们的平均伪对数似然值（average pseudo-log-likelihood），也就是对qs的每一个设定值做一个fitting model。它接受spatialCluster()的任何参数。
• qPlot()：绘制q值的伪对数似然值，帮助选择合适的q值。

npcs = 15
sce <- qTune(sce, qs=seq(2, 15), platform="Visium", d=npcs)
p <- qPlot(sce)
p
ggsave("Elbowplot_bc1.pdf", p, width = 8, height = 6)

bayes聚类

spatialCluster()功能对spots进行聚类，并向SingleCellExperiment对象增加预测的cluster标签。通常，我们建议至少10,000次迭代（nrep=10000），但是这个时间通常会很长，所以为了演示，这里只做1000次迭代。（注意：为了结果的一致性，这一步需要设置随机种子。）

### 贝叶斯聚类
sce <- readRDS("sce.bc1a.rds")
# 设置聚类参数
npcs = 30; q = 15 ##会聚出15个cluster
# bayes聚类(nrep: 马尔可夫链蒙特卡罗算法MCMC迭代次数，burn.in:需要排除的MCMC磨合期迭代次数)
set.seed(149)
sce <- spatialCluster(sce, q=q, platform="Visium", d=npcs, nrep=1000, burn.in=100, 
                        init.method="mclust", model="t", gamma=3)

在这一步之后，mclust初始化(cluster.init)和贝叶斯cluster分配(spatial.cluster)的结果都已经得到并储存在SingleCellExperiment的colData中。

head(colData(sce))

空间聚类可视化

clusterPlot()的结果是基于ggplot2的，因此可以使用ggplot2的语法对结果图进行优化。

### 结果展示
p1 <- clusterPlot(sce, label = "cluster.init") + ggtitle("mclust") + theme(plot.title = element_text(hjust=0.5))
p2 <- clusterPlot(sce, label = "spatial.cluster") + ggtitle("bayes") + theme(plot.title = element_text(hjust=0.5))
p3=clusterPlot(sce, color="black") +
    theme_bw() +
    xlab("Column") +
    ylab("Row") +
    labs(fill="BayesSpace\ncluster", title="Spatial clustering of BC1")
p = p1|p2|p3
p
ggsave(paste0("cluster_mclust_bayes_bc1_q", q, ".pdf"), p, width = 15, height = 5)

# 保存结果
write.csv(bayes_cluster, "bayes_cluster_bc1.csv")
saveRDS(sce, file = paste0("sce.bc1b.q", q,".rds"))

增强分辨率

Clustering at enhanced resolution

spatialEnhance()功能将增强主成分的分辨率，并将这些PC以及subspot分辨率的预测聚类标签添加到新的SingleCellExperiment对象。和前面的spatialCluster一样，为了演示这里的迭代也是使用了1000，但是建议使用10000或更大的迭代值。

sce.enhanced <- spatialEnhance(sce, platform="Visium", d=npcs, q=q, model="t", 
                                   gamma=3, nrep=1000, burn.in=10, save.chain=FALSE)

增强的SingleCellExperiment包括了原始sce中父类spot的指标(spot.idx)，也包含了subspot的指标(subspot.idx)。它将offsets添加在原始点坐标上，并提供了增强的聚类坐标(spatial.cluster)。

head(colData(sce.enhanced))

可视化

p <- clusterPlot(sce.enhanced)
p
ggsave("enhanced_resolution_bc1.pdf", p, width = 8.5, height = 6.5)
saveRDS(sce.enhanced, file = "sce.enhanced.rds")

Enhancing the resolution of gene expression

BayesSpace对基因表达矩阵的主成分进行运算，因此spatialEnhanced()计算的是增强分辨率的PC向量。增强的基因表达并不是直接计算的，而是通过回归算法计算的。对于每个基因，使用每个点的PC向量训练了模型来预测spot-level的基因表达，并使用拟合的模型来预测来自subspot PC的subspot表达。
基因表达增强是使用enhancedFeature()函数实现的。BayesSpace默认使用xgboost预测基因表达，但也可以通过模型参数使用线性和 Dirichlet狄利克雷回归。使用xgboost时，我们建议将nrounds参数设置为0来自动调整它，尽管这大大增加了运行时间（比预先指定的nrounds慢了约4倍）。
EnnhancedFeatures()可用于计算所有基因或感兴趣基因集的subspot表达水平。在这里演示了四个常见基因。

markers <- c("LYZ","S100A9","ACTA2",'COL1A1')
sce.enhanced <- enhanceFeatures(sce.enhanced, sce.ref = sce, nrounds = 0,
                                  feature_names = markers)

默认情况下，计算的是log标准化的expression(logcounts(sce))，尽管也可以指定其他assays或者arbitrary features矩阵。

logcounts(sce.enhanced)[markers, 1:5]

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

rowData(sce.enhanced)[markers, ]

增强基因表达的可视化

p1 <- lapply(markers, function(x){
    featurePlot(sce.enhanced, x) + ggtitle(paste0("imputed_",x)) + theme(plot.title = element_text(hjust=0.5))
  })
  p2 <- lapply(markers, function(x){
    featurePlot(sce, x) + ggtitle(x) + theme(plot.title = element_text(hjust=0.5))
  })
  plist <- c(p1,p2)
  p <- wrap_plots(plist, ncol = 4)
  ggsave("markers_expression_enhanced.pdf", p, width = 15, height = 6.5)  
  
  saveRDS(sce.enhanced, file = "sce_enhanced.rds")

这个增强效果好棒哦

参考：https://edward130603.github.io/BayesSpace/articles/BayesSpace.html

---

主要函数	含义
spatialPreprocess	SingleCellExperiment对象的预处理
qTune/qPlot	帮助选择q值
spatialCluster	空间聚类
clusterPlot	绘图函数，展示聚类结果
spatialEnhance	增强聚类的分辨率
enhanceFeatures	增强基因表达的分辨率
featurePlot	绘图函数，展示基因表达