3  构建模型

3.1 模型、反应和代谢物

这个简单示例演示了如何创建模型，创建反应，然后将反应添加到模型中。

我们将使用 STM_1.0 模型中的“3OAS140”反应：

1.0 malACP[c] + 1.0 h[c] + 1.0 ddcaACP[c] → 1.0 co2[c] + 1.0 ACP[c] + 1.0 3omrsACP[c]

首先，创建模型和反应。

from cobra import Model, Reaction, Metabolite

model = Model('example_model')

reaction = Reaction('R_3OAS140')
reaction.name = '3 oxoacyl acyl carrier protein synthase n C140 '
reaction.subsystem = 'Cell Envelope Biosynthesis'
reaction.lower_bound = 0.  # This is the default
reaction.upper_bound = 1000.  # This is the default

我们也需要创造代谢物。如果我们使用现有模型“Model.get_by_id”来获得适当的代谢物对象。

ACP_c = Metabolite(
    'ACP_c',
    formula='C11H21N2O7PRS',
    name='acyl-carrier-protein',
    compartment='c')
omrsACP_c = Metabolite(
    'M3omrsACP_c',
    formula='C25H45N2O9PRS',
    name='3-Oxotetradecanoyl-acyl-carrier-protein',
    compartment='c')
co2_c = Metabolite('co2_c', formula='CO2', name='CO2', compartment='c')
malACP_c = Metabolite(
    'malACP_c',
    formula='C14H22N2O10PRS',
    name='Malonyl-acyl-carrier-protein',
    compartment='c')
h_c = Metabolite('h_c', formula='H', name='H', compartment='c')
ddcaACP_c = Metabolite(
    'ddcaACP_c',
    formula='C23H43N2O8PRS',
    name='Dodecanoyl-ACP-n-C120ACP',
    compartment='c')

旁注：SId

强烈建议反应、代谢物和基因的 ID 是有效的 SBML 标识符 （SId）。SId 是从基本 XML 类型字符串派生的数据类型，但对允许的字符以及这些字符可能出现的顺序有限制 letter ::= ’a’..’z’,’A’..’Z’ digit ::= ’0’..’9’ idChar ::= letter | digit | ’’ SId ::= ( letter | ’’ ) idCharidChar 主要限制是 id 不能以数字开头。使用 SIds 允许序列化到 SBML。此外，通过点语法进行代码补全和对象访问等功能将在 cobrapy 中工作。

使用代谢物及其化学计量系数的字典可将代谢物添加到反应中。一组代谢物可以一次全部添加，也可以一次添加一个。

reaction.add_metabolites({
    malACP_c: -1.0,
    h_c: -1.0,
    ddcaACP_c: -1.0,
    co2_c: 1.0,
    ACP_c: 1.0,
    omrsACP_c: 1.0
})

reaction.reaction  # This gives a string representation of the reaction

'2.0 ddcaACP_c + 2.0 h_c + 2.0 malACP_c --> 2.0 ACP_c + 2.0 M3omrsACP_c + 2.0 co2_c'

gene_reaction_rule是该反应活跃的基因要求的布尔表示，如 Schellenberger 等人 的研究Schellenberger et al 2011 Nature Protocols 6(9):1290-307. 我们将分配基因反应规则字符串，该字符串将自动创建相应的基因对象。

reaction.gene_reaction_rule = '( STM2378 or STM1197 )'
reaction.genes

frozenset({<Gene STM1197 at 0x1a2c2619f90>, <Gene STM2378 at 0x1a2c261a150>})

此时，模型仍为空

print(f'{len(model.reactions)} reactions initially')
print(f'{len(model.metabolites)} metabolites initially')
print(f'{len(model.genes)} genes initially')

3 reactions initially
8 metabolites initially
2 genes initially

我们将反应添加到模型中，该模型还将添加所有相关的代谢物和基因

model.add_reactions([reaction])

# The objects have been added to the model
print(f'{len(model.reactions)} reactions')
print(f'{len(model.metabolites)} metabolites')
print(f'{len(model.genes)} genes')

Ignoring reaction 'R_3OAS140' since it already exists.

3 reactions
8 metabolites
2 genes

我们可以遍历模型对象来观察内容

# Iterate through the the objects in the model
print("Reactions")
print("---------")
for x in model.reactions:
    print("%s : %s" % (x.id, x.reaction))

print("")
print("Metabolites")
print("-----------")
for x in model.metabolites:
    print('%9s : %s' % (x.id, x.formula))

print("")
print("Genes")
print("-----")
for x in model.genes:
    associated_ids = (i.id for i in x.reactions)
    print("%s is associated with reactions: %s" %
          (x.id, "{" + ", ".join(associated_ids) + "}"))

Reactions
---------
R_3OAS140 : 2.0 ddcaACP_c + 2.0 h_c + 2.0 malACP_c --> 2.0 ACP_c + 2.0 M3omrsACP_c + 2.0 co2_c
EX_co2_e : co2_e <=> 
SK_glycogen_c : glycogen_c <=> 

Metabolites
-----------
 malACP_c : C14H22N2O10PRS
      h_c : H
ddcaACP_c : C23H43N2O8PRS
    co2_c : CO2
    ACP_c : C11H21N2O7PRS
M3omrsACP_c : C25H45N2O9PRS
glycogen_c : None
    co2_e : None

Genes
-----
STM1197 is associated with reactions: {R_3OAS140}
STM2378 is associated with reactions: {R_3OAS140}

3.2 目标

最后，我们需要设定模型的目标。在这里，我们只希望这是我们添加的单个反应中通量的最大化，我们通过将反应的标识符分配给模型的“目标”属性来做到这一点。

model.objective = 'R_3OAS140'

创建的目标是一个符号代数表达式，我们可以通过打印来检查它

print(model.objective.expression)
print(model.objective.direction)

1.0*R_3OAS140 - 1.0*R_3OAS140_reverse_60acb
max

这里表明求解器将在正向上最大化通量。

3.3 模型验证

为了与其他工具进行交换，您可以验证模型并将其导出到 SBML。 有关序列化和可用格式的详细信息，请参阅“读取和写入模型”部分

import tempfile
from pprint import pprint
from cobra.io import write_sbml_model, validate_sbml_model
with tempfile.NamedTemporaryFile(suffix='.xml') as f_sbml:
    write_sbml_model(model, filename=f_sbml.name)
    report = validate_sbml_model(filename=f_sbml.name)

pprint(report)

(None,
 {'COBRA_CHECK': [],
  'COBRA_ERROR': ['No SBML model detected in file.'],
  'COBRA_FATAL': [],
  'COBRA_WARNING': [],
  'SBML_ERROR': ['E0 (Error): Operating system (core, L1); File unreadable; '
                 'File unreadable. '
                 'C:\\Users\\ADMINI~1\\AppData\\Local\\Temp\\tmpc4e0sees.xml\n'],
  'SBML_FATAL': [],
  'SBML_SCHEMA_ERROR': [],
  'SBML_WARNING': []})

该模型有效，没有 COBRA 或 SBML 错误或警告。

3.4 交换（exchange）、汇（sinks）和需求（demands）

可以使用模型的“add_boundary”方法添加边界反应 有三种不同类型的预定义边界反应：交换反应、需求反应和汇反应。它们都是不平衡的伪反应，这意味着它们通过在模型系统中添加或去除代谢物来实现建模功能，但不是基于真实的生物学。交换反应是一种可逆反应，在细胞外区室中增加或去除细胞外代谢物。需求反应是消耗细胞内代谢物的不可逆反应。汇类似于交换，但专门用于细胞内代谢物，即添加或去除细胞内代谢物的可逆反应。

print("exchanges", model.exchanges)
print("demands", model.demands)
print("sinks", model.sinks)

exchanges [<Reaction EX_co2_e at 0x1a2c266f190>]
demands []
sinks [<Reaction SK_glycogen_c at 0x1a2c264de50>]

边界反应在代谢物上定义。首先，我们在模型中添加两种代谢物，然后 我们定义边界反应。我们将糖原添加到胞质区室“c”中，将 CO2 添加到外部区室“e”。

model.add_metabolites([
    Metabolite(
    'glycogen_c',
    name='glycogen',
    compartment='c'
    ),
    Metabolite(
    'co2_e',
    name='CO2',
    compartment='e'
    ),
])

# create exchange reaction
model.add_boundary(model.metabolites.get_by_id("co2_e"), type="exchange")

ValueError: Boundary reaction 'EX_co2_e' already exists.

# create exchange reaction
model.add_boundary(model.metabolites.get_by_id("glycogen_c"), type="sink")

ValueError: Boundary reaction 'SK_glycogen_c' already exists.

# Now we have an additional exchange and sink reaction in the model
print("exchanges", model.exchanges)
print("sinks", model.sinks)
print("demands", model.demands)

exchanges [<Reaction EX_co2_e at 0x1a2c266f190>]
sinks [<Reaction SK_glycogen_c at 0x1a2c264de50>]
demands []

若要创建需求反应而不是接收器，请使用类型“demand”而不是“sink”。

有关所有边界反应的信息都可以通过模型的属性“边界”获得。

# boundary reactions
model.boundary

[<Reaction EX_co2_e at 0x1a2c266f190>,
 <Reaction SK_glycogen_c at 0x1a2c264de50>]

获得所有代谢反应的一个巧妙技巧是

# metabolic reactions
set(model.reactions) - set(model.boundary)

{<Reaction R_3OAS140 at 0x1a2c2602b90>}