python script, extracting stoichiometries form SBML file, printing as LP

print_stoichiometries.py

# -*- coding: utf-8 -*-
"""
Basic example program reading an SBML file going through the reactions 
and printing the stoichiometries as constraints in LP format. 

Created on Thu Apr 14 19:12:41 2016

@author: fbergmann
"""

import sys

from  libsbml import *

def printStoichiometry(stoichiometries):
  """ This functions prints the stoichiometries found as LP constraints"""
  print '/* constraints */'
  for item in stoichiometries:
      print '{0} = 0;'.format(stoichiometries[item])

def addSpeciesReference(stoichiometries, reactionId, reference, operation):
  """Adds the given species reference to the list of stoichiometries
  
     Parameters:   
     
        stoichiometries: dictionary of speciesId with the stoichiometry
        reactionId:      the current reaction id to add
        reference:       the species reference to add
        operation:       string, '-' in case of reactant, '+' in case of product
  
  """
  species = reference.getSpecies()
  stoichiometry = reference.getStoichiometry()
  current = ''
  if species in stoichiometries: 
    current = stoichiometries[species]
  if stoichiometry == 1: 
    stoichiometries[species] = '{0} {1}{2}'.format(current,operation,reactionId)
  else: 
    stoichiometries[species] = '{0} {1}{2} {3}'.format(current, operation, stoichiometry, reactionId)

def printStoichiometryFor(file):
  """Reads the SBML file, extracts the stoichiometry, and prints it as constraint"""
  doc = readSBMLFromFile(file)
  model = doc.getModel()
  stoichiometries = {}

  for i in range(0, model.getNumReactions()): 
      reaction = model.getReaction(i)
      reactionId = reaction.getId();
      
      for j in range (0, reaction.getNumReactants()):
          reactant = reaction.getReactant(j);
          addSpeciesReference(stoichiometries, reactionId, reactant, '-')

      for j in range (0, reaction.getNumProducts()):
          product= reaction.getProduct(j);
          addSpeciesReference(stoichiometries, reactionId, product, '+')
  
  printStoichiometry(stoichiometries)
  pass

if __name__ == '__main__':
    if len(sys.argv) > 1: 
      printStoichiometryFor(sys.argv[1])
    else: 
      printStoichiometryFor('BorisEJB.xml')