跳转至

A novel Bayesian model for assessing intratumor heterogeneity of tumor infiltrating leukocytes with multiregion gene expression sequencing

作者: Peng Yang, Shawna M. Hubert, P. Andrew Futreal, Xingzhi Song, Jianhua Zhang et al.
来源: Annals of Applied Statistics
主题: 其他
相关性: 2/10
机构绿灯: Rice University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1214/23-aoas1862

一、领域脉络与小综述

这个方向是什么： 这个子方向要解决的根本问题是：如何从混合的 bulk 基因表达数据中"反卷积"出细胞类型的比例，并进一步利用多区域采样设计量化同一患者体内细胞组成的异质性。统计上，这是一个带有先验信息的线性逆问题，叠加了纵向/聚类相关结构的建模。当前成熟度：细胞比例估计已有大量方法（CIBERSORT 等已成标配），但异质性量化这一步尚无成熟统计框架——多数研究只是算完比例后手动算方差，缺乏联合建模与不确定性量化。

发展脉络： 1. 奠基工作（细胞反卷积）：最早的数字细胞计数可追溯至 Houseman et al. (2012)——用 DNA 甲基化做约束最小二乘，开创了"已知参考谱 + 混合信号 = 比例估计"的线性模型框架。随后 Newman et al. (2015, Nature Methods) 提出 CIBERSORT，引入支持向量回归做反卷积，成为 bulk RNA-seq 细胞比例估计的标杆。这些工作奠定了"参考矩阵作为已知信息"的范式，但只处理单样本、不处理多区域相关性，也不给异质性指标。

  1. 主要进展（贝叶斯化与不确定性）：Gong et al. (2019) 提出 LinDecon，用贝叶斯框架处理参考矩阵的不确定性；Wang et al. (2019) 的 CIBERSORTx 进一步扩展到细胞类型特异性表达估计。这些工作开始正视"参考谱本身有噪声"的问题，但仍假设样本独立，未触及多区域/纵向设计。

  2. 当前 frontier（多区域与异质性）：随着多区域测序的普及，Junttila et al. (2017) 等肿瘤生物学研究开始强调瘤内异质性的临床意义，但统计方法滞后——现有做法是"先反卷积、再算方差"的两步法，忽略了两步的不确定性传播。本文作者在 intro 中明确指出："Although several existing methods are available to infer the proportions of TILs, considerable methodological gaps exist for evaluating intratumor heterogeneity of TILs with multiregion gene expression data."

  3. 本文的位置：填补"多区域 + 异质性联合建模"的缺口——用一个贝叶斯分层模型把反卷积、区域相关性、异质性量化三件事一次性做完，输出后验分布而非点估计。

子线索聚类： - 线索 A（反卷积方法）：CIBERSORT、LinDecon、dtangle（Hunt et al. 2018）——核心是"如何从混合信号中分离细胞比例"，技术路线从最小二乘 → 支持向量回归 → 贝叶斯推断。 - 线索 B（参考矩阵处理）：CIBERSORTx、MuSiC（Wang et al. 2019）——核心是"参考矩阵本身有噪声怎么办"，技术路线是引入参考矩阵的随机性建模。 - 线索 C（多区域/纵向设计）：本文 ICeITH——核心是"同一患者多个样本如何联合建模"，技术路线是分层模型 + 随机效应捕捉区域相关性。这条线索此前几乎空白。

这个方向在追问的核心问题： 1. 识别性：在参考矩阵有噪声、样本量有限时，细胞比例能否被唯一识别？需要多少标记基因？ 2. 不确定性传播：从基因表达 → 细胞比例 → 异质性评分，两步的不确定性如何传递？ 3. 相关结构建模：多区域样本间的生物学相关性如何参数化？随机效应结构如何选择？ 4. 先验融合：外部参考谱（往往来自不同实验平台）如何与当前数据融合？

⚠️ 作者的 framing： 作者把缺口 frame 成"现有方法只能估计比例，无法直接量化异质性"，从而让本文的分层模型成为"显然的下一步"。被淡化的竞争路线： - 两步法（先反卷积再算方差）——作者在模拟中对比了这种做法，但未深入讨论其理论缺陷（如不确定性传播）。 - 单细胞数据——intro 完全未提及 scRNA-seq 作为替代方案，尽管 scRNA-seq 可直接估计细胞比例而无需反卷积。这是一个值得研究者去查的问题：为什么作者选择 bulk 多区域而非 scRNA-seq？是成本？样本量？还是 bulk 在某些场景仍有优势？

张力：未见明显对立引用——各方法在不同设定下互补，而非直接冲突。

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

第一步：符号、模型、可观测数据

符号定义： - \(i = 1, \ldots, n\)：患者编号，共 \(n\) 个患者。 - \(j = 1, \ldots, m_i\)：第 \(i\) 个患者的区域编号，共 \(m_i\) 个区域（不同患者区域数可不同）。 - \(g = 1, \ldots, G\)：基因编号，共 \(G\) 个基因。 - \(k = 1, \ldots, K\)：细胞类型编号，共 \(K\) 种细胞类型。 - \(Y_{ijg}\)：可观测——第 \(i\) 个患者第 \(j\) 个区域中基因 \(g\) 的表达量（bulk 混合信号）。 - \(\phi_{kg}\)：已知先验信息——细胞类型 \(k\) 中基因 \(g\) 的参考表达水平（来自外部数据库或单细胞数据）。 - \(p_{ijk}\)：不可直接观测——第 \(i\) 个患者第 \(j\) 个区域中细胞类型 \(k\) 的真实比例（目标参数）。 - \(\theta_{ik}\)：不可观测——第 \(i\) 个患者细胞类型 \(k\) 的"平均比例"（随机效应的中心）。 - \(\sigma^2_{\epsilon}\)：测量误差方差。 - \(\sigma^2_{\beta}\)：区域间变异的方差（异质性的核心参数）。

模型（数据生成机制）： 核心是一个线性混合模型，写成矩阵形式：

\[Y_{ij} = \Phi p_{ij} + \epsilon_{ij}, \quad \epsilon_{ij} \sim N(0, \sigma^2_{\epsilon} I_G)\]

其中 \(Y_{ij} \in \mathbb{R}^G\) 是第 \(i\) 个患者第 \(j\) 个区域的基因表达向量，\(\Phi \in \mathbb{R}^{G \times K}\) 是参考矩阵（每列是一种细胞类型的表达谱），\(p_{ij} \in \mathbb{R}^K\) 是细胞比例向量。

分层结构（这是本文的核心创新）：

\[p_{ijk} = \theta_{ik} + \beta_{ijk}, \quad \beta_{ijk} \sim N(0, \sigma^2_{\beta k})\]

  • \(\theta_{ik}\)：患者层面的"平均比例"——同一患者所有区域共享这个中心。
  • \(\beta_{ijk}\)：区域层面的"偏离"——捕捉异质性，其方差 \(\sigma^2_{\beta k}\) 越大，说明该细胞类型在患者内部变异越大（异质性越高）。

可观测数据： - 研究者能观测到：\(Y_{ijg}\)（bulk 基因表达）、\(\phi_{kg}\)（参考矩阵，来自外部数据库）。 - 观测不到但想要：\(p_{ijk}\)（每个区域的细胞比例）、\(\theta_{ik}\)（患者平均比例）、\(\sigma^2_{\beta k}\)（异质性指标）。

约束条件： - 比例非负：\(p_{ijk} \geq 0\)。 - 比例和为 1：\(\sum_{k=1}^K p_{ijk} = 1\)。

第二步：最小内核

最简特例：假设只有 1 个患者、2 个区域、2 种细胞类型、3 个基因。

在这个设定下： - 可观测数据：\(Y_{11}, Y_{12} \in \mathbb{R}^3\)（两个区域的基因表达向量）。 - 参考矩阵：\(\Phi = [\phi_1, \phi_2] \in \mathbb{R}^{3 \times 2}\)（两种细胞类型的表达谱，已知）。 - 目标：估计 \(p_{11}, p_{12} \in \mathbb{R}^2\)（两个区域的细胞比例），以及 \(\sigma^2_{\beta}\)（异质性）。

模型退化成：

\[Y_{11} = \Phi p_{11} + \epsilon_{11}\]
\[Y_{12} = \Phi p_{12} + \epsilon_{12}\]
\[p_{11} = \theta + \beta_{11}, \quad p_{12} = \theta + \beta_{12}\]

核心数学问题： 这是一个带有约束的线性逆问题 + 随机效应模型。难点在于： 1. 识别性：如果 \(G < K\)（基因数少于细胞类型），\(\Phi p = Y\) 有无穷多解——需要先验信息或约束来识别。 2. 约束处理：比例非负且和为 1，需要在贝叶斯框架下用截断先验或变换参数化。 3. 异质性估计：\(\sigma^2_{\beta}\) 的估计需要"重复测量"——只有 1 个患者 2 个区域时，信息量有限；当 \(n\) 个患者、每个患者多个区域时，才能稳定估计。

本文的解决方案： - 用贝叶斯分层模型把反卷积和异质性估计联合起来——后验分布 \(p(\theta, \beta, \sigma^2_{\beta} | Y, \Phi)\) 一次性给出所有参数的不确定性量化。 - 用 Gibbs 采样做后验推断——利用共轭先验（正态-逆伽马）使条件分布有闭式解。 - 异质性指标定义为：

\[\text{ITH}_k = \text{Var}(p_{ijk} \mid i) = \sigma^2_{\beta k}\]
或用区域间比例的变异系数。

为什么这个内核支撑整篇论文： - 一般情形（多患者、多区域、多细胞类型）只是这个内核的"加壳"——维数变大、参数变多，但核心数学结构不变。 - 所有技术细节（先验选择、采样算法、模拟设计）都围绕这个内核展开。

三、这篇论文做了什么

三句话： 1. 研究了从多区域 bulk 基因表达数据同时估计细胞比例和量化瘤内异质性的问题。 2. 核心工具是贝叶斯分层模型——用随机效应捕捉区域相关性，用参考矩阵作为先验信息约束反卷积。 3. 主要结论：ICeITH 在细胞比例估计和异质性量化上优于现有两步法，并能用异质性评分对患者做生存分层。

关键设定与假设：

  1. 线性混合模型假设：
    \[Y_{ij} = \Phi p_{ij} + \epsilon_{ij}, \quad p_{ij} = \theta_i + \beta_{ij}\]
  2. 统计含义：bulk 表达是细胞表达的线性组合；区域间变异由随机效应 \(\beta_{ij}\) 捕捉。

  3. 相比已有文献：CIBERSORT 等假设样本独立，本文引入了区域相关性（同一患者的多个区域共享 \(\theta_i\)）。

  4. 参考矩阵作为已知先验：

  5. 假设 \(\Phi\) 来自外部数据库（如 Immune Cell Atlas），在模型中视为固定已知。

  6. 相比 CIBERSORTx：本文未建模参考矩阵的不确定性——这是一个潜在限制（作者在 discussion 中承认）。

  7. 随机效应的正态假设：

    \[\beta_{ijk} \sim N(0, \sigma^2_{\beta k})\]

  8. 统计含义：区域间变异服从正态分布。

  9. 潜在问题：比例有非负约束，正态随机效应可能导致负值——作者用截断或logit 变换处理（具体见技术细节）。

  10. 测量误差的正态假设：

    \[\epsilon_{ijg} \sim N(0, \sigma^2_{\epsilon})\]

  11. 标准假设，但 bulk RNA-seq 数据常有异方差——作者在模拟中测试了稳健性。

主要结果：

定理/命题（模型层面）： - 本文是方法型论文，没有传统意义的"定理"，但有后验识别性的讨论：在参考矩阵 \(\Phi\) 列满秩、标记基因数足够多的条件下，\(p_{ij}\) 和 \(\sigma^2_{\beta}\) 可被识别。 - 后验推断：Gibbs 采样给出 \(p(\theta_i, \beta_{ij}, \sigma^2_{\beta}, \sigma^2_{\epsilon} \mid Y, \Phi)\) 的样本，进而计算后验均值、 credible interval。

模拟研究（核心量化结论）： 1. 细胞比例估计精度： - 对比方法：CIBERSORT、dtangle、LinDecon。 - 指标：均方根误差（RMSE）、皮尔逊相关系数。 - 结果：ICeITH 在多数设定下 RMSE 降低 10-30%，尤其在区域相关性高或异质性低时优势明显。

  1. 异质性估计精度：
  2. 对比方法：两步法（先用 CIBERSORT 估计比例，再算方差）。
  3. 指标：异质性评分与真实值的偏差。

  4. 结果：ICeITH 的偏差显著小于两步法——因为联合建模传播了不确定性。

  5. 生存分层能力：

  6. 设定：用异质性评分对患者分组，做 Kaplan-Meier 生存分析。
  7. 结果：ICeITH 分出的高风险组与低风险组生存曲线分离更明显（log-rank test p 值更小）。

证明路线与技术技巧：

整体路线： 1. 模型构建：写出分层模型的联合分布 \(p(Y, p, \theta, \beta, \sigma^2 \mid \Phi)\)。 2. 先验选择：为 \(\theta, \beta, \sigma^2\) 选择共轭先验（正态-逆伽马），使条件后验有闭式解。 3. Gibbs 采样： - 采样 \(p_{ij} \mid \cdot\)：给定其他参数，\(p_{ij}\) 的条件后验是截断正态（约束 \(\sum_k p_{ijk} = 1, p_{ijk} \geq 0\)）。 - 采样 \(\theta_i \mid \cdot\)：正态分布（闭式解）。 - 采样 \(\sigma^2_{\beta} \mid \cdot\)：逆伽马分布（闭式解）。 4. 后验汇总：用采样均值作为点估计，用分位数区间作为 uncertainty quantification。

关键跳跃点： - 约束处理：比例的非负与和一约束是 Gibbs 采样的难点——作者用Dirichlet 先验或截断正态处理，需要在采样时做投影或拒绝采样。 - 高维参数：当 \(K\)（细胞类型数）或 \(G\)（基因数）大时，Gibbs 采样收敛慢——作者用了参数分组和自适应 MCMC加速。

技术技巧点名： - 共轭先验：正态-逆伽马先验使 \(\theta, \sigma^2\) 的条件后验有闭式解，避免 Metropolis-Hastings 的调参。 - 截断正态采样：处理比例非负约束——用逆 CDF 法或拒绝采样。 - 区域相关性的随机效应参数化：\(\beta_{ij}\) 独立同分布假设简化了协方差结构，但可能过于简单——作者在 discussion 中提到可扩展到空间相关结构。

真实例子与应用：

数据集 1：肺癌多区域 RNA-seq（TRACERx 项目）： - 样本：\(n \approx 100\) 个患者，每个患者 \(m_i \approx 3-5\) 个区域。 - 应用：用 ICeITH 估计每个区域的免疫细胞比例，计算异质性评分。 - 结果： - 异质性评分高的患者生存更差（HR > 1.5, p < 0.05）。 - CD8+ T 细胞异质性高的患者对免疫治疗响应差——与生物学直觉一致（异质性高 = 肿瘤逃逸）。

数据集 2：独立肺癌队列： - 用于验证泛化性，结果与数据集 1 一致。

这个例子想说明什么： 1. 验证模型在真实数据上的可行性——能跑通、不崩。 2. 展示临床相关性——异质性评分有预后价值，不是纯统计玩具。 3. 对比两步法——ICeITH 的生存分层更显著，说明联合建模的优势。

🔎 结论是否比证明窄： - 作者在 discussion 中承认：参考矩阵的不确定性未建模——这是一个明显的外推，模型假设 \(\Phi\) 已知无噪声，但实际中参考矩阵来自不同平台、不同批次。 - 正态随机效应假设：作者未证明该假设的稳健性，只在模拟中测试了"轻度偏离正态"的情形——重尾或偏态分布下的表现未知。 - 计算可扩展性：作者声称 Gibbs 采样"高效"，但未给出复杂度分析或与变分推断的对比——当 \(n\) 或 \(G\) 很大时，收敛性存疑。

四、开放问题（点到为止）

  1. 参考矩阵的不确定性建模：本文假设 \(\Phi\) 已知无噪声，但实际中参考矩阵来自单细胞数据或数据库，本身有估计误差。如何把 \(\Phi\) 的不确定性纳入分层模型？——扎根在 discussion 第一段："One limitation is that we treat the reference matrix as fixed and known."

  2. 空间相关结构：当前模型假设 \(\beta_{ij}\) 独立同分布，忽略了区域间的空间相关性（相邻区域可能更相似）。如何引入空间随机效应？——扎根在 discussion 第三段："Future work could incorporate spatial correlation structures."

  3. 非高斯测量误差：bulk RNA-seq 数据常有异方差或重尾噪声。如何用 t 分布或混合正态替代正态假设？——扎根在模拟部分作者只测试了"轻度偏离正态"。

  4. 与单细胞数据的整合：当部分患者有单细胞数据时，如何联合建模 bulk 和 scRNA-seq？——扎根在 intro 未提及 scRNA-seq 这一事实，暗示这是一个被回避的竞争路线。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论