当前位置：首页 > news >正文

ARGEN：单细胞因果基因网络重建方法解析

news 2026/6/5 8:22:46

1. ARGEN：基于Perturb-seq数据的因果基因网络重建方法解析

在单细胞生物学领域，理解基因间的调控关系是揭示细胞功能机制的核心挑战。传统方法如基因共表达网络只能反映相关性，而因果基因网络（Causal Gene Networks, CGNs）通过有向无环图（DAGs）揭示基因间的定向调控关系，为机制研究提供更直接的证据。Perturb-seq技术的出现为因果推断带来了新机遇——它结合单细胞CRISPR扰动与转录组测序，能在数千个基因敲除/激活的背景下获取单细胞分辨率的功能响应数据。

然而，现有DAG学习方法（如PC算法、GES、NOTEARS等）存在两大局限：一是假设所有混杂因素可观测，二是要求干预完美实施。这与实际Perturb-seq数据特性严重不符：细胞周期、染色质状态等潜在混杂因素无法完全测量；实验设计上往往只干预关键基因；数据存在高维度、样本不平衡等技术噪声。这些因素导致传统方法在Perturb-seq场景下会产生严重偏差。

ARGEN（Arbitrary-confounder Robust causal GEne Network）的创新在于：

工具变量策略：利用gRNA结合信号作为工具变量，通过代理基因表达抵消未观测混杂的影响
扩展的结构方程模型：在经典scRNA-seq测量模型基础上嵌入扰动项，显式建模CRISPR干预效应
两阶段回归框架：先估计代理表达量，再通过QMLE识别父节点，配合Murphy-Topel标准误校正确保有限样本下的有效性

2. 核心模型与算法实现

2.1 数据生成模型与识别理论

ARGEN的输入数据包括：

单细胞UMI计数矩阵Y∈ℝ^(N×p)
gRNA扰动指示矩阵D∈{0,1}^(N×p)
技术协变量X∈ℝ^(N×J)

其核心测量模型延续scRNA-seq的经典假设：

Y_{ij} | λ_{ij}, ℓ_i ∼ Poisson(ℓ_iλ_{ij})

其中λ_{ij}代表潜在真实表达水平，ℓ_i为测序深度。关键创新在于对λ_{ij}的建模：

log μ_{ij} = θ_{j0} + ∑_{k∈pa(j)} log μ_{ik}θ_{jk} + τ_j D_{ij} + β_j^T X_i + γ_j^T U_i + ε_{ij}

这里μ_{ij}通过SEM编码了基因调控网络结构，D_{ij}反映CRISPRi干预效应，U_i代表未观测混杂。

识别性证明（Theorem 1-2）：

在"非退化直接效应"假设（θ_{jk}≠0, τ_j≠0）下，通过"排他性有向路径"概念证明祖先/后代节点可识别
当未观测混杂U与X独立或线性相关时，父节点系数θ_{jk}可通过代理变量η_k(D,X)=log E[Y_k|D,X,ℓ] - log ℓ识别

2.2 两阶段估计流程

阶段一：后代集合估计

通过检验假设：

H_0^{jk}: E[Y_k|D_j=1,X,ℓ,D_{-j}=0] = E[Y_k|D_j=0,X,ℓ,D_{-j}=0]

具体步骤：

对每个基因j，拟合GLM模型：
```
Y_k ~ offset(log ℓ) + 1 + D_j + X
```
使用受j扰动或未受任何干预的细胞数据
计算Wald统计量z_{τ*_jk} = τ*jk / √(V{τ*_jk})，获取p值
通过BH法控制FDR，确定显著后代集合des(j)

阶段二：父节点识别

代理构造：对每个k∈anc(j)，拟合：
```
η_k(D,X) = (1, D_k, D_{anc(k)}, X^T)^T ξ_k
```

QMLE回归：

Y_j ~ offset(log ℓ) + X + D_j + ∑_{k∈anc(j)} η_k(D,X)

标准误校正：采用Murphy-Topel方差估计：
```
V_MT = A_θθ^{-1}(B_θθ + A_θξ V_ξξ A_θξ^T)A_θθ^{-1}
```
解决代理变量估计误差传播问题

2.3 DAG搜索算法

算法1通过干预评分ν_j=|des(j)|实现：

按ν_j升序排列节点，从无后代节点开始构建
对每个节点π_{j+1}：
- 在anc(π_{j+1})中估计η_k
- 通过在线FDR控制确定父节点集合
- 添加π_{j+1}到图中并更新干预评分
保证输出为有效DAG（无环性）

3. 实际应用与验证

3.1 模拟研究：抗混淆能力验证

在p=8基因网络中测试ARGEN对未观测混杂和基因遗漏的鲁棒性：

系数估计偏差：当故意遗漏Y7,Y8时，ARGEN估计的θ_{jk}与全基因分析结果无显著差异（Wilcoxon p>0.05），而INSPRE和Naive GLM出现显著偏移
结构恢复：仅ARGEN能正确识别子图结构（图2a），其他方法产生虚假边
计算效率：8000细胞规模下可在2小时内完成分析

3.2 K562细胞数据应用

3.2.1 染色体内调控验证

对23条染色体分别构建DAG，平均每染色体40基因含100条边。关键发现：

方向性验证：98%的边满足CRISPRi预期——当父基因被抑制时：
- θ_{jk}<0对应子基因表达上升
- θ_{jk}>0对应子基因表达下降（图3d-e）
3D基因组支持：
- 调控基因对更可能位于同一TAD内（OR=2.1, p<0.001）
- Hi-C接触频率与ARGEN边强度正相关（r=0.33, p=3e-6）
- B区室（低表达区域）的调控信号更强（p<5e-3）
表观遗传证据：
- 6/8测试的TF（如SPI1、GATA2）在ARGEN边两端基因的启动子区显示共富集
- H3K27ac信号在调控对中显著增强（log2OR=1.8）

3.2.2 跨染色体调控网络

分析874个必需基因构建的全基因组网络显示：

核心调控枢纽：染色质修饰基因（如EP300）和转录因子（如MYC）处于网络顶层
通路模块性：核糖体蛋白基因形成紧密子网，与翻译起始因子相连
新型调控关系：发现线粒体基因TFAM对细胞周期基因CDK1的跨区室调控

4. 实操建议与注意事项

4.1 数据预处理要点

细胞质量控制：
- 保留UMI>500且<20%线粒体基因表达的细胞
- 每个扰动至少需要50个细胞以保证检验效能
技术协变量选择：
- 必须包含：批次、测序深度、细胞周期评分
- 建议用PCA检测潜在技术变异源
基因筛选：
- 在>10%细胞中表达的基因
- 扰动效率通过τ_j的显著性评估（FDR<0.1）

4.2 参数调优经验

代理模型选择：
- 高稀疏数据：负二项GLM
- 高测序深度：泊松GLM
- 过分散明显时：考虑零膨胀模型
FDR控制策略：
- 后代检测用BH法（α=0.05）
- 父节点识别用在线FDR（Zrnic方法），适应动态测试量
计算加速技巧：
- 对大型网络（>100基因）：
  - 先用Louvain聚类分模块
  - 并行处理各模块
- 内存优化：将表达式矩阵转为稀疏格式

4.3 结果解释陷阱

因果时效性：
- Perturb-seq捕获的是小时级响应
- 间接效应可能反映次级调控而非直接作用
细胞类型混杂：
- 未标注的亚群会导致虚假边
- 建议先用UMAP检查聚类，必要时加入聚类标签作为协变量
网络动态性：
- 静态DAG无法反映条件特异性调控
- 可分层分析（如按细胞周期阶段）

5. 方法比较与扩展方向

5.1 与传统方法对比

特性	ARGEN	INSPRE	NOTEARS
未观测混杂鲁棒性	✓✓✓	✓	×
干预数据利用	工具变量	干预目标	未使用
单细胞数据适配	测量模型	忽略技术噪声	高斯假设
计算复杂度	O(Np^2)	O(Np^3)	O(p^4)