modelmage.generator

31 """ 32 This class is the core of model generation in modelMaGe. It creates a graph of 33 from the model and provides methods to generate candidate models from this. 34 """ 35 DG = NW.XDiGraph() 36 # Dicts to store all the graphs and models produced by removing different nodes 37 __changedDG = {} 38 # Dicts to store different attributes of the model 39 kinetics = {} 40 __stoich = {} 41 __graphHistory = {} # stores which are the parent reactions of a combined reaction 42 possibleKinetics = {} 43 44

45 - def __init__(self, model, functionsIdName, verbose=None):

46 """ 47 Init method of Generator class. 48 49 @type model: Model 50 @param model: the master model that will be used to generate new models 51 52 @type verbose: str 53 @param verbose: switch that defines if the method returns vrbose output or not. 54 """ 55 56 self.__model = model 57 self.__stoich = model.getStoich() 58 self.verbose = verbose 59 # self.lawFormularModifier = {} # stores the law function for copasi file remove 60 self.functionsIdName = functionsIdName 61 62 # content defined in def __lawtoSBML() 63 # used in def assignFunctionNames() 64 self.nameKineticNameMap = {} # stores the changed names of kinetics 65 66 self.idNameReactions = {} # stores the id-name map of reactions 67 self.definedKinetic = {} # stores the previous defined kinetics 68 # this pat used in __updateReaction, importTask, assignFunctionNames 69 self.pat = re.compile('reaction_re') 70 self.REVERSEMACROS = '' #store aliases for modelnames

71

72 - def getGraphHistory(self, nodes=None):

73 """ 74 Returns the history of the nodes in the generated models. The key of a 75 newly generated reaction is the combined key of the reactions it is combined of. 76 77 @rtype: list 78 @return: history of nodes 79 """ 80 if nodes: 81 return self.__graphHistory[nodes] 82 else: 83 return self.__graphHistory

84

85 - def setGraphHistory(self, removedNodes, history):

86 """ 87 Set the history for a model. The model is defined by the nodes which are 88 removed. 89 90 @type removedNodes: str 91 @param removedNodes: Removed nodes in the model. This serves as a key in the history dict. 92 93 @type history: dict 94 @param history : the dictionary containing the history of the newly added nodes 95 """ 96 self.__graphHistory[removedNodes] = history

97

98 - def getStoich(self, pair):

99 if self.__stoich.has_key(pair): 100 return self.__stoich[pair] 101 else: 102 return False

103

104 - def setStoich(self, pair, stoich):

105 if self.__stoich.has_key(pair): 106 self.__stoich[pair] += stoich 107 else: 108 self.__stoich[pair] = stoich

109

110 - def getModel(self):

111 return self.__model

112

113 - def __checkSelfLoop(self, removed, predecessor, successor):

114 """ 115 This method checks for resulting self-loops in case a species is removed. 116 Self-loops are removed from the graph/model. 117 """ 118 suc_products = [x[1] for x in removed.out_edges(successor) if x[2] == 'product'] 119 suc_reactants = [x[0] for x in removed.in_edges(successor) if x[2] == 'reactant'] 120 for val in [x for x in suc_products if x in suc_reactants]: 121 suc_products.remove(val) 122 123 pred_products = [x[1] for x in removed.out_edges(predecessor) if x[2] == 'product'] 124 pred_reactants = [x[0] for x in removed.in_edges(predecessor) if x[2] == 'reactant'] 125 for val in [x for x in pred_products if x in pred_reactants]: 126 pred_reactants.remove(val) 127 128 pred_reactants.sort() 129 suc_products.sort() 130 131 if pred_reactants == suc_products: 132 # self.__updateSelfLoops(products) 133 if self.verbose: 134 print removed.predecessors(successor), "=", removed.successors(predecessor), ": don't connect", successor, 'and', predecessor 135 return 1 136 else: 137 return 0

138

139 - def getChangedDG(self):

140 return self.__changedDG

141 142

143 - def printReactionsSpecies(self, remove=None):

144 """ 145 Prints lists of reactions and species of an SBML model to stdout. 146 147 @param remove: the model to get reactions and species from 148 @type remove: libsbml.model 149 150 @return: none 151 @rtype: none 152 """ 153 154 removeElements = [] 155 if remove != None: 156 for ree in remove: 157 removeElements.extend(ree) 158 species = self.__model.getModelElement("Species", ["Id", "Name"]) 159 reactions = self.__model.getModelElement("Reactions", ["Id", "Name"]) 160 reacModifiersMap = {} 161 162 graph = self.__model.getGraph() 163 164 print '\n\nModel species:' 165 for s in species: 166 if (species.index(s) > 0) & (species.index(s) % 3 == 0) : 167 if (s[0] in removeElements) | (s[1] in removeElements): 168 print "%s(remove), " % s 169 else: 170 print "%s, " % s 171 elif species.index(s) == len(species) - 1: 172 if (s[0] in removeElements) | (s[1] in removeElements): 173 print "%s(remove)" % s 174 else: 175 print s 176 else: 177 if (s[0] in removeElements) | (s[1] in removeElements): 178 print "['%s','%s'](remove), " % (s[0], s[1]), 179 else: 180 print "['%s','%s'], " % (s[0], s[1]), 181 182 print "\nModel reactions:" 183 for reaction in reactions: 184 reactants = [r[0] for r in graph.in_edges(reaction) if r[2] == 'reactant'] 185 modifiers = [r[0] for r in graph.in_edges(reaction) if r[2] == 'modifier'] 186 products = [r[1] for r in graph.out_edges(reaction) if r[2] == 'product'] 187 188 try: 189 reactants.sort() 190 modifiers.sort() 191 products.sort() 192 except: 193 pass 194 195 if (reaction[0] in removeElements) | (reaction[1] in removeElements): 196 print "%s(remove)" % reaction 197 else: 198 print reaction 199 200 if reactants: 201 print "\tsubstrates: %s"%reactants 202 if products: 203 print "\tproducts : %s"%products 204 if modifiers: 205 print "\tmodifiers : %s"%modifiers

206 207

208 - def removeNode(self, nodes):

209 """ 210 Removes the nodes v from the graph DG. In case the node is a species, it combines the reactions 211 which are connected to v in an intelligent way. This always leads to connected graphs 212 if species are removed, but might lead to unconnected graphs in case reactions are removed. 213 214 @type nodes: [string] 215 @param nodes: list of nodes to remove from the model 216 217 @return: removed 218 @rtype: graph 219 """ 220 DG = self.__model.getGraph() 221 removed = DG.copy() 222 223 #dummyhistory for each structure 224 nodeHistory = {} 225 self.setGraphHistory(str(nodes), nodeHistory) 226 227 # order nodes to be removed according to reaction, modifiers,species 228 species = self.__model.getModelElement("Species", ["Id", "Name"]) 229 species_ids = [x[0] for x in species] 230 reactions = self.__model.getModelElement("Reactions", ["Id", "Name"]) 231 reaction_ids = [x[0] for x in reactions] 232 233 # first remove reactions 234 for node in nodes: 235 if node in reaction_ids: 236 self.__removeReaction(removed, node, nodes) 237 # them remove modifiers 238 for node in nodes: 239 if node.find(':') != -1: 240 item = node.split(':') 241 reaction = item[0] 242 modifier = item[1] 243 self.__removeModifier(removed, reaction, modifier) 244 # finally remove species 245 for node in nodes: 246 if node in species_ids: 247 nodeHistory = self.__removeSpecies(removed, node, nodeHistory) 248 249 #check for reactions that occur twice but with different names and delete them: 250 reactions = [] 251 for node in removed.nodes(): 252 if [x for x in removed.in_edges(node) if x[2] == 'reactant'] or [x for x in removed.out_edges(node) if x[2] == 'product']: 253 reactions.append(node) 254 for node in reactions: 255 involvedSpecies = ([x[0] for x in removed.in_edges(node) if x[2] == 'reactant'], \ 256 [x[0] for x in removed.in_edges(node) if x[2] == 'modifier'], \ 257 [x[1] for x in removed.out_edges(node)]) 258 for otherNode in [x for x in reactions if x != node]: 259 otherInvolvedSpecies = ([x[0] for x in removed.in_edges(otherNode) if x[2] == 'reactant'], \ 260 [x[0] for x in removed.in_edges(otherNode) if x[2] == 'modifier'], \ 261 [x[1] for x in removed.out_edges(otherNode)]) 262 if involvedSpecies == otherInvolvedSpecies: 263 try: 264 removed.delete_node(otherNode) 265 except: 266 pass 267 if self.verbose: 268 print 'Deleting \"%s\" beause it is the same as \"%s\".' % (otherNode, node) 269 270 return removed

271

272 - def __removeModifier(self, removed, reaction, modifier):

273 """ 274 Removes a modifier from a reaction. 275 276 @param removed: graph that will be the result of removed node 277 @type removed: networkx.DG 278 279 @param reaction: the id of the reaction to remove the modifier from 280 @type reaction: str 281 282 @param modifier: id of the species to be removed as a modifier from reaction 283 @type modifier: str 284 285 """ 286 287 print 'Warning: %s acts as a modifier in %s and will be removed.' %(modifier,reaction) 288 # modifiers is a true modifier 289 if removed.has_edge(modifier,reaction,'modifier'): 290 if self.verbose == 1: 291 print "Removing %s as modifier in %s..." % (modifier , reaction) 292 removed.delete_edge(modifier, reaction,'modifier') 293 # modifier is substrate and product 294 elif removed.has_edge(modifier,reaction,'reactant') and removed.has_edge(reaction,modifier,'product'): 295 if self.verbose == 1: 296 print "Removing %s as modifier in %s..." % (modifier , reaction) 297 removed.delete_edge(modifier,reaction,'reactant') 298 removed.delete_edge(reaction,modifier,'product') 299 else: 300 print "%s is not a modifier of reaction %s. Skipping ..."%(modifier,reaction)

301 302

303 - def __removeReaction(self, removed , reaction, nodes):

304 """ 305 Removes a reaction from a given graph. 306 307 @param removed: graph that will be the result of remove node 308 @type removed: networkx.DG 309 310 @param reaction: the id of the reaction to remove 311 @type reaction: str 312 313 @return: none 314 @rtype: none 315 """ 316 try: 317 if self.verbose: 318 print 'removing reaction %s...' % reaction 319 removed.delete_node(reaction) 320 if len(NW.connected_components(removed.to_undirected())) > 1 and self.verbose: 321 print "Warning: removing %s creates unconnected subnetworks!" % reaction 322 323 except NW.exception.NetworkXError, e: 324 try: 325 newNodes = [x[0] for x in self.__graphHistory[str(nodes)].items() if reaction in x[1]] 326 for node in newNodes: 327 removed.delete_node(node) 328 if len(NW.connected_components(removed.to_undirected())) > 1 and self.verbose: 329 print "Warning: removing %s creates unconnected subnetworks!" % node 330 except: 331 raise generatorError(NW.exception.NetworkXError, e)

332

333 - def __removeSpecies(self, removed, species, nodeHistory):

334 """ 335 Removes a species from a graph. 336 337 @param removed: graph that will be the result of remove node 338 @type removed: networkx.DG 339 340 @param species: the id of the species to be removed 341 @type species: str 342 343 @param nodeHistory: dictionary that stores the history of newly created reactions 344 @type nodeHistory: dict 345 346 @return: nodeHistory 347 @rtype: dict 348 349 """ 350 #remove one node after another from the graph "removed" 351 #nodeHistory = {} #dict to save the reactions that are combined to a new reaction 352 # list for reactions that can be removed afterwards 353 remove = [] 354 if self.verbose: 355 print 'removing species %s...' % species 356 try: 357 predecessors = removed.predecessors(species) 358 successors = [x[1] for x in removed.out_edges(species) if x[2] == "reactant"] 359 modified = [x[1] for x in removed.out_edges(species) if x[2] == "modifier"] 360 except NW.exception.NetworkXError, e: 361 raise generatorError(NW.exception.NetworkXError, e) 362 363 # if species does neither have predecessors nor successors it can be deleted, because 364 # there is no connection between other nodes 365 if not predecessors or not successors: 366 remove.append(species) 367 368 # first remove the species from all reactions, 369 # where it is a modifier. 370 for i in modified: 371 self.__removeModifier(removed, i, species) 372 373 # and for reactions where it is substrate as well as product 374 r_remove=[] 375 for i in predecessors: 376 for j in successors: 377 if i==j: 378 self.__removeModifier(removed, i, species) 379 r_remove.append(i) 380 381 for i in r_remove: 382 predecessors.remove(i) 383 successors.remove(i) 384 385 # now check for self-loops, otherwise combine appropriate reactions 386 for i in predecessors: 387 for j in successors: 388 #detect and leave out selfloops 389 if self.__checkSelfLoop(removed, i, j): 390 remove.extend([i, j]) 391 continue 392 #combine reactions 393 else: 394 newNode = 'reaction_%s_%s' % (i, j) 395 self.__combine(removed, i, j, species, newNode) 396 #save the parent nodes of the combined node in dict 397 nodeHistory[newNode] = [] 398 399 if nodeHistory.has_key(i): 400 nodeHistory[newNode].extend(nodeHistory[i]) #copy previous history 401 if nodeHistory.has_key(j): 402 nodeHistory[newNode].extend(nodeHistory[j]) #copy previous history 403 nodeHistory[newNode].extend([i, j]) 404 #self.history[str(species)] = nodeHistory 405 remove.extend([i, j]) 406 407 # Mark reactions which have no substrates or products left for removal 408 for i in predecessors: 409 if not [x for x in removed.in_edges(i) if x[2] == 'reactant'] \ 410 and not [x for x in removed.out_edges(i) if x[2] == 'product']: 411 remove.extend(i) 412 if self.verbose: 413 print 'Removing Reaction %s' % i 414 for j in successors: 415 if removed.degree(j) <= 1: 416 remove.extend(j) 417 if self.verbose: 418 print 'Removing Reaction %s' % j 419 420 try: 421 removed.delete_node(species) 422 except: 423 pass 424 #remove reactions that are no longer needed 425 for node in remove: 426 try: 427 removed.delete_node(node) 428 except: 429 pass 430 431 return nodeHistory

432

433 - def __combine(self, DG, i, j, v, newNode):

434 """ 435 Combines two reactions if they have more than one reactant or product. 436 Therefore it sets the products and reactants of reaction 'j' directly to 437 products and reactants of reaction 'i' and then renames 'i' to 'reaction_i_j'. 438 Species that are modifiers to reaction i as well as to reaction j are kept only once. 439 440 @type DG: networkx.XDiGraph 441 @param DG: the graph in which the reactions will be combined(should be the graph where v is removed) 442 443 @type i: string 444 @param i: first reaction for combine 445 446 @type j: string 447 @param j: second reaction for combine 448 449 @type v: string 450 @param v: removed species 451 """ 452 #combine i and j 453 products = [x[1] for x in DG.out_edges(i) if x[2] == 'product'] 454 reactants = [x[0] for x in DG.in_edges(i) if x[2] == 'reactant'] 455 for edge in DG.in_edges(i): 456 if edge[0] != v and edge[0] not in products: 457 DG.add_edge(edge[0] , newNode, edge[2]) 458 try: 459 self.__stoich[(edge[0], newNode)] = self.__stoich[(edge[0], i)] 460 except: 461 self.__stoich[(edge[0], newNode)] = 1 462 # check for modifiers that are declared as substrate as well as product 463 # and re-define them as true modifiers 464 elif edge[0] != v and edge[0] in products: 465 DG.add_edge(edge[0] , newNode, 'modifier') 466 for edge in DG.out_edges(i): 467 if edge[1] != v and edge[1] not in reactants: # modifiers have already been detected 468 DG.add_edge(newNode, edge[1], edge[2]) 469 try: 470 self.__stoich[(newNode, edge[1])] = self.__stoich[(i, edge[1])] 471 except: 472 self.__stoich[(newNode, edge[1])] = 1 473 474 products = [x[1] for x in DG.out_edges(j) if x[2] == 'product'] 475 reactants = [x[0] for x in DG.in_edges(j) if x[2] == 'reactant'] 476 for edge in DG.in_edges(j): 477 if edge[0] != v and edge[0] not in products: 478 DG.add_edge(edge[0] , newNode, edge[2]) 479 try: 480 self.__stoich[(edge[0], newNode)] = self.__stoich[(edge[0], j)] 481 except: 482 self.__stoich[(edge[0], newNode)] = 1 483 # check for modifiers that are declared as substrate as well as product 484 # and re-define them as true modifiers 485 elif edge[0] != v and edge[0] in products: 486 DG.add_edge(edge[0] , newNode, 'modifier') 487 for edge in DG.out_edges(j): 488 if edge[1] != v and edge[1] not in reactants: # modifiers have already been detected 489 DG.add_edge(newNode, edge[1], edge[2]) 490 try: 491 self.__stoich[(newNode, edge[1])] = self.__stoich[(j, edge[1])] 492 except: 493 self.__stoich[(newNode, edge[1])] = 1 494 #adjust stoichiometric values for non-modifiers 495 modifiers = [x[0] for x in DG.in_edges(newNode) if x[2] == 'modifier'] 496 equalProducts = [x for x in DG.successors(j) if x != v and x in DG.successors(i) and x not in modifiers]#intersection of lists 497 equalReactants = [x for x in DG.predecessors(j) if x != v and x in DG.predecessors(i) and x not in modifiers] 498 for p in equalProducts: 499 self.__stoich[(newNode, p)] += 1 500 501 for p in equalReactants: 502 self.__stoich[(p, newNode)] += 1 503 504 #remove duplicates in modifiers: 505 tmp = [] 506 for i in modifiers: 507 if i in tmp: 508 DG.delete_edge(i,newNode,'modifier') 509 tmp.append(i) 510 if self.verbose: 511 print "combined %s and %s" % (i, j)

512 513

514 - def generate(self, remove=[], kinetics=False, REVERSEMACROS=None, verbose=False):

515 ''' 516 Returns model objects that are generated according to give kinetics and remove 517 commands. Further the method returns dictionaries containing Maps of funcitons 518 to reactions and ids to names. 519 ''' 520 #print '\n------------------------' 521 #print "Generating structural alternatives:" 522 #print '----------------------\n' 523 524 # save the functionsname from copasi 525 self.functionsName = [] 526 remove.sort() 527 for nodes in remove: 528 try: 529 if verbose: 530 print "\nCreating structure %s:" % nodes 531 self.__changedDG[str(nodes)] = self.removeNode(nodes) 532 except generatorError, e: 533 print e.message 534 print 'Species or reaction is not part of the model please choose from the lists below:\n' 535 self.printReactionsSpecies() 536 exit(1) 537 if not remove: 538 self.setGraphHistory('', {}) 539 self.__changedDG[''] = self.__model.getGraph() 540 541 self.REVERSEMACROS = REVERSEMACROS 542 if kinetics: 543 self.kinetics.update(kinetics) 544 newModels = self.__createKineticSBML() 545 546 return (newModels, self.nameKineticNameMap, self.modelFunctionsIdName, self.idNameReactions)

547 548

549 - def __findStructureKinetics(self):

550 """ 551 Determines the possible kinetics of all generated Structures and returns the kinetics in a dictionary. 552 553 @rtype: {str:{str:[str]}} 554 @return: a dictionary containing the possible kinetics for each reaction in every structure. 555 {structure:{reaction:[kinetics]}} 556 """ 557 558 masterGraph = self.__model.getGraph() 559 kineticsForStructure = {} 560 species = self.__model.getModelElement("Species", "id") 561 reactions = self.__model.getModelElement("Reactions", ["Id", "Name"]) 562 563 # define map with key reaction-id and reaction-name 564 # key: reaction_1 value: v1 565 for r in reactions: 566 rId = r[0] 567 rName = r[1] 568 self.idNameReactions[rId] = rName 569 570 ReverseMacros = self.REVERSEMACROS[0] 571 ReverseKineticMacros = self.REVERSEMACROS[1] 572 #for all changed structures create a dictionary containing the kinetic alternatives that apply to them 573 for id, graph in self.__changedDG.items(): 574 assignModel = '' 575 if id in ReverseMacros.keys(): 576 assignModel = ReverseMacros[id] 577 578 idKin = {} #alternative kinetics for this structure 579 for key in ReverseKineticMacros.keys(): 580 revReact, revKin, revMa = key.split('.') 581 if revMa == assignModel: 582 if revReact in idKin.keys(): 583 idKin[revReact].append(revKin) 584 else: 585 idKin[revReact] = [revKin] 586 587 combineTable = {} 588 589 #iterate over all nodes in the graph that are reactions 590 for reaction in [x for x in graph.nodes() if x not in species]: 591 possibleKinetics = [] 592 593 #number of modifiers and reactants of the current reaction 594 reactants = [x[0] for x in graph.in_edges(reaction) if x[2] == 'reactant'] 595 modifiers = [x[0] for x in graph.in_edges(reaction) if x[2] == 'modifier'] 596 products = [x[1] for x in graph.out_edges(reaction) if x[2] == 'product'] 597 598 for val in [x for x in products if x in reactants]: 599 modifiers.append(val) 600 products.remove(val) 601 reactants.remove(val) 602 try: 603 reactants.sort() 604 modifiers.sort() 605 products.sort() 606 except: 607 pass 608 609 # self.graphHistory key: ['species_4'] value: {'reaction_reaction_3_reaction_4': ['reaction_3', 'reaction_4']} 610 # if reaction is not combination of reactions but might be modified (modifier removed), 611 # check whether an alternative kinetic is given or orignal kinetic still applies. 612 if not self.getGraphHistory(id).has_key(reaction): #did the reaction exist in the master model? 613 # Check if the original kinetic law still fits the reaction. If it does, insert it 614 # first check whether there is an alternative kinetic specified 615 # self.kinetics key: reaction_10 value: ['mMA'] 616 if reaction in idKin.keys(): 617 possibleKinetics.extend(idKin[reaction]) 618 else: 619 masterReactants = [x[0] for x in masterGraph.in_edges(reaction) if x[2] == 'reactant'] 620 masterModifiers = [x[0] for x in masterGraph.in_edges(reaction) if x[2] == 'modifier'] 621 masterProducts = [x[1] for x in masterGraph.out_edges(reaction) if x[2] == 'product'] 622 623 masterProducts.sort() 624 masterModifiers.sort() 625 masterReactants.sort() 626 #check for substrates that are products as well and count as modifiers 627 for val in [x for x in masterProducts if x in masterReactants]: 628 masterModifiers.append(val) 629 masterProducts.remove(val) 630 masterReactants.remove(val) 631 # 632 try: 633 masterReactants.sort() 634 masterModifiers.sort() 635 masterProducts.sort() 636 except: 637 pass 638 639 if masterModifiers == modifiers and \ 640 masterReactants == reactants and \ 641 masterProducts == products: 642 possibleKinetics.append(self.__model.reactionId[reaction].getKineticLaw()) 643 644 # if the the reaction is a combination of two or more reactions check if there are 645 # alternative kinetics for the parents 646 elif self.getGraphHistory(id).has_key(reaction): 647 combinationOf = self.getGraphHistory(id)[reaction] 648 for parentReaction in [x for x in combinationOf if x in idKin.keys()]: 649 possibleKinetics.extend(idKin[parentReaction]) 650 651 # mass action kinetics can always be inserted if there is nothing else 652 if possibleKinetics == []: 653 654 if len(reactants) == 0: 655 defaultKinetic = "CF" 656 else: 657 defaultKinetic = "MA" 658 for m in modifiers: 659 defaultKinetic = "m" + defaultKinetic 660 possibleKinetics.append(defaultKinetic) 661 662 #remove duplicates: 663 possibleKinetics.sort() 664 new = [] 665 for x in possibleKinetics: 666 if x not in new: 667 new.append(x) 668 possibleKinetics = new 669 670 # print 'reaction: %s, kin: %s' % (reaction, possibleKinetics) 671 combineTable[reaction] = possibleKinetics 672 673 #the possible kinetics for all structures are returned 674 kineticsForStructure[id] = combineTable 675 676 return kineticsForStructure

677

678 - def __checkConsistentReverseMacros(self, ReverseMacros, ReverseKineticMacros):

679 """ 680 For case ca1 = species_2 | species_5, there will be ca1.cps, ca1K1.cps, ca1K2.cps 681 the keys in 2 ReverseMacros should be correspondingly changed 682 ReverseMacros: {"['species_2', 'species_5']": 'ca1', "['species_5']": 'ca1', "['species_2']": 'ca1'} 683 ReverseKineticMacros: {'reaction_6:MM:ca1': 'ca1', 'reaction_5:mMM:ca1': 'ca1'} 684 """ 685 usedName = {} 686 addReverseKineticitems = {} 687 for key, value in ReverseMacros.items(): 688 if value in usedName.keys(): 689 del ReverseMacros[key] 690 changed = value + 'K' + str(usedName[value]) 691 ReverseMacros[key] = changed 692 va = usedName[value] 693 694 for kineticKey in ReverseKineticMacros.keys(): 695 if value in ReverseKineticMacros.values(): 696 newKineticKey = re.sub(value, changed, kineticKey) 697 #ReverseKineticMacros[newKineticKey] = changed 698 addReverseKineticitems[newKineticKey] = changed 699 usedName[value] = va + 1 700 else: 701 usedName[value] = 1 702 703 if len(addReverseKineticitems) > 0: 704 for addKey, addValue in addReverseKineticitems.items(): 705 ReverseKineticMacros[addKey] = addValue 706 707 return (ReverseMacros, ReverseKineticMacros)

708 709

710 - def __createKineticSBML(self):

711 """ 712 Creates SBML models from structure graphs and kinetics. 713 714 @rtype: list 715 @return: a list of generated models that are using different kinetics 716 """ 717 718 kineticsForStructure = self.__findStructureKinetics() 719 structureModels = self.__createStructureSBML() 720 721 newModels = {} 722 ReverseMacros = self.REVERSEMACROS[0] 723 ReverseKineticMacros = self.REVERSEMACROS[1] 724 ReverseMacros, ReverseKineticMacros = self.__checkConsistentReverseMacros(ReverseMacros, ReverseKineticMacros) 725 #print "ReverseMacros: %s\tReverseKineticMacros: %s"%(ReverseMacros, ReverseKineticMacros) 726 # stores function id name for each model 727 self.modelFunctionsIdName = {} 728 729 # initialise the index 730 indexFuDe = len(self.functionsIdName.items()) 731 732 # use id as model alias name, in case when ReverseMarcros not available 733 # it means there is no user defined model alias like in initial file 734 useIdasName = False 735 if (len(structureModels) > 1) & (len(ReverseMacros) == 0): 736 useIdasName = True 737 738 for id in structureModels.keys(): 739 # number of laws per reaction: 740 numberOfLaws = [] 741 kinetics = {} 742 numberOfModels = 1 743 self.indexFuDe = indexFuDe + 1 744 # get model 745 sbmlModel = structureModels[id].getModel() 746 747 # get graph 748 sbmlGraph = self.__changedDG[id] 749 750 # Get the total number of models that will be created for the structure 751 listOfReactions = kineticsForStructure[id].keys() 752 753 # make sure if all reations is taking place in the same compartment 754 self.inTheSameCompartment = self.checkCompartment(sbmlModel) 755 756 # get the assigned model name for current structure 757 assignModel = '' 758 if id in ReverseMacros.keys(): 759 assignModel = ReverseMacros[id] 760 # default value used, doesn't mean anything 761 self.modelFunctionsIdName[assignModel] = [assignModel] 762 763 # want to match file like M1Rereaction_7K0.cps M1Rereaction_7reaction_1K0.cps 764 if useIdasName: 765 idTrueNames = re.findall('\'(.*?)\'', id) 766 for idT in idTrueNames: 767 assignModel = assignModel + idT 768 self.modelFunctionsIdName[assignModel] = [assignModel] 769 770 # initialize a new functionsIdName, because of each sbmlModel 771 # needs his own. 772 fIdName = {} 773 # copy the standard functionsIdName for master model 774 for ttid, ttname in self.functionsIdName.items(): 775 fIdName[ttid] = ttname 776 777 for reaction in listOfReactions: 778 779 # change id to name for reaction 780 if reaction not in self.idNameReactions.keys(): 781 self.idNameReactions[reaction] = reaction 782 783 # get possible kinetics 784 listOfLaws = kineticsForStructure[id][reaction] 785 numberOfModels *= len(listOfLaws) 786 numberOfLaws.append(len(listOfLaws)) 787 SBMLListOfLaws = [] 788 789 for law in listOfLaws: 790 if (type(law) == str): 791 SBMLlaw, functionDefinition = self.__lawToSBML(law, reaction, sbmlGraph, assignModel, fIdName) 792 SBMLListOfLaws.append(SBMLlaw) 793 if functionDefinition.getMath() != None: 794 sbmlModel.addFunctionDefinition(functionDefinition) 795 else: 796 SBMLListOfLaws.append(law) 797 kinetics[reaction] = SBMLListOfLaws 798 799 self.modelFunctionsIdName[assignModel] = fIdName 800 801 oldDocument = lSBML.SBMLDocument() # only needed for transcription to xmlobject 802 oldDocument.setModel(sbmlModel) 803 804 for modelNumber in range(numberOfModels): 805 newDocument = lSBML.SBMLDocument() 806 newSbmlModel = sbmlModel.clone() 807 newSbmlModel.setSBMLDocument(newDocument) 808 newDocument.setModel(newSbmlModel) 809 newModel = model.Model(newDocument) 810 811 #Copy Annotation from old model and add the information which nodes have been removed 812 #use xmlobject to manipulate the xml inside the annotation tag 813 #writer = lSBML.SBMLWriter() 814 raw = '<?xml version="1.0" encoding="UTF-8"?>\n%s' % (oldDocument.toSBML()) 815 rootObject = XMLFile(raw=raw) 816 try: 817 annotationObject = rootObject.sbml.model.annotation 818 except Exception, e: 819 rootObject.sbml.model._addNode('annotation') 820 annotationObject = rootObject.sbml.model.annotation 821 raise generatorError(e, 'Model annotation in rootfile is missing or not valid.') 822 try: 823 annotationObject.modelMaGe 824 except: 825 annotationObject._addNode('modelMaGe') 826 annotationObject.modelMaGe.xmlns = "http://www.modelmage.org/" 827 annotationObject = rootObject.sbml.model.annotation 828 number = annotationObject.modelMaGe._addNode('kineticNumber') 829 number.id = str(modelNumber) 830 831 newSbmlModel.setAnnotation(str(annotationObject._toxml())) 832 833 #compute which kinetic law has to be used 834 factorA = float(modelNumber) 835 for reaction in listOfReactions: 836 factorB = 1. 837 reactionIndex = listOfReactions.index(reaction) 838 for number in numberOfLaws[reactionIndex + 1:]: 839 factorB = factorB * number 840 if reactionIndex == len(listOfReactions) - 1: 841 factorB = 1. 842 kineticLawIndex = int(math.floor((factorA / factorB) % numberOfLaws[reactionIndex])) 843 844 reactionInModel = newSbmlModel.getReaction(reaction) 845 reactionInModel.unsetKineticLaw() 846 #print "reaction: %s\tkinetic: %s"%(reaction, kinetics[reaction][kineticLawIndex].getFormula()) 847 reactionInModel.setKineticLaw(kinetics[reaction][kineticLawIndex]) 848 newSbmlModel.setName(sbmlModel.getName()) 849 newDocument.setModel(newSbmlModel) 850 if id in newModels.keys(): 851 newModels[id].append(newModel) 852 else: 853 newModels[id] = [newModel] 854 return newModels

855 856

857 - def checkCompartment(self, sbmlModel):

858 ''' 859 check if there is more than one compartment of current system 860 861 @type sbmlModel: Model 862 @param sbmlModel: a model object 863 864 @rtype: boolean 865 @return: 1 if one compartment, 0 if more compartments 866 ''' 867 compartments = sbmlModel.getListOfCompartments() 868 #print "compartment: %s\n"%sbmlModel.getCompartment(0).getId() 869 if len(compartments) == 1: 870 self.thisCompartment = compartments[0].getId() 871 return 1 872 else: 873 return 0

874 875

876 - def __lawToSBML(self, law, reaction, graph, assignModel, fIdName):

877 """ 878 Constructs an SBML KineticLaw object out of an XML Node 879 880 @type law: minidom:Dom Element 881 @param law: A kinetic law element from the input XML-file which will be translated to SBML 882 883 @type reaction: str 884 @param reaction: The reaction this law applies to 885 886 @type graph: networkx.XDiGraph 887 @param graph: the structure graph the reaction belongs to 888 889 @type assignModel: str 890 @param assignModel: the name of current candidate model 891 892 @type fIdName: list 893 @param fIdName: {function_id:function_name} 894 895 @return: Kinetic law that can be inserted into a reaction and functionDefinition 896 @rtype: libSBML:KineticLaw 897 898 """ 899 #get the species which belong to the reaction 900 products = [x[1] for x in graph.out_edges(reaction)] 901 #take the first element from the in_edges which are reactants XOR modifiers 902 reactants = [x[0] for x in graph.in_edges(reaction) if x[2] == 'reactant'] 903 modifiers = [x[0] for x in graph.in_edges(reaction) if x[2] == 'modifier'] 904 905 # check wehther there are modifiers that are declared as substrate as well as products 906 for val in [x for x in products if x in reactants]: 907 modifiers.append(val) 908 products.remove(val) 909 reactants.remove(val) 910 911 tmp = [] 912 #print "lawtoSBML: %s"%reaction 913 914 # get the function id-Name Map 915 #withIIHM = 0 # if the assigned kinetic has specious modification 916 mathFormula = '' 917 formula = '' 918 name = '' 919 suffix = '' 920 count = 0 # count how many modifiers a reaction has 921 amountModi = [] # initial value of assigned kinetic 922 923 functionDefinition = lSBML.FunctionDefinition() 924 functionDefinition.setId('function_' + str(self.indexFuDe)) 925 functionVariables = '' 926 # create a new kinetic law for xml 927 kineticLaw = lSBML.KineticLaw() 928 929 # if law is defined by user, 930 if (type(law) == str): 931 amountModi = re.findall('ih|i|m', law) 932 if len(amountModi) != len(modifiers): 933 print "The assigned kinetic %s with '%s' modifiers is not equal with present %s modifiers in %s"%(law, len(amountModi), len(modifiers), reaction) 934 exit(1) 935 936 #constant flux irrevesible 937 # kinetic: compartment_1 * function_6(v) 938 #if len(reactants) == 0 : 939 if re.findall('CF$', law): 940 if len(reactants) > 0: 941 print "\nThe assigned kinetic %s for %s is not appropriate,"%(law, reaction) 942 print "%s has one or more reactants. Exiting ..."%reaction 943 exit(1) 944 suffix = 'Constant flux (irreversible)' 945 mathFormula = 'k' 946 formula = '(k' + ',' 947 functionVariables = 'k' + ',' 948 #add parameters to the reaction 949 lawParameterV = lSBML.Parameter() 950 lawParameterV.setId('k') 951 lawParameterV.setValue(0.1) 952 kineticLaw.addParameter(lawParameterV) 953 954 #michaelis menten kinetics 955 elif re.findall('MM$', law): 956 if len(reactants) != 1: 957 print "\nThe assigned kinetic %s for %s is not appropriate,"%(law, reaction) 958 print "%s has not exactly one reactant. Exiting ..."%reaction 959 exit(1) 960 suffix = 'Michaelis-Menten' 961 mathFormula = 'vmax*S/(Km+S)' 962 formula = '(vmax,' + reactants[0] + ',Km' + ',' 963 functionVariables = 'vmax,' + 'S,' + 'Km' + ',' 964 #add parameters to the reaction 965 lawParameterV = lSBML.Parameter() 966 lawParameterV.setId('vmax') 967 lawParameterV.setValue(0.1) 968 kineticLaw.addParameter(lawParameterV) 969 970 lawParameterKM = lSBML.Parameter() 971 lawParameterKM.setId('Km') 972 lawParameterKM.setValue(0.1) 973 kineticLaw.addParameter(lawParameterKM) 974 975 #mass action kinetics: 976 elif re.findall('MA$', law): 977 suffix = 'Mass action' 978 #string to insert mass action kinetics into the reaction 979 mathFormula = 'k'+ '*' 980 formula = '(k' + ',' 981 functionVariables = functionVariables + 'k' + ',' 982 reacIndex = 1 983 for reactant in reactants: 984 for i in range(int(self.__stoich[(reactant, reaction)])): 985 mathFormula = mathFormula + 'S' + str(reacIndex) + '*' 986 formula = formula + reactant + ',' 987 functionVariables = functionVariables + 'S' + str(reacIndex) + ',' 988 reacIndex += 1 989 mathFormula = mathFormula.strip('*') 990 #add parameters to the reaction 991 lawParameter = lSBML.Parameter() 992 lawParameter.setId('k') 993 lawParameter.setValue(0.1) 994 kineticLaw.addParameter(lawParameter) 995 996 # currently only the mass action and his modified kinetics and michaelis menten and his 997 # modified kinetics are allowed, otherwise will appear the following info, then determine 998 else: 999 print "The kinetic %s you have specified for the %s is not allowed" % (law, reaction) 1000 sys.exit(1) 1001 1002 # if law is defined by libsbml, then it is not a string 1003 # the following check would be superfluence 1004 if (type(law) == str): 1005 #if re.findall('ih|i|m', law): 1006 if len(amountModi) > 0: 1007 # filter out duplicates: 1008 m = [] 1009 for x in modifiers: 1010 if x not in m: 1011 m.append(x) 1012 modifiers = m 1013 addParameterH = False 1014 preFactors = re.findall('ih|i|m', law) 1015 for num, value in enumerate(preFactors): 1016 mIndex = str(count + 1) 1017 if value == 'i': 1018 name = 'inhibited' 1019 try: 1020 #mathFormula = mathFormula + '*1/(1+' + modifiers[num] + ')' 1021 formula = formula + modifiers[num] + ',' 1022 mathFormula = mathFormula + '*1/(1+M' + str(mIndex) + ')' 1023 except Exception, e: 1024 print '\nError: \n%s\nFor %s there are less modifiers than user defined' % (e, reaction) 1025 exit(1) 1026 1027 if value == 'ih': 1028 name = 'inhibited' 1029 try: 1030 #mathFormula = mathFormula + '*1/(1+' + modifiers[num] + '**h)' 1031 mathFormula = mathFormula + '*1/(1+M' + str(mIndex) + '^h' + str(mIndex) + ')' 1032 formula = formula + modifiers[num] + ',h' + str(mIndex) + ',' 1033 addParameterH = True 1034 except Exception, e: 1035 print '\nError: \n%s\nFor %s there are less modifiers than user defined' % (e, reaction) 1036 exit(1) 1037 1038 if value == 'm': 1039 name = 'modified' 1040 try: 1041 #mathFormula = mathFormula + '*' + modifiers[num] 1042 mathFormula = mathFormula + '*M' + str(mIndex) 1043 formula = formula + modifiers[num] + ',' 1044 except Exception, e: 1045 print '\nError: \n%s\nFor %s there are less modifiers than user defined' % (e, reaction) 1046 exit(1) 1047 count += 1 1048 functionVariables = functionVariables + 'M' + str(count) + ',' 1049 if addParameterH: 1050 lawParameterH = lSBML.Parameter() 1051 lawParameterH.setId('h' + str(mIndex)) 1052 lawParameterH.setValue(2) 1053 kineticLaw.addParameter(lawParameterH) 1054 functionVariables = functionVariables + 'h' + str(mIndex) + ',' 1055 addParameterH = False 1056 1057 # doesn't match any patter for kinetic definition 1058 if suffix == '': 1059 return (law, functionDefinition) 1060 1061 else: 1062 formula = formula.strip(',') + ')' 1063 1064 if count == 0: 1065 value = name + ' ' + suffix 1066 else: 1067 if name != '': 1068 value = name + '(' + str(count) + ') ' + suffix 1069 else: 1070 value = 'modified(' + str(count) + ') ' + suffix 1071 1072 function_index = 'existIndex' 1073 passFunction = False 1074 for id, va in fIdName.items(): 1075 if va == value: 1076 passFunction = True 1077 function_index = id 1078 1079 if passFunction: 1080 pass 1081 else: 1082 # lambda(v, species_3, species_4, v * species_3 * species_4) 1083 # for modifier(2) Constant flux 1084 # 1085 # lambda(k1,S1,M1,h, k1*S1*1/(1+M1**h)) 1086 # for inhibited(1) Mass action 1087 # 1088 myFkey = 'function_' + str(self.indexFuDe) 1089 fIdName[myFkey] = value 1090 functionDefinition.setName(value) 1091 functionFormula = 'lambda(%s, %s)' % (functionVariables.strip(','), mathFormula) 1092 functionDefinition.setMath(lSBML.parseFormula(functionFormula)) 1093 1094 # important debug output 1095 #print "reaction: %s\tfunctiondefinition: %s\t fuId: %s"%(reaction, value, myFkey) 1096 #print "functionFormula: %s"%functionFormula 1097 #print "math: %s"%functionDefinition.getMath() 1098 #print "value: %s\n"%value 1099 1100 # keep the unique relation between law and assignModel 1101 if assignModel != '': 1102 key = self.idNameReactions[reaction] + ':' + assignModel 1103 else: 1104 if reaction in self.idNameReactions.keys(): 1105 key = self.idNameReactions[reaction] 1106 else: 1107 key = '' 1108 #functionDefinition = None 1109 1110 # if all reaction is in the same compartment 1111 if self.inTheSameCompartment: 1112 1113 # if the function not exists 1114 if function_index == 'existIndex': 1115 formula = self.thisCompartment + '*function_' + str(self.indexFuDe) + formula 1116 self.indexFuDe += 1 1117 1118 # if the function already exists in muster cps file 1119 else: 1120 formula = self.thisCompartment + '*' + function_index + formula 1121 1122 # if in different compartments 1123 else: 1124 if function_index == 'existIndex': 1125 formula = 'function_' + str(self.indexFuDe) + formula 1126 self.indexFuDe += 1 1127 else: 1128 formula = function_index + formula 1129 if key != '': 1130 # if entry parameter has reaction_6, law: MM, assignModel: ca1 1131 # key: v6:ca1 value: Michaelis-Menten 1132 # this essential dict will be used in def assignFunctionNames 1133 self.nameKineticNameMap[key] = value 1134 1135 #self.__model.reactionId[reaction].setName(value) 1136 kineticLaw.setFormula(formula) 1137 return (kineticLaw, functionDefinition)

1138 1139

1140 - def __createStructureSBML(self):

1141 """ 1142 Creates new models from given graphs(usually processed by Generator.removeNode()) 1143 and the old model. 1144 1145 @rtype: {str:libSBML.SBMLDocument} 1146 @return: a dict of SBMLDocuments containing the created models 1147 """ 1148 changedModels = {} 1149 sbmlModel = self.__model.getSbmlDocument().getModel() 1150 for removedNodes, graph in self.__changedDG.items(): 1151 newDocument = lSBML.SBMLDocument() 1152 # Workaround to change the namespaces in the SBML tag, Celldesigner inserts its namespace into the SBML tag 1153 # although this is not possible with libSBML, so we have to convert the model to minidom, which does the trick. 1154 newDocumentDom = minidom.parseString(unicode('<?xml version="1.0" encoding="UTF-8"?>\n%s' % (newDocument.toSBML()))) 1155 1156 sbml = newDocumentDom.getElementsByTagName('sbml')[0] 1157 oldNamespaces = self.__model.getSbmlDocument().getNamespaces() 1158 # newNamespaces = newDocument.getNamespaces() 1159 1160 for namespace in range(oldNamespaces.getLength()): 1161 prefix = oldNamespaces.getPrefix(namespace) 1162 uri = oldNamespaces.getURI(namespace) 1163 if prefix: 1164 sbml.setAttribute('xmlns:' + str(prefix), uri) 1165 1166 #can't overwrite libsbml.Document object due to error in libsbml so we need a new variable 1167 reader = lSBML.SBMLReader() 1168 newDocumentWNS = reader.readSBMLFromString(str(newDocumentDom.toxml())) 1169 1170 newModel = newDocumentWNS.createModel() 1171 newModel.setName(sbmlModel.getName()) 1172 1173 #Copy Annotation from old model and add the information which nodes have been removed 1174 removedNodesList = removedNodes.replace('[', '').replace(']', '').replace('\'', '').replace(' ', '').split(',') 1175 1176 #use xmlobject to manipulate the xml inside the annotation tag and to add namespaces 1177 writer = lSBML.SBMLWriter() 1178 rootObject = XMLFile(raw=writer.writeToString(newDocumentWNS)) 1179 1180 try: 1181 annotationObject = rootObject.sbml.model.annotation 1182 except Exception, e: 1183 rootObject.sbml.model._addNode('annotation') 1184 annotationObject = rootObject.sbml.model.annotation 1185 #raise generatorError(e, 'Model annotation in rootfile is missing or not valid.') 1186 try: 1187 annotationObject.modelMaGe 1188 except: 1189 annotationObject._addNode('modelMaGe') 1190 annotationObject.modelMaGe.xmlns = "http://www.modelmage.org/" 1191 try: 1192 #if there is already a node listOfRemovedNodes 1193 annotationObject.modelMaGe.listOfRemovedNodes 1194 except Exception, e: 1195 annotationObject.modelMaGe._addNode('listOfRemovedNodes') 1196 #raise generatorError(e, 'ModelMage annotation in rootfile is missing or not valid.') 1197 1198 noderecorde = '' 1199 for removed in removedNodesList: 1200 node = annotationObject.modelMaGe.listOfRemovedNodes._addNode('node') 1201 node.id = removed 1202 noderecorde += '<p>node id= ' + removed + '</p>' 1203 #annotationStream = lSBML.XMLInputStream(str(annotationObject._toxml()),0) 1204 #annotation = lSBML.XMLNode(annotationStream) 1205 1206 newModel.setAnnotation(str(annotationObject._toxml())) 1207 1208 ###Copy elements that are never altered by the class:### 1209 #Add notes from old model 1210 notes = sbmlModel.getNotes() 1211 notes = "<notes>\ 1212 <html xmlns='http://www.w3.org/1999/xhtml'>\ 1213 <head> <meta name='qrichtext' content='1'/></head>\ 1214 <body style='font-size:8pt;font-family:MS Shell Dlg'>\ 1215 <p>listOfRemoveNodes:</p>" 1216 notes = notes + noderecorde 1217 notes = notes + "</body> </html> </notes>" 1218 newModel.setNotes(notes) 1219 1220 #Add compartments from old model 1221 compartments = sbmlModel.getListOfCompartments() 1222 for compartment in compartments: 1223 newModel.addCompartment(compartment) 1224 try: 1225 #Add FormulaUnitsDatas from old model 1226 FormulaUnitsDatas = sbmlModel.getListFormulaUnitsData() 1227 for FormulaUnitsData in FormulaUnitsDatas: 1228 newModel.addFormulaUnitsData(FormulaUnitsData) 1229 except: 1230 pass 1231 1232 #Add constraints from old model 1233 constraints = sbmlModel.getListOfConstraints() 1234 for constraint in constraints: 1235 newModel.addConstraint(constraint) 1236 1237 #Add InitialAssingmentss from old model 1238 InitialAssignments = sbmlModel.getListOfInitialAssignments() 1239 for InitialAssignment in InitialAssignments: 1240 newModel.addInitialAssignment(InitialAssignment) 1241 1242 #Add UnitDefinitions from old model 1243 UnitDefinitions = sbmlModel.getListOfUnitDefinitions() 1244 for UnitDefinition in UnitDefinitions: 1245 newModel.addUnitDefinition(UnitDefinition) 1246 1247 #Add parameters from old model 1248 parameters = sbmlModel.getListOfParameters() 1249 for parameter in parameters: 1250 newModel.addParameter(parameter) 1251 1252 #Add events from old model: 1253 events = sbmlModel.getListOfEvents() 1254 for event in events: 1255 newModel.addEvent(event) 1256 1257 #Add rules from old model: 1258 rules = sbmlModel.getListOfRules() 1259 for rule in rules: 1260 if rule.getVariable() in graph.nodes(): 1261 newModel.addRule(rule) 1262 1263 #Add functions from old model: 1264 functionDefinitions = sbmlModel.getListOfFunctionDefinitions() 1265 for functionDefinition in functionDefinitions: 1266 newModel.addFunctionDefinition(functionDefinition) 1267 1268 ###Species and reactions are more complicated, because they can be subject to change.### 1269 #Add species and reactions from old model: 1270 nodes = graph.nodes() 1271 nodes.sort() 1272 for node in nodes: 1273 #Node is a species: 1274 if (self.__model.speciesId.has_key(node)): 1275 newModel.addSpecies(self.__model.speciesId[node]) 1276 1277 #Node is a reaction: 1278 if (self.__model.reactionId.has_key(node) or \ 1279 node not in self.__model.reactionId.keys() and node not in self.__model.speciesId.keys()): #Update reactions or add reactions created by removeNode() 1280 reactants = [] 1281 products = graph.successors(node) 1282 modifiers = [] 1283 for edge in graph.in_edges(node): 1284 if edge[2] == 'reactant': 1285 reactants.append(edge[0]) 1286 elif edge[2] == 'modifier': 1287 modifiers.append(edge[0]) 1288 reactants.sort() 1289 products.sort() 1290 modifiers.sort() 1291 1292 if node not in self.__model.reactionId.keys(): 1293 self.__updateReaction(node, reactants, products, modifiers, newModel, 1) 1294 else: 1295 self.__updateReaction(node, reactants, products, modifiers, newModel, 0) 1296 1297 #Add created Model to dict 1298 changedModels[removedNodes] = newDocumentWNS 1299 1300 return changedModels

1301

1302 - def __updateReaction(self, node, reactants, products, modifiers, model, add):

1303 """ 1304 Decides if a reaction must be updated or added according to the changes in the graph or not. 1305 If an update is neccessary the given arguments are passed to the reaction as products 1306 modifiers and reactants. 1307 1308 @type node: string 1309 @param node: new/changed reaction in the graph 1310 1311 @type reactants: [string] 1312 @param reactants: list of reactants of the reaction 1313 1314 @type products: [string] 1315 @param products: list of products of the reaction 1316 1317 @type modifiers: [string] 1318 @param modifiers: list of modifiers of the reaction 1319 1320 @type model: libSBML.SBMLModel 1321 @param model: the model the reactions has to be added to 1322 1323 @type add: bool 1324 @param add: 1=add reaction, 0=update reaction 1325 """ 1326 if not add and \ 1327 [reac.getSpecies() for reac in self.__model.reactionId[node].getListOfReactants()] == reactants and \ 1328 [prod.getSpecies() for prod in self.__model.reactionId[node].getListOfProducts()] == products and \ 1329 [mod.getSpecies() for mod in self.__model.reactionId[node].getListOfModifiers()] == modifiers: 1330 model.addReaction(self.__model.reactionId[node]) 1331 else: 1332 reaction = model.createReaction() 1333 reaction.setId(node) 1334 if self.__model.reactionId.has_key(node): 1335 nName = self.__model.reactionId[node].getName() 1336 #newNname = self.pat.sub('v_reaction',nName) 1337 newNname = self.pat.sub('v_re', nName) 1338 reaction.setName(newNname) 1339 ##reaction.setName(self.__model.reactionId[node].getName()) 1340 else: 1341 #node: reaction_reaction_3_reaction_4 1342 #changed to v_reaction_3_reaction_4 1343 #newNname = self.pat.sub('v_reaction',node) 1344 newNname = self.pat.sub('v_re', node) 1345 reaction.setName(newNname) 1346 #reaction.setName(node) 1347 reaction.setReversible(0) 1348 1349 #string to insert mass action kinetics into the reaction 1350 kineticLaw = '' 1351 for product in products: 1352 newProduct = model.createProduct() 1353 newProduct.setSpecies(product) 1354 newProduct.setStoichiometry(self.__stoich[(node, product)]) 1355 for reactant in reactants: 1356 newReactant = model.createReactant() 1357 newReactant.setSpecies(reactant) 1358 newReactant.setStoichiometry(self.__stoich[(reactant, node)]) 1359 for i in range(int(self.__stoich[(reactant, node)])): 1360 kineticLaw = kineticLaw + reactant + '*' 1361 1362 # filter out duplicates: 1363 m = [] 1364 for x in modifiers: 1365 if x not in m: 1366 m.append(x) 1367 modifiers = m 1368 1369 for modifier in modifiers: 1370 model.createModifier().setSpecies(modifier) 1371 kineticLaw = kineticLaw + modifier + '*' 1372 #remove trailing '*' from formula 1373 kineticLaw = kineticLaw.strip('*') 1374 1375 #insert kinetic law into reaction 1376 model.createKineticLaw().setFormula(kineticLaw)

1379 """ 1380 Holds a Copasifile and has methods to manipulate it. 1381 The class is used to merge the root copasifile with the submodels. 1382 """ 1383

1384 - def __init__(self, filename, history, map=None):

1385 """ 1386 1387 @type filename: str 1388 @param filename: path to the copasifile 1389 1390 @type history: dict 1391 @param history: the node history of the candidate models 1392 """ 1393 self.cpsFile = minidom.parse(filename) 1394 self.filename = filename 1395 self.graphHistory = history 1396 self.idNameMap = map 1397 # used in __updateReaction, importTask, assignFunctionNames 1398 self.pat = re.compile('reaction_re')

1399

1400 - def __addFitItem(self, importFile, problem, model, ReactionName, parameter):

1401 fitItem = minidom.parseString(self.__fitItem(model, ReactionName, parameter)) 1402 fitItem = fitItem.firstChild.firstChild 1403 groups = problem.getElementsByTagName("ParameterGroup") 1404 for group in groups: 1405 if group.attributes['name'].value == u"OptimizationItemList": 1406 group.appendChild(fitItem)

1407

1408 - def __removeExperimentalPar(self, valueDict, parameter):

1409 valueDict['Metabolite'] = valueDict['Vector'].split('[')[1][:-1] # this is the metabolite which belongs to the experimental data 1410 if valueDict['Metabolite'] in self.__getMetaboliteNames(): 1411 return 0 1412 elif valueDict['Metabolite'] in self.__getModelValues(): 1413 return 0 1414 else: 1415 for child in parameter.parentNode.childNodes: 1416 if child.nodeType == 1 and child.attributes['value'].value == u'2': 1417 child.attributes['value'].value = u'0' 1418 parameter.parentNode.removeChild(parameter) 1419 return 1

1420 1421

1422 - def __getTask(self, type):

1423 """ 1424 Returns all tasks of type 'type'. For example: 1425 1426 >>>self.getTask('parameterFitting') 1427 1428 would return a minidom object containing the corresponding task. 1429 1430 @type type: unicode 1431 @param type: type of the task to return 1432 1433 @rtype: minidom.Node 1434 @return: task of type type 1435 """ 1436 tasks = self.cpsFile.getElementsByTagName('Task') 1437 task = [x for x in tasks if x.attributes['type'].value == type] 1438 if task: 1439 return task[0] 1440 else: 1441 return 0

1442

1443 - def __getReport(self, type):

1444 """ 1445 Returns all tasks of type 'type'. 1446 For example: 1447 1448 >>>self.getReport('parameterFitting') 1449 1450 would return a minidom object containing the corresponding report. 1451 1452 @type type: unicode 1453 @param type: type of the report to return 1454 1455 @rtype: minidom.Node 1456 @return: report of type type 1457 """ 1458 reports = self.cpsFile.getElementsByTagName('Report') 1459 for report in reports: 1460 if report.attributes['taskType'].value == type: return report 1461 else: 1462 return False

1463 1464

1465 - def __getReactionNames(self):

1466 """ 1467 Return the names of all reactions of the model in a list. 1468 1469 @rtype: [str] 1470 @return list of reaction names 1471 """ 1472 return [str(x.attributes['name'].value) for x in self.cpsFile.getElementsByTagName('Reaction')]

1473 1474 1475

1476 - def __getReactionIds(self):

1477 """ 1478 Return the ids of all reactions of the model in a list. 1479 1480 @rtype: [str] 1481 @return list of reaction ids 1482 """ 1483 return [str(x.attributes['key'].value) for x in self.cpsFile.getElementsByTagName('Reaction')]

1484 1485

1486 - def __getMetaboliteNames(self):

1487 """ 1488 Return the names of all metabolites/species of the model in a list. 1489 1490 @rtype: [str] 1491 @return list of metabolite names 1492 """ 1493 return [str(x.attributes['name'].value) for x in self.cpsFile.getElementsByTagName('Metabolite')]

1494

1495 - def __getMetaboliteIds(self):

1496 """ 1497 Return the ids of all metabolites/species of the model in a list. 1498 1499 @rtype: [str] 1500 @return list of metabolite ids 1501 """ 1502 return [str(x.attributes['key'].value) for x in self.cpsFile.getElementsByTagName('Metabolite')]

1503

1504 - def __getModelValues(self):

1505 """ 1506 Return the names of all global quantities of the model in a list. 1507 1508 @rtype: [str] 1509 @return list of global quantities 1510 """ 1511 return [str(x.attributes['name'].value) for x in self.cpsFile.getElementsByTagName('ModelValue')]

1512

1513 - def __getReactionParameters(self, reaction):

1514 """ 1515 Returns all parameter names that belong to a specific reaction. 1516 1517 @type reaction: str 1518 @param reaction: a reaction belonging to the model 1519 @rtype: [str] 1520 @return: a list of parameternames 1521 """ 1522 reactions = self.cpsFile.getElementsByTagName('Reaction') 1523 found = 0 1524 for x in reactions: 1525 if x.attributes['name'].value == reaction: 1526 found = x 1527 if found == 0: 1528 return [] # if the reaction isn't found return an empty list 1529 1530 return [str(x.attributes['name'].value) for x in found.getElementsByTagName('Constant')]

1531

1532 - def __getGlobalParameters(self):

1533 """ 1534 Returns all parameter names that belong to a specific reaction. 1535 1536 @rtype: [str] 1537 @return: a list of parameternames 1538 """ 1539 values = self.cpsFile.getElementsByTagName('ModelValue') 1540 return [str(x.attributes['name'].value) for x in values]

1541

1542 - def __fitItem(self, model, ReactionName, parameter):

1543 fitItem = """<document><ParameterGroup name="FitItem"> 1544 <Parameter name="ObjectCN" type="cn" value="CN=Root,Model=%s,Vector=Reactions[%s],ParameterGroup=Parameters,Parameter=%s,Reference=Value"/> 1545 <Parameter name="LowerBound" type="cn" value="0.00001"/> 1546 <Parameter name="UpperBound" type="cn" value="100"/> 1547 <Parameter name="StartValue" type="float" value="1"/> 1548 <ParameterGroup name="Affected Experiments"> 1549 </ParameterGroup> 1550 </ParameterGroup></document>""" % (model, ReactionName, parameter) 1551 return fitItem

1552

1553 - def getExperimentalDataFiles(self):

1554 """ 1555 Gets the paths to experimental data files specified for a parameter estimation and returns them in a list. 1556 1557 @rtype: list 1558 @return: ist of paths to the specified experimental data files 1559 """ 1560 estimationTask = self.__getTask('parameterFitting') 1561 if estimationTask: 1562 filename = [x for x in estimationTask.getElementsByTagName('Parameter') if x.attributes['name'].value == 'File Name'] 1563 else: 1564 raise generatorError('', "No existing parameter estimation task in {root}.cps") 1565 return [x.attributes['value'].value for x in filename]

1566

1567 - def __getFitItems(self, task):

1568 1569 problem = task.getElementsByTagName("Problem")[0] 1570 parameters = problem.getElementsByTagName("Parameter") 1571 return [x for x in parameters if x.attributes['name'].value == u"ObjectCN"]

1572

1573 - def __getExperimentalData(self, task):

1574 1575 problem = task.getElementsByTagName("Problem")[0] 1576 parameters = problem.getElementsByTagName("Parameter") 1577 return [x for x in parameters if x.attributes['name'].value == u"Object CN"]

1578

1579 - def __createValueDict(self, longString):

1580 items = str(longString).split(',') 1581 valueDict = {} 1582 for pair in items: 1583 key, value = pair.split('=') 1584 valueDict[key] = value 1585 1586 vectorValue = valueDict['Vector'].split('[')[0] # determines if the parameter is a reaction or not 1587 ###Check what kind of parameter(global/local) will be estimated: 1588 #if a parameter is local(reactions) we have to check if it is still there 1589 if vectorValue == 'Reactions': 1590 #reaction the fitItem belongs to: 1591 valueDict['Reaction'] = valueDict['Vector'].split('[')[1][:-1] 1592 1593 valueDict['Reaction'] = valueDict['Reaction'].replace("\\", "") # this corrects a parsing error with non ASCII characters 1594 1595 #if the parameter is global(value) than it has to stay in the task: 1596 elif vectorValue == 'Values': 1597 valueDict['Reaction'] = None 1598 return valueDict

1599 1600

1601 - def assignFunctionNames(self, functionReactionMap, functionsIdName, modelName, KineticNameMap):

1602 ''' 1603 gives kinetics with a changed formula the corresponding name 1604 1605 @type functionReactionMap: map ( produced in def getFunction_Id_Name_XML in modelmage.py) 1606 @param functionReactionMap: {function_id : reaction_id} 1607 1608 @type functionsIdName: map (produced in def getMasterSBMLFunctions(input) in modelmage.py) 1609 @param functionsIdName: {function_id : function_name} 1610 1611 @type modelName: str 1612 @param modelName: current name of candidate model 1613 1614 @type KineticNameMap: map (produced in def __lawtoSBML()) 1615 @param KineticNameMap: {function_id: kineticName} 1616 1617 ''' 1618 1619 # print "\nmodelName: %s\nfunctionReactionMap: %s"%(modelName, functionReactionMap) 1620 # corresponding substitution in __updateReaction 1621 for key, value in functionReactionMap.items(): 1622 #print "key: %s\tvalue: %s"%(key,value) 1623 #if 'reaction_reaction' in key: 1624 if 'reaction_re' in key: 1625 newKey = self.pat.sub('v_re', key) 1626 functionReactionMap[newKey] = value 1627 del functionReactionMap[key] 1628 1629 # print "\nfunctionsIdName:" 1630 # for key, value in functionsIdName.items(): 1631 # print "key: %s\tvalue: %s"%(key,value) 1632 # print "\nkineticNameMap:" 1633 # for key, value in KineticNameMap.items(): 1634 # print "key: %s\tvalue: %s"%(key,value) 1635 1636 # get all functions from the .cps file 1637 functions = self.cpsFile.getElementsByTagName('Function') 1638 dummyIndex = 1 1639 for f in functions: 1640 newF = f 1641 functionName = f.attributes['name'].value 1642 reff = functionName.split('_') 1643 #print "TTTTT function Name: %s"%functionName 1644 vReaction = '_'.join(reff[2:]) 1645 #print "TTTTT vReaction: %s"%vReaction 1646 if vReaction in functionReactionMap.keys(): 1647 func = functionReactionMap[vReaction] 1648 funcName = functionsIdName[func] 1649 #print "funcName: %s"%funcName 1650 newF.attributes['name'].value = funcName + '(' + vReaction + ')' 1651 f.parentNode.replaceChild(newF, f)

1652 1653

1654 - def importTask(self, file, taskType):

1655 """ 1656 Adds a task to the cps file and overwrites tasks of the same type. 1657 The method automatically sets attributes that relate to other elements 1658 in the file(e.g. report reference) to the correct values. 1659 1660 Filenames for resultfiles are automatically set to: filname_est.txt 1661 1662 @type file: str 1663 @param file: the file from which the task will be imported 1664 1665 @type taskType: str 1666 @param taskType: the type of task that shall be imported 1667 """ 1668 #open the file as a new CpsHandler 1669 removed = [] 1670 importFile = CpsHandler(file, {}) 1671 task = importFile.__getTask(taskType) 1672 model = self.cpsFile.getElementsByTagName('Model')[0].attributes['name'].value 1673 1674 #set the correct new key: 1675 task.attributes['key'].value = self.__getTask(taskType).attributes['key'].value 1676 #Tag elements: scheduled, updateModel should always be assigned with true 1677 task.attributes['scheduled'].value = 'true' 1678 task.attributes['updateModel'].value = 'true' 1679 1680 #set the report reference: 1681 reportReference = self.__getTask(taskType).getElementsByTagName('Report')[0].attributes['reference'].value 1682 task.getElementsByTagName('Report')[0].attributes['reference'].value = reportReference 1683 #Give a new name for the file that holds reports because the otherwise the file would be overwritten by 1684 #each of the models. 1685 filename = os.path.split(self.filename)[-1] 1686 filename = filename.split('.')[0] + '_est.txt' 1687 task.getElementsByTagName('Report')[0].attributes['target'].value = filename 1688 1689 #set the references to timecourse and steady state: 1690 steadyState = '' 1691 timeCourse = '' 1692 for parameter in self.__getTask(taskType).getElementsByTagName('Parameter'): 1693 if parameter.attributes['name'].value == u"Steady-State": 1694 steadyState = parameter.attributes['value'].value 1695 if parameter.attributes['name'].value == u"Time-Course": 1696 timeCourse = parameter.attributes['value'].value 1697 1698 for parameter in task.getElementsByTagName('Parameter'): 1699 if parameter.attributes['name'].value == u"Steady-State": 1700 parameter.attributes['value'].value = steadyState 1701 if parameter.attributes['name'].value == u"Time-Course": 1702 parameter.attributes['value'].value = timeCourse 1703 if parameter.attributes['name'].value == 'File Name': 1704 expDataName = parameter.attributes['value'].value 1705 if expDataName.rfind('/') != -1: 1706 index = expDataName.rfind('/') 1707 newDataName = expDataName[(index + 1):] 1708 parameter.attributes['value'].value = newDataName 1709 1710 #get the parameter tags of the tasks problem definition 1711 problem = task.getElementsByTagName("Problem")[0] 1712 #parameters = problem.getElementsByTagName("Parameter") 1713 1714 # Update the fitItem parameters 1715 importFitItems = copy.deepcopy(self.__getFitItems(task)) 1716 for fitItem in self.__getFitItems(task): 1717 #split the long identifierstring substrings and create a dict of these 1718 valueDict = self.__createValueDict(fitItem.attributes["value"].value) 1719 if valueDict.has_key('Reaction'): 1720 if valueDict['Reaction'] == None: 1721 continue 1722 1723 # check if the reaction and it's parameter still exist in the cps file if they both do continue 1724 if valueDict['Reaction'] in self.__getReactionNames() and \ 1725 valueDict['Parameter'] in self.__getReactionParameters(valueDict['Reaction']): 1726 continue 1727 1728 # remove parameter if it is not in the new file 1729 elif valueDict['Reaction'] in self.__getReactionNames() and not \ 1730 valueDict['Parameter'] in self.__getReactionParameters(valueDict['Reaction']): 1731 for importParam in importFile.__getReactionParameters(valueDict['Reaction']): 1732 if importParam not in self.__getReactionParameters(valueDict['Reaction']): 1733 try: 1734 fitItem.parentNode.parentNode.removeChild(fitItem.parentNode) 1735 except: 1736 #print 'could not remove fitItem %s' % importParam 1737 pass 1738 1739 # remove the fititem: 1740 else: 1741 fitItem.parentNode.parentNode.removeChild(fitItem.parentNode) 1742 removed.append((valueDict['Reaction'], valueDict['Parameter'])) 1743 1744 # Update the Experimental Data parameters for local parameters 1745 for expData in self.__getExperimentalData(task): 1746 valueDict = self.__createValueDict(expData.attributes["value"].value) 1747 if valueDict.has_key('Reaction'): 1748 if expData.attributes['name'].value == u"Object CN": 1749 self.__removeExperimentalPar(valueDict, expData) 1750 1751 # add parameters of reactions, which have an ancestor in the imported task, to the new task 1752 newReactions = [x for x in self.__getReactionNames() if x not in importFile.__getReactionNames()] 1753 1754 # get the history of reactions and reformat it 1755 flatHistory = [] 1756 for item in [x.items() for x in self.graphHistory.values()]: 1757 flatHistory.extend(item) 1758 1759 history = {} 1760 for tupel in flatHistory: 1761 history[tupel[0]] = [] 1762 for elem in tupel[1]: 1763 try: 1764 history[tupel[0]].append(self.idNameMap[elem]) 1765 except KeyError: 1766 history[tupel[0]].append(elem) 1767 1768 # get all reactions which have fitted parameters and put them in a list 1769 fittedReactions = [] 1770 for item in importFitItems: 1771 valueDict = self.__createValueDict(item.attributes["value"].value) 1772 if valueDict.has_key('Reaction'): 1773 fittedReactions.append(valueDict['Reaction']) 1774 addParametersOf = [] 1775 1776 # corresponding to the substitution in function __updateReaction() 1777 for key, value in history.items(): 1778 if 'reaction_reaction' in key: 1779 # newKey = self.pat.sub('v_reaction',key) 1780 newKey = self.pat.sub('v_re', key) 1781 history[newKey] = value 1782 del history[key] 1783 else: 1784 continue 1785 #print "key: %s\tvalue: %s"%(key, value) 1786 # find reactions that have ancestors that were estimated in the master model 1787 # and add all of their parameters 1788 for reaction in newReactions: 1789 for ancestor in history[reaction]: 1790 if ancestor in fittedReactions: 1791 if reaction not in addParametersOf: 1792 addParametersOf.append(reaction) 1793 1794 # if a reaction has new parameters, add them to the estimation: 1795 for reactionName in self.__getReactionNames(): 1796 if reactionName in fittedReactions: 1797 for parameter in [x for x in self.__getReactionParameters(reactionName) if x not in importFile.__getReactionParameters(reactionName)]: 1798 # if parameter not in self.__getGlobalParameters(): 1799 self.__addFitItem(importFile, problem, model, reactionName, parameter) 1800 1801 # add fitItem parameters for newly generated reactions to the task: 1802 for reactionName in addParametersOf: 1803 ReactionParameters = self.__getReactionParameters(reactionName) 1804 for parameter in ReactionParameters: 1805 self.__addFitItem(importFile, problem, model, reactionName, parameter) 1806 1807 #update parameter for Evolutionary Programming, Genetic Algorithm, Genetic Algorithm SR, Evolution Strategy (SRES) 1808 EvolutionGeneticMethods = ('Evolutionary Programming', 'Evolutionary Strategy (SRES)', 'Genetic Algorithm', 'Genetic Algorithm SR') 1809 numberOfFitItems = len(self.__getFitItems(task)) 1810 #print "number of parameter: %s\n"%(numberOfFitItems) 1811 for methode in task.getElementsByTagName('Method'): 1812 #print "Method name: %s type: %s\n"%(methode.attributes['name'].value, methode.attributes['type'].value) 1813 methodName = methode.attributes['name'].value 1814 if methodName in EvolutionGeneticMethods: 1815 for parameter in reversed(task.getElementsByTagName('Method')[0].getElementsByTagName('Parameter')): 1816 #print "parameter: %s\tvalue: %s\n"%(parameter.attributes['name'].value, parameter.attributes['value'].value) 1817 if parameter.attributes['name'].value == 'Population Size': 1818 PopulationSizeValue = 10 * numberOfFitItems 1819 parameter.attributes['value'].value = str(PopulationSizeValue) 1820 if parameter.attributes['name'].value == 'Number of Generations': 1821 parameter.attributes['value'].value = str(10 * PopulationSizeValue) 1822 if parameter.attributes['name'].value == 'Seed': 1823 parameter.attributes['value'].value = '1' 1824 1825 1826 #insert updated task into file: 1827 for oldTask in self.cpsFile.getElementsByTagName("Task"): 1828 if oldTask.attributes["type"].value == taskType: 1829 oldTask.parentNode.replaceChild(task, oldTask) 1830 1831 ###Add corresponding report to the file. 1832 taskReport = '' 1833 for report in importFile.cpsFile.getElementsByTagName('ListOfReports')[0].getElementsByTagName('Report'): 1834 if report.attributes['taskType'].value == taskType: 1835 taskReport = report 1836 for report in self.cpsFile.getElementsByTagName('ListOfReports')[0].getElementsByTagName('Report'): 1837 if report.attributes['taskType'].value == taskType: 1838 report.parentNode.replaceChild(taskReport, report) 1839 1840 return removed

1841

1842 - def __encodeAttributes(self, string):

1843 """ 1844 Converts encoding of attributes. There are some characters which Copasi interprets as 1845 html encoding, these characters are switched to this encoding. 1846 1847 @type string: str 1848 @param string: the string that will be changed 1849 1850 @rtype: str 1851 @return: a new string with exchanged special characters 1852 """ 1853 newString = r'' 1854 isAttribute = False 1855 #every special character inside an attribute is switched converted to its html representation 1856 for char in string: 1857 if not isAttribute and char == '"': 1858 isAttribute = True 1859 newString += char 1860 continue 1861 if isAttribute and char == '"': 1862 isAttribute = False 1863 newString += char 1864 continue 1865 if isAttribute: 1866 newString += { '&': "&", 1867 '<': "<", 1868 '\"': """, 1869 '\t': "	", 1870 '\n': "
", 1871 }.get(char, char) #char is a default if it is not found in dict 1872 else: 1873 newString += char 1874 1875 return newString

1876

1877 - def writeCPS(self, filename):

1878 """ 1879 Writes the file to the path filename. 1880 1881 @type filename: str 1882 @param filename: path where the file will be written 1883 1884 @rtype: bool 1885 @return: True if file was written, else False 1886 """ 1887 1888 try: 1889 self.cpsFile.normalize() 1890 f = open(filename, 'w') 1891 output = self.__encodeAttributes(self.cpsFile.toxml()) 1892 1893 if not 'name="Seperator" type="string" value="	"' in output: 1894 output = self.cpsFile.toxml().encode('utf8', 'ignore') 1895 1896 f.write(output) 1897 f.close() 1898 return True 1899 except Exception, e: 1900 print e 1901 return False

Source Code for Module modelmage.generator