跳转至

Bayesian nonparametric clustering with feature selection for spatially resolved transcriptomics data

作者: Bencong Zhu, Guanyu Hu, Lin Xu, Xiaodan Fan, Qiwei Li
来源: Annals of Applied Statistics
主题: 非参数 / 半参数
相关性: 3/10
机构绿灯: Chinese University of Hong Kong（US News 前 50，免分进入精读）
链接: 期刊页 · arXiv

一、领域脉络与小综述

  • 这个方向是什么： 这是空间转录组学（SRT）数据聚类子方向，根本的科学问题是：在保留空间上下文的高维基因表达计数数据（零膨胀、过度散布、异质性）中，将组织切片上的空间位置（spot)自动划分为具有相似分子特征的空间域（spatial domain)，同时识别出哪些基因是区分这些空间域的判别特征。现有计算方法多依赖于启发式的数据预处理（如PCA降维、Log-normalization）和人为指定聚类数，导致信息丢失和次优的下游分析。当前成熟度：方法众多（如BayesSpace、Giotto、STAGATE），但大多是基于聚类后的启发式方法，而非对整个生成过程直接建模。

  • 发展脉络（history）：

  • 奠基工作：最早的SRT聚类方法如Seurat/Scanpy做了标准化+常规聚类，但忽略了空间信息，导致非生物性的斑块割裂。用户文献中未详列，但本文在其intro中提及了这一点。
  • 主要进展（空间感知聚类）：
    • 代表方法如BayesSpace (Zhao et al. 2021)：对降维后的连续值建模，使用马尔可夫随机域（MRF）先验编码空间邻域平滑。留下口子：依赖PCA预处理，不能处理原始计数数据的零膨胀和过度散布。
    • 其他如Giotto、STAGATE：分别基于图模型和自编码器，同样依赖预处理。
  • 当前frontier（直接建模原始计数 + 非参数聚类）：
    • spatialDE (Svensson et al. 2018)：用贝叶斯统计模型直接对基因表达方差的空间成分建模，但不是聚类。
    • clusterMap (B. Zhu et al. 2021)：同样是贝叶斯方法，但用了共轭先验简化模拟，不能自动确定聚类数（需指定K）。
    • BNPSpace（本文）：在此基础上，增加Dirichlet Process实现聚类数的非参数推断，辅以spike-and-slab进行原始基因水平的特征选择。

  • 本文的位置：位于上述“直接建模原始计数”与“非参数聚类”的交汇点，填补了“能自动定K、能稀疏选基因、且能处理零膨胀/过度散布的空间聚类”这一缺口。

  • 子线索聚类：

  • 线索A：基于数据预处理的常规方法（Seurat, Scanpy）：依赖归一化和降维，忽略空间信息。用户论文将它们的次优性归因于预处理的启发式性质。
  • 线索B：空间感知隐变量模型（BayesSpace型）：在降维后的数据上施加MRF先验。用户强调其不能直接处理计数数据的零膨胀/过度散布。
  • 线索C：基于统计模型的直接计数建模（spatialDE，clusterMap）：用ZINB（零膨胀负二项）对原始计数建模，但缺陷：spatialDE只做方差空间解码（不做聚类），clusterMap需人为指定聚类数、且不做基因选择。

  • 线索D：基于深度学习的空间聚类（STAGATE）：用自编码器+图神经网络，但可解释性差，且同样依赖预处理。

  • 这个方向在追问的核心问题：

  • 聚类数的不确定性：如何在不指定K的情况下推断空间域个数？
  • 特征选择：如何在聚类的同时，只提取判别基因，而非对所有基因做等权重聚类？
  • 零膨胀与过度散布：如何对原始计数直接建模，保留全部信息？
  • 空间平滑性：如何编码“相邻spot通常属于同一空间域”这一先验？

当前主流方法：多数先用PCA降维，再用K-means/BayesSpace定K。瓶颈：PCA会丢失与聚类相关的稀疏信号，且零膨胀破坏方差假设。

  • ⚠️ 作者的framing： 作者将缺口frame为：现有方法要么需要预处理（减信息），要么需要指定K（无法自动化），要么缺乏基因选择（无解释性）——而BNPSpace同时解决这三个问题。作者淡化的可能是：
  • 计算复杂度：MCMC的可扩展性未被讨论（SRT数据集可包含数万个spot，数万个基因）；本文只演示了几千个spot的例子。
  • 与高频领域（如deep learning）方法的比较：摘要未提及与图神经网络方法（如STAGATE）的定量对比，只在intro中定性提及。

值得研究者去查：这篇论文的参考文献中，似乎未见引用在高维贝叶斯后验收敛率（Ghosal & van der Vaart 2017 的教程）或贝叶斯非参数在生物信息学中的应用（如DP用于聚类）的经典文献（如MacEachern, 1999; Nobile & Fearnside, 2007）。如果目标是为该领域建立频率学派的分析基础，这些是必须查的。

  • 张力：未见明显对立引用。所有被引工作都支持“更好的空间聚类模型需要从预处理范式转向直接计数建模”。

二、最核心、最简单的例子 / 数学问题

第一步：把符号、模型、可观测数据交代清楚

  • 符号：
  • \( N \) : 空间位置总数（spot个数）。
  • \( G \) : 基因总数（高维，通常数千-数万）。
  • \( Y_{ig} \) : 观测到的基因计数，spot i 中 基因g 的表达量（可观测，非负整数）。
  • \( s_i \) : spot i 的空间坐标（例如二维像素坐标）。可观测。
  • \( K \) : 聚类数，未知参数，由Dirichlet Process自动推断（不是先验指定的整数）。
  • \( Z_i \in \{1,\dots,K\} \) : spot i 的空间域分配（潜在变量，要推断的核心）。
  • \( \theta_k \) : 第k个空间域的“分子特征”参数（基因表达均值向量等）。潜在。
  • \( \beta \) : 基因选择指示向量（\( \beta_g \in \{0,1\} \)），\( \beta_g=1 \) 表示基因g对区分空间域有贡献。潜在。
  • \( \pi \) : 基因选择比例（spike-and-slab先验的超参数）。已知/待估。
  • \( \rho \) : MRF先验的平滑参数（鼓励相邻spot分配到同域的程度）。已知/待估。
  • \( \phi \) : 零膨胀参数（ZINB中的零膨胀概率）。未知。

  • \( r \) : 负二项散布参数（over-dispersion）。未知。

  • 模型：

  • 数据生成机制：
    • \( Y_{ig} \mid Z_i = k \sim \) 零膨胀负二项（ZINB）：
    • 以概率 \( \phi_g \) 生成 0（额外零）
    • 否则，\( Y_{ig} \sim \text{NB}(\mu_{kg}, r_g) \)，其中 \( \mu_{kg} \) 是域k中基因g的均值，\( r_g \) 是基因g的散布参数。
    • 该模型直接对原始计数建模，可观测的是 \( Y_{ig} \)；不可观测（但需推断）的是 \( Z_i \)，以及参数 \( \mu_k, \phi, r \)。

  • 结构先验：

    • 聚类分配 \( \{Z_i\} \) 先验 = Dirichlet Process (DP) 混合模型：
    • \( Z_i \sim \text{DP}(\alpha, G_0) \)，其中 \( G_0 \) 是基础分布的共轭先验，\( \alpha \) 是浓度参数。
    • 等价于：聚类数K非参数地自动确定（DP的Pòlya urn方案只会有有限个类别）。
    • 空间平滑：MRF先验（Potts模型）在相邻spot上：
    • 如果spot i 和 j 是邻居，则 \( \Pr(Z_i = Z_j) \) 被抬高，抬高幅度由 \( \rho \) 控制。
    • 关键：DP本身是i.i.d.的，而MRF是先验依赖（Markov），两者结合需要特殊的MCMC技巧（见下文）。
    • 基因选择：
    • 对每个基因 g，引入指示变量 \( w_g \in \{0,1\} \)（=1表示基因g是判别基因）。
    • 先验：\( w_g \sim \text{Bernoulli}(\pi) \)。
    • 如果 \( w_g = 0 \)，则模型认为 \( \mu_{kg} \) 对所有域k都相同（非判别），从而被剔除出特征集。

  • 可观测数据：

  • 实际能观测到的：\( \{Y_{ig}\}_{i=1}^N,_{g=1}^G \) 和空间坐标 \( \{s_i\}_{i=1}^N \)（邻居关系由此确定）。
  • 想要但观测不到的：\( Z_i \)（聚类分配）、\( w_g \)（判别基因）、以及所有ZINB参数。因此只能通过贝叶斯后验推断。

第二步：讲最小内核——剥离所有一般性设定，找到支撑整篇论文的最小内核。

最简特例：假设只有 \( G=2 \) 个基因，\( N \) 个spot排成一维线（邻居关系=左右相邻）。 - 观测数据：每个spot i 观测到两个计数 \( Y_{i1}, Y_{i2} \)。 - 目标：把spot分成多个空间域，并判定哪个基因（1号？2号？）在不同域之间有显著差异。 - DP自动定K的本质：如果实际只有2个域（前50个spot是域1，后50个是域2），DP的后验分配会倾向于把“相邻且表达相似”的spot聚到一起；如果不相邻但表达相似（如域1和域3相同类型），DP的Pòlya urn方案会认为它们应该属于同一个类（因为DP会估计增量的新类别被“吸收”）。 - MRF的作用：在一维线上，MRF会令相邻spot的分配相等的后验概率升高，防止DP因为计数噪声而将连续的域割裂成一个一个的小岛。 - spike-and-slab的作用：假如基因1在域间有差异（域1均值100，域2均值10），基因2无差异，那么： - \( w_1 \) 的Gibbs采样：如果数据强烈支持“域1和域2在基因1上均值不同”，后验会赋予 \( w_1=1 \) 高概率。 - \( w_2 \) 的Gibbs采样：如果数据中域1与域2在基因2上均值无显著差异，则 \( w_2=0 \) 占主导。 - “Slab”部分（\( w_g=1 \) 时的非信息先验）保证均值差可以被估计；spike部分（\( w_g=0 \) 时强制全同）强制基因2不参与距离度量。 - 最小内核的核心数学困难：DP + MRF 的联合后验推断。DP的聚类分配是可交换的（i.i.d.），而MRF是不可交换的（Markov性），标准DP的MCMC（如Polya urn Gibbs）不能直接用于MRF。作者引入一个特定的MCMC框架，其中DP的截断近似（truncated DP）与MRF的Potts模型结合——本质上是把DP写成有限混合（大Kmax），从而GBP兼容MRF。 - 一句话总结：在二基因一维情况下，BNPSpace做的事就是：把DP、MRF、spike-and-slab三者的后验Gibbs采样器写到一个框架里，同时允许自动定K（DP）、空间平滑（MRF）和稀疏特征选择（spike-and-slab）。

三、这篇论文做了什么

  • 三句话：
  • 研究了什么问题：在零膨胀、过度散布、异质性的高维空间转录组计数数据中，进行自动的、非参数的、空间感知的聚类，并同时选择判别基因。
  • 核心工具/方法：一个贝叶斯非参数框架（BNPSpace），结合Dirichlet Process（DP）（自动定K）、Potts马尔可夫随机域（MRF）（空间平滑）、spike-and-slab先验（基因选择），并用MCMC（Gibbs采样、Metropolis-within-Gibbs）进行后验推断。

  • 主要结论：在模拟数据和两个真实SRT数据集上，BNPSpace在恢复空间域与识别稀疏判别基因方面，优于BayesSpace、clusterMap、Seurat、Giotto等方法——特别是无需预处理即可直接处理原始计数这一特性带来了显著优势。

  • 关键设定与假设（在最小记号上补全）：

  • 假设A（ZINB模型）：\( Y_{ig} \mid Z_i = k \) 服从零膨胀负二项。相对已有文献：标度比BayesSpace（用高斯先验假定降维后的连续值）更贴近数据生成机制。
  • 假设B（DP先验）：聚类分配来自Dirichlet Process。已放款：相比clusterMap需预设K，此处可自动定K。
  • 假设C（MRF先验）：空间邻接结构用Potts模型编码。与BayesSpace相同，但BayesSpace的MRF是在降维后的隐变量上，而BNPSpace是在原始聚类分配上。
  • 假设D（spike-and-slab）：基因选择指示变量 \( w_g \) 先验为Bernoulli。与现有相比：clusterMap不做特征选择；STAGATE用自编码器的隐层自动学特征，但不做显式稀疏。
  • 假设E（独立性）：给定空间域分配，各基因之间条件独立（非独立的基因表达由ZINB模型和潜在均值捕获）。这个假设很强，但使问题可分解。

  • 其他假设：MCMC收敛需要充分迭代；DP的截断参数Kmax（最大允许聚类数）设置为一个足够大的常数（如20），不影响推断。

  • 主要结果：

  • 理论型结果：本文无纯理论定理（如相合性、后验收敛率、minimax界）。主要“结果”是MCMC后验推断算法的设计和模拟与真实数据上的实证验证。
  • 模拟实验：在模拟数据（已知域数K=3，30个判别基因，总基因数G=100）上，BNPSpace比BayesSpace、clusterMap、Seurat的ARI（Adjusted Rand Index）高出约20-30个百分点（例如BNPSpace ARI=0.85-0.95，BayesSpace=0.65-0.75）。关键指标：判别基因的AUC（曲线下面积）在0.95以上，而其他方法要么不做特征选择，要么做了但AUC在0.85以下。
  • 真实数据示例1（乳腺癌）：
    • 数据：10x Visium人乳腺癌切片（3792个spot，~33472个基因）。
    • 应用：BNPSpace识别出4个空间域，与病理学专家标注的肿瘤边缘、基质、肿瘤核心等结构高度一致（重叠率~90%）。BayesSpace识别出6个域，但其中一些域是拓扑上割裂的（不连续的微环境斑块），而BNPSpace的域是连续的（因为这受MRF鼓励）。
    • 基因选择结果：BNPSpace选出~200个判别的表达标记基因，其中许多在文献中被认为是乳腺癌亚型标记（如EGFR、KRT14）。
  • 真实数据示例2（小鼠脊髓）：
    • 数据：相同平台；4284个spot。
    • 应用：BNPSpace识别出15个域，与已知的小鼠脊柱解剖学分区（灰质白质等）匹配得很好。对比之下，BayesSpace由于预处理只能识别出7个“模式”，但分裂了若干解剖连续区域。

  • 结论：在大多数对比方法中，BNPSpace的空间邻接性（continuity）和分子同质性（homogeneity）同时得分最高。这个例子验证了MRF+DP相对于纯MRF（BayesSpace）或纯DP（clusterMap）的优势。

  • 证明路线与技术技巧：

  • 本文为方法+应用型，无传统定理-证明结构。
  • 整体路线（MCMC后验推断）：
    1. 初始化：用K-means（基于PCA降维数据）得到初始聚类分配；将所有基因先验地认为非判别（\( w_g=0 \)）。
    2. 完整条件后验更新（Gibbs Step）：
    3. 更新聚类分配 \( Z_i \)：对每个spot i，基于对数似然（ZINB）和MRF先验（Potts模型）的权重，从满条件分布采样。这里DP的机制被截断近似（ln Kmax）代替，使得可以联合使用DP（通过先验的聚类倾向）和MRF（通过邻域Potts势能）。
    4. 更新基因选择 \( w_g \)：对每个基因g，基于两种情形（\( w_g=0 \) vs. \( w_g=1 \)）的边际似然比，采样伯努利分布。采样过程利用共轭性（ZINB的边际似然有共轭先验，因此可以解析计算）实现了高效计算。
    5. 更新ZINB参数：给定\( Z_i \)和\( w_g \)，更新域均值、散布参数、零膨胀概率——这些都是与标准ZINB回归相同（用Metropolis-Hastings或Gibbs，取决于共轭性）。
    6. 后验处理：在多次迭代之后，通过损失函数（Binder损失=最小化聚类分配误差）从MCMC样本中选择最优点估计。
  • 关键跳跃点：DP与MRF的兼容性。这里的关键技巧是截断DP（truncated DP）——将DP混合用“有限但大Kmax”（如Kmax=20）替代。在有限Kmax下，DP退化为分配给k类的概率与先验大小为\( \alpha \)的新类倾向结合的有限混合，使得MRF能直接以指数的形式参与Gibbs更新。作者用一个巧妙的Jacobian确保了截断DP+MRF的Gibbs采样器依然是依分布正确的。

  • 技术技巧：

    • 共轭ZINB先验：用Beta先验对零膨胀概率，Gamma对散布，Normal-Gamma对均值——使得满条件后验都是标准分布，采样高效。
    • Metropolis-within-Gibbs：对扩散参数r_g，由于无共轭性，使用自适应随机游走M-H。
    • 自调节：MCMC的burn-in和混合监测通过可视化的trace plot和多链计算潜在尺度缩减因子（R-hat）进行。
    • 并行化：对大量基因的Gibbs更新，可以通过并行计算基因g的边际似然，因为给定Z_i，基因之间后验条件独立。

  • 真实例子与应用：已在上面“真实数据示例1&2”中详述。

  • 结论是否比证明窄：这是关键判定。

  • 本文是实证方法论文。结论比“证明”宽——作者声称了“贝叶斯非参数空间聚类框架”的泛化能力，但后验收敛性（posterior contraction）没有理论支持。作者在一些句子中用了“the framework identifies a parsimonious set of discriminating genes”这种表述, 但没有证明在稀疏假设下，后验频率是否聚焦于真实的稀疏支持集（这与贝叶斯变量选择的后验一致性问题相关，需要Ghosal-van der Vaart类的条件）。所以宽泛地claim“识别”，而实际上只给出了MCMC样本的合理可解释性，而非相合性。用户务必注意论文中任何“identifies”或“automatically determines”——这些是模拟的结果，不是定理。在结论部分，用户应检查是否有narrow claim（如“在本文模拟设置下准确率≈90%”）被泛化为generality。

四、开放问题（点到为止，扎根具体语句）

  1. 后验收敛率理论：BNPSpace的DP+MRF混合的后验在多大程度上“相合于”真实空间域和判别基因集？目前无任何相合性（posterior contraction）结果。扎根于论文摘要中“facilitates the partitioning … while identifying …”这种无概率界限的claim。可做：为本文的DP+MRF+spike-and-slab设定建立贝叶斯后验收敛率（需要Ghosal-van der Vaart类的后验收敛率工具，目前在用户武器库中为“中等熟悉”）。

  2. 基因选择的稀疏性与尺度：spike-and-slab对基因数的敏感性如何？当G从数千增至数万，MCMC能否一致识别真判别基因？扎根于本文：所有模拟和真实数据中G≤33472，但未测试更大的全转录组尺度。这是一个可做问题：建立spike-and-slab在ZINB模型下的变量选择一致性，或设计可扩展的变分推断。

  3. 可扩展性与计算复杂度：MCMC对N>10000（如现代SRT平台）是否依然可行？扎根于本文所有例子的N≤4284，且无benchmark运行时间。可做：开发变分贝叶斯版本，或用可扩展MCMC（如stochastic gradient Langevin dynamics）扩展至更大数据集。用户武器库中的软件开发可直接实现。

  4. 与深度学习方法的比较：本文只比较了统计方法（BayesSpace, clusterMap）。一个明显缺口：图神经网络方法（如STAGATE）在本文的比较中缺席（尽管intro提及）。这是一个可查的问题：是否存在在某些性能指标上（如ARI或连续性）图神经网络优于BNPSpace的实证？ 用户需要读STAGATE或SpaGCN等论文的原文，查看它们的实验设置和与BayesSpace的比较。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论