第4章：SPARK算法原理（中）：零膨胀负二项模型处理稀疏计数

上一章我们聊了空间变异检测的基本框架，今天要深入一个非常实际的问题——如何处理空间转录组数据中的零值。为什么单独拿一章讲这个？因为零值太多了，多到如果不妥善处理，整个分析结果都会跑偏。

我刚开始接触空间转录组数据时，第一反应是“这数据也太稀疏了吧”。一个典型的10x Visium切片，高达70%-90%的spot-基因对都是零表达。如果你直接拿这些零去拟合标准模型，模型会非常困惑——它分不清这个零是因为基因在这个位置本来就不表达，还是表达了但没被检测到。

为什么标准模型搞不定稀疏数据

先看看传统方法怎么做。大多数RNA-seq差异性分析工具假设表达量服从负二项分布（Negative Binomial，NB）。这个模型有两个参数：均值 μ 和离散参数 φ。在处理bulk RNA-seq时它表现不错，因为每个样本的测序深度高，零值比例可控。

但空间转录组不一样。每个spot的UMI count本来就低，再加上测序灵敏度限制，很多低表达基因会出现大量“假阴性”零。标准负二项模型把所有的零都当作“真零”来处理，结果就是它严重低估了数据的离散程度，导致空间变异检测失去统计效力。

你想想看，如果一个基因在组织切片上其实存在空间模式，但某些区域因为捕获效率低出现了零值，标准模型会告诉你“这个基因不表达，没有模式”。这就是典型的假阴性。

零膨胀负二项模型的核心思想

SPARK采用的方案很聪明——它引入了一个混合模型，把零值拆成两个来源：

• 结构零（structural zero）：基因在这个spot真的不表达，生物学意义上的零
• 采样零（sampling zero）：基因本来有表达，但没被检测到，技术缺陷造成的零

ZINB模型用一个两部分的机制来处理：第一部分用一个逻辑回归模型判断零的来源；第二部分用负二项模型处理非零和“真零”。这就是所谓的“零膨胀”。

我习惯把这个过程理解成一个两步决策：首先是“这个spot到底表达不表达？”如果是，才进入“表达多少？”的问题。空间转录组的数据生成过程，本质上就是这么个层次结构。

模型组成：数学上的拆解

我们一步步看ZINB的数学形式。对于基因 g 在 spot i 的表达量 Y_gi：

Y_gi ~ ZINB(π_i, μ_i, φ)

其中：

P(Y_gi = 0) = π_i + (1-π_i) * NB(0 | μ_i, φ)
P(Y_gi = y) = (1-π_i) * NB(y | μ_i, φ),  y > 0

NB(y | μ, φ) = Γ(y+φ) / (Γ(φ) * y!) * (φ/(μ+φ))^φ * (μ/(μ+φ))^y

这里面三个参数我分别解释下：

π_i（零膨胀概率）：这个参数决定了spot i 产生结构零的概率。SPARK会通过一个logistic链接函数把这个概率和空间位置相关的协变量联系起来。π_i 越大，说明这个spot越可能产生结构零。

μ_i（均值参数）：给定基因在spot i 实际表达的情况下，表达量的期望值。这里SPARK会通过一个对数线性模型，把μ_i分解成空间模式和非空间基线的组合。

φ（离散参数）：控制负二项分布的方差，决定了数据相对于均值的离散程度。φ 越小，数据越容易产生极端高表达值——这在空间转录组里很常见，有些hotspot的表达量会高出平均水平一个数量级。

说实话，我自己早期写ZINB的似然函数时踩过不少坑。最大的一个就是参数初始化不对，导致EM算法收敛到局部极值。后来我养成了一个习惯：先用标准负二项模型跑一遍，把估计出来的μ和φ作为ZINB的初始值，π初始设为0.2。这个小技巧帮我避免了很多收敛问题。

ZINB vs 标准NB：在空间转录组中有多重要？

直接看对比数据。我拿一个真实的空间转录组数据集做过测试，选了1000个基因同时用NB和ZINB做拟合：

注意看AIC那一行。AIC越小模型越好，ZINB的AIC比NB低了122——这在模型选择里是一个巨大的差距。而且ZINB检出了更多空间变异基因，同时假阳性率降低了74%。

为什么会差这么多？因为标准NB模型把那些因为spatial dropout造成的零值当成了真实的生物学零，模型被迫用更大的离散参数去“解释”那些零，结果就是统计检验的能力大打折扣。ZINB通过引入π参数，把零的来源分开建模，离散参数φ就能更准确地反映真实的生物学变异性。

参数估计和收敛性

SPARK使用EM算法来估计ZINB的三个参数。具体的更新公式我不展开讲（论文里有详细推导），但有几个实战中的坑我需要提醒你：

第一个坑：EM收敛慢

EM算法在空间转录组数据上收敛确实慢。我遇到过跑40个迭代还在小幅震荡的情况。后来我加了一个early stopping：当对数似然的提升小于0.01%时提前终止。一般10-15个迭代就能收敛。

第二个坑：φ参数的边界问题

离散参数φ如果太小，模型会变得不稳定。我一般把φ的下界设为0.001，上界设为100。如果模型拟合后φ卡在下界上，说明这个基因的过离散程度比模型能处理的还严重，这时候我会换一个更灵活的模型——SPARK后面会用到广义Pareto分布来处理这种情况。

第三个坑：π参数的可识别性

如果基因的表达量整体很低（比如平均每个spot不到0.1个count），ZINB很难把结构零和采样零区分开。我定义了一个过滤器：如果一个基因在超过95%的spot上都是零，我会直接把它标记为“非表达基因”，不参与ZINB拟合。这样既省计算时间，也避免了模型在数据不足时的胡乱估计。

ZINB在SPARK中的具体实现

SPARK的源码里，ZINB的拟合是由一个Rcpp函数完成的。这里我把核心的计算逻辑摘出来，虽然不是完整的代码，但你可以看到关键步骤：

// 核心：ZINB的EM算法单次迭代
// 输入：counts, spatial_coords, current_params
// 输出：更新后的params

Step 1: E-step
  for each spot i:
    // 计算后验概率：该spot属于“零膨胀组”的概率
    posterior_zero[i] = pi[i] / (pi[i] + (1 - pi[i]) * NB(0 | mu[i], phi))
    
Step 2: M-step
  // 更新pi参数：用逻辑回归拟合空间协变量
  pi_new = logistic_regression(spatial_coords, posterior_zero)
  
  // 更新mu参数：用加权泊松回归拟合空间模式
  weights = 1 - posterior_zero
  mu_new = weighted_poisson_regression(spatial_coords, counts, weights)
  
  // 更新phi参数：通过一维优化最大化似然
  phi_new = optimize_over_dispersion(counts, mu_new, pi_new, phi_old)

Step 3: 检查收敛
  if |logLik(new) - logLik(old)| < 1e-4:
    break
  else:
    repeat Step 1-2

这个EM算法的复杂度是O(n·k·t)，其中n是spot数，k是空间基函数数量，t是迭代次数。对于10x Visium的4000个spot，一次完整的ZINB拟合大概需要0.5-2秒——这个速度是可以接受的。

SVG：ZINB模型的概率生成过程

下面这张图展示了ZINB如何将稀疏计数数据分解为两个生成路径：

一个实战场景：你该什么时候用ZINB？

如果有人问我“是不是所有空间转录组数据都该用ZINB”，我会说“不一定”。下面是我总结的几个判断标准：

• 零值比例 > 60%：建议用ZINB。这个阈值是经验的，我用过多个数据集验证过。低于这个比例，标准NB模型的误差可以接受。
• 检测到的基因数量 > 5000：每个基因都跑ZINB会比较慢。我会先用一个快速筛选（比如简单的Moran's I）把明显没有空间模式的基因过滤掉，只对候选基因做ZINB拟合。
• 组织异质性高：像脑、肿瘤这种有明确空间结构的组织，ZINB的π参数能捕捉到不同区域的技术差异，效果显著提升。
• 高分辨率数据：Slide-seq、MERFISH这类高分辨率数据往往更稀疏，零值比例可以到90%以上。这种情况下ZINB几乎是必须的选择。

我记得有一次帮一个同事处理他的MERFISH数据，他用标准方法只检出了15个空间变异基因，换用ZINB后检出了47个。他跑来问我“是不是我的数据有问题？”我说不是，“是你的模型选错了”。

参数设置的工程经验

在SPARK的实现中，有几个参数你可以根据数据特点调整。我列一个我自己用的推荐配置：

最后说一句，ZINB不是银弹。如果一个基因的表达模式根本不符合负二项分布（比如有多个离散的表达簇），ZINB的拟合也会有问题。但作为空间转录组数据分析的基础模型，它已经能够处理绝大多数实际数据中的零膨胀问题。SPARK后面会引入更复杂的空间协方差结构来提升检测能力，但ZINB是地基——地基打好了，上层建筑才能稳。