从单细胞基因表达数据推断细胞特异性基因调控网络

编译| 林荣鑫 审稿| 程昭龙，王静

本文介绍由佐治亚理工学院计算科学与工程系的Xiuwei Zhang等人的研究成果。基因调控网络（GRN）可以被视为细胞的另一个特征，有助于发现每个细胞的独特性。然而，目前仍然缺少重建细胞特异性GRN的方法。作者提出了一种从单细胞基因表达数据推断细胞特异性GRN的方法（简写为CeSpGRN）。CeSpGRN使用高斯加权核，从发育过程中的细胞以及该细胞上游和下游细胞的基因表达谱中构建给定细胞的GRN。CeSpGRN可用于推断任何轨迹或簇结构的细胞群中的细胞特异性GRN，并且不需要额外输入细胞的时间信息。经实验证明，CeSpGRN在重建每个细胞的GRN以及检测细胞间的相互调节作用方面性能优越。

01简介

基因调控网络（GRN）表示的是基因在生物过程中如何相互调控。从基因表达数据推断GRN是一个具有挑战性的问题。单细胞基因表达数据已被用于推断GRN，其中每个细胞作为一个样本。虽然已经开发了许多方法，但这些方法的目的是从多个单细胞的基因表达数据中构建一个GRN。然而，GRN是高度动态的，并且拓扑结构会随着时间发生变化。有研究者开发了在特定时间点使用微阵列或RNA-seq数据构建时变GRN的方法。然而，由于这些数据集中的时间点数量通常较少，并且每个时间点只能测量一组细胞的批量基因表达数据，因此可能无法检测到某些网络重构活动。

当在动态过程中分析scRNA-seq数据时，每个细胞被认为是动态过程中的一个独特阶段，其基因表达谱代表相应阶段的细胞状态。因此，不同的细胞可以具有不同的GRN。本文作者提出的CeSpGRN可以推断每个单细胞的GRN，并且获得的细胞特异性GRN对研究GRN在动态过程中或跨越空间情况下的变化具有重要意义。

在GRN的推断过程中，同时考虑所有的基因具有较大的挑战。经典的GRN推断方法需要足够数量的样本，这使得为每个细胞推断出细胞特异性GRN似乎成为不可能。作者受到现有方法KELLER的启发，该方法使用加权内核来推断时变GRN，该内核允许从附近的时间点“借用”信息。在CeSpGRN中，作者假设经历发育或分化过程的细胞的GRN沿细胞轨迹平稳变化，该假设通过加权内核实现了跨不同细胞的“信息共享”。但是，CeSpGRN和KELLER在方法上存在实质性差异：（1）KELLER使用马尔可夫随机场（MRF）对GRN进行建模，该GRN采用二进制基因表达数据。为了避免二值化过程中发生的信息丢失，从而更好地适应单细胞基因表达数据的分布，作者在CeSpGRN中使用高斯Copula图模型（GCGM）；（2）KELLER需要每个样本的时间信息，并使用每个样本的一维时间信息构建内核。而CeSpGRN不需要单个细胞的已知时间信息，并且使用细胞的高维基因表达谱构建内核；（3）由于KELLER使用细胞的时间信息设计内核，因此无法区分处于相同发育时间但类型不同的细胞。相比之下，CeSpGRN中使用的高维核可以处理任何轨迹结构的细胞数据集，也可以处理不同簇的细胞数据集。

高斯图模型（GGM）被用来表示生物网络或其他网络。TREEGL中使用GGM来构建细胞分化过程中的多个GRN，并通过总变差（TV）惩罚来约束GRN，其目的是使构建的GRN之间的差异最小化。与KELLER和CeSpGRN中采用的加权核方法相比，TV惩罚框架的计算成本更高。在给定大量单细胞的情况下推断细胞特异性GRN时，这种时间复杂度的差异会更加明显。此外，TREEGL是以树作为输入，在使用scRNA-seq数据的情况下，需要推断表示细胞轨迹的树。相比之下，CeSpGRN可以处理细胞形成的任何图结构轨迹的数据集，并且不需要图或树信息作为输入，因为图结构可以在高维基因表达数据构建的核函数中表达。

最近，GCGM已被用于从单细胞基因表达数据，尤其是scRNA-seq数据中构建GRN。与GGM相比，GCGM解释了单细胞基因表达数据的非高斯性。然而，现有使用GCGM的方法不能推断细胞特异性网络，但可以推断出一组细胞的整个GRN。

通过使用GCGM模型，CeSpGRN不局限于二进制或高斯分布的基因表达数据；通过使用高维加权核，CeSpGRN可以为数据集中的每个细胞推断出一个GRN，其中细胞可以形成任何轨迹或簇结构。作者还将CeSpGRN扩展到一个名为CeSpGRN-TF的版本中，该版本可以利用部分基因是转录因子（TF）的先验信息。作者对来自两个不同模拟器的模拟数据集和两个真实数据集进行实验，并将它们与CSN和最先进的GRN推理方法进行比较。

02结果

作者在来自两个不同模拟器的模拟数据集和两个真实数据集上测试CeSpGRN的性能。首先基于模拟器为特定细胞生成GRN，进而为每个给定GRN的细胞生成基因表达数据。对于后一步，第一个模拟器生成多变量高斯（基因表达）数据，第二个模拟器是从BoolODE修改而来，使用带有非线性希尔函数的微分方程生成基因表达数据。第一个模拟器上的测试用于验证该方法的性能，第二个模拟器上的测试用于证明该方法在实际单细胞基因表达数据集上的适用性和准确性。在测试中，作者考虑两种类型的细胞轨迹，一种是“线性”，细胞形成线性轨迹；另一种是“分叉”，细胞分化成两种不同的细胞类型。

作者将CeSpGRN与基线方法CSN、GENIE3、SCODE相比较。GENIE3和SCODE是在scRNA-seq数据上表现最好的GRN推理方法，GENIE3是一种基于回归树的方法，而SCODE是一种基于常微分方程的方法。由于GENIE3和SCODE是从所有细胞中推断出静态GRN，为了获得细胞特异性GRN与CeSpGRN进行比较，作者设计了它们的“动态”版本GENIE3-Dyn 和 SCODE-Dyn。动态版本的工作方式是沿着细胞轨迹，将细胞分成100个细胞片段，然后在每个片段上运行GENIE3或SCODE，推断出的细胞特异性GRN就是相应片段的GRN。

作者在同一细胞中计算真实GRN和预测GRN邻接矩阵间的Pearson相关性和Spearman相关性，用以评估推断的细胞特异性GRN的质量。并且为了评估这些方法检测不同细胞GRN之间差异（即变化的边）的能力，作者还使用了四种相同的测量方法：Pearson相关性、Spearman相关性、AUPRC和早期精度分数。

2.1 多元高斯数据检验

作者用一系列超参数运行CeSpGRN，并使用PCA可视化一个模拟数据集和CeSpGRN推断出的相应GRN（图1a），结果表明推断的GRN沿着与基因表达数据一致的连续分叉轨迹演化。

图1 GRN推理方法的比较

由于GRN中存在正负权重的边，因此在计算AUPRC和早期精度分数时，作者区分了正边和负边：先计算正边的AUPRC和早期精度分数，然后再计算负边，最后计算正边和负边分数的平均值，并分别称为“AUPRC（有符号）”和“早期精度（有符号）”分数。因为CSN输出的二进制网络不会产生精确召回曲线，所以CSN的AUPRC分数不可用。GENIE3和CSN无法区分权重为正的边和权重为负的边，它们只返回非负值的边。在计算GENIE3和CSN的评估分数时，作者首先获取真实网络中边的绝对值，然后进行比较。实验结果表明，在不同的参数和评估指标下，与基线方法相比，CeSpGRN推断出的细胞特异性GRN更准确（图1b、c）。

为进一步评估不同方法检测不同细胞GRN之间差异的能力（图1d、e）。作者通过测量两个推断图之间的变化，并将它们与两个真实图中相应的变化边进行比较。在图2中，作者可视化了模拟时间为0.50的细胞i到模拟时间为0.54的细胞j的边的变化。结果发现CeSpGRN可以正确检测到大部分变化的边。作者还研究了超参数对CeSpGRN性能的影响，并为参数选择提供了指导。

图2 从细胞i到细胞j边的变化

总体来说，CeSpGRN比基线方法性能更加优越。这可能得益于连续变化的建模，也可能得益于生成数据的高斯分布。接下来，作者将用非高斯数据集测试CeSpGRN。

2.2 BoolODE模型模拟数据测试

作者使用模拟器生成了5个具有线性轨迹的数据集和5个具有分叉轨迹的数据集，其中包含不同的随机种子。每个数据集有1000个细胞和20个基因。

图3 BoolODE模型模拟数据的测试实验

为评估推断细胞特异性GRN的准确性，作者在模拟数据集上测试了CeSpGRN和基线方法。实验结果表明，从CeSpGRN推断的GRN中可以清楚地检测到预期的轨迹（图3a、b）。与其他方法相比，CeSpGRN推断的细胞特异性GRN准确性更高（图3c）。作者还通过合并TF信息来测试CeSpGRN-TF的性能，图中可以看出合并TF信息进一步提高了CeSpGRN的性能（图3d）。

2.3 人类THP-1细胞数据集测试

作者通过CeSpGRN来获取该数据集的细胞特异性GRN。虽然真实数据的细胞特异性GRN无法获得，但THP-1细胞的45个 TF存在一个实验测量的GRN，作者将之称为扰动矩阵。由于该GRN不是细胞特异性GRN，而是对应于细胞群的GRN，因此作者从CeSpGRN结果中获得了群体水平的GRN，用以与扰动矩阵进行比较。作者还总结了CeSpGRN推断的所有细胞的GRN的邻接矩阵，以获得一个网络。通过与GENIE3、SCODE和CSN推断出的群体水平网络进行比较，发现CeSpGRN推断GRN的准确度得分更高（图4）。

图4推断的THP-1单细胞群体水平GRN的准确性

为了分析与变化边相关的基因，作者选取一组随着时间推移获得相互作用的基因，以及一组随着时间推移失去相互作用的基因，获得相互作用的一组基因是：BCL6、ETS1、PPARD、PPARG、SNAI1、SNAI3、SP3、SPIB、TCF3，失去相互作用的基因是：FOSB、MYEF2、TCFL5、NFATC2、NFATC1、UHRF1、JUN。基于群体水平上通过实验测量的定向GRN中的基因位置，作者将45个基因分为上游、中游和下游基因。实验结果可观察到获得相互作用的基因大多是中下游基因，而失去相互作用的基因大多是中上游基因（图5）。这表明，在对应于早期时间点的细胞特异性GRN中，群体水平GRN中上游基因的相互作用比下游基因的相互作用更活跃，而在晚期时间点的GRN中，下游基因的相互作用更活跃。

图5 THP-1数据集中获得相互作用和失去相互作用的基因表

2.4 小鼠胚胎干细胞数据集测试

作者进一步将CeSpGRN应用于小鼠胚胎干细胞（mESC）数据集。使用CeSpGRN-TF推断细胞特异性GRN，并将CeSpGRN-TF和CSN推断的细胞特异性网络进行比较。由于AUPRC比率不能用于CSN，因此只能比较两种方法的早期精度比率（图6）。又因为CeSpGRN-TF推断的细胞特异性GRN能够分析数据集中GRN的动态过程，作者为此计算了每个推断的GRN的密度（图7）。结果表明，在分化过程开始时，调节活动较多，随着分化的进行，需要的调节活动也会减少。进一步实验发现，关键的调节活动不需要在整个分化过程中都保持活跃，并且在某些阶段，它们的强度可能较低。

图6 CeSpGRN-TF和CSN的早期精度比曲线

图7 对GRN进行UMAP可视化

03总结

作者提出的CeSpGRN，它从单细胞基因表达数据中推断加权、有符号和无向细胞特异性GRN。CeSpGRN可以应用到任何轨迹结构和簇结构的数据集中，因为其设计的核函数是基于细胞之间的成对流形距离，这使得CeSpGRN比需要时间序列数据的方法更适用。高斯Copula图模型的使用使CeSpGRN能够处理具有分布特征的数据，该分布可以通过单变量单调函数转换为多变量高斯分布。作为第一种推断细胞特异性GRN的方法，CeSpGRN可以在以下方向进一步改进。首先，如果存在在细胞间不会发生变化的相互调节作用的先验信息，可以将这些信息结合起来，以提高推断GRN和检测变化边的准确性。其次，CeSpGRN将GRN建模为无向图，未来的工作可以将GRN建模为有向图。此外，随着单细胞多组学技术的发展，将多模态数据建模到网络推理框架中也是一个很有前途的发展方向。

参考资料

Ziqi, Z., Jongseok, H., Le, S. et al. Inferring cell-specific gene regulatory networks from single cell gene expression data. bioRxiv (2022).

https://doi.org/10.1101/2022.03.03.482887

数据链接：

https://github.com/PeterZZQ/CeSpGRN/tree/master/data

代码链接：

https://github.com/PeterZZQ/CeSpGRN