HALO:用于单细胞多组学数据的分层因果建模.

0. TL; DR

HALO (Hierarchical cAusal modeLing)是一个从因果视角出发的单细胞多组学建模框架,核心假设是染色质可及性的变化是基因表达变化的原因,但这种因果关系存在两种模式:耦合(coupled)解耦(decoupled)

在多个真实的多组学数据集上,HALO不仅在分离耦合/解耦表征方面优于现有SOTA方法(如MIRA, GLUE, MultiVI),还成功地利用解耦表征识别出了具有“染色质潜能”的干细胞亚群、揭示了在疾病(如肺纤维化)和正常发育中,不同细胞谱系特异性的表观遗传调控因子、发现了与T细胞分化和功能相关的、新的远距离cis调控关系。

1. 背景介绍

在基因调控时,染色质可及性(ATAC)和基因表达(RNA)通常被认为具有同步耦合关系:染色质开放,转录随之启动。大多数单细胞多组学整合方法,如MultiVI、scMVP等,致力于寻找和对齐两种模态共享的、同步变化的信息。

然而,真实的生物学过程远比这复杂。

现有的整合方法大多忽略了这种异步解耦调控模式,限制了捕捉细胞状态转换、谱系决定等关键动态过程的能力。HALO 引入了因果建模的思想,将多组学数据分解为耦合与解耦两个部分,并分层级地进行深入探索,旨在揭示ATAC-RNA关系背后更深层次的调控逻辑。

2. HALO 模型

HALO框架通过一种自顶向下的方式,在表征层面(representation level)基因层面(individual gene level)ATAC-RNA的因果关系进行建模。

2.1 表征层面:基于因果约束的潜在空间分解

在宏观的表征层面,HALO的核心是一个经过特殊因果约束的变分自编码器(VAE)框架。其目标是将scATAC-seq数据 $X_A$ 和scRNA-seq数据 $X_R$ 分别分解为耦合与解耦两部分潜在表征。ATAC的潜在表征 $Z_A$ 被分解为:$Z_A = [Z_A^d, Z_A^c]$,RNA的潜在表征 $Z_R$ 被分解为:$Z_R = [Z_R^d, Z_R^c]$,其中上标 $d$ 代表解耦(decoupled),$c$ 代表耦合(coupled)

HALO引入了基于时间(或伪时间)的因果约束(Causal Constraints)。其基本假设是ATAC的变化是RNA变化的原因,但这种因果关系在耦合与解耦部分表现不同:

为了在模型中实现这些约束,HALO引入了一个基于核函数嵌入的统计量 $\Delta$,用于衡量两个分布随时间变化的依赖程度:

\[\Delta_{\hat{Z}_A \to Z_R} = \frac{\text{tr}(M_{Z_A} H M_{Z_R \mid Z_A} H)}{\sqrt{\text{tr}(M_{Z_A} H) \text{tr}(M_{Z_R \mid Z_A} H)}}\]

其中 $M$ 是格拉姆矩阵(Gram matrix)。这个 $\Delta$ 值越小,表示两个分布的变化越独立。因此因果约束可以被形式化为以下损失项:

通过将这些约束作为正则项加入到VAEELBO损失函数中,模型被强迫去学习满足上述因果假设的、解耦的潜在表征。

HALO还设计了一个可解释的解码器,它能将最终的重构结果分解为来自每个潜在维度(无论是耦合还是解耦)的加性贡献,从而能够直观地理解每个潜在维度所代表的具体生物学功能(如富集了哪些基因或通路)。

2.2 基因层面:耦合/解耦分数的计算与因果推断

在微观的基因层面,HALO旨在为每一个基因-peak对量化其耦合或解耦的程度,并对解耦现象进行因果推断。

首先,对于每个基因,使用负二项回归将其表达量与邻近一定距离内(如上下游600kb)的所有peaks的开放度进行关联,并根据回归系数和基因组距离进行加权,从而为每个基因匹配上一组最相关的局部peaks

对于匹配好的基因-peak对,再次使用上一节提到的依赖性度量 $\Delta$,来计算它们的表达谱($\tilde{X}_R(t)$)和聚合后的peak信号($\tilde{X}_A(t)$)在时间(或伪时间)上的变化依赖性。

\[\text{Decouple Score} = \text{sign}(\alpha - \Delta_{\tilde{X}_A \to \tilde{X}_R})(\Delta_{\tilde{X}_A \to \tilde{X}_R} - \Delta_{\tilde{X}_R \to \tilde{X}_A})\] \[\text{Couple Score} = -\text{Decouple Score}\]

当解耦分数为正时,我们认为该基因-peak对是解耦的。当耦合分数为正时,我们认为它们是耦合的。

对于那些表现出明显解耦行为的基因(即其局部peaks开放度随时间显著增加,但自身表达量却保持不变),一个自然的假设是这些开放的peaks可能不是在调控这个基因,而是在远距离调控(distal regulation)其他基因。为了验证这一假设,HALO引入了经典的格兰杰E因果检验(Granger Causality Test)。该检验的核心思想是:如果一个时间序列 $A$ 的过去值能够帮助预测另一个时间序列 $B$ 的未来值(在已知$B$自身过去值的条件下),那么我们就说$A$是$B$的“格兰杰原因”。

在我们的场景中,$A$是某个解耦peak的时间序列信号,$B$是其附近其他基因的表达时间序列。通过格兰杰因果检验可以系统性地筛选出那些由解耦peaks驱动的、潜在的远距离调控关系。

3. 实验分析

作者在多个具有复杂动态过程的真实多组学数据集上,对HALO的性能和生物学发现能力进行了深入的验证。

3.1 实验一:表征层面的耦合/解耦分析

作者首先在小鼠皮肤毛囊发育的SHARE-seq数据集上测试了HALO的表征分解能力。

HALO成功地将ATACRNA数据分解为耦合与解耦两部分。UMAP图显示,耦合表征($Z_A^c, Z_R^c$)在两种模态间展现出高度相似的结构,反映了共享的发育进程(图B, F)。而解耦表征($Z_A^d, Z_R^d$)则呈现出截然不同的模式,特别是在ATAC的解耦表征中,一个独特的细胞群被凸显出来(图C, G)。

深入分析发现,这个被ATAC解耦表征凸显出的细胞群,正是在原研究中被定义为具有高染色质潜能的新型根细胞(novel root cells)。这些细胞的染色质已经为分化做好了准备,但转录尚未启动,这正是典型的解耦现象。HALO通过其解耦表征,精准地捕捉到了这一关键的、具有谱系可塑性的干细胞亚群。

MIRA, GLUE, MultiVI等方法相比,HALO学习到的耦合表征相关性最高,而解耦表征相关性最低,定量地证明了其在分离两种信息上的优越性(图M)。

3.2 实验二:基因层面的耦合/解耦分析

作者进一步在基因层面,分析了不同谱系分化过程中的耦合与解耦基因。

HALO成功地识别出了在毛囊分化不同分支(皮层、内根鞘、髓质)中表现出特异性解耦行为的基因。例如,关键发育调控因子Dlx3在所有三个分支中都表现出解耦行为:其邻近的peaks开放度随时间持续增加,但其自身表达量却保持稳定(图E-G)。

对于Dlx3这种解耦基因,作者提出了一个假设:它邻近的那些开放的peaks可能不是在调控Dlx3,而是在调控其他基因。通过格兰杰因果检验,作者发现,这些位于Dlx3超级增强子区域内的peaks,其开放动态与邻近的另一个基因Itga3的表达动态存在显著的因果关系。这一发现揭示了一个新的、由超级增强子介导的远距离调控网络,而这是传统方法无法发现的。

作者发现在多个数据集中,处于分化过渡状态或干细胞状态的细胞(如毛囊中的TAC细胞、大脑中的V-SVZ细胞),其解耦基因的比例显著更高。而处于终末分化状态的细胞,其耦合基因的比例则更高。这一普适性规律表明,解耦是细胞可塑性的一个重要分子特征。

3.3 实验三:揭示疾病中的表观遗传调控

作者将HALO应用于系统性硬化症相关间质性肺病(SSc-ILD)的多组学数据,以探究疾病状态下的调控失常。

SSc-ILD中,肺泡II型(AT2)细胞的再生能力受损。HALOATAC解耦表征成功地将AT2细胞分为了两个亚群:一个仍能分化为AT1细胞(正常再生),另一个则异常地分化为分泌型细胞(纤维化相关)。进一步分析发现,后一个亚群的表观遗传状态(如TF基序活性、超级增强子富集)已经偏向于气道上皮细胞,揭示了其谱系失调的表观遗传基础。

作者发现,EMT主调节因子SOX4SSc-ILD患者的AT1细胞分化过程中表现出异常的解耦行为。通过格兰杰因果检验,进一步发现SOX4的局部peaks与一个已知的纤维化相关长非编码RNA CASC15的表达存在因果关系。这一发现为SSc-ILD的发病机制提供了一个全新的、表观遗传驱动的解释。