Structure learning for zero-inflated counts with an application to single-cell RNA sequencing data¶

作者: Thi Kim Hue Nguyen, Koen van den Berge, Monica Chiogna, Davide Risso
来源: Annals of Applied Statistics
主题: 其他
相关性: 6/10
链接: 期刊页 · arXiv

一、领域脉络与小综述¶

这个方向是什么¶

本论文研究的子方向是从高维稀疏计数数据中推断变量间的条件独立结构（即无向图结构学习）。在单细胞RNA测序（scRNA-seq）应用中，观测数据是一个n×p的计数矩阵（n个细胞，p个基因），每个元素记录单个细胞中某个基因的转录本数目（UMI计数），具有“高维、高方差、大量零值”的特点。问题核心是：在非高斯、零膨胀的数据特征下，如何恢复基因之间的局部相互作用网络（无向图），使得边代表给定其他基因后的条件依赖关系。目前该领域的成熟度属于“方法已有大量积累，但尚未就scRNA-seq的最佳建模策略达成一致”。

发展脉络¶

以下基于论文引用句与被引摘要，梳理成时间线与关键节点：

奠基工作（2000年代初）：图模型的结构学习最初面向高斯连续数据或分类数据。Spirtes等（2000）提出PC算法（constraint-based结构学习），Csiszár & Talata（2004）给出马尔可夫随机场邻域一致估计。Drton & Maathuis（2016）的综述系统总结了基于惩罚似然（graphical lasso）和邻域选择（neighborhood selection）的框架。这些方法假设节点条件分布来自多元高斯或有限状态空间，不适用于高方差计数数据。
向非高斯模型扩展（2012–2015）：Yang et al.（2013）提出一类“节点条件分布取自单变量指数族”的图模型（如泊松、负二项），用节点wise的M估计加ℓ1惩罚联合恢复图结构。该工作为计数数据图模型提供了理论支撑（sparsistency），但泊松不能处理过离散，负二项不能处理零膨胀。几乎同时，Chandrasekaran et al.（2010）引入潜变量高斯图模型，可间接处理观测变量中的异常相关性，但仍是高斯框架。
单细胞数据特化与零膨胀争议（2017–2021）：单细胞技术导致数据中出现大量结构零（dropout）。Risso et al.（2017）提出ZINB-WaVE，用零膨胀负二项模型为scRNA-seq做降维与归一化，获广泛使用。但Townes et al.（2019）、Sarkar & Stephens（2020）基于UMI计数数据论证“不存在真正的零膨胀，负二项就够了”，引发激烈争论。本文作者在引用句中明确提及这一争议（“the need for modeling zero-inflation … is at the core of an ongoing debate”），并选择站在ZINB一边。
当前frontier与本文位置：截至论文发表，面向scRNA-seq的图结构学习贫乏直接方法。一方面，通用指数族图模型（Yang et al. 2013，Gallopin et al. 2013的Poisson log-normal）可以处理计数，但未被针对零膨胀设计；另一方面，零膨胀模型已有降维应用（ZINB-WaVE），但未经改造用于图学习。本文的定位是：将ZINB分布嵌入节点条件模型，以此构建专为scRNA-seq设计的图结构学习方法。

子线索聚类¶

线索A：基于节点wise指数族回归的图学习（neighborhood selection）：Yang et al.（2013）为理论代表；PC-stable算法（Colombo & Maathuis 2012）的邻域搜索也可归入此路。本文采用类似思路：对每个节点用ZINB回归拟合邻居，再对称化得到图。
线索B：基于惩罚似然的联合估计：经典graphical lasso（高斯）或泊松图（如Gallopin 2013的Hierarchical Poisson log-normal），试图一次性优化整个图。本文似乎未采用这一路线（从Abstract看是“learning the structure of a graph”，但方法细节未给出）。
线索C：零膨胀建模与降维：Risso et al.（2017）的ZINB-WaVE将ZINB用作因子模型；本文将其改为图结构学习任务中的节点条件分布。
线索D：模型选择与稳定性：Liu et al.（2010）的StARS被用于选择正则化参数，本文也使用了StARS（见“In fact, when the model is correctly specified, the F1-scores”的引用语境）。

这个方向在追问的核心问题¶

零膨胀是否是scRNA-seq图学习必须建模的特征？
若ZINB比NB更准，则假阴性边（因dropout被误判为独立）可以减少。若ZINB过度拟合（参数多），可能引入虚假边。
高维p≫n下，ZINB图学习的估计一致性需要什么条件？
节点wise回归要求每个节点的回归稀疏；参数估计的不稳定性（尤其是零膨胀概率π）是否影响图选择？
如何选择合适的正则化强度才能平衡边检测的灵敏度和特异度？
StARS提供了保守选择，但本文未给出关于error rate的理论保证。
图的因果方向能否恢复？
本文是无向图（条件独立线图），不提供方向；但许多下游生物学分析需要方向（regulatory network）。

⚠️ 作者的framing¶

“作者的缺口”：现有图学习方法专为连续或简单计数变量设计，不能恰当处理scRNA-seq的零膨胀与高方差；本文提供专为ZINB设计的结构学习框架是“显然的下一步”。
被淡化/回避的竞争路线：
- Townes（2019）、Sarkar & Stephens（2020）的“零膨胀是多余”论据，几乎未被正面回应（仅在引用句中提及争议，未见模拟对比负二项与ZINB的图恢复效果）。
- 单纯的负二项图模型（如Yang 2013中的NB特例）作为baseline，但作者未在Abstract中报告其相对表现。
- 降维后图学习（先PCA再高斯图）与现成工具（如WGCNA）被回避：它们不基于概率模型，但简单。

值得查的遗漏：没有引用和讨论基于copula的图模型（如非参秩相关方法）或稀疏广义线性模型中的edge selection方法（如Meinshausen & Bühlmann 2006的高斯邻域选择在非高斯下的扩展）。另外，作者本可以引用Sarkar & Stephens（2021）的完整论文（只有预印本？实际上被引条目17是2021的？），但似乎并未深入讨论其技术含义。

张力¶

被引工作中存在明显张力：Risso et al.（2017）坚定使用ZINB vs. Townes et al.（2019）、Sarkar & Stephens（2020）认为零膨胀不是真正的零膨胀。该张力在引用句中被直接点出。但本文选择ZINB后，并未安排针对不同的零假设（NB vs ZINB）进行比较验证（根据Abstract，模拟只说明了在“多种设定”下能恢复，但未说是何种设定）。这一张力是未来研究者值得深挖的切入点。

二、最核心、最简单的例子/数学问题（先把符号/模型/可观测数据交代清楚）¶

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

n: 细胞个数（样本量）
p: 基因个数（变量维数）
Yₙₓₚ: n×p的观测矩阵，条目yᵢⱼ表示第i个细胞中第j个基因的UMI计数（非负整数）
G = (V, E): 目标无向图，V = {1,…,p}，E⊆V×V表示条件依赖关系。邻接矩阵A，Aⱼₖ=1当且仅当(j,k)有边。
节点j的邻居集: ne(j) = {k∈V: (j,k)∈E}

模型（ZINB分布，每个节点j独立给定其邻居）：

对每个节点j，给定其邻居的计数yᵢₙₑ₍ⱼ₎，yᵢⱼ的条件分布为ZINB(μⱼ,θⱼ,πⱼ)，其中：
- μⱼ = g⁻¹(ηⱼ)（均值），ηⱼ = βⱼ₀ + Σₖ₋ₙₑ₍ⱼ₎ βⱼₖ× (某种链接函数下的邻居值)
- θⱼ: 离散参数（负二项的形状），πⱼ: 零膨胀概率。
实际参数化：
ZINB的pmf:
P(Y=0) = πⱼ + (1-πⱼ) (θⱼ/(θⱼ+μⱼ))^θⱼ
P(Y=y) = (1-πⱼ) × Γ(y+θⱼ)/(Γ(θⱼ) y!) (θⱼ/(θⱼ+μⱼ))^θⱼ (μⱼ/(θⱼ+μⱼ))^y, y>0.

可观测数据：Y（n×p计数矩阵）。研究者能观测到的只有yᵢⱼ。
想要但观测不到的：图G（邻接矩阵）、每个节点的回归系数βⱼ、离散参数θⱼ、零膨胀概率πⱼ、以及所有潜在（latent）的细胞状态或批次效应（虽然ZINB-WaVE包含潜变量，但本文的图模型是否纳入？未明说）。

核心识别假设：图G是马尔可夫网——给定邻居集，节点j与其他非邻居节点条件独立；同时假设分布是ZINB节点条件同构的（类似Yang 2013的“节点wise条件模型确定联合分布”），需满足某种一致性条件（pseudolikelihood或compatibility）。

第二步：最小内核¶

最简特例：假设p=2，即只有两个基因。目标：判断它们之间是否有边（即是否条件独立）。由于没有其他节点，“条件独立”退化为“独立”。此时：

观测到n对计数 (aᵢ,bᵢ) ∈ ℕ²。若边存在，则联合分布必须是可表示成在某个隐变量（可能无）下的某种依赖——但用ZINB节点条件无法直接定义二元联合。最简单的方法：拟合两个模型：
- 模型0（独立）：Y₁~ZINB, Y₂~ZINB，相互独立。
- 模型1（相依）：用节点wise回归思想：Y₁ = ZINB₁(μ₁=exp(β₁₀+β₁₂b))，Y₂ = ZINB₂(μ₂=exp(β₂₀+β₂₁a))，两式联合最大化某种损失。

对于p>2的一般情况，此处的最小内核其实是邻域选择（Meinshausen & Bühlmann 2006）在ZINB节点回归下的应用：
- 对每个节点j，以yⱼ为响应，其余p-1个变量为候选预测变量，拟合带ℓ₁惩罚的ZINB回归。被选中的预测变量构成ne(j)的候选。
- 对称化（AND/OR规则）：若j在k的邻居集中且k在j的邻居集中，则确定边(j,k)。

因此，本文的核心技术是将节点wise ZINB回归 + ℓ₁惩罚 + 对称化应用于图学习——这是Yang et al.（2013）框架的特例，只是把指数族从NB换成ZINB。

最小问题：给定n个独立同分布观测(yᵢ₁,…,yᵢₚ)，假设真实图G稀疏且回归系数稀疏，能否用ℓ₁正则化ZINB节点回归一致地恢复边？关键困难：ZINB的似然非凸（因零膨胀概率π与θ、β耦合），且ℓ₁正则化只能处理凸损失（NB似然是凸的，加上π后不是）。本文如何绕开？很可能采用矩估计或分步估计（如先估计π，再估计β,θ），或者采用惩罚似然的近似（如基于负对数似然的凸上界）。论文摘要提到“penalized likelihood or moment estimation methods”，但具体未明。

三、这篇论文做了什么（本次重心，务必讲透）¶

由于只看到Abstract和被引论文摘要，无法精确复述每项技术细节，以下用合理推断并标注“推断”部分；若有实际全文，应以真实内容为准。

三句话¶

研究了什么问题：如何利用零膨胀负二项分布，从单细胞RNA-seq计数数据中学习无向条件独立图。
核心工具/方法：以ZINB为每个节点的条件分布，采用节点wise的惩罚似然或矩估计（如method-of-moments）与ℓ₁正则化，结合StARS选择正则化参数，并使用PC-stable算法的对称化变体得到最终图。
主要结论：在模拟设定下，该方法（记为PC-zinb1及其变体PC-zinb0, PC-nb, PC-pois）能有效恢复图结构（F1-score高），并在真实嗅觉干细胞分化和p53数据集上得到生物可解释的模块。

关键设定与假设¶

SUTVA？ 不适用，非因果。
零膨胀假设：每个基因的表达在给定邻居基因的表达后，观测到额外零的概率由独立参数πⱼ控制。该假设是模型核心，也是争议焦点。
稀疏性：真实图边数远小于p²；与大多数高维图学习相同。
节点条件分布的一致性（compatibility）：所有ZINB节点条件构成一个有效联合分布——这需要验证。不同于高斯、Ising，ZINB节点条件不一定能嵌入多元分布（有潜在可加性冲突）。方法实际是在估计一个伪似然或结点wise回归，并不保证存在一个联合分布的边际条件与模型一致。作者可能使用伪似然框架（Besag 1975），并在渐近论证中引用Yang et al.（2013）中类似处理。
正则化选择：采用StARS而非CV，以获得更稳定的稀疏解。
PC-stable算法的使用：与Colombo & Maathuis（2012）相同，用条件独立性检验来导图，但检验统计量基于ZINB模型的偏差或似然比——这要求条件独立性检验是可计算的（如通过拟合ZINB回归并检验回归系数为零）。作者将节点wise回归的邻域选择作为条件独立性检验的替代（类似“neighborhood selection corresponds to conditional independence testing in exponential families”的论点，见Yang 2013）。

主要结果¶

模拟结果（来自引用语境：“the PC-zinb1 algorithm and its variants … are consistent in terms of reconstructing the structure from given data. In fact, when the model is correctly specified, the F1-scores”）：
- PC-zinb1（用ZINB + ℓ₁ + StARS）在模型正确时F1高；
- PC-zinb0（π=0退化为NB）、PC-nb、PC-pois作为对照表现递减；
- 高维场景下仍能恢复。
真实数据：
- 嗅觉干细胞分化（Fletcher et al. 2017）：分析242个“stem cell differentiation”基因（注释自MGD），在激活的HBC细胞类型中的表达，构建图后识别出已知调控模块（如Sox9, Ptn等）。
- p53网络：Marin Navarro et al.（2020）的数据，展示p53与p63靶点基因网络的结构。

证明路线与技术技巧¶

由于缺少全文，以下基于此类论文的标准路线推断（标记为[推断]）：

整体路线
对每个节点j，用带ℓ₁惩罚的ZINB极大似然（或惩罚矩估计）拟合与其余所有变量Y₋ⱼ的回归。
通过交叉验证或StARS确定λ，得到节点j的邻居候选je。
对所有j重复，再应用对称规则（AND/OR）得到最终图。
用StARS稳定选择（多次subsampling）筛选边。
关键跳跃点
ZINB似然的非凸性：若同时优化π,θ,β，似然非凸。作者可能采用交替优化：固定π，对β,θ做凸优化轮换；或者用矩估计（如基于经验零率与NB比率）先估计π，再套用NB惩罚回归。
一致性条件：需要证明节点wise回归恢复的邻域以高概率等于真实图，即sparsistency。这通常需加条件如最小系数β_min足够大（beta-min条件）、维数log(p)/n→0、设计矩阵的受限特征值（Restricted Eigenvalue）条件——但ZINB的Fisher信息矩阵的谱将与μ与θ相关。
模型错误指定：若真实分布是NB、ZINB是过度拟合，差异如何影响图选择（专门讨论有限）。
技术技巧
节点wise pseudo-likelihood（避免定义多元分布，只对各节点条件做MLE+lasso）
StARS（Liu et al. 2010）用于正则化选择
PC-stable算法（Colombo & Maathuis 2012）的条件独立性检验替代为lasso选择的邻域（实际该算法本来的独立性检验是离散/高斯z检验；此处用ZINB的似然比检验或系数检验，但文中提及PC-stable时称“to estimate the neighborhood of each node, we employ the PC-stable algorithm”，可能意味着以条件独立检验的p值为依据逐步删除边，而非lasso——矛盾待核实）
[推断] 可能使用赤池信息或BIC变体选择模型进行小规模的模型比较。

真实例子与应用¶

嗅觉干细胞（Fletcher et al. 2017）：基础基因集242个，检测激活HBC细胞中的表达。结果识别的基因模块包含已知干性因子（Ptn, Sox9），方法重建的模块支持“stem cell lineage trajectories”的文献结论（但论文未声称发现新边，而是验证已知调控关系）。
p53干细胞（Marin Navarro et al. 2020）：分析人类神经干细胞，展示p53缺失扰动后的基因网络，发现OXPHOS下调相关边连通性变化。意义：方法能复现已知p53靶基因互作，同时发现意外连接。
这些例子旨在说明：方法的输出具有生物学可解释性，并不显示盲选的结果。

🔎 结论是否比证明窄¶

[推测] 文章在摘要中声称“retrieve the structure of a graph in a variety of settings”，但未提供关于图估计的渐近一致收敛速度或error control的理论。只通过模拟展示了F1分数——它小于100%，说明在稀疏或噪声大的设置下恢复并不完美。此外，没有给出“零膨胀不必要”场景下的基准表现；结论的普遍性被限制在ZINB模型是正确的设定下，但现实scRNA-seq数据模型的正确性正是辩论焦点。论文可能在讨论中承认此点，应核实。

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

零膨胀假设的必要性检验
请读Townes et al.（2019）和Sarkar & Stephens（2020）的完整论文，用本文方法在确实只有NB（非零膨胀）的UMI数据上比较ZINB vs NB图恢复的F1差；若ZINB并不更好，则本文所用模型的识别优势源于过拟合。这篇论文的模拟仅针对模型正确设定（自然偏向ZINB）——这是引言中“ongoing debate”的落脚点。
图估计的一致性理论
本文未在文中给出渐近一致性的数学证明。可追问：在p随n增长、真实图稀疏且ZINB节点条件正确设定下，ℓ₁惩罚的ZINB回归是否达成sparsistency？需要什么条件（如beta-min, 受限特征值在东不同θ,π下的界）？这与Yang et al.（2013）的结果直接相接，但ZINB的非凸性使得该问题不平凡。
伪似然与联合分布的存在性
文中是否确认了节点wise ZINB条件对应的多元分布存在且唯一？在没有可加性结构时，伪似然框架可能不保证图参数与真实图的一一对应。可参考Besag（1974）和Drton（2016）的权威讨论。
扩展至因果图（有向无环图）
本文的无向图不能直接用于解析调控方向。若将ZINB节点分布推广至有向模型（类似贝叶斯网，但有向边不能出现在循环中），评分函数的可分解性仍是开放问题。这依赖于ZINB在因果推断中作为conditional分布的使用——这与您（陈星宇）的因果推断兴趣恰好吻合。
计算效率与可扩展性
单细胞数据可达数万细胞、数万基因；节点wise ZINB回归需要估计每个节点的p个系数(β)和1个θ、1个π，复杂度高。可否用二阶矩估计（如截面矩条件）代替MLE，实现比lasso更快的求解？有没有利用您熟悉的高阶U统计计算（如树分解）来加速邻域搜索的可能？——这直接衔接您武器库中的“high-order U-statistics computation via tensor contraction”。

Maintained by 陈星宇 · Homepage · Source on GitHub