植物抗旱性途径的贝叶斯模型

背景

植物是无梗生物，在极端环境条件下无法迁移到有利的位置。因此，他们别无选择，只能适应并最终适应恶劣的条件，以确保他们的生存。随着增强植物抵御压力环境的传统方法达到其生物极限，我们需要更新的方法来加强植物的内部防御机制。这些因素促使我们研究植物的遗传网络。WRKY转录因子因其在植物抵御生物胁迫中的作用而闻名，但最近的研究揭示了它们对干旱等非生物胁迫的作用。我们利用贝叶斯网络建立了WRKY转录因子网络模型，并应用推理算法寻找干旱响应的最佳调控因子。通过生物干预(激活/抑制)这些调控因子可以增强植物对干旱的防御反应。

结果

我们使用NCBI GEO数据库中的真实数据和贝叶斯网络中依赖关系生成的合成数据来学习网络参数。这些参数是用贝叶斯和频率论方法估计的。将这两组参数应用于基于效用的推理算法，以确定WRKY转录因子网络中植物干旱响应的最佳调控因子。

结论

我们的分析表明，在WRKY转录因子网络的所有其他元素中，激活WRKY18具有诱导干旱响应的最高可能性。我们的观察也得到了生物学文献的支持，因为WRKY18已知能积极调节干旱响应基因。我们还发现，激活WRKY60-60蛋白复合物具有诱导干旱防御反应的第二高可能性。与现有的生物学文献一致，我们还发现转录因子WRKY40和蛋白质复合物WRKY40-40可以抑制干旱反应。

到2050年，全球人口将增长34%，提高作物产量以确保粮食安全已成为一项重大挑战[1］．全球变暖导致的全球气温上升增加了干旱影响作物产量的风险，使这一挑战进一步复杂化。研究表明，到2050年，全球受干旱影响的地区将显著增加，并伴随着作物产量的急剧下降[2］．世界人口的空前增长伴随着对作物供应需求的增加，而与此同时，使作物产量最大化的传统方法正走向其生物极限。因此，开发抗旱作物已成为确保粮食安全的全球优先事项。幸运的是，植物有多种天生的应激感应机制，可以检测到环境中的不利变化，并部署适当的防御反应。因此，了解植物防御机制背后的遗传网络，以增加其遗传产量潜力，同时降低其对恶劣条件的敏感性，是非常有意义的。

脱落酸(ABA)是一种众所周知的植物激素，在干旱胁迫条件下被诱导，并通过转录因子的作用调节植物的基因表达[3.，4］．WRKY转录因子家族传统上与植物对抗病原体的防御机制有关;然而，最近许多研究强调了WRKY在非生物应激反应中的作用[5- - - - - -7］．由于WRKY是植物中最大的转录因子家族之一，在植物防御机制中发挥着如此多样化的作用，因此对WRKY信号通路的各个组分之间的相互作用进行建模是可行的，以获得对这些相互作用的有价值的见解[8］．在本文中，我们使用贝叶斯网络(BN)来模拟aba诱导的WRKY转录因子网络。然后，我们应用基于效用的推理技术来确定BN中干旱胁迫响应基因的重要调控因子。这种方法允许我们将现有的生物学知识整合到我们的模型中。

生物学背景回顾

类似于肾上腺素在动物体内作为应激激素的作用方式，植物通过分泌植物激素(如ABA、细胞分裂素、水杨酸和乙烯)来应对恶劣的环境变化、病原体的攻击或伤害，从而触发自身的防御机制。在植物的环境中，干旱的特点是无法获得水，这可以阻止植物进行基本的生存过程，如光合作用。当植物面临缺水条件时，它可以通过避免或容忍的过程来保护自己。在避免的情况下，植物可以在雨季完成它的生命周期。而在耐受性的情况下，植物可以通过改变其基因表达将可逆的变化引入其生理，从而初步适应条件的变化;然而，如果干旱条件仍然存在，那么植物就会把改变过的基因传递给下一代，这样这些新一代的植物就已经适应了干旱条件。9］．

为了采用这两种防御机制中的任何一种，当植物得到干旱的初始信号时，例如外质体中水势的下降和离子浓度的上升，它必须经历一个信号转导过程[9］．所有这些信号和许多其他信号一起导致植物激素ABA水平的迅速上升，ABA起着压力传感器的作用，随后它激活了次级信使，例如C一个^{2 +}活性氧(ROS)和环磷酸腺苷(cAMP)。这些二级信使启动各自的信号通路(例如MAPK, CDPK)，其中蛋白质磷酸化(添加磷酸盐)\ ({PO_ {4} ^ {3}} \))，而去磷酸化可分别通过激酶(酶)和磷酸酶(酶)的作用进行[10］．在激酶和磷酸酶的信号传导作用下，转录因子被激活或失活以调节下游基因的表达[10］．转录因子是与基因的特定DNA序列结合以激活或使其失活的蛋白质。最后，转录因子直接负责启动应激反应基因并关闭任何其他非必需基因。

WRKY、bZIP、NAC等转录因子家族都调控大量基因。因此，了解转录因子的活性对于了解植物的胁迫响应机制至关重要。WRKY是一个转录因子家族，在植物抵御非生物和生物胁迫的防御机制中发挥作用。直到最近，WRKY在处理非生物胁迫中的作用还没有像在生物胁迫中那样被广泛地探索。由于这些原因，缺乏可用的实验数据[3.］．在这篇论文中，我们感兴趣的是研究WRKY转录因子信号通路各成员之间的相互作用(图。1)在干旱胁迫下ABA快速诱导。了解这些相互作用可以让我们更深刻地了解这一途径的功能，这可以帮助我们制定培育抗旱植物的干预策略。

研究表明，ABA在水分亏缺和盐胁迫条件下诱导转录因子WRKY18、WRKY40和WRKY60 (Chen et al.) [11］．此外，也有报道WRKY18和WRKY60在抑制种子萌发、根系生长和增强植物对水分亏缺胁迫敏感性方面对ABA具有正敏感性;相反，WRK40可以拮抗WRKY18和WRKY60，从而影响植物的ABA敏感性和非生物胁迫反应(Chen et al.)。用WRKY18和WRKY40缺失突变体进行实验，结果表明WRKY60的表达可以忽略不计。这意味着WRKY18和WRKY40通过识别WRKY60启动子中的一簇W-BOX序列直接诱导WRKY60 (Chen et al.)。除了这三种WRKY转录因子的各种调控行为外，人们注意到这三种转录因子不仅相互作用形成三种同质复合物，而且相互作用形成异质复合物[12］．

贝叶斯网络建模

生物网络本质上是曲折和随机的。通常很难解释网络不同组成部分之间的多元相互作用。BN是一个有向无环图，它决定了网络中一组随机变量联合概率分布的条件分解，从而简化了其联合概率分布的计算(Sinoquet and Mourad) [13］．因此，我们对使用神经网络来模拟生物网络中的相互作用很感兴趣，因为它们提供了一个清晰而紧凑的框架来表示联合概率分布，并从这些网络中得出推论[14］．检查bn可以帮助增强我们对网络中不同元素之间关系的信念，并提供对网络因果关系的洞察。

在本文中，我们建立了WRKY信号通路的模型。1)参与模式植物拟南芥的干旱胁迫响应。基于图中所示的信号转导途径。1，我们构造了一个BN，如图所示。2．每个圆形节点(A,B,C，…，H)代表一个基因、转录因子或蛋白质复合体，节点之间的每条有向边代表存在于WRKY信号通路中的因果关系。每个节点上都有一个矩形，表示与该节点相关的参数或局部概率模型。例如,θ_C|一个，B表示节点C在父节点A和b的情况下的条件概率密度。这些参数可以从数据中学习到，对于理解整体图结构很重要。

参数估计

根据给定应用程序中先验知识的可用性，可以使用频率论方法或贝叶斯方法来学习BN的参数。频率论方法，如最大似然估计(MLE)假设所学习的参数是固定的，并在不考虑先验信息的情况下产生点估计。另一方面，贝叶斯估计将参数视为随机变量，利用参数的数据和先验分布来获得参数的后验分布。此外，贝叶斯方法考虑了零概率估计的问题，这可能会影响学习算法。即使先验信息遵循均匀分布(非信息先验)，贝叶斯方法也提供了非零概率估计。这是因为后验信念是由数据和先验知识共同支配的，因此对概率的零估计只与事件的不发生有关。然而，贝叶斯估计过程在计算上具有挑战性，因为它需要执行积分以获得证据(数据)的概率。由于这个原因，我们在贝叶斯估计过程中对给定似然函数使用共轭先验的概念[15］．

在图的BN中。2，我们假设BN中的每个节点X只能获得二进制值，X = 0或X = 1。当一个节点X =1时，表示该节点所代表的基因、转录因子或蛋白质复合体被激活，当X =0时，则相反(基因、转录因子或蛋白质复合体被抑制)。这个公式允许我们使用伯努利分布对网络中每个节点的状态进行建模，给定其父节点的状态。现在，考虑一个有N个节点的BNθ_X是X =1(成功)和1-的概率θ_X为X =0(失败)的概率。假设我们对每个节点的状态进行n(>0)次观察，并设k(≤n)为节点状态为1的次数。我们进一步假设对每个节点进行n次观测后得到的随机变量X1, X2，…，Xn序列是独立且同分布的。因此，给定父节点(P_一个（X)符合二项分布，由:

为了估计后验分布，我们需要定义参数的先验θ_X对于我们的模型。由于与我们的模型相关的似然函数是二项的，我们选择先验分布来遵循带有一些形状参数的Beta分布(α_X，β_X)，结果为:

由于将先验建模为二项似然下的beta分布，从共轭族的性质可以得出，后验也将遵循具有形状参数的beta分布\(\左(\alpha ^{'}_{X}，\beta ^{'}_{X}\右)\)［15］．在我们的模型中，参数的后验分布θ_X由:

在哪里\(\alpha ^{'}_{X} =\)（α_X+ k)和，\(\beta ^{'}_{X} =\)（β_X+ n - k)，该分布的期望值为:

我们可以利用实验数据进行迭代更新α_X而且β_X得到后验分布。随着数据的增多，后验分布将向实际的后验分布收敛。我们将先验建模为二项式似然下的Beta，这使我们能够获得后验的封闭形式解。其他的非共轭先验也可以使用;但是，封闭形式的解决方案可能无法保证。注意，这种方法给出了与每个节点相关的边际和条件后验分布，而不是它们的概率(θ_X，θ_Y|X)．在本文中，为了学习这些概率，我们用期望值(Eq。5)对应节点的后验分布。此外，我们还使用MLE的频率方法学习概率，以便比较使用两种方法得到的最终结果。理想情况下，当数据丰富时，贝叶斯方法和MLE估计收敛于同一点[16］．二元随机变量的边际概率和条件概率可以用等式估计。6而且7分别。

其中M [X¹]是随机变量X = 1的次数，M[X⁰为X为0的次数，M [Y¹，X⁰是X = 0 Y = 1的次数，M[Y⁰，X⁰]是X等于0 Y等于0的次数。我们在BN建模中所做的一个关键假设是，节点集的联合分布根据图中的BN进行因式分解。2．这一假设意味着，生物结构中的依赖性反映在我们从中学习网络参数的数据中。人们可以使用基于约束或基于分数的学习技术从数据中推导出图结构，然后随后学习网络参数[17］．然而，在本文中，我们避免学习图结构，因为公开的WRKY转录因子在非生物胁迫下的实验数据非常有限，而且我们的网络中包含蛋白质复合物(节点C、D、E和G)，其表达数据不与基因表达数据(节点A、B、F和H)同时存在。通常，包含基因表达数据的数据集不包含蛋白质复合物的表达数据，反之亦然。因此，利用BN中的依赖关系和网络中其他非蛋白质复合物节点的实验数据生成了蛋白质复合物的合成数据。

基于效用的贝叶斯网络推理

从数据中推断出网络参数或与每个节点相关的局部概率后，BN就有了足够的信息进行推理。我们的目标是在WRKY BN中找到一个能最大限度上调抗旱基因下游表达的节点。换句话说，我们感兴趣的是在BN中找到一个单一的节点(节点a - g)，当上调或下调时，我们的应激反应基因(节点H)上调的机会最大化。

在BN中执行推理有多种方法。每当我们有一个单连通图时，Pearl的消息传递算法就会受到青睐，因为它允许我们执行精确的推断[18］．然而，这里考虑的BN不是单连通的，也有循环，这是Pearl算法无法处理的。其他基于非精确采样的技术需要大量的数据来提供可靠的推断。因此，在本文中，我们考虑了另一种类型的近似推理技术，该技术根据在特定节点上采取的操作计算一个分数，通常称为预期效用。效用衡量的是行动的效果。为了在BN中实现实用程序，我们首先需要理解贝叶斯决策网络的概念，以及如何从BN中创建一个。为了说明这些概念，考虑如下图所示的一个简单BN的例子。3.．

图中的BN。3.有三个节点基因A，基因B和基因C，我们假设每个节点的二进制值为0(抑制)或1(激活)。基因A和基因B是基因C的父节点，具有与之相关的边际概率，如图所示。3.．同时，我们假设当基因A活跃时，它会激活基因C，当基因B活跃时，它会抑制基因C。基于这个BN，我们构建了一个贝叶斯决策网络，如图所示。4．矩形节点作为决策(行动)节点，菱形节点作为效用节点，圆形节点表示机会(自然或概率)节点。在这个例子中，我们感兴趣的是让基因C的值为1，这是预期效用所衡量的。在这种情况下，我们可以选择在机会节点基因A或基因b上采取行动。一旦我们决定在机会节点上采取行动，它就不再是一个机会节点，而是成为一个确定性节点。

根据所采取的行动，预期效用可由式计算。8：

在P (O_我|一个)表示结果的概率(O_我)，与动作A一致，U(O_我)表示动作a下该结果的效用值。效用表在table中定义1，其中第一行表示基因A活跃而基因B被抑制的最佳情况，最后一行表示基因A被抑制而基因B活跃的最坏情况。其余的行表示其他可能的场景。分配的效用分数相对于最佳(最高效用)和最差(最低效用)情况场景，这些值不是唯一的，可以以不同的方式重新定义，但是，分数必须反映决策网络中描述的场景。

用情商。8我们首先计算对基因A采取行动的期望效用如下:

案例1:动作:基因A被激活(A = 1)。

案例2:所采取的行动:基因A被抑制(A = 0)。

因此，当基因A =1或激活时，预期效用更大。同样，让我们计算一下对基因B采取行动的预期效用。

案例1:采取的行动:基因B被激活(B = 1)。

案例2:所采取的行动:基因B被抑制(B = 0)。

因此，当基因B被抑制时，预期效用更大。然而，基因A被激活的效用大于基因B被抑制的效用。因此，我们必须选择激活基因A而不是抑制基因B，以最大限度地提高C基因被激活的几率。

为了估算这些参数并在BN中进行效用计算，我们需要获得干旱胁迫条件下WRKY转录因子的数据。由于WRKY转录因子最近才被发现在非生物应激反应中起作用，因此很难获得公开的大规模数据。然而，我们能够从数据集GSE46365、GSE65046和GSE76827中获得BN中所有基因和转录因子(节点A、B、F和H)的真实微阵列基因表达数据，这些数据集可从NCBI GEO数据库公开获得[19- - - - - -21］．这些数据集被单独归一化、二值化并聚合成一个复合数据集，每个非蛋白质复合体节点(基因和转录因子)包含116个数据点。一旦真实数据的关键,他们一起使用的依赖BN为蛋白复合物生成数据用节点C, D, E, g .例如,为了为节点生成数据集D的表达值节点和节点H观察,即所有的父母和孩子节点D .如果两个节点A股和H观察调节(状态= 1)节点D分配确定性是调节(状态= 1),如果节点A和H都被观察到下调(状态=0)，则节点D被确定性地分配为下调(状态=0)。这是因为我们从生物学文献(Chen等)中了解到，节点A上调节点D，节点D又上调节点H。如果节点A的表达状态上调(state =1)，节点H的表达状态下调(state =0)，则节点D被赋值为1(上调)，概率大于0.5。这是因为如果节点A被上调，节点D也很可能被上调，但不足以抵消节点E对节点H的下调效应，这可能导致节点H被下调。节点D被上调的概率是从[0.6,0.7,0.8,0.9和1]的一组离散概率值中随机选择的，其中每个值被选中的概率相等。类似地，当节点A =0且节点H =1时，节点D被概率地赋值为0(下调)。以这种方式也生成了节点C、G和E的数据，因此这些合成数据反映了网络依赖关系和非蛋白质复合体节点的真实数据。表2，3.，4，5，下面给出的伪代码1以及附加在附加文件部分的R代码(附加文件1，2，3.，4，5，6，7，8，9而且10)进一步解释了节点D, E, G和c的数据生成。一般来说，基因表达数据集不包含蛋白质-蛋白质相互作用数据，这是我们网络中避免生成蛋白质复合物合成数据所需要的。虽然我们可以找到蛋白质-蛋白质相互作用数据集，但这些数据集不包含基因表达数据，因此为了避免这个问题，我们考虑为BN中的蛋白质复合物生成合成数据。

一旦生成了所有蛋白质复合体节点的合成数据，它们就会与真实世界的数据一起聚合在一个数据集中。该数据集用于估计网络参数的目的，使用贝叶斯方法和参数估计部分概述的最大似然方法。对于贝叶斯方法，每个节点的先验首先被初始化为beta(1,1)分布，这是一个在区间[0,1]上的均匀分布。使用数据和Eq。4更新每个节点的后验分布，并使用Eq计算期望值。5．期望值近似为节点的条件概率。使用MLE方法以及使用等式学习了一组单独的参数。6而且7．

然后，计算单点干预的效用。将效用节点设置在节点H，进行效用分析，找出单个干预点上调节点H6)基于贝叶斯决策网络进行定义。节点H的效用直接依赖于节点D、G和E。最好的情况是节点D和G上调，节点E下调;最坏的情况是节点D和G下调，节点E上调。这些场景代表了网络中实际的生物过程，效用分数是相对于这些场景定义的，效用分数高是有利的。模拟使用R软件[22]，效用计算使用Netica [23］．

图中的分割条形图。5表示预处理数据后BN中每个基因的激活和抑制状态。表格7显示使用贝叶斯和MLE方法估计的条件概率。使用贝叶斯和MLE方法得到的参数得到的最大期望效用如图所示。6而且7分别。

从使用贝叶斯和MLE方法的效用分析中，我们发现WRKY18的激活效用得分最高。这意味着WRKY18的上调是导致下游干旱胁迫反应基因上调最有效的单基因干预。这一结果与生物学文献一致，表明WRKY18在干旱胁迫条件下对ABA具有正向敏感性，并在下游基因表达上调中起关键作用。我们也可以从图中的条形图中看到。6而且7第二个最佳的干预点是蛋白质复合体WRKY60-60，它也上调了下游干旱应激反应基因的表达。此外，与文献一致，我们看到WRKY40和WRKY40-40都具有很高的抑制效用，因为它们负责下调下游干旱胁迫反应基因。我们还看到，我们在贝叶斯和MLE方法中的效用分数是可比的，这是由于估计的概率(表7)使用这两种方法是非常相似的。

在本文中，我们介绍了WRKY信号通路，传统上它与植物对生物胁迫的防御反应有关，但最近它被证明在植物对非生物胁迫(如干旱)的防御反应中发挥着重要作用。由于其在植物防御中的多种作用，这是一个有趣的途径选择。我们使用BN对WRKY通路进行建模，其中网络中的每个节点都代表通路中的一个基因、转录因子或蛋白质复合体，节点之间的每条边都代表通路中存在的因果关系。与每个节点相关联的是一个条件或边际概率，它表示节点被激活或被抑制的概率。对于我们的分析，我们假设网络中的节点只能取1(激活)或0(抑制)的二进制值。由于BN可以捕捉因果生物学关系和生物途径的概率性质，因此它是建模目的的理想选择。

为了学习网络中的参数，我们使用了真实的基因表达数据，并生成了反映网络依赖关系的合成数据。为了估计局部条件概率和边际概率，使用了贝叶斯方法和频率论方法。在贝叶斯方法中，我们假设每个节点都有Beta(1,1)的先验分布(在范围[0,1]内均匀分布)，这意味着我们对我们的模型没有先验知识。由于我们的似然服从二项分布，我们能够通过共轭族的性质使用封闭形式公式来得到每个节点的后验分布。后验分布的期望值被用作局部概率的估计。我们选择共轭族是为了简化我们的计算，并得到一个封闭形式的解，然而，选择共轭先验并不总是最好的选择。如果有足够的信息可用，先验可以使用非共轭族分布建模，后先验可以使用(Markov-Chain-Monte-Carlo) MCMC技术估计，尽管这可能需要大量的计算。在频率论方法中，我们简单地使用最大似然估计来获得局部概率。使用这两种方法得到的概率被发现非常相似。

一旦学习了每种方法的参数，就可以使用效用的概念来进行推断干预的最佳节点的任务。在我们的模型中，我们使用了一种非精确推断技术，因为我们不能使用精确的技术，例如Pearl的消息传递算法在我们的贝叶斯网络中，因为前者只适用于单连接和环路较少的网络。此外，数据点的数量非常有限，这使得用于推理目的的实用程序的选择相对于基于数据密集型采样的推理技术非常合理。此外，基于效用的推理可以很容易地应用于更大的BN。使用贝叶斯和MLE方法的参数进行的效用分析显示，上调时WRKY18是最有效的干预节点。因此，上调WRKY18会进一步上调WRKY信号通路下游的应激反应基因。这一结果与生物学文献一致，即WRKY18在干旱条件下积极上调干旱反应基因的基因表达。我们在这项研究中的下一步将是在贝叶斯参数估计方法中探索和实现更多的信息先验，而不是使用我们在这里使用的非信息Beta(1,1)先验，因为我们的研究受到缺乏关于网络的先验知识的限制。我们还想研究与植物对干旱的防御反应有关的其他信号通路，并找到这些网络中的关键调控因子，并将其功效与WRKY18进行比较。我们的研究范围仅限于单节点干预，未来我们还希望将研究扩展到多节点干预，并研究这些干预如何调节多个下游基因。

我们利用贝叶斯网络建立了植物干旱响应单通路模型，并应用了基于效用的推理算法，结果显示WRKY18激活后最有可能激活下游抗旱基因。这一结果与生物学文献以及基于效用分析的其余结果一致。未来，我们计划利用CRISPR-CAS9在田间植物中激活WRKY18，以检测WRKY18抗旱效果。这种使用贝叶斯网络寻找干旱响应调控因子的过程，可以应用于寻找其他植物网络中的关键调控因子，这有助于在未来培育出健壮而有价值的作物。

阿坝:

脱落酸

BN:

贝叶斯网络

ROS:

活性氧

阵营:

环磷酸腺苷

大中型企业:

极大似然估计

并:

最大期望效用

密度:

马尔科夫链蒙特卡洛

我们感谢BMC植物生物学的审稿人和编辑对改进这篇手稿的建设性反馈。我们也要感谢德州农工大学写作中心在编辑这篇手稿方面的帮助。

资金

这项工作部分得到了TEES-AgriLife生物信息学和基因组系统工程中心(CBGSE)启动资金的支持，德克萨斯州A&M X-Grant计划，部分由国家科学基金会资助ECCS-1404314。资助机构在研究设计、数据收集、分析和解释以及撰写手稿方面没有发挥任何作用。

数据和材料的可用性

本研究中使用的数据集在NCBI公开提供，登录号为GSE46365、GSE65046和GSE76827。从这三个数据集中提取的支持本文结论的数据子集包含在本文中。

从属关系

1. 德州农工大学电气与计算机工程系，188 Bizzell St, College Station, 77843，美国

   Aditya Lahiri, Priyadharshini S. Venkatasubramani和Aniruddha Datta

贡献

AL发现了这个问题。PSV和AD对建模和推理的方法论做出了贡献。AL实施了该方法并起草了手稿。所有三位作者都阅读并批准了手稿。

相应的作者

对应到Aditya Lahiri．

伦理批准并同意参与

不适用。

发表同意书

不适用。

相互竞争的利益

作者宣称他们之间没有利益冲突。

出版商的注意

施普林格自然对出版的地图和机构从属关系中的管辖权主张保持中立。

作者的信息

不适用。

开放获取本文根据创作共用属性4.0国际许可协议(http://creativecommons.org/licenses/by/4.0/)，允许在任何媒介上不受限制地使用、分发和复制，前提是您对原作者和来源给予适当的赞扬，提供到创作共用许可证的链接，并注明是否进行了更改。创作共用公共领域奉献弃权书(http://creativecommons.org/publicdomain/zero/1.0/)除另有说明外，适用于本条所提供的资料。

转载及权限