from cobra.io import load_model
= load_model("textbook") model
5 使用FBA进行模拟
使用磁通平衡分析的仿真可以使用 Model.optimize() 进行求解。这将最大化或最小化(最大化是默认设置)通过客观反应的通量。
5.1 运行 FBA
= model.optimize()
solution print(solution)
<Solution 0.874 at 0x195d0067050>
Model.optimize()
方法会给出一个解决方案(Solution)对象。这个对象包含了几个关键属性:
objective_value
:表示目标函数的优化值,也就是模型追求最大化或最小化的那个数值结果。status
:显示线性规划求解器给出的解决状态,告诉你问题是否得到解决、有没有特殊情况需要留意。fluxes
:这是一个按照反应标识符来组织的pandas序列,里面记录了每个反应的通量值。通量是指正向反应速率与逆向反应速率之差,直观反映了物质在该反应中的净流动情况。shadow_prices
:同样是一个pandas序列,不过这次是依据代谢物标识符来排列的,展示了每个代谢物的影子价格。影子价格揭示了在保持目标函数最优的前提下,如果系统可利用的某种代谢物增多或减少一个单位,目标函数值的变化量,有助于理解代谢物的重要性及限制条件。
例如,在最后一次调用“model.optimize()”之后,如果优化成功,则其状态将为最佳。如果模型不可行,则会引发错误。
solution.objective_value
0.8739215069684279
可与 cobrapy 一起使用的求解器速度非常快,以至于对于许多中小型模型来说,计算求解的速度甚至比从求解器收集值并将 python 对象转换为它们所需的速度还要快。使用“model.optimize”,我们收集所有反应和代谢物的值,如果重复进行,这可能需要大量时间。如果我们只对单个反应或目标的通量值感兴趣,那么使用“model.slim_optimize”会更快,它只执行优化并返回目标值,由你来获取你可能需要的其他值。
%%time
model.optimize().objective_value
CPU times: total: 0 ns
Wall time: 1e+03 µs
0.8739215069684305
%%time
model.slim_optimize()
CPU times: total: 0 ns
Wall time: 1e+03 µs
0.8739215069684305
5.1.1 分析 FBA solutions
使用FBA求解的模型可以通过使用摘要方法进一步分析,这些方法输出打印文本以快速表示模型行为。对整个模型调用 summary 方法将显示有关模型的输入和输出行为的信息,以及优化的目标。
model.summary()
Objective
1.0 Biomass_Ecoli_core = 0.8739215069684297
Uptake
Metabolite | Reaction | Flux | C-Number | C-Flux |
---|---|---|---|---|
glc__D_e | EX_glc__D_e | 10 | 6 | 100.00% |
nh4_e | EX_nh4_e | 4.765 | 0 | 0.00% |
o2_e | EX_o2_e | 21.8 | 0 | 0.00% |
pi_e | EX_pi_e | 3.215 | 0 | 0.00% |
5.1.1.1 Secretion
Metabolite | Reaction | Flux | C-Number | C-Flux |
---|---|---|---|---|
co2_e | EX_co2_e | -22.81 | 1 | 100.00% |
h2o_e | EX_h2o_e | -29.18 | 0 | 0.00% |
h_e | EX_h_e | -17.53 | 0 | 0.00% |
此外,还可以使用汇总方法检查单个代谢物的输入输出行为。例如,以下命令可用于检查模型的整体氧化还原平衡
model.metabolites.nadh_c.summary()
nadh_c
C21H27N7O14P2
Producing Reactions
Percent | Flux | Reaction | Definition |
---|---|---|---|
13.14% | 5.064 | AKGDH | akg_c + coa_c + nad_c --> co2_c + nadh_c + succoa_c |
8.04% | 3.1 | Biomass_Ecoli_core | 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.361 e4p_c + 0.0709 f6p_c + 0.129 g3p_c + 0.205 g6p_c + 0.2557 gln__L_c + 4.9414 glu__L_c + 59.81 h2o_c + 3.547 nad_c + 13.0279 nadph_c + 1.7867 oaa_c + 0.5191 pep_c + 2.8328 pyr_c + 0.8977 r5p_c --> 59.81 adp_c + 4.1182 akg_c + 3.7478 coa_c + 59.81 h_c + 3.547 nadh_c + 13.0279 nadp_c + 59.81 pi_c |
41.58% | 16.02 | GAPD | g3p_c + nad_c + pi_c <=> 13dpg_c + h_c + nadh_c |
13.14% | 5.064 | MDH | mal__L_c + nad_c <=> h_c + nadh_c + oaa_c |
24.09% | 9.283 | PDH | coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c |
5.1.1.2 Consuming Reactions
Percent | Flux | Reaction | Definition |
---|---|---|---|
100.00% | -38.53 | NADH16 | 4.0 h_c + nadh_c + q8_c --> 3.0 h_e + nad_c + q8h2_c |
或者了解主要的能源生产和消费反应
model.metabolites.atp_c.summary()
atp_c
C10H12N5O13P3
Producing Reactions
Percent | Flux | Reaction | Definition |
---|---|---|---|
66.58% | 45.51 | ATPS4r | adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c |
23.44% | 16.02 | PGK | 3pg_c + atp_c <=> 13dpg_c + adp_c |
2.57% | 1.758 | PYK | adp_c + h_c + pep_c --> atp_c + pyr_c |
7.41% | 5.064 | SUCOAS | atp_c + coa_c + succ_c <=> adp_c + pi_c + succoa_c |
5.1.1.3 Consuming Reactions
Percent | Flux | Reaction | Definition |
---|---|---|---|
12.27% | -8.39 | ATPM | atp_c + h2o_c --> adp_c + h_c + pi_c |
76.46% | -52.27 | Biomass_Ecoli_core | 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.361 e4p_c + 0.0709 f6p_c + 0.129 g3p_c + 0.205 g6p_c + 0.2557 gln__L_c + 4.9414 glu__L_c + 59.81 h2o_c + 3.547 nad_c + 13.0279 nadph_c + 1.7867 oaa_c + 0.5191 pep_c + 2.8328 pyr_c + 0.8977 r5p_c --> 59.81 adp_c + 4.1182 akg_c + 3.7478 coa_c + 59.81 h_c + 3.547 nadh_c + 13.0279 nadp_c + 59.81 pi_c |
0.33% | -0.2235 | GLNS | atp_c + glu__L_c + nh4_c --> adp_c + gln__L_c + h_c + pi_c |
10.94% | -7.477 | PFK | atp_c + f6p_c --> adp_c + fdp_c + h_c |
5.2 更改目标
目标函数由客观反应的objective_coefficient属性确定。通常,使用描述构成细胞的代谢物组成的“生物量”函数。
= model.reactions.get_by_id("Biomass_Ecoli_core") biomass_rxn
目前在模型中,目标中只有一个反应(生物质反应),线性系数为1。
from cobra.util.solver import linear_reaction_coefficients
linear_reaction_coefficients(model)
{<Reaction Biomass_Ecoli_core at 0x195cff33490>: 1.0}
可以通过分配 Model.objective 来更改目标函数,Model.objective 可以是反应对象(或只是它的名称),也可以是 {Reaction: objective_coefficient} 的字典。
# 将目标改为 ATPM
= "ATPM"
model.objective
# 上界应该是1000,这样我们才能得到实际的最佳值
"ATPM").upper_bound = 1000.
model.reactions.get_by_id( linear_reaction_coefficients(model)
{<Reaction ATPM at 0x195cff09450>: 1.0}
model.optimize().objective_value
175.00000000000006
我们还可以有更复杂的目标,比如二次项。
5.3 运行 FVA
FBA不会总是给出独特的解决方案,因为多个通量状态可以达到相同的最佳状态。FVA(或通量变异性分析)找到每个代谢通量处于最佳状态的范围。
import cobra
from cobra.flux_analysis import flux_variability_analysis
10]) flux_variability_analysis(model, model.reactions[:
minimum | maximum | |
---|---|---|
ACALD | -9.375513e-15 | 0.000000e+00 |
ACALDt | -1.198454e-14 | 0.000000e+00 |
ACKr | -1.832929e-14 | 1.039260e-14 |
ACONTa | 2.000000e+01 | 2.000000e+01 |
ACONTb | 2.000000e+01 | 2.000000e+01 |
ACt2r | -3.615752e-15 | 0.000000e+00 |
ADK1 | 0.000000e+00 | 1.924093e-13 |
AKGDH | 2.000000e+01 | 2.000000e+01 |
AKGt2r | -9.821844e-15 | 0.000000e+00 |
ALCD2x | 9.994214e-15 | 0.000000e+00 |
设置参数“fraction_of_optimium=0.90”将给出 90% 最佳反应的通量范围。
cobra.flux_analysis.flux_variability_analysis(10], fraction_of_optimum=0.9) model, model.reactions[:
minimum | maximum | |
---|---|---|
ACALD | -2.692308 | 0.000000e+00 |
ACALDt | -2.692308 | 0.000000e+00 |
ACKr | -4.117647 | 1.039260e-14 |
ACONTa | 8.461538 | 2.000000e+01 |
ACONTb | 8.461538 | 2.000000e+01 |
ACt2r | -4.117647 | 0.000000e+00 |
ADK1 | 0.000000 | 1.750000e+01 |
AKGDH | 2.500000 | 2.000000e+01 |
AKGt2r | -1.489362 | 0.000000e+00 |
ALCD2x | -2.333333 | 0.000000e+00 |
标准FVA可能包含环路,即高绝对通量值,只有在允许它们参与环路时才能很高(一种在体内不会发生的数学伪影)。使用“loopless”参数来避免此类循环。接下来,我们可以在环中看到 FRD7 和 SUCDi 反应,与此同时可使用无环 FVA 来避免这种情况。
= [model.reactions.FRD7, model.reactions.SUCDi]
loop_reactions =loop_reactions, loopless=False) flux_variability_analysis(model, reaction_list
minimum | maximum | |
---|---|---|
FRD7 | 0.0 | 980.0 |
SUCDi | 20.0 | 1000.0 |
=loop_reactions, loopless=True) flux_variability_analysis(model, reaction_list
minimum | maximum | |
---|---|---|
FRD7 | 0.0 | 0.0 |
SUCDi | 20.0 | 20.0 |
5.3.1 在摘要方法中运行 FVA
通量变异性分析也可以嵌入到对摘要方法的调用中。例如,可以通过以下方式快速找到基材消耗和产物形成的预期变化
model.optimize()=0.95) model.summary(fva
Objective
1.0 ATPM = 175.0000000000001
Uptake
Metabolite | Reaction | Flux | Range | C-Number | C-Flux |
---|---|---|---|---|---|
glc__D_e | EX_glc__D_e | 10 | [9.5; 10] | 6 | 100.00% |
o2_e | EX_o2_e | 60 | [55.88; 60] | 0 | 0.00% |
5.3.1.1 Secretion
Metabolite | Reaction | Flux | Range | C-Number | C-Flux |
---|---|---|---|---|---|
ac_e | EX_ac_e | 0 | [-2.059; 0] | 2 | 0.00% |
acald_e | EX_acald_e | 0 | [-1.346; 0] | 2 | 0.00% |
akg_e | EX_akg_e | 0 | [-0.7447; 0] | 5 | 0.00% |
co2_e | EX_co2_e | -60 | [-60; -54.17] | 1 | 100.00% |
etoh_e | EX_etoh_e | 0 | [-1.167; 0] | 2 | 0.00% |
for_e | EX_for_e | 0 | [-5.833; 0] | 1 | 0.00% |
glu__L_e | EX_glu__L_e | 0 | [-0.6731; 0] | 5 | 0.00% |
h2o_e | EX_h2o_e | -60 | [-60; -54.17] | 0 | 0.00% |
h_e | EX_h_e | 0 | [-5.833; 0] | 0 | 0.00% |
lac__D_e | EX_lac__D_e | 0 | [-1.129; 0] | 3 | 0.00% |
nh4_e | EX_nh4_e | 0 | [0; 0.6731] | 0 | 0.00% |
pi_e | EX_pi_e | 0 | [0; 0.171] | 0 | 0.00% |
pyr_e | EX_pyr_e | 0 | [-1.346; 0] | 3 | 0.00% |
succ_e | EX_succ_e | 0 | [-0.875; 0] | 4 | 0.00% |
同样,代谢物质量平衡的变异性也可以通过通量变异性分析来检查。
=0.95) model.metabolites.pyr_c.summary(fva
pyr_c
C3H3O3
Producing Reactions
Percent | Flux | Range | Reaction | Definition |
---|---|---|---|---|
50.00% | 10 | [9.5; 10] | GLCpts | glc__D_e + pep_c --> g6p_c + pyr_c |
0.00% | 0 | [-1.129; 0] | LDH_D | lac__D_c + nad_c <=> h_c + nadh_c + pyr_c |
0.00% | 0 | [0; 8.75] | ME1 | mal__L_c + nad_c --> co2_c + nadh_c + pyr_c |
0.00% | 0 | [0; 8.75] | ME2 | mal__L_c + nadp_c --> co2_c + nadph_c + pyr_c |
50.00% | 10 | [1.25; 18.75] | PYK | adp_c + h_c + pep_c --> atp_c + pyr_c |
0.00% | 0 | [-1.346; 0] | PYRt2 | h_e + pyr_e <=> h_c + pyr_c |
5.3.1.2 Consuming Reactions
Percent | Flux | Range | Reaction | Definition |
---|---|---|---|---|
0.00% | 0 | [-0.1316; 0] | Biomass_Ecoli_core | 1.496 3pg_c + 3.7478 accoa_c + 59.81 atp_c + 0.361 e4p_c + 0.0709 f6p_c + 0.129 g3p_c + 0.205 g6p_c + 0.2557 gln__L_c + 4.9414 glu__L_c + 59.81 h2o_c + 3.547 nad_c + 13.0279 nadph_c + 1.7867 oaa_c + 0.5191 pep_c + 2.8328 pyr_c + 0.8977 r5p_c --> 59.81 adp_c + 4.1182 akg_c + 3.7478 coa_c + 59.81 h_c + 3.547 nadh_c + 13.0279 nadp_c + 59.81 pi_c |
100.00% | -20 | [-28.75; -13] | PDH | coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c |
0.00% | 0 | [-5.833; 0] | PFL | coa_c + pyr_c --> accoa_c + for_c |
0.00% | 0 | [-8.75; 0] | PPS | atp_c + h2o_c + pyr_c --> amp_c + 2.0 h_c + pep_c + pi_c |
在这些汇总方法中,这些值报告为中心点 +/- FVA 解的范围,由最大值和最小值计算得出。
5.4 运行 pFBA
简约的 FBA(通常写成 pFBA)找到一个通量分布,该分布可提供最佳增长率,但使通量总和最小化。这涉及求解两个顺序线性程序,但由 cobrapy 透明地处理。有关pFBA的更多详细信息,请参阅 Lewis et al. (2010).
= 'Biomass_Ecoli_core'
model.objective = model.optimize()
fba_solution = cobra.flux_analysis.pfba(model) pfba_solution
这些函数可能会给出完全不同的目标值,因为 pFBA 显示的目标值定义为 ‘sum(abs(pfba_solution.fluxes.values))’,而标准 FBA 的目标值定义为通过优化反应的加权通量(例如 ‘fba_solution.fluxes[“Biomass_Ecoli_core”]’)。
pFBA 和 FBA 都应在求解器容差范围内为要优化的目标返回相同的结果。例如,对于使反应“Biomass_Ecoli_core”最大化的 FBA 问题:
abs(fba_solution.fluxes["Biomass_Ecoli_core"] - pfba_solution.fluxes[
"Biomass_Ecoli_core"])
5.551115123125783e-16
import numpy as np
np.isclose("Biomass_Ecoli_core"],
fba_solution.fluxes["Biomass_Ecoli_core"]
pfba_solution.fluxes[ )
True
5.5 运行几何FBA
几何 FBA 找到一个独特的最佳通量分布,这是可能通量范围的核心。有关几何FBA的更多详细信息,请参阅 K Smallbone, E Simeonidis (2009).
= cobra.flux_analysis.geometric_fba(model)
geometric_fba_sol geometric_fba_sol
fluxes | reduced_costs | |
---|---|---|
ACALD | 0.000000e+00 | 0.0 |
ACALDt | 0.000000e+00 | 0.0 |
ACKr | 1.260427e-13 | 0.0 |
ACONTa | 6.007250e+00 | 0.0 |
ACONTb | 6.007250e+00 | 0.0 |
... | ... | ... |
TALA | 1.496984e+00 | 0.0 |
THD2 | 0.000000e+00 | 0.0 |
TKT1 | 1.496984e+00 | 0.0 |
TKT2 | 1.181498e+00 | 0.0 |
TPI | 7.477382e+00 | 0.0 |
95 rows × 2 columns