MIRA：单细胞多模态表达和染色质可及性的联合调控建模.

0. TL; DR

MIRA (Probabilistic Multimodal Models for Integrated Regulatory Analysis)包含了两个相互关联的概率模型，从全局和局部两个层面系统性地整合多组学信息：

  1. 细胞层面（全局）- 主题模型 (Topic Modeling): 利用一个基于VAE的主题模型从scRNA-seq和scATAC-seq数据中分别学习到定义细胞身份的表达主题和可及性主题。
  2. 基因层面（局部）- 调控潜能模型 (Regulatory Potential, RP Modeling): 针对每个基因座构建一个调控潜能模型量化基因周围的cis-调控元件对其表达的影响。

基于RP模型，MIRA能够区分两类基因：

  1. LITE基因 (Local accessibility-influenced Transcriptional Expression): 其表达主要受局部染色质可及性紧密调控的基因。
  2. NITE基因 (Non-local accessibility-influenced Transcriptional Expression): 其表达更多地受到非局部信号（如远端增强子或转录因子丰度）的影响，而与局部染色质状态相对解耦。

作者将MIRA应用于表皮维持分化和胚胎大脑发育两个系统。分析发现，在这些发育过程中，早期发育基因倾向于被LITE机制紧密调控，而终末命运基因则倾向于被NITE机制调控。

1. 背景介绍

在同一个单细胞中同时分析基因表达（expression）和染色质可及性（accessibility），为我们理解驱动细胞沿发育轨迹分化的转录和表观遗传机制的相互作用，提供了前所未有的机会。

近年来，涌现出许多计算方法来联合分析这些多组学数据。它们大多通过将细胞投影到一个基于两种模态的共享潜在空间（shared latent space），来更好地捕捉数据的生物学结构。然而，该领域仍然缺乏能够超越可视化和聚类，去严谨地对比这两种模态，以揭示驱动发育命运决定的复杂调控回路（regulatory circuitry）的工具。

一个理想的整合分析框架应该能回答更深层次的问题：

2. MIRA 方法

MIRA通过主题模型（topic modeling）和调控潜能（Regulatory Potential, RP）建模，旨在系统性地对比转录和可及性信息，以确定驱动细胞沿发育连续体分化的调控回路。其核心是其独特的双层概率模型，它分别在细胞层面和基因层面对多组学数据进行建模和整合，从而实现从宏观到微观的全面分析。

2.1 细胞层面：基于主题模型的联合表征

在细胞层面，MIRA的目标是学习一个能够捕捉细胞身份和状态的、可解释的低维表征。它引入了概率主题模型（Probabilistic Topic Modeling）的思想。

主题模型源于自然语言处理，它旨在从一堆文档中发现抽象的主题。在单细胞领域，一个细胞可以被看作一篇文档，一个基因或一个peak可以被看作一个单词。MIRA利用一个基于变分自编码器的方法，从scRNA-seq和scATAC-seq数据中，分别学习出$K$个表达主题和$K$个可及性主题。每个主题都是一组共调控的基因或共开放的cis-调控元件。每个细胞的分子谱则被表示为这些主题的一个稀疏组合。

MIRA为两种模态设计了不同的生成分布，以适应它们各自的统计特性。例如，为scRNA-seq的过离散计数和scATAC-seq的稀疏数据设计不同的分布）。模型中加入了一个稀疏性约束，以确保每个细胞只由少数几个主题构成，从而保证了主题的可解释性。模型通过贝叶斯超参数优化来自动学习最合适的主题数量。

学习完成后，MIRA将表达主题和可及性主题的细胞组成（cell-topic compositions）进行整合，形成一个联合表征。这个联合表征被用于构建一个K-近邻图（KNN graph），并进一步用于细胞可视化和聚类、构建高保真度的谱系轨迹（lineage trees）、在轨迹的分叉点进行主题分析，以确定驱动命运决定的关键调控因子。

2.2 基因层面：调控潜能 (RP) 模型

在基因层面，MIRA的目标是定量地评估每个基因周围的局部染色质状态对其表达的影响。

作者借鉴了调控潜能（Regulatory Potential, RP）建模的思想。对于一个给定的基因，其RP被定义为其附近所有cis-调控元件（peaks）可及性的加权总和。MIRA的RP模型的一个关键创新是，它假设调控元件的影响力会随着其与基因转录起始位点（TSS）的基因组距离呈指数衰减。模型为上游和下游学习独立的衰减速率，并包含额外的参数来权衡启动子区（promoter）、上游和下游区域的贡献。

为了量化局部染色质对基因表达的调控程度，MIRA构建了两个嵌套的预测模型：

通过比较这两个模型的预测能力，MIRA可以将基因分为两类：

MIRA定义了一个染色质差异(Chromatin Differential)指标，即LITE模型预测值与NITE模型预测值之间的差异。这个指标量化了在每个细胞中，一个基因的表达与其局部染色质状态解耦的程度。

3. 实验分析

作者将MIRA应用于两个来自不同平台（SHARE-seq和10x Genomics）的真实发育系统：小鼠毛囊分化和胚胎大脑发育。

3.1 实验一：毛囊分化（SHARE-seq数据）

MIRA的联合主题表征构建的UMAP图，其潜在结构与毛囊的真实空间布局高度吻合。推断出的谱系树准确地重建了从外根鞘细胞到基质祖细胞，再到内根鞘、髓质和皮质的祖先层级关系。

通过在谱系树上可视化主题的动态变化（stream graphs），MIRA成功地识别出了驱动髓质（medulla）和皮质（cortex）命运决定的关键调控模块。例如，髓质特异性的可及性主题a5富集了Smad家族的基序，而皮质特异性的表达主题e6则与Notch信号通路相关，这与已知的生物学知识完全一致。

作者发现一个有趣的现象：在髓质和皮质命运分叉点，存在一个短暂的、由Wnt信号驱动的可及性主题a4。这个主题的出现，早于任何一个谱系特异性表达主题的激活。这表明，在这些祖细胞中，全基因组范围的染色质重塑，先于决定命运的转录程序的启动。

3.2 实验二：LITE与NITE基因的区分

MIRA的RP模型成功地将毛囊中的基因分为了LITE和NITE两类。例如，Braf是一个典型的LITE基因，其表达与局部染色质可及性高度同步；而Krt23则是一个典型的NITE基因，其表达与局部可及性解耦，但与一个全基因组范围的髓质特异性可及性主题a5高度相关。

作者发现，许多终末分化基因都属于NITE基因，并表现出一种分支点启动（branch-primed）的模式：在命运决定之前的分叉点，它们周围的染色质就已经变得开放（即被启动）；而在命运决定之后，它们只在特定的谱系中被高表达（被激活）。

MIRA的pISD分析进一步推断出，皮质谱系中启动基因的激活，是由Notch信号通路的效应因子Rbpj驱动的；而髓质谱系则是由Bmp/Tgf-β信号通路的效应因子Smad家族驱动的。

3.3 实验三：胚胎大脑发育（10x Multiome数据）

作者进一步将MIRA应用于胚胎大脑发育数据，以探究星形胶质细胞、兴奋性神经元和抑制性神经元这三种命运的调控机制。

MIRA再次成功地重建了发育轨迹，并识别出驱动不同命运的关键TF。例如，Notch信号驱动星形胶质细胞命运，Neurog2和Neurod1驱动兴奋性神经元命运，而Ascl1则驱动抑制性神经元命运。

同样，作者也观察到了LITE和NITE基因的差异分布。Pax6等早期命运决定TF属于LITE基因，受局部染色质紧密调控；而神经递质和离子通道等在终末神经元功能中起作用的基因则属于NITE基因，其表达可能需要对外界信号做出快速响应，因此与局部染色质状态解耦。