TRANSWESD：用传递约简推断细胞网络

摘要

动机：区分直接影响和间接影响是生物网络逆向工程中的一个核心问题，因为它有助于检测和去除假阳性边缘。传递约简是一种消除反映间接影响的边的方法，但在重建具有真正冗余结构的循环交互图时，它的使用存在问题。

结果：我们提出了TRANSWESD，这是一种针对加权有符号有向图的TRANSive约简的详细变体，它克服了现有版本的概念问题。主要变化和改进涉及：（i）从系统扰动实验生成高质量扰动图的新统计方法；（ii）使用边缘权重（关联强度）识别真正的冗余结构；（iii）周期的因果解释；（iv）放宽了传递归约的定义；（v）大型网络的近似算法。通过使用标准化的基准测试，我们证明了我们的方法优于现有的传递约简变体，并且尽管概念简单，但与其他逆向工程方法相比具有很强的竞争力。

联系人： klamt@mpi-magdeburg.mpg.de

补充信息： 补充数据可在生物信息学在线。

1简介

蜂窝网络的逆向工程已经成为分析和利用由组学技术生成的越来越多数据的关键方法（Gardner和Faith，2005；赫克等。,2009；Markowetz和Spang，2007). 然而，代谢反应网络的结构可以主要根据基因组信息重建，对于许多生物体来说，这是非常详细的（Oberhardt等。,2009)在许多情况下，对调控和信号转导网络拓扑结构的了解仍然不完整，即使是“典型信号通路”的接线图在不同的细胞系中也可能不同（Saez-Rodriguez等。,2009). 反向工程方法的最终目标是通过分析系统和受控扰动实验的数据来识别相关参与者（基因、蛋白质等）之间的相互作用。其结果是一个网络，在许多情况下表示为一个图，它可以是有向的或无向的，并且在其边缘可能有符号和/或权重。一些算法提供了精细表示，如布尔网络（Akutsu等。,2003；萨伊兹·罗德里格斯等。,2009)反应网络（Durzinsky等。,2008)或微分方程（Nelander等。,2008)但主要的结果仍然是底层的网络拓扑。

一种简单而智能的逆向工程方法基于过渡还原，一种图论方法（Aho等。,1972)瓦格纳首先认识到了其重建监管网络的潜力(2001). 基本思路如下：重建监管网络n个节点，首先测量野生类型中节点的状态，然后至少执行n个微扰实验：实验中我节点我感到不安，而所有其他n个−1个节点与野生型节点相比，是否改变了其状态。如果发生扰动我影响j个，来自节点的定向边我到j个，表示为我→j个绘制。在所有微扰实验中观察到的这些效应的完整集合产生了扰动图摄动图中的每条边反映了一个节点对另一个节点的直接或间接影响。下一步处理网络重建中的一个中心问题，即识别和删除表示间接影响的边缘。Wagner使用的过渡还原(2001)目的是找到能够解释实验中所有效应的最小（最简约）子图。最一般形式的传递约简允许删除和添加边以找到最小图（Aho等。,1972). 然而，在网络重建的背景下，人们通常关注的是只能删除边缘的特殊情况，即搜索最小值子图解释扰动图[也称为最小等价图问题（伯曼等。,2009；莫尔斯和汤普森，1969)]. 其中，我们只考虑基于边缘去除的传递约简。瓦格纳(2001)通过去除所有边从扰动图中确定最小子图我→j个从中开始的（简单）路径我结束于j个（未使用我→j个)可以找到，假设我在j个是间接的，因此可以通过路径来解释。结果图是可传递性约简扰动图的。下面描述了一个简单的示例图1a.每个非循环图都有一个唯一的传递约简（具有最小数量的边；Aho等。,1972)解释所有测量的扰动效应。

Wagner提出的方法(2001)很容易实现，但也有一些缺点，这可能是它很少应用的原因。首先，如上所述的传递性约简并没有考虑摄动实验的全部信息，即使只考虑定性观测。如果一个节点对扰动显示出显著的响应，那么至少可以将测量到的效果分为“向上”或“向下”。通过向扰动图中的每条边添加一个符号标签，可以考虑到该信息，扰动图变成有符号有向图（请参见图1b-a签名版本图1a） ●●●●。然后可以以类似的方式执行传递性还原：边缘我→j个仅当存在来自的路径时才删除我到j个其整体符号（相关边符号的乘积）与该边的符号相对应。如中的示例所示图1b、 这可能会保存未签名版本中错误删除的边。原始传递约简方法的第二个缺点是，即使在有符号扰动图中，也存在删除真边的风险。传递约简的激进剪枝策略旨在最小化重构网络中的假阳性边缘，但它可能会导致大量的假阴性。这种效应在由许多（相干的）前馈环路组成的网络中变得明显，其中一个节点可能通过相同符号的直接（边缘）和间接（路径）链接影响另一个节点。由于前馈环路已被证明在基因调控网络中频繁发生（Shen-Orr等。,2002)，此属性可能会成为该方法的严重限制。第三个缺点是，扰动图是非循环的前提条件，这一条件在现实的生物网络中通常不满足。如果扰动图是循环的，传递约简的解通常不是唯一的。正如我们将看到的，符号扰动图中的负圈可能会给传递约简带来更多的复杂性。

在本文中，我们将介绍TRANSWESD（加权有符号有向图中的传递归约），它是传递归约的一个新变体，旨在克服这些问题。其他作者（Albert等。,2007；特里西等。,2007). 然而，我们的方法结合并扩展了现有的变体，并且在几个关键方面有所不同（例如，使用加权扰动图、负循环的处理以及初等与非初等路径的处理）。我们还讨论了与识别显著扰动效应相关的问题，这是生成扰动图的关键步骤。通过使用标准化的基准测试，我们证明了我们的方法优于现有的传递约简变体，并且尽管概念简单，但与其他逆向工程方法相比具有很强的竞争力。

2方法

2.1定义

我们从图论中总结了一些将在本文中使用的标准术语和符号。A类图形G= (V（V）,E类)由一组V（V）节点（或顶点）和集合E类连接成对节点的边。在我们的例子中，节点可能代表基因、mRNA、蛋白质等，而边缘相应地代表物理或有影响的节点间相互作用。我们只关心定向的图形(二合字母)边是定向的（也称为弧），即。e（电子）∈E类是有序对e（电子）= (u个,v（v）)不同节点的u个,v（v）∈V（V），也表示为u个→v（v），其中u个是开始节点v（v）结束节点。A类有符号有向图G= (V（V）,E类，φ）还包含一个符号映射φ：E类→ {−，+}表示每条边的开始节点对其结束节点是否有促进或抑制作用。定向边缘(u个,v（v）)带有标志秒表示为u个→_秒v（v）.A型加权有符号有向图G= (V（V）,E类，φ，γ）还包含一个权重映射γ：E类→ ℜ_≥0它为每条边分配一个权重，我们在这里假设它是非负的。边缘(u个,v（v）)带有标志秒和重量w个表示为u个→_秒,w个v（v）.

A类步行有向图中是节点和边的交替序列v（v）₀,e（电子）₁,v（v）₁,e（电子）₂,…e（电子）_n个,v（v）_n个以节点开始和结束v（v）₀和v（v）_n个，分别满足节点v（v）_我−1和v（v）_我由边缘连接e（电子）_我。在我们的术语中路径是具有附加条件的行走，即没有节点出现两次，即路径不包含循环。后一个属性有时通过将路径称为“简单”或“基本”来强调，行走有时被称为非基本路径。最后，一个（基本）周期是一个没有重复节点的封闭行走，除了第一个和最后一个节点重合。具体的路径或循环被写成节点和箭头（边）的交替序列，例如。u个→v（v）→w个，如果所涉及的节点之间不存在平行边，则它为路径或循环提供唯一标识符。具有起始节点的路径u个和结束节点v（v）表示为u个⇒v（v）也可以由单个边缘组成u个→v（v）.

这个长度路径/循环的权重是根据相关边的权重计算出来的，例如，将它们相加（∑-公制）。我们还需要另一个变量，称为MAX-metric，其中路径的长度是其所有边的最大权重。路径/循环的符号是通过边缘符号相乘获得的（因此，有符号有向图不等价于具有正负边缘权重的加权有向图）。来自的路径u个到v（v）带有整体标志秒表示为u个⇒_秒v（v）。如果此路径有长度d日然后我们写u个⇒_秒,d日v（v）.

2.2整个过程的工作流程

我们首先概述我们的算法。

第1步：如引言所述，起点是一个野生型实验加n个每个扰动实验中n个节点受到扰动，并测量其他节点的响应，无论是在瞬态还是稳态（如果没有另行说明，我们假设后者）。我们用表示野生型状态x个⁰(x个_我⁰表示的野生类型状态我-th节点），并假设每个物种的测量值均归一化为各个物种所观察到的最大值，即。x个_我⁰∈ [0, 1]. 在k个-第次扰动实验（其中节点k个扰动）表示为x个^k个，即。x个_我^k个是的状态我-实验中的第th个节点k个。使用第页^k个∈{−1，+1}我们表示插管是否k个强迫节点数量/活动减少（-1）（例如通过淘汰或淘汰）或增加（+1）（例如过度表达）k个.

第2步(第2.3节)：对于每个节点，我们比较未扰动状态(x个_我⁰)扰动实验中的测量状态(x个_我^k个). 使用适当的阈值策略，显著的变化被识别并包括为有符号边k个→我在得到的扰动图中。

步骤3(第2.4节)：扰动图中的每个识别边都有一个从测量值中提取的权重，表示两个连接节点之间的关联/交互强度。

第4步(第2.5节)：最后一步是使用我们新的TRANSWESD算法计算传递性约简，该算法可以处理也可能包含循环的加权有符号有向图。（请注意，原则上，TRANSWESD可以接受任何扰动图，即使生成图的方式与步骤1-3不同。）

2.3生成扰动图：阈值

引入阈值ϑ来设置搜索边时所需的总体最小变化量。对于所有节点对，它都保持不变。第二阈值β_我考虑了每个节点的单独动态特性，并从节点的方差中导出我不包括节点扰动的整个扰动剖面我自身。

因此，我们从节点引入边k个到我仅当两个条件（i）|Δ_我^k个|>ϑ和（ii）|Δ_我^k个| > β_我，满足。我们计算β_我:Γσ_我具有方差标度因子Γ和SDσ_我基因的x个_我注意，取决于ϑ和Γ的选定值以及节点的波动x个_我，我们要么有ϑ≥β_我或ϑ<β_我.基准测试(第3节)表明条件（i）或（ii）单独导致预测性能较弱（结果未显示）。仅使用条件（i）会忽略单个节点的动力学。例如，由于其他节点的抑制而具有较小绝对变化的节点边缘可能会丢失。仅条件（ii）容易出现稠密图中的实验数据错误，这增加了测量噪声的概率。参数ϑ和Γ的幅值可以根据已知相互作用图的扰动数据进行估计，其功能应接近所研究的系统。或者，如果噪声分布函数已知，则可以直接计算给定阈值P（P）-价值。

从含噪数据中获取高质量扰动图的合适阈值策略是一个重要步骤，因为每个图都有一个临界边密度，传递约简相关算法在修剪结果和计算时间方面都能很好地达到该临界边密度。虽然反射间接效果的边缘可能会在稍后阶段被TRANSWESD过滤（见下文），但既不表示直接影响也不表示间接影响（因此噪声）的边缘无法纠正，并将导致重建错误。另一方面，FN的数量也应尽量减少，因为它们无法通过传递还原恢复。示例中给出了我们的阈值策略补充材料我们注意到，我们的方法与叶提出的噪声学习模型有一些相似之处等。(2010)用于过滤非确定性影响。主要区别在于，我们的方法不假设某种噪声分布函数。

2.4量化关联强度

对于传递约简的变体，我们需要为扰动图中的符号边分配权重，以量化有向关系的强度。因此，对于每个有序的节点对(u个,v（v）)，我们确定两两条件相关（ρ_{u个,v（v）})来自u个-th和v（v）-被测状态向量的第个元素x个⁰,x个¹,x个²,…,x个^v（v）−1,x个^v（v）+1,…,x个^n个（参见Rice等。,2005). ρ_{u个,v（v）}计算为线性相关系数，其中我们从v（v）-因为我们想要ρ_{u个,v（v）}量化v（v）在u个和外部扰动v（v）无法解释u个相应地，ρ_{u个,v（v）}不是对称的。

对于每个边缘u个→_秒v（v）摄动图中捕获的P（P）派生于第2.3节，我们将其权重指定为1−|ρ_{u个,v（v）}|即权重越小关联越高。这种加权方案，其中边缘权重表示两个节点行为之间的“距离”，与其他权重较大通常表示高度关联的工作有些相反。然而，我们需要这种表示，因为我们将使用最短路径计算来找到具有最高总体关联性（最低总体权重）的路径。

有时会这样秒≠sgn（ρ_{u个,v（v）})表明v（v）扰动时u个不反映相关性分析得出的符号。尽管如此，保持这一边缘似乎很有用，但权重较高（接近最大可能权重1）表明相关性较弱。我们还测试了其他权重，例如v（v）当扰动时u个, |x个_v（v）^u个−x个_v（v）⁰|，但事实证明，该算法对ρ的性能更好_{u个,v（v）}，只是因为它评估了许多(n个)实验。另一方面，用于确定边缘u个→_秒v（v）存在于P（P）（为了固定其符号）u个导入节点v（v）似乎比条件相关性更适合(第2.3节).

2.5使用TRANSWESD进行过渡还原

在这个阶段，我们得到了一个加权的、有符号的和有向的扰动图P（P）= (V（V）,E类，φ，γ），其中我们假设其边显示直接或间接关系。保持边缘捕获间接影响将导致FP预测，传递性约简试图删除FP以获得真阴性（TN），但存在删除真阳性（TP）预测导致FN的风险。从简单的扰动图开始，以循环扰动图结束，我们逐步推广了传递约简的思想，并解释了我们在TRANSWESD中的扩展，这些扩展寻求最小化先前变量的缺点。

2.5.1有符号非循环图的传递约简

瓦格纳(2001)使用传递约简来修剪无符号非循环扰动图。将此过程推广到有符号非循环扰动图是很简单的P（P）（此时我们忽略了重量）。基本思想是检查每条边u个→_秒v（v）在里面P（P）是否存在基本路径u个⇒_秒v（v）不涉及该边缘，这可以被视为对观察到的影响的解释u个→_秒v（v）让一个人来移除这个边缘。为此，在第一步中，我们计算每对节点(u个,v（v）)最短的正路径和最短的负路径告诉我们从u个到v（v）根本不存在。由于我们只对路径的存在感兴趣，因此可以使用任意边权重，例如将所有边权重设置为一，以及任意度量。我们可以使用双标签算法，这是Dijkstra算法的推广版本，用于计算∑度量中的最短正/负路径。如果有符号图是非循环的，它可以在多项式时间内提供精确的结果（Hansen，1984；克拉姆特和冯·坎普，2009). 我们将最短正负路径的长度存储在矩阵中S公司⁺和S公司⁻分别是。例如，S公司⁺(u个,v（v）)存储最短正路径的长度u个到v（v）.无限长(inf公司)如果不存在路径，则存储。

在非循环有符号图中，借助于S公司⁺和S公司⁻：我们检查每个边缘u个→_秒v（v）我们能否找到继任者z（z）≠v（v）属于u个这样的边缘u个→_q个z（z）和一条路z（z）⇒_t吨v（v）满足符号条件的存在q个·t吨=秒（如果S公司^t吨(z（z）,v（v）) <inf公司). 如果是这样，我们可以得出以下结论：u个→_秒v（v）可以用扩展路径来解释u个→_q个z（z）⇒_t吨v（v），这是基本的，因为我们有一个非循环图。因此，我们删除u个→_秒v（v）然后继续下一条边。请注意，没有必要重新计算最短路径长度S公司⁺和S公司⁻边缘移除后u个→_秒v（v）：在使用此边的所有路径中，我们可以将后者替换为u个→_q个z（z）⇒_t吨v（v）因为，在非循环图中，可以确保生成的路径仍然是基本的，因此是有效的解释。消除所有可移动边，得到唯一的最小等价图P（P）_信托收据产生与原始图形相同的扰动效果P（P）。无符号图中的传递约简使用相同的算法，但忽略了符号条件。

我们认为条件(3)可以放松到(2)，因为在我们应用传递约简时，没有必要保留路径u个⇒_t吨v（v）两个节点之间u个和v（v）如果没有边缘u个→_t吨v（v）（即既不是直接影响也不是间接影响u个在v（v）)可以从实验中推断出来。然而，只要我们考虑非循环图，两种定义都会导致相同的结果，因为(3)以下为(2).

中的示例图1b表明，考虑边缘标志可以避免删除无法解释的边缘：与之相反图1a（无符号扰动图）边a→₋保留B是因为路径A→₊C类→₊B无法解释这条边的负号。

2.5.2有符号图和加权非循环图中的传递约简

为了实现(4)，传递约简步骤必须修改如下：我们删除一条边u个→_秒,w个v（v）如果我们能找到继任者z（z）≠v（v）属于u个这样的边缘u个→_{q个,c（c）}z（z）和一条路z（z）⇒_t吨,d日>v（v）满足符号条件的存在q个·t吨=秒现在再加上重量条件max(c（c）,d日) < α ·w个.正因子α控制路径必须具有的整体关联强度，以解释给定边。通常情况下，人们会选择一个接近于1的值（我们使用0.95），但人们也可能更喜欢较小的值，这就要求在路径的所有边上都有明显更大的关联来解释边。在极端情况下α=0，我们有P（P）_信托收据=P（P）如果α>1，则允许路径中的边与路径解释的边具有更低的关联。带有α=inf公司条件(4)符合条件(2)因此，我们回到了未加权图的传递性约简。同样，具有非循环图首先确保了扩展路径u个→_{q个,c（c）}z（z）⇒_t吨,d日v（v）顺从的u个⇒_{秒，最大值(c（c）,d日)}v（v）是基本的，即。z（z）⇒_t吨,d日v（v）不包含边u个→_{q个,c（c）}z（z），因此是对u个→_秒,w个v（v）其次，我们不需要重新计算S公司⁺和S公司⁻去除边缘后。因此，与上一节类似，如果条件(4)对于移除的边，也将对所有其他边执行此操作。

图1c证明，如果替代路径不能解释其高关联强度，则边缘保持不变。与…对比图1b、 A类→_−,0.3保留D是因为路径A→_−,0.6B→_+，0.5D的长度为0.6，因此在选择α<1时不是一个有效的解释（但在α>2时可能是这样）。

我们注意到Rice中提出的三角形约简方案等。(2005)使用本节所述程序的类似版本；然而，该方案仅适用于三角形，即边u个→v（v）仅当两条连续边u个→z（z）→v（v）解释一下。

2.5.3有符号和加权循环图中的传递约简

我们现在讨论扰动图可能包含循环的最一般情况，这是许多蜂窝网络的固有特性。反馈不仅会导致复杂的动态网络行为，还会阻碍因果关系的推断。因此，传递约简变得更加复杂也就不足为奇了，不仅在算法的结构方面，而且在计算复杂性方面。

与非循环情况一样，我们的过程TRANSWESD从计算最短路径长度开始S公司⁺和S公司⁻这里，我们面临一个内在的算法问题：在含有负圈的图中，这个问题对于基本路径来说是NP-完全的（Lapaugh和Papadimitriou，1984). 幸运的是，人们可以用较低的计算需求检查负循环是否存在。如果没有，我们可以再次使用双标签算法在多项式时间内计算精确结果。事实证明，即使存在负循环，通过使用深度优先搜索或其特殊变体（Klamt和von Kamp，2009). 后一篇文章还描述了在大规模网络中生成合理近似值的多项式算法。

第二个技术问题涉及负循环因果关系的解释。在图1d、 我们看到一个包含负循环C的扰动图的小例子→_+,0.3D类→_−,0.4C.关键问题是我们是否考虑负非初等路径（行走）A→_+,0.2C类→_+,0.3D类→_−,0.4C类→_+,0.3D类→_+,0.3B作为对负面影响a的有效解释→_−,0.6当扰动A时，我们观察到。当α<1时，这个行走的符号和长度实际上允许这样做。特里西等。(2007)认为步行是可能的解释，尽管阿尔伯特等。(2007)没有考虑权重，他们的方法也是基于这种解释的。这带来了一个优点，即只需要计算最短的正/负行走，这在计算上很容易[例如，通过一种改进的Floyd-Warshall算法（Albert等。,2007；特里西等。,2007)]与最短的相比初级的路径。然而，我们认为应该保持A和B之间的负边，原因如下：我们假设网络在A中受扰动时处于稳定状态（在不损失通用性的情况下，我们假设A中存在过表达）。扰动图中的负边图1d表示我们测量到B的激活水平降低。从系统理论（Maurya等。,2003)，可以证明，如果我们测量初始响应或稳态响应在B中，网络中的初始响应由基本路径的符号控制，并且由于边缘a的移除→_−,0.6B意味着只有从a到B的正基本路径仍然存在，B中的初始响应将是正的（简单地说，当观察B中的最初响应时，正路径的影响不能被该路径所诱导的负反馈的影响所超越）。同样在稳定状态下，如果从a到B的负边缘被移除，B不会表现出活性降低（与未受干扰的野生型相比）。如果只有从A到B的正基本路径存在，则在稳态下，负反馈只能与包括正反馈在内的其他结构要求一起诱导相反的效果（Maurya等。,2003). 尽管B中的负效应可能会暂时观察到，但我们通常认为包含负循环的非初等路径不足以解释边；只接受具有适当符号和权重的基本路径。因此，从A到B的负边保持在图1d。

循环图中可能出现的第三个问题是非唯一性。我们的方法的优点是边权重消除了许多可能的非唯一性来源，特别是那些与正循环有关的来源。图1e表示包含正循环的未加权扰动图。从A到B的正边可以用正路径A来解释→₊C类→₊B.另一方面，从A到C的正边可以用正路径A来解释→₊B→₊C.基于Albert中未加权扰动图的方法等。(2007)因此，将删除其中一条边并保留另一条边。选择取决于边缘处理顺序。有了关于关联强度（边权重）的附加信息，通常可以找到α<1的唯一解决方案，如所示图1f： 我们将去除从A到B的边缘，这可以通过从A到C经由B的正路径来解释。

然而，即使使用边权重，也可能出现非唯一性，如所示图1g.在第一步中，我们可以删除边缘a→_+,0.8C（α=0.95，可通过路径A解释→_+,0.5B类→_+,0.6C或A→_+,0.35D类→_+,0.4B→_+,0.6C） ●●●●。在第二步中，我们可以删除边a→_+,0.5B（可由A解释→_+,0.35D类→_+,0.4B） 或边缘D→_+，0.4B（可由D解释→_−,0.3E类→_+,0.2C类→_−,0.25B） ●●●●。我们只能删除其中一个，然后必须停止修剪，否则无法解释删除的边A→_+,0.8C将留在网络中，从而违反条件(4). 因此，我们可能会得到重建图的两个可能的最小解。通常，只有当对于给定的边存在至少两个解释路径，并且网络包含负循环时，才会发生这种情况。在我们的算法中，我们使用贪婪策略，即在每次迭代中，我们尝试删除满足条件的权重最大（关联强度最低）的可解释边(4).

因此，我们进行如下操作：计算后S公司⁺和S公司⁻，我们使用这些矩阵来检测潜在的可解释边。潜在的可解释边缘u个→_秒,w个v（v）是我们可以找到继任者的地方z（z）≠v（v）属于u个这样的边缘u个→_{q个,c（c）}z（z）和一条路z（z）→_t吨,d日v（v）满足符号条件的存在q个·t吨=秒以及最大重量条件(c（c）,d日) < α ·w个与非循环网络相比，可能会出现扩展路径u个→_{q个,c（c）}z（z）→_t吨,d日v（v）不是基本的，因为路径z（z）→_t吨,d日v（v）可能会溢出u个从而在中引入循环u个。移除该边缘后重新计算路径长度时，将看到候选边缘是否真的可以解释（见下文）。所有可能解释的边都是根据它们的权重排序的（首先是最高的），现在一个边以降序迭代这些边。因此，在图1g、 我们将首先删除A→_+,0.8C然后A→_+,0.5B，我们必须保留D→_+,0.4B.这个例子还说明了我们在循环扰动图中必须考虑的第四个问题：原则上，D→_+,0.4B可以用D来解释→_−,0.3E类→_+,0.2C类→_−，0.25B.然而，如前所述，不允许移除，因为这样会影响A→_+,0.8C最初包含在P（P）在修剪后的图和条件中无法解释(4)会被侵犯。因此，在循环图中，如果一条边可以用一条路径来解释，我们就不能期望在删除这条边时所有其他的基本路径都保持不变。在我们的示例中，原始边A→_+,0.8C已经不能用路径A来解释了→_+,0.35D类→_+,0.4B→_+,0.6C如果我们移除边缘D→_+,0.4B.后者的解释，D→_−,0.3E类→_+,0.2C类→_−,0.25B、 无法嵌入路径A中→_+,0.35D类→_+，0.4B→_+,0.6C，因为生成的路径不是基本路径。同样，负循环导致如此复杂的结构。因此，在切割可解释的边缘之前，我们必须检查移除后的最短路径长度S公司_我⁺和S公司_我⁻仍满足条件(4)在生成的中间图中P（P）_我因此，我们必须重新计算最短的有符号路径。由于这是整个算法中最耗时的部分，我们可以尝试简化这一步，例如通过快速计算上述最短路径长度的近似值。此外，可以完全忽略重新计算步骤（即使用原始S公司⁺和S公司⁻并在整个过程结束后检查条件(4)在剪枝图中被违反。事实上，正如我们将在第4节，在许多实际应用中，如果α<1且(4)因此成立。如果不是，则可以接受少量错误或重新引入节点对的边(u个,v（v）)违反(4).

使用精确的算法，得到的修剪图满足条件(4)但不一定是唯一的，也不一定是最小的。然而，在补充材料证明了剪枝图对于排序的边权重列表上的偏序通常是唯一的和最小的。此外，TRANSWESD还包括前面章节中讨论的非循环加权/未加权扰动图的特殊情况：如果没有（实验推导的）权重可用，只需将所有边权重设置为1和α=inf公司在这些非循环情况下，解决方案对于传递约简的原始定义中的边总数也是最小的。

TRANSWESD算法的伪代码版本在补充材料一个实现已经作为API函数集成在我们的MATLAB工具箱中CellNet分析仪（克拉姆特等。,2007).

3结果

使用本文提出的已开发的逆向工程方法，我们参加了逆向工程评估与方法对话（DREAM4）的第四项挑战生物信息学基因网络重建。DREAM倡议提供了一个平台，根据生物信息学数据为高通量基因表达谱分析和基因调控网络重建提供了现实场景（马尔巴赫等。,2009；斯托洛维茨基等。,2007,2009). 从DREAM4挑战中，我们展示了我们的数据集方法的结果Insilico_Size_100号子挑战，可从DREAM网站下载(http://wiki.c2b2.columbia.edu/dream/index.php/The_dream_Project网站). 基于从基因网络中抽取的100个节点的5个子网络大肠杆菌生成了具有随机选择参数的真实酵母动力学模型，并用GeneNetWeaver（Marbach等。,2009)使用随机微分方程。为了重建这些网络，生物信息学测量数据包括野生型和单基因敲除和敲除实验的噪声稳态mRNA表达水平以及时间序列数据。五个网络的黄金标准是在公布所有提交结果后提供的，因此我们可以将计算结果与实际网络进行比较。

对于每个网络，我们首先使用野生型和淘汰稳态数据生成如上所述的扰动图。所需的两个参数是从DREAM3挑战中训练出来的。边缘权重根据淘汰和淘汰数据计算为条件相关系数。仅使用淘汰数据时，结果非常相似。提供的时间序列数据根本没有使用。然后，我们使用TRANSWESD对生成的扰动图进行传递约简，得到最终的重构图。根据性能分析所需的权重对发现的边缘进行排序[与AUROC（接收器-操作员特征曲线下的面积）和AUPR（精度-召回曲线下的区域）值相关；见下文]。为了与阿尔伯特的方法进行比较等。(2007)，我们使用了网络同步（Kachalo等。,2008)在扰动图的未加权版本中应用他们的传递约简算法。我们还实施了Tresch提出的SOS（保存我们的标志）修剪程序等。(2007). 这种传递约简方法也适用于加权图，但与TRANSWESD在两个关键方面有所不同：（i）如第2节与NET-SYNTHESIS类似，这种剪枝方法接受非基本路径（包含负循环）来解释边，并且（ii）路径的长度作为边概率的乘积进行计算。基于此度量，边缘z（z）→v（v）如果存在（基本或非基本）路径，则删除z（z）⇒v（v）总概率大于边缘概率的同一符号z（z）→v（v）实际上，SOS剪枝区分了两种概率，一种是具有正边缘(第页₊)一个是因为有负边(第页₋)和标志秒边缘的符号为第页₊−第页₋对于DREAM设置，我们使用了与TRANSWESD相同的扰动图，并将条件相关系数分配为第页₊（如果秒为正）或第页₋（如果秒为负），而另一个概率设置为零。

表1总结了所有五个网络的结果，并指出了每个推断网络的TP/FP/TN/FN边数、计算时间和评估重建网络质量的标准统计指标。后者由DREAM评估脚本确定，包括AUROC、AUPR以及pAUROC和pAUPR（通过随机网络链路排列获得给定或更大AUROC/AUPR值的概率，如从10万次运行中估计的；参见Stolovitzky等。,2009). 由于AUPR值在稀疏网络中更为敏感，因此评估重建基因调控网络的质量特别有用。

这个P（P）-数值表明，我们的方法产生的结果明显优于随机选择的网络。这甚至仅适用于扰动图，表明将观察到的变化分类为相关或不相关的适当阈值策略可以提供大量有意义的信息。请注意，即使通过我们的阈值策略获得的原始扰动图也比提交给DREAM4挑战的许多解决方案要好。在五个网络中的四个网络中，将TRANSWESD应用于扰动图增加了AUPR值，而在一个网络中(2)AUPR略有下降。正如预期的那样，通过过渡还原法去除FP不可避免地伴随着一些TP的去除。这通常意味着AUROC值降低。然而，在大多数情况下，这种减少比AUPR的改善低一个数量级。当观察P（P）-值以及true和FP的数量。例如，在网络3中，FP的数量可以减少98个（从291个减少到193个），只牺牲3个TPs（从85个减少到82个），这导致AUPR从0.309适度增加到0.326，而P（P）-值从6.23e-111到1.97e-116。如果我们采用相同的数据，但基于无任何噪声的确定性模拟（这些数据是在发布挑战结果时提供的），则TRANSWESD的效果会更加显著。为了便于说明，表1显示为网络5获取无噪声数据时的结果。如果没有噪声，我们可以选择较小的阈值，例如ϑ=0.005和Γ=0。正如预期的那样，摄动图中的AUROC和AUPR已经高得多。但使用TRANSWESD，我们可以将AUPR（pAUPR）值从0.442（6.82e-102）进一步增加到0.567（9.11e-132）。FP的数量从476减少到110，而FN的数量从33适度增加到61。因此，高质量的扰动图可以提高传递约简的有效性。

在网络2中，我们观察到应用TRANSWESD时AUPR测量值减少，尽管FP的数量比额外的12个FN减少了71个。引人注目的是，在宣布DREAM4挑战的结果时，有人提到网络2出现了振荡，因此，提供的淘汰数据代表瞬态数据而非稳态数据。在这些条件下，扰动图中的一些边可能对应于非基本影响路径（具有负循环）的作用导致TRANSWESD较高错误率的扰动响应。

正如预期的那样，生成扰动图的计算成本一直很低。相反，应用完整TRANSWESD算法的运行时间因不同网络而异，并且可能会变得广泛。在两个网络（1和4）中，我们能够在合理的时间内将精确算法应用于摄动图。对于其他网络，我们在5小时后中断了精确算法，并使用了近似变量（近似最短路径计算，并且在删除边后不重新计算路径），在所有网络中，该变量需要<5秒。检查近似结果中的错误[删除边违反了条件(4)]，我们看到简化算法在五个网络中没有引入任何错误，这表明简化算法在实际基因调控网络中提供了合理的近似值。然而，当选择α>1时，错误数可能会迅速增加，尤其是当网络包含正循环时（未显示数据）。

我们还使用Albert的方法计算了传递归约等。(2007)在NET-SYNTHESIS中实现，输入与TRANSWESD相同的扰动图（无权重）。由于计算结果是非唯一的，并且对边缘排序非常敏感，我们对10种不同（随机）排序重复了计算，并获得了最佳结果[就AUPR值而言；为了公平比较，我们还根据摄动图中的边缘权重对最终图中的（剩余）边缘进行了排序]。该算法比完整的TRANSWESD快得多（但与近似版本相比：在所有网络中小于5 s），并且删除了更多的边缘，导致FN显著减少。然而，这是以消除相对较多的TPs为代价的。由此产生的AUPR值始终显著低于TRANSWESD获得的网络。结果甚至表明，通过网络同步获得的所有网络的AUROC和AUPR值（部分相当低）都低于生成它们的扰动图。这些结果表明，基于边缘权重的TRANSWESD衰减剪枝策略对于重构调控网络非常有利。也许正是出于这个原因，Tresch的SOS修剪策略等。(2007)-它也对加权图进行操作，产生了比NET-SYNTHESIS更好的结果。然而，AUPR值的改善（如果它增加了摄动图的AUPR）明显低于TRANSWESD，但非平稳网络2的问题情况除外，其中SOS剪枝表现稍好，但也导致AUPR比摄动图降低。通常，SOS剪枝删除的边缘明显少于TRANSWESD，这可能是由于SOS剪除用于量化路径长度的非常保守的度量（概率乘法）。

作为原理证明，结果如下表1，我们的算法在DREAM4中排名第3位（共19次提交）Insilico_Size_100号次级挑战。这一结果令人鼓舞，尤其是因为我们的方法只需要提供的模拟数据的一小部分。

4讨论和结论

在这项工作中，我们提出了TRANSWESD，它是传递约简的一个详细变体，适用于一类扩展的扰动图，即循环有向图、有符号有向图和加权有向图。主要变化和改进涉及：（i）生成加权扰动图和符号扰动图的新统计方法；（ii）使用边缘权重（关联强度）识别真正的冗余结构；（iii）周期的因果解释；（iv）放宽了传递归约的定义；（v）大型网络的近似算法。

传递约简的成功在很大程度上取决于扰动图的质量，因此取决于所选的阈值方法以及（间接）可用数据的类型和质量。虽然质量主要取决于信噪比，但数据类型（如基因表达、蛋白质水平、蛋白质磷酸化水平等）可能会对可观察到的扰动效应产生深远影响。

我们提出了一个生成扰动图的模块化过程，为FP约简方法提供了基础。此过程的主要任务是从噪声中过滤可解释的（直接和间接）影响。虽然这是FP减少方法（如TRANSWESD）的主要任务，但在此过程中也可以过滤间接影响。我们生成扰动图的工作流程由三个顺序模块组成：（i）规划和执行扰动实验；（ii）从实验数据生成有符号扰动图；以及（iii）从相关性度量中分配边缘权重（反映关联强度）。每个模块都可以进行交换或调整，例如，如果有其他类型的数据可用。例如，某些相互作用可能无法从单个扰动或/和稳态数据中推断出来，可能需要模块（i）中的特殊扰动策略。例如，如果通过布尔OR逻辑将此影响与其他影响结合，则只有多个淘汰将检测到一个节点对另一个节点的积极影响。在推导模块（ii）中的扰动图时，很容易集成单个和多个扰动的信息。此外，在生成扰动图时，还可以考虑瞬态响应阶段的数据以及模块（ii）和（iii）中的适当数据分析。请注意，根据具体的扰动数据（瞬态、稳态、时间进程），其他可能的非线性相关度量（如相互信息）可能更适合量化关联强度（Daub等。,2004)，尽管如果可以假设单调依赖性（唯一边缘符号），则线性度量似乎是合适的。也许生成扰动图的最关键步骤是将扰动效果分类为显著（引入边）或不显著。与相关性度量相比，该分类基于一个单一值。因此，实验重复将有助于为边缘提供更高的置信水平。一般来说，我们生成扰动图的方法所需的数据量相当大，因为所有节点都必须单独扰动，并且必须测量所有其他节点中的相应响应。然而，TRANSWESD的核心过程与推导扰动图的方法无关，因此，对于稀疏数据集，可以使用其他方法。

基准测试表明，我们用于生成扰动图的双阈值策略提供的网络本身已经具有相当高的重建质量。叶也做了类似的观察等。(2010)：从相关扰动效应中过滤噪声的简单噪声模型比更精细的微分方程模型具有更高的精度（尽管两者结合可以略微改善结果）。作者没有提到从间接影响中去除边缘，但他们提到这是一种改进结果的潜在方法。事实上，基准测试表明，我们的TRANSWESD算法可以通过仔细删除可能是FP的边缘来显著提高重建质量。

我们已经说明了有符号非循环图的传递约简在算法上相当简单。然而，即使在非循环图中，使用边权重也可能非常有益，因为它有助于避免消除真正的冗余结构。修剪循环图带来了几个问题，其中许多问题在TRANSWESD中通过使用边权重来解决。我们进一步说明了解释负循环结构所产生的因果关系的问题，并建议使用基本路径来解决更多计算时间的开销。通过这些特性，TRANSWESD在实际和客观基准测试中优于现有的传递约简算法。TRANSWESD的另一个优点是，如果所有边权重都是不同的，那么将提供一个按部分顺序（定义在边权重上）的唯一图形。将TRANSWESD与其他修剪策略进行比较是很有意思的，例如基于部分相关性（de la Fuente等。,2004)或数据处理不等式（Margolin等。,2006). 后两种方法最初是为无向图开发的，但可能会修改为有向图。

在具有多个循环的大型网络中，TRANSWESD可能会变得效率低下，因为它需要计算最短签名路径，这是一个NP完全问题。我们提出了两种计算近似解的方法，一种是通过近似计算最短有符号路径，另一种是在删除边时放弃重新计算最短路径。在基准中，近似解与精确解相同。

与阿尔伯特的方法类似等。(2007)，TRANSWESD可以通过为已知交互分配0的权重来轻松解释先验知识。如果路径用于解释边，则还可以为路径可能包含的边数引入上边界。到目前为止，TRANSWESD仅限于边缘去除。当集成多组扰动数据时，将TRANSWESD扩展到传递约简的一般情况可能是有益的，也就是说，允许以智能、数据驱动的方式插入边。这需要进一步调查。

总之，我们提出的重建工作流需要简单的数据，并以相对较高的概率提供边缘候选。边缘通过权重、符号和方向进行识别，为设计新实验和测试新假设提供了额外的关键信息。TRANSWESD是该工作流的一个重要组成部分，但它可以与其他逆向工程方法相结合，单独用作一般FP减少方法。

致谢

我们感谢Regina Samaga和Axel von Kamp对手稿的评论，感谢DREAM4组织者创建了生物信息学挑战并提供数据和黄金标准。

基金：德国联邦教育和研究部[HepatoSys，Virtual Liver；FORSYS-Centre MaCS（马格德堡系统生物学中心）]；萨克森大学教育和研究部（研究中心“过程工程和生物医学中的动力系统”）。

利益冲突：未声明。

参考文献

作者注释