libTI2.model.model

11 - def __make_sbmldoc_available__(self, force=False):

12 """Currently not available, since parallel computing is disabled""" 13 if not self.use_pp and not force: return 14 self.sbmldoc = libsbml.readSBMLFromString(self.__sbmldoc_string__) 15 del self.__sbmldoc_string__

16 - def __hide_sbmldoc__(self, force=False):

17 """Currently not available, since parallel computing is disabled""" 18 if not self.use_pp and not force: return 19 self.__sbmldoc_string__ = libsbml.writeSBMLToString(self.sbmldoc) 20 del self.sbmldoc

21 - def load_sbml(self,filename):

22 self.filename = filename 23 path = filename.split(os.sep) 24 if len(path) == 1: path.insert(0,".") 25 self.path_to_xmlfile = os.sep.join( path[:-1] ) 26 self.path_to_localfiles = self.path_to_xmlfile + os.sep + ".".join( filename.split(os.sep)[-1].split(".")[:-1] ) 27 self.aih = additional_information_handler.additional_information_handler(self) 28 if self.aih.exists: 29 filename = self.aih.filename 30 print "Reloading cached file" 31 else: 32 try: 33 os.mkdir(self.path_to_localfiles) 34 except OSError: 35 # not necessary to create a new directory 36 pass 37 r = libsbml.SBMLReader() 38 try: 39 d = r.readSBML(filename) 40 except IOError: 41 print "Could not find file",filename 42 return None 43 return d

44 - def __get_working_filename__():

45 return self.filename+".working.xml"

46 - def __get_sbml_element_from_list_by_name__(self,name,list):

47 for e in list: 48 if e.getId() == name: 49 return e 50 return None

51 - def __add_reaction_math_to_differential_equation__(self,math,species,dictionary,is_positive_summand=True):

52 # test for constant 53 species_reference = self.__get_sbml_element_from_list_by_name__(species,self.sbmldoc.getModel().getListOfSpecies()) 54 if species_reference.getConstant() or species_reference.getBoundaryCondition(): return dictionary 55 if species not in dictionary: 56 dictionary[ species ] = "0" 57 if is_positive_summand: 58 premath = " + " 59 else: 60 premath = " - " 61 dictionary[ species ] = dictionary[ species ] + premath + math 62 return dictionary

63 - def mathml2scipy(self,formula,replace_dictionary={}):

64 formula = " "+formula+" " 65 formula = self.function_definition_replace(formula) 66 formula = self.assignment_rule_replace(formula) 67 for key in rre_generator.scipy_replaces: 68 formula = formula.replace( " "+key+" ", " scipy."+rre_generator.scipy_replaces[key]+" " ) 69 for key in replace_dictionary: 70 formula = formula.replace( " "+key+" ", " "+str(replace_dictionary[key])+" " ) 71 return formula

72 - def assignment_rule_replace(self,formula):

73 for ar in [ x for x in self.sbmldoc.getModel().getListOfRules() if x.isAssignment()]: 74 if " "+ar.getVariable()+" " in formula: 75 formula = formula.replace(" "+ar.getVariable()+" "," ( "+function_manipulation.include_spaces_in_formula(ar.getFormula())+" ) ") 76 return formula

77 - def function_definition_replace( self,formula ):

78 for fd in self.sbmldoc.getModel().getListOfFunctionDefinitions(): 79 if " "+fd.getId()+" " in formula: 80 paramlist = [] 81 for i in range(fd.getNumArguments()): 82 paramlist.append( fd.getArgument(i).getName() ) 83 matchstring = "\s*"+fd.getId()+"\s*\("+ ",".join([".*"]*len(paramlist)) +"\)\s*" 84 match = re.search(matchstring,formula) 85 fd_formula = function_manipulation.include_spaces_in_formula(libsbml.formulaToString(fd.getMath().getRightChild())) 86 while match: 87 matching = match.group(0) 88 vars_in_fd = [x.strip() for x in matching.split("(")[1].strip().strip(")").split(",")] 89 this_fd_formula = copy.deepcopy(fd_formula) 90 for paramindex in range(len(paramlist)): 91 this_fd_formula = this_fd_formula.replace(" "+paramlist[paramindex]+" "," "+vars_in_fd[paramindex]+" ") 92 formula = formula.replace(matching,this_fd_formula) 93 94 match = re.search(matchstring,formula) 95 96 return formula

97 - def load_species(self):

98 self.species = [] 99 self.species_values = {} 100 for sbmlspecies in self.sbmldoc.getModel().getListOfSpecies(): 101 if sbmlspecies.getConstant() or sbmlspecies.getBoundaryCondition(): continue 102 if sbmlspecies.isSetInitialConcentration(): 103 self.species_values[ sbmlspecies.getId() ] = sbmlspecies.getInitialConcentration() 104 elif sbmlspecies.isSetInitialAmount(): 105 self.species_values[ sbmlspecies.getId() ] = sbmlspecies.getInitialAmount() 106 else: 107 self.species_values[ sbmlspecies.getId() ] = 0 108 self.species.append( sbmlspecies.getId() ) 109 self.species.sort()

110 - def load_parameters(self):

111 self.parameter_values = {} 112 self.parameters = [] 113 for sbmlparameter in self.sbmldoc.getModel().getListOfParameters(): 114 self.parameter_values[ sbmlparameter.getId() ] = sbmlparameter.getValue() 115 self.parameters.append( sbmlparameter.getId() ) 116 for compartment in self.sbmldoc.getModel().getListOfCompartments(): 117 self.parameter_values[ compartment.getId() ] = compartment.getSize() 118 self.parameters.append( compartment.getId() ) 119 for reaction in self.sbmldoc.getModel().getListOfReactions(): 120 rname = reaction.getId() 121 for parameter in reaction.getKineticLaw().getListOfParameters(): 122 newname = rname+"__"+parameter.getId() 123 self.parameter_values[newname] = parameter.getValue() 124 self.parameters.append( newname ) 125 for sbmlspecies in self.sbmldoc.getModel().getListOfSpecies(): 126 if sbmlspecies.getConstant() or sbmlspecies.getBoundaryCondition(): 127 if sbmlspecies.isSetInitialConcentration(): 128 self.parameter_values[ sbmlspecies.getId() ] = sbmlspecies.getInitialConcentration() 129 elif sbmlspecies.isSetInitialAmount(): 130 self.parameter_values[ sbmlspecies.getId() ] = sbmlspecies.getInitialAmount() 131 else: 132 self.parameter_values[ sbmlspecies.getId() ] = 0 133 self.parameters.append( sbmlspecies.getId() ) 134 self.parameters.sort()

135 - def load_assignment_rules(self):

136 self.assignment_rules = [] # has to be a list since it is ordered 137 for ar in [x for x in self.sbmldoc.getModel().getListOfRules() if x.isAssignment()]: 138 self.assignment_rules.append( (ar.getVariable(), function_manipulation.clean_up_formula( ar.getFormula() )) )

139 - def __add_stoichiometric_coefficient__(self,stoi_dic,reaction,substrate,coeff):

140 if reaction not in stoi_dic: 141 stoi_dic[reaction] = {} 142 if substrate not in stoi_dic[reaction]: 143 stoi_dic[reaction][substrate] = 0 144 stoi_dic[reaction][substrate] = stoi_dic[reaction][substrate] + coeff

145 - def load_differential_equations(self):

146 self.load_species() 147 self.load_parameters() 148 self.load_assignment_rules() 149 self.differential_equations = {} 150 for species in self.species: 151 self.differential_equations[ species ] = "0" 152 self.reaction_equations = {} 153 self.reactions = [] 154 self.stoichiometry_dic = {} 155 self.reaction_number = {} 156 self.reaction2substrate2stoich = {} 157 rnumber = 0 158 for reaction in self.sbmldoc.getModel().getListOfReactions(): 159 self.reaction_number[reaction.getId()] = rnumber 160 rnumber = rnumber + 1 161 # get stoichiometry for participants 162 s2s,m2s,p2s = {},{},{} 163 for sr in reaction.getListOfReactants(): 164 if sr.getSpecies() in self.species: 165 s2s[ sr.getSpecies() ] = sr.getStoichiometry() 166 for sr in reaction.getListOfModifiers(): 167 if sr.getSpecies() in self.species: 168 m2s[ sr.getSpecies() ] = 1 169 for sr in reaction.getListOfProducts(): 170 if sr.getSpecies() in self.species: 171 p2s[ sr.getSpecies() ] = sr.getStoichiometry() 172 self.reaction2substrate2stoich[ reaction.getId() ] = [ s2s, m2s, p2s ] 173 if not reaction.isSetKineticLaw(): 174 print reaction.getId(),"has no Kinetic Law!" 175 continue 176 # acquiring parameters 177 params = {} 178 for p in reaction.getKineticLaw().getListOfParameters(): 179 ps = p.getId() 180 params[ps] = reaction.getId()+"__"+ps 181 # getting the formula 182 f = reaction.getKineticLaw().getFormula() 183 f = function_manipulation.include_spaces_in_formula( f ) 184 f = self.mathml2scipy( f, params ) 185 # adding it to the reaction equations 186 self.reactions.append( reaction.getId() ) 187 self.reaction_equations[ reaction.getId() ] = f 188 # appending it 189 for s in reaction.getListOfReactants(): 190 self.differential_equations = self.__add_reaction_math_to_differential_equation__( str(s.getStoichiometry())+" * ( "+f+" ) ", s.getSpecies(), self.differential_equations, False ) 191 self.__add_stoichiometric_coefficient__(self.stoichiometry_dic, reaction.getId() , s.getSpecies(), - s.getStoichiometry() ) 192 for s in reaction.getListOfProducts(): 193 self.differential_equations = self.__add_reaction_math_to_differential_equation__( str(s.getStoichiometry())+" * ( "+f+" ) ", s.getSpecies(), self.differential_equations ) 194 self.__add_stoichiometric_coefficient__(self.stoichiometry_dic, reaction.getId() , s.getSpecies(), s.getStoichiometry() ) 195 self.reactions.sort() 196 N = [] 197 for species in self.species: 198 line = [] 199 for reaction in self.reactions: 200 if species in self.stoichiometry_dic[reaction]: 201 coeff = self.stoichiometry_dic[reaction][species] 202 else: 203 coeff = 0 204 line.append( coeff ) 205 N.append(line) 206 self.N = scipy.matrix(N) 207 self.find_conservation_relations()

208

209 - def find_conservation_relations(self):

210 "docstring for find_conservation_relations" 211 # try to import a correct library for singular value decomposition 212 self.conservation_relations = [] 213 if len(self.species) == 0 or len(self.reactions) == 0: return 214 try: 215 import scipy 216 import scipy.linalg 217 svdfunc = scipy.linalg.svd 218 except ImportError: 219 import numpy.linalg 220 svdfunc = numpy.linalg.svd 221 (u,s,v) = svdfunc(self.N) 222 u = scipy.matrix(u) 223 s = scipy.matrix(s) 224 v = scipy.matrix(v) 225 cutoff = max(self.N.shape)*s.max()*rre_generator.precision()() 226 dimnull = len(self.species)-[x>cutoff for x in s.tolist()[0]].count(True) 227 if dimnull == 0: return 228 # else we have conservation relations 229 conservation_rows = u[:,-dimnull:].round(int(-scipy.log10(cutoff))) 230 for j in range(dimnull): 231 conservationstring = " " 232 for i in range(len(self.species)): 233 if conservation_rows[i,j] != 0.0: 234 conservationstring += "+ ( "+str(conservation_rows[i,j])+ " * "+self.species[i]+" ) " 235 result = function_manipulation.evaluate_clean_formula(conservationstring,self.species_values) 236 conservationstring = " ( "+str(-result)+" ) "+conservationstring 237 self.conservation_relations.append( conservationstring )

238 239

240 - def evaluate_formula(self,formula):

241 # test for new module 242 return function_manipulation.evaluate_general_formula(formula,[self.parameter_values,self.species_values]) 243 # end of new module test 244 this_formula = copy.deepcopy(formula) 245 for p in self.parameters: 246 this_formula = this_formula.replace(" "+p+" "," "+str( self.parameter_values[p] )+" ") 247 for s in self.species: 248 this_formula = this_formula.replace(" "+s+" "," "+str( self.species_values[s] )+" ") 249 return eval(this_formula)

250

251 - def __elasticities_for_list__(self, elasticity_dictionary, differential_equation_dictionary, value_dictionary, first_keys, second_keys):

252 for key1 in first_keys: 253 if key1 not in elasticity_dictionary: elasticity_dictionary[key1] = {} 254 normal_change = self.evaluate_formula( differential_equation_dictionary[ key1 ] ) 255 for key2 in second_keys: 256 value_dictionary[key2] = value_dictionary[key2] + self.small_value 257 new_change = self.evaluate_formula( differential_equation_dictionary[ key1 ] ) 258 value_dictionary[key2] = value_dictionary[key2] - self.small_value 259 elasticity_dictionary[key1][key2] = (new_change - normal_change) / self.small_value 260 return elasticity_dictionary

261

262 - def compute_elasticities(self):

263 self.elasticities_dic = {} 264 self.elasticities_dic = self.__elasticities_for_list__(self.elasticities_dic, 265 self.differential_equations, self.parameter_values, 266 self.species, self.parameters) 267 self.elasticities_dic = self.__elasticities_for_list__(self.elasticities_dic, 268 self.differential_equations, self.species_values, 269 self.species, self.species) 270 self.elasticities_dic = self.__elasticities_for_list__(self.elasticities_dic, 271 self.reaction_equations, self.parameter_values, 272 self.reactions, self.parameters) 273 self.elasticities_dic = self.__elasticities_for_list__(self.elasticities_dic, 274 self.reaction_equations, self.species_values, 275 self.reactions, self.species) 276 substrate_elasticities = [] 277 parameter_elasticities = [] 278 for reaction in self.reactions: 279 line = [] 280 for species in self.species: 281 line.append( self.elasticities_dic[reaction][species] ) 282 substrate_elasticities.append(line) 283 line = [] 284 for parameter in self.parameters: 285 line.append( self.elasticities_dic[reaction][parameter] ) 286 parameter_elasticities.append(line) 287 self.substrate_elasticities = scipy.matrix( substrate_elasticities ) 288 self.parameter_elasticities = scipy.matrix( parameter_elasticities )

289

290 - def __mca_to_dic__(self,dic,matrix,rownames,colnames):

291 for i in range(len(rownames)): 292 rowname = rownames[i] 293 if rowname not in dic: 294 dic[rowname] = {} 295 for j in range(len(colnames)): 296 colname = colnames[j] 297 dic[rowname][colname] = matrix[i,j] 298 return dic

299

300 - def mca(self):

301 self.compute_elasticities() 302 self.mca_coefficients = {} 303 Jacobian = self.N*self.substrate_elasticities 304 try: 305 Inverse_Jacobian = Jacobian.getI() 306 except Exception,inst: # catching numpy and scipy singular matrix errors 307 if "ingular" in str(inst) and "atrix" in str(inst): 308 # singular matrix error 309 Inverse_Jacobian = numpy.linalg.linalg.pinv(Jacobian) 310 elif "NaN" in str(inst): 311 Inverse_Jacobian = scipy.nan * Jacobian 312 else: 313 # some other error 314 raise 315 self.mca_substrate_control = -Inverse_Jacobian * self.N 316 self.mca_substrate_response = self.mca_substrate_control * self.parameter_elasticities 317 self.mca_flux_control = scipy.eye(len(self.reactions)) - self.substrate_elasticities * Inverse_Jacobian * self.N 318 self.mca_flux_response = self.mca_flux_control * self.parameter_elasticities 319 320 self.mca_coefficients = self.__mca_to_dic__( self.mca_coefficients, self.mca_substrate_control, self.species, self.reactions ) 321 self.mca_coefficients = self.__mca_to_dic__( self.mca_coefficients, self.mca_substrate_response, self.species, self.parameters ) 322 self.mca_coefficients = self.__mca_to_dic__( self.mca_coefficients, self.mca_flux_control, self.reactions, self.reactions ) 323 self.mca_coefficients = self.__mca_to_dic__( self.mca_coefficients, self.mca_flux_response, self.reactions, self.parameters )

324

325 - def toPDF(self):

326 if bool( self.global_options.get_property("write_pdfs") ): 327 import draw_elasticities 328 draw_elasticities.draw_elasticities( self.elasticities_dic, self.filename[:-4]+os.path.sep+"elasticities" ) 329 draw_elasticities.draw_elasticities( self.mca_coefficients, self.filename[:-4]+os.path.sep+"mca" )

330

331 - def solve(self):

332 self.solver = libTI2.solvers.solver.solver(self,signalling=self.signalling) 333 self.species_values = self.solver.steady_state 334 self.steady_state_fluxes = self.solver.steady_state_fluxes 335 if self.is_root_model_instance: 336 self.save_root_results()

337

338 - def dont_solve(self,parentModel):

339 #self.solver = solver.solver(self,just_set_fluxes=True) 340 #self.steady_state_fluxes = self.solver.steady_state_fluxes 341 self.steady_state_fluxes = parentModel.steady_state_fluxes

342

343 - def set_parent_characteristics(self):

344 if self.signalling: 345 self.parent_tau = self.solver.tau 346 self.parent_sigma = self.solver.sigma 347 self.parent_amplitude = self.solver.amplitude 348 self.parent_auc = self.solver.auc

349

350 - def write_ai_to_tempfile(self):

351 self.aih.write_file(self.species_values)

352

353 - def write_ai_to_temptempfile(self):

354 self.aih.write_file(self.species_values,"temp.xml")

355

356 - def mca_ss(self,similarity,parameters_to_be_modified=None):

357 if not parameters_to_be_modified: parameters_to_be_modified = self.parameters 358 self.mca2_coefficients = {} 359 self.mca2_coefficients_normalized = {} 360 for species in self.species+self.reactions: 361 self.mca2_coefficients[species] = {} 362 self.mca2_coefficients_normalized[species] = {} 363 fast_estimated_parameter_number = len([x for x in parameters_to_be_modified if "__i" in x]) 364 for parameter,index in zip(parameters_to_be_modified,range(len(parameters_to_be_modified))): 365 # REMOVE THIS LATER 366 shortparameter = parameter.split("__")[1] 367 if shortparameter != "i": continue 368 # REMOVE UNTIL HERE 369 self.parameter_values[parameter] = self.parameter_values[parameter] + similarity 370 newsolver = solver.solver(self,signalling = self.signalling) 371 new_species_values = newsolver.steady_state 372 self.parameter_values[parameter] = self.parameter_values[parameter] - similarity 373 # WARNING IF CHANGING THE WAY DIFFERENCES ARE RECORDED ALSO CHANGE IN THE TITRATION RECORDS 374 # SEARCH FOR: TITRATION RECORDS 375 for species in self.species: 376 self.mca2_coefficients[species][parameter] = (new_species_values[species] - self.root_steady_state[species]) #/ similarity / self.shift_factor 377 try: 378 self.mca2_coefficients_normalized[species][parameter] = self.mca2_coefficients[species][parameter] * self.parameter_values[parameter] / self.species_values[species] 379 except ZeroDivisionError: 380 self.mca2_coefficients_normalized[species][parameter] = scipy.inf 381 # setting the signal characteristics 382 if self.signalling: 383 self.mca2_coefficients[species]["tau_"+parameter] = newsolver.tau[species] - self.root_tau[species] 384 self.mca2_coefficients[species]["sigma_"+parameter] = newsolver.sigma[species] - self.root_sigma[species] 385 self.mca2_coefficients[species]["amplitude_"+parameter] = newsolver.amplitude[species] - self.root_amplitude[species] 386 self.mca2_coefficients[species]["auc_"+parameter] = newsolver.auc[species] - self.root_auc[species] 387 # evaluate the reaction velocity at this point 388 for reaction in self.reactions: 389 self.mca2_coefficients[reaction][parameter] = (newsolver.steady_state_fluxes[reaction] - self.root_steady_state_fluxes[reaction]) #/ similarity / self.shift_factor 390 # normalize them 391 try: 392 self.mca2_coefficients_normalized[reaction][parameter] = self.mca2_coefficients[reaction][parameter] * self.parameter_values[parameter] / self.steady_state_fluxes[reaction] 393 except ZeroDivisionError: 394 self.mca2_coefficients_normalized[reaction][parameter] = scipy.inf 395 self.global_options.status_bar("simdifparams",index,fast_estimated_parameter_number)

396

397 - def target_identification_for_single_reaction(self,reaction,rre,original_formula,inhibition_formula):

398 """TI for a single reaction. commands are mca.""" 399 ti_new_parameters = [] 400 # get the parameters in the original formula 401 parameters_for_orig_formula = {} 402 compval = rre.pi.compatible_values 403 for sbmlp in reaction.getKineticLaw().getListOfParameters(): 404 sbmlname = sbmlp.getId() 405 val = sbmlp.getValue() 406 try: 407 parameters_for_orig_formula[ compval[ reaction.getId()+"__"+sbmlname ][0] ] = val 408 except KeyError: 409 # this might happen if the kinetic law contains a parameter which is not in the formula 410 pass 411 # replace species names in the inhibition formula 412 my_inhibition_formula = function_manipulation.clean_up_formula(inhibition_formula.formula) 413 for species in compval: 414 if species not in [x.getId() for x in self.sbmldoc.getModel().getListOfSpecies()]: 415 continue 416 my_inhibition_formula = my_inhibition_formula.replace( " "+compval[species][0]+" "," "+species+"_replaced " ) 417 # load data 418 d = libsbml.readSBMLFromString( libsbml.writeSBMLToString(self.sbmldoc) ) 419 if not self.signalling: 420 # write the old steady state to the document in order to accelerate the finding of the steady state 421 self.aih.__write_steady_state_values_to_document__(d,self.species_values) 422 for nsr in d.getModel().getListOfReactions(): 423 if nsr.getId() == reaction.getId(): 424 newsbmlreaction = nsr 425 break 426 # renaming of original parameters 427 for p in newsbmlreaction.getKineticLaw().getListOfParameters(): 428 pid = newsbmlreaction.getId()+"__"+p.getId() 429 p.setId( compval[pid][0] ) 430 # add sbml stuff 431 newparameters = list([(x,inhibition_formula.parameters[x]) for x in inhibition_formula.parameters.keys() if x not in original_formula.parameters]) 432 for key,val in newparameters: 433 p=newsbmlreaction.getKineticLaw().createParameter() 434 p.setId(key) 435 p.setValue(val) 436 # ti_new_parameters.append(newsbmlreaction.getId()+"__"+key) # Testing to remove the real parameters from the list of changeable parameters 437 # create a compartment for the modifiers 438 compartmentname = "TIdes_new_compartment" 439 c=d.getModel().createCompartment() 440 c.setId(compartmentname) 441 c.setVolume(1) 442 # same for the modifiers 443 newmodifiers = list([(x,inhibition_formula.species[1][x]) for x in inhibition_formula.species[1].keys() if x not in original_formula.species[1]]) 444 for key,val in newmodifiers: 445 s=d.getModel().createSpecies() 446 newname = newsbmlreaction.getId()+"__"+key 447 s.setId(newname) 448 s.setInitialAmount(0) 449 s.setCompartment(compartmentname) 450 s.setBoundaryCondition(True) 451 s.setConstant(True) 452 sbmlmod = newsbmlreaction.createModifier( ) 453 sbmlmod.setSpecies(newname) 454 my_inhibition_formula = my_inhibition_formula.replace(" "+key+" "," "+newname+"_replaced ") 455 ti_new_parameters.append(newname) 456 # including the volume of compartments 457 my_inhibition_formula = function_manipulation.clean_up_formula(my_inhibition_formula) 458 compartmentnames = [x.getId() for x in self.sbmldoc.getModel().getListOfCompartments()] 459 for key,lis in compval.iteritems(): 460 if key not in compartmentnames: continue 461 nameinformula = lis[0] 462 my_inhibition_formula = my_inhibition_formula.replace(" "+nameinformula+" "," "+key+"_replaced ") 463 # rename global parameters in the new formula 464 for p in d.getModel().getListOfParameters(): 465 if p.getId() in compval: 466 my_inhibition_formula = my_inhibition_formula.replace(" "+compval[p.getId()][0]+" "," "+p.getId()+"_replaced ") 467 # finally setting the formula 468 # perhaps replace the inhibitor 469 # TODO: what has to be done to include a different inhibitor? ################################################### 470 #if "inhibitor_for_target_identification" in dir(self): 471 # my_inhibition_formula = my_inhibition_formula.replace(" i "," "+self.inhibitor_for_target_identification+"_replaced ") 472 473 474 my_inhibition_formula = function_manipulation.scipy2mathml( my_inhibition_formula ) 475 476 # delete all the replaced marks from the actual formula 477 my_inhibition_formula = my_inhibition_formula.replace("_replaced "," ") 478 479 newsbmlreaction.getKineticLaw().setFormula( my_inhibition_formula ) 480 481 # this looks a bit weird 482 # but perhaps I want to create a new file that should remain in the file system 483 tempdir = "temp"+"".join( random.Random().sample(string.letters+string.digits,12) ) 484 tempfile = tempdir+".xml" 485 w = libsbml.writeSBML(d,tempfile) 486 newModel = Model(tempfile,False,level=self.level-1,reaction_leave_out_list=self.reaction_leave_out_list+[reaction.getId()], signalling = self.signalling, shiftlist=self.shiftlist,is_root_model_instance=False,investigate_reactions=self.investigate_reactions,estimate_inhibitor_concentration=None if not self.estimate_inhibitor_concentration else self.reduction_for_estimation,use_titration=self.use_titration,write_internal_sbml_files = self.write_internal_sbml_files) 487 os.remove(tempfile) 488 os.removedirs(tempdir) 489 self.copy_aim(newModel) 490 self.copy_params(newModel) 491 self.copy_signalling_params(newModel) 492 newModel.after_copy_functions() 493 newModel.dont_solve(self) 494 newModel.mca() 495 newModel.ti_new_parameters = ti_new_parameters 496 return newModel

497 498

499 - def copy_aim(self,newModel):

500 # copy all variables that start with aim_formula to the new Model 501 for var in [x for x in dir(self) if x.startswith("aim_formula")]: 502 newModel.__dict__[var] = self.__dict__[var]

503

504 - def copy_signalling_params(self,newModel):

505 if self.signalling: 506 copyparams = ["parent_tau","parent_sigma","parent_amplitude","parent_auc"] 507 for var in copyparams: 508 newModel.__dict__[var] = self.__dict__[var]

509

510 - def copy_params(self,newModel):

511 """Copies attributes from this model instance to a new Model. 512 513 This should include all important parameters for the target identification. """ 514 copyparams = ["shift_factor","signalling","shiftlist","root_results_filename"] 515 for var in copyparams: 516 newModel.__dict__[var] = self.__dict__[var]

517

518 - def after_copy_functions(self):

519 if not self.is_root_model_instance: 520 self.load_root_results()

521 522

523 - def set_aim_for_target_identification(self,formula):

524 """Currently only substrates and reactions are possible.""" 525 numbers = "1234567890" 526 self.aim_formula = function_manipulation.include_spaces_in_formula( formula ) 527 # identify variables in the formula 528 variables = [] 529 for possible_variable in self.aim_formula.split(): 530 # sort out non-variables 531 if possible_variable[0] in numbers or possible_variable[0] in rre_generator.spacechars: continue 532 if possible_variable not in variables: variables.append(possible_variable) 533 self.aim_formula_variables = variables

534

535 - def set_inhibitor_for_target_identification(self,inhibitor):

536 self.inhibitor_for_target_identification = inhibitor

537

538 - def compute_aim_fulfillment(self,doubledic):

539 returndic = {} 540 for var in self.aim_formula_variables: 541 if var not in doubledic: 542 import sys 543 print "The variable",var,"defined in the AIM formula is not a valid variable." 544 print "Please restrict yourself to the following set:" 545 print "\t"+"\n\t".join( doubledic.keys() ) 546 sys.exit(1) 547 # keys we want to add to the output 548 # it is for adding the tau, sigma, amplitude and auc to the output 549 keys = [] 550 for k in doubledic[self.aim_formula_variables[0]]: 551 keys.append(k) 552 553 # really do something 554 for key in keys: 555 if key not in self.ti_new_parameters: 556 continue 557 vals_for_this_key = {} 558 for var in self.aim_formula_variables: 559 vals_for_this_key[var] = doubledic[var][key] 560 try: 561 returndic[key] = function_manipulation.evaluate_general_formula( self.aim_formula, vals_for_this_key ) 562 except ZeroDivisionError: 563 returndic[key] = scipy.inf 564 return returndic

565

566 - def __identify_reaction__(self,safe_reaction,asTuple=False):

567 """Function was included in order to test parallel programming on the reaction identification part.""" 568 self.__make_sbmldoc_available__() 569 # if the reaction is safely pickled, unpickle it 570 try: 571 reaction = safe_reaction.unpickle() 572 except AttributeError: 573 reaction = safe_reaction 574 # identify the reaction 575 rg = rre_generator.rre_generator(self,reaction) 576 rg.find_formula() 577 self.global_options.condprint("Identified "+reaction.getId()+" as "+rg.name,"identified_reaction") 578 if rg.pi.identical: 579 self.global_options.condprint(rg.old_formula+" corresponds to","identified_reaction") 580 self.global_options.condprint(rg.corresponding_formula,"identified_reaction") 581 if rg.name == 'unknown': 582 # no reaction fits 583 rg.add_formula_and_artificial_inhibition_to_known_formulas() 584 rid = reaction.getId() 585 rg.prepare_to_pickle() 586 self.__hide_sbmldoc__() 587 if asTuple: 588 return (rid,(rg,self.pdf)) 589 else: 590 return rg

591

592 - def compute_signal_characteristics(self,doubledic):

593 reordered_doubledic = {} 594 for key1,single in doubledic.iteritems(): 595 for key2,val in single.iteritems(): 596 for prepend in ["tau","sigma","amplitude","auc"]: 597 if key2.startswith(prepend+"_"): 598 if key2 not in reordered_doubledic: 599 reordered_doubledic[key2] = {} 600 reordered_doubledic[key2][key1] = val 601 aim_fulfillment = {} 602 for parameter in reordered_doubledic: 603 local_parameters = {} 604 for species in reordered_doubledic[parameter]: 605 local_parameters[species] = reordered_doubledic[parameter][species] 606 aim_fulfillment[parameter] = function_manipulation.evaluate_general_formula( self.aim_formula, local_parameters ) 607 return aim_fulfillment

608

609 - def __get_copasi_reactions_description__(self,newModel):

610 """Should be called with the newest Model instance as a parameter from the unmodified model.""" 611 # get the inhibitor species 612 inhibitor_name = newModel.ti_new_parameters[0] 613 # set up the reaction variable 614 reactions = [] 615 for r in newModel.reactions: 616 [s2s,m2s,p2s] = self.reaction2substrate2stoich[r] 617 m2s = copy.deepcopy(m2s) 618 formula = copy.deepcopy(function_manipulation.include_spaces_in_formula(newModel.reaction_equations[r])) 619 # add the inhibitor and possible other modifiers as a modifier if it appears in the formula 620 for modifier in newModel.species+[inhibitor_name]: 621 if " "+modifier+" " in formula: 622 m2s[modifier] = 1 623 for p,v in newModel.parameter_values.iteritems(): 624 # forget the inhibitor which does in copasi account as a metabolite 625 if p == inhibitor_name: continue 626 formula = formula.replace(" "+p+" "," "+str(v)+" ") 627 reactions.append( [ r, s2s, p2s, m2s, formula, True ] ) 628 629 # add the inhibitor as a species into the copasi file 630 copasi_species_values = copy.deepcopy( newModel.species_values ) 631 copasi_species_values[ inhibitor_name ] = 0 632 633 copasi_species = copy.deepcopy(newModel.species)+[inhibitor_name] 634 635 return [copasi_species,copasi_species_values,reactions,inhibitor_name]

636

637 - def __copy_global_options_module_to_global_namespace__(self):

638 if self.use_pp: 639 # the borg pattern global options is not useful for parallelpython things 640 import global_options 641 go = libTI2.global_options.global_options() 642 go.initialize_from_object(self.global_options)

643 644

645 - def __estimate_inhibitor_concentration__(self,newModel,inhibition_type):

646 import copasisolver,copy,function_manipulation 647 648 self.__copy_global_options_module_to_global_namespace__() 649 650 [copasi_species,copasi_species_values,reactions,inhibitor_name] = self.__get_copasi_reactions_description__(newModel) 651 652 # get the correct value for the aim 653 current_val_dict = {} 654 for d in [self.species_values,self.steady_state_fluxes]: 655 for key in d: 656 current_val_dict[key] = d[key] 657 val = eval(self.aim_formula,current_val_dict) 658 aim_value = self.reduction_for_estimation * val 659 quality_factor = 0.001 660 max_quality = quality_factor * abs(1.0-self.reduction_for_estimation) 661 662 if self.use_titration: 663 if self.global_options.get_property("solvertype") == "copasi": 664 ce = copasisolver.copasititrator(copasi_species,copasi_species_values,reactions,False) 665 elif self.global_options.get_property("solvertype") == "soslib": 666 import soslib 667 ce = soslib.SOStitrator(copasi_species,copasi_species_values,reactions,False) 668 elif self.global_options.get_property("solvertype") in bestfortransolver.knownsolvers_and_choice: 669 solvername = self.global_options.get_property("solvertype") 670 solvermodule = __import__(solvername+"solver") 671 titrationconstructor = eval("solvermodule."+solvername+"titrator") 672 ce = titrationconstructor(copasi_species,copasi_species_values,reactions,False) 673 else: 674 print "Titration only available for solvers:" 675 print bestfortransolver.knownsolvers_and_choice 676 print "not for",self.global_options.get_property("solvertype") 677 import sys 678 sys.exit(1) 679 else: 680 ce = copasisolver.copasiestimator(copasi_species,copasi_species_values,reactions,False) 681 result_of_estimation = ce.estimate_inhibitor_concentration(self.aim_formula.strip(),aim_value,inhibitor_name,max_quality) 682 683 # add this to the sdh, if a solution is found 684 newsdh = libTI2.results.simulation_data_handler.simulation_data_handler(None) # this one is constructed explicitly without a model since it has to be pickled 685 686 if abs(result_of_estimation)>1e-10: 687 newsdh.add_data(inhibitor_name,"estimation_"+inhibition_type,0,self.reduction_for_estimation,result_of_estimation) 688 689 return newsdh

690 691

692 - def __simulate_for_shift__(self,newModel,shift,inhibition_type,of,number):

693 self.__make_sbmldoc_available__() 694 newModel.__make_sbmldoc_available__() 695 self.global_options.status_bar("simdifshifts",of,number) 696 697 self.__copy_global_options_module_to_global_namespace__() 698 699 newsdh = libTI2.results.simulation_data_handler.simulation_data_handler(None) # this one is constructed explicitly without a model since it has to be pickled 700 if not self.estimate_inhibitor_concentration: 701 # if the titration is used 702 if self.use_titration: 703 import copasisolver 704 # get species,species_values,and reactions instances 705 [copasi_species,copasi_species_values,reactions,inhibitor_name] = self.__get_copasi_reactions_description__(newModel) 706 if self.global_options.get_property("solvertype") == "copasi": 707 ct = copasisolver.copasititrator(copasi_species,copasi_species_values,reactions,self.signalling,inhibition_type=inhibition_type) 708 elif self.global_options.get_property("solvertype") == "soslib": 709 import soslib 710 ct = soslib.SOStitrator(copasi_species,copasi_species_values,reactions,self.signalling,inhibition_type=inhibition_type) 711 elif self.global_options.get_property("solvertype") in bestfortransolver.knownsolvers_and_choice: 712 solvername = self.global_options.get_property("solvertype") 713 solvermodule = __import__(solvername+"solver") 714 titrationconstructor = eval("solvermodule."+solvername+"titrator") 715 ct = titrationconstructor(copasi_species,copasi_species_values,reactions,self.signalling,inhibition_type=inhibition_type) 716 else: 717 print "Titration only available for solvers:" 718 print bestfortransolver.knownsolvers_and_choice 719 print "not for",self.global_options.get_property("solvertype") 720 import sys 721 sys.exit(1) 722 ct.get_titration_data(inhibitor_name,self.aim_formula.strip(),shift) 723 # TITRATION RECORDS 724 # modify them to contain the difference rather than the absolute values 725 new_values = {} 726 for dic in [newModel.root_steady_state,newModel.root_steady_state_fluxes]: 727 for key in dic: 728 new_values[key] = dic[key] 729 root_aim_value = eval(self.aim_formula,new_values) 730 # substract all the values from the simulation data handler 731 for sd in ct.sdh.data: 732 sd.set_aim_sensitivity( sd.get_aim_sensitivity() - root_aim_value ) 733 newsdh.integrate_sdh( ct.sdh ) 734 else: 735 # do the MCA analysis on this model 736 newModel.mca_ss(shift,newModel.ti_new_parameters) 737 af = newModel.compute_aim_fulfillment(newModel.mca2_coefficients) 738 for p in newModel.ti_new_parameters: 739 if p in af: 740 newsdh.add_data(p,inhibition_type,0,shift,af[p]) 741 742 if self.signalling: 743 af2 = self.compute_signal_characteristics(newModel.mca2_coefficients) 744 for p in af2: 745 newsdh.add_data( p, inhibition_type, 0, shift, af2[p] ) 746 747 if newModel.level > 0 or self.write_internal_sbml_files: 748 # in this if statement the temporary sbml file is constructed and written to disk 749 possible_new_filename = self.path_to_localfiles+os.sep+newModel.ti_new_parameters[0]+"."+str(shift).replace(".","_")+"."+inhibition_type+".xml" 750 # changing the parameters in the libsbml document 751 parameterids = [x.getId() for x in newModel.sbmldoc.getModel().getListOfParameters()] 752 speciesids = [x.getId() for x in newModel.sbmldoc.getModel().getListOfSpecies()] 753 for p in newModel.ti_new_parameters: 754 if p in parameterids: 755 sbmlelem = newModel.sbmldoc.getModel().getParameter(p) 756 getfunc = sbmlelem.getValue 757 setfunc = sbmlelem.setValue 758 elif p in speciesids: 759 sbmlelem = newModel.sbmldoc.getModel().getSpecies(p) 760 if sbmlelem.isSetInitialAmount(): 761 getfunc = sbmlelem.getInitialAmount 762 setfunc = sbmlelem.setInitialAmount 763 else: 764 getfunc = sbmlelem.getInitialConcentration 765 setfunc = sbmlelem.getInitialConcentration 766 setfunc( getfunc() + shift ) 767 # writing out the file 768 libsbml.writeSBML(newModel.sbmldoc,possible_new_filename) 769 # and changing them back 770 for p in newModel.ti_new_parameters: 771 if p in parameterids: 772 sbmlelem = newModel.sbmldoc.getModel().getParameter(p) 773 getfunc = sbmlelem.getValue 774 setfunc = sbmlelem.setValue 775 elif p in speciesids: 776 sbmlelem = newModel.sbmldoc.getModel().getSpecies(p) 777 if sbmlelem.isSetInitialAmount(): 778 getfunc = sbmlelem.getInitialAmount 779 setfunc = sbmlelem.setInitialAmount 780 else: 781 getfunc = sbmlelem.getInitialConcentration 782 setfunc = sbmlelem.getInitialConcentration 783 setfunc( getfunc() - shift ) 784 if newModel.level > 0: 785 # the level information and the reaction_leave_out_list are already updated during the creation of newModel 786 levelBelowModel = Model(possible_new_filename,True,newModel.level,newModel.reaction_leave_out_list,signalling = self.signalling, shiftlist = self.shiftlist, is_root_model_instance=False, investigate_reactions=self.investigate_reactions,estimate_inhibitor_concentration=None if not self.estimate_inhibitor_concentration else self.reduction_for_estimation,use_titration=self.use_titration,write_internal_sbml_files=self.write_internal_sbml_files) 787 levelBelowModel.multiply_shift_factor(shift) 788 levelBelowModel.dont_solve(newModel) 789 newModel.copy_params(levelBelowModel) 790 newModel.copy_signalling_params(levelBelowModel) 791 levelBelowModel.after_copy_functions() 792 #levelBelowModel.write_ai_to_tempfile() # do not save the steady state 793 levelBelowModel.set_aim_for_target_identification(self.aim_formula) 794 levelBelowModel.target_identification() 795 levelBelowModel.write_ti_output() 796 levelBelowModel.drop_sdh() 797 798 799 self.__hide_sbmldoc__ () 800 newModel.__hide_sbmldoc__() 801 return newsdh

802

803 - def target_identification(self):

804 import rre_generator,sboreader,pre_defined_formulas,sbmlpickle,datetime 805 """Runs through all reactions, replaces them with inhibitory kinetics.""" 806 self.pdf = pre_defined_formulas.pre_defined_formulas() 807 self.sdh = libTI2.results.simulation_data_handler.simulation_data_handler(self.filename,reload_file=self.path_to_localfiles+os.sep+"sdh") 808 rglist = [] 809 810 811 if self.use_pp: 812 import pp 813 # get pp servers 814 servers = [] 815 f = file("servers","r") 816 for l in f: 817 if l.startswith("#"): continue 818 servers.append(l.split(",")[1]) 819 f.close() 820 # read secret file 821 try: 822 f = open("secret","r") 823 secret = f.readlines()[0] 824 f.close() 825 except: 826 import sys 827 print "Please create a file calles \"secret\" with one long string in the first line." 828 sys.exit(1) 829 830 print "starting servers",servers 831 job_server = pp.Server(ppservers=tuple(servers),secret=secret) 832 #job_server = pp.Server() 833 safe_reactions = [sbmlpickle.sbmlpickle(reaction) for reaction in self.sbmldoc.getModel().getListOfReactions()] 834 self.__hide_sbmldoc__() 835 for reaction in safe_reactions: 836 # REMOVED UNTIL FIXED 837 #if self.use_pp: 838 # rglist.append( job_server.submit( self.__identify_reaction__ , ( reaction, True ),globals=globals() ) ) 839 #else: 840 rglist.append(self.__identify_reaction__(reaction,True)) 841 self.__make_sbmldoc_available__() 842 # REMOVE UNTIL FIXED 843 #if self.use_pp: 844 # rglist = [job() for job in rglist] 845 rid2rg=dict(rglist) 846 rglist = [] 847 848 # modify the shiftlist for the titration experiment 849 if self.use_titration: 850 shiftlist = [(max(self.shiftlist),0)] if not self.estimate_inhibitor_concentration else [] 851 else: 852 shiftlist = zip(self.shiftlist,range(len(self.shiftlist))) 853 854 len_shiftlist = len(shiftlist) 855 # in the next step we want to consider only reaction with a higher reaction number than the ones already tried out 856 minimum_reaction_number = max([self.reaction_number[rid] for rid in self.reaction_leave_out_list]+[-1]) 857 for rid in [r.getId() for r in self.sbmldoc.getModel().getListOfReactions()]: 858 ###if rid in self.reaction_leave_out_list: continue # leave out this reaction 859 if self.reaction_number[rid] <= minimum_reaction_number: continue # leave out reactions which have a lower id than the ones already tried out 860 if len(self.investigate_reactions)>0: 861 if rid not in self.investigate_reactions: 862 continue 863 (rg,pdf) = rid2rg[rid] 864 self.pdf.integrate(pdf) 865 for inhibition_type in pre_defined_formulas.inhibition_types: 866 inhibition_formula = self.pdf.get_inhibition_formula_for_name(rg.name,inhibition_type) 867 # do we know an inhibition 868 if not inhibition_formula: continue 869 # ok, we know one 870 self.global_options.condprint("Found "+inhibition_type+" inhibition for "+rid,"identified_reaction") 871 original_formula = self.pdf.find_formula_with_name(rg.name) 872 reaction = self.sbmldoc.getModel().getReaction(rid) 873 newModel = self.target_identification_for_single_reaction(reaction,rg,original_formula,inhibition_formula) 874 if not self.estimate_inhibitor_concentration: 875 af = newModel.compute_aim_fulfillment(newModel.mca_coefficients) 876 for p in newModel.ti_new_parameters: 877 if p in af: 878 self.sdh.add_data(p,inhibition_type,0.0,0.0,af[p]) 879 self.__hide_sbmldoc__() 880 newModel.__hide_sbmldoc__() 881 if self.estimate_inhibitor_concentration: 882 if self.use_pp: 883 rglist.append( job_server.submit( self.__estimate_inhibitor_concentration__, (newModel,inhibition_type), modules=("global_options",),globals=globals() ) ) 884 else: 885 rglist.append( self.__estimate_inhibitor_concentration__(newModel,inhibition_type) ) 886 for shift,index in shiftlist: 887 if self.use_pp: 888 rglist.append( job_server.submit( self.__simulate_for_shift__, (newModel,shift,inhibition_type,index,len_shiftlist), modules=("global_options",),globals=globals() ) ) 889 else: 890 rglist.append( self.__simulate_for_shift__(newModel,shift,inhibition_type,index,len_shiftlist) ) 891 newModel.__make_sbmldoc_available__() 892 self.__make_sbmldoc_available__() 893 if self.use_pp: 894 newrglist = [] 895 rglen = len(rglist) 896 index = 0 897 for job in rglist: 898 index = index + 1 899 newrglist.append( job() ) 900 self.global_options.status_bar("ppjobid",index,rglen) 901 job_server.print_stats() 902 rglist = newrglist 903 for sdh in rglist: 904 self.sdh.integrate_sdh(sdh)

905 906 907 908

909 - def write_ti_output(self):

910 # repeat writing the pdfs in case of an error 911 try: 912 if self.global_options.get_property("write_pdfs"): 913 self.sdh.write_pdfs(self.path_to_localfiles+os.sep+"global.pdf") 914 except IOError,inst: 915 print inst 916 print inst.args 917 print "Problem with the pdf writing." 918 self.global_options.unset_property("write_pdfs") 919 try: 920 if not self.signalling: 921 self.sdh.write_svg(self.path_to_localfiles+os.sep+"global.svg") 922 except IOError: 923 print "Problem with the svg writing" 924 except (ValueError, TypeError): 925 print "Problem with the svg writing. perhaps a strange value in the SDH."

926

927 - def drop_sdh(self):

928 self.sdh.drop(self.path_to_localfiles+os.sep+"sdh")

929

930 - def load_sdh(self):

931 self.sdh = libTI2.results.simulation_data_handler.load_sdh(self.path_to_localfiles+os.sep+"sdh")

932

933 - def multiply_shift_factor(self,val):

934 self.shift_factor = self.shift_factor * val

935

936 - def save_root_results(self):

937 root_res = libTI2.results.simulation_data_handler.root_results(self.signalling) 938 root_res.save(self) 939 self.root_results_filename = root_res.filename

940

941 - def load_root_results(self):

942 root_res = libTI2.results.simulation_data_handler.root_results(self.signalling) 943 root_res.load(self,self.root_results_filename)

944

945 - def __init__(self,filename=None,reload=True,level=1,reaction_leave_out_list=[],signalling=False,shiftlist=[0.1],is_root_model_instance=True,investigate_reactions=[],estimate_inhibitor_concentration=None,use_titration=False,write_internal_sbml_files=False):

946 """Filename can also be a libsbml.Document().""" 947 self.level = level 948 self.reaction_leave_out_list = reaction_leave_out_list 949 self.global_options = libTI2.global_options.global_options() 950 self.use_pp = bool( self.global_options.get_property( "use_parallel" ) ) 951 self.reload_cached_file = reload 952 self.shift_factor = 1.0 953 self.use_titration = use_titration 954 self.signalling = signalling 955 self.shiftlist = shiftlist 956 self.is_root_model_instance = is_root_model_instance 957 self.investigate_reactions = investigate_reactions 958 self.estimate_inhibitor_concentration = (estimate_inhibitor_concentration != None) 959 self.write_internal_sbml_files = write_internal_sbml_files 960 if self.estimate_inhibitor_concentration: 961 self.reduction_for_estimation = float(estimate_inhibitor_concentration) 962 # forcing other variables 963 self.level = 1 964 self.shiftlist = [] 965 # if the filename is a string containing xml then convert it to a libsbml document 966 if type(filename) == type("") and filename.startswith( "<?xml" ): 967 filename = libsbml.readSBMLFromString(filename) 968 if not filename: 969 self.sbmldoc = libsbml.SBMLDocument() 970 self.filename = "test.xml" 971 self.sbmldoc.createModel() 972 elif type(filename) == type(""): 973 self.sbmldoc = self.load_sbml(filename) 974 else: 975 # filename is a libsbml instance or something 976 self.sbmldoc = filename.getSBMLDocument() 977 self.filename = False 978 self.small_value = 10e-6 979 self.load_differential_equations()

Source Code for Module libTI2.model.model