import cobra
from cobra.io import load_model
# "iJO1366" and "salmonella" are also valid arguments
model = load_model("textbook")1 快速上手指南
1.1 加载模型并检查它
我们首先从两个最简单的内置模型开始。它们是沙门氏菌的模型，大肠杆菌的模型，以及一个叫做“教材”（textbook）的模型。textbook 模型描述了大肠杆菌的核心代谢。
cobrapy模型里的反应、代谢物和基因是一种叫做 “cobra.DictList”的特殊列表, 并且每一个分别以 cobra.Reaction, cobra.Metabolite 和cobra.Gene为对象。
print(len(model.reactions))
print(len(model.metabolites))
print(len(model.genes))95
72
137当我们使用Jupyter notebook时，这些信息将会以表格形式呈现。
model| Name | e_coli_core |
| Memory address | 124eb3a3050 |
| Number of metabolites | 72 |
| Number of reactions | 95 |
| Number of genes | 137 |
| Number of groups | 0 |
| Objective expression | 1.0*Biomass_Ecoli_core - 1.0*Biomass_Ecoli_core_reverse_2cdba |
| Compartments | cytosol, extracellular |
就像常规列表一样， DictList中的对象可以通过索引检索。例如，要获得该模型中的第30个反应（ (at index 29 because of 0-indexing):
model.reactions[29]| Reaction identifier | EX_glu__L_e |
| Name | L-Glutamate exchange |
| Memory address | 0x124ecc53fd0 |
| Stoichiometry | glu__L_e --> L-Glutamate --> |
| GPR | |
| Lower bound | 0.0 |
| Upper bound | 1000.0 |
此外，可以使用DictList.get_by_id() 函数按项目的id检索项目。例如，要获得胞质atp代谢产物对象（id为“atp_c”），我们可以执行以下操作：
model.metabolites.get_by_id("atp_c")| Metabolite identifier | atp_c |
| Name | ATP |
| Memory address | 0x124ecd18c90 |
| Formula | C10H12N5O13P3 |
| Compartment | c |
| In 13 reaction(s) | PPS, PYK, Biomass_Ecoli_core, PPCK, GLNS, PFK, ATPS4r, PGK, SUCOAS, GLNabc, ACKr, ATPM, ADK1 |
作为额外的收获，具有交互式shell（如IPython）的用户将能够完成制表，列出列表中的元素。虽然这不是大多数代码的推荐行为，因为id中可能存在“-”等字符，但在交互式提示中这非常有用：
model.reactions.EX_glc__D_e.bounds(-10.0, 1000.0)1.2 反应
我们将以葡萄糖-6-磷酸异构酶参与的葡萄糖 6-磷酸和果糖 6-磷酸相互转化反应作为演示对象。在我们的测试的模型中，该反应的反应 ID 是 PGI。
pgi = model.reactions.get_by_id("PGI")
pgi| Reaction identifier | PGI |
| Name | glucose-6-phosphate isomerase |
| Memory address | 0x124ecca7c10 |
| Stoichiometry | g6p_c <=> f6p_c D-Glucose 6-phosphate <=> D-Fructose 6-phosphate |
| GPR | b4025 |
| Lower bound | -1000.0 |
| Upper bound | 1000.0 |
我们可以将全名和催化的反应视为字符串
print(pgi.name)
print(pgi.reaction)glucose-6-phosphate isomerase
g6p_c <=> f6p_c因为 pgi 是可逆的，所以我们可以通过pgi.lower_bound < 0 和 pgi.upper_bound> 0来查看该反应的下限和上限
print(pgi.lower_bound, "< pgi <", pgi.upper_bound)
print(pgi.reversibility)-1000.0 < pgi < 1000.0
True反应的下限和上限被修改时，其可逆性属性也会自动更新。操作边界的首选方法是使用reaction.bounds，例如
old_bounds = pgi.bounds
pgi.bounds = (0, 1000.0)
print(pgi.lower_bound, "< pgi <", pgi.upper_bound)
print("Reversibility after modification:", pgi.reversibility)
pgi.bounds = old_bounds
print("Reversibility after resetting:", pgi.reversibility)0 < pgi < 1000.0
Reversibility after modification: False
Reversibility after resetting: True还可以使用reaction.lower_bound 或 reaction.upper_bound一次修改一个边界。这种方法没有使用 reaction.bounds同时设置两个边界的方法实用，因为用户可能会不小心设置下限高于上限（反之亦然）。如果出现下限高于上限（反之亦然），则会导致错误。
old_bounds = pgi.bounds
print('Upper bound prior to setting new lower bound:', pgi.upper_bound)
pgi.lower_bound = 1100
print('Upper bound after setting new lower bound:', pgi.upper_bound)
pgi.bounds = old_boundsUpper bound prior to setting new lower bound: 1000.0ValueError: The lower bound must be less than or equal to the upper bound (1100 <= 1000.0).我们需要确保反应是质量平衡的。如果反应不满足质量平衡，则该函数会返回元素；如果它返回为空，则表示反应是质量平衡的。
pgi.check_mass_balance(){}为了添加代谢物，我们传入一个带有代谢物对象及其系数的 dict
pgi.add_metabolites({model.metabolites.get_by_id("h_c"): -1})
pgi.reaction'g6p_c + h_c <=> f6p_c'反应不再保持质量平衡
pgi.check_mass_balance(){'charge': -1.0, 'H': -1.0}我们可以去除代谢物，反应将再次平衡。
pgi.subtract_metabolites({model.metabolites.get_by_id("h_c"): -1})
print(pgi.reaction)
print(pgi.check_mass_balance())g6p_c <=> f6p_c
{}也可以从字符串构建反应。但是，在 执行此操作时必须小心，以确保反应 ID 与模型中的反应 ID 匹配。箭头的方向也用于更新上限和下限。
pgi.reaction = "g6p_c --> f6p_c + h_c + green_eggs + ham"unknown metabolite 'green_eggs' created
unknown metabolite 'ham' createdpgi.reaction'g6p_c --> f6p_c + green_eggs + h_c + ham'pgi.reaction = "g6p_c <=> f6p_c"
pgi.reaction'g6p_c <=> f6p_c'1.3 代谢产物
我们将以胞质 atp 作为我们的代谢物，它在我们的测试模型中的 id 为"atp_c"。
atp = model.metabolites.get_by_id("atp_c")
atp| Metabolite identifier | atp_c |
| Name | ATP |
| Memory address | 0x124ecd18c90 |
| Formula | C10H12N5O13P3 |
| Compartment | c |
| In 13 reaction(s) | PPS, PYK, Biomass_Ecoli_core, PPCK, GLNS, PFK, ATPS4r, PGK, SUCOAS, GLNabc, ACKr, ATPM, ADK1 |
我们可以直接把代谢物名称和隔室（在本例中为胞质溶胶）以字符串的形式打印出来。
print(atp.name)
print(atp.compartment)ATP
c我们可以看到 ATP 在我们的模型中是一种带电分子。
atp.charge-4我们也可以看到代谢物的化学式。
print(atp.formula)C10H12N5O13P3反应属性给出了使用给定代谢物的所有反应的frozenset 。我们可以用它来计算使用 atp 的反应数量。
len(atp.reactions)13葡萄糖 6-磷酸等代谢物将参与较少的反应。
model.metabolites.get_by_id("g6p_c").reactionsfrozenset({<Reaction Biomass_Ecoli_core at 0x124ecca4a50>,
<Reaction G6PDH2r at 0x124ecc50090>,
<Reaction GLCpts at 0x124ece28190>,
<Reaction PGI at 0x124ecca7c10>})1.4 基因
gene_reaction_rule是该反应活跃的基因要求的布尔值表示，Schellenberger et al 2011 Nature Protocols 6(9):1290-307.GPR作为GPR类存储在反应的GPR中。它的字符串表示形式存储为 Reaction 对象的gene_reaction_rule。
gpr = pgi.gpr
print(gpr)
gpr_string = pgi.gene_reaction_rule
print(gpr_string)b4025
b4025相应的基因对象也存在。这些对象由反应本身以及模型追踪。
pgi.genesfrozenset({<Gene b4025 at 0x124eccbd1d0>})pgi_gene = model.genes.get_by_id("b4025")
pgi_gene| Gene identifier | b4025 |
| Name | pgi |
| Memory address | 0x124eccbd1d0 |
| Functional | True |
| In 1 reaction(s) | PGI |
每个基因都追踪它催化的反应。
pgi_gene.reactionsfrozenset({<Reaction PGI at 0x124ecca7c10>})如有必要，改变gene_reaction_rule将创建新的基因对象并更新所有关系。
pgi.gene_reaction_rule = "(spam or eggs)"
pgi.genesfrozenset({<Gene eggs at 0x124eea1bed0>, <Gene spam at 0x124eea1bfd0>})pgi_gene.reactionsfrozenset()新创建的基因也能添加到模型中
model.genes.get_by_id("spam")| Gene identifier | spam |
| Name | |
| Memory address | 0x124eea1bfd0 |
| Functional | True |
| In 1 reaction(s) | PGI |
knock_out_model_genes函数能够评估 GPR，并在反应被敲除时将上限和下限设置为 0。
cobra.manipulation.knock_out_model_genes(
model, ["spam"])
print("after 1 KO: %4d < flux_PGI < %4d" % (pgi.lower_bound, pgi.upper_bound))
cobra.manipulation.knock_out_model_genes(
model, ["eggs"])
print("after 2 KO: %4d < flux_PGI < %4d" % (pgi.lower_bound, pgi.upper_bound))after 1 KO: -1000 < flux_PGI < 1000
after 2 KO: 0 < flux_PGI < 0当在上下文中敲除模型基因时，它在离开上下文时会相反。
很多时候，人们希望对模型进行一些小的更改并评估这些更改的影响。例如，我们可能希望按顺序敲除所有反应，看看这对目标函数有什么影响。一种方法是在每次敲除之前使用 model.copy() 创建模型的新副本。然而，即使使用小型模型，这也是一种非常缓慢的方法，因为模型是相当复杂的对象。更好的做法是进行敲除、优化，然后手动重置反应边界，然后再进行下一个反应。然而，由于这是一个非常常见的场景，cobrapy 允许我们将模型用作上下文，以自动恢复更改。
model = load_model('textbook')
for reaction in model.reactions[:5]:
with model as model:
reaction.knock_out()
model.optimize()
print('%s blocked (bounds: %s), new growth rate %f' %
(reaction.id, str(reaction.bounds), model.objective.value))ACALD blocked (bounds: (0, 0)), new growth rate 0.873922
ACALDt blocked (bounds: (0, 0)), new growth rate 0.873922
ACKr blocked (bounds: (0, 0)), new growth rate 0.873922
ACONTa blocked (bounds: (0, 0)), new growth rate -0.000000
ACONTb blocked (bounds: (0, 0)), new growth rate -0.000000如果我们看一下那些被敲击的反应，就会发现它们的界限都被恢复了。
[reaction.bounds for reaction in model.reactions[:5]][(-1000.0, 1000.0),
(-1000.0, 1000.0),
(-1000.0, 1000.0),
(-1000.0, 1000.0),
(-1000.0, 1000.0)]上下文也支持被嵌套
print('original objective: ', model.objective.expression)
with model:
model.objective = 'ATPM'
print('print objective in first context:', model.objective.expression)
with model:
model.objective = 'ACALD'
print('print objective in second context:', model.objective.expression)
print('objective after exiting second context:',
model.objective.expression)
print('back to original objective:', model.objective.expression)original objective: 1.0*Biomass_Ecoli_core - 1.0*Biomass_Ecoli_core_reverse_2cdba
print objective in first context: 1.0*ATPM - 1.0*ATPM_reverse_5b752
print objective in second context: 1.0*ACALD - 1.0*ACALD_reverse_fda2b
objective after exiting second context: 1.0*ATPM - 1.0*ATPM_reverse_5b752
back to original objective: 1.0*Biomass_Ecoli_core - 1.0*Biomass_Ecoli_core_reverse_2cdba大多数修改模型的方法都像这样支持，包括添加和去除反应和代谢物以及设定目标。支持的方法和函数在相应的文档中提到了这一点。
虽然它没有任何实际效果，但为了语法上的方便，也可以使用与上下文外部不同的名称来引用模型。如
with model as inner:
inner.reactions.PFK.knock_out