Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Parse SBML stoichiometry matrix
from __future__ import print_function
import libsbml
import argparse
def _parser():
parser = argparse.ArgumentParser(description="Parse stoichiometry matrix of SBML file")
parser.add_argument('file', metavar="filename", type=argparse.FileType('r'),
help="Filename of SBML file to parse")
return parser
def _main():
parser = _parser()
args = parser.parse_args()
file_ = args.file
species, reactions, stoichiometry_matrix = parse_file(file_)
_print_stoichiometry_matrix(species, reactions, stoichiometry_matrix)
def parse_file(open_file_):
reader = libsbml.SBMLReader()
document = reader.readSBML(open_file_.name)
sbml_model = document.getModel()
species = [s.getName() for s in sbml_model.getListOfSpecies()]
reactions = [r.getId() for r in sbml_model.getListOfReactions()]
stoichiometry_matrix = {}
for reaction_index, reaction in enumerate(sbml_model.getListOfReactions()):
reactants = {r.getSpecies(): r.getStoichiometry() for r in reaction.getListOfReactants()}
products = {p.getSpecies(): p.getStoichiometry() for p in reaction.getListOfProducts()}
for species_index, species_node in enumerate(sbml_model.getListOfSpecies()):
species_id = species_node.getId()
net_stoichiometry = int(reactants.get(species_id, 0)) - int(products.get(species_id, 0))
stoichiometry_matrix[species_index, reaction_index] = net_stoichiometry
return species, reactions, stoichiometry_matrix
def _print_stoichiometry_matrix(species, reactions, stoichiometry_matrix):
print('\t'.join(['---'] + reactions))
for species_ix, species_label in enumerate(species):
to_print = [species_label]
for reaction_ix in range(len(reactions)):
stoichiometry = stoichiometry_matrix[species_ix, reaction_ix]
to_print.append(stoichiometry)
print('\t'.join(map(str, to_print)))
if __name__ == '__main__':
_main()
@qacwnfq

This comment has been minimized.

Copy link

qacwnfq commented Sep 25, 2019

Concering following statement
net_stoichiometry = int(reactants.get(species_id, 0)) - int(products.get(species_id, 0))

This only works for systems where there are no averaged biomass equations, because these averaged are usually not integers.

@lukauskas

This comment has been minimized.

Copy link
Owner Author

lukauskas commented Sep 27, 2019

@qacwnfq It's been a while since I worked on this so not sure how to fix it.
Would you suggest making it into a float instead?

@qacwnfq

This comment has been minimized.

Copy link

qacwnfq commented Sep 28, 2019

@lukauskas
Yes, that could work :)
The main reason behind my comment was just documentation for any one else using this gist, because it worked very well on some sbml files and not on others for me.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.