Editing an existing model#

Curation is mostly editing: a stoichiometry that was wrong, a gene association that has been superseded, a reaction that should never have been there. Each of those edits has a safe procedure in both toolboxes.

Functions#

MATLAB

Python

Purpose

changeRxns

change_reaction_equations

replace a reaction's stoichiometry

changeGrRules

change_gene_reaction_rules

set or extend gene associations

standardizeGrRules

gpr_to_dnf

normalise a GPR

setParam

Reaction.bounds cobrapy

change bounds, objective

removeReactions

Model.remove_reactions cobrapy

delete reactions

removeMets

remove_metabolites

delete metabolites

removeGenes, deleteUnusedGenes

remove_genes

delete genes

Setup#

smallYeast.yml from docs/data/.

model = readYAMLmodel('smallYeast.yml');

from raven_toolbox.io import read_yaml_model

model = read_yaml_model("smallYeast.yml")

Fixing a stoichiometry#

Give the whole equation, not a coefficient: the toolbox rewrites the reaction's column in the stoichiometric matrix, which keeps everything consistent.

model = changeRxns(model, {'PGI'}, {'G6P_c <=> F6P_c'});
idx = getIndexes(model, 'PGI', 'rxns');
eqn = constructEquations(model, model.rxns(idx));
fprintf('%s\n', eqn{1});

alpha-D-glucose 6-phosphate[c] <=> beta-D-fructofuranose 6-phosphate[c]

from raven_toolbox.manipulation import change_reaction_equations

change_reaction_equations(model, {"PGI": "G6P_c <=> F6P_c"})
print(model.reactions.get_by_id("PGI").reaction)

G6P_c <=> F6P_c

changeRxns takes two matching lists, reactions and equations, and change_reaction_equations a mapping of reaction to equation, so a batch of curated reactions can be applied in one call, which is what a curation spreadsheet turns into. Every id must exist in the model; these functions change reactions and do not add them.

Bounds are left unchanged. Rewriting an equation from => to <=> does not make the reaction reversible, because lb and ub set the direction and the arrow does not; Changing bounds and the objective sets those separately. changeRxns takes the eqnType, compartment and allowNewMets options of addRxns, and change_reaction_equations the mets_by, compartment, allow_new_mets and new_met_prefix options of add_reactions_from_equations; both match metabolites by id by default. allowNewMets defaults to false and allow_new_mets to True.

Changing a gene association#

model = changeGrRules(model, {'PGI'}, {'YBR196C or YLR354C'});
disp(model.grRules{idx});

YBR196C or YLR354C

from raven_toolbox.manipulation import change_gene_reaction_rules

change_gene_reaction_rules(model, {"PGI": "YBR196C or YLR354C"})
pgi = model.reactions.get_by_id("PGI")
print(pgi.gene_reaction_rule)
print("genes now:", sorted(gene.id for gene in pgi.genes))

YBR196C or YLR354C
genes now: ['YBR196C', 'YLR354C']

A gene that is not in the model is created. Pass replace=False to append an isozyme to the rule, giving (old) or (new), which records an additional catalyst and not a correction; the default overwrites the rule. A reaction with no rule gets the new rule alone, without empty brackets.

Normalising a GPR#

Rules that mix and and or in nested brackets are hard to compare and hard to score. Both toolboxes rewrite them into disjunctive normal form, a list of alternative complexes. is_dnf takes a rule string; gpr_to_dnf takes cobrapy's parsed GPR object and returns the complexes as lists.

grRules = standardizeGrRules(model);
fprintf('%s\n', grRules{idx});

YBR196C or YLR354C

standardizeGrRules returns the rules and a matching rxnGeneMat, not a model, so the results have to be assigned back into the struct to persist. Disjunctive normal form is a flat or of and groups, one group per alternative complex, which is the form the scoring in Context-specific models expects.

from cobra.core.gene import GPR

from raven_toolbox.manipulation import gpr_to_dnf
from raven_toolbox.utils import is_dnf

rule = "YBR196C and (YLR354C or YGR192C)"
print("already DNF:", is_dnf(rule))
print(gpr_to_dnf(GPR.from_string(rule)))

already DNF: False
[['YBR196C', 'YLR354C'], ['YBR196C', 'YGR192C']]

Changing bounds and the objective#

model = setParam(model, 'lb', 'PGI', 0);        % make it irreversible
model = setParam(model, 'obj', 'biomassOUT', 1);
fprintf('bounds: [%g %g]\n', model.lb(idx), model.ub(idx));

bounds: [0 1000]

model.reactions.get_by_id("PGI").bounds = (0, 1000)
model.objective = "biomassOUT"
print("bounds:", model.reactions.get_by_id("PGI").bounds)

bounds: (0, 1000)

setParam takes 'lb', 'ub', 'eq', 'obj', 'rev', 'var' and 'unc', and accepts a list of reactions with a list of values, or one value for all, applying them in one call.

Deleting metabolites, reactions and genes#

Deletion is where the two designs differ most. RAVEN must remove the matching row or column from every field, so deletion has its own functions; cobrapy's objects record what they are attached to.

before = numel(model.rxns);
reduced = removeReactions(model, {'ACO'}, ...
    'removeUnusedMets', true, 'removeUnusedGenes', true);
fprintf('%d -> %d reactions, %d -> %d genes\n', before, numel(reduced.rxns), ...
    numel(model.genes), numel(reduced.genes));

53 -> 52 reactions, 61 -> 60 genes

Both flags default to false, so a plain removeReactions leaves metabolites and genes that nothing refers to, which later appear as gaps. removeUnusedComps does the same for compartments. deleteUnusedGenes performs the gene half on its own.

before = len(model.reactions), len(model.genes)

with model:
    model.remove_reactions([model.reactions.get_by_id("ACO")], remove_orphans=True)
    print(f"{before[0]} -> {len(model.reactions)} reactions, "
          f"{before[1]} -> {len(model.genes)} genes")

53 -> 52 reactions, 61 -> 60 genes

remove_orphans=True is cobrapy's equivalent of those two flags. Deletions inside with model: are rolled back on exit, so the effect of a deletion can be measured without keeping it.

Deleting a gene#

Deleting a gene differs from deleting its reactions. Removing a gene takes the gene and its occurrences in the reaction rules out of the model, and the reactions stay with rewritten GPRs. With the default policy for reactions whose last catalyst is removed, the flux state equals that of a knockout (bounds (0, 0)); the model differs in that the gene is removed.

reduced = removeGenes(model, {'YBR196C'});
fprintf('%d -> %d genes\n', numel(model.genes), numel(reduced.genes));
disp(reduced.grRules{getIndexes(reduced, 'PGI', 'rxns')});

61 -> 60 genes
YLR354C

from raven_toolbox.manipulation import remove_genes

with model:
    remove_genes(model, ["YBR196C"])
    print(len(model.genes), "genes")
    print("PGI rule:", repr(model.reactions.get_by_id("PGI").gene_reaction_rule))

60 genes
PGI rule: 'YLR354C'

removeGenes and remove_genes rewrite each affected GPR without the removed gene. A reaction whose last catalyst is removed stays in the model with its rule dropped and its bounds set to (0, 0): this is the default of remove_genes (blocked_reactions="constrain"), while removeGenes does it when removeBlockedRxns is false, its default, and it is the bound state that knock_out_model_genes leaves. With blocked_reactions="keep" the blocked reactions keep their bounds; setting removeBlockedRxns to true, or blocked_reactions="remove", deletes them instead, which turns a gene deletion into a reaction deletion; Deletions and essentiality covers when that is the question being asked. When at least one requested gene is in the model, genes that no remaining reaction used are removed as well, also when they were not requested; when none is in the model, no gene is removed. In removeGenes, standardizeRules, on by default, normalises the rewritten rules.

Troubleshooting

  • Editing model.S by hand. It leaves model.rev, the bounds and the gene matrix describing a different model than the matrix does. Use changeRxns.

  • Deleting a reaction and leaving its metabolites without a reaction. Without the orphan flags, metabolites that nothing produces or consumes remain in the model and appear as gaps. See Quality control.

  • Overwriting a GPR that was meant to be extended. changeGrRules and change_gene_reaction_rules replace by default.

  • Python edits in place. MATLAB copies the struct on assignment, so modelB = removeReactions(modelA, ...) leaves modelA unchanged. In Python, use model.copy() or with model:.

See also#