2. BayesSpace算法原理与数学模型归纳总结

上一章我们聊了空间转录组学的基本背景，以及为什么需要专门的算法来分析那些带坐标的基因表达数据。这一章我们正式进入 BayesSpace 的核心。我个人习惯先讲清楚一个算法“凭什么能工作”，再去看它的代码怎么写。所以，我们先从它的数学引擎说起。

你可能会问：“空间转录组数据不就是一堆点吗？为什么非要搞贝叶斯？” 打个比方：普通聚类就像在地图上画圈，每个点只能属于一个圈；但空间数据里相邻的spot（捕获点）本来就该互相影响——比如一个spot的基因表达大概率跟它周围的几个spot相似。BayesSpace把这个直觉用数学语言写死了，而且加上了统计推断的“概率”维度。

2.1 模型框架：隐马尔可夫随机场与空间先验

BayesSpace本质上是一个 隐马尔可夫随机场（HMRF）模型。这里的“隐”指每个spot真正的“组织类型”（也就是潜变量）我们看不到，只能通过观测到的基因表达来推测。而“随机场”意味着相邻元素在空间上存在依赖关系。

形式上，我们把每一个spot i 的基因表达向量记为 \(y_i\)（高维，比如200个高变基因），对应的潜变量 \(z_i\) 表示该spot所属的空间域（如1到K类）。模型的核心假设是：

• 观测层：给定 \(z_i = k\)，\(y_i\) 服从一个多元正态分布，均值和协方差矩阵由第k个域的参数决定。
• 空间先验层：潜变量 \(z_i\) 的联合概率由一个 Potts模型（一种马尔可夫随机场）控制——两个相邻的spot越倾向于被分配到相同域，能量越低。

用数学写出来就是：

p(z | β) ∝ exp{ β * Σ_{(i,j)∈E} I(z_i = z_j) }

其中：
- β 是空间平滑参数（越大越鼓励相邻spot同域）
- E 是邻接边集合（比如基于visium的六边形网格）
- I(·) 是指示函数

注意这个β在经典Potts模型里可能是固定值，BayesSpace把它也当作一个随机参数来推断。那我之前做过一个项目——用的还是旧版BayesSpace——没有去推断β，结果聚类结果对β取值特别敏感，换了不同组织切片就得手工调。后来换成带β推断的新版本，省了很多麻烦。

2.2 先验设定与变分推断

纯贝叶斯公式长这样：后验 ∝ 似然 × 先验。这里我们需要给每个高斯分量（每个域）的均值μ_k、协方差Σ_k，以及空间参数β都配上先验。

• μ_k：用Normal先验，均值全局均值，协方差大一点表示不确定。
• Σ_k：Inverse Wishart先验，提供对协方差结构的基本控制。
• β：Gamma先验（因为β>0），形状和速率参数可自由设置。

后验分布没有解析解。直接跑MCMC（马尔可夫链蒙特卡洛）在几万个spot的数据上会慢到令人抓狂。所以BayesSpace用了 平均场变分贝叶斯——假设潜变量之间在变分后验下相互独立，然后用坐标上升法迭代更新。我当年第一次跑这个算法时，看见它几十秒就收敛了，对比MCMC动辄几小时，感慨变分真的是空间转录组分析的“工业级”解决方案。

变分迭代的核心步骤简化如下（伪代码只展示思想，实际paper中有详细公式推导）：

初始化：随机分配每个spot到某域
循环直到收敛：
  1. 固定变分参数，更新每个域的μ, Σ（用加权最大似然）
  2. 固定μ, Σ和β，更新每个spot的分属概率q(z_i)（利用空间邻接信息）
  3. 更新β的后验（用Gamma分布拟合）
检查ELBO（证据下界）是否增长

这里有一个关键操作：空间平滑项进入每个spot的后验概率更新公式。简单说，如果某个spot周围邻居大部分属于域k，那么它被分到域k的概率会叠加一个“奖励因子”。这个奖励正比于β。

2.3 核心参数：空间平滑参数β与邻接矩阵

我把这两个放在一起讲，因为它们是BayesSpace与普通GMM（高斯混合模型）最大的区别。

2.4 算法流程：从输入到输出

整体流程我用一个SVG图来归纳，你看一眼就能记住：

从左到右主路径很清晰：你先准备好标准化后的高变基因矩阵（通常2000个基因以内），然后根据坐标建邻接图，接着用变分贝叶斯迭代求解，最后得到每个spot的域归属。右上角的虚线框表示BayesSpace还提供一个增强功能——spatialEnhance()可以把spot级别的聚类结果“插值”到亚spot网格上，生成类似增强分辨率的热图。不过这个功能不是必须的，我一般只在做高级可视化时才调用。

2.5 变分更新的数学细节（选读）

你可能不想看满屏的推导，但理解关键公式有助于调参。以下是最核心的 q(z_i) 更新方程：

对于 spot i，其属于域 k 的变分概率 q(z_i=k) 正比于：

exp {  log_pi_k  +  log_N(y_i; μ_k, Σ_k)  +  β * Σ_{j∈N(i)} q(z_j=k)  }

其中：
- log_pi_k：域k的混合权重对数
- log_N(·)：多元正态的对数密度
- 最后一项是空间“奖励”：β乘以邻居中属于域k的期望数量

你看，只有最后一项跟邻居有关。当β=0时，BayesSpace退化成一个标准的GMM，完全忽略空间信息。所以β的大小直接决定了空间聚类的“块状程度”。我建议新手先试β=1，然后在结果图上观察——如果边界区域出现了许多孤立的“小斑点”（斑点噪声），说明β偏小；如果边界被过度模糊、不同区域连成一片，说明β偏大。

2.6 与SpaGCN的区别——先打个底

这一节我们不深入SpaGCN，但为了教程体系连贯，先记下BayesSpace的核心优势：它是一个纯概率模型，输出的是每个spot属于每个域的概率分布，而不是硬分类。这意味着你可以用这些概率做后续的空间差异表达分析（比如检测两个域之间哪些基因差异表达），而且不确定性被量化了。SpaGCN则更偏图神经网络+聚类，本质是一个判别式模型。我们到第10章再专门对比。