Package SloppyCell :: Package ReactionNetworks :: Module SBMLInterface
[hide private]

Source Code for Module SloppyCell.ReactionNetworks.SBMLInterface

  1  """ 
  2  Methods for loading from and saving to SBML files. 
  3  """ 
  4  __docformat__ = "restructuredtext en" 
  5   
  6  import os 
  7  import sys 
  8  import logging 
  9  logger = logging.getLogger('RxnNets.SBMLInterface') 
 10   
 11  import libsbml 
 12   
 13  import Network_mod 
 14  import SloppyCell.ExprManip as ExprManip 
 15  import SloppyCell.KeyedList_mod 
 16  KeyedList = SloppyCell.KeyedList_mod.KeyedList 
 17   
 18  import SloppyCell.ExprManip as ExprManip 
 19   
20 -def toSBMLFile(net, fileName):
21 sbmlStr = toSBMLString(net) 22 f = file(fileName, 'w') 23 f.write(sbmlStr) 24 f.close()
25
26 -def SBMLtoDOT(sbmlFileName, dotFileName):
27 raise DeprecationWarning, 'SBMLtoDOT has been deprecated. Instead, use IO.net_DOT_file(net, filename)'
28
29 -def rxn_add_stoich(srxn, rid, stoich, is_product=True):
30 sr = libsbml.SpeciesReference(rid) 31 try: 32 stoich = float(stoich) 33 if stoich < 0: 34 sr.setStoichiometry(-stoich) 35 srxn.addReactant(sr) 36 if stoich > 0: 37 sr.setStoichiometry(stoich) 38 srxn.addProduct(sr) 39 if stoich == 0: 40 sr = libsbml.ModifierSpeciesReference(rid) 41 srxn.addModifier(sr) 42 except ValueError: 43 formula = stoich.replace('**', '^') 44 math_ast = libsbml.parseFormula(formula) 45 smath = libsbml.StoichiometryMath(math_ast) 46 sr.setStoichiometryMath(smath) 47 if is_product == True: 48 srxn.addProduct(sr) 49 else: 50 srxn.addReactant(sr)
51 52 53 54
55 -def toSBMLString(net):
56 m = libsbml.Model(net.id) 57 m.setName(net.name) 58 59 for id, fd in net.functionDefinitions.items(): 60 sfd = libsbml.FunctionDefinition(id) 61 sfd.setName(fd.name) 62 formula = fd.math 63 formula = formula.replace('**', '^') 64 formula = 'lambda(%s, %s)' % (','.join(fd.variables), formula) 65 sfd.setMath(libsbml.parseFormula(formula)) 66 m.addFunctionDefinition(sfd) 67 68 for id, c in net.compartments.items(): 69 sc = libsbml.Compartment(id) 70 sc.setName(c.name) 71 sc.setConstant(c.is_constant) 72 sc.setSize(c.initialValue) 73 m.addCompartment(sc) 74 75 for id, s in net.species.items(): 76 ss = libsbml.Species(id) 77 ss.setName(s.name) 78 ss.setCompartment(s.compartment) 79 if s.initialValue is not None and not isinstance(s.initialValue, str): 80 ss.setInitialConcentration(s.initialValue) 81 ss.setBoundaryCondition(s.is_boundary_condition) 82 m.addSpecies(ss) 83 84 for id, p in net.parameters.items(): 85 sp = libsbml.Parameter(id) 86 sp.setName(p.name) 87 if p.initialValue is not None: 88 sp.setValue(p.initialValue) 89 sp.setConstant(p.is_constant) 90 m.addParameter(sp) 91 92 for id, r in net.rateRules.items(): 93 sr = libsbml.RateRule() 94 sr.setVariable(id) 95 formula = r.replace('**', '^') 96 sr.setMath(libsbml.parseFormula(formula)) 97 m.addRule(sr) 98 99 for id, r in net.assignmentRules.items(): 100 sr = libsbml.AssignmentRule() 101 sr.setVariable(id) 102 formula = r.replace('**', '^') 103 sr.setMath(libsbml.parseFormula(formula)) 104 m.addRule(sr) 105 106 for r, r in net.algebraicRules.items(): 107 sr = libsbml.AlgebraicRule() 108 formula = r.replace('**', '^') 109 sr.setMath(libsbml.parseFormula(formula)) 110 m.addRule(sr) 111 112 for id, rxn in net.reactions.items(): 113 srxn = libsbml.Reaction(id) 114 srxn.setName(rxn.name) 115 # Handle the case where the model was originally read in from an 116 # SBML file, so that the reactants and products of the Reaction 117 # object are explicitly set. 118 if rxn.reactant_stoichiometry != None and \ 119 rxn.product_stoichiometry != None: 120 for rid, stoich_list in rxn.reactant_stoichiometry.items(): 121 for stoich in stoich_list: 122 rxn_add_stoich(rxn, rid, stoich, is_product=False) 123 for rid, stoich_list in rxn.product_stoichiometry.items(): 124 for stoich in stoich_list: 125 rxn_add_stoich(rxn, rid, stoich, is_product=True) 126 # Handle the case where the model was created using the SloppyCell 127 # API, in which case reactants and products are inferred from their 128 # stoichiometries 129 else: 130 for rid, stoich in rxn.stoichiometry.items(): 131 rxn_add_stoich(srxn, rid, stoich) 132 formula = rxn.kineticLaw.replace('**', '^') 133 kl = libsbml.KineticLaw(formula) 134 srxn.setKineticLaw(kl) 135 m.addReaction(srxn) 136 137 for id, e in net.events.items(): 138 se = libsbml.Event(id) 139 se.setName(e.name) 140 formula = e.trigger.replace('**', '^') 141 ast = libsbml.parseFormula(formula) 142 if ast is None: 143 raise ValueError('Problem parsing event trigger: %s. Problem may ' 144 'be use of relational operators (< and >) rather ' 145 'than libsbml-friendly functions lt and gt.' 146 % formula) 147 try: 148 se.setTrigger(ast) 149 except TypeError: 150 se.setTrigger(libsbml.Trigger(ast)) 151 152 formula = str(e.delay).replace('**', '^') 153 try: 154 se.setDelay(libsbml.parseFormula(formula)) 155 except TypeError: 156 se.setDelay(libsbml.Delay(libsbml.parseFormula(formula))) 157 for varId, formula in e.event_assignments.items(): 158 sea = libsbml.EventAssignment() 159 sea.setVariable(varId) 160 formula = str(formula).replace('**', '^') 161 ast = libsbml.parseFormula(formula) 162 sea.setMath(ast) 163 se.addEventAssignment(sea) 164 m.addEvent(se) 165 166 for id, con in net.constraints.items(): 167 scon = libsbml.Constraint() 168 scon.setId(con.id) 169 scon.setName(con.name) 170 formula = con.trigger.replace('**', '^') 171 ast = libsbml.parseFormula(formula) 172 if ast is None: 173 raise ValueError('Problem parsing constraint math: %s. Problem may ' 174 'be use of relational operators (< and >) rather ' 175 'than libsbml-friendly functions lt and gt.' 176 % formula) 177 #try: 178 scon.setMath(ast) 179 #except TypeError: 180 # scon.setMath(libsbml.Trigger(ast)) 181 182 m.addConstraint(scon) 183 184 185 d = libsbml.SBMLDocument() 186 d.setModel(m) 187 sbmlStr = libsbml.writeSBMLToString(d) 188 189 return sbmlStr
190
191 -def _print_sbml_fatals(doc):
192 for ii in range(doc.getNumFatals()): 193 logger.critical(d.getFatal(ii).getMessage())
194
195 -def to_SBML_l2v1(from_name, to_name):
196 """ 197 Convert an SBML file to level 2, version 1 using lisbml. 198 199 from_name Name of file to read SBML from. 200 to_name File to output SBML to. 201 """ 202 doc = libsbml.readSBML(os.path.abspath(from_name)) 203 if isinstance(doc, int): 204 logger.critical('Fatal Errors reading SBML from file %s!' % from_name) 205 _print_sbml_fatals(doc) 206 errors = doc.setLevel(2) 207 if errors: 208 logger.critical('Fatal Errors converting %f to level 2!' % from_name) 209 _print_sbml_fatals(doc) 210 errors = doc.setVersion(1) 211 if errors: 212 logger.critical('Fatal Errors converting %f to level 2, version 1!' 213 % from_name) 214 _print_sbml_fatals(doc) 215 success = libsbml.writeSBML(doc, to_name) 216 if not success: 217 logger.critical('Error writing to %s' % to_name)
218
219 -def fromSBMLFile(fileName, id = None, duplicate_rxn_params=False):
220 f = file(fileName, 'r') 221 net = fromSBMLString(f.read(), id, duplicate_rxn_params) 222 f.close() 223 return net
224
225 -def stoichToString(species, stoich):
226 if stoich is None: 227 stoich = str(species.getStoichiometry()) 228 elif hasattr(stoich, 'getMath'): # libsbml > 3.0 229 stoich = libsbml.formulaToString(stoich.getMath()) 230 else: # libsbml 2.3.4 231 stoich = libsbml.formulaToString(stoich) 232 return stoich
233
234 -def fromSBMLString(sbmlStr, id = None, duplicate_rxn_params=False):
235 r = libsbml.SBMLReader() 236 d = r.readSBMLFromString(sbmlStr) 237 if d.getNumErrors(): 238 message = 'libSBML reported errors in SBML file. Try running file '\ 239 'through the online validator: '\ 240 'http://www.sbml.org/Facilities/Validator . Specific errors '\ 241 'noted are: ' 242 errors = [] 243 for ii in range(d.getNumErrors()): 244 pm = d.getError(ii) 245 errors.append(pm.getMessage()) 246 raise ValueError(message + '; '.join(errors)) 247 248 m = d.getModel() 249 250 modelId = m.getId() 251 if (id == None) and (modelId == ''): 252 raise ValueError('Network id not specified in SBML or passed in.') 253 elif id is not None: 254 modelId = id 255 256 rn = Network_mod.Network(id = modelId, name = m.getName()) 257 258 for f in m.getListOfFunctionDefinitions(): 259 id, name = f.getId(), f.getName() 260 math = f.getMath() 261 variables = [] 262 for ii in range(math.getNumChildren() - 1): 263 variables.append(libsbml.formulaToString(math.getChild(ii))) 264 265 math = libsbml.formulaToString(math.getRightChild()) 266 267 rn.addFunctionDefinition(id, variables, math) 268 269 for c in m.getListOfCompartments(): 270 id, name = c.getId(), c.getName() 271 size = c.getSize() 272 isConstant = c.getConstant() 273 274 rn.addCompartment(id = id, size = size, 275 isConstant = isConstant, 276 name = name) 277 278 for s in m.getListOfSpecies(): 279 id, name = s.getId(), s.getName() 280 compartment = s.getCompartment() 281 if s.isSetInitialConcentration(): 282 iC = s.getInitialConcentration() 283 elif s.isSetInitialAmount(): 284 iC = s.getInitialAmount() 285 else: 286 iC = 1 287 isBC, isConstant = s.getBoundaryCondition(), s.getConstant() 288 289 xml_text = s.toSBML() 290 uniprot_ids = set([entry[1:].split('"')[0] 291 for entry in xml_text.split('uniprot')[1:]]) 292 293 rn.addSpecies(id = id, compartment = compartment, 294 initialConcentration = iC, 295 isConstant = isConstant, 296 is_boundary_condition = isBC, 297 name = name, uniprot_ids = uniprot_ids) 298 299 for p in m.getListOfParameters(): 300 parameter = createNetworkParameter(p) 301 rn.addVariable(parameter) 302 303 for rxn in m.getListOfReactions(): 304 id, name = rxn.getId(), rxn.getName() 305 kL = rxn.getKineticLaw() 306 kLFormula = kL.getFormula() 307 308 substitution_dict = {} 309 # Deal with parameters defined within reactions 310 for p in kL.getListOfParameters(): 311 parameter = createNetworkParameter(p) 312 # If a parameter with this name already exists, **and it has a 313 # different value than this parameter** we rename this parameter 314 # instance by prefixing it with the rxn name so there isn't a 315 # clash. 316 if parameter.id in rn.variables.keys(): 317 logger.warn('Parameter %s appears in two different reactions ' 318 'in SBML file.' % parameter.id) 319 if parameter.value != rn.variables.get(parameter.id).value or\ 320 duplicate_rxn_params: 321 oldId = parameter.id 322 parameter.id = id + '_' + parameter.id 323 substitution_dict[oldId] = parameter.id 324 logger.warn('It has different values in the two positions ' 325 'so we are creating a new parameter %s.' 326 % (parameter.id)) 327 else: 328 logger.warn('It has the same value in the two positions ' 329 'so we are only defining one parameter %s. ' 330 'This behavior can be changed with the option ' 331 'duplicate_rxn_params = True' % (parameter.id)) 332 333 if parameter.id not in rn.variables.keys(): 334 rn.addVariable(parameter) 335 kLFormula = ExprManip.sub_for_vars(kLFormula, substitution_dict) 336 337 # Assemble the stoichiometry. SBML has the annoying trait that 338 # species can appear as both products and reactants and 'cancel out' 339 # For each species appearing in the reaction, we build up a string 340 # representing the stoichiometry. Then we'll simplify that string and 341 # see whether we ended up with a float value in the end. 342 stoichiometry = {} 343 reactant_stoichiometry = {} 344 product_stoichiometry = {} 345 for reactant in rxn.getListOfReactants(): 346 species = reactant.getSpecies() 347 stoichiometry.setdefault(species, '0') 348 stoich = reactant.getStoichiometryMath() 349 stoich = stoichToString(reactant, stoich) 350 stoichiometry[species] += '-(%s)' % stoich 351 if species in reactant_stoichiometry: 352 reactant_stoichiometry[species].append(stoich) 353 else: 354 reactant_stoichiometry[species] = [stoich] 355 356 for product in rxn.getListOfProducts(): 357 species = product.getSpecies() 358 stoichiometry.setdefault(species, '0') 359 stoich = product.getStoichiometryMath() 360 stoich = stoichToString(product, stoich) 361 stoichiometry[species] += '+(%s)' % stoich 362 if species in product_stoichiometry: 363 product_stoichiometry[species].append(stoich) 364 else: 365 product_stoichiometry[species] = [stoich] 366 367 for species, stoich in stoichiometry.items(): 368 stoich = ExprManip.simplify_expr(stoich) 369 try: 370 # Try converting the string to a float. 371 stoich = float(stoich) 372 except ValueError: 373 pass 374 stoichiometry[species] = stoich 375 376 for modifier in rxn.getListOfModifiers(): 377 stoichiometry.setdefault(modifier.getSpecies(), 0) 378 379 rn.addReaction(id = id, stoichiometry = stoichiometry, 380 kineticLaw = kLFormula, 381 reactant_stoichiometry = reactant_stoichiometry, 382 product_stoichiometry = product_stoichiometry) 383 384 for ii, r in enumerate(m.getListOfRules()): 385 if r.getTypeCode() == libsbml.SBML_ALGEBRAIC_RULE: 386 math = libsbml.formulaToString(r.getMath()) 387 rn.add_algebraic_rule(math) 388 else: 389 variable = r.getVariable() 390 math = libsbml.formulaToString(r.getMath()) 391 if r.getTypeCode() == libsbml.SBML_ASSIGNMENT_RULE: 392 rn.addAssignmentRule(variable, math) 393 elif r.getTypeCode() == libsbml.SBML_RATE_RULE: 394 rn.addRateRule(variable, math) 395 396 397 398 for ii, e in enumerate(m.getListOfEvents()): 399 id, name = e.getId(), e.getName() 400 if id == '': 401 id = 'event%i' % ii 402 403 try: 404 # For libSBML 3.0 405 trigger_math = e.getTrigger().getMath() 406 except AttributeError: 407 # For older versions 408 trigger_math = e.getTrigger() 409 trigger = libsbml.formulaToString(trigger_math) 410 411 if e.getDelay() is not None: 412 try: 413 # For libSBML 3.0 414 delay_math = e.getDelay().getMath() 415 except AttributeError: 416 # For older versions 417 delay_math = e.getDelay() 418 delay = libsbml.formulaToString(delay_math) 419 else: 420 delay = 0 421 422 timeUnits = e.getTimeUnits() 423 eaDict = KeyedList() 424 for ea in e.getListOfEventAssignments(): 425 eaDict.set(ea.getVariable(), libsbml.formulaToString(ea.getMath())) 426 427 rn.addEvent(id = id, trigger = trigger, eventAssignments = eaDict, 428 delay = delay, name = name) 429 430 431 for ii, con in enumerate(m.getListOfConstraints()): 432 id, name = con.getId(), con.getName() 433 if id == '': 434 id = 'constraint%i' % ii 435 436 trigger_math = con.getMath() 437 438 trigger = libsbml.formulaToString(trigger_math) 439 440 if con.isSetMessage(): 441 message = con.getMessage() 442 else: 443 message = None 444 445 rn.addConstraint(id = id, trigger = trigger, message = message, 446 name = name) 447 448 return rn
449
450 -def createNetworkParameter(p):
451 id, name = p.getId(), p.getName() 452 v = p.getValue() 453 isConstant = p.getConstant() 454 455 parameter = Network_mod.Parameter(id = id, value = v, is_constant = isConstant, 456 name = name, typical_value = None, is_optimizable = True) 457 # optimizable by default 458 459 return parameter
460