探究基因调控网络中的转录因子子集

《Algorithms for Molecular Biology》：Probing transcription factor subsets in gene regulatory networks

摘要
转录因子（TFs）在基因表达的调控中起着至关重要的作用，因此在疾病研究以及特定细胞功能或途径的操控方面一直备受关注。TFs与基因之间的复杂相互作用可以被建模为一个网络，其中每条边代表一种调控作用。在这样的网络中，识别出具有最大调控作用的TFs集合是一个挑战，这在扰动研究的背景下尤为重要。我们将寻找一组能够最大化对一组基因影响的TFs的任务建模为二分图中的探测问题。我们在模拟图上测试了不同的自适应和非自适应算法，并讨论了它们的特性和适应性差距。然后，我们将这些算法应用于实际数据，以找到参与T细胞介导的免疫和淋巴细胞白血病的调控基因的TFs。该方法对数据的需求最小，并且可以很容易地应用于所有其他类型的二分网络。

引言
转录因子（TFs）是通过以序列特异性方式结合DNA来影响基因表达的调控蛋白。TFs特异性地结合染色质的可访问区域，这些区域也被称为调控元件（REs）[1,2,3]。由于TFs在众多生物过程中的核心作用以及它们与疾病的关联，它们一直是研究的重点。要理解TF的功能，了解其靶基因是必不可少的。有多种方法可以将TFs映射到基因上，例如依赖于收集TF-基因相互作用的数据库[4,5,6,7]，在一系列样本或单个细胞中关联基因和TFs的表达[8, 9]，比较不同条件下的变化[10, 11]，观察TFs扰动的后果[12,13,14,15]，或者通过注释来自实验数据（如ChIP-seq [16, 17]）或基于已知序列结合偏好的预测的转录因子结合位点（TFBS）[2]。由于一个基因不仅受到其附近TFs结合的影响，因此应该将基因的所有REs都视为潜在的TFBS位置[3]。基因和TFs之间的关系通常通过网络来建模，其中每条边代表一种调控作用[9, 18,19,20,21]。

在这样的网络中，识别出调控一组感兴趣的基因（例如，特定途径或疾病的基因）的TFs是一个挑战。通常，TFs的重要性是逐一评估的[22,23,24,25]。然而，已知TFs会协同作用，且途径可以由一组TFs调控[1,2,3, 26]。一个著名的例子是四个Yamanaka因子Oct3/4、Sox2、c-Myc和Klf4[27]。在网络中识别一组TFs，无论是为了提出关于途径如何被调控的假设，还是为了寻找扰动实验的潜在目标，都可能成为一个组合问题。例如，从400个TFs中挑选出三个TFs会产生超过1000万种可能的组合。很少有工作明确旨在寻找TFs的集合。Gross和Blüthgen[28]旨在通过将其描述为一个最大流问题来识别能够最大化对生物网络内关系理解的扰动组合。类似地，MEED工具建议进行能够最大化对调控关系理解的扰动实验，使用预测逻辑模型[29]。Wang等人[30]根据扰动数据集为每个基因构建了TF调控模型，随后使用这些模型来定义最有助于理解基因调控的扰动集合。虽然这些方法有助于设计验证网络的实验，但它们并不旨在找到在网络中覆盖基因最多的TFs集合。

我们通过解决一个随机探测问题来解决寻找一组能够最大化目标基因数量的TFs的任务，即被TFs调控的基因。探测问题要求从给定的基础集合中选择一组元素，以最大化给定的随机目标函数[31,32,33]。形式上，对于给定的n个随机变量集合\(A = \{X_1,\ldots ,X_n\}\)和一个随机函数\(f:2^A \rightarrow \mathbb {R}\)，目标是选择一个子集\(S \subseteq A\)，使其满足某些约束条件，以便f(S)达到最大值（期望值）。通常，选定的子集受到大小的限制[31, 32]。这个问题有两种类型的算法：非自适应算法事先决定要探测哪些元素；而更通用的自适应算法可以在做出下一个探测决策时考虑最近探测的结果。这使得自适应算法通常比非自适应算法更强，但往往也更加复杂和昂贵。为此，研究适应性差距（即最优自适应算法与（次）最优非自适应算法的预期目标值之间的最坏情况比率）是有趣的。

在这项工作中，我们首先模拟了具有不同设置的TF-基因网络作为加权二分图，并测试了自适应和非自适应算法来选择一组TFs（图1）。我们使用不同的目标函数评估了算法的性能，讨论了它们的特性，并分析了它们的适应性差距。尽管计算最优解的成本非常高，但我们考虑了高效的非自适应算法。它们的结果与更昂贵的自适应算法相比表现良好。换句话说，我们提出了具有接近最优解的高效探测算法，克服了测试大量组合的问题。我们的模型本质上等同于经典的随机探测模型[31, 34, 35]，因为任何目标函数也可以用二分图来定义。然而，现有的工作仅将目标函数形式化在一组随机变量上。我们专门化了模型，以便于更清晰和直观的表述，这在使用边不独立的目标准时尤其相关。

在应用于实际数据时，我们构建了涉及T细胞介导的免疫和淋巴细胞白血病相关基因的网络。我们的目标是找到具有最佳网络覆盖率的TFs三元组。更准确地说，我们希望找到网络中总边权重最高的TFs组合，而每个基因只能考虑一个TF。识别出的TFs对基因集具有特异性，在三元组内表现出显著的共表达，并且我们可以用文献证据来支持它们的相关性。我们的网络构建方法可以轻松应用于不同的数据，因为唯一需要的数据类型是特定样本的REs活性测量，例如ATAC-seq、DNase1-seq或ChIP-seq测量的RE相关组蛋白修饰。据我们所知，没有其他工作优化了选择TFs子集以最大化调控基因覆盖率的策略。此外，探测框架也可以应用于从任何其他可用方法派生的TF-基因网络。原则上，它也可以用于任何其他类型的二分网络，例如REs与其调控基因之间的相互作用或micro RNAs与其mRNA目标之间的相互作用。代码可在GitHub上获取。

图1
此图像的替代文本可能是使用AI生成的。

研究概述
转录因子（TFs）及其靶基因被建模为二分网络。在模拟图中测试了不同的算法、目标函数和图设置。在应用于实际数据时，基于T细胞的H3K27ac ChIP-seq数据构建了两个网络。TFs和基因之间的边依赖于基于基序的转录因子结合位点（TFBS），考虑了基因的所有预测调控元件（REs）[36]。目标是找到在T细胞介导的免疫和淋巴细胞白血病中具有最佳基因覆盖率的TFs集合。

方法
随机图模型
我们考虑使用图论方法来表达引言中描述的调控因子探测问题。在我们的模型中，有一个加权的完全二分图\(G=(A \cup B, A \times B)\)，其中包含不同的节点集\(\mathcal {A}\)和\(\mathcal {B}\)。我们将A的元素称为调控因子，将B的元素称为位置。我们使用\(n_{\mathcal {A}} :=|\mathcal {A}|\)和\(n_{\mathcal {B}} :=|\mathcal {B}|\)来表示它们的数量。在生物数据的应用中，调控因子代表TFs，位置代表基因。

一般来说，每个调控因子\(a \in \mathcal {A}\)可以影响每个位置\(b \in \mathcal {B}\)。这种影响的程度由一个强度或边权重\(w_{a,b} \in [V] = \{0,1,\ldots ,V\}\)表示，其中\(V \in \mathbb {N}\)。这里\(w_{a,b} = 0\)表示完全没有影响。为了捕捉影响的不确定性，每条边\((a,b) \in \mathcal {A} \times \mathcal {B}\)都伴随着一个在\(\mathcal {V}\)上的离散分布\(D_{a,b}\)。为了简化，我们假设随机权重\(w_{a,b}\)是根据\(D_{a,b}\)独立抽取的。最初，我们知道\(D_{a,b}\)，但不知道实现值\(w_{a,b}\)。

探测和目标
对于我们的探测问题，我们给定了两个整数参数k和\(\ell\)，其中\(\ell \le k \le n_\mathcal {A}\)。我们的目标是选择一个最多包含\(\ell\)个调控因子的子集\(S \subseteq \mathcal {A}\)。我们通过一个目标函数\(f:2^A \rightarrow \mathbb {R}_{\ge 0}\)来评估这个子集S的质量。这个函数取决于所选调控因子相关的实际边权重。因此，我们的目标是找到一个最大化该函数的子集S。我们说有一个\(\ell\)的基数约束，因为只有最多包含\(\ell\)个调控因子的子集是可行的。

为了确定一个好的调控因子子集，我们允许探索或探测最多k个调控因子，其中\(\ell \le k \le n_\mathcal {A}\)。在探测一个调控因子a后，所有相关边的实际权重都会被揭示，即我们可以看到每个\(b \in \mathcal {B}\)的\(w_{a,b}\)。给定\(k \ge \ell\)个被探测的调控因子，我们可以选择一个在f(S)方面尽可能好的\(\ell\)个被探测调控因子的子集S。

形式上，我们在目标函数的扩展定义中表达了后一个选择步骤，该定义直接以所有k个被探测调控因子的集合\(S'\)作为参数。定义1：对于给定的目标函数\(f:2^{\mathcal {A}} \rightarrow \mathbb {R}_{\ge 0}\)和基数约束\(\ell\)，我们定义其父函数为$$f^*:2^{\mathcal {A}} \rightarrow \mathbb {R}_{\ge 0}, \qquad S' \mapsto \max _{S \subseteq S':|S| \le \ell } f\left( S\right) 。$$ 因此，我们的目标是选择一个大小为k的子集S'来进行探测，以最大化\(f^*(S')\)。

为了模拟真实的基因调控网络，我们选择了一个覆盖函数作为目标函数。定义2：覆盖函数定义为$$f_{cov}(S) = \sum _{b \in \mathcal {B}} \max _{a \in S} w_{a,b}\hspace{5.0pt}。$$ 这个名称来源于（加权的）MaximumCoverage问题[37, 38]：在探测了一组调控因子后，我们必须选择一个子集，以最大化总覆盖权重，其中每个位置最多只能被一个调控因子覆盖。因此，所选的调节因子之间是相互依赖的，这正是在TF基因网络中所期望的。我们探索了两个额外的目标函数\(f_{max}\)和\(f_{sum}\)，这两个函数将调节因子视为完全独立的。这些函数是一大类调节因子加性可分函数的例子，我们在补充材料中详细定义和分析了这些函数。定义3：我们定义$$f_{max}(S) = \sum _{a \in S}\max _{b \in \mathcal {B}}w_{a,b} \quad \text {和} \quad f_{sum}(S) = \sum _{a \in S}\sum _{b \in \mathcal {B}}w_{a,b}$$，分别表示集合S中所有调节因子的最大值和所有相邻边权重的总和。

我们的探测问题可以很容易地被构建为一个有限马尔可夫决策过程[39]。因此，我们区分了两类探测算法（通常也称为探测策略或政策）：自适应和非自适应的。自适应策略会依次做出k次探测决策，并在每一步中使用之前步骤中实现的边权重的反馈。相比之下，非自适应策略会预先决定完整的k次探测集合\(S'\)。注意，在这两种情况下，最终从探测元素集合中选择\(\ell\)个元素都可以基于在探测过程中观察到的所有实现的边权重。

任何自适应策略都可以用决策树来描述。由于我们的问题是一个有限马尔可夫决策过程，因此存在一个最优的自适应策略\(\pi ^*\)。然而，这样的策略通常计算成本非常高。相反，我们遵循关于探测的文献，并采用概念上和计算上更简单的非自适应策略。我们通过标准的适应度差距方法来衡量目标值的恶化。定义4：设\(\pi ^*\)是一个最优（自适应）策略，\(\pi '\)是一个最优的非自适应策略。假设\(S'_{\pi ^*}\)和\(S'_{\pi '}\)分别是\(\pi ^*\)和\(\pi '\)执行的k次探测的（随机）集合。适应度差距由预期目标函数值的比率\(\mathbb {E}[f^*(S'_{\pi ^*})] \, / \, \mathbb {E}[f^*(S'_{\pi '})]\)给出。在本文的其余部分，我们将策略称为算法。

在过去二十年里，随机探测问题中的适应度差距一直备受关注[40,41,42]。特别是，单调次模函数\(f^*\)是一个重要的研究领域，受到了很多关注。定义5：一个非负的集合函数\(f:2^{\mathcal {A}} \rightarrow \mathbb {R}_{\ge 0}\)被称为单调的，如果对于任意两个子集\(S \subseteq T \subseteq \mathcal {A}\)，我们有\(f(S) \le f(T)\)。如果对于任意两个子集\(S \subseteq T \subseteq \mathcal {A}\)和\(j \in \mathcal {A}\)，我们有\(f(T \cup \{j\}) - f(T) \le f(S \cup \{j\}) - f(S)\)，即向集合中添加一个元素的边际价值随着集合大小的增加而减少，那么这个函数被称为次模的。次模性的另一个定义是$$f(S \cup T) + f(S \cap T) \le f(S) + f(T) \qquad \forall S,T \subseteq A.$$值得注意的是，在一项开创性的贡献中，Asadpour和Nazerzadeh [31]证明了如果\(f^*\)是单调次模的，那么适应度差距被限制在\(\frac{e}{e - 1} \approx 1.58\)以内。在某些情况下，这个界限是紧的。他们还提供了最优的多项式时间非自适应算法，可以将最优（自适应）算法的质量近似到接近\(\frac{e}{e-1}\)的因子。

在我们的探测问题中，我们感兴趣的是选择k次探测来最大化父函数\(f^*\)。显然，如果底层函数f是单调的，那么\(f^*\)也是单调的。与单调性不同，这种继承性对于次模性并不成立。我们在补充材料中展示了对于每一个调节因子加性可分函数\(f_{ras}\)，父函数是次模的。相比之下，覆盖函数\(f_{cov}\)是次模的，但其父函数不是次模的。这激发了以下定义。定义6：如果父函数\(f^*\)是次模的，那么函数f被称为父次模的。

我们的目标是设计好的（非）自适应算法来进行调节因子探测。我们的方法与Asadpour和Nazerzadeh [31]类似，因为我们使用贪婪算法来优化不同的目标函数。我们提供了我们算法的概述，应用它们时面临的挑战，以及一些调整和对它们分析的影响。我们所有算法的伪代码以及关于它们的近似保证的新结果可以在补充材料中找到。

Asadpour和Nazerzadeh [31]分析了基本的自适应（Amp）和非自适应（Namp）算法，这些算法基于之前选择的调节因子集合来短视地选择下一个要探测的候选者。在Amp中，第i步选择的调节因子是那个能够最大化目标函数预期边际增加的调节因子（即，从之前探测到的最佳\(\ell\)个调节因子中选择的那个以及第i步中的那个）。由于Amp是自适应的，它可以考虑之前探测的确切实现，但它依赖于第i步选择的那个的预期值。相反，Namp是这个算法的非自适应版本，在这个版本中，预期边际增加是在所有之前选择的调节因子集合以及第i步中的那个上计算的。

对于\(f_{max}^*\)和\(f_{sum}^*\)，我们在补充材料中指定的一个小预计算/简化之后，直接应用Amp而不做任何修改。对于Namp和\(f^*_{sum}\)（或\(f_{max}^*\)，我们根据预期边值的总和（或最大值）对调节因子进行排序，并分别探测预期值最高的k个调节因子。这简化了Namp原始定义中的边际增加规则。

在应用这些算法进行目标为\(f^*_{cov}\)的调节因子探测时，我们面临多个挑战。值得注意的是，我们观察到\(f^*_{cov}\)不满足次模性。因此，[31]中关于适应度差距的小常数保证并不适用。命题1：\(f_{cov}\)是单调次模的。对于\(k > \ell + 1\)，它不是父次模的。证明：\(f_{cov}\)和\(f_{cov}^*\)的单调性是显而易见的，因为添加到集合S中的每个调节因子都只能改善S当前获得的值。现在设\(S, T \subseteq \mathcal {A}\)。那么，$$\begin{aligned} f_{cov}\text {是次模的} &\Leftrightarrow f_{cov}\left( S \cup T\right) + f_{cov}\left( S \cap T\right) \le f_{cov}\left( S\right) + f_{cov}\left( T\right) \\&\Leftrightarrow \sum _{b \in \mathcal {B}}\max _{a \in S \cup T}w_{a,b} + \sum _{b \in \mathcal {B}}\max _{a \in S \cap T}w_{a,b} \le \sum _{b \in \mathcal {B}}\max _{a \in S}w_{a,b} + \sum _{b \in B}\max _{a \in T}w_{a,b} \\&\Leftrightarrow \sum _{b \in \mathcal {B}}\left( \max _{a \in S \cup T}w_{a,b} + \max _{a \in S \cap T}w_{a,b}\right) \le \sum _{b \in \mathcal {B}}\left( \max _{a \in S}w_{a,b} + \max _{a \in T}w_{a,b}\right) \\&\Leftrightarrow \forall b \in \mathcal {B} :\max _{a \in S \cup T}w_{a,b} + \max _{a \in S \cap T}w_{a,b} \le \max _{a \in S}w_{a,b} + \max _{a \in T}w_{a,b}。\end{aligned}$$现在固定\(b \in \mathcal {B}\)。显然，\(\max _{a \in S}w_{a,b}\)以及\(\max _{a \in T}w_{a,b}\)必须大于或等于\(\max _{a \in S \cap T}w_{a,b}\)。此外，\(\max _{a \in S}w_{a,b}\)和\(\max _{a \in T}w_{a,b}\)中至少有一个必须等于\(\max _{a \in S \cup T}w_{a,b}\)，并且显然\(\max _{a \in S \cup T}w_{a,b} \ge \max _{a \in S \cap T}w_{a,b}\)。因此，最后一个不等式成立，\(f_{cov}\)是次模的。为了证明如果\(k > \ell + 1\)，\(f_{cov}\)不是父次模的，考虑以下问题实例：我们有\(\mathcal {A} :=\{1, 2, 3, 4\}, \mathcal {B} :=\{1, 2, 3, 4\}, \mathcal {V} :=\{0,1\}, k :=4, \ell :=2\)。每个\(D_{a,b}\)都是确定性的，意味着每个概率要么是1要么是0，所以我们总是得到实现$$\(w_{a,b}) :=\left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 \end{array}\)。$$现在考虑子集\(S \subseteq T \subseteq \mathcal {A}, j \in \mathcal {A}\)，其中\(S = \{1,2\}, T = \{1,2,3\}\)且\(j = 4\)。那么，$$\begin{aligned}&f_{cov}^*\text {不是次模的}\\&\Longleftrightarrow \quad f_{cov}^*(T \cup \{j\}) - f_{cov}^*(T)> f_{cov}^*(S \cup \{j\}) - f_{cov}^*(S) \\&\Longleftrightarrow \quad f_{cov}^*(\{1,2,3,4\}) - f_{cov}^*(\{1,2,3\})> f_{cov}^*(\{1,2,4\}) - f_{cov}^*(\{1,2\}) \\&\Longleftrightarrow \quad 4 - 2> 3 - 2 \\&\Longleftrightarrow \quad 2 > 1 \end{aligned}$$由于我们可以轻松地缩放这样的实例，对于每个\(k > \ell + 1\)，\(f_{cov}^*\)不是次模的。此外，如上所述，选择调节因子来最大化\(f_{cov}\)是（加权）MaximumCoverage问题[37, 38]的泛化。这个问题是NP难的，在计算复杂性的标准假设下，无法在\(\frac{e - 1}{e} - o(1)\)的因子内进行近似。因此，即使我们已经完全探测了k个元素（非自适应或自适应地），我们也面临着计算\(\ell\)个元素的最优选择是NP难的问题。此外，在算法的每次迭代中，我们需要评估从扩展的探测集合中选择\(\ell\)个调节因子的预期边际增加。这再次需要解决同一个问题的实例。

我们的（离线）调节因子选择问题可以通过一个通用的贪婪算法实现\(\frac{e - 1}{e}\)的近似比率，该算法用于具有基数约束的次模函数的最大化[38]。我们称这个算法为Greedy。我们的算法依赖于这个程序来估计探测的预期边际收益，以及从探测到的调节因子集合中进行最终选择。更具体地说，在Amp中，我们使用Greedy从我们已经探测到的\(\ell\)个调节因子开始，计算\(f^*_{cov}\)的预期增加。为此，我们遍历所有未探测到的调节因子，将它们视为已探测到的，并将它们所有的相邻边权重设置为各自的预期值，然后对现在“未覆盖”的子图运行Greedy来估计预期边际增加。

修改Namp并不简单，因为原始算法依赖于在每次迭代中重新计算每个位置的概率分布。我们精确计算了前\(\ell\)轮的预期增加，因为我们知道直到这一点所有调节因子也将在最优解中被选中。然而，这种方法在后续步骤中由于计算难度而无法保持高效。在后面的轮次中，我们改为选择具有最高预期相邻边权重总和的下一个调节因子。这种调整可能与Namp的结构有显著不同。不幸的是，由于底层问题的难度，展示这种（或任何其他简单）修改的通用近似因子是非常非平凡的。我们在补充材料中展示了按度参数化的近似因子，但它们严重低估了实际观察到的性能。在实验中，我们发现这个算法的结果与自适应变体的结果非常接近。

为了强调我们的修改，我们将修改后的算法分别称为AmpCov和NampCov。对于\(f_{max}\)和\(f_{sum}\)，我们仍然称这些算法为Amp和Namp，因为我们只是稍微修改了非自适应算法，而没有修改自适应算法。

最后，为了在实验中准确评估这些算法对于\(f_{cov}\)的效果，还需要计算完全揭示的图的最优（离线）调节因子选择。同样，这在多项式时间内是不可行的。相反，我们制定了一个整数程序（IP）（见补充图10），可以正确解决这样的实例。不幸的是，即使是对这个公式的线性规划（LP）松弛，也可以在多项式时间内解决，但对于网络设置中的更大实例来说，这需要太长的时间（见下文）。原始的IP本身对于较小的实例就已经是不可行的。因此，我们使用LP-Relaxation来上界估计最优离线值（从而下界估计我们算法的近似比率）。对于其他情况，我们使用Greedy来近似最优离线值，并将这种应用称为Off。

为了探索上述算法的实用性，我们在随机图的实例上运行了自适应短视算法（Amp/AmpCov）以及非自适应短视算法（Namp/NampCov）。我们关注了两种情况：1. 均匀设置：我们设置\(n_{\mathcal {A}} = n_{\mathcal {B}} = |\mathcal {V}| = 16 \cdot z\)，对于\(1 \le z \le 20\)。边缘权重分布\(D_{a,b}\)是通过在[0, 1]范围内为边缘权重可能取的每个值随机抽取一个权重（然后对所有权重进行归一化）而产生的随机权重分布。2. 网络设置：我们设置\(n_{\mathcal {A}} = 16 \cdot z, n_{\mathcal {B}} = 400 \cdot z, |\mathcal {V}| = 10\)，其中\(1 \le z \le 20\)。边缘权重模拟具有参数\(\lambda \sim [0.5,2.5]\)的泊松分布随机变量。这种设置旨在模拟基因调控网络。对于每个z值，我们生成了100个随机图模型的实例，并在每个图实例上对每对\((k,\ell )\)运行两种算法，其中\(k \in \left\{ \frac{1}{4}n_{\mathcal {A}}, \frac{2}{4}n_{\mathcal {A}}, \frac{3}{4}n_{\mathcal {A}}\right\} , \qquad \ell \in \left\{ \frac{1}{4}k, \frac{2}{4}k, \frac{3}{4}k, \frac{4}{4}k\right\}\)。然后使用通过选择\(\ell\)个调控因子可以达到的最优离线值\(OPT_O\)来评估这些算法，即最优离线算法Opt。对于\(f_{cov}\)，我们在（完整的）均匀设置以及\(n_{\mathcal {A}} \le 160\)的网络设置中使用了LP-Relaxation。对于网络设置中的更大实例，我们使用了离线近似方法Off。算法的伪代码可以在补充材料中找到。

计算资源
所有算法和实验都是用Rust编程语言实现的。代码可以在GitHub上找到：https://github.com/lukasgeis/BipartiteRegulatorProbing。所有实验都在一台配备AMD EPYC 7702P 64核处理器和512GB RAM的服务器上运行。为了加快实验速度，实例生成和评估可以同时使用多达64个线程进行。

构建TF-基因网络
作为这些算法的一个实际应用，我们构建了一个由TFs（调控因子）和基因（位置）作为节点的网络，边缘权重表示TF对基因的调控活性。这种调控活性\(w_{TF,G}\)是通过计算基因中所有含有TF结合位点的REs的比例来估计的：
$$\begin{aligned} w_{TF,G}=\frac{|c_{TF}| }{|C_G|} * x {\left\{ \begin{array}{ll} x=1,\ \text {如果\ \text {p值} \le 0.05 \\ x=0,\ \text {否则} \end{array}\right. }。 \end{aligned}$$
其中\(C_G\)是基因的所有REs，\(c_{TF} \subseteq C_G\)是包含TF预测结合位点的REs的子集。我们要求基因中具有TFBS的REs数量比所有REs的数量更多（在Fisher精确检验中p值\(\le\) 0.05，\(H_a~=~更大\)）。如果不满足这个要求，\(w_{TF,G}\)将变为零。富集检验的列联表如下：
与基因GREs相关的REs 与GREs无关的REs 没有TFBS_{TF}的REs
富集要求的背后原因是为了抵消TFs之间结合位点数量的不平衡。否则，即使对于含有平均数量motif的基因，具有大量motif出现的TFs也会有较高的边缘权重。因此，所有非零边缘权重的边都被认为代表了TF和基因之间的调控相互作用。换句话说，在我们的网络中，一个基因预计会被所有与该基因相连的TFs调控。

我们使用国际人类表观基因组联盟（IHEC）EpiATLAS（https://ihec-epigenomes.org/epiatlas/data/，元数据版本1.1）中的健康T细胞样本（IHECRE00000187）构建了该网络。为了计算边缘权重，需要两种类型的信息：哪些REs调控哪些基因，以及哪些REs具有哪些TF结合位点。作为候选REs，我们使用了H3K27ac ChIP-seq信号（MACS2 narrowPeak文件[43]）中标记的峰值。候选REs与基因之间的相互作用是用改进的Activity-by-Contact模型[36]预测的，使用MACS2信号值作为活性，并使用35个生物样本的平均Hi-C矩阵作为接触数据[44]。排除了已知会积累异常信号的REs重叠区域[45, 46]。RE与转录起始位点（TSS）之间的最大距离被设置为2.5 MB，gABC-score的截止值设置为0.02。相互作用是针对GENCODE v38注释中的所有基因预测的[47]：
此图像的替代文本可能是使用AI生成的。

TFBS的注释使用了FIMO [48]（v.5.4.1），它结合了来自JASPAR 2022 [49]、HOCOMOCO v11 [50]以及[51]的工作（可在https://github.com/SchulzLab/STARE/blob/main/PWMs/2.2/Jaspar_Hocomoco_Kellis_human_transfac.txt）的非冗余motif集合。只有在其启动子（\(\pm ~200~bp\)范围内有H3K27ac ChIP-seq峰值的TFs被考虑在内（\(n = 420\)）。对于TF二聚体，所有组成TF的启动子都必须有峰值。RE的平均碱基含量被用作FIMO的背景。对于基因组坐标的处理，我们使用了pybedtools（v0.9.0）[52, 53]。MyGene.info的API被用来将基因名称与Ensembl IDs匹配[54,55,56]。

我们定义了两组基因，这些基因比所有注释的基因更可能被相似的TFs调控。首先，我们从GSEA（生物过程本体GO:0002456）[57, 58]中选取了所有已知参与T细胞介导的免疫的基因。作为第二组基因，我们收集了DISGENET（curated v24.3）[59]中与以下疾病之一相关的基因：“Leukemia, T Cell”、“human T cell leukemia”、“Early T-Cell Precursor Lymphoblastic Leukemia”、“Lymphoid Leukemias”。我们只能使用在基因注释中找到的基因，因此T细胞介导的免疫有131个基因，淋巴细胞白血病有83个基因。

为了找到覆盖T细胞网络的TFs，我们使用了\(f_{cov}\)，并对每个\(k \in \{30, 40, 50\}\)进行了500次迭代，总共进行了1,500次运行。选定的调控因子数量保持不变，为\(\ell =3\)。换句话说，目标是找到具有最佳网络覆盖率的TF三元组。对于每次运行，都选取了\(f_{cov}\)返回值最高的TF三元组。为了测试其他TF集合大小，我们对每个\(3\le \ell \le 10\)的情况进行了10次迭代，同时设置\(k=50\)。由于对于较大的\(\ell\)，OFF的运行时间不可行，因此没有对其进行测试。

为了测量TF三元组中的共表达，我们使用了IHEC EpiATLAS中的表达矩阵’genes_TPM.csv.gz’，并使用皮尔逊相关系数来评估基因之间的相关性。对于三元组中的每对可能的TF，计算了相关系数，然后取平均值以获得每个三元组的值。对于TF二聚体，首先对组成TF的表达值进行了平均。

在模拟图上的性能
近似因子
我们创建了多个具有不同参数的图实例，以评估Amp/AmpCov和Namp/NampCov在网络和均匀设置中的近似因子，同时改变了目标函数。在具有\(f_{cov}\)的网络设置中，NampCov在所有情况下都几乎达到了与AmpCov相同的近似因子（图2）。只有在\(k = \ell\)的情况下，我们有近似保证时，自适应算法的表现略优于非自适应算法。平均近似因子范围从0.94到0.99。对于较大的\(n_{\mathcal {A}}\)，我们观察到了更好的近似因子，因为单个调控因子对总值的影响减小了。我们还看到了一个明显的趋势，即更多的探测和更少的探测会产生更高的回报，因为我们有更多的机会探测到优秀的调控因子，而不会因为探测到不良调控因子而损失价值。从使用低估近似因子的LP-Relaxation切换到高估近似因子的离线贪婪近似是相当明显的。然而，对于较小的实例，以及\(n_{\mathcal {A}}\)的增加趋势（接近\(n_{\mathcal {A}} = 160\)的步骤），表明无论实例大小如何，我们的算法至少达到了最优离线值的\(94\%\)（这本身是最优自适应策略值的上限）。

对于均匀设置和其他两个目标函数\(f_{max}\)和\(f_{sum}\)（补充图1+2），与\(f_{cov}\)相比，\(f_{max}\)和\(f_{sum}\)的近似比率更高，因为单条边的值空间更小。

图2
此图像的替代文本可能是使用AI生成的。

网络设置中AmpCov/NampCov的近似比率与Off相比，作为\(n_{\mathcal {A}}\)和\((k,\ell ) \in \left\{ \frac{1}{4}n_A, \frac{2}{4}n_A, \frac{3}{4}n_A\right\} \times \left\{ \frac{1}{4}k, \frac{2}{4}k, \frac{3}{4}k, \frac{4}{4}k\right\}\)的函数。垂直线是Off从IP切换到Greedy的截止点。

运行时间
我们比较了算法、目标函数和图设置在模拟图上的运行时间（图3，补充图3–5）。总体而言，对于\(\ell < k\)，NampCov总是比AmpCov快。由于在自适应算法中，每次探测后需要多次运行Nemhauser等人[38]的贪婪算法，因此在\(\ell < k\)且\(f_{cov}\)的网络设置中，AmpCov要慢得多。然而，对于\(\ell = k\)，AmpCov在两种设置中都比Namp快，因为我们不再使用贪婪近似作为子程序。在截止点之后的轻微下降可能是由于在服务器上不同的调度造成的，因为我们在这两部分实验中运行的时间不同。然而，这并不影响我们的观察结果。

对于其他两个目标函数\(f_{max}\)和\(f_{sum}\）（补充图3+4），每个算法的运行时间都不到百分之一秒，这导致时间测量中有很多噪声，这可能进一步被引入的并行性增强了。

图3
此图像的替代文本可能是使用AI生成的。

网络设置中AmpCov/NampCov的执行时间，作为\(n_{\mathcal {A}}\)和\((k,\ell ) \in \left\{ \frac{1}{4}n_A, \frac{2}{4}n_A, \frac{3}{4}n_A\right\} \times \left\{ \frac{1}{4}k, \frac{2}{4}k, \frac{3}{4}k, \frac{4}{4}k\}\)的函数。垂直线是Off从IP切换到Greedy的截止点。

识别覆盖T细胞网络的TFs
对于一个生物学应用，我们在由TFs和基因组成的网络上使用了探测算法，这些数据来源于T细胞样本。为两组基因构建了网络：一组是与T细胞介导的免疫相关的基因，另一组是报告参与淋巴细胞白血病的基因。在我们的设置中，重复运行算法以选择三个TF的集合来最大化目标函数\(f_{cov}\)。我们还测试了更大的集合，但主要关注TF三元组。正如预期的那样，最常选择的TF三元组比随机组装的TF三元组调控了更多的网络基因（图4a，补充图6a）。然而，对于一部分基因，没有选择到调控TF。特别是那些只与很少TF相连的基因很少被包括在内，这表明即使在预定义的基因集合中，调控程序也是多样的，因此仅用三个TF很难捕捉到。通过增加TF集合的大小，可以覆盖更多的基因，但即使是大小为10的集合也只覆盖了\(53\%\)的T细胞介导的免疫基因和\(65\%\)的淋巴细胞白血病基因（补充图6b+c）。选择的大小为三的集合主要由相似的TF组成，在排名靠前的TF三元组之间变化很小（补充图6d, e）。最常选择的TF三元组在Off和NampCov之间相似，而AmpCov选择了一些不同的TF三元组（图4b, c）。对于k（可以探测的TF数量）的三种测试选择，结果相似，特别是AmpCov和NampCov产生了几乎相同的结果（补充图6f）。我们还测试了网络构建的不同阈值，即当在基因和TF之间绘制边时的阈值。正如预期的那样，更宽松的阈值导致更多基因被TF靶向（补充图9a–d）。由于阈值的变化改变了TF在连接基因数量方面的排名（补充图9b, e），探测算法的优先级也相应地改变了（补充图9c+f）。然而，除非移除了太多TF-基因边，否则结果保持稳定。

图4
此图像的替代文本可能是使用AI生成的。

全尺寸图像

应用算法在T细胞网络中寻找转录因子（TFs）。a. NampCov覆盖的基因在T细胞介导的免疫网络中的覆盖率，与三个最常选择的TF三元组相比，以及三个随机选择的TF三元组。从排名前十的转录因子（TF）中选取了最佳的三个TF组合，以减少单个TF的冗余性。每一列代表网络中的一个基因。第一行显示了与该基因有非零边连接的TF数量，即那些在基因的调控元件（REs）中转录因子结合位点（TFBS）富集的TF（p值≤0.05）。基因的调控元件中TFBS的比例代表了网络的边权重。关于淋巴细胞白血病网络的覆盖情况以及完整网络的可视化，请参见补充图6a和补充图7+8。b、c展示了三种算法选出的十个最常见TF组合的重叠情况。d、e显示了所有被TF调控的基因组范围基因的百分比（背景调控）与相应基因集中被调控基因的百分比。仅考虑了至少有一个RE的基因。颜色表示NampCov选择该TF的次数。根据NampCov选出的TF组合中的出现次数，对前五个TF进行了排名，并附上了支持这些TF在T细胞介导的免疫（f）和淋巴细胞白血病（g）中功能的文献参考。

为了探讨TF组合内的关系，我们基于IHEC EpiATLAS的RNA-seq数据（1555个样本）研究了TF之间的共表达情况。对于所有算法和两个网络，最常被选出的TF组合中的TF之间的共表达显著高于随机TF组合（p值≤0.05，双侧Mann–Whitney U检验，补充图6g、h）。这表明这些TF组合中的TF具有相似的调控程序。

为了评估单个TF，我们根据它们在选出的TF组合中的出现频率对它们进行了排名。我们重点关注了NampCov算法，因为该算法在模拟网络上仍然接近最优解的同时效率最高。我们通过比较TF在网络中连接的基因集与基因组范围背景中的基因百分比，来评估最常出现的TF的特异性。背景包括所有至少有一个RE的注释基因。对于两个测试的网络，NampCov选择的TF调节的基因集成员相对于背景更多（图4d+e）。这种效应是数据驱动的，因为算法仅在特定基因集的网络上运行，并没有使用任何关于背景的知识。尽管如此，它表明了优先考虑的TF的特异性，并暗示了基因集的不同调控模块。对于许多选出的TF，它们在T细胞介导的免疫和淋巴细胞白血病中的参与已经在文献中有所报道（图4f、g）。一个经常被选出的TF是ZBTB17，除了其对早期胚胎发育的重要性外，还发现它参与了淋巴细胞的发育和T细胞表型的指定[60, 61]。IRF1被报道对于特定T细胞亚型的分化和功能是必需的[62,63,64]。另一个经常被选出的与T细胞介导的免疫相关的TF是PATZ1，多项研究表明它在不同T细胞群体的发育中起重要作用[61, 65, 66]。从TCGA图谱的淋巴细胞白血病基因网络中反复选出的多个TF在各种癌症样本中表现出差异表达[67]，包括ZNF354A、ZNF530和ZNF320。ZNF320还被发现与肝细胞癌中的免疫浸润有关[68]。ZSCAN22被预测是与急性髓系白血病患者生存相关的micro RNA的目标[69]。另一个经常被选出的TF是KLF10，它被反复发现是多种癌症类型的肿瘤抑制因子[70,71,72]，并调节T细胞亚型的分化[73]。对于一些报告的TF，文献支持较少或相对间接，这可能是因为许多TF研究不足，因此缺乏已知的功能。

总体而言，我们能够将这种方法应用于生物数据集，并找到了可能调节T细胞中具有特定功能的基因集的有意义的TF组合。

讨论与结论
我们提出了在调控网络中找到一组TF的挑战，这些TF能够最大化对受影响基因的影响。我们评估了不同的目标函数，并在不同的情况下展示了不同的近似保证。在模拟图上，我们发现非自适应贪婪算法的性能与自适应贪婪算法相似，而且通常计算速度更快。虽然我们发现像\(f_{max}\)和\(f_{sum}\)这样的简单目标函数容易满足Asadpour和Nazerzadeh [31]的近似保证，但更复杂和现实的目标函数如\(f_{cov}\)则未能达到必要的要求。即使对于规模合理的实例，使用LP-Relaxation解决离线问题也很快变得不可行。因此，我们使用了在实践中表现非常好的修改后的贪婪算法来处理\(f_{cov}\)。我们还证明了基于最大度的算法未修改版本的性能保证。更具体地说，我们展示了将基因视为独立的通用目标函数在应用Asadpour和Nazerzadeh的算法时，会导致与最大基因节点度成比例的损失。对于\(f_{cov}\)，我们还展示了这种损失与网络中的全局最大度（TF或基因）成正比。

在将算法应用于基于T细胞数据构建的网络时，我们能够找到在调节相应基因集中具有已知功能的有意义的TF组合。与其他检查TF-基因关系的方法不同，我们考虑了基因的所有潜在REs，而不仅仅局限于启动子区域。关于网络的构建，需要考虑一些注意事项。我们使用位置权重矩阵来找到TFBS，这是一个简化的TF结合模型，它假设各个结合位置是独立的。此外，TFBS被二值化了，这忽略了低亲和力结合位点的潜在影响[2]。基于基序的TFBS预测也引入了比较具有冗余结合偏好的TF的挑战。此外，我们通过TF结合的REs的比例来建模TF和基因之间的调控强度，要求结合的REs必须富集。这种调控强度的估计不太可能完全捕捉TF-基因关系，因为TF结合的频率并不期望线性转化为表达调控。可能还涉及许多其他未考虑的因素，如共因子的存在、TF之间的不同合作模式、翻译后修饰或染色质修饰[2, 74, 75]。重要的是，本工作中使用的网络作为应用探测算法的示例。它们可以被用任何可用的方法或数据库构建的其他TF-基因网络替换[4,5,6,7,76]。唯一的要求是它们被构建为二分图，即TF和基因之间只有边，并且边权重是正的。在设计扰动实验的背景下，我们选择的TF子集是针对特定基因集优化的，而不需要了解整个基因组网络。因此，TF不必特定于基因集，但也可能调节其他基因。尽管如此，在我们的应用中，选出的TF对特定基因集表现出特异性，这可能是由于具有相似功能的基因更可能被类似的TF靶向。另一个需要考虑的点是，我们设置中的所有基因都被赋予了相同的权重。在某些情况下，可能希望不同地加权基因，例如对于路径中的基因，位于中心位置的基因应该比下游的基因更重要。这可以通过根据基因的重要性来调整边权重来实现。

我们方法的最小数据要求，即RE活性的测量，使其能够轻松应用于其他数据。识别在TF-基因网络中覆盖率最高的TF集合可以作为寻找扰动研究候选者的起点。TF集合的大小是灵活的，可以根据手头的问题来选择。未来的一个有趣方向是将这些算法应用于其他类型的网络，例如micro RNA与其mRNA目标之间的相互作用，或配体-受体网络。