多组学背景下的基因调控网络推断

染色质、转录因子和基因之间的相互作用产生了复杂的调节回路，可以表示为基因调节网络（GRNs）。GRNs的研究有助于了解疾病中细胞身份是如何建立、维持和破坏的。GRN可以从实验数据——历史上的大量组学数据——或文献中推断出来。单细胞多组学技术的出现引领了新的计算方法发展，这些方法利用转录组，染色质可及性信息，蛋白质组，以单细胞分辨率推断GRN。这里主要回顾了从转录组学和染色质可及性数据推断包含转录因子-基因相互作用的GRN。主要强调了GRN推断中的挑战，特别是在基准测试方面，以及使用额外数据模态的未来进一步发展。

来自：Gene regulatory network inference in the era of single-cell multi-omics

---

TF（转录因子）是一类能够结合到DNA上调控基因表达的蛋白质（可以促进也可以抑制），而motif（基序）则是指在DNA序列中存在的具有特定模式的短序列，通常与TF的结合位点相关联。研究转录因子结合位点的基序可以帮助我们理解具体基因调控机制的物理过程（可以想象为建模配体-受体的相互作用）。

从scRNA-seq中获取TF：使用基因注释数据库将scRNA-seq数据中的基因与已知的转录因子基因进行比对和注释。基因注释数据库提供了基因的功能注释和分类信息，其中包括转录因子的注释。

---

背景概述

细胞调节基因转录以协调细胞对细胞内和细胞外信号的反应。转录在很大程度上受转录因子（TF，transcription factors）的调节，转录因子是一种与DNA特定序列（DNA结合位点，DNA binding sites）结合的蛋白质，可对靶基因的转录率产生积极或消极影响。DNA与结构蛋白紧密堆积在称为核小体（nucleosomes）的复合体中，核小体是染色质（chromatin）的基本单元。为了实现转录，需要通过置换紧密堆积的核小体来暴露基因转录起始位点附近的区域，即启动子（promoter）。DNA可及性的变化可以由所谓的pioneer TFs的结合触发。其他TF可以与DNA的末端顺式调节元件（CRE，cis-regulatory elements）结合，并与辅因子和其他蛋白质一起，协同实现从gene body DNA合成mRNA的RNA聚合酶-蛋白质复合物的募集和稳定（图1a）。

• 图1a：基因调控及其关键要素。转录因子（TF）与启动子区域和顺式调节元件（CRE）结合，取代核小体并使转录起始位点可接近。转录因子、辅因子和其他蛋白质之间的合作允许RNA聚合酶-蛋白质复合物的募集和稳定，该复合物从gene body DNA合成mRNA。

基因调控网络（GRNs，Gene regulatory networks）是以网络表示的基因调节模型，在数学上被定义为graph。基因调控的多种成分，如转录因子、剪接因子、长非编码RNA、微小RNA和代谢产物，都可以整合到GRN中。在这里，主要关注它们最简单的表示，即只捕捉转录因子和靶基因之间的相互作用，其中GRN的节点由基因组成，其中一些用基因代表转录因子，GRN的边代表基因之间的调节相互作用（图1b）。

• 图1b：基因调控网络（GRN）可以从测量的组学数据中推断出来，并通过对TF结合预测或染色质可及性等额外信息的建模，可以进行细化，以更好地近似真正的GRN。GRN的节点是TF及其调控基因，节点之间的边表示调控方式（激活或抑制，activation or inhibition）。

理解GRNs是生物学中一项长期的探索，如果有足够的转录组学数据可用，可以推断GRN以拟合手头的生物学问题。然而，转录组学数据并没有直接捕捉到许多潜在的调控机制，如TF蛋白丰度和DNA结合事件、TF和辅因子的合作、替代转录物剪接、翻译后蛋白修饰事件以及基因组的可及性。考虑这些其他方面有可能产生更好地代表体内基因调控的GRN（图1b）。例如，纳入染色质可及性数据可以通过考虑基因是否开放来修正TF-靶基因的连接。

此外，单细胞技术允许推断不同细胞类型、分化轨迹下的GRN（图1c）。因此，随着多模态分析技术的引入，最近出现了新的GRN推理方法。这篇综述概述了GRN推理的一般原理及其潜在的局限性。此外，描述了如何利用多模态信息来推断更准确的GRN，并对为此任务开发的几种新工具进行了分类和简要描述。

• 图1c：由单细胞多组学数据生成的GRN可以理解细胞类型，解释动态轨迹的进展，并识别条件之间的差异。

推断GRN

GRN推断是指使用计算方法从数据中将基因调控（一个高度复杂和动态的过程）总结为可解释的网络结构的过程。这是基于这样的假设：

• 可以在分子数据中观察和测量真正的潜在GRN的影响（图1b）。

GRN中的相互作用可以是有向的或无向的（分别表示基因之间的因果关系或缺乏因果关系）、有符号的（表示调节模式，阳性或阴性）和加权的（表示相互作用的强度）。

基于转录组数据

这一类方法适用于基于其他基因表达解释感兴趣观测基因表达差异的模型。加权基因共表达网络分析（WGCNA）是最简单和最流行的方法之一。它在整个转录组中进行成对相关性分析，以识别共表达基因。由此产生的网络通常被称为基因共表达网络，由于相关性计算是对称的，因此edge代表的相互作用是无方向的。尽管这种策略有助于以无监督的方式识别基因组合，但缺乏因果调控的联系阻碍了其可解释性，导致graph中存在大量的假阳性关联。为了解决其局限性，GENIE3及其更快的实现GRNBoost2等方法首先基于文献报道的调节活性（regulatory activity）将转录因子与靶基因区分开，然后训练仅基于转录因子表达预测靶基因表达的模型，这显著减少了需要考虑的相互作用的数量。通过这样做，无方向的相互作用变成了有方向的联系，从而引入了假定的因果关系。然而，仅从转录组数据推断依然会引入假阳性，因为许多其他参与基因调控的机制，如染色质可及性都被忽视了。此外，由于编码TF的mRNA转录物需要许多过程才能成为功能蛋白，因此仅转录物水平可能没有足够的信息量。这些限制都会阻碍推断GRN，总体而言，这些方法在准确推断GRNs方面只取得了适度的成功。

基于TF结合数据或染色质可及性

染色质免疫共沉淀测序（ChIP-seq，Chromatin immunoprecipitation followed by sequencing）和CUT&Tag（cleavage under targets and tagmentation，新一代ChIP-seq）使TF binding数据能够在整个基因组中进行测量。这些信息可以通过将TF结合位点分配给假定的靶基因来直接用于构建GRN（TF binding信息可以直接构建GRN）。然而，尽管有高通量的方案，TF binding的分析仍然是昂贵的，并且仅限于部分抗体检测的TFs（因为在ChIP-seq实验中，常用抗体来识别和捕获与转录因子结合的DNA片段）。此外，单独使用TF binding数据通常需要通过最接近的基因组将结合的TF分配给靶基因，这忽略已知与基因调控相关的可能的远端相互作用事件。

另一种方法是使用染色质可及性数据来推断转录因子可能靶向的基因调控元件。由于其简单且相对便宜的方案，最常用的技术是转座酶可及染色质测序分析（ATAC-seq），当然也存在其他技术，如DNase-seq和NOME-seq。利用染色质可及性数据的方法将GRN推断分为两个步骤：

• 首先，将转录因子分配给基因调控元件（开放染色质区域open chromatin regions，通常称为peaks）；
• 第二，将这些调控元件分配给基因（见图2）。

• 图2：基于染色质可及性的GRN推断流程：GRN推断的方法涉及不同的步骤。上述例子更符合配对数据场景。
• 转录组学数据首先经过预处理和标准化，以构建基因表达矩阵，包含样本或细胞中每个基因的转录水平。已知转录因子（TF）基因的列表是从其他来源获得的，以区分具有调节能力的基因（用基因代表TF）。然后通过构建模型来推断TF和靶基因之间的相互作用，该模型试图从TF对应的基因表达量预测矩阵中的其他基因表达，从而产生TF-基因关联。最后，所获得的交互作用被聚合并表示为GRN。
• 对于染色质可及性数据，首先对染色质可及性数据进行预处理，并用峰值来构建峰值可及性矩阵，该矩阵包含关于顺式调节元件（CRE）在样品或细胞中的开放性的二进制信息。CRE与基于基因组距离限制的基因相关，并且使用TF结合基序数据库和基序匹配算法预测TF与CRE结合。
• 这些信息一起用于获得"TF–CRE–基因"三元组。最后，这些相互作用被简化为"TF-基因"对，并聚集成GRN。当通过转录组学和染色质可及性（多组学数据）对样本进行分析时，需要对每种模态进行预处理，如果需要，还要对未配对的模态进行整合。

对于第一步，利用TF结合基序数据库（TF binding motif databases）和基序匹配器算法（motif matcher algorithms）对可访问CRE上的TF进行结合预测（方框1）。对于第二步，将可访问的CRE与一定基因组距离内的基因联系起来。距离截止是基于这样的假设：远端CRE（如增强子enhancers）通常以一定的距离与基因的启动子区域相互作用。简单来说，就是利用TF binding motif databases考虑染色质可及性信息筛选TF-CRE，再根据CRE获得附近基因，从而得到TF-gene。

这种基于ATAC推理的方法包括ATAC2GRN、LISA和SPIDER。这些方法假设：如果一个基因的启动子区域是可访问的，那么该基因正在被转录，但实际情况可能并不总是如此。

• 方框1：生成多种转录因子（TF）的全基因组结合数据需要很费力的实验，因此GRN推断方法倾向于基于数据库统计的TF结合事件。这些信息来自大量的TF–DNA结合测定，如ChIP-seq实验，可用于提取给定TF特异性结合的最可能的基因组序列。方框中的几个数据库已经收集了这样的测定，并生成了TF binding motif集合。由于数据库之间的覆盖范围可能不同，因此可以将它们合并以增加GRN推理过程中考虑的TF数量。
• 此外，已经开发了几种利用TF binding motif来预测结合事件的计算算法，称为基序匹配算法。所有这些算法都是基于从基序序列中计算TF结合事件的概率并过滤得到重要事件。由于不同的方法对TF结合的建模不同，因此它们之间的结果可能不同，在GRN推理过程中应该考虑这些结果。

关于ATAC中的二进制

二进制矩阵是指矩阵中的元素只有两种取值，通常是0和1。在scATAC-seq中，每个细胞的染色质开放区域可以被编码成一个二进制向量，其中1表示该位置在染色质中是开放的，0表示该位置是关闭的。

而有些scATAC-seq数据可能不是以二进制矩阵的形式存储的，可能是采用其他的数据表示方式，例如以整数形式表示染色质的开放程度或者使用浮点数来表示染色质的峰值强度等。这样的数据表示形式可以提供更多的细节信息，但可能需要更多的存储空间和计算资源。

基于单细胞转录组

使用大量组学数据的GRN推断方法能够表征全基因组调控事件，但在组织混合样本的情况下，难以根据数据捕捉GRNs，以及细胞类型或状态特异性。此外，GRN推理方法需要大量样本，这在bulk分析中可能会变得非常昂贵。

随着单细胞技术的出现，特别是单细胞RNA测序（scRNA-seq），GRN重建方法已被用于推断细胞类型特异的TF-基因相互作用，以及这些GRN在发育和条件下发生的动态变化（图1c）。第一种专门针对scRNA-seq数据的GRN推断方法是SCENIC，这是GRNBoost2方法的扩展，该方法通过利用TF-基因共表达产生细胞类型特异的GRN，此外，基于基因启动子区域上的TF binding motif富集来修剪GRN的边。

单细胞测量分辨率的提高也使得能够识别可能不容易分化为不同组的动态细胞状态及其转变，例如发育、细胞分化或疾病进展过程。Pseudotime的排序表征了这些连续的变化，并且可用于GRN推断。由此产生的GRN为关键命运决策所涉及的复杂过程提供了有价值的见解。LEAP和SINCERITIES是利用伪时间排序来做GRN推断的例子。

使用差异测试（differential testing）后获得的对比度统计数据是识别疾病之间差异的有效手段，例如健康个体和疾病患者队列之间的差异。

单细胞染色质可及性分析（scATAC-seq）的最新进展可以与单细胞转录组学一起进行，这使得GRN重建的精细化达到了无与伦比的定义。一些早期工作从未配对的多组学数据推断出GRN，以研究人类骨髓细胞分化、小鼠胚胎发育和树突状细胞的HIV感染。然而，他们并没有将自己的方法作为工具提供给其他人使用。随后，利用scRNA-seq和scATAC-seq的GRN推断新方法激增（表1和图2）。

• 表1：利用scRNA-seq和scATAC-seq的GRN推断方法。

• 表1-1：表1内容的补充。

如果两个测量值来自同一个细胞，则用于GRN推断的多模态数据是成对的，如果它们来自不同的细胞，则是不成对的（未配对）。一些方法不需要匹配每个细胞的染色质可及性和基因表达，因为它们要么总结跨细胞组的读数，要么为每个模态独立构建GRN，然后进行合并。相比之下，其他方法同时对同一细胞中的两种模态进行建模。在这些"同时"方法中，如果使用整合方法匹配两种模态，则仍然可以对未配对的数据进行建模。为了便于使用，其中一些方法（例如，DeepMAPS、FigR、GLUE、scAI和SOMatic）在GRN推断过程中实现了自己的整合方法。

对于单细胞数据，多模态GRN推理方法使用了一个可扩展的框架来重构GRN。具体而言，从TF基因表达预测基因表达，使用结合基序信息将TF分配给可访问的CRE，并将CRE与受基因组距离限制的靶基因关联（见图2）。GRN推断方法使用不同的基因组距离截断来将开放染色质区域分配给靶基因。一些人认为近距离高达10kb，另一些人认为中距离高达100kb，其他人认为远端效应高达1000kb，其他人则没有在原始出版物或源代码中指定距离截止值（见表1）。

在进行上述步骤后，多模态GRN推断方法生成由与连接到靶基因的CRE相关的TF的三元组组成的候选网络。为了生成最终的GRN，使用了不同的数学策略。其中一些策略假设TF、CRE和基因之间存在线性关系，而另一些策略则假设存在非线性关系（见表1）。线性模型假设一个变量，例如基因转录物，与另一个变量（例如TF转录物或CRE开放度）成正比变化。相比之下，非线性建模可以适应变量之间更复杂的相互作用，例如协同效应。尽管人们普遍认为基因表达是一个非线性过程，但GRN的线性建模往往是优选的，因为它在公式和解释方面很简单。独立于所使用的建模策略，可以使用频率论或贝叶斯概率统计框架来评估所获得的调节相互作用的重要性（见表1）。频率论方法将事件的概率定义为事件在大量相同实验中发生的次数比例，而贝叶斯概率将其定义为基于观测数据和先前信息对所述事件发生的置信度的度量。贝叶斯方法可以考虑可用的先验知识，但它们通常比频率论方法需要更大的计算资源，这在用大规模单细胞数据推断全基因组GRN时可能是一个限制。此外，贝叶斯推理的成功与否取决于所使用的先验知识的质量。因此，当没有可用的先验信息或怀疑其不准确时，频率推断可能更准确。

多模态GRN推理方法可以根据其建模策略和接受的输入类型的组合进行分组（表1）。大多数方法旨在通过频率回归对不同群体（通常是细胞类型）的GRN进行建模。FigR和GRaNIE等使用了频率线性回归；DIRECT-NET和SCENIC+使用频率非线性回归（随机森林）；PECA和Symphony使用贝叶斯建模。相比之下，CellOracle、Infereator 3.0和Pando为用户提供了多种建模策略。在由于数据的连续性而无法从数据中定义不同的组的情况下，例如在细胞发育中，scMEGA和IReNA利用轨迹分别线性和非线性地推断GRN。此外，Dictys、scMTNI和TimeReg使用细胞类型分组和轨迹数据的组合来为GRN建模提供信息，而CellOracle和SCENIC+使用后者来进行下游分析。

GRN的下游分析

GRN支持各种下游分析，以提供新的生物学见解。

拓扑分析

尽管GRN是简单且可解释的基因调控模型，但graph中仍然可以包含大量基因以及它们之间更大数量的相互作用。网络中心性度量可以帮助确定哪些TF或基因对网络的连通性或信息流更重要（图3a）。网络中心性度量的一些例子包括度中心性、贴近度中心性，介数中心性和特征向量中心性。这些度量有助于识别在不同生物环境中驱动细胞命运变化的转录因子。

表征GRN拓扑的另一种方法是使用基于谱图理论的方法，该方法探索网络在表示为矩阵时的性质。例如，应用于GRN的邻接矩阵的非负矩阵分解已经确定了协同驱动小鼠胚胎干细胞谱系转变的TF组。类似地，GRN拓扑结构的聚类确定了人类造血细胞分化和巨噬细胞对干扰素-γ反应中的已知调节因子。然后可以对获得的基因调控模块进行基因集富集，以表征其潜在的生物功能。

• 图3a：拓扑分析。网络中心性度量可用于识别转录因子（TF）中枢或基因调控网络（GRN）内高度连接的基因。基于节点连接性的节点聚类产生了可以与生物功能相关联的子网络模块。

比较分析

GRN的比较分析可以揭示导致细胞类型、细胞状态、疾病状态、治疗方法和生物体之间差异的事件（图3b）。比较分析最简单的方法是成对减去GRN之间的TF-基因相互作用。该方法已经确定了淋巴细胞白血病患者B细胞亚群中的关键调节因子、用于将成纤维细胞转分化为不同人类细胞类型的TF组、候选阿尔茨海默病特异性转调节因子和人类T细胞中的细胞状态特异性调节因子。它还被用于评估TF-基因相互作用的进化保守性和跨物种转录调节的适应性。然而，由于GRN的稀疏性和噪声性，TF-基因相互作用的直接比较往往不够稳健。

主题建模策略（Topic modelling strategies），如潜在狄利克雷分配（latent Dirichlet allocation，来自：Inference of Population Structure Using Multilocus Genotype Data），一种最初为自然语言处理开发的无监督贝叶斯模型，允许生成密集的低维表示，过滤GRN结构中的噪声，从而更稳健地捕捉调控关系中的差异。这一策略对预测癌症患者的生存率和识别人类造血系统中的rewriting events非常有用。

• 图3b：比较分析。通过成对减去GRN之间的TF-基因相互作用来比较不同GRN中的连接性，可以深入了解不同细胞类型、个体、条件或生物体之间的基因调控rewriting。

TF活性推断

GRN可以与富集方法相结合，从转录组数据推断TF活性（比如：decoupleR: ensemble of computational methods to infer biological activities from omics data）。这种方法允许将观察到的基因表达与GRN拓扑结构整合，以提取哪些转录因子在某些情况下可能具有相关作用（图3c）。常见的富集方法包括GSEA、AUCell和VIPER等。在大量研究中，通过富集方法推断TF活性使得例如能够鉴定可药用癌蛋白、对药物治疗的细胞系分层以及鉴定作为乳腺癌转移启动子的主调节因子。

在单细胞研究中，富集方法已经确定了人类T细胞免疫疗法耐药性的新机制、少突胶质瘤的调节因子和诱导因子，以及新冠肺炎患者病理性成纤维细胞的潜在药物靶点（来自：A molecular single-cell lung atlas of lethal COVID-19）。这些方法最近也被应用于空间分辨转录组学数据，例如，提出参与心肌细胞在人类心肌梗死中缺血病变周围边界区的功能转换的调节因子（来自： Spatial multi-omic map of human myocardial infarction）。

• 图3c：TF活性推断。GRN可以与富集方法相结合，从转录组学数据中推断哪些转录因子可能具有功能活性。从多组学数据推断出的GRN可以用于推断其他情况下的TF活性，如空间转录组或bulk转录组。

细胞命运扰动预测

GRN可用于通过以迭代方式将TF表达传播到靶基因来模拟随时间的基因表达值。有了这个框架，可以通过改变候选TF的表达，然后观察它在给定次数的迭代后如何影响最终的转录组来进行计算机扰动（图3d）。之后，可以将模拟值与局部邻近细胞的基因表达进行比较，以类似于RNA速度分析来估计细胞身份转换概率。该策略首先由CellOracle引入，提出了Zfp57在产生和维持小鼠诱导的内胚层祖细胞中的作用，后来通过体外扰动实验进行了实验验证。SCENIC+使用类似的策略将RUNX3鉴定为黑色素细胞到间充质黑色素瘤细胞的潜在驱动因素，展示了GRN捕捉和模拟复杂调节事件的能力。

• 图3d：在计算机扰动实验中。GRN可以通过在短迭代时间内通过网络传播基因表达的变化来模拟扰动。然后可以使用所获得的模拟基因表达谱来推断细胞命运发展。

多组学GRN推断的未来

转录组学和染色质可及性数据的配对分析使GRN的推断可能更准确，但仍然是一种成本高昂的测定方法，限制了其广泛使用。ISSAAC-seq等较新的替代品能够以比商业10×Multiome试剂盒低得多的成本进行多组学测量。尽管如此，单独的scRNA-seq和scATAC-seq联合数据可能无法提供足够的信息来完全表征基因调控。在这种情况下，包括更多数据模态的单细胞多组学分析技术的进步将是至关重要的。

在这些有前景的技术中，NEAT-seq可以同时分析核内蛋白、染色质可及性和基因表达，通过包括TF蛋白丰度，可以在GRN建模中丢弃可能的假阳性。另一个例子是scChaRM-seq，它同时描述DNA甲基化、染色质可及性和基因表达。利用它们的联合分析允许根据其甲基化状态对CRE的TF分配进行微调。此外，ATAC-STARR-seq可以同时进行大规模基因测定和染色质可及性测定。

非靶向单细胞蛋白质组学和磷酸蛋白质组学的进展使功能活性TFs的分析成为可能。一个例子是Phospho-seq，这是一种在单细胞水平上分析染色质可及性和磷酸化蛋白质的新技术。众所周知，遗传信息在个体群体中是异质的，但大多数方法都假设它们共享相同的基因组。