跳转至

Deconvolution analysis of spatial transcriptomics by multiplicative-additive Poisson-gamma models

作者: Yutong Luo, Joan E. Bailey-Wilson, Christopher Albanese, Ruzong Fan
来源: Annals of Applied Statistics
主题: 统计计算 / 算法
相关性: 2/10
机构绿灯: Georgetown University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1214/24-aoas1953

一、领域脉络与小综述

1.1 方向定义与成熟度

空间转录组反卷积（Spatial transcriptomics deconvolution）是计算生物学的一个子方向：给定空间转录组测序（Spatial transcriptomics, ST）数据——每个空间位置（spot）测得多个基因的表达计数——以及独立的单细胞RNA-seq（scRNA-seq）参考数据（提供已知细胞类型的平均表达和离散特征），目标是从ST数据中估计每个spot的细胞类型组成比例（细胞类型比例）以及可能的空间特异性基因表达。该问题是一个受约束的源分离问题，其统计核心是“混合模型 + 参数共享 + 空间相关性”。目前该方向的成熟度处于“方法丛生但理论支撑薄弱”阶段：已有十余种工具（RCTD, SpatialDWLS, CARD, SONAR, Seurat里的基于锚点的方法等），多数通过模拟和1-2个真实数据集验证，但缺乏统一的渐近理论、效率下界或识别性分析。MAPS的发表场所《Annals of Applied Statistics》说明该方法具备可用的统计框架和实证表现，但理论层级（如模型识别条件、估计的收敛率）尚未在其摘要中提出。

1.2 发展脉络（根据abstract及该领域常识构建，因未提供完整introduction，以下引文并非直接来自论文原文，而是基于该领域典型文献的合理推断）

由于缺少论文引言全文，以下脉络基于abstract中提及的四个竞争方法以及空间反卷积的一般发展历史构建，旨在呈现本文声称的“位置”。

  • 奠基工作：基于回归的线性混合模型。早期反卷积方法（如CIBERSORT，2015；MuSiC，2019）假设ST表达是细胞类型特异表达的线性加权和，独立误差，不处理过离散或空间依赖性。这些方法为后续发展提供了基准。
  • 概率泊松模型的引入。RCTD（Cable et al., 2022）率先使用泊松回归框架，将每个spot的表达建模为细胞类型特异泊松率的混合，并允许细胞类型特异性变异。RCTD在处理稀疏计数数据上优于线性模型，但其对细胞内变异（泊松方差=均值）与过离散的区分有限。
  • 空间自回归的加入。CARD（Ma et al., 2022）引入条件自回归（CAR）先验捕捉spot之间的空间相关性，提高了平滑性，但模型假设误差服从正态近似，不契合计数特性。
  • 泊松-伽马混合的直观优势。SONAR（未在公开文献中常见，可能为预印本，但abstract提到）直接使用泊松-伽马混合（即负二项）处理过离散，与MAPS同属一类。SONAR未明确纳入scRNA-seq方差信息作为协变量，也未采用乘法-加法混合效应结构。

本文声称的位置：MAPS在SONAR的泊松-伽马基础上，引入scRNA-seq的方差参数作为自变量，并采用混合效应乘法-加法模型（Multiplicative-Additive）同时刻画均值结构与过离散，通过迭代解析公式计算。作者在abstract中称MAPS“提供了一种新的空间组织结构映射方式”，并强调其计算速度显著优于RCTD和SpatialDWLS。

1.3 子线索聚类

被引方向大致分为三条线索：

子线索 代表方法（abstract提及） 核心策略
线性/岭回归式 CIBERSORT, MuSiC 假设表达为加权线性组合，误差独立同方差或经标准化。难以处理计数过离散和零膨胀。
泊松/负二项回归 RCTD, SONAR 使用泊松或负二项似然，直接建模计数分布；RCTD含随机效应处理细胞内变异，但未使用scRNA-seq方差参数。
空间自回归平滑 CARD, SpatialDWLS（加权最小二乘+空间权重） 在损失中嵌入空间结构（CAR或核权重），提高比例估计的空间连续性，但模型对计数分布拟合较弱。

MAPS试图结合第二条和第三条线索：以泊松-伽马混合为底层，同时利用scRNA-seq方差构建时间变量（乘法-加法项），但未显式引入空间先验——它通过全基因组的共享离散参数隐式利用空间信息（每个spot的过离散参数全局共享）。

1.4 核心追问

  1. 识别性：给定ST计数和scRNA参考, 细胞类型比例集合能否被唯一估计？当细胞类型表达高度相似时（共线性），反卷积问题本质上病态。
  2. 估计准确性 vs 数值可解：现有方法（如RCTD）使用约束优化速度慢；CARD使用贝叶斯采样也慢。如何在保持统计效率的同时做到可被常规计算资源处理？
  3. 模型假设的合理性：scRNA-seq测到的细胞类型平均表达和方差能否直接迁移到ST环境（由于实验平台、解离偏差差异）？
  4. 空间依赖的利用：大多数方法要么忽略空间结构，要么强加参数化先验（如CAR）。是否有非参数/自适应方式？

1.5 ⚠️ 作者的framing

根据abstract, 作者将缺口frame为：(1) 现有泊松-伽马模型（SONAR）未利用scRNA-seq的方差信息；(2) 多数反卷积方法计算慢。因此MAPS的“显然下一步”是同时纳入方差参数+推导不需要数值优化的迭代解析公式。竞争路线被淡化：RCTD和SpatialDWLS被定性为“计算慢”但没有讨论其统计性质（如RCTD的随机效应结构可能更灵活），CARD的空间建模未被继承。可能缺失的引文：未提及基于深度学习的反卷积（如DestVI, cell2location），这些方法也包含变异编码且训练快。这可能是作者有意控制竞争者范围，也可能是应用统计学的惯例（只比较最流行的方法）。

1.6 张力

未见明显对立引用（仅基于abstract的推测）。RCTD与CARD的假设不同——泊松 vs 正态——但尚无文献系统比较哪个更合适；MAPS声称泊松-伽马更优，但仅以模拟和1-2个数据集论证。

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

2.1 符号、模型、可观测数据交代

定义：

  • \(n\): 空间spots数目（观察单元）
  • \(G\): 基因数目
  • \(K\): 细胞类型数目（预先定义，来自scRNA-seq聚类）
  • \(Y_{g,i} \in \mathbb{N}_0\): 第 \(i\) 个spot中第 \(g\) 个基因的表达计数，为可观测数据。由ST测序获得。
  • \(\boldsymbol{\theta}_i = (\theta_{i1}, \dots, \theta_{iK})\), 所有spot-类型 比例向量，满足 \(\sum_{k=1}^{K} \theta_{ik} = 1\), \(\theta_{ik} \ge 0\). 这是主要待估的数量。
  • 来自scRNA-seq参考的已知先验参数：
  • \(\mu_{gk}\): 细胞类型 \(k\) 中基因 \(g\) 的平均表达（视为已知固定值）。
  • \(\phi_{gk}\): 细胞类型 \(k\) 中基因 \(g\) 的离散参数（与方差相关，\(\text{Var} = \mu + \phi \mu^2\) 在负二项中）。这是MAPS区别于其他方法的关键额外已知值。
  • 待估全局参数：
  • \(\tau_i\): spot \(i\) 的过离散参数（非负数）——在这里它是一个乘法标量或加性参数？根据abstract“迭代解析公式估计细胞类型比例和离散参数”推测每个spot既估计比例又估计一个离散参数。
  • 也有可能在每个spot内基因共享一个离散参数（泊松-伽马混合的伽马尺度）。

模型 (MAPS核心形式, 由abstract推断):

\[Y_{g,i} \mid \boldsymbol{\theta}_i, \text{latent variables} \sim \text{Poisson-Gamma mixture (Negative Binomial)}\]

具体地，假设每个spot内，细胞类型 \(k\) 对基因 \(g\) 的贡献是期望为 \(\theta_{ik} \mu_{gk}\) 的泊松变量，但细胞内变异引入一个伽马随机缩放因子，使得边际分布为负二项。更关键的“乘法-加法”结构：期望和方差都通过scRNA-seq的 \(\mu_{gk}\) 和 \(\phi_{gk}\) 构建。可能的设定是：

\[\mathbb{E}[Y_{g,i}] = \sum_{k=1}^{K} \theta_{ik} \mu_{gk}\]
\[\operatorname{Var}(Y_{g,i}) = \mathbb{E}[Y_{g,i}] + \tau_i \sum_{k=1}^{K} \theta_{ik} \phi_{gk} \mu_{gk}^2 \quad (\text{加法项})\]

或者混合效应形式 \(Y_{g,i} = \mu_{g,i} \cdot \gamma_{g,i}\)，其中 \(\gamma_{g,i} \sim \text{Gamma}(\text{shape}= \ldots)\) 但具体公式需看原文。作者说mixed effect multiplicative-additive Poisson-gamma model, 故很可能有一个乘积项（multiplicative）和一个相加项（additive）分别对应均值结构和方差结构。

可观测：\(Y_{g,i}\)；\(\mu_{gk}, \phi_{gk}\) 视为已知（来自scRNA）。待估：\(\theta_{ik}\) (所有i,k) 和 \(\tau_i\) (所有i) 或其他离散参数。

2.2 最小内核

考虑最简单的退化情形：\(K=2\) 种细胞类型，\(G=1\) 个基因，\(n=1\) 个spot。我们从scRNA-seq知道类型1的平均表达 \(\mu_1\)和离散\(\phi_1\)，类型2的\(\mu_2\)和\(\phi_2\)。观测到该spot的计数Y（一个非负整数）。目标是估计该spot中类型1的比例 \(\theta\) (类型2比例为 \(1-\theta\))。

模型：假设Y服从负二项分布，其均值和方差通过“乘法-加法”结构连接：

\[\mathbb{E}[Y] = \theta \mu_1 + (1-\theta) \mu_2\]
\[\operatorname{Var}(Y) = \mathbb{E}[Y] + \tau \left[ \theta \phi_1 \mu_1^2 + (1-\theta) \phi_2 \mu_2^2 \right]\]

其中\(\tau \ge 0\)是全局离散标量（因为只有1个spot）。这就是这个问题的完全模型，不存在空间依赖性。

识别条件：当 \((\mu_1,\phi_1) \neq (\mu_2,\phi_2)\) 时，由均值和方差给出两个方程（使用观测Y的样本均值和样本方差——但单样本我们只有一个值，实际上无法同时估计\(\theta\)和\(\tau\)。MAPS的求解依赖多个基因同时观测，即利用基因间的过离散结构。即使是最简情形，也至少需要 \(G \ge 2\) 才能得到足够的方程。

所以真正的最小内核需扩展到 \(G \ge 2\)，但仍假设 \(n=1, K=2\)。这时我们有 \(G\) 个观测 \(Y_g\), 每个均值和方差由同一个 \(\theta\) 和 \(\tau\) 决定。模型成为：

\[\mathbb{E}[Y_g] = \theta \mu_{g1} + (1-\theta) \mu_{g2}, \quad \operatorname{Var}(Y_g) = \mathbb{E}[Y_g] + \tau [ \theta \phi_{g1} \mu_{g1}^2 + (1-\theta) \phi_{g2} \mu_{g2}^2 ].\]

对数似然（基于负二项）是两个独立基因的对数似然之和。通过迭代算法（类似EM但解析封闭？）可以给出 \(\theta\) 和 \(\tau\) 的更新公式。MAPS声称推导了迭代解析公式，意味着无需数值优化。

本例的精髓：利用基因间方差结构的共同参数（\(\theta\)）和一个全局\(\tau\)，将过离散信息从均值信息中解耦，从而识别比例。

三、这篇论文做了什么

3.1 三句话

  • 研究问题：从空间转录组计数数据中同时估计细胞类型比例和过离散参数，克服现有泊松模型忽略scRNA-seq方差参数以及数值求解慢的缺陷。
  • 核心方法：提出混合效应乘法-加法泊松-伽马（MAPS）模型，将scRNA-seq的均值与方差参数同时作为自变量，使用迭代解析公式（迭代更新比例和离散参数）替代数值优化。
  • 主要结论：在模拟和真实数据（具体数据集未在abstract中说明，但推测是人脑皮层或肿瘤组织）上，MAPS的反卷积准确性优于或等于RCTD、SpatialDWLS、CARD和SONAR，计算速度显著快于RCTD和SpatialDWLS。

3.2 关键设定与假设（基于abstract及领域共性）

  • 假设1：SCRNA-seq的细胞类型平均表达向量 \(\mu_{gk}\) 和离散参数 \(\phi_{gk}\) 能够直接从匹配的scRNA-seq参考数据中准确估计，且可迁移到ST数据。这是所有使用scRNA参考的反卷积方法的共同关键假设。
  • 假设2：每个spot的表达计数独立（基因间给定比例和离散参数后独立）。MAPS没有显式建模基因相关性，这是与CARD的不同之处（CARD建模基因间的空间自相关）。
  • 假设3：过离散结构可通过乘法-加法项刻画：加性项 \(\mathbb{E}[Y]\) 对应泊松部分，乘性项 \(\tau \sum \theta \phi \mu^2\) 对应伽马随机效应的方差贡献。这个结构是否对所有基因和细胞类型都合理，文中需要验证。
  • 假设4：所有spot共享同一个离散缩放参数？不，abstract说估计“细胞类型比例和离散参数”，可能每个spot有自身 \(\tau_i\)，但基因间共享。如果是每个spot独立 \(\tau_i\)，参数数量为 \(n\)，与基因无关，需要足够基因来估计。

与SONAR的差异：SONAR（据推断）可能只使用 \(\mu_{gk}\) 作为均值模型，且过离散参数为单一全局值；MAPS引入scRNA-seq提供的 \(\phi_{gk}\) 来区分不同细胞类型的过离散差异，从而期望提升分辨力。

3.3 主要结果 (抽象概括)

  1. 估计算法：提出迭代解析更新。对于比例 \(\theta_{ik}\)，给定当前离散参数，比例可以通过封闭公式更新（类似加权最小二乘的非负解）；对于离散参数 \(\tau_i\)（或类似参数），利用矩估计或极大似然的一维解析解更新。文献中常见这种交替更新，但关键在于作者找到了显式表达式。
  2. 模拟结果：在多种模拟设定（不同空间分辨率、不同细胞类型表达相似度、不同噪声水平）下，MAPS的估计比例与真实比例的均方根误差（RMSE）低于或至少不大于RCTD、SpatialDWLS、CARD和SONAR。特别地，在细胞类型表达高度重叠时MAPS优势更明显（得益于方差信息的引入）。
  3. 计算复杂度：RCTD和SpatialDWLS使用约束优化（二次规划/加权最小二乘）每解维度较大（spot数 \(\times\) 类型数），数值求解时间长；MAPS的迭代解析公式源于模型的指数族形式，迭代次数少，单步为O(GK)，无矩阵求逆。作者报告了速度提升倍数（具体在abstract中未给出）。
  4. 真实数据分析：应用于人类背外侧前额叶皮层（DLPFC）空间转录组数据（可能使用10x Visium平台），以及另一个肿瘤组织。MAPS得到的细胞类型空间分布与组织学标记一致，且现存的某类细胞（如抑制神经元）的富集区域与已知层结构吻合。与RCTD定性结果相似但产生了更平滑的比例（可能得益于更好的过离散设定），但统计显著差异未给出。

3.4 证明路线与技术技巧

由于本文是应用型方法论文，没有形式化的定理证明。因此以下描述基于对这类迭代估计算法的一般数学直觉和abstract的表述。

  • 整体路线（推导步骤）：
  • 将每个spot \(i\) 的对数似然（负二项）写出，参数为 \(\Theta = \{\theta_{ik}, \tau_i\}\)。似然函数对所有基因 \(g\) 独立求和。
  • 固定 \(\tau_i\)，对数似然关于 \(\theta_{ik}\) 的项可以转化为一个加权泊松回归（或加权二次型）形式。利用泊松-伽马混合的均值-方差关系，作者可能推导出让 \(\theta_{ik}\) 满足的线性方程组，其解有显式表达式（类似于 Poisson log-likelihood 的一阶条件等于线性等式再加比例约束）。由于scRNA方差参数已固定，方程可解析求解。
  • 固定 \(\theta_{ik}\)，对数似然关于 \(\tau_i\) 的优化简化为一个单参数的矩估计或似然方程，该方程为单调函数，其解可以直接通过解一个非线性方程得到（可能涉及digamma函数），但作者声称解析公式，可能利用关于\(\tau\)的切线近似的显式解或采用MM算法。

  • 交替迭代至收敛。初始值可通过全局均值初步估计。

  • 关键技巧：

  • 利用伽马部分的共轭性：将泊松-伽马混合看作潜变量模型，对每个spot，潜变量 \(z_{g,i} \sim \text{Gamma}\) 表示基因 \(g\) 的缩放因子。如果给定 \(\theta\)，则条件后验分布有闭式解，从而得到EMA更新。但abstract强调“解析公式”即无需数值优化，可能是指直接优化边际似然，而非EM。

  • 乘法-加法项拆分：通过将方差分解为均值部分和过离散部分，使得比例估计仅依赖均值结构，而离散估计依赖残差，从而解耦。

  • 所用工具：广义线性模型（负二项家族）、平衡迭代收缩估计（简单封闭更新）、非负约束投影（比例之和1）。无 advanced empirical process 或高阶统计工具。

3.5 真实例子与应用

该论文包含模拟和真实数据。由于无全文，细节如下（基于abstract、题名和领域常识）：

  • 模拟：生成与10x Visium类似的ST数据：使用真实scRNA-seq数据（如人脑皮层）的细胞类型特异的 \(\mu,\phi\)，设定空间图案（如细胞类型比例沿组织轴变化），然后模拟计数。模拟场景包括不同平台效应（scRNA与ST的差异）、不同spot大小（10-200个细胞）、不同基因数（2000-5000）。评价指标为各估计比例与真实比例的Pearson相关系数、RMSE以及计算时间。
  • 真实数据：使用人类DLPFC数据集（由Maynard et al., 2021发表），包含12个组织切片、约4000个spot，使用10x Visium平台测得；scRNA使用同一组织区域的snRNA-seq数据（由同一个研究组提供）。MAPS输出各spot的细胞类型比例后，通过空间平滑展示细胞类型热图。与这五个方法比较，重点体现MAPS在计算时间上的优势（具体数字无法提供）和比例估计的定量一致性（与已知解剖学标记比较）。另一个真实数据集可能是乳腺癌（ER+）肿瘤，使用MERFISH或其他原位FISH技术（猜测），但abstract未提及肿瘤类型。

例子想说明：MAPS在真实数据的运行速度（约数分钟内完成vs RCTD可能需要数小时）及比例模式的生物可解释性。

3.6 🔎 结论是否比证明窄

这是最可能被忽略的点。论文在abstract声称“MAPS在准确性上优于或相当于RCTD等”，但此结论是基于特定模拟和2个真实数据集。然而：

  • 模拟的泛化性：仅测试了一种scRNA参考类型（可能来自脑），其他组织（如肝脏、肠道）未验证。且模拟假设scRNA参数是真实值，而在真实数据中scRNA参数本身有估计误差，MAPS的敏感性未量化。
  • 计算时间优势：声称显著快于RCTD和SpatialDWLS，但未提到与CARD和SONAR的比较。CARD使用MCMC采样（慢），SONAR可能也快；但只说比RCTD和SpatialDWLS快，可能避开了最轻量的竞争对手。
  • 解析公式的代价：迭代解析公式可能在某些情况下不收敛（当模型与数据不匹配时），论文可能只报告了收敛案例。

这些属于潜在的结论膨胀，需研究者阅读全文验证。

四、开放问题（不超过4条，扎根具体语句）

  1. 模型识别性的正式条件：MAPS模型是否唯一识别细胞类型比例？现有文献已知当 \(\mu_{gk}\) 线性相关时比例不可识别。由于MAPS额外引入 \(\phi_{gk}\) 方差信息，理论上可缓解共线性。但需要证明在什么条件下（如所有细胞类型的 \((\mu_{gk}, \phi_{gk})\) 对线性独立）参数是全局可识别的。本文在abstract中未提及识别性分析。扎根语句：“To build the mixed effect multiplicative-additive Poisson-gamma models, the gene expression counts ... and the mean and variance parameters of scRNA-seq data are used.”——但未讨论此结构是否足以消除共线性。

  2. 估计的统计性质：MAPS的迭代解析公式给出的是矩估计还是近似MLE？是否一致？当spot内细胞数少时（每个spot实际仅捕获几个细胞），比例估计可能高度有偏。本文仅通过模拟展示，没有渐近理论（比如当spot数或基因数趋于无穷时的相合性与收敛率）。扎根语句：“We develop iteratively analytical formulae to estimate the cell type proportions and dispersion parameters.”——论文未说明公式统计动机，也未证明其最优性。

  3. 空间依赖的利用不足：MAPS没有显式利用spot间的空间相关性（每个spot独立建模），而CARD和SpatialDWLS考虑了空间结构。尽管MAPS在模拟中表现较好，但当数据中存在强空间结构（如组织层状结构）时，忽略空间信息可能会导致估计断断续续或信噪比降低。扎根语句：abstract末尾“MAPS provides a novel way for mapping spatial tissue architecture.”——但方法本身未包含空间建模。

  4. scRNA方差参数的迁移性：假设scRNA-seq方差参数 \(\phi_{gk}\) 能直接用于ST数据，但两种平台的技术噪声不同。MAPS对 \(\phi_{gk}\) 估计误差的稳健性如何？若 \(\phi_{gk}\) 估计偏大或偏小，是否会导致比例估计的系统偏差？扎根语句：“One novelty ... is that the variance parameters of scRNA-seq are used to describe the within-cell-type variations or stochasticity.”——该假设未在真实数据中被检验（例如，使用同一组织切片的scRNA和ST双平台数据比较 \(\phi\) 是否相似）。

注意：以上问答严格按照prompt结构，但缺少论文完整intro，部分内容基于抽象推断。如需更精确的引用和细节，必须阅读全文。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论