预测影响特定代谢产物产量关键反应的通量域比较方法

本申请得到国家自然科学基金（NO.21106095）和天津师范大学引进人才基金课题的资助。 

技术领域

本发明属于生物信息学技术领域，涉及预测影响特定代谢产物产量关键反应的通量域比较方法。 

背景技术

目前，国际国内的研究中，已经构建的基因组尺度代谢网络模型超过100个，现有的计算机模拟方法主要针对野生型物种通过最优化的方法（如通量平衡分析）进行模拟和分析。分析的方法包括： 

（1）特定产物产量最大化时物种内部各反应的通量分布，即以特定产物的合成反应为目标函数，在给定的培养基条件下，计算代谢网络中各个反应的通量值；

（2）通过单基因敲除预测必需基因，即在其他参数设定不变的情况下，将单个基因对应的反应通量设定为0，计算目标函数的值，目标函数值为0则该单个基因为必需基因，目标函数不为0则该单个基因为非必需基因；

（3）    鲁棒性分析，即分析单个反应的流量变化对目标函数最优值的影响。

（4）通量可变性分析,通量可变性分析在给定的稳态空间中，对每个反应都做两次优化，分别求出最大值和最小值。用最大和最小值之间的区间大小来度量反应的可变性。 

然而，在很多应用最优化算法的方法中，由于最优化算法本身的结果不具有唯一性，因而根据一次计算的结果对反应通量进行分析具有很大的特殊性和片面性。其次，现有方法仅局限于针对野生型物种的分析，没有将野生型物种与优化菌种的模拟相结合，因而难以产生对物种改造策略有实际指导意义的预测结果。 

本发明人于2012年申请了（申请号：201210099321.2） 基于基因组尺度代谢网络模型的代谢工程设计预测方法专利，随着实验的不断进行发现，此专利申请仅考虑了单次通量平衡分析的结果，而没有对通量域内的其他最优解进行分析，并且没有对预测结果进行排序，因而本专利利用通量域的比较进一步完善了上一个专利的内容，并且对预测结果进行了排序，从而对生物实验提供更具有指导意义的模拟结果。 

发明内容

本发明提供的方法引入了通量域比较的方式，并将野生型物种的生长状态与优化型工程物种相结合，能够针对特定代谢产物，得到与提高其产量密切相关的关键反应，从而对生物实验及工业生产的实验设计产生实际的指导价值。该方法将代谢网络模型用于对物种改进的计算和预测，对于生物实验和工业生产具有重要的指导意义。 

为实现上述目的，本发明公开了如下的技术内容： 

一种预测影响特定代谢产物产量关键反应的通量域比较方法，其特征在于该方法由两次计算，一次比较、一次排序组成，按如下的步骤进行：

（1）计算野生型目标物种自然生长状态的通量域：

首先采用通量平衡分析的方法，利用已知的基因组尺度代谢网络模型，计算生物量合成速率的最大值，计算时参数设置方法如下：根据已知的目标物种的培养条件及表型，将模型中培养基组成、溶氧、目标产物产量、其他副产物产量按照实验条件进行相同的设置；目标函数设定为生物量合成反应，计算得到生物量合成速率最优值f；

之后，采用通量可变性分析的方法，得到野生型物种的代谢网络通量域分布，计算时参数设置方法如下：将生物量合成反应的速率设定为上一步骤中得到的最大值f，其它条件与上一步骤相同，计算得到两个向量V_{w_min}和V_{w_max}，其中V_{w_min}表示代谢网络中各反应通量值的最小值（即通量下限）组成的向量，V_{w_max}表示代谢网络中各反应通量值的最大值（即通量上限）组成的向量；计算所得的通量下限与通量上限组成的区域即为野生型物种的代谢通量域分布；

（2）模拟优化型物种的生长状态：

第一，采用通量平衡分析的方法，利用已知的基因组尺度代谢网络模型，计算生物量合成速率的最大值，计算时参数设置方法如下：根据已知的目标物种的培养条件及表型，将模型中培养基组成、溶氧按照实验条件进行相同的设置，而对目标产物产量和其他副产物的产量不作限制，生物量合成反应的通量设置为步骤（1）中生物量合成反应的最优值f，目标函数设定为特定代谢产物的产量，进行第二次模拟，计算得到特定产物产量的最大值p；

第二，采用通量可变性分析的方法，得到优化型物种的代谢网络通量域分布，计算时参数设置方法如下：将特定代谢产物的产量设定为上一步骤中得到的最大值p，其它条件与上一步骤相同，计算得到两个向量V_{e_min}和V_{e_max}，其中V_{e_min}表示代谢网络中各反应通量值的最小值（即通量下限）组成的向量，V_{e_max}表示代谢网络中各反应通量值的最大值（即通量上限）组成的向量；计算所得的通量下限与通量上限组成的区域即为优化型物种的代谢通量域分布；

（3）一次比较：将两次计算的通量域进行比较，确定优化型与野生型之间反应通量域的不同，从而确定将实验型改造成优化型所需进行的代谢工程改造，进而制定相应的湿实验策略，比较及预测方法如下：

对于代谢网络中的每一个反应，

其中，v_{w_max}表示野生型物种模拟得到的通量上限，v_{w_min}表示野生型物种模拟得到的通量下限，v_{e_max}表示优化型物种模拟得到的通量上限，v_{e_min}表示优化型物种模拟得到的通量下限；

上述预测结果有一个特殊情况，当(v_{e_max}-v_{w_max})/ v_{w_max}与(v_{e_max}-v_{w_max})/ v_{w_max}计算得到的结果符号相反时，预测结果为未知，即无法得到准确的预测结果；

（4）一次排序：对于每一个预测得到的反应，将其代谢通量分别设置为野生型通量的上限或下限，这样有两种设置方式，以v_l表示反应通量的下限，v_u表示反应通量的上限，则这两种设置方式可以表示为：

v_l=v_u=v_{w_min}   和  v_l=v_u= v_{w_max}

其他参数的设置与步骤（2）相同，每次设定后都以特定产物的产量为目标函数进行计算，计算方法为通量平衡分析，两种条件的设定各得到一个计算结果，取两个计算结果的最大值，记为q，计算该反应通量改变后特定产物产量的改变量C，即

将所有预测反应的C值按照从大到小的顺序进行排列，即为预测反应对特定产物产量影响大小的排序，根据该排序对代谢工程基因操作的实验策略进行指导。

利用将野生型物种与优化型物种计算得到的通量域进行比较，对提高特定代谢产物合成量的关键反应及其对产物产量的影响进行计算机预测，从而有针对性的指导生物实验的进行。需要说明的是：首先保证作为基础的全基因组尺度代谢网络模型具有较高的质量。其次实施者需具备利用COBRA工具包进行编程计算的能力，同时具有对大量模拟结果通过计算机编程进行比较的能力，具备重构基因组尺度代谢网络模型。 

本发明更进一步公开了预测影响特定代谢产物产量关键反应的通量域比较方法在制备枯草芽孢杆菌在最小培养基上进行核黄素生产中的应用。 

申请号：201210099321.2也用到了此方法，两个方法应用于同一个模型，本例中采用了通量域比较法，因而预测的结果更加精确，从数量上看也可以看出来，上个专利预测了127个，本次预测了112个反应，排除了一些价值不大的预测结果，另外，本方法对预测的反应进行了影响力的排序，因而具有更好地指导效果，详细方法如下： 

将该方法用在枯草芽孢杆菌最小培养基上进行的核黄素生产上。所采用的模型是已发表的野生型枯草芽孢杆菌基因组尺度代谢网络模型，实验菌生长条件的数据来自史硕博等2009年在ME上发表论文^[1]中Bacillus subtilis168的实验数据，该菌种在最小培养基上进行培养，最小培养基包含的元素组成及其在模型中通量的上下限如表1所示。各通量交换反应上下的设定主要依据以下几点：

(1) K+，Na+，Mg²⁺，Ca²⁺，Fe³⁺，CO₂，H₂O和H+可以自由地进入和离开网络，因此它们通量反应的上限为1000 mmol g DW-1 h-1，下限为-1000 mmol g DW-1 h-1；

(2) 由于枯草芽孢杆菌为需氧菌，因此设氧气可以自由进入网络，却不可能离开网络（即生成氧气），因此设氧气通量交换反应的下限为-1000 mmol g DW-1 h-1，上限为0；

(3) 碳源、氮源、硫源、磷源作为从培养基中获得的有限底物，其最大吸收速率根据实验设定为5 mmol g DW-1 h-1，即其通量交换反应的下限为-5 mmol g DW-1 h-1，上限为1000 mmol g DW-1 h-1。其中，根据实验条件，菌体利用的碳源为α-D-葡萄糖。

表1 最小培养及条件及通量设置 

根据上文介绍的方法，步骤如下：

（1）按照表1给出的培养基组成和通量上下限对模型中相应组分交换反应的通量上下限进行设置，代谢产物（包括乙酸和核黄素）交换反应上下限的设置按照史硕博等2009年在ME上发表论文^[1]中的数据进行设置，即乙酸交换反应上下限均为1.63 mmol g DW^-1 h^-1，核黄素交换反应上下限均为0。以生物量合成反应为目标函数，进行第一次模拟计算，得到生物量合成量的最大值f=0.334gh^-1, 再将生物量合成反应的上下限同时设定为0.334 gh^-1，进行FVA分析，得到野生菌的通量域分布。

（2）在第二次模拟计算中，生物量合成反应的上下限同时设定为0.334 gh^-1，乙酸交换反应上下限分别设定为默认值，即下限为-1000 mmol g DW^-1 h^-1，上限为1000 mmol g DW^-1 h^-1，其他限制条件与第一次计算相同，以核黄素交换反应为目标函数进行计算，得到核黄素的最大合成量为0.109 mmol g DW^-1 h^-1。再将核黄素的交换反应通量上下限同时设定为0.109 mmol g DW^-1 h^-1，进行FVA分析，得到优化菌的通量域分布。 

（3）将两次计算得到的反应通量域进行比较，当一个反应在优化型中的通量域范围高于野生型1%时，表明该反应需要通过过表达（overexpress）进行优化；当该反应在优化型中的通量域范围低于野生型1%时，表明该反应需要通过低表达（underexpress）进行优化；当该反应在优化型中通量域的上下限均为0，而在野生型中通量域上限或下限有一个不为0时，表明该反应需要通过敲除（knockout）进行优化。结果表明共有112个反应需进行改造，其中有80个需高表达的反应，27个需低表达的反应，另有5个反应由于其通量域的上限和下限改变的方向相反，因而不可预测。这些反应在不同代谢系统中的分布如图1所示。需高表达的反应多分布在碳水化合物和核酸代谢中，这与核黄素的两个前体（Ribulose-5-phosphate 和 GTP）分别分布在碳水化合物和核酸代谢中有关。其中prs^[1],purFMNHD^[2] 和 ribABGHT^[3,4]基因对应反应的高表达预测在之前的研究中得到了证实。 

（4）对上述112个预测反应，在优化型模型中，将其代谢通量分别设置为野生型通量域的上限（v_{w_max}）或下限（v_{w_min}），每次设定后都以核黄素的产量为目标函数进行计算。对于每个反应，取其min和max条件下计算得到结果的最大值，计算该反应通量改变后核黄素产量的改变量C，将所有预测反应的C值按照从大到小的顺序进行排列，即为预测反应对核黄素影响的排序。在反应对核黄素产量影响的排序中，由PRPP合成GTP的途径中包含的其他pur基因（purQLCEK）和ndk基因被预测为过表达，且对核黄素产量的影响最大，另外，一些合成谷氨酸盐或甘氨酸反应的过表达对核黄素的影响非常大，而消耗谷氨酸盐或甘氨酸反应的低表达同样对核黄素的产量非常大，这一点表明了谷氨酸盐和甘氨酸在核黄素合成过程中有非常重要的作用。 

本发明公开的预测影响特定代谢产物产量关键反应通量域比较方法所具有的优点在于： 

（1）通过通量域的比较，完整的覆盖了所有可能的通量分布情况。

（2）结合野生型和优化型物种的计算结果进行预测，对反应的过表达和低表达情况进行了预测，体现了优化型物种代谢通量的变化情况。 

（3）通过计算对预测得到的反应进行影响力大小的排序，从而对湿实验及生产产生具有实际应用价值的指导。本发明可应用于任何具有基因组尺度代谢网络的物种，以及模拟预测网络模型计算能力范围内的任何产品，特别对于基因序列尚不明确的代谢工程菌具有很好的指导意义。 

为了验证该方法在代谢工程菌上的应用性，按照文献“Transcriptome analysis guided metabolic engineering of Bacillus subtilis for riboflavin production”和文献“Increased production of riboflavin by metabolic engineering of the purine pathway in Bacillus subtilis”两篇文章的实验条件分别模拟预测了RH33和B.subtilis PK两株工程菌的核黄素生产情况，结果表明在实验给定的最小培养基条件下，对RH33的模拟得到了118个预测反应，其中77个为过表达，40个为低表达，1个为未知。对B.subtilis PK的模拟得到了102个预测反应，其中65个为过表达，35个为低表达，2个为未知。该方法对RH33和B.subtilis PK的模拟结果进一步印证了生物湿实验中prs和purFMNHD，以及ribABGHT的过表达对提高核黄素产量的作用。同时，预测结果表明由PRPP合成GTP的途径中包含的其他pur基因（purQLCEK）和ndk基因被预测为过表达，其对应反应对核黄素产量的影响最大。除此之外，大量未经实验验证的预测反应为进一步的实验提供了由针对性的指导。RH33和B.subtilis PK都是经过多次进化实验和得到的菌型，它们的基因组序列尚不清晰。上述经过证实的基因预测表明了本发明对于工程菌代谢工程设计预测的准确性。这两个实验表明该方法对于序列不清楚的工程菌具有很好的应用性，这一点很好地解决了野生型菌株基因组尺度代谢网络模型难以应用到工程菌实验指导方面的问题。 

附图说明

图1 核黄素产量提高需改造的预测反应在不同代谢系统中的分布；其中carbon：碳水化合物代谢；Amino：氨基酸代谢；Cofact：辅因子及维生素代谢；Fatty：脂肪酸代谢； Nucleo：核酸及核苷酸代谢； 

图2为预测影响特定代谢产物产量关键反应的通量域比较方法的流程图。

具体实施方式

下面结合实施例说明本发明，这里所述实施例的方案不限制本发明，本领域的专业人员按照本发明的精神可以对其进行改进和变化，所述的这些改进和变化都应视为在本发明的范围内，本发明的范围和实质由权利要求来限定。本发明所用到的原料、试剂均有市售。枯草芽孢杆菌可以在中国微生物菌种保藏管理委员会普通微生物保藏中心获得，编号:1.3376。 

实施例1

一种预测影响特定代谢产物产量关键反应的通量域比较方法，该方法由两次计算，一次比较、一次排序组成，按如下的步骤进行：

（1）计算野生型目标物种自然生长状态的通量域：

首先采用通量平衡分析的方法，利用已知的基因组尺度代谢网络模型，计算生物量合成速率的最大值，计算时参数设置方法如下：根据已知的目标物种的培养条件及表型，将模型中培养基组成、溶氧、目标产物产量、其他副产物产量按照实验条件进行相同的设置；目标函数设定为生物量合成反应，计算得到生物量合成速率最优值f；之后，采用通量可变性分析的方法，得到野生型物种的代谢网络通量域分布，计算时参数设置方法如下：将生物量合成反应的速率设定为上一步骤中得到的最大值f，其它条件与上一步骤相同，计算得到两个向量V_{w_min}和V_{w_max}，其中V_{w_min}表示代谢网络中各反应通量值的最小值（即通量下限）组成的向量，V_{w_max}表示代谢网络中各反应通量值的最大值（即通量上限）组成的向量；计算所得的通量下限与通量上限组成的区域即为野生型物种的代谢通量域分布；

（2）模拟优化型物种的生长状态：

第一，采用通量平衡分析的方法，利用已知的基因组尺度代谢网络模型，计算生物量合成速率的最大值，计算时参数设置方法如下：根据已知的目标物种的培养条件及表型，将模型中培养基组成、溶氧按照实验条件进行相同的设置，而对目标产物产量和其他副产物的产量不作限制，生物量合成反应的通量设置为步骤（1）中生物量合成反应的最优值f，目标函数设定为特定代谢产物的产量，进行第二次模拟，计算得到特定产物产量的最大值p。

第二，采用通量可变性分析的方法，得到优化型物种的代谢网络通量域分布，计算时参数设置方法如下：将特定代谢产物的产量设定为上一步骤中得到的最大值p，其它条件与上一步骤相同，计算得到两个向量V_{e_min}和V_{e_max}，其中V_{e_min}表示代谢网络中各反应通量值的最小值（即通量下限）组成的向量，V_{e_max}表示代谢网络中各反应通量值的最大值（即通量上限）组成的向量；计算所得的通量下限与通量上限组成的区域即为优化型物种的代谢通量域分布。 

（3）一次比较：将两次计算的通量域进行比较，确定优化型与野生型之间反应通量域的不同，从而确定将实验型改造成优化型所需进行的代谢工程改造，进而制定相应的湿实验策略：比较及预测方法如下： 

对于代谢网络中的每一个反应，

其中，v_{w_max}表示野生型物种模拟得到的通量上限，v_{w_min}表示野生型物种模拟得到的通量下限，v_{e_max}表示优化型物种模拟得到的通量上限，v_{e_min}表示优化型物种模拟得到的通量下限；

上述预测结果有一个特殊情况，当(v_{e_max}-v_{w_max})/v_{w_max}与(v_{e_max}-v_{w_max})/v_{w_max}计算得到的结果符号相反时，预测结果为未知，即无法得到准确的预测结果；

（4）一次排序：对于每一个预测得到的反应，将其代谢通量分别设置为野生型通量的上限或下限，这样有两种设置方式，以v_l表示反应通量的下限，v_u表示反应通量的上限，则这两种设置方式可以表示为：

v_l=v_u=v_{w_min}   和  v_l=v_u= v_{w_max}

其他参数的设置与步骤（2）相同，每次设定后都以特定产物的产量为目标函数进行计算，计算方法为通量平衡分析，两种条件的设定各得到一个计算结果，取两个计算结果的最大值，记为q，计算该反应通量改变后特定产物产量的改变量C，即

将所有预测反应的C值按照从大到小的顺序进行排列，即为预测反应对特定产物产量影响大小的排序，根据该排序对代谢工程基因操作的实验策略进行指导。

实施例2 

对生物湿实验的实际指导：

1.对枯草芽孢杆菌生产核黄素预测的结果表明共有112个反应需进行改造，其中有80个需高表达的反应，27个需低表达的反应，另有5个反应由于其通量域的上限和下限改变的方向相反，因而不可预测。其中prs^[1],purFMNHD^[2] 和 ribABGHT^[3,4]基因对应反应的高表达预测在之前的研究中得到了证实。而其他尚未经实验证实的高表达及地表达反应可作为进一步实验的研究对象。

2.在反应对核黄素产量影响的排序中，由PRPP合成GTP的途径中包含的其他pur基因（purQLCEK）和ndk基因被预测为过表达，且对核黄素产量的影响最大，另外，一些合成谷氨酸盐或甘氨酸反应的过表达对核黄素的影响非常大，而消耗谷氨酸盐或甘氨酸反应的低表达同样对核黄素的产量非常大，这一点表明了谷氨酸盐和甘氨酸在核黄素合成过程中有非常重要的作用。在湿实验研究中可根据排序结果按照从高到低的顺序针对关键反应进行相关的实验设计。 

实施例3 

对比试验

常规的针对野生型物种通过最优化的方法（如通量平衡分析）进行模拟和分析，

方法： 通量平衡分析

步骤：利用已知的基因组尺度代谢网络模型，计算生物量合成速率的最大值。计算时参数设置方法如下：根据已知的目标物种的培养条件及表型，将模型中培养基组成、溶氧、目标产物产量、其他副产物产量及生物量合成速率按照实验条件进行相同的设置；目标函数设定为特定代谢产物的产量。

结果：得到野生型物种在实验条件下的通量分布，没有与优化型物种的比较，因而对实验设计参考价值不大。 

本发明的通量域比较方法： 

方法： 通量域比较方法同实施例1。

结果：得到可以通过实验进行改进以提高目标代谢产物产量的反应列表、改进方法，以及这些反应对目标代谢产物产量的影响力排序。 

实施例4 

本发明预测方法的验证情况：

为了验证该方法在代谢工程菌上的应用按照按照文献“Transcriptome analysis guided metabolic engineering of Bacillus subtilis for riboflavin production”和文献“Increased production of riboflavin by metabolic engineering of the purine pathway in Bacillus subtilis”两篇文章的实验条件，并按照本专利介绍的方法，模拟预测了RH33和B.subtilis PK两株工程菌核黄素的生产情况，结果表明在实验给定的最小培养基条件下，对RH33的模拟得到了118个预测反应，其中77个为过表达，40个为低表达，1个为未知。对B.subtilis PK的模拟得到了102个预测反应，其中65个为过表达，35个为低表达，2个为未知。该方法对RH33和B.subtilis PK的模拟结果进一步印证了生物湿实验中prs和purFMNHD，以及ribABGHT的过表达对提高核黄素产量的作用。同时，预测结果表明由PRPP合成GTP的途径中包含的其他pur基因（purQLCEK）和ndk基因被预测为过表达，其对应反应对核黄素产量的影响最大。除此之外，大量未经实验验证的预测反应为进一步的实验提供了由针对性的指导。RH33和B.subtilis PK都是经过多次进化实验和得到的菌型，它们的基因组序列尚不清晰。上述经过证实的基因预测表明了本发明对于工程菌代谢工程设计预测的准确性。这两个实验表明该方法对于序列不清楚的工程菌具有很好的应用性，这一点很好地解决了野生型菌株基因组尺度代谢网络模型难以应用到工程菌实验指导方面的问题。

参考文献

[1] Shi, S., Chen, T., Zhang, Z., Chen, X., Zhao, X., 2009a. Transcriptome analysis guided metabolic engineering of Bacillus subtilis for riboflavin production. Metab Eng. 11, 243-52.

[2] Shi, S., Zhang, Z., Chen, X., Chen, T., Zhao, X., 2009b. Increased production of riboflavin by metabolic engineering of the purine pathway in Bacillus subtilis. Biochemical Engineering Journal. 46, 28-33.

[3] Hümbelin M, Griesser V, Keller T, Schurter W, Haiker M, Hohmann HP, Ritz H, Richter G, Bacher A, 1999. GTP cyclohydrolase II and 3,4-dihydroxy-2-butanone 4-phosphate synthase are ratelimiting enzymes in riboflavin synthesis of an industrial Bacillus subtilis strain used for riboflavin production. J Ind Microbiol Biotech. 22, 1-7.

[4] Li XJ, Zhou SQ, Chen T, Chen X, Zhao XM, 2005. Effect of transcriptional modified riboflavin operon on riboflavin biosynthetic ability in Bacillus subtilis. 2nd Chinese National Chemical and Biochemical Engineering Annual Meeting. Beijing, China.。