通过从单细胞转录组数据中推断基因调控网络来揭示基因扰动效应

《PLOS Computational Biology》：Unveiling gene perturbation effects through gene regulatory networks inference from single-cell transcriptomic data

　　摘要：生理和病理过程受到称为基因调控网络（GRNs）的基因网络的控制。通过重构GRNs，我们可以准确模拟细胞在自然状态下的行为，并预测基因变化将如何影响它们。现在，已经获得了多种物种中多种细胞过程的单细胞转录组数据。因此，一种无需任何额外先验知识即可从单细胞RNA测序（scRNA-seq）数据构建预测性GRNs的方法，可能会对我们的生物学过程及其关键基因的理解产生重大影响。为此，我们开发了IGNITE（利用逆动力学理论和实验推断基因网络），这是一种无监督的机器学习框架，旨在直接从未受干扰的单细胞RNA测序数据中推断出有向、加权和有符号的GRNs。IGNITE利用这些GRNs来生成在单个或多个基因扰动下的基因表达数据。IGNITE基于动力学Ising模型的逆问题，这是一种已经成功应用于生物学网络的统计物理模型。我们在两个互补的多能干细胞（PSCs）系统上测试了IGNITE：从小鼠PSCs（从幼稚状态转变为形成状态）和人类PSCs（分化为定型内胚层）。这些数据集在物种、发育轨迹和单细胞技术（10X与Fluidigm C1）方面存在差异，从而提供了对通用性的严格测试。仅使用未受干扰的scRNA-seq数据，IGNITE模拟了单个和多个基因敲除（KOs），并产生了与独立实验观察结果一致的预测。在小鼠PSCs中，IGNITE生成的野生型数据与实验高度相关，并准确预测了Rbpj、Etv5和三重KOs的效果；而在人类PSCs中，它正确预测了促进分化和阻止分化的扰动，这与已发表的研究结果一致。这些结果表明，IGNITE能够跨物种和技术稳健地捕捉基因交互逻辑，从而可以直接从scRNA-seq数据中进行稳健的计算机模拟扰动分析。

作者总结：生理和病理过程依赖于基因之间的复杂相互作用，这些基因组织成基因调控网络，控制细胞的行为和对变化的响应。理解这些网络对于预测基因扰动如何影响细胞状态至关重要，但实验测试基因功能通常耗时、成本高昂或技术上受到限制。计算方法可以帮助应对这一挑战，尽管许多现有方法需要额外的实验数据或先验生物学知识。在这里，我们介绍了IGNITE，这是一种无监督的机器学习框架，它可以直接从未受干扰的单细胞基因表达数据中推断出基因调控网络。该方法无需依赖先验信息，即可重建有效的基因相互作用，并利用它们来模拟细胞对基因扰动（如基因敲除）的反应。我们将这种方法应用于两个不同的生物学系统：从小鼠幼稚多能干细胞转变为形成状态的系统，以及人类多能干细胞向内胚层分化的系统。尽管在物种、发育轨迹和实验技术方面存在差异，该模型仍准确地再现了观察到的细胞状态，并预测了与独立实验结果一致的基因扰动效应。总体而言，这项工作介绍了一个通用且可扩展的机器学习框架，用于直接从单细胞数据构建基因调控的预测模型。

1 引言：细胞过程的正确运作需要多个组分的受调控相互作用。通过单独研究每个组成部分，我们无法理解参与特定过程的各个组分的组织如何导致复杂的行为。然而，研究整个系统仍然是一项具有挑战性的任务。单细胞RNA测序（scRNA-seq）数据[1,2]具有解开不同生物学过程的潜力[3]，并将其作为复杂系统进行建模可以使用统计物理描述来提供对其潜在机制的更深入的见解[4–6]。在细胞的所有基本组成部分中，基因是通过相互作用和组织成基因调控网络（GRN）来调节生物学过程的关键要素。为了研究基因之间的相互作用以及它们如何塑造细胞，可以通过实验干扰个别基因的调控，强制它们的激活或失活。一种通过失活特定基因来干扰GRN的实验技术是基因敲除。该技术允许测量单个或多个基因失活对其他基因调控的影响，以及这如何影响细胞的表观动态。然而，由于可以关注的基因数量众多，实验方法将过于耗时且成本高昂。计算机模型代表了一个可行的替代方案，因为它们可以推断复杂的基因相互作用并预测基因扰动的结果，最终揭示关键的调控机制。新兴的计算策略旨在预测扰动如何影响细胞表型。为了有效应用这些方法，理解它们的能力和限制是至关重要的。例如，依赖于实验扰动获得的数据进行训练步骤限制了这些模型的适用性[7]。此外，一些方法需要将多组学数据作为GRN推断的先验知识[8,9]。这意味着需要对同一生物系统进行不同的测量以获取这些数据。或者，一些方法仅可以从未受干扰的单细胞RNA测序数据中推断GRN相互作用[10]。然而，它们仅关注某些特定方面，例如GRN推断、生成野生型（WT）数据或扰动模拟。在GRN推断方法中，RE:IN[11–14]能够使用布尔网络建模方法推断众多基因之间的相互作用，并模拟一个或多个基因的敲除。然而，RE:IN依赖于敲除实验和大量转录组数据，缺乏单细胞分辨率，并且需要包括其他类型的数据（ChIP测序数据，称为ChIP-seq，以及文献整理的相互作用）来约束模型。此外，RE:IN需要探索的可能网络数量随着基因数量的增加而呈指数级增长，使得推断任务在计算上非常耗时。另一种推断GRN的方法是MaxEnt[15]。它利用最大熵原理从基因表达模式中推断基因相互作用网络。通过在满足数据给出的特定约束的同时最大化熵来获得GRN[16,17]。然而，MaxEnt只能推断无向相互作用，且不适合生成扰动数据。SCODE[18]是用于从scRNA-seq数据开始推断有向加权网络的金标准模型之一，它使用一组耦合的线性常微分方程来描述基因表达动态。此外，它可以生成新数据，但尚未验证其用于模拟基因扰动的能力。CellOracle[8]是一种基于网络的机器学习框架。它使用来自不同细胞簇的scRNA-seq数据，并将活跃的调控连接映射到使用scATAC-sec数据（或ATAC-seq数据）计算的GRN框架上。它为数据集的每个簇生成GRN，并模拟具有已知结合基序的转录因子的扰动效应。CellOracle关注GRN扰动在修改细胞表型中的作用，但没有描述GRN如何生成未受扰动的表型。因此，所有这些方法仅研究了一些类似于真实系统的特定特征，我们在表1中进行了总结。然而，开发只需要一种类型的未受干扰数据（例如单细胞转录组学）的可解释模型是至关重要的，这简化了算法的使用，并提供了包括GRN推断和生成规则的描述，解释了GRN如何产生细胞表型，以及模拟基因扰动效应的能力。

下载：PNG（更大图像）| TIFF（原始图像）

表1. 不同计算方法用于GRN推断的比较，突出显示了它们的输入数据类型、对先验信息的需求、生成新数据的能力、实现多个基因敲除（KO）模拟的能力以及推断有向相互作用的能力。星号（*）表示该方法理论上支持该特征，但在参考出版物中未实现和讨论。

为了应对这些建模挑战，我们提出了IGNITE（利用逆动力学理论和实验推断基因网络）。这是一种基于统计物理概念的无监督计算方法，它从未受干扰的scRNA-seq数据开始，推断出能够模拟未受干扰或受扰动细胞表型的GRN模型。特别是，它依赖于动力学Ising模型的逆问题[19]，遵循Zeng等人在[20]中提出的方法。Ising模型最初是为描述铁磁性而开发的，通过使用逆方法，已经应用于包括生物网络和神经网络在内的多个多体系统的建模[15,21,22]。在这里，我们利用其多功能性开发了一个计算模型，利用表1中报告的所有特征来描述GRN。我们将IGNITE应用于研究多能干细胞（PSCs）的发育进展。多能性是指产生成人中发现的所有细胞的潜力，这在植入前胚胎的幼稚外胚层中被观察到。胚胎植入后，外胚层过渡到一个不同的多能性阶段，称为形成阶段，其特征是不同的基因表达谱、表观遗传状态和发育特征。在这个过渡过程中，幼稚多能性标记被下调，表观遗传修饰因子和形成标记被激活，细胞变得能够形成胚层并指定原始生殖细胞（E6.5–7.5）[14]。早期胚胎的多能细胞可以在体外捕获为幼稚胚胎干细胞（ESC），它们保持多能性并具有从幼稚状态向形成状态发展的能力。已经使用了几种遗传方法来识别对维持幼稚状态或向形成状态过渡至关重要的基因。然而，对这种过渡如何被调控的全面理解仍然缺失。以前的方法主要集中在表征幼稚多能性的状态以及阐明控制和维持它的分子组分之间的功能相互作用[23–26]。此外，还有努力致力于在不同细胞培养条件下建模这些核心基因及其动态。例如，Walker等人[23]通过实验推导出小鼠ESC 在多能状态下的可能GRNs，而其他方法依赖于已建立的GRNs，使用微分方程表征了少数相互作用基因的动态[24–26]。然而，尽管这些方法具有信息性和精确性，但它们仅限于描述有限数量基因的相互作用。布尔网络方法不仅应用于批量RNA-seq数据（如RE:IN所做的那样），还应用于单细胞数据[27–30]。这些方法可以获得单细胞分辨率数据，但仍然需要整合不同类型的数据（例如ChIP-Seq和敲除实验），或者需要探索随着基因数量增加而呈指数级增长的可能网络数量，使得推断任务在计算上非常耗时。其他计算策略，包括基于树的模型、相关性和部分相关性分析、信息论和微分方程，已被用于从未受干扰的单细胞多组学数据中建模分化[10]。然而，即使在关于分化的研究中，这些方法也主要侧重于推断GRNs，如前文所述，而没有阐明这些GRNs如何生成分化细胞或检查扰动在分化中的作用。其他方法使用单细胞多组学数据来推断GRN并检查可能的扰动效应，但没有理论描述可以生成推断GRNs的模型[8,9]。因此，研究这个系统以获得一个可靠的模型是至关重要的，该模型能够生成所有观察到的细胞表型，并能够模拟敲除以研究基因在过渡过程中的功能。我们将IGNITE应用于两个互补的多能干细胞系统。首先，使用一个新颖的10X scRNA-seq数据集，该数据集包含从小鼠ESC从幼稚状态转变为形成状态的过程，IGNITE推断出了GRN。从GRN出发，IGNITE生成了新的合成数据，而不依赖于生物学先验知识，例如不同的数据类型（例如ATACseq数据）或基于文献的信息（例如已建立的相互作用）。值得注意的是，生成的野生型数据与实验测量结果在统计学上相当。IGNITE还模拟了三个基因Rbpj、Etv5和Tcf7l1的敲除，无论是单独还是同时进行。这些模拟在很大程度上得到了实验数据的支持。值得注意的是，Tcf7l1是一个在蛋白质稳定性而非转录水平上受到调控的转录因子；IGNITE和其他方法对其失活的模拟并不准确，揭示了仅使用转录组学方法的局限性。其次，为了进一步评估IGNITE的通用性，我们将其应用于Chu等人[31]提供的独立scRNA-seq数据集，该数据集描述了向定型内胚层分化的人类PSCs。这个数据集之前已经被其他团队用于生成规则基因网络（GRN）的推断[10,18]，它是使用Fluidigm C1平台生成的，该平台与我们在小鼠数据中使用的10X技术不同，并且捕获的细胞数量较少。除了代表一个独特的发育轨迹外，这个数据集还为跨物种、细胞数量和实验技术的IGNITE测试提供了严格的基准。IGNITE推断出了一个生物学上一致的GRN，再现了野生型细胞的分化动态，并正确预测了促进和阻断分化的扰动，这与已发表的研究结果一致。总之，IGNITE代表了一种基于模型的计算方法，用于生成新的单细胞表型，作为基因激活/失活的模式，并通过推断所研究的生物系统的GRN来模拟扰动实验，其性能优于当前的黄金标准方法[8,18]。

**2 结果**

**2.1 IGNITE算法**
IGNITE算法包括三个主要步骤：预处理、建模和数据生成（图1A）。
（i）预处理：scRNA-seq输入数据经过质量控制和对数标准化处理，随后进行伪时间（PST）计算，并应用Mini-Bulk（MB）方法（详见方法部分）。PST计算需要数据降维和聚类步骤。PST算法Slingshot[32]将单个细胞沿着伪时间轨迹排序，为每个细胞分配一个代表其在该轨迹上位置的伪时间值[33]。

**2.2 通过scRNA-seq解码多能干细胞动态**
我们分析了小鼠（mPSC）和人类多能干细胞（hPSC），以表征细胞命运转换背后的调控动态。在本节中，我们重点关注mPSC，稍后我们将分析扩展到hPSC，以测试IGNITE在跨物种和实验平台上的性能。我们的mPSC数据集捕捉了在2iL withdrew后从幼稚状态到形成状态的转换过程[35]。这种化学定义的培养条件确保了干细胞的幼稚多能状态[36,37]。在2iL withdrawn后，我们通过scRNA-seq在多个时间点分析基因表达，以监测分化过程中的基因活动动态（图1B）。我们重点研究了文献中已知在PSC分化决策过程中起驱动作用的24个基因[11,14,38,39]，分析了mPSC发育的前48小时。这个时间窗口允许我们研究指导从幼稚细胞到形成细胞转换的相互作用，共分析了9894个细胞。基因表达谱在伪时间上的分布显示了四种时间模式（图1C）：幼稚状态下的基因活性较高，随后下降；形成早期基因在2iL withdraw后不久被激活；形成晚期基因在后期被激活；以及一个没有单调趋势的混合组，如先前报道的[14]。为了计算PST，我们首先对选定的基因集进行了主成分分析（PCA）（见S1A和S1B图，详见方法部分），然后使用t分布随机邻域嵌入（t-SNE）[40]和均匀流形近似与投影（UMAP）[41]（见S1C图）对数据进行投影，并对数据进行了聚类（见S1D图，详见方法部分）。为了验证所得到的簇的生物学意义，我们调查了每个簇内基因的表达分布（参见S1E图中的三个示例）。此外，由于簇主要由2Cell-like细胞（2CLCs）组成（通过全面的差异基因表达分析确定[42]），我们将簇1从后续分析中排除。我们发现上调的标记基因包括Eif1a样基因和Zscan4。由于这个簇偏离了我们分析的主要范围，即从幼稚状态到形成状态的转换，因此将其从进一步的研究中省略。在计算PST之后，我们检查了每个时间点的PST值分布（见S1F图），观察到采样时间与PST排序之间的一致性。即使在PST排序之后，基因的预期时间表达模式也得到了保持。在低维UMAP空间（图1D）中，可以观察到分配给每个细胞的PST值。PST值较低的细胞表现出较高的幼稚基因表达（例如图1E中的Klf4），而在伪时间值较高的细胞中，形成基因表现出增加的活性（例如图1E中的Otx2和Utf1）。经过PST和MB处理后，基因表达动态得到了保留。实际上，其随时间的变化（图1F和S1G）与图1C中描述的结果相似。最后，我们计算了基因活性（GA），仍然观察到了从幼稚状态到形成状态的预期转换（见S1H图）。

**2.3 小鼠GRN推断和野生型数据生成：IGNITE和SCODE**
IGNITE仅使用scRNA-seq数据作为输入来推断一组250个GRN（图1A）。它利用生成数据的相关矩阵及其时间导数来推断最准确和有效的GRN。通过比较输入数据（经过PST和MB处理的scRNA-seq数据）与生成数据集的相关矩阵，并通过欧几里得距离来量化它们的差异来实现这一点。我们将这个量定义为相关矩阵距离（CMD，详见方法部分）。选择了CMD最低（0.43）的GRN（图2A–2C），该GRN最好地再现了初始数据的相关性。通过双样本t检验进行统计验证，得到了t统计量和p值，确认了IGNITE模型生成的CMD与空模型CMD之间存在显著差异。

**2.4 使用scRNA-seq解码多能干细胞动态**
我们分析了小鼠（mPSC）和人类多能干细胞（hPSC），以表征细胞命运转换背后的调控动态。在本节中，我们主要关注mPSC，稍后我们将分析扩展到hPSC，以测试IGNITE在跨物种和实验平台上的性能。我们的mPSC数据集捕捉了在2iL withdraw后从幼稚状态到形成状态的转换[35]。这种化学定义的培养条件确保了干细胞的幼稚多能状态[36,37]。在2iL withdraw后，我们通过scRNA-seq在多个时间点分析基因表达，以监测分化过程中的基因活动动态（图1B）。我们重点研究了24个在文献中已知在PSC分化决策过程中起驱动作用的基因[11,14,38,39]，分析了mPSC发育的最初48小时。这个时间窗口允许我们研究指导从幼稚细胞到形成细胞转换的相互作用，共分析了9894个细胞。基因表达谱在伪时间上的分布显示了四种时间模式（图1C）：幼稚状态下的基因活性较高，随后下降；形成早期基因在2iL withdraw后不久被激活；形成晚期基因在后期被激活；以及一个没有单调趋势的混合组，如先前报告的[14]。为了计算PST，我们首先对选定的基因集进行了主成分分析（PCA）（见方法部分），然后使用t分布随机邻域嵌入（t-SNE）[40]和均匀流形近似与投影（UMAP）[41]对数据进行了投影（见S1C图），并对数据进行了聚类（见S1D图，详见方法部分）。为了验证所得簇的生物学意义，我们调查了每个簇内基因的表达分布（参见S1E图中的三个示例）。此外，由于簇1主要由2Cell-like细胞（2CLCs）组成，我们将其从后续分析中排除，这是通过全面的差异基因表达分析确定的[42]。我们发现上调的标记基因包括Eif1a样基因和Zscan4。由于该簇偏离了我们分析的主要范围，即从幼稚状态到形成状态的转换，因此将其从进一步的研究中省略。在计算PST之后，我们检查了每个时间点的PST值分布（见S1F图），观察到采样时间与PST排序之间的一致性。我们验证了即使在PST排序之后，基因的预期时间表达模式也得到了保持。在低维UMAP空间（图1D）中，可以观察到分配给每个细胞的PST值。PST值较低的细胞表现出较高的幼稚基因表达（例如图1E中的Klf4），而形成基因在伪时间值较高的细胞中表现出增强的活性（例如图1E中的Otx2和Utf1）。经过PST和MB处理后，基因表达动态得到了保留。事实上，其随时间的变化（图1F和S1G）与图1C中描述的情况类似。最后，我们计算了基因活性（GA），仍然观察到了从幼稚状态到形成状态的预期转换（见S1H图）。

**2.5 小鼠GRN推断和WT数据生成：IGNITE和SCODE**
IGNITE仅使用scRNA-seq数据作为输入来推断250个GRN（图1A）。它利用生成数据的相关矩阵及其时间导数来推断最准确和有效的GRN。通过比较输入数据（经过PST和MB处理的scRNA-seq数据）与生成数据集的相关矩阵，并通过欧几里得距离来量化它们的差异来实现这一点。我们将这个量定义为相关矩阵距离（CMD，详见方法部分）。选择了CMD最低（0.43）的GRN（图2A–2C），该GRN最好地再现了初始数据的相关性。统计验证通过双样本t检验得到了t统计量和p值，确认了IGNITE模型生成的CMD与空模型CMD之间存在显著差异。对于Rbpj、Etv5和Tcf7l1的三重敲除实验，模拟数据和实验数据均显示天真基因（na?ve genes）显著增加，而形成基因（formative genes）显著减少。为了量化这些变化，我们计算了敲除（KO）组和野生型（WT）组之间的基因活性/表达差异，并将其缩放到最大绝对值（缩放后的KO-WT差异，见图3A和3B）。然后，我们将生成的数据的缩放后KO-WT差异与从文献[38,39]中获得的实验敲除数据中的差异进行了比较。实验中的缩放后KO-WT差异是通过缩放单个敲除数据[39]和三重敲除数据[38]的对数2倍荧光强度变化（log2FC）得到的。利用这些值，我们使用Spearman相关性（Spearman等级相关测试）和一致性分数（FoA，见方法部分）来量化预测结果与实验结果之间的一致性。IGNITE模拟结果显示，在Rbpj（，）、Etv5（，）和三重敲除（，）条件下，预测结果与实验结果之间存在显著相关性，而Tcf7l1则没有相关性（，），这与该因子的已知转录后调控机制一致[43]（见图3A和3B，表2）。下载：PNG更大图像；TIFF原始图像。

表2：每种方法预测的与实验测量的平均基因活性变化之间的Spearman相关性（Spearman等级相关测试）。括号中的数值表示以科学记数法表示的关联p值。统计上显著的相关性（p < 0.05）以粗体显示。

除了以上内容外，我们还使用一致性分数（FoA）来量化预测基因表达变化与实验结果之间的一致性，FoA定义为预测变化和实验变化方向相同基因的比例（见表3和方法部分）。对于Rbpj KO，FoA为0.65；对于Etv5 KO，FoA为0.74，这两个值与空模型相比都具有统计学意义（二项式检验；见方法部分）。相比之下，Tcf7l1 KO未观察到显著的一致性。

此外，据报道，在实验条件下，所考虑的三重敲除能够扩增细胞，同时大多数细胞仍保持天真标志物的表达[38]。由于IGNITE能够模拟野生型分化细胞的细胞组成（见图2E），我们研究了它是否也能模拟三重敲除细胞的组成（见图3C）。大多数细胞（79%）表现出所有天真基因的激活，剩余21%的细胞表现出一些形成基因（Utf1、Etv4和Tcf15）的活性，并且只有两个天真基因（Klf4、Stat3）被失活。将IGNITE生成的三重敲除基因活性的PCA图与输入的MB基因活性数据的PCA图进行比较（见图3D），发现三重敲除模拟结果位于输入数据早期时间点的区域，该区域以天真状态细胞为主。这些观察结果与先前关于三重敲除的实验发现一致[38]。为了评估IGNITE的结果，我们还使用SCODE和CellOracle模拟了单个和三重基因扰动[8]。我们计算了SCODE和CellOracle模拟的缩放后KO-WT差异，并将其与实验测量的log2FC进行了比较（见图3A和3B）。在所有测试的敲除条件下，IGNITE得到的Spearman相关系数和FoA值都高于SCODE和CellOracle：对于SCODE，Spearman相关系数介于[值范围]到[值范围]之间，而对于CellOracle，相关系数介于[值范围]到[值范围]之间，且FoA值始终低于IGNITE得到的结果（见表3）。此外，SCODE生成的敲除数据并未表现出天真基因的上调，天真细胞的比例仅从WT条件下的21%（S2D图）增加到敲除模拟条件下的28%（S2F图）。对于CellOracle，WT模拟产生的天真细胞比例较低（18%），其特征是Klf4、Klf2和Nanog等标记物的高表达，这与分化中的野生型细胞的预期组成一致（S2G图）。然而，在三重敲除模拟中，CellOracle产生的天真样细胞比例较高（48%），表明方向发生了改变。不过，这一比例仍低于实验报告的结果，其中大多数三重敲除细胞保持天真状态（S2H图；[38]）。对于IGNITE目前获得的结果，我们选择了CMD值最低的有效基因调控网络（GRN）（见图1A）。作为进一步分析，我们验证了IGNITE在模拟生物系统扰动方面的性能是否对于具有低标准化距离的任何GRN都保持稳定。为此，我们选择了CMD值最低的10个GRN，并计算了模拟结果与实验结果之间的Spearman相关性（见图3E）。这些结果表明，IGNITE的预测结果对不同表现最佳的GRN模型之间的变化具有鲁棒性。

接下来，我们分析了IGNITE推断出的GRN是否包含了之前使用独立技术验证过的基因相互作用。这些GRN被可视化为相互作用矩阵，其中基因j（列）与基因i（行）之间的相互作用由矩阵中的对应条目（i,j）表示。值得注意的是，正如先前的报道，天真组和形成组中普遍存在正相互作用。而在天真基因和形成基因之间观察到了大量的负相互作用（见图4A）。

为了进行比较，我们使用三种额外的方法SCODE、CellOracle和MaxEnt也推断出了GRN。MaxEnt的GRN相互作用矩阵显示天真基因之间存在正相互作用，某些形成基因之间也存在正相互作用（见S3A图）。然而，在IGNITE的GRN中观察到的天真基因和形成基因之间的负相互作用模式并不明显。我们还使用SCODE和CellOracle计算了相互作用矩阵（见S3B和S3C图）。SCODE的矩阵没有显示出明显的相互作用模式，因为没有明显的正负相互作用块。CellOracle推断出了五个GRN，每个输入的scRNA-seq数据集中对应一个GRN。这些CellOracle GRN的相互作用模式与IGNITE相似，天真基因和形成基因之间都存在正相互作用。只是不同簇之间的相互作用强度发生了变化，但模式保持不变。为了定量比较这些推断的GRN，我们确定了文献中18个经过实验验证的相互作用[38,44–49]。我们认为，当某个因子被证明可以直接与目标基因的启动子相互作用，并且其遗传失活或过度激活导致其直接目标的水平发生显著变化时，这种相互作用就被认为是经过实验验证的。在图4C中，我们展示了仅包含已知相互作用的IGNITE子网络。其链接代表了正确推断出的已知相互作用。重要的是，对空模型进行的统计测试（排列测试，见方法部分）证实这些相互作用不可能仅仅是随机发生的。从这些已知的相互作用中，我们计算了正确推断的相互作用比例（FCI），即与已知相互作用符号相同的相互作用比例（见表4）。由于MaxEnt GRN的无向性质，我们没有考虑已知相互作用的方向性。所有方法的FCI值是可比的，表明它们在识别相互作用方面具有相似的能力。

除了FCI之外，我们还计算了IGNITE推断的GRN中的相互作用强度与其他方法推断的相互作用强度之间的Spearman相关性，以评估整体相互作用结构的相似性。虽然MaxEnt和CellOracle与IGNITE显示出中等程度的正相关性，但SCODE没有显示出显著相关性（见表4）。

此外，我们使用一套指标（S1表）评估了GRN推断、WT生成、针对已知相互作用的验证以及KO扰动分析的效果。接下来，我们利用这些指标来评估先验知识整合和缩放属性的影响。IGNITE仅以scRNA-seq数据为输入，并且不需要先验知识（见图1A），因为它通过最小化CMD来选择GRN。然而，其他方法（如CellOracle）会整合额外的实验数据集（例如ATAC-seq）作为先验知识。因此，我们研究了将经过实验验证的相互作用作为先验知识是否可以提高IGNITE的性能。为此，我们寻找了最大化FCI的GRN。我们找到了十个包含18个实验验证相互作用中的14个的GRN（FCI = 0.78）。为了研究这两种方法之间的关系，我们计算了IGNITE生成的所有GRN的CMD和FCI（见图5A）。最大化FCI的10个GRN（橙色点）与我们之前分析的最小化CMD的GRN（红色点）不同。然而，这两个参数呈负相关（Spearman相关性，p值=，Spearman等级相关测试），这可能表明这两种方法可能识别出具有相似属性的GRN。在小鼠基因网络（GRN）推断中，CMD和FCI选择标准的比较。A. 使用IGNITE从输入数据集（经过对数标准化处理的scRNA-seq，包含PST和MB）推断出的250个GRN中，正确推断的相互作用比例（FCI）与相关矩阵距离（CMD）之间的关系。每个点代表一个GRN，对应于一组特定的超参数。CMD最低的10个模型以红色突出显示，FCI最高的10个模型以橙色显示。https://doi.org/10.1371/journal.pcbi.1014067.g005随后，我们评估了这两个GRN集合的KO预测性能。表5报告了每种选择策略下，十个GRN的模拟和实验KO-WT差异之间的平均Spearman相关性（SEM）。通过最小化CMD选择的GRN在Rbpj、Etv5和三重KO方面显示出一致且统计上显著的正面相关性，而Tcf7l1的预测仅较弱（见方法部分）。相比之下，通过最大化FCI选择的GRN表现较低且变化较大，在FDR校正后相关性不显著。这些结果表明，将实验验证的相互作用作为先验知识并不会提高IGNITE的推断质量，CMD仍然是最可靠的选择标准。下载：PNG更大图像TIFF原始图像表5.最低CMD或最高FCI的十个GRN的模拟和实验KO-WT差异之间的Spearman相关性（平均值±SEM）。P值来自针对零的单侧t检验，并在四种KO条件下进行了FDR校正。显著结果用星号标记（* p<0.05，** p<0.01，*** p<0.001）。https://doi.org/10.1371/journal.pcbi.1014067.t005接下来，我们评估了该方法在基因和细胞数量方面的可扩展性。我们使用IGNITE根据输入的大型scRNA-seq数据，生成了从24到200个基因的GRN。我们在标准笔记本电脑（配备Apple M1 Pro处理器和16GB RAM的MacBook Pro）上进行了这项分析并测量了计算时间。计算时间与GRN大小呈二次方趋势（二次曲线拟合的R2=0.9997），对于最大的GRN，计算时间是可以接受的。为了评估细胞数量方面的可扩展性，我们随机抽取了500个细胞25次（占数据集的5％），并在每批中保留了最佳的GRN。如表6所示，仅从500个细胞推断出的最佳GRN仍然能够高精度地再现Rbpj、Etv5和三重KO的实验KO-WT差异，而Tcf7l1的预测效果较弱（见方法部分）。这些结果表明，即使输入数据有限，IGNITE也能保持KO预测的准确性，同时还能减少运行时间，从而支持其实际的可扩展性。受此可扩展性的鼓舞，我们接下来将IGNITE应用于一个人类多能干细胞（PSCs）的独立数据集，以评估其在不同物种和实验平台上的通用性。下载：PNG更大图像TIFF原始图像表6.从500个随机抽取的细胞中推断出的最佳GRN的模拟和实验KO-WT差异之间的Spearman相关性（平均值±SEM）。P值来自针对零的单侧t检验，并在四种KO条件下进行了FDR校正。https://doi.org/10.1371/journal.pcbi.1014067.t006

2.7 人类GRN推断和WT数据生成
Chu等人[31]的人类PSCs scRNA-seq数据集捕捉了细胞从多能性到中内胚层（mesendoderm）再到终末内胚层（DE）的分化过程（图6A）。该数据集包含在六个时间点（0、12、24、36、72和96小时）收集的758个细胞（S5A图）。选择这个数据集有三个主要原因：i) 它已经被广泛研究，并且之前由Matsumoto等人和Pratapa等人[10]用于GRN推断；ii) 它使用的是Fluidigm C1系统生成的，这是一种与小鼠数据集中使用的10X平台不同的单细胞转录组技术。这个系统捕获的细胞较少，使我们能够用更有限的输入数据来评估IGNITE的性能；iii) 它代表了一个与之前分析的小鼠模型不同的生物学系统。下载：PNG更大图像TIFF原始图像图6. 人类PSC数据集预处理。A. 研究的生物学系统[31]：人类PSC分化的时间动态，涵盖多能性、中内胚层和终末内胚层。该数据集基于Smart-seq scRNA-seq。细胞在六个时间点被采样：0小时（n=92）、12小时（n=102）、24小时（n=66）、36小时（n=172）、72小时（n=138）和96小时（n=188）。B. 对数标准化数据集的PCA。每个点代表一个细胞，按聚类标签着色。C. 数据集的PCA，细胞按伪时间着色。伪时间的范围从0到183.58。D. 对数标准化scRNA-seq数据集中三个关键基因（Pou5f1、T和Gata4）的PCA可视化。每个点代表一个个体细胞，按相应基因的标准化表达值着色。E. 三个代表性基因Pou5f1、T和Gata4的伪时间中的迷你批量基因表达模式。https://doi.org/10.1371/journal.pcbi.1014067.g006hPSC分析遵循应用于mPSC数据集的相同基准框架，包括GRN推断、WT数据生成、针对已知相互作用的验证以及必要的数据集特定指标的KO扰动分析（S1表）。根据IGNITE的要求，我们使用对数标准化方法对数据进行了标准化，将数据投影到低维PCA空间，并对细胞进行聚类以计算伪时间（见方法部分）。我们使用Slingshot算法[32]计算了PST：PST值与实验时间点对齐（图6C），并且每个采样时间内的PST分布确认了预期的顺序（S5B图）。为了评估识别的簇的生物学一致性，我们检查了随时间变化的簇组成（S5C图）以及簇内代表性基因的表达分布（S5D图），这与Chu等人[31]报告的趋势一致。PST排序后，关键的时间动态得到了保留（图6D和6E）：PST较低的细胞显示出较高的多能性标记物表达（例如POU5F1），中内胚层标记物在中间PST时达到峰值（例如T），而DE标记物在较高PST时增加（例如GATA4）。作为下一步，我们定义了GRN节点。我们选择了数据集中所有标记为转录因子（TFs）[50]的基因，保留了在至少一个簇中差异表达显著的基因。这种选择产生了一组92个基因，作为GRN推断的基础（见方法部分），另外还有6个基因用于验证。然后我们通过沿PST聚合相邻细胞来构建迷你批量轮廓。生成的 heatmap（S5E图）突出了不同的时间模式：一些基因在早期达到峰值，而大多数基因在中期或晚期被激活。最后，我们将表达二值化以获得用于GRN推断的基因活动数据集（S5F图）。

2.8 使用GA作为输入，IGNITE推断出250个基因调控网络（GRNs），并计算了每个网络的相关矩阵距离（CMD）（S6A图）。我们选择了最能再现基因之间相关性的GRN（图7A和7B），其CMD为0.35。使用双样本t检验的统计验证确认了IGNITE模型生成的CMD值与零模型生成的CMD值之间存在显著差异（t统计量：，p值：）。推断出的GRN表示为一个交互矩阵，其中每个条目（i,j）代表从基因j（列）到基因i（行）的交互强度（S6B图）。交互矩阵显示了具有正负调控关系的异质交互模式，以及潜在的调控模块和交互模式。我们选择了22个在人类PSCs分化为终末内胚层过程中起作用的基因之间的交互（S6C图）[31,51]。交互矩阵显示，在不同分化阶段的活跃基因组中，正交互占主导地位。相比之下，负交互主要观察到在不同分化阶段表达的基因之间。具体来说，幼稚的多能性基因倾向于抑制中内胚层标记物，例如T的抑制和终末内胚层基因的抑制，例如HNF4A和SOX17的抑制。相反，晚期基因负调控早期表达的因子，例如NANOG和POU5F1的抑制。下载：PNG更大图像TIFF原始图像图7. IGNITE在hPSCs中的推断和验证：GRN和生成的数据。A. 输入数据集（人类PSCs分化scRNA-seq数据，包含LogNorm、PST和MB）中的基因活动皮尔逊相关矩阵。B. 使用IGNITE生成的基因活动的皮尔逊相关矩阵。C. 三个数据集的基因活动PCA：输入数据集（含LogNorm、PST和MB），IGNITE生成的野生型（中间），以及IGNITE生成的POU5F1敲除型（右侧）。颜色渐变代表每个正方形区域内的细胞密度。虚线轮廓表示在输入数据上定义的簇边界，并重新用于生成的数据集。D. 分配给簇（1-5）或排除区域的细胞比例，对于C中显示的相同三个数据集。对于生成的数据，条形图显示了150个模拟数据集的均值（误差条表示SEM）。条形图的颜色与C中的簇相匹配。https://doi.org/10.1371/journal.pcbi.1014067.g007为了验证推断出的GRN，我们选择了之前由SCODE用于验证的TRRUST调控网络中的基因[52]。IGNITE成功地推断出了我们GRN和TRRUST网络之间共享的7个相互作用中的6个（S6D图）。此外，IGNITE生成了与输入数据集维度相同的150个基因活动数据集。我们通过将IGNITE生成的WT数据与PCA空间中的WT基因活动数据进行比较来评估IGNITE生成数据的质量。IGNITE生成的每个二进制数据集都被投影到WT基因活动数据的PCA空间中（图7C）。为了定量评估数据集的相似性，我们使用输入WT数据定义了这个PCA空间中先前识别的簇区域（见方法部分）。簇5和簇6被合并，因为簇6的区域完全包含在从GA数据计算出的簇5的区域内。这个结果是预期的，因为这些簇对应于PST值较高的细胞（S6E图），这些细胞对应于较晚的时间点（72和96小时），在Chu等人[31]中已经观察到单个细胞之间的重叠域。通过利用定义的簇区域，我们验证了IGNITE生成的细胞在不需要引入原始数据中不存在的细胞类型的情况下，再现了相同的激活模式。此外，我们确认IGNITE准确再现了每个簇中的细胞比例，如图7D所示，正确推断了数据集在WT状态下的细胞组成模式。

2.9 人类PSC KO预测
为了评估IGNITE预测人类PSCs中基因敲除后细胞模式变化的能力，我们为GRN中的每个基因生成了敲除条件下的基因活动谱。基因KO的效果体现在WT和KO条件下的基因活动模式的变化中。为了评估这些变化，我们检查了PCA空间中KO生成细胞的分布（图8A），并分析了每个簇中的细胞比例与WT条件的差异（图8B）。下载：PNG更大图像TIFF原始图像图8. 评估IGNITE对hPSC分化中基因敲除的预测。A. WT生成数据和KO生成数据之间的簇组成差异，平均值取自150次重复实验。行对应于模拟的单基因KO，列对应于簇（1-5）或排除的细胞。行按层次聚类排序，这也识别出了11个组，颜色表示出来。红色表示KO后的增加，蓝色表示减少。B. 不同条件下的簇组成PCA：WT输入、WT生成的以及每个单基因KO。每个点对应一个条件，并根据簇（1-5）和排除的细胞中的细胞比例进行定位（生成条件的比例取自150次重复实验的平均值）。颜色表示图8A中定义的组。图中显示的基因对应的KO条件用C中标出。C. 选择与[31]中分析的基因重叠的选定KO基因的分化得分，计算方法如方法部分所述。轴断裂突出了POU5F1的异常值幅度。https://doi.org/10.1371/journal.pcbi.1014067.g008为了验证这些KO扰动效应，我们将IGNITE的模拟结果与Chu等人的分析[31]进行了比较。在这种分析中观察到的预期行为是，T表达在分化早期短暂达到峰值，而CXCR4在后期升高，标志着从中内胚层到DE的转变。Chu等人定义了分化得分（DS）作为CXCR4+细胞百分比与T+（EGFP+）细胞百分比的比率，通过在分化第2天通过流式细胞术测量，并根据WT中的相同比率进行归一化（见方法部分）。得分低于1表示分化受损，而得分高于1则表明分化过程有所改善。为了将IGNITE的预测结果与Chu等人的研究进行比较，我们计算了每个簇中表达CXCR4和T的细胞比例，这些数据包括了野生型（WT）条件（S7A图）和选定基因敲除（KO）条件下的情况。在WT条件下，第4簇中的T+细胞数量减少，而CXCR4细胞数量增加，这一趋势与实验数据中分化第2天的观察结果相符。因此，我们接着计算了DS，即KO条件下第4簇中CXCR4+细胞与T+细胞的比例，并将其标准化为WT条件下的同一比例。这样得到了图8C中展示的各种KO模拟的DS值。与Chu等人的实验数据一致，只有POU5F1的KO显著促进了分化，其DS值大于1。相反，SOX17的KO导致得分接近1，表明对分化没有显著影响，这是预期的结果。所有其他基因的KO都损害了分化，其中KLF8的KO显示出最低的DS值。这些结果表明，IGNITE正确预测了基因KO对分化的总体影响。此外，我们系统地分析了GRN所有成分的KO扰动效应。对于每个模拟的KO，我们计算了WT条件和KO条件下各个簇之间细胞比例的差异。正差异表示KO条件下某个簇内细胞比例的增加，而负差异表示减少。我们对每个条件下（基因KO或WT）每个簇中的细胞比例进行了层次聚类，得到了11种不同的行为模式。促进分化的基因KO组包括第1组（例如包含POU5F1）和第7组，其中的基因敲除导致细胞在后期簇（分别为第5簇和第4簇）中的比例增加。这表明这些基因通常具有分化抑制作用，它们的缺失促进了向更成熟细胞状态的进展。这一效应在POU5F1中尤为明显，因为已知其KO会驱动细胞进入后期发育阶段。相反，第3、4、5和6组包含的基因KO损害了分化（例如EOMEs、ID1、ID2），导致细胞在早期簇（第1、2或3簇）中积累。第8组和第9组（例如KLF8、CXCR4、T）表现出对后期分化的抑制作用，细胞在中间簇中积累，而不是完全停留在最早的分化阶段。第10组和第11组（例如SOX17）包含的基因KO减少了后期簇中的细胞比例，同时略微增加了早期簇中的细胞比例。这表明这些基因在稳定内胚层命运方面起着关键作用。然后，我们将我们的分析结果与已知调控这些发育阶段的实验数据进行了比较。POU5F1的KO促进了分化，减少了早期簇中的细胞比例，同时增加了第5簇中的细胞比例。这一结果与实验观察一致，因为POU5F1的KO已被证明会丧失多能性[31]并加速分化。我们还检查了另一项独立研究[51]的实验数据，该研究观察到多能基因减少而内胚层基因表达增加（S7B图）。NANOG的KO导致早期分化簇的比例增加（S7C图）。这与先前研究[51]的数据一致，表明早期内胚层标记物的表达增加。正向诱导内胚层命运的基因（例如HNF4A、CXCR4、KLF8和SOX17）的KO导致早期簇中的细胞比例增加，这与后期分化阶段相关。总体而言，IGNITE正确预测了基因KO对分化的效应，区分了促进和抑制分化的基因KO。

3 讨论
计算机模拟方法可以用来建模和预测细胞受到干扰（如基因敲除、强制基因激活或外部刺激）的影响。一个高效的算法将能够检查系统上的大量干扰，优先考虑那些最有可能对实验验证产生影响的干扰。许多方法专注于GRN推断[10]或预测干扰效应[7,8]。然而，之前的干扰方法需要在干扰数据上进行模型训练，从而限制了它们的应用范围[7]。此外，各种算法依赖不同类型的数据作为输入，这取决于数据的可用性和兼容性，从而影响了它们的适用性。例如，RE:IN[11]算法需要ChIP-Seq和批量转录组数据，并且需要敲除实验来构建GRN。如何将细胞表型与推断出的GRN联系起来，以及不同条件如何影响细胞行为，仍然是一个挑战。为了解决这一基本的生物学问题，我们开发了IGNITE这种方法，它利用GRN推断来模拟基因干扰。IGNITE的一个重要优势是它能够研究基因相互作用在塑造细胞表型中的作用。我们在两个互补的多能干细胞（PSCs）系统中展示了其能力：从小鼠PSCs从幼稚状态向形成状态转变，以及人类PSCs向内胚层分化[31]。这些数据集在物种、发育轨迹和单细胞平台（10X vs. Fluidigm C1）上有所不同，提供了对通用性的严格测试。IGNITE具有几个优势：(1) 它能够生成野生型数据，复制每个细胞内的基因激活模式；(2) 模拟GRN干扰（如基因敲除），并准确预测其结果，预测KO如何改变GRN中的基因活性；(3) 从未受干扰的单细胞RNA测序数据开始，因此不需要额外的输入（如ATAC-seq数据）；(4) 通过无监督机器学习方法，IGNITE推断出最能描述观察到的细胞表型的定向、加权且有符号的GRN，其中基因组之间的激活模式是可辨别的。这种预测能力为深入理解控制基因相互作用的生物学机制提供了基础。尽管我们没有明确模拟强制基因激活，但IGNITE可以扩展以包含此类干扰。IGNITE的这些优势在应用于分化中的PSCs时得到了体现，它成功预测了野生型或干扰条件下的细胞表型，阐明了PSCs在分化过程中的行为。
为了对IGNITE进行基准测试，我们适配并测试了其他算法以产生可比较的输出。对于我们的分化数据，MaxEnt[16,17]的信息量不足，因为它无法推断网络的方向性。SCODE[18]旨在从分化细胞的scRNA-seq数据中进行GRN推断，但它没有模拟数据中看到的内在变异，因为它使用确定性模型进行推断。此外，SCODE没有利用模型生成系统干扰模拟的能力，即使这是可行的。我们探索了SCODE未测试的生成单基因和三基因敲除模拟的能力。CellOracle[8]使用从单细胞多组学数据推断出的GRN来模拟转录因子的计算机模拟干扰，基于未受干扰的野生型数据模拟细胞身份的变化。然而，这种方法没有描述潜在的过程，因为它只关注GRN中干扰的传播。我们没有使用CellOracle通过网络推断来解析细胞身份，而是将其应用限制在测量KO干扰的效果上，与IGNITE的方法类似。我们包括了三基因敲除的模拟，这是原始工作中未研究的情景。尽管所有方法都能正确推断出GRN中的已知相互作用，但IGNITE在生成野生型和干扰数据方面表现最佳。这是因为IGNITE利用了其生成性的非线性动态来学习潜在的有效GRN属性，而不是依赖于数据约束的GRN推断或KO后的基因表达变化规则。为了进一步评估IGNITE的鲁棒性和通用性，并测试为mPSCs建立的基准框架是否适用于其他物种和实验平台，我们将其应用于一个先前发布的捕捉人类PSCs分化的scRNA-seq数据集[31]。这个数据集在生物学和技术上与原始的小鼠模型不同：它涉及的人类细胞经历了不同的发育轨迹，并使用Fluidigm C1平台生成，该平台捕获的细胞数量（758个）少于10X Genomics平台（9894个）。尽管细胞数量较少且缺乏基因选择的先验知识，IGNITE仍然成功地推断出了一个生物学上有意义的基因调控网络，该网络由从scRNA-seq数据中无偏选择的最具变异性的98个基因组成。推断出的GRN揭示了连贯的调控模式，包括基因组之间的阶段特异性正负相互作用，与已知的分化动态一致。此外，IGNITE有效地再现了野生型基因活性分布，捕捉到了分化过程的全局结构和细胞在各个簇中的相对丰度。值得注意的是，IGNITE预测了在不同阶段阻断分化的KO干扰。一些KO被预测会增加多能细胞的数量，而其他KO则会增加前内胚层细胞的数量，表明不同的基因在分化的不同阶段促进分化。IGNITE正确预测了基因敲除的表型效应，这与两项独立的研究结果[31,51]一致。mPSC和hPSC的分析共同表明，IGNITE可以有效地应用于来自两个不同物种的各种单细胞数据集，即使先前的生物学信息有限或实验设计有所不同。这种灵活性对于研究那些扰动实验不可行且 kurated 知识有限的系统至关重要。然而，IGNITE的性能取决于一些关键的建模假设，这些假设对其使用设定了限制。首先，IGNITE需要正确的伪时间排序。这是其他推断算法（如SCODE）也需要的重要步骤，因为在PST排序后基因表达的时间模式仍然一致。PST和MB允许我们避免使用可能无法精确捕捉实际生物变异性的替代方法，例如假设零基因表达值的插补。其次，IGNITE采用Ising模型将基因表达简化为二元状态，通过将基因近似为自旋。这种抽象足以捕捉与特定细胞行为相关的GRN中的功能有效组织。这些优势在IGNITE应用于分化中的PSCs时得到了体现，它推断出的有效GRN成功预测了野生型或干扰条件下的细胞表型，阐明了PSCs在分化过程中的行为。这种稳定性至关重要，因为它意味着IGNITE能够找到具有生物学意义的真正有信息量的基因 Regulatory Networks（GRNs）。总体而言，IGNITE是一种研究多能干细胞（PSC）分化的宝贵机器学习方法。其适用性可以推广到表征不同的细胞系统及其GRNs，而无需进行实验性干扰，这与其他几个模型[54]相似。相反，IGNITE不仅依赖于未受干扰的数据，还能够对系统干扰进行可靠的计算机预测，即使对于那些结合位点未知的基因也是如此。通过推断即使表达水平较低的基因的调控相互作用，它也解决了单细胞RNA测序（scRNA-seq）数据集固有的数据稀疏性问题。这使得IGNITE成为在数据稀疏或表达水平低的情境下导航基因调控网络复杂性的强大工具。此外，IGNITE提供了一个强大且可扩展的框架，用于揭示基因相互作用如何从单细胞数据中塑造细胞谱型，为预测复杂生物系统中的细胞命运决策提供了途径。

**4 方法**

**4.1 小鼠scRNA-seq数据加工与分析**
- **4.1.1 实验程序**
- mESCs的维护：小鼠胚胎干细胞（mESCs，R1）的培养方法如先前所述[55]。简而言之，基础培养基是HyClone DMEM/F12基础培养基（不加Hepes，Cytiva），含有4mg/mL AlbuMAX II富脂牛血清白蛋白（GIBCO），1× MACS NeuroBrew-21含维生素A（Miltenyi Biotec），1× MEM NEAA（GIBCO），50U/mL青霉素-链霉素（GIBCO），1mM丙酮酸钠（GIBCO），以及1×2-巯基乙醇（GIBCO）。基础培养基中还添加了3.3mM CHIR-99021（Selleckchem），0.8mM PD0325901（Selleckchem）和10ng/mL hLIF（由VBCF Protein Technologies Facility提供），以配制2iL自我更新培养基。细胞在Greiner Bio-One CELLSTAR聚苯乙烯6孔细胞培养多孔板上培养，该板预先用poly-L-鸟氨酸盐酸盐（6g/mL，溶于1xPBS，37°C下孵育1小时，Sigma-Aldrich）处理，随后用Engelbreth-Holm-Swarm鼠肉瘤基底膜蛋白（1.2 mg/mL，溶于1xPBS，37°C下孵育1小时，SigmaAldrich）处理。常规细胞传代比例为1:6，每2天进行一次，使用1×胰蛋白酶-EDTA溶液（Sigma），在37°C下操作，并在基础培养基中加入10%胎牛血清（FBS，Sigma）终止反应。

**4.2 10X基因组scRNA-seq数据与MULTI-Seq条形码标记及验证**
- 10X基因组scRNA-seq的分步操作如[55]所示。细胞通过37°C下胰蛋白酶处理10分钟进行收集。胰蛋白酶处理后使用10% FBS在基础培养基中终止。随后将细胞沉淀重悬于PBS中以去除剩余的FCS，并使用CASY细胞计数器（Biovendis）进行计数。细胞用MULTI-Seq条形码标记的过程如下：每个样本重悬0.5百万个细胞，并与条形码和脂质锚定剂混合物在冰上孵育5分钟。然后引入共锚定剂再孵育5分钟，最后用1% BSA/PBS终止反应。细胞用1% BSA/PBS洗涤两次，之后将细胞重悬于0.04% BSA中。合并所有样本后，通过FACS过滤器套管过滤，并将浓度调整至1百万细胞/mL。根据制造商的说明，只有存活率超过80%的细胞池被用于后续的10X文库制备。MULTI-Seq的鉴定与在Illumina NextSeq550或NovaSeq平台上测序的cDNA文库同时进行。为了确认MULTI-Seq的特异性，使用了ATTO-488和ATTO-590标记的条形码，具体过程见[56]。细胞条形码标记遵循上述程序，并通过FACS分析验证标记的有效性。

**4.1.2 小鼠scRNA-seq数据集**
- cDNA读取结果使用Cellranger映射到mm10参考基因组，并使用CITE-Seq-Count进行多路复用分割到各个滴液中，然后使用Seurat的“Demultiplexing with Hashtag oligos”功能进行细胞分类。数据集中的细胞来自5个时间点（0小时、6小时、12小时、24小时和48小时）。

**4.1.3 质量控制与细胞过滤**
- 为了确保数据质量，我们移除了重复纳尔和阴性细胞。使用scater R包中的perCellQCMetrics函数计算了每个细胞的质量控制指标。使用quickPerCellQC R函数，以中位数的3倍绝对偏差（MAD）作为阈值，识别并移除了异常细胞。此外，表达水平低（在少于20个细胞中观察到的）的基因也被排除在分析之外。完成这些步骤后，数据集包含Ngenes = 13833个基因和Ncells = 9894个独立细胞。

**4.1.4 标准化**
- 我们通过将每个细胞的 total 计数除以相应的大小因子（该细胞的计数总和）来标准化原始基因表达计数。随后，通过添加PseudoCount = 1并进行log2转换来对标准化后的表达值进行对数变换。因此，表示细胞i中的原始基因表达为，我们得到其标准化后的基因表达为。

**4.1.5 降维**
- 为了进行聚类和伪时间推断，同时保留分化过程的时间结构，我们对小鼠PSC数据集执行了降维处理。所有降维步骤都使用相同的输入特征集：一组2318个先前被确定为与此分化过程相关的基因[14]，这些基因是通过基于Pearson相关性的层次聚类识别的，如Carbognin等人所述。我们最初评估了这个基因集的主成分分析（PCA）。PCA只能捕获有限的总方差，并且没有清楚地区分不同时间点采样的细胞（见S1A和S1B图），这与密集采样和分化轨迹的近似线性性质一致。因此，我们采用了基于t-SNE[40]的非线性降维策略，随后使用UMAP[41]，这更有效地捕获了细胞之间的局部邻域关系和时间顺序。具体来说，我们计算了一个二维的t-SNE嵌入，并将其投影到一个二维的UMAP空间中（见S1C图）。UMAP分析显示，从t = 0小时的细胞样本到t = 48小时的样本，细胞沿着UMAP的第一维呈现出一致的时间进展。这种进展是明显的，不同时间点的细胞有明显的分离，形成了从左到右的有序序列。这一观察结果通过检查UMAP空间内24个选定基因的基因表达模式得到了进一步支持（例如，见图1E）。

**4.1.6 聚类**
- 我们在UMAP转换后的数据上使用了基于图的方法进行聚类分析。首先，我们构建了一个包含30个最近邻居的共享最近邻居图（k = 30）。接着，我们使用igraph R包中的buildSNNGraph函数计算了聚类。然后，我们应用Walktrap算法合并较小的聚类以形成较大的聚类，利用时间样本标签来理解在同一时间点生成的细胞组成的聚类。这一过程增强了我们对细胞亚群演变的分析和可解释性。为了验证结果聚类的生物学意义，我们调查了每个聚类中24个相关基因的表达分布（例如，见S1E图中的三个基因）。我们的探索发现了不同的模式： naive基因主要在左侧聚类中表现出较高的活性（标记为2和3的聚类），而形成性基因则主要在UMAP空间右侧的聚类中表现出较高的活性（聚类4、5和6）。

**4.1.7 2Cell-Like细胞聚类**
- 我们调查了聚类1，由于该聚类中的198个细胞的24个选定基因的表达情况，我们无法确定它是一个以naive基因为主还是以早期形成性或晚期形成性基因为主的聚类。我们使用scran R包的findMarkers函数计算了该聚类中上调幅度最大的前100个基因和下调幅度最大的前100个基因。值得注意的是，我们观察到上调基因集中包含了与2Cell-Like阶段相关的基因[42]。这一发现使我们能够确信将这个聚类中的细胞归类为2CLCs。从这一点开始，我们将这个聚类从我们的研究数据集中排除，因为我们对此细胞状态不感兴趣。

**4.1.8 伪时间**
- 我们使用了Slingshot算法[32]来计算细胞的伪时间（PST）。该算法从低维空间开始，并需要每个细胞的聚类标签。为此，我们利用了UMAP空间数据以及聚类标签。我们设置了slingshot函数的输入：我们将聚类“2”指定为起始聚类，聚类“6”指定为结束聚类。这个选择是基于UMAP空间中每个选定基因的高基因表达区域来决定的。此外，我们通过分别生成每个时间点细胞的伪时间值的小提琴图（见S1F图），发现采样时间和PST进展是一致的。

**4.1.9 基因选择**
- 伪时间处理后的数据集包含13833个基因和9696个细胞。我们只关注了文献中认为与PSCs从naive到形成性转变相关且有信息量的24个基因[11,14,38,39]。

**4.1.10 Mini-Bulk**
- 我们将计算数据集中细胞移动平均值的步骤称为Mini-Bulk（MB）：我们采用了一个窗口大小为150个细胞、步长为一的窗口来计算移动平均值。选择这个步长确保了结果数据集的维度与原始数据集相当。我们将MB程序应用于我们的单细胞数据集，并按照伪时间值对细胞进行了排序。最终，我们得到了一个包含24个基因和9547个细胞的数据集（见S1G图）。这个数据集将成为所有GRN构建方法的输入，除了不需要PST和MB的CellOracle方法。

**4.2 人类scRNA-seq数据加工与分析**
- **4.2.1 质量控制与标准化**
- Chu等人[31]提供的人类scRNA-Seq时间序列数据集（0小时、12小时、24小时、36小时、72小时和96小时）包含了758个来自人类胚胎干细胞（ES细胞）分化为定型内胚层（DE）的细胞。原始作者对这些细胞进行了预过滤，保留了那些至少有5000个基因且TPM > 1的细胞，并且表达值已经使用中位数比率标准化进行了标准化。作为额外的质量控制步骤，我们移除了在少于20个细胞中表达的4348个基因。然后，我们进行了对数变换。

- **4.2.2 降维**
- 为了减少人类scRNA-seq数据的维度，我们使用scran R包中的modelGeneVar函数识别了表达变化最大的前200个基因（HVGs）。然后对这些HVGs的对数转换后的表达值进行了主成分分析（PCA），保留了前10个主成分。S5A图中的PCA坐标显示了一个清晰的时间结构，表明主要轴捕捉到了分化过程中的生物学上有意义的进展。

- **4.2.3 聚类**
- 我们使用基于图的方法对PCA降维后的数据进行了聚类。首先，使用scran R包中的buildSNNGraph函数构建了一个共享最近邻居（SNN）图，设置邻居数为k = 35。然后使用igraph包中的Walktrap算法识别了聚类。S5C图显示了聚类标签与时间点之间的对应关系。得到的结构反映了清晰的时间进展，特定聚类在连续时间点占据主导地位，支持了聚类的生物学一致性。

- **4.2.4 伪时间**
- 我们使用Slingshot算法[32]来计算细胞的伪时间（PST）。PCA坐标和聚类标签被用作分析的输入。根据聚类与时间点之间的对应关系，我们将聚类“1”设定为起点，聚类“6”设定为终点。伪时间的进展与实际采样时间一致，这一点通过分别生成每个时间点细胞的伪时间值的小提琴图得到了验证（见S5B图）。

- **4.2.5 基因选择**
- 伪时间推断后的数据集包含14,841个基因和758个细胞。我们专注于与分化过程相关的转录因子。首先，我们从数据集中选择了Vaquerizas等人[50]编制的已知人类TFs列表中找到的基因。接下来，我们使用scran R包中的findMarkers函数进行了差异表达分析，并对各个簇应用了Wilcoxon秩和检验。通过识别在至少一项与其他簇的成对比较中显示出显著上调的转录因子，为每个簇选取了标记物。具体来说，我们保留了在每个簇与所有其他簇的比较中满足以下阈值的基因：FDR<0.05和AUC>0.92。这确保了每个被选中的基因不仅在统计学上显著，而且在至少一项比较中显示出强烈的区分能力。最终选出的转录因子共有92个。由于SALL1、CXCR4、GLIS1、HNF4A、FOXA2和SOX17基因在胚胎内胚层分化过程中具有已知的生物学相关性（参见参考文献[31,51]），因此我们将它们手动添加到最终列表中。因此，最终的转录因子集合共有98个。

4.2.6 Mini-Bulk策略。我们将之前用于小鼠数据集的Mini-Bulk（MB）策略应用到了人类scRNA-seq数据上。在这里，我们将窗口大小设置为50，步长设为1。由此产生的数据集包含了98个转录调节因子和709个细胞（见S5E图），并作为输入用于GNR推断，使用的是IGNITE工具。

4.3 IGNITE假设基因的活性可以通过一个动态Ising模型来有效描述，该模型由一种称为Glauber动力学的简单自旋翻转动态所控制[19]。为了推断GRN交互作用，我们评估了这个模型的逆问题。遵循Zeng等人的工作[20]，我们还假设自旋是异步更新的，并且整个GRN是一个完全连接的网络，存在可能的交互作用。

4.3.1 基因活性：输入数据的二值化。我们将基因表示为一组自旋，其中每个si是一个二进制变量，可以取值+1或-1。因此，为了从Ising模型中推断GRN并运行IGNITE，首先需要将转录组数据二值化。我们从输入数据集（经过PST和MB处理的对数标准化scRNA-seq矩阵）开始。如果用xi(t)表示具有PST值t的细胞的基因i的表达量，那么xi(t)的二值化方法如下：(2)其中si(t) = +1表示基因i是活跃的，-1表示基因i不活跃。根据基因表达量，这种方法计算基因活性（GA）作为活性的二进制度量，不活跃的基因赋值为-1，活跃的基因赋值为+1。我们通过验证这一点来测试这种方法：在PST值较低时，惰性基因主要是活跃的；而在伪时间时间演化的后期步骤中，形成基因主要是活跃的，这是通过计算GA的热图来确定的（见S1H图）。

4.3.2 非对称动力学Ising模型。动力学Ising模型由一个耦合矩阵Jij指定，该矩阵量化了基因间的相互作用并描述了GRN，以及一组外部场hi，这些场决定了给定基因是活跃（hi>0）还是不活跃（hi<0）。IGNITE的目标是从时间序列si(t)中推断出Jij和hi，这两个值被认为在时间上是恒定的。hi和Jij都以温度单位表示。考虑时间离散化，步长为Δt。在模型的时间演化过程中，自旋si是异步更新的。在每个时间步骤t，以一定的概率随机选择一个自旋i，并根据转移概率更新其值：(3)其中π_{ij}(t|t?1)是给定时间t-1的自旋值i′导致自旋i在时间t取新值j的概率，而h_i(t)是作用在时间t的第i个自旋上的有效局部场。

4.3.3 基于数据的模型重建。原则上，为了进行推断，我们需要知道每个时间的自旋值以及给定时间间隔T内的更新时间。然后，从自旋历史开始，我们的目标是重建模型参数hi和Jij。我们可以通过最大化似然来确定这些参数的学习规则[20]。似然可以表示为：(4)对于每个自旋i，更新时间是一个离散的泊松过程。这意味着每个时间步骤t都是更新时间集合的元素，其概率为Δt。因此，更新时间概率与我们要推断的模型参数无关。为了简化，我们假设Δt是已知的，这样在推断过程中就不需要确定它。因此，似然最大化可以在Δt上直接进行。忽略更新时间后，我们可以得到第i个自旋的更新概率的简化形式，这正是Glauber动力学更新规则：(5)然后，似然的对数，即对数似然，可以表示为：(6)其中N_{ii}是基因i的更新次数。我们强调，最大化对数似然等同于最大化似然本身。对方程6求偏导数，我们得到了交互矩阵的学习规则：(7)重要的是，如果我们定义Δt和s0(t)=1，方程7还包括了场hi的学习规则[20]。关键的是，方程(7)需要知道自旋的更新时间，而这通常是无法获得的。实际上，在每个时间间隔中，只有一个是随机选择的自旋进行更新，但它不一定会发生翻转。因此，通常无法仅从时间序列中测量更新时间。为了克服这个问题，我们可以按照Zeng等人的方法进行[20]。我们定义了一个相关函数ρ_ij，其中ρ_ij是随机动力学的所有实现的平均值。相关函数的时间导数为：(8)注意到我们可以区分有自旋翻转和无自旋翻转的情况，我们可以写成：无翻转项消失是因为自旋没有更新，因此ρ_ij=0。由于Δt是实验轨迹中最小的时间间隔，我们可以写成：(9)这样(10)更新规则（方程7）可以通过将方程右侧的第一项替换为方程10来重写：(11)其中方程11右侧的第二项是对所有时间点的平均值。我们注意到这一项等同于方程7中的平均值，因为Δt不影响是否进行了更新。由于时间间隔的选择对推断无关紧要，我们可以不失一般性地设置Δt=1。虽然不能保证收敛到全局最大值，但方程(11)中的更新规则将达到一个局部最大值，该局部最大值对应于代表功能GRN的耦合矩阵。我们注意到交互矩阵J包括所有基因的自耦合。然而，测量mRNA并不能区分目标基因和调节基因。因此，自耦合缺乏直接的生物学意义。然而，这种计算上的包含虽然不反映直接的生物学重要性，但通过引入一个有信息的偏差来稳定推断过程，这可能补偿了模型未能直接捕获的特定动态或相互作用。基于这个原因，我们推断出自耦合，然后从GRN中移除它们。

4.3.4 优化算法。为了找到方程11的最大值，我们实现了两种优化方法：Momentum梯度上升（MGA）算法（实现了随机梯度下降算法[57]的最大化版本）和Nesterov加速自适应矩估计算法（NADAM）[58]。通常，似然的最大化是从Jij和hi的估计值开始的，并沿着梯度方向用一个小值更新这个估计值。这一过程重复进行，直到收敛。学习规则的每次循环称为一个 epoch。对于MGA算法，我们采用了步长衰减的学习率和L1正则化。第n个epoch的学习率为：(12)其中θ_0是初始学习率，δ_r是每ndr epoch学习率的减少量。场的第n次优化步骤可以表示为：(13)(14)耦合的第n次优化步骤为：(15)(16)方程13和15中的参数λ是动量参数，α是正则化参数。对于NADAM算法，我们实现了步长衰减的学习率（方程12）和L2正则化。为了计算耦合和场的 gradient 的一阶矩m和二阶矩v，我们可以将第n次优化步骤写为：(17)(18)其中g_i是1或-1（见方程11）。为了执行学习步骤，我们需要引入校正偏差的矩：(19)(20)以便：(21)(22)其中λ是L2正则化参数，η和α是NADAM算法的固定参数。

4.3.5 优化算法的超参数。为了从方程11重建GRN，我们需要在MGA和NADAM优化算法中设置以下超参数。对于MGA优化器：步长衰减的学习率参数θ_0和δ_r、δ_n，以及L1惩罚项λ。对于NADAM优化器：步长衰减的学习率参数θ_0和δ_r、δ_n，以及L2惩罚项λ。我们还注意到，优化问题不是凸的，因为我们预期方程6会有许多局部最大值。为了避免次优解，我们在Ntrial=250组随机选定的超参数中进行了随机搜索。考虑的超参数可能值显示在表7中。最佳超参数集可以在有或没有先验知识的情况下定义。有先验知识时：我们知道18个交互作用，并选择FCI（拟合系数）最大的超参数集。如果有多个模型具有相同的FCI，我们随机选择一个。没有先验知识时：我们选择CMD（一致性矩阵）最小的超参数集。如果有多个模型具有相同的CMD，我们随机选择一个。

4.3.6 与空模型的比较。我们进行了排列测试来评估推断出的交互作用的统计显著性。排列测试涉及将基因表达数据的行和列索引洗牌Ntest=50次，然后为每个测试数据集评估Nsets=50组超参数，以确保对不同超参数配置进行彻底测试。使用IGNITE推断出的最终空GRN数量为2500个。我们研究了使用IGNITE推断出的GRN的统计显著性。我们假设随机生成的网络将表现出与空模型相似的交互作用值。我们通过将IGNITE预测的这些交互作用的值与这些空模型GRN的概率分布进行比较，来评估IGNITE正确推断出的12个已知交互作用的显著性。具体来说，我们检查了IGNITE预测的这些交互作用的值是否位于洗牌数据的第五到第九十五百分位数范围内。由于只有三个交互作用与空模型的预期一致（Nanog-Otx2、Stat3-Gbx2和Tcf7l1-Esrrb），我们拒绝了零假设，支持了IGNITE GRN的统计显著性。

4.4 MaxEnt重建。MaxEnt交互矩阵中基因i和j之间的元素可以计算为：(23)其中Φ是通过对我们的转录组数据集计算得出的协方差矩阵的逆矩阵，如[16,17]所述。交互矩阵Φ是对称的。

4.5 SCODE重建 SCODE [18]是一种从scRNA-Seq数据中推断调控网络的方法，专为区分PSCs而设计。它利用线性常微分方程（ODEs）和线性回归来捕捉观察到的基因表达动态。我们选择这种方法作为与IGNITE可比的黄金标准方法，因为正如[10]中强调的，它具有与IGNITE相似的特点：(i)它需要时间有序的细胞来捕捉网络动态，(ii)推断出的交互作用是有方向的（不对称交互矩阵），(iii)是带符号的，以及(iv)可以在野生型或受扰动条件下生成新数据。在[10]中评估的其他模型中，SCODE是具有这些特征且在实验性单细胞RNA测序（scRNA-seq）数据集上准确度最高的模型。SCODE的效率基于一个基本概念：表达动态的模式是有限的，并且可以使用有限数量的这些模式来准确重建表达动态。因此，表达动态向量应该从Ngenes个基因的长度减少到D长度。D是一个模型参数，其值应根据作者的建议进行优化。为此，我们将数据集分为训练集和测试集（占总数据集细胞的20%）。我们将SCODE应用于训练数据，并通过计算测试集的残差平方和（RSS）来评估优化后的模型的有效性，对不同的D值（D=2、4、6和8）进行了测试。对于每个D值，像在原始工作中一样，我们独立执行了100次SCODE。我们观察到RSS值的中位数在D=6时几乎达到了饱和。按照作者的建议，我们检查了推断出的交互矩阵的可重复性。使用每个D的测试数据，计算了50次重复实验中RSS值最低的优化交互矩阵之间的相关系数。直到D=6，重复实验之间的相关性很高，因此优化后的交互矩阵在D=6时仍然稳定。考虑到测试数据的RSS值饱和以及D=6时估计的交互矩阵的稳定性，我们选择了这个值用于我们的分析。最终的推断出的基因网络（GRN）交互矩阵是前50次重复实验交互矩阵的逐元素平均值。

4.6 CellOracle应用程序
CellOracle [8]是一个机器学习框架，它利用从单细胞多组学数据中推断出的GRN来模拟基因扰动及其对细胞身份的影响。模拟扰动分为四个步骤：（1）通过对多组学数据进行聚类正则化线性回归，生成特定于细胞类型或细胞状态的GRN。（2）使用GRN模型模拟这些变化导致的级联效应，估计目标基因表达的变化。CellOracle的扰动建模是确定性的。（3）通过比较基因表达的变化与相邻细胞的变化来计算细胞身份转换概率。（4）将这些概率转换为加权向量，将复杂的基因表达变化简化为二维表示，以预测TF扰动后的细胞状态转换。为了将CellOracle与IGNITE和SCODE进行比较，我们专注于CellOracle的前两个步骤，得到的输出是推断出的GRN和KO扰动后的模拟基因表达。CellOracle的设计仅能生成扰动后的基因表达数据，而不能生成野生型（WT）数据。CellOracle需要输入已进行降维和聚类的处理过的scRNA-seq数据集。我们使用的输入数据是按照5.1节中的描述处理过的数据，直到去除了2CLCs簇。处理后的数据集包含9696个细胞和2078个基因。这些基因由Carbognin等人[14]鉴定为在分化过程中具有显著性，由于它们与研究过程的相关性已得到证实，因此选择了这些基因而不是高度变异的基因。我们使用一个全连接网络作为基础GRN，并用CellOracle推断这些连接的权重（不包括自环，因为CellOracle不推断自环）。我们没有使用作者建议的预构建GRN，因为在这个网络中大多数已知的相关性都不存在（78%）。然后，我们使用CellOracle Python包中的simulate_shift函数生成了扰动后的KO数据。在GRN推断和数据生成之后，我们仅关注了对我们的研究至关重要的24个基因及其相应的GRN。

4.7 正确推断交互的比例（FCI）
正确推断交互的比例（FCI）是根据推断出的交互值计算得出的。FCI的计算公式为：
分子是通过统计推断出的交互值的符号与文献中已知交互符号匹配的实例数量来确定的。分母等于18。

4.8 WT数据的生成与评估
4.8.1 使用IGNITE生成基因活性数据
我们使用IGNITE方法推断出的交互矩阵生成了基因活性数据。这个数据集中的基因，与IGNITE的输入一样，被表示为自旋（spins）。为了生成一个新的数据集，我们利用了伊辛模型参数、推断出的交互矩阵（J）和外部场（h），通过Glauber自旋更新规则来确定不同细胞的基因状态。结果数据集的结构与原始数据相同，行代表基因，列代表对应于不同时间步长的不同细胞。数据维度与输入数据相同。为了生成数据，我们为网络中的每个基因初始化了随机的自旋值。在每个时间步长，我们根据每个基因与其他基因的相互作用独立更新自旋的值。根据Glauber动力学，每个自旋都有概率翻转。

4.8.2 使用SCODE生成基因表达数据
我们使用了SCODE作者提供的代码来重建基因表达。我们生成了一个新的数据集，捕捉了24个选定基因在100个时间步长内的动态变化，如作者建议的那样。每个基因的初始基因表达值设置为原始对数标准化scRNA-seq数据集前1000个细胞的平均表达值。由于数据生成过程是确定性的，我们只生成了一个数据集，以确保所有生成的数据集都是相同的。

4.8.3 相关矩阵距离
我们按照以下方法评估了使用IGNITE或SCODE生成的WT数据集与输入数据集之间的相似性：
i. 使用IGNITE生成了Ntrial=250个基因活性数据集，其维度与原始基因活性数据集相同（Ngenes=24，Ncells=9547）。由于SCODE是确定性的，因此我们只生成了一个细胞数量符合作者建议的基因表达数据集（Ngenes=24，Ncells=100）。
ii. 我们计算了输入数据集和每个模拟数据集的皮尔逊相关矩阵。
iii. 仅对于IGNITE，我们计算了模拟数据中Ntrial个相关矩阵的逐元素平均值，得到了一个平均相关矩阵。
iv. 我们按照以下方式计算了输入数据集和生成数据集之间的距离：
（25）
其中xij表示输入基因活性/表达数据的相关矩阵中的一个元素，表示模拟数据的相关矩阵（针对IGNITE平均后的）中的一个元素。
v. 为了缩放距离度量d，我们通过随机打乱输入基因活性数据和基因表达数据集来创建了Nshuffle=250个随机数据集，如之前为IGNITE空模型所描述的那样。然后，我们通过将其除以输入数据的相关矩阵与空模型数据集的平均距离来缩放d。我们称缩放后的d量为相关矩阵距离（CMD）。一个好的推断方法将产生CMD小于1的模拟数据。CMD=1的值意味着模拟数据的相关矩阵与偶然发生的情况相当。为了量化IGNITE和SCODE结果的质量，我们进行了两个统计测试：
（i）对IGNITE生成的d值与空模型生成的d值分布进行双样本t检验。
（ii）使用z分数比较SCODE获得的单个d值与空模型d值的分布。

4.9 数据的聚类
我们通过根据细胞模式的相似性对数据集中的细胞进行分组来进行聚类。这使我们能够探索基因表达或基因活性数据集的潜在结构。我们首先使用SciPy库中的scipy.cluster.hierarchy.linkage函数对数据集进行了链接分析。我们采用了Ward方差最小化算法进行链接分析。这是一种分层聚类方法，旨在最小化每个簇内的方差。随后，使用 dendrogram和热图可视化了链接分析的结果。热图直观地展示了整个细胞群体中的基因表达/活性模式。dendrogram表示了簇的分层结构，说明了单个细胞或细胞群体如何根据它们的基因表达/活性相似性合并成更大的簇（图2D为IGNITE的输入基因活性，图2E为IGNITE生成的基因活性，S2C图为SCODE的输入基因表达，S2D图为SCODE生成的基因表达）。

4.10 系统的基因敲除（KO）扰动
我们通过实施单个基因（针对Rbpj、Etv5或Tcf7l1）和三个基因（同时针对Rbpj、Etv5和Tcf7l1）的敲除模拟来扰动系统，以评估IGNITE、SCODE和CellOracle算法的预测能力。

4.10.1 KO数据的生成
使用IGNITE：为了在单个或三个基因KO后生成扰动数据，我们从推断出的交互矩阵中移除了KO基因的行和列。然后我们像在4.8.1节中描述的那样生成了Ntrial=250个基因活性数据集。这些数据集包含NKN个KO基因。细胞数量（Ncells=9547）与输入数据集相同。
使用SCODE：我们通过从SCODE推断出的交互矩阵中移除KO基因来生成模拟的KO。然后我们按照4.8.2节中描述的程序为每个KO生成了一个基因表达数据集。细胞数量为Ncells=100，这是SCODE作者建议的。由于SCODE的动力学是确定性的，我们只生成了一个数据集。
使用CellOracle：我们使用CellOracle Python包中的simulate_shift函数生成了模拟的KO数据。与SCODE一样，对于每个CellOracle KO，我们也生成了一个基因表达数据集。

4.10.2 测量KO基因对GRN的影响
为了评估基因敲除对基因活性或表达模式的影响，我们比较了WT数据集和KO数据集，考虑了单个和三个基因的KO。对于IGNITE，我们生成了250个数据集，如4.8.1节中详细描述的。对于SCODE，我们生成了一个WT数据集，如4.8.2节中描述的。对于CellOracle，WT数据集是用于推断GRN的输入数据集，即经过KNN插补处理的原始数据。作为输入和生成的数据，我们在考虑IGNITE数据时使用了基因活性数据，在考虑SCODE和CellOracle数据时使用了基因表达数据。基因活性可与用于IGNITE的输入数据相比较：经过对数转换的scRNA-seq数据，应用了伪时间（pseudotime）、Mini-Bulk和二值化处理。相反，使用SCODE和CellOracle生成的基因表达数据可与这些方法所需的输入数据相比较：对于SCODE，使用了经过对数转换的scRNA-seq数据，应用了伪时间和Mini-Bulk；对于CellOracle，使用了经过对数转换的scRNA-seq数据，并应用了伪时间。为了将WT数据集与KO数据集进行比较，我们移除了扰动基因的数据。这种影响的量化方法如下：首先，我们计算了所有细胞中每个基因（以及在所有IGNITE生成的数据集中）的基因活性/表达水平的平均值。接下来，我们计算了每个基因i的KO和WT基因活性/表达之间的差异：
（26）
其中表示KO数据集中基因i的平均基因活性/表达，表示WT数据集中的相同量。负值表示KO扰动导致基因i的平均基因活性/表达降低，而正值表明KO导致该基因的基因活性/表达增加。相反，表示KO没有影响基因i的基因活性/表达。实验数据。考虑到对单个基因j或三个基因{j,k,l}的扰动，为了测量这种扰动对实验数据的影响，我们采用了[39]中的log2FC指标。负值表示由于扰动导致基因i的表达减少，而正值表示表达增加。接近零的log2FC值意味着对基因i的影响可以忽略不计。对于三基因敲除（KO）的实验数据，我们使用了[38]中的数据集。为了量化扰动对基因表达的影响，我们计算了每个基因的两个log2FC值之间的差异：一个是三基因KO与对照（2i）条件之间的差异，另一个是去除2i后72小时与0小时时WT（野生型）基因表达之间的差异。这仍然是一个log2FC指标，它能够更好地反映KO对基因表达影响的相对程度，相对于WT条件下的基线表达。4.10.3 比较KO指标。为了确保模拟数据的值与实验数据的log2FC值之间有意义的比较，我们对这两个量进行了归一化。对于生成的数据，所有从四个KO模拟中获得的值都除以了所有敲除实验中值的最大绝对值。同样，单基因KO和三基因KO的Log2FC值也分别除以了它们的最大绝对值。这种归一化程序确保了所有值都是可比的，并且落在一定的范围内。值得注意的是，我们避免了将值和log2FC值都除以它们的共同绝对最大值，因为生成这两个数据集的步骤不同，直接比较可能不合适。接下来，为了清晰起见，我们将这两个归一化后的值（和log2FC）称为缩放的KO-WT差异。KO实验的同意度（FoA）。为了评估模拟和实验KO-WT变异指标之间的一致性，我们定义了同意度（FoA）。为了计算FoA，我们设定了一个阈值thr = 0.05。随后，我们将差异值分为三类：如果则为正，如果则为负，如果落在范围内则为零。然后我们统计了每个基因i的生成值和实验值在分类范围内匹配的次数（正数、负数或零）。这个计数然后除以基因的总数，得到FoA值（对于单基因KO模拟为23，对于三基因KO为21）。为了评估观察到的FoA值的统计显著性，我们使用二项分布作为空白模型。试验次数为Ntrials。与之前一样，我们将生成值和实验值落在同一分类范围内的情况定义为成功。在任何给定试验中成功的概率是p = 1/3。为了确定FoA值的显著性阈值，我们计算了基于二项分布观察到等于或更高FoA值的累积概率小于或等于0.05所需的最小FoA值。这个阈值定义了被认为是统计上显著的FoA值的下限。对于单基因KO和三基因KO情况，确定的阈值为FoAthr = 0.48。KO预测性能的统计分析。对于每个KO条件，计算了评估中的一组GRN（无论是具有最低CMD的10个GRN、具有最高FCI的10个GRN，还是来自25个子集样本的最佳GRN）上模拟和实验KO-WT差异之间的Spearman相关性。通过单侧t检验对相关性与零进行统计显著性评估。P值使用Benjamini–Hochberg FDR程序进行了多重检验调整。4.11 人类系统的扰动：基因敲除分析 4.11.1 PCA空间中的簇区域估计。为了计算PCA空间中每个簇对应的区域（图7C和7D，以及图8A），我们通过计算Delaunay三角剖分的局部边连接来估计平滑的簇边界，并只保留那些短于固定alpha阈值（我们设为α）的边[59]。得到的边集用于定义一个外边界，通过保留段的凸包来近似。为了考虑边界位置的不确定性——由于局部稀疏或投影噪声——我们通过应用一个固定宽度的几何缓冲区（在PCA空间中为0.5单位）来扩展每个alpha形状。这个缓冲区是通过原始多边形和半径等于缓冲区大小的圆的Minkowski和计算得出的，它在几何上增厚了边界，确保边细胞包含在簇区域内。缓冲区宽度的选择是基于两个因素的平衡：(i) 包含有意义的边界细胞，以及 (ii) 最小化相邻簇之间的重叠。在缓冲区区域重叠的情况下，使用预定义的优先顺序来解决冲突。未分配到任何区域的细胞被标记为“排除”。4.11.2 生成GA数据中每个簇的细胞比例。我们生成了150个WT基因活动数据集，其维度与输入的WT基因活动数据相同。每个数据集都独立处理，计算分配给每个簇的细胞比例。最终簇组成报告为所有生成数据集的平均值及其标准误（SEM）（图7C和7D）。4.11.3 与Chu等人的KO行为比较。为了评估基因敲除（KOs）对分化动态的影响，我们采用了Chu等人最初引入的一个指标，称为分化评分（DS）。该评分使用了KO和WT条件下分化开始两天后（大约45-48小时）的CXCR4+和T+细胞的比例。其定义为：(27)其中和分别代表KO条件下45-48小时的CXCR4+和T+细胞的比例，和是WT条件下的相应比例。DS > 1表示KO后分化加速，而DS < 1表明分化延迟或受损。我们将这个指标应用于我们生成的IGNITE KO数据。为此，我们首先确定了数据集中与Chu等人工作中第二天的分化最接近的簇。我们比较了各簇中表达T或CXCR4（基因活性为1）的细胞比例（S7A图）。在原始研究中，32.9%的WT细胞是CXCR4+，6.44%是T+。在我们的数据中，簇4最符合这一特征。因此，我们使用分配给簇4的WT和KO细胞计算了分化评分。4.11.4 与Wang等人的POU5F1 KO行为比较。为了验证我们的KO模拟结果，我们重新分析了Wang等人[51]的微阵列基因表达数据，该数据可在GEO访问号GSE34921下找到。这个数据集包括了人类多能因子敲除（shRNA介导）后PSCs的表达谱。在可用的KO条件下，我们重点比较了POU5F1 KO，因为这个基因既存在于我们推断的GRN中，又与实验数据集中的明显分化表型相关。我们关注了一个包含多能性和分化相关基因的精选列表，对应于S6C图中的基因。比较使用了从分化第1天到第8天的8个重复实验，包括WT和POU5F1 KO条件。我们计算了每个基因在重复实验中的平均表达值，并用它来计算S7B图中的z分数。支持信息 参考文献