On the Recoverability of Causal Relations from Bulk Gene Expression Data¶
会议: NeurIPS 2026
arXiv: 2606.00568
领域: 其他
关键词: Causal Discovery, Bulk Gene Expression, Structural Causal Model, Identifiability, Transcriptomics
一句话总结¶
系统形式化了单细胞混合聚合对因果发现的破坏机制，证明了函数形式一致性与条件独立一致性在非退化线性聚合下可恢复的充要条件必须要求结构方程为仿射变换，并在真实转录组数据中实证了线性假设的不成立。
研究背景与动机¶
从基因表达数据中推断基因之间的因果调控网络（Gene Regulatory Networks, GRN）是计算生物学与系统生物学的核心挑战之一。尽管单细胞转录组测序技术（scRNA-seq）在近年来得到了长足发展，群体混合转录组（Bulk Gene Expression Profiling，包括 Microarray 和 Bulk RNA-seq）因其具有测序通量高、批次与测量噪声低、检测灵敏度高以及建库成本低廉等优势，依然在大规模临床队列和药物分子筛选中扮演着不可替代的角色。受此驱动，过去数十年间学者们提出了大量基于信息论模型、贝叶斯网络和深度神经网络的计算发现算法，尝试直接利用组织水平的 Bulk 表达矩阵重构基因间的线性和非线性因果关系。
然而，现有方法普遍存在一个被长期忽略的根本性前提缺陷：算法几乎无一例外地将 Bulk 表达量视为单个同质因果系统的直接观测变量，默认宏观混合测量能够保留微观细胞层面的因果语义。事实上，Bulk 表达测量是将生物样本中成千上万个细胞的 RNA 物理混合后进行检测的有损（lossy）、不可逆（non-invertible）粗粒化聚合（coarsening）过程。这种跨细胞尺度的聚合操作会严重扭曲统计相关性、函数依赖结构以及条件独立性（Conditional Independence, CI）。例如，在微观真实调控链 \(X \to Y \to Z\) 中本来满足条件独立 \(X \perp\!\!\!\perp Z \mid Y\)，但在宏观观测其聚合总量 \(\bar{X}, \bar{Y}, \bar{Z}\) 时，由于以 \(\bar{Y}\) 为条件无法完全消除底层潜变量 \(Y\) 的残留变异协调效应，往往导致条件独立性破坏（\(\bar{X} \not\perp\!\!\!\perp \bar{Z} \mid \bar{Y}\)）。由此引发了一个悬而未决的基础理论问题：从 Bulk 表达数据中推断出的因果关系究竟在何种条件下能够真实还原底层的基因调控网络？
针对这一关键空白，本文在因果发现与聚合因果理论框架下展开严格的形式化推导，通过函数形式一致性（Functional-form Consistency）与条件独立一致性（Conditional-independence Consistency）给出了因果可恢复性的充要条件。核心 idea：形式化单细胞聚合过程为有损投影，从数学上严格证明了在常见的线性混合聚合下，只有当底层基因结构方程完全退化为仿射关系（Affine structural equations）时，因果关系与条件独立结构才能在聚合数据中被忠实保留；进一步通过真实转录组统计诊断证明基因间普遍呈强非线性调控，从而警示在缺乏强辅助先验时不可直接从 Bulk 表达数据恢复因果图。
方法详解¶
整体框架¶
本文构建了一套用于评估聚合系统因果可恢复性的形式化理论框架，其核心逻辑分为三步：首先严格刻画微阵列（Microarray）分位数归一化与 Bulk RNA-seq（DESeq2）标准归一化流水线，确立宏观测量本质上是底层单细胞转录本的加权线性求和；随后建立两套正交的一致性准则（函数形式一致性与条件独立一致性），推导可恢复性的充要代数条件；最后设计两类统计度量用于量化真实转录组数据与理论充要条件的契合程度。
%%{init: {'theme': 'base', 'themeVariables': {'primaryColor': '#EEF2FF', 'edgeLabelBackground':'#FFFFFF', 'tertiaryColor': '#F9FAFB'}}}%%
flowchart TD
subgraph MicroLevel ["微观潜变量层（Single-Cell Micro-Level）"]
SC["单细胞表达状态 V ~ P_V<br/>结构因果模型 SCM: Y = f(X, ε_Y)"]
CI_true["微观条件独立性<br/>Chain/Fork: X ⟂⟂ Z | Y<br/>Collider: X !⟂⟂ Z | Y"]
end
subgraph Aggregation ["跨细胞物理聚合与归一化（Lossy Coarsening）"]
Sum["细胞群物理池化（RNA Pooling）<br/>x_bar = g({x_j}_j∈K) = Σ w_j x_j"]
Norm["标准生物信息归一化<br/>Microarray: Quantile f_M<br/>RNA-seq: DESeq2 Size Factor f_R"]
end
subgraph MacroLevel ["宏观观测层（Bulk Macro-Level）"]
MacroObs["Bulk 观测变量: X_bar, Y_bar, Z_bar"]
MacroCI["宏观关联与条件依赖检验<br/>X_bar ⟂⟂? Z_bar | Y_bar"]
end
subgraph ConsistencyTheorems ["理论充要条件证明（Theoretical Guarantees）"]
T1["定理 4.3 & 推论 4.4/4.5<br/>函数形式一致性: g(F(x_K)) = h(g(x_K))<br/>ANM / PNL 保持可恢复 ⇔ 结构方程 f 为仿射映射"]
T2["定理 4.6 & 推论 4.7<br/>条件独立一致性: Cov_Y|Y_bar(P(X_bar|Y), P(Z_bar|Y)) = 0<br/>Fork/Chain CI 保持 ⇔ 线性潜变量系统"]
end
subgraph RealValidation ["转录组非线性实证诊断（Empirical Diagnostics）"]
GAM["广义加性模型 GAM: y = μ + s(x) + η<br/>反向有效自由度: L^edf = 1 / max{edf, 1}"]
AffDev["最优仿射偏差归一化评分: D^aff<br/>评估单细胞与 Bulk 基因对偏离仿射的程度"]
end
SC --> Sum
CI_true --> Sum
Sum --> Norm
Norm --> MacroObs
MacroObs --> MacroCI
SC -.推演约束.-> ConsistencyTheorems
MacroCI -.推演约束.-> ConsistencyTheorems
ConsistencyTheorems --> RealValidation关键设计¶
1. 归一化保留聚合本性（Aggregation Invariance under Normalization）
在推导因果理论之前，必须明确宏观观测数据的物理生成机制。对于 Microarray 技术，荧光强度 \(D_{ij}^{mi} = S_{ij} + N_{ij}\) 经分位数归一化 \(f_M\) 与探针集汇总后，基因表达估计值可形式化为： $\(\tilde{\theta}_j^k = f_M\left(\sum_{i \in P_K} D_{ij}^{mi}\right) = f_M\left(\sum_{i \in P_K} \sum_{c=1}^n \theta_{kc}^{ij}\right)\)$ 对于 Bulk RNA-seq，测序读段数 \(D_{kj}^{se} \sim \mathrm{NB}(\mu_{kj}, \gamma_k)\)，经 DESeq2 大小因子 \(s_j = \mathrm{median}_k (D_{kj}^{se} / G_k)\) 缩放与广义线性模型拟合后，归一化表达值为 \(\tilde{o}_j^k = f_M\left(\sum_{c=1}^n o_{kc}^j\right)\)。作者明确指出，常规归一化旨在校正文库大小和批次效应等技术变异，但完全无法逆转底层成千上万个细胞转录本的物理求和操作。宏观变量本质上可统一建模为微观潜变量集上的加权线性映射 \(g(x_K) = \sum_{j=1}^m w_j x_{k_j}\)。
2. 函数形式一致性与柯西函数方程约束（Functional-Form Consistency via Cauchy Equation）
函数形式一致性要求聚合后的宏观变量仍然遵从与微观相同的结构因果模型（SCM）函数族，从而使得基于函数不对称性的因果发现算法（如 ANM 和 PNL）成立。设微观加性噪声模型为 \(Y = f(X) + \varepsilon_Y\)，聚合变量满足 \(\bar{Y} = g(y_K) = g(F(x_K)) + g(\varepsilon_{y_K})\)。由于外生噪声独立于原因变量，聚合后的模型仍为 ANM 的充要条件是存在单变量函数 \(h\) 使得： $\(g(F(x_K)) = h(g(x_K)) \implies \sum_{j=1}^m w_j f(x_{k_j}) = h\left(\sum_{j=1}^m w_j x_{k_j}\right)\)$ 对于非退化线性聚合（存在至少两个非零权重），令非主导坐标置零，并定义偏置消除变换 \(\phi(x) = f(x) - f(0)\) 与 \(H(t) = h(t) - (\sum w_j)f(0)\)，该条件退化为经典的柯西函数方程： $\(H(s + t) = H(s) + H(t), \quad \forall s, t \in \mathbb{R}\)$ 在可测性条件下，柯西方程的唯一解为齐次线性函数 \(H(t) = ct\)。由此严格导出推论 4.4：微观函数 \(f(X)\) 必须严格为仿射函数 \(f(X) = \alpha X + \beta\)。对于后非线性模型（PNL）\(Y = f_2(f_1(X) + \varepsilon_Y)\)，两步推导同样强力限制外层可逆变换 \(f_2\) 与内层函数 \(f_1\) 必须同时为仿射形式。
3. 条件独立性保持与残余协同变异定理（CI Preservation & Zero Residual Co-Variation）
在多变量图模型（如 PC 算法）中，因果骨架依赖于条件独立测试。论文深入剖析了宏观聚合对马尔可夫独立性的冲击。在对撞结构（Collider \(X \to Y \leftarrow Z\)）下，以宏观后代 \(\bar{Y}\) 为条件自然打开路径，保留 \(\bar{X} \not\perp\!\!\!\perp \bar{Z} \mid \bar{Y}\)；然而在链式（Chain \(X \to Y \to Z\)）与叉式（Fork \(X \leftarrow Y \rightarrow Z\)）结构中，微观上 \(X \perp\!\!\!\perp Z \mid Y\)，但宏观以 \(\bar{Y}\) 为条件并不等同于固定了微观状态向量 \(\mathbf{Y}\)。 作者通过全概率公式边缘化微观状态，推导出了宏观条件独立的充要条件（定理 4.6）： $\(\operatorname{Cov}_{\mathbf{Y} \mid \bar{Y}}\left(P(\bar{X} = \bar{x} \mid \mathbf{Y}), P(\bar{Z} = \bar{z} \mid \mathbf{Y})\right) = 0\)$ 该协方差项度量了微观未知潜状态在被宏观均值约束下的“残余协同变异”。只有当不同微观组合不会系统性地同时拉高或拉低 \(\bar{X}\) 与 \(\bar{Z}\) 的条件概率时，独立性才得以维持。在推论 4.7 中，作者证明了当潜变量完全由线性系统生成且噪声独立时，残余变异被退化隔离，此时宏观条件独立性严格保持；但一旦存在非线性关系，协方差普遍非零，必然引入伪因果关联或漏检真实调控边。
损失函数 / 训练策略¶
作为因果可恢复性理论推导与诊断工作，本文不涉及传统意义上的神经网络反向传播损失函数，而是构建了严格的非线性诊断度量与统计检验流水线：
1. 广义加性模型反向有效自由度评分（Inverse-edf Linearity Score）
为了经验性检验基因调控函数是否满足定理要求的仿射性，对标准化表达矩阵中的任意有序基因对 \((X, Y)\)，拟合单变量样条广义加性模型（GAM）： $\(y_i = \mu + s(x_i) + \eta_i\)$ 其中 \(s(\cdot)\) 为三次平滑样条函数。提取拟合平滑项的有效自由度（Effective Degrees of Freedom, \(\widehat{\mathrm{edf}}_{j \to l}\)）。\(\mathrm{edf}\) 越接近 1 说明关系越趋向线性，自由度越大说明非线性曲率越高。定义反向自由度指标： $\(L_{j \to l}^{\mathrm{edf}} = \frac{1}{\max\{\widehat{\mathrm{edf}}_{j \to l}, 1\}} \in (0, 1]\)$
2. 归一化最优仿射偏差评分（Normalized Affine-Deviation Score）
直接从几何距离角度量化非线性残差。令 \(\hat{f}_i = \hat{\mu} + \hat{s}(x_i)\) 为 GAM 的预测值，通过普通最小二乘法求解最佳仿射逼近参数： $\((\hat{\alpha}, \hat{\beta}) = \arg\min_{\alpha, \beta} \sum_{i=1}^n \left[\hat{f}_i - (\alpha x_i + \beta)\right]^2\)$ 计算预测函数偏离最优仿射超平面的归一化标准差： $\(D_{j \to l}^{\mathrm{aff}} = \left( \frac{\sum_{i=1}^n \left[\hat{f}_i - (\hat{\alpha} x_i + \hat{\beta})\right]^2}{\sum_{i=1}^n \left[\hat{f}_i - \bar{f}\right]^2} \right)^{1/2}, \quad \bar{f} = \frac{1}{n}\sum_{i=1}^n \hat{f}_i\)$ 若调控函数严格为仿射形式，则 \(D_{j \to l}^{\mathrm{aff}} = 0\)；偏离程度越大，该指标越接近 1。
实验评估¶
数据集与实验设置¶
- 可控合成数据集（Synthetic Benchmarks）：构建微观结构因果模型，涵盖加性噪声（ANM）与后非线性模型（PNL），因果函数包含线性（\(ax+b\)）以及由 \(\{\exp(x), x^p, \mathrm{ReLU}(wx+b), \sin(x)\}\) 组成的强非线性集合；后处理变换 \(f_2 \in \{\exp(x), x^3, \sinh(x)\}\)。聚合机制分别对比线性求和/均值（Mean）与非线性次序聚合（Median）。针对 Fork 与 Collider 模体进行 1,000 次独立重复实验，条件独立性判定阈值取 \(p = 0.05\)。
- 真实转录组数据集（Real-world Transcriptomics）：
- 4 个大规模 Bulk 表达数据集：微阵列结肠癌数据集 GSE39582（585 样本）、乳腺癌数据集 GSE20142（1,240 样本）；Bulk RNA-seq 数据集 GSE57945（260 样本）、GSE245006（518 样本）。
- 4 个高通量单细胞数据集（Perturb-seq / CRISPRi）：Adamson、Norman、Srivatsan、Tian 单细胞文库。
- 数据预处理严格保留四分位距（IQR）排名前 5,000 的高变基因进行全配对检验。
实验结果与对比分析¶
1. 合成数据因果可恢复性基准评测¶
下表展示了函数型方法（ANM, PNL）与约束型方法（Fork-UI/CI, Collider-UI/CI）在聚合前后对因果方向与条件独立性判定的精确率（Precision，基于 1,000 次蒙特卡洛试验）：
| 生成因果机制 | 聚合算子类型 | ANM (原始) | ANM (聚合) | PNL (原始) | PNL (聚合) | Fork-CI (原始) | Fork-CI (聚合) | Collider-CI (原始) | Collider-CI (聚合) |
|---|---|---|---|---|---|---|---|---|---|
| 线性 (Linear) | 线性均值 (Mean) | 1.000 | 0.830 | 1.000 | 0.912 | 0.943 | 0.932 | 1.000 | 0.999 |
| 线性 (Linear) | 非线性中位数 (Median) | 1.000 | 0.613 | 1.000 | 0.350 | 0.960 | 0.256 | 1.000 | 0.930 |
| 非线性 (Nonlinear) | 线性均值 (Mean) | 0.958 | 0.450 | 0.838 | 0.554 | 0.846 | 0.402 | 0.932 | 0.508 |
| 非线性 (Nonlinear) | 非线性中位数 (Median) | 0.967 | 0.591 | 0.840 | 0.460 | 0.865 | 0.617 | 0.953 | 0.432 |
从表 1 数据可以清晰观察到： - 只有在“线性因果 + 线性聚合”的严格理想设定下，ANM、PNL 及 Fork-CI 的准确率在聚合后依然能够保持在 0.830 至 0.932 的高位水平，因果语义基本完好； - 一旦因果函数引入非线性（如第二行至第四行），即使用最简单的线性均值聚合，ANM 的可恢复精度断崖式下跌至 0.450，Fork-CI 暴跌至 0.402，Collider-CI 从 0.932 跌至 0.508； - 若聚合操作本身呈现非线性（如中位数），因果发现精度出现系统性瓦解。这完整印证了定理 4.3 与定理 4.6 的必要性判定。
2. 真实转录组数据线性度统计诊断¶
下表汇总了 4 个 Bulk 转录组与 4 个单细胞转录组在高变基因对上的线性诊断统计指标分布特征：
| 数据集类型 | 数据集名称 | 测序技术平台 | 样本量 / 细胞数 | 反向自由度中位数 \(1/\mathrm{edf}\) | 仿射偏差中位数 \(D^{\mathrm{aff}}\) | 近似线性基因对占比 (\(L^{\mathrm{edf}} > 0.8\)) |
|---|---|---|---|---|---|---|
| Bulk 表达谱 | GSE39582 | Microarray (Affymetrix) | 585 | 0.284 | 0.621 | < 6.8% |
| Bulk 表达谱 | GSE20142 | Microarray (Illumina) | 1,240 | 0.312 | 0.589 | < 7.9% |
| Bulk 表达谱 | GSE57945 | Bulk RNA-seq | 260 | 0.245 | 0.674 | < 5.2% |
| Bulk 表达谱 | GSE245006 | Bulk RNA-seq | 518 | 0.263 | 0.648 | < 5.9% |
| 单细胞表达谱 | Adamson | Perturb-seq (10x) | 87,000+ | 0.118 | 0.842 | < 1.4% |
| 单细胞表达谱 | Norman | Perturb-seq (10x) | 110,000+ | 0.105 | 0.867 | < 1.1% |
| 单细胞表达谱 | Srivatsan | sci-Plex | 45,000+ | 0.132 | 0.825 | < 1.7% |
| 单细胞表达谱 | Tian | CRISPRi CROP-seq | 32,000+ | 0.126 | 0.839 | < 1.5% |
诊断结果表明： - 在单细胞真实分辨率下，绝大多数基因调控对的 \(1/\mathrm{edf}\) 聚集在 0.1 左右，\(D^{\mathrm{aff}} > 0.82\)，表明真实的细胞内生化分子调控具有极强的饱和、门控和协同非线性效应； - 经过 Bulk 混合物理聚合后，根据多元中心极限定理，混合变量分布向高斯分布收敛，使得 Bulk 数据的平滑度略微提高（\(1/\mathrm{edf}\) 提升至 0.25~0.31），但这仅仅是统计平均带来的“表观线性”，实际偏离仿射的程度依然极其严重（近线性基因对比例不足 8%）； - 实验结论无可辩驳地表明，真实基因表达谱根本无法满足可恢复性定理所要求的仿射假设。
总结与思考¶
本文是一篇针对计算生物学因果发现研究具有重要“清醒剂”意义的奠基性理论工作。长期以来，生物信息学界大量借鉴计算机视觉与纯统计领域的因果网络学习算法，直接套用在组织样本 Bulk 表达矩阵上，并根据算法输出的 DAG 图做出生物调控假说。本文明确揭示了这种做法在数学基础上的内在脆弱性：微观异质细胞在物理混合聚合过程中，必然产生不可逆的因果信息丢失与虚假条件关联。
从理论深度来看，论文运用柯西函数方程精巧证明了保持 ANM/PNL 函数形式一致性时对仿射关系的绝对依赖，并通过残余条件协方差刻画了宏观条件独立的破坏机制，论证逻辑严丝合缝。从实用角度看，尽管文章指出了纯 Bulk 表达量因果图推断的理论极限，但为未来提供了清晰的研究指引：要实现可靠的因果发现，计算生物学家必须引入外生干预先验（如基因敲除与化学扰动数据）、时序动态测量（Time-series RNA-seq）、单细胞参考图谱解卷积（Deconvolution），或者将研究目标从还原微观细胞因果网络调整为推断宏观组织级表型因果模块。