• Main
• About
• FBA
• Prerequisites
• Guide
• Group

---

In this tutorial we make basic operations with Genome Scale Metabolic (GSM) models, chemical reactions and the underlying network topology

---

• The programming examples are based on Python3 programming language.

from numpy import *

Import the COBRApy (Constraint-Based Reconstruction and Analysis in Python) package to load and explore SBML (Systems Biology Markup Language) models.

• For instuction to install the package read: https://pypi.org/project/cobra/
• For documentation read: https://cobrapy.readthedocs.io/en/latest/

import cobra

Import the NetworkX and igraph packages for statistics on the underlying network topologies:

• For instuctions to install the packages read: https://networkx.org/documentation/stable/install.html and https://igraph.org/python/
• For documentations read: https://networkx.org/documentation/stable/index.html and https://igraph.org/python/doc/igraph.Graph-class.html

import networkx
import igraph

Import RDKit and its necessary packages to visualize chemical elements and reactions.

• For instalation read: https://depth-first.com/articles/2020/08/17/getting-started-rdkit-and-jupyter/ (Note here that this can be done only by creating a virtual environment by Anaconda https://docs.anaconda.com/anaconda/install/linux/)
• For documentation read: https://www.rdkit.org/docs/GettingStartedInPython.html and https://iwatobipen.wordpress.com/2020/01/10/edit-reaction-image-and-draw-it-rdkit-memo/

from rdkit import Chem
from rdkit.Chem import rdChemReactions as Reactions
from rdkit.Chem.Draw import IPythonConsole

Load the model downloaded from the BiGG database (http://bigg.ucsd.edu/models) and check some basic properties. More information and useful tutorials can be found on https://github.com/webermarcolivier/metabolic_modelling_jupyter_tutorial

model = cobra.io.load_json_model('iJO1366.json')

print('Nb of reactions: ',len(model.reactions))
print('Nb of metabolites: ',len(model.metabolites))
print('Nb of genes: ',len(model.genes))

Nb of reactions:  2583
Nb of metabolites:  1805
Nb of genes:  1367

Here we choose a custom reaction, check some of its properties and visualize it by RDKit

rexpl = 603
reaction = model.reactions[rexpl]

for r in reaction.reactants:
    print("NAME =", r.name, "| FORMULA =", r.formula)

NAME = L-Alanine | FORMULA = C3H7NO2
NAME = Pyridoxal 5'-phosphate | FORMULA = C8H8NO6P

for p in reaction.products:
    print("NAME =", p.name, "| FORMULA =", p.formula)

NAME = Pyridoxamine 5'-phosphate | FORMULA = C8H12N2O5P
NAME = Pyruvate | FORMULA = C3H3O3

print("Name of reaction:",reaction.name)

Name of reaction: Alanine transaminase

Converting reaction names to SMILES (Symplified Molecular Input Line Entry System, https://archive.epa.gov/med/med_archive_03/web/html/smiles.html, https://www.daylight.com/dayhtml_tutorials/languages/smiles/index.html) used by RDKit. One can use for example this online automatic converter: https://opsin.ch.cam.ac.uk/

name2smiles = {
               "L-Alanine"                 : "N[C@@H](C)C(=O)O",
               "Pyridoxal 5'-phosphate"    : "P(=O)(O)(O)OCC=1C(=C(C(=NC1)C)O)C=O",
               "Pyridoxamine 5'-phosphate" : "P(=O)(O)(O)OCC=1C(=C(C(=NC1)C)O)CN",
               "Pyruvate"                  : "C(C(=O)C)(=O)[O-]"
              }

Visualize a metabolite:

mol = Chem.MolFromSmiles(name2smiles["Pyridoxamine 5'-phosphate"])
mol

Visualize a reaction:

rxn = Reactions.ReactionFromSmarts(name2smiles["L-Alanine"]+"."+name2smiles["Pyridoxal 5'-phosphate"]+">>"+
                                   name2smiles["Pyridoxamine 5'-phosphate"]+"."+
                                   name2smiles["Pyruvate"], useSmiles=True)
rxn

Convert SBML model to NetworkX graph object:

• enzymatic reaction were ommited (fully catalized system)
• directions were considered
• $(A + B \rightarrow C + D) \equiv (A \rightarrow C; A \rightarrow D; B \rightarrow C; B \rightarrow D)$

def SBML2NetworkX(model):
    graph = networkx.DiGraph()
    for i in range(len(model.reactions)):
        Ri = model.reactions[i]
        if Ri.reversibility == False:
            Rir = Ri.reactants
            Rip = Ri.products
            for j in range(len(Rir)):
                for k in range(len(Rip)):
                    if (Rir[j].id[-1:] != 'e') and (Rip[k].id[-1:] != 'e'):
                        graph.add_edge(Rir[j].id, Rip[k].id)
        if Ri.reversibility == True:
            Rir = Ri.reactants
            Rip = Ri.products
            for j in range(len(Rir)):
                for k in range(len(Rip)):
                    if (Rir[j].id[-1:] != 'e') and (Rip[k].id[-1:] != 'e'): 
                        graph.add_edge(Rir[j].id, Rip[k].id)
                        graph.add_edge(Rip[k].id, Rir[j].id)
    return graph

gnx = SBML2NetworkX(model)

Extract the edge-list and save it for further use:

gnx = networkx.convert_node_labels_to_integers(gnx)
EL = array(gnx.edges).astype("int32")
EL.tofile("expl_metabolic.bin")

Read back edge-list and construct igraph Graph object:

EL = fromfile("expl_metabolic.bin", dtype = int32).reshape((-1, 2)) # read in binary file
EL = EL[where((EL[:, 0] == EL[:, 1]) == 0)[0]] # delete self-loops
UEL = unique(EL)
EL = vectorize(dict(zip(UEL, arange(len(UEL)))).__getitem__)(EL) # make indexing after 0 for igraph

gig = igraph.Graph()
gig.add_vertices(UEL)
gig.add_edges(EL)

print("Number of nodes: ", gig.vcount()) # number of vertices
print("Number of edges: ", gig.ecount()) # number of edges

Number of nodes:  1481
Number of edges:  7694

Visualization of the graph:

• For a better interpretability and to avoid overcrowding of the figure we use only a small part of the network

pEL = EL[:1000]
pUEL = unique(pEL)
pEL = vectorize(dict(zip(pUEL, arange(len(pUEL)))).__getitem__)(pEL) # make indexing after 0 for igraph
gig = igraph.Graph()
gig.add_vertices(pUEL)
gig.add_edges(pEL)

out_fig_name = "graph.eps"

visual_style = {}
# Set vertex size
visual_style["vertex_size"] = 5
# Don't curve the edges
visual_style["edge_curved"] = False
# Set the layout
my_layout = gig.layout_fruchterman_reingold()
visual_style["layout"] = my_layout

# Plot the graph
igraph.plot(gig, out_fig_name, **visual_style)

For a more detailed graph investigations check: https://github.com/MateJozsaPhys/CNDinvestigation/tree/main/measure-codes