基于单细胞RNA 测序数据的基因调控网络推断算法综述

张少强，潘镜伊

（天津师范大学计算机与信息工程学院，天津 300387）

基因调控网络（gene regulatory networks，GRN）[1]由影响生物体生物过程的调控因子间的相互作用关系组成，是不同转录因子（transcription factors，TF）或辅因子（co-factors）与其靶基因或靶转录本之间相互作用的图谱.GRN 通常表示为连接调控组件与目标组件的图，通过网络建模，能够表示各组件如何相互作用以及它们可能产生的行为.虽然基因间相互作用是动态的，但当前许多生物学研究侧重于基因的定性分析或基因对（gene pair）间的相互作用，即研究一个基因扰动如何影响另一个基因表达.因此，通过基因表达的变化可以推断GRN；通过分析不同细胞类型在不同生物过程中的基因表达关系可以研究基因间的动态相互作用.了解每个基因的作用及其与其他基因的关系是在分子水平上解释生物过程的关键.因此，重构和分析GRN 已成为研究基因型驱动表型潜在机制的有效方法，并有望应用于医学和药物设计等领域[2].

GRN 具有2 个重要的特征[3]：①网络拓扑结构，它表示网络连接模式.某个基因的逻辑结构可以表示为元素是0 或1 的向量，元素的取值代表该基因与其他基因的关系.②组件之间的调控作用，即相互作用的强度和类型，该作用由一组权值来描述，权值为正表示基因的正调控（基因被激活），权值为负表示基因的负调控（基因被抑制），权值为0 表示无调控关系.

单细胞RNA 测序（single-cell RNA sequencing，scRNA-seq）为推断细胞周期或分化等时间依赖性生物过程的GRN 提供了新的可能，并且基于scRNAseq 产生了大量的GRN 推断算法[4].与批量RNA 测序（bulk RNA-seq）数据不同，scRNA-seq 没有混合所有细胞样本的基因表达而掩盖特定细胞的生物信号；另外，scRNA-seq 不仅能够提供静态的细胞快照数据，也能够提供细胞不同时序（time-series）的动态数据.因此scRNA-seq 数据更适合于GRN 推断[4].本文综合分析了26 种基于scRNA-seq 数据的GRN 推断算法，其中包括3 种早期基于批量RNA 测序数据的设计方法，但也可以用于scRNA-seq 数据.根据不同的网络构造方式和模型，另外23 种GRN 推断算法大致分为6类：①基于布尔网络（Boolean network）的算法2 种[5-6]；②基于微分方程（differential equations）的算法3 种[7-9]；③基于伪时序（pseudo-time-series）基因相关性集成策略的算法5 种[10-14]；④基于共表达（co-expression）基因的算法4 种[15-18]；⑤基于细胞特异性（cell specific）的算法3 种[19-21]；⑥基于深度学习的算法6 种[22-27].详细描述了每类算法的工作机制或方法原理、适用性、主要优势及局限性等，对相关的算法比较研究进行分析，并利用scRNA-seq 数据简单验证了这些推断算法重建网络的精确度，最后分析了当前GRN 推断算法面临的主要局限和突出挑战，并建议新的研究思路.

1 GRN 推断算法

1.1 早期基于批量RNA 测序数据的算法

基于批量测序数据重构GRN 的算法主要包括3种，分别为WGCNA（weighted correlation network analysis）[28]、GENIE3（GEne network inference with ensemble of trees）[29]和GRNBoost2[30].WGCNA 使用层次聚类和动态树分割识别高度相关（共表达）的基因簇（模块），利用模块特征基因或模块内的中心基因表征这些模块，此外，将模块内基因相连和因果关系与外部信息（如基因功能语义、功能富集和蛋白质组学等）进行关联来寻找模块中的关键驱动基因.GENIE3 是一种基于特征选择和回归树集成（tree-based ensemble）的GRN推断算法.GENIE3 算法的基本思想是将p 个基因的GRN 问题分解为p 个不同的回归问题，目标是确定网络中目标基因的调控因子，在每个回归问题中，通过树集成的方法（如随机森林），基于所有其他基因的表达模式预测得到目标基因的表达模式.GRNBoost2是基于GENIE3 架构的GRN 推断算法，是GENIE3的快速替代方案，可扩展至数万样本规模的数据集.与GENIE3 类似，GRNBoost2 训练一个回归模型为数据集中的每个基因选择最重要的调控因子，并利用启发式早停正则化的随机梯度提升机（gradient boosting machine，GBM）回归来推断网络.

WGCNA 利用数据探索性工具或基因筛选（排序）构造加权和未加权的无向网络，GENIE3 和GRNBoost2 构造允许存在反馈回路的有向网络.WGCNA 通过增加可测试的假设（如某模块与某疾病相关），以便在独立的数据集中进行验证，而GENIE3 没有对基因调控的性质做任何假设，可以潜在地处理组合调控和非线性关系.GENIE3 和GRNBoost2 在已知大量基因表达的情况下运行良好，计算速度快，且可扩展，WGCNA 没有考虑输入数据是否被正确预处理或规范化，所以有可能导致输出结果存在偏差或无效.此外，WGCNA 虽然实现了几种共表达模块的检测方法，但没有明确其中的最好方法.GENIE3 在基因排序中没有考虑树的质量，所有树模型的权重取值相同.这3种基于批量测序数据的GRN 推断算法也可以用于单细胞测序数据，但在各种单细胞实验中性能表现不佳，其主要原因是单细胞基因表达数据存在大量缺失（dropout）零值，零值率大大高于批量样本数据.相关实验结果也表明，在没有dropout 的模拟单细胞数据的实验中，GENIE3 和GRNBoost2 的表现优于很多通用方法，但在有dropout 的模拟数据的实验中，其表现较差[31].

1.2 基于scRNA-seq 数据的算法

单细胞测序技术能够在细胞分化过程中追踪细胞谱系并识别新的细胞类型.此外，单细胞测序技术使得控制细胞分化并驱动其细胞类型转变的GRN 成为可能.但与批量数据相比，scRNA-seq 较大的噪声、细胞异质性、细胞间测序深度的差异、dropout 引起的基因表达的高稀疏性和细胞周期相关效应等都为GRN 推断带来新的挑战.

基于scRNA-seq 数据的GRN 推断算法的总体工作流程如图1 所示.输入为某scRNA-seq 数据集的基因表达矩阵，行和列分别代表基因和细胞.首先进行基因过滤，将分析范围缩小到具有高度可变性的基因或用户感兴趣的基因（预定义基因）；然后根据数据假设建模和构造中间数据（图1 给出了6 类建模方法的示意图）；最后推断网络，输出结果可以是无向的共表达网络，也可以是具有基因间调控关系的有向网络.

图1 基于scRNA-seq 数据的GRN 推断算法的整体工作流程Fig.1 Overall workflow of GRN inference algorithm based on scRNA-seq data

1.2.1 基于布尔网络的算法

布尔网络模型是最简单的网络推断方法之一，一个布尔网络模型由n 个基因x1，…，xn和n 个基因关联更新函数f1，…，fn：{0，1}n→{0，1}构成.每个基因取一个布尔值x∈{0，1}表示基因表达是否存在；每个更新函数fi用布尔逻辑表示，使用布尔算子AND、OR 或NOT 指定基因x1，…，xn间的关系.布尔网络的GRN推断算法一般分为3 步：①对基因表达数据进行二值化，生成初始布尔状态；②用布尔状态优化器进行模型状态优化；③输出具有激活边（edge）、抑制边和一组布尔函数的GRN.应用布尔网络推断基因间关系的算法主要有SCNS（single cell network synthesis）[5]和BTR（BoolTraineR）[6].

SCNS 以跨时间段（时序）单细胞基因表达数据作为输入，通过构建布尔网络模型搜索从初始细胞状态到后期细胞状态的进展和转化的逻辑规则.SCNS 产生的逻辑模型有助于预测基因扰动（如基因敲除或基因过表达）对特定谱系的影响.对于几千个细胞的数据集，SCNS 在普通计算机上几分钟内即可重建布尔网络模型.对于较大规模的数据集，SCNS 支持云计算和并行计算，能够轻松部署到云或高性能计算服务器上.BTR 算法是一种异步布尔模型，它使用新型布尔状态空间评分函数训练单细胞表达数据.BTR 通过改进模型预测与表达式数据的匹配来细化现有的布尔模型，并重构新的布尔模型.BTR 对细胞间表达状态的关系不做任何假设，这有助于在条件尽可能少的情况下重构GRN.BTR 也可用于优化现有假定的调控网络.

BTR 将网络作为图对象输出，而SCNS 只输出优化后的布尔规则.SCNS 依赖已知细胞状态的一般轨迹，这要求单细胞表达数据至少来自2 种具有已知关系的细胞类型，而BTR 不需要细胞状态的轨迹信息来推断网络结构和布尔规则.布尔网络对dropout 具有一定的鲁棒性，但是布尔网络将基因表达信息二值化会降低网络预测的准确性.此外，这类算法需要通过状态转换优化布尔规则，分析中包含的基因越多，构建连接状态转换图所需的细胞就越多.

1.2.2 基于微分方程的算法

基于微分方程的算法将基因表达动力学描述为与其他基因表达和环境因子相关的函数.在给定一组参数的情况下，微分方程可以表达数据随时间的变化率.使用微分方程进行网络推断需要细胞的时间顺序，因此这类算法适用于推断细胞分化期scRNA-seq 数据的GRN.这类算法主要包括Inference Snapshot[7]、SCODE[8]和SCOUP（single-cell expression data during differentiation with Ornstein-Uhlenbeck process）[9].基于微分方程算法的工作流程如图2 所示，分为4 步：①利用外部软件或内嵌函数推断细胞的伪时序；②使用微分方程描述基因与时序间的关系；③使用不同的优化技术估计参数；④利用优化的参数通过微分方程推断基因间的关系，最后输出基因亲和度矩阵作为GRN.

图2 基于微分方程的GRN 推断算法的整体工作流程Fig.2 Overall workflow of GRN inference algorithm based on differential equations

Inference Snapshot 利用模型选择与参数优化相结合的方法构建解释伪时序的常微分方程.在构建常微分方程模型的过程中，为减少模型数量，首先用GENIE3 生成粗糙的GRN 结构作为先验知识，然后重建分化通路中的基因表达动态，并推断出关键的GRN 结构.SCODE 以微分方程对表达数据的最优拟合为目标，优化转录因子间调控关系的矩阵，通过整合线性微分方程变换和线性回归实现矩阵优化.SCODE 只使用少量转录因子即可重建表达动力学，从而显著降低了时间复杂度.SCOUP 使用随机微分方程直接描述细胞分化程度和细胞生命整个分化过程中的基因表达动力学.SCOUP 首先构造从初始细胞到所有其他细胞的最小生成树以获取细胞的伪时序；然后将每个基因的表达随伪时序的分化建模为一个连续的随机扩散过程，即Ornstein-Uhlenbeck（OU）过程，其中单个基因在某一时间点的表达通过当前OU 过程的正态分布来估计；最后通过计算所有细胞的Z 值得到基因间的相关性.

这3 种基于微分方程的GRN 推断算法存在一些技术差异.在数据降维和构建伪时序的过程中，Inference Snapshot 通过非线性降维方法diffusion map[32]降维输入数据后调用Wanderlust 算法[33]生成细胞伪时序；SCODE 使用轨迹推断方法Monocle[34]降维数据并寻找最小生成树轨迹来构造细胞伪时序；SCOUP 使用主成分分析方法进行数据降维并采用Prim 算法[35]寻找最小生成树的最短路径作为细胞伪时序.另外，这3 种算法中只有SCODE 假设基因间的依赖关系是线性的.在算法输出结果中，SCODE 和SCOUP 输出代表基因关系的相关矩阵，并可以调用第三方网络可视化应用以实现可视化，而Inference Snapshot 只输出模型的估计参数.

由于微分方程可以模拟RNA、蛋白质等在GRN中与时间有关的各种相互作用的参数变化，并可以推断基因间的某些因果关系，因此比其他方法更接近实际网络和细节.然而，由于微分方程的推断算法需要输入大量数据来估计参数且计算量较大，因此只能应用于有限数量的基因调控关系.

1.2.3 基于伪时序基因相关性集成策略的算法

与基于微分方程的算法一样，基于伪时序基因相关性集成策略的推断算法也需要推断细胞的伪时序或者需要用户提供细胞时序，通过时序基因表达数据推断出动态GRN[36].这类算法假设基因间的关系随着细胞的不同发育阶段而变化，并将时序数据分成较短的时间窗口，然后计算每个时间窗口的基因相关性，最后运用集成策略将多个相关矩阵聚合成单一基因邻接矩阵.这类算法主要包括LEAP（lag-based expression association for pseudo-time-series）[10]、AR1MA1-VBEM（first-order autoregressive moving-average）[11]、SINCERITIES[12]、SCIMITAR（single-cell inference of morphing trajectories and their associated regulation）[13]和SINGE[14].

LEAP 以伪时序数据为输入，首先计算不同时滞的固定时长窗口读数（read counts）的皮尔逊相关性；然后将每对基因间的相关强度设置为所有时滞中皮尔逊相关性绝对值的最大值；最后通过置换检验估计错误发现率得到最终的基因相关矩阵.由于计算的相关性不是对称的，因此该算法可以输出有向网络.AR1MA1-VBEM 的输入是单细胞时序数据或使用模块最优叶序算法得到的伪时序，然后构造一阶自回归移动平均模型来拟合基因表达数据，并在变分贝叶斯（variational-Bayesian）框架内推断GRN.SINCERITIES利用时序基因表达分布的时间变化，通过正则化线性回归（岭回归）恢复基因间的定向调控关系，同时通过成对基因间的偏相关分析预测基因调控的激活和抑制模式.该算法的计算复杂度较低，能够进行大规模基因的预测.SCIMITAR 首先通过在高维曲线中使用连续参数化（时间点作为参数）的渐变高斯混合模型（morphing Gaussian mixture model）推断静态scRNA-seq数据的伪时序轨迹，同时考虑了数据中的异方差噪声，并加入了检测轨迹进展中相关基因的统计能力.SCIMITAR 从数据中标记出与时序轨迹进展相关的共表达模式基因，通过计算协方差矩阵间的距离获得基因相关矩阵.针对伪时序的不均匀分布，SINGE 使用基于核的Granger 因果回归来平滑不规则的伪时序和补全缺失的表达值来提高GRN 推断的鲁棒性.SINGE输入伪时序的单细胞基因表达数据，输出预测的调控因子和靶基因相互作用的排序列表，利用不同的超参组合的多个广义Lasso-Granger 检验来控制网络的稀疏性、核平滑和伪时序分辨率等，最后将所有广义Lasso-Granger 检验得到的部分网络聚合成整体GRN.

在这5 种算法中，SCIMITAR 和AR1MA1-VBEM可以从基因表达数据中预测时序，其他3 种算法均需要输入单细胞时序数据.这类GRN 推断算法将不同细胞命运决定的多分支细胞轨迹强制合并为一条线性发展轨迹，可能会影响网络的准确性.此外，这些算法利用皮尔逊相关性或Granger 因果关系来推断基因间的关系，但由于生物系统的复杂性，很可能产生非线性的基因相互作用，这也会影响网络的准确性.在5 种算法中，LEAP 因计算相关性比其他算法更加简单而具有明显的优势.

与基于微分方程的算法一样，这类算法的性能严重依赖于细胞时序的准确性.时序信息通常需要数据提供，或者通过伪时序轨迹推断方法得到.大多数伪时序轨迹推断方法需要提供某些先验信息，如起始细胞、结束细胞和分支数，若无法获得这些信息，则很难建立可靠的时间轨迹，所以这类推断算法不适于处理多分支细胞轨迹.

1.2.4 基于共表达基因的算法

基于共表达的GRN 推断算法需要调用某种相关性度量来计算两两基因间的关系，其整体工作流程如图3 所示.首先通过计算每对基因的表达相关性初始化边的权值；然后通过假设检验估计每条边的显著性，使用预定义的阈值去除被认为不显著的边；最后输出最大连通分支.这类算法主要包括NLNET[15]、PIDC（PID and context）[16]、SINCERA（single-cell RNA-seq profiling analysis）[17]和SCENIC（single-cell regulatory network inference and clustering）[18].

图3 基于基因共表达的GRN 推断算法的整体工作流程Fig.3 Overall workflow of GRN inference algorithm based on gene co-expression

细胞之间存在着非线性的复杂关系，因此可以利用非线性依赖关系的方法推断网络.NLNET 是一种基于非线性关联的网络重建算法，利用基于条件有序列表的距离[37]（distance based on conditional ordered list，DCOL）计算每对基因间关系，具有较高的灵敏度和计算效率.PIDC 使用信息论中的部分信息分解（partial information decomposition，PID）来确定基因间的关系.PIDC 将一对基因X 和Y 提供的关于目标基因Z 的信息划分为冗余信息、协同信息和独有信息.首先计算提供给Z 的X 独有信息和X 与Y 的互信息（mutual information，MI）的比率；然后计算为所有其他基因提供的X 和Y 的独有信息比率的平均值，该平均值被称为比率独有贡献（proportionaluniquecontribution，PUC）.为量化PUC 值的置信度，对每个基因估计其零假设下PUC 值的分布，然后计算其置信分数.由于X 和Y 的PUC 值是对称的，因此产生的网络是无向的.SINCERA通过给定其他基因的变量，利用一阶条件相关性度量一对基因调控相关的重要性，即给定第3 个基因Z，基因X 和Y 的相关性被定义为X 和Z 的线性回归与Y 和Z 的线性回归所产生的残差之间的相关性，然后使用最小二乘法计算回归的目标基因和条件基因的系数，并使用单样本t 检验检测基因X 和Y 关系的显著性.SINCERA 使用的是线性相关性方法，因此不能获得基因间的非线性关系.SCENIC 首先使用可获得基因间非线性关系的GENIE3 或GRNBoost2 算法推断转录因子和候选靶基因之间的共表达模块；然后通过CisTarget（https：//resources.aert-slab.org/cis-target/）识别出调控因子结合模体（binding motif）在靶基因中的显著富集模块，并构建调控子；最后通过计算细胞中所有基因表达排序的恢复曲线下的面积（area under curve，AUC），对每个细胞中的每个调控子的活性进行评分，从而生成一个二值化的活性矩阵.SCENIC 最初在R/Bioconductor 中实现，后来被移植到python，命名为pySCENIC[38]，pySCENIC通过软件容器和并行化大幅缩短了程序运算时间，SCENIC/pySCENIC 的主要局限是CisTarget 步骤无法预测具有未知模体的转录因子调控子.

这类算法需要某种度量方法计算每对基因间的关系.除SINCERA 外，其他算法均能获得基因间的非线性关系.当样本量不够大时，NLNET 检测线性关系的能力低于基于线性关系的推断算法，因此NLNET为了获得非线性关系可能会失去一些弱线性关系，但是可以将其与线性方法相结合，使结果更全面.此外，这类算法一般要求数据尽量同构或尽量少的细胞类型，对于异构数据（不同发育阶段或不同细胞类型），算法可能忽略基因间的变化关系.这类算法的一个常见缺点是不能推断出基因间的具体调控关系（如激活或抑制），但是SCENIC 和SINCERA 可以利用现有数据库中的已知调控网络进行富集分析来弥补这一缺点.

1.2.5 基于细胞特异性的算法

传统构建GRN 的算法都是基于所有细胞的相似性而忽略了细胞间的异质性.在单细胞的基础上，细胞特异性网络将数据从“不稳定”的基因表达形式转化为“稳定”的基因关联形式.这类算法的输入是细胞的原始基因表达矩阵，对于每类细胞，输出一个细胞特异性网络，即一对基因可能在某类细胞中有关联，但在其他类细胞中没有关联.这类算法主要包括CSN（cell specific network）[19]、CCSN（conditional cell specific network）[20]和CeSpGRN（cell specific GRN）[21].

CSN 设计了每个细胞中每对基因独立性的统计量，以统计独立性确定该对基因的关联，该统计量能够很好地区分独立关系和非独立关系的基因对.如果一对基因相互独立，则归一化后的统计量遵循标准正态分布.CSN 作为一种无监督算法，直接由基因表达矩阵构建而成，不需要预先知道细胞类型或对细胞聚类.CCSN 是为了解决CSN 对间接效应的过高估计而提出的条件细胞特异性网络算法.与CSN 不同，CCSN构造2 个基因在第3 个基因条件下的条件独立性统计量.CCSN 对每个条件基因都构造该细胞的特异性网络，然后将所有基因的特异性网络整合为最终的细胞特异性网络.虽然CCSN 的计算复杂度高于CSN，但是它不仅可以通过消除基因间的间接关联来测量基因间的直接关联，还可以基于单细胞网络进行细胞聚类和降维.另外，通过定义每个细胞的基因间关联数据的网络流熵（network flow entropy，NFE），CCSN可以估计单个细胞的分化潜能，从而掌握细胞分化的谱系动力学性质.CCSN 的缺点是无法构建因果关系的基因关联.CeSpGRN 使用高斯Copula 图模型对正调控边和负调控边进行推断，而CSN 和CCSN 算法无法推断正负调控关系.CeSpGRN 可以在任何轨迹结构或没有轨迹信息的数据集上工作，不需要预先知道细胞的时间信息.该算法根据基因表达相似性设计高维加权核获得细胞轨迹平滑变化，通过计算每对细胞在K 最近邻图中测地距离的高斯核函数得到这对细胞的核权重，最后通过交替方向乘子法优化模型得到每个细胞的基因邻接矩阵.CeSpGRN 的主要缺点是没有将GRN 建模为有向图.

基于细胞特异性网络推断算法可以识别基因间的相互作用，并在单细胞数据上描绘网络的异质性.这类模型有助于发现新的细胞类型，并能够揭示某些基因在特定细胞类型中的重要作用.

1.2.6 基于深度学习的算法

非深度学习的推断算法大多是无监督的，不能在标记好的数据上进行训练，属于特定假设下的特定方法.这些算法虽然在某些情况下表现较好，但是容易导致更多的假阳性或假阴性结果.而深度学习模型一般不需要对输入数据进行模型或分布假设，因此能够模拟细胞异质性下的基因相互作用并克服噪声和减少错误.这类算法主要包括CNNC（convolutional neural network for co-expression）[22]、DeepDRIM（deep learning-based direct regulatory interaction model）[23]、Deep-SEM[24]、scGeneRAI（single-cell gene regulatory network prediction by explainable AI）[25]、STGRNS[26]和GENELink（gene regulatory network via link prediction）[27].

CNNC 运用卷积神经网络（convolutional neural network，CNN）进行有监督的基因关系推断.该算法首先将所有细胞的每对基因汇总成共现（co-occurrence）直方图，然后将直方图转换为类图像形式（经典概率函数矩阵）以便于进行神经网络训练.CNNC 以少量带标记的基因对集合作为训练集，可以学习区分各类基因关系.CNNC 的输出层可以根据应用需要输出多种类型的关系，如对于因果推断关系，CNNC 可以输出两基因无交互的概率和一个基因调控另一个基因的概率等.DeepDRIM 是构造细胞类型特异性（cell type specific）GRN 的有监督深度神经网络模型，与CNNC不同，DeepDRIM 主要针对转录因子基因对（TF-gene pair）进行预测.DeepDRIM 模型的结构如图4 所示.

图4 DeepDRIM 模型的结构Fig.4 Structure of DeepDRIM

DeepDRIM 首先将转录因子基因对的联合表达二维直方图转化为主图像，并将与该基因对具有强正协方差的基因对图像作为主图像的最近邻，以消除传递交互引起的错误.然后将2 个堆叠卷积层嵌入结构网络A 和B，分别用于处理主图像和最近邻图像.网络A 类似于VGGnet[39]，包含堆叠卷积层和最大池化层.网络B 的结构与A 相似，是A 的类孪生网络，用于处理多个最近邻图像，所有子网络共享权值.网络B 由已知转录因子的基因对相互作用训练，这些相互作用来自公开可用的细胞类型特异性染色质免疫共沉淀测序（ChIP-seq）数据.最后DeepDRIM 将每个图像对应的向量串联成高维向量，经过2 个堆叠全连接层压缩，并使用Sigmoid 函数进行二分类，得到有或无交互2 种结果.虽然DeepDRIM 使用的神经网络模型比CNNC 复杂，但是CNNC 能够得到一个基因调控另一基因的概率，而DeepDRIM 只能得到简单的二分类.

DeepSEM 的算法架构是一个beta 变分自动编码器（beta variational autoencoder，beta-VAE），包括编码器、GRN 层、用于矩阵逆运算的逆GRN 层和解码器.DeepSEM 的编码器每次输入一个基因的表达，不同基因共享beta-VAE 的权重.GRN 层和逆GRN 层均为基因交互矩阵，不同基因之间条件依赖关系的邻接矩阵使用结构方程模型表示，GRN 层用来拟合基因表达矩阵和邻接矩阵之间的非线性关系.由于矩阵逆运算的复杂性，该算法的运行时间会随基因数的增加而增加.scGeneRAI 使用分层相关传播（layer-wise relevance propagation，LRP）估计每个基因的相关性，从静态scRNA-seq 数据中推断GRN，以LRP 得分衡量目标基因与所有预测基因之间相互作用的强度.STGRNS的模型结构包括4 个模块，分别为GEM（gene expression motif）模块、位置编码层、Transformer 编码器和分类层.GEM 模块将基因对转化为可输入Transformer编码器的格式；位置编码层用于捕获位置或时间信息；Transformer 编码器用于计算不同子向量（特别是关键子向量）的相关性；分类层输出最终的分类结果.对于静态数据、伪时序数据或时序数据，STGRNS 都可以根据已知的基因间关系推断GRN.GENELink 是基于图注意力网络（graph attention network，GAT）的有监督模型，该算法首先聚合条件特异性基因表达谱和基于知识的基因相互作用矩阵，用以学习基因的低维向量化表示，然后通过优化嵌入空间，并学习特定的基因特征，用于下游相似性测量或基因调控的因果推理.

在基于深度学习的算法中，CNNC 能够区分相互作用、因果对、负对，或任何可以定义的基因关系；DeepDRIM 只可以识别TF 和基因的相互作用及因果关系；DeepSEM 不仅可以利用GRN 作为“桥梁”，整合不同的单细胞形态，构建共同的潜在空间，还可以集成其他分子相互作用网络，如蛋白和蛋白相互作用网络，开放染色质数据、DNA 的结合模体和遗传相互作用网络；STGRNS 可以处理静态和时序数据，根据已知基因关系训练Transformer 进行GRN 推断，基于Transformer 的方法在适用性和性能方面优于其他深度学习模型的算法[26].基于深度学习的算法相较于传统算法有3 个显著优势：①scRNA-seq 数据隐含着生物实体间和过程间复杂的相互依赖关系，这些实体和过程通常具有内在的噪声，并且发生在多个尺度上，深度学习算法能够检测生物信息数据中的复杂模式；②神经网络可以自然地扩展到包括表观遗传和序列信息在内的互补数据；③深度学习算法适用于高通量的测序数据，随着测序深度和规模的提升以及训练样本的增加，深度学习算法的性能会得到提升.

2 GRN 推断算法综合比较

2.1 各算法的相关信息和应用场景

本文综述的GRN 推断算法都有自身的优势和局限性.表1 给出了这些算法的相关信息，包括算法工具、使用语言和网络类型等，表2 给出了这些算法的主要优势和局限性.在实际应用中，可根据不同的应用场景和具体研究问题选择GRN 推断算法.基于单细胞数据研究细胞发育、细胞分化、疾病发展或疾病疗程，则需要考虑基于时序数据或者基于伪时序轨迹推断的算法，如SCOUP、SCODE、SCENIC、LEAP、PIDC和SINGE 等.对于需要有向网络的研究，可以选择SCOUP、SCODE、LEAP 和SINGE 等算法.非线性关系的基因网络比线性关系更贴合实际的基因调控关系，因此SCOUP 和SINGE 更适合于时序数据的非线性有向网络的预测.对于仅需了解与某些基因相关的关键调控因子，则应选择能够严格界定基因关系的二值化网络，如SCNS、BTR 和DeepDRIM.对于已得到足够多的已知基因调控关系的数据（包括ChIP-seq 或ACTC-seq 等），则可以选择基于深度学习的算法进行训练和预测.

表1 GRN 推断算法的相关信息Tab.1 Relevant information of GRN inference algorithms

表2 GRN 推断算法的应用场景、优势和局限Tab.2 Application scenarios，advantages and limitations of GRN inference algorithms

2.2 算法主要评价指标

GRN 推断算法的主要评价指标包括AUROC（area under the ROC）[41]、AUPRC（area under the precisionrecall curve）[42]、早期精确度（early precision，EP）[31]和早期精确比（early precision ratio，EPR）[31].

ROC（receiver operating characteristic）曲线的横轴和纵轴分别为假阳性率和真阳性率，AUROC 定义为ROC 曲线下面积.相对于准确率和召回率等指标，ROC 能够反映真假阳性在不同阈值下的分数变化，而AUROC 比ROC 能更直观地反映算法的性能，因此大部分GRN 推断算法使用AUROC 代替ROC，AUROC的数值越大，说明算法性能越好.AUPRC 也常用来评价GRN 推断算法性能，该指标定义为精确率-召回率（precision-recall，PR）曲线下面积，同AUROC 一样，AUPRC 的数值越大，说明算法性能越好.EP 为预测网络前k 条边中真阳性的比例，其中k 为基准真实（ground true）网络的边数.EPR 为模型随机预测的早期精度与模型早期精度的比值，其中随机预测的精度为基准真实网络的边密度.

2.3 推断算法比较研究

2018 年，Chen 等[43]使用GeneNetWeaver[44]生成的两组大肠杆菌单细胞仿真数据和2 组分别来自老鼠胚胎干细胞（ESC）和造血干细胞（HSC）样本的真实单细胞数据评估了8 种早期的GRN 推断算法.该研究分别绘制了算法在仿真数据集和真实数据集上的PR 曲线和ROC 曲线，发现8 种算法的性能都不高（SCODE在仿真数据集上的AUROC 值最高，但没超过0.65），分析表明，基于scRNA-seq 数据的GRN 推断算法并不一定比基于批量RNA 测序的推断算法性能更好.

2020 年，Pratapa 等[31]分析了GeneNetWeaver 的缺点，设计了一个基于布尔模型的仿真数据生成工具BoolODE，使用BoolODE 生成了400 多个包含6 个合成网络和4 个精选布尔模型的仿真数据集，然后在生成的仿真数据和5 个实验性人类和小鼠scRNA-seq数据集上评估了12 种非监督、无附加信息和非细胞特异性的GRN 推断算法，构建了基准数据集综合评估框架BEELINE.该研究结果表明，PIDC、GENIE3 和GRNBoost2 在精确度方面领先其他算法；GENIE3 和PIDC 在静态细胞数据集上多次运行可表现出更好的稳定性，而GRNBoost2 对dropout 的灵敏度较低；在给定较好伪时序数据的情况下，SINCERITIES 具有较好的性能和稳定性.同年，Dai 等[45]简单评估了19 种GRN推断算法，发现非线性关系的推断算法通常比线性关系的算法能更准确地预测基因调控关系.

2021 年，Nguyen 等[46]使用GeneNetWeaver 工具生成不同基因规模、不同细胞数量和不同dropout 率的仿真scRNA-seq 数据集，评估了15 种GRN 推断算法，考察其重构参考网络的AUROC 值、不同dropout率和稀疏水平的敏感度以及算法的时间复杂性.该研究发现，大部分算法只能推断小规模基因的网络，如BTR 要求基因数小于30，SCNS 要求基因数小于64，SCINCERA 要求基因数不超过500，而SCODE、NLNET、SCENIC、LEAP 等算法能够分析所有基因规模的数据；SCENIC 在大部分模拟数据上AUROC 的中位值最高；LEAP 和NLNET 是运行最快的算法，即使是3 000 个基因、1 000 个细胞的数据集也可以在数分钟内完成一次分析；在不同dropout 率和稀疏水平下，SCODE、SCOUP 和SCENIC 的AUROC 值始终高于0.5，其中SCOUP 的表现最稳定.

2023 年，Xu 等[26]在48 个涉及不同物种、谱系、网络类型和网络规模的基准数据集上进行实验，发现STGRNS 优于有监督算法CNNC 以及4 种无监督算法PIDC、LEAP、SINGE 和DeepSEM，并且STGRNS 在“TF 基因对”预测上具有一定的可移植性，且有更好的超参数鲁棒性.

本文在人类胚胎干细胞（hESC）数据集上简单评估以上26 种算法，并计算其AUROC 值，结果如图5所示，可以发现DeepDRIM、GENELink 和STGRNS 的表现优于其他算法.

图5 GRN 推断算法的性能评估Fig.5 Performance evaluation of GRN inference algorithms

3 GRN 推断算法面临的机遇与挑战

由于scRNA-seq 数据可能无法为可靠的网络推断提供足够的分辨率，因此，当前的推断算法难以预测真实的GRN 或推测性能较低.一般来说，GRN 推断算法包括3 个关键组成部分：数据预处理、特征提取和潜在基因调控模式的推断.数据预处理的目的是转换数据，使其更易于分析，如对细胞进行伪时序排序和降维等.即使单细胞数据不构成实际的时间序列，也可以根据伪时序索引对其进行排序.特征提取旨在选取特定基因和提取基因调控信息，如基因表达均值和基因互信息等，但是在很多情况下，统计量没有完全刻画基因间的关系，表达模式间的统计关系对应于调控相互作用的假设也可能存在固有缺陷.对潜在基因调控模式的推断是在所有可能的GRN 空间中找到与数据统计最匹配的GRN，但是基因调控的随机性引起的数据内在噪声会降低预测的准确性.

单细胞数据本身的生物异质性和技术异质性也为GRN 推断带来挑战.生物异质性主要由基因随机突变或内部细胞相互作用和环境变化引起的表型变化导致.即使是批次相同的细胞类型，由于噪声和化学批次不同等因素，单个细胞类型的细胞数量和转录状也可能会有显著差异.GRN 主要代表了所有细胞群的平均活动，这些细胞群可以由最丰富的细胞类型所控制.因此，组织网络无法捕捉单个细胞群的独特行为，以及细胞如何相互作用来执行更高层次的组织功能.尽管一些GRN 推断算法声称能够模拟基因表达的随机部分，但是区分不同异质性仍是一个重大挑战.此外，GRN 推断算法的验证和性能评估是最具挑战性的问题，目前仍缺乏高标准的指标评估这些数学模型的性能.针对不同的应用场景，整合其他类型的数据用于训练和验证非常重要，如运用已知的转录因子结合位点信息的ChIP-seq 数据，或者依靠染色质可及性和基因甲基化等多组学测序数据.

随着单细胞测序技术的发展，单细胞数据GRN推断算法的进步将推动细胞功能的研究，在细胞水平上发现缺失的GRN 成分，可应用于识别疾病的生物标志物和信号通路，从而为疾病精准医疗和设计药物提供支持.