Getting started#

Load a model, see how large it is, and inspect a single reaction, metabolite and gene. Every later page assumes these operations.

Functions#

MATLAB

Python

Purpose

readYAMLmodel

read_yaml_model

read a RAVEN YAML model

importModel

read_sbml_model cobrapy

read an SBML model

printModelStats

model attributes cobrapy

how big the model is

getIndexes

get_by_id cobrapy

look something up by identifier

constructEquations

Reaction.reaction cobrapy

a reaction as a readable string

getElementalBalance

check_mass_balance cobrapy

is a reaction balanced

Setup#

The examples use smallYeast.yml, a model of central carbon metabolism in yeast that ships with RAVEN. Download it from docs/data/smallYeast.yml and run everything from the directory that contains it.

Install one of the toolboxes first (Installation). Every page lists the files it needs under its Setup heading; they are all in docs/data/, or clone the repository and work in raven-docs/docs/data. The blocks on a page build on each other: run them in order in one MATLAB or Python session.

Loading the model#

model = readYAMLmodel('smallYeast.yml');
disp(model.id);

smallYeast

A RAVEN model is a MATLAB struct of parallel arrays. model.rxns, model.mets and model.genes hold the identifiers, and every other field lines up with one of them by position: model.lb(4) is the lower bound of the reaction named in model.rxns{4}. The stoichiometric matrix model.S has one row per metabolite and one column per reaction, in those same orders. Working with the model means working with indices into these arrays, so looking an identifier up comes first.

from raven_toolbox.io import read_yaml_model

model = read_yaml_model("smallYeast.yml")
print(model.id)

smallYeast

In Python a RAVEN model is a cobra.Model. raven-toolbox defines no model class of its own, so every cobrapy tool, and anything else in the Python COBRA ecosystem, operates on it directly. Reactions, metabolites and genes are objects held in collections. Each object records what it is connected to, so there are no indices to keep aligned.

For SBML use importModel in MATLAB and cobrapy's read_sbml_model in Python; Reading and writing models covers every supported format.

Model size#

printModelStats(model);

Network statistics for smallYeast: Central carbon metabolism for yeast
Genes*				61
	cytosol	52
	mitochondria	17

Reactions*			53
	cytosol	45
	mitochondria	19
Unique reactions**	53

Metabolites			52
	cytosol	35
	mitochondria	17
Unique metabolites	45

* Genes and reactions are counted for each compartment if any of the corresponding metabolites are in that compartment. The sum may therefore not add up to the total number.
** Unique reactions are defined as being biochemically unique (no compartmentalization)

The per-compartment counts do not sum to the total, because a reaction or gene is counted in every compartment in which one of its metabolites is located; a transport reaction appears under both. Unique means biochemically unique, ignoring compartments: 52 metabolites reduce to 45 because seven of them exist in both the cytosol and the mitochondrion. The same 45 is what Combining and simplifying models arrives at when it flattens the compartments.

print(len(model.reactions), len(model.metabolites), len(model.genes))

53 52 61

There is no equivalent of printModelStats. The collections are attributes, so their lengths are the totals, and per-compartment counts come from filtering them. After a solve, model.summary() prints the exchange fluxes and the objective.

Inspecting a reaction#

Glucose-6-phosphate isomerase, PGI, shows most of what a reaction carries: it is reversible, it has a gene association, and it should be mass balanced.

idx = getIndexes(model, 'PGI', 'rxns');
disp(model.rxnNames{idx});
eqn = constructEquations(model, model.rxns(idx));
fprintf('%s\n', eqn{1});
fprintf('bounds: [%g %g]\n', model.lb(idx), model.ub(idx));
disp(model.grRules{idx});

Glucose-6-phosphate isomerase
alpha-D-glucose 6-phosphate[c] <=> beta-D-fructofuranose 6-phosphate[c]
bounds: [-1000 1000]
YBR196C

getIndexes turns an identifier into the position used by every other field. constructEquations builds the readable equation from the stoichiometric column, substituting metabolite names for identifiers, so the equation reads in words.

pgi = model.reactions.get_by_id("PGI")
print(pgi.name)
print(pgi.reaction)
print("bounds:", pgi.bounds, "reversible:", pgi.reversibility)
print(pgi.gene_reaction_rule)

Glucose-6-phosphate isomerase
G6P_c <=> F6P_c
bounds: (-1000.0, 1000.0) reversible: True
YBR196C

Reaction.reaction substitutes metabolite identifiers, so the same reaction reads differently from the MATLAB equation above. Both describe the same stoichiometry.

The bounds carry the direction. A lower bound of −1000 means the reaction may run backwards; a lower bound of zero means it may not. cobrapy derives reversibility from the bounds each time it is asked, so the two cannot disagree. RAVEN stores an explicit model.rev field alongside the bounds, so one can be set without the other; Model structure and identifiers covers what to do about that.

Inspecting a metabolite#

idx = getIndexes(model, 'G6P_c', 'mets');
fprintf('%s (%s) in compartment %s\n', model.metNames{idx}, ...
    model.metFormulas{idx}, model.comps{model.metComps(idx)});
fprintf('takes part in %d reactions\n', nnz(model.S(idx, :)));

alpha-D-glucose 6-phosphate (C6H13O9P) in compartment c
takes part in 4 reactions

The reaction count comes from counting non-zero entries in the metabolite's row of model.S. model.metComps holds an index into model.comps, not the compartment letter itself, so the compartment needs two lookups.

g6p = model.metabolites.get_by_id("G6P_c")
print(f"{g6p.name} ({g6p.formula}) in compartment {g6p.compartment}")
print(f"takes part in {len(g6p.reactions)} reactions")

alpha-D-glucose 6-phosphate (C6H13O9P) in compartment c
takes part in 4 reactions

Metabolite.reactions is maintained by the model as reactions are added and removed, so it needs no matrix lookup.

The formula makes a mass-balance check possible. A metabolite without one leaves every reaction it appears in undecidable, which is the distinction drawn in Mass balance of a reaction.

Inspecting a gene#

gi = find(strcmp(model.genes, 'YBR196C'));
rxnIds = model.rxns(model.rxnGeneMat(:, gi) ~= 0);
fprintf('%s -> %s\n', model.geneShortNames{gi}, strjoin(rxnIds, ', '));

PGI1 -> PGI

model.rxnGeneMat is a reactions-by-genes incidence matrix; a non-zero entry means that gene appears in that reaction's rule. It records which genes are involved; how they combine is recorded only in model.grRules, as and and or in text.

gene = model.genes.get_by_id("YBR196C")
print(gene.name, "->", [rxn.id for rxn in gene.reactions])

PGI1 -> ['PGI']

Both toolboxes distinguish the identifier from the name. YBR196C is the systematic identifier and PGI1 the standard gene name; the model is keyed on the former.

Mass balance of a reaction#

Draft models routinely contain reactions that do not balance. Checking one is the same operation in both toolboxes; doing it for a whole model is covered in Quality control.

balance = getElementalBalance(model, 'rxns', {'PGI'});
fprintf('elemental %d, charge %d\n', balance.balanceStatus, balance.chargeStatus);

elemental 1, charge -1

The two statuses are reported separately, and each has three values: 1 balanced, 0 unbalanced, -1 undecidable. PGI balances elementally and its charge is undecidable: the metabolites in this model carry formulas but no charges, so there is nothing to sum. That is a different statement from "the charges do not balance": one says the check could not run, the other says it failed.

print(pgi.check_mass_balance())

{}

An empty dictionary means the reaction balances. Anything listed is an element, or charge, mapped to the amount by which the two sides differ, so the result says both whether and by how much. cobrapy reports charge in the same dictionary, not separately, and omits what it cannot decide, so an empty result covers both "balanced" and "nothing to check".

Troubleshooting

  • A lookup fails. Identifiers are case-sensitive and carry the compartment suffix (G6P_c, not G6P). Neither toolbox returns a sentinel for a name it cannot find: getIndexes errors with Could not find object 'X' in the model and get_by_id raises KeyError, so a typo stops the script; nothing else is indexed without warning.

  • The model loads but nothing grows. In smallYeast.yml every uptake reaction is closed (glcIN has bounds [0 0]). Opening a medium is covered in Growth media and conditions.

  • Gene identifiers differ between model and FASTA. Systematic names (YBR196C) and standard names (PGI1) are not interchangeable; RAVEN stores the systematic name as the identifier and the standard name as the gene name.

See also#