Source code for bionetgen.atomizer.sbml2bngl

# -*- coding: utf-8 -*-
"""
Created on Tue Dec  6 17:42:31 2011

@author: proto
"""
from copy import deepcopy, copy
from bionetgen.atomizer.writer import bnglWriter as writer

log = {"species": [], "reactions": []}
import re, sympy
from collections import Counter
from collections import defaultdict
import math as pymath
from bionetgen.atomizer.utils.util import logMess, TranslationException
import libsbml
from bionetgen.atomizer.bngModel import bngModel

import sympy
from sympy import Function
from sympy.abc import _clash
from sympy.printing.str import StrPrinter
from sympy.core.sympify import SympifyError

# Define 2 and 3 argument functions
# for sympy parsing
[docs]class sympyPiece(Function): nargs = (3, 4, 5)
[docs]class sympyIF(Function): nargs = 3
[docs]class sympyGT(Function): nargs = 2
[docs]class sympyLT(Function): nargs = 2
[docs]class sympyGEQ(Function): nargs = 2
[docs]class sympyLEQ(Function): nargs = 2
[docs]class sympyAnd(Function): nargs = (2, 3, 4, 5)
[docs]class sympyOr(Function): nargs = (2, 3, 4, 5)
[docs]class sympyNot(Function): nargs = 1
[docs]def factorial(x): temp = x acc = 1 while temp > 0: acc *= temp temp -= 1 return acc
[docs]def comb(x, y, exact=True): return factorial(x) / (factorial(y) * factorial(x - y))
bioqual = [ "BQB_IS", "BQB_HAS_PART", "BQB_IS_PART_OF", "BQB_IS_VERSION_OF", "BQB_HAS_VERSION", "BQB_IS_HOMOLOG_TO", "BQB_IS_DESCRIBED_BY", "BQB_IS_ENCODED_BY", "BQB_ENCODES", "BQB_OCCURS_IN", "BQB_HAS_PROPERTY", "BQB_IS_PROPERTY_OF", "BQB_HAS_TAXON", "BQB_UNKNOWN", ] modqual = [ "BQM_IS", "BQM_IS_DESCRIBED_BY", "BQM_IS_DERIVED_FROM", "BQM_IS_INSTANCE_OF", "BQM_HAS_INSTANCE", "BQM_UNKNOWN", ] annotationHeader = {"BQB": "bqbiol", "BQM": "bmbiol"}
[docs]def unrollSBMLFunction(function, sbmlFunctions): """ remove calls to functions inside functions """ modificationFlag = True recursionIndex = 0 # remove calls to other sbml functions while modificationFlag and recursionIndex < 20: modificationFlag = False for sbml in sbmlFunctions: if sbml in function: temp = writer.extendFunction(function, sbml, sbmlFunctions[sbml]) if temp != function: function = temp modificationFlag = True recursionIndex += 1 break return function
[docs]class SBML2BNGL: """ contains methods for extracting and formatting those sbml elements that are translatable into bngl """ def __init__(self, model, useID=True, replaceLocParams=True, obs_map_file=None): self.bngModel = bngModel() self.bngModel.useID = useID self.bngModel.replaceLocParams = replaceLocParams self.bngModel.obs_map_file = obs_map_file self.useID = useID self.replaceLocParams = replaceLocParams self.model = model self.tags = {} self.speciesUnits = {} self.isConversion = True self.boundaryConditionVariables = [] self.speciesDictionary = {} self.speciesMemory = [] self.speciesAnnotationDict = None self.reactionDictionary = {} self.speciesAnnotation = None self.speciesCompartments = None self.unitDefinitions = self.getUnitDefinitions() self.convertSubstanceUnits = False # keep track of obs names self.obs_names = [] self.obs_map = {} self.param_repl = {} # ASS - I think there should be a check for compartments right here # to determine if a) any compartment is actually used and # b) remove unused ones and c) if only one is used with volume # 1.0, remove that and write a normal model without compartments # The following will give access to all methods to this information self.noCompartment = None # multi compartment warning flag self.multi_comp_warn = False self.check_noCompartment() self.getSpecies() # We want to re-use all the symbols # generated by Sympy and then access it later on self.all_syms = {} self.all_syms.update(_clash) self.all_syms.update({"pow": pow}) # self.all_syms["__epsilon__"] = sympy.symbols("__epsilon__") # Do we want this? we'll have to handle it downstream # handling a lot of functions via a dummy function # that can take 1,2,3 variables self.all_syms["Piecewise"] = sympyPiece self.all_syms["piecewise"] = sympyPiece self.all_syms["gt"] = sympyGT self.all_syms["lt"] = sympyLT self.all_syms["if"] = sympyIF self.all_syms["sympyAnd"] = sympyAnd self.all_syms["sympyOr"] = sympyOr self.all_syms["sympyNot"] = sympyNot # We are trying to replace things that we know # are only in assignment rules in functions self.only_assignment_dict = {} if not hasattr(self, "used_molecules"): self.used_molecules = [] if not hasattr(self, "used_symbols"): self.used_symbols = [] if not hasattr(self, "arule_map"): self.arule_map = {} # Only write epsilon if we must self.write_epsilon = False
[docs] def setConversion(self, conversion): self.isConversion = conversion
[docs] def reset(self): self.tags = {} self.boundaryConditionVariables = [] self.speciesDictionary = {} self.speciesMemory = [] self.getSpecies() self.reactionDictionary = {}
[docs] def static_var(varname, value): def decorate(func): setattr(func, varname, value) return func return decorate
[docs] def extractModelAnnotation(self): metaInformation = {} annotation = self.model.getAnnotation() lista = libsbml.CVTermList() libsbml.RDFAnnotationParser.parseRDFAnnotation(annotation, lista) for idx in range(lista.getSize()): # biol,qual = lista.get(idx).getBiologicalQualifierType(), lista.get(idx).getModelQualifierType() qualifierType = lista.get(idx).getQualifierType() qualifierDescription = ( bioqual[lista.get(idx).getBiologicalQualifierType()] if qualifierType else modqual[lista.get(idx).getModelQualifierType()] ) if qualifierDescription not in metaInformation: metaInformation[qualifierDescription] = set([]) for idx2 in range(0, lista.get(idx).getResources().getLength()): resource = lista.get(idx).getResources().getValue(idx2) metaInformation[qualifierDescription].add(resource) return metaInformation
[docs] def getMetaInformation(self, additionalNotes): # get unit information unitList = self.getUnitDefinitions() metaInformation = self.extractModelAnnotation() modelHistory = self.model.getModelHistory() if modelHistory: try: tmp = libsbml.ModelHistory.getCreator( self.model.getModelHistory(), 0 ).getFamilyName() tmp += ( " " + libsbml.ModelHistory.getCreator( self.model.getModelHistory(), 0 ).getGivenName() ) metaInformation["creatorEmail"] = ( "'" + libsbml.ModelHistory.getCreator( self.model.getModelHistory(), 0 ).getEmail() + "'" ) metaInformation["creatorName"] = "'" + tmp + "'" except: metaInformation["creatorEmail"] = "''" metaInformation["creatorName"] = "''" metaInformation.update(additionalNotes) metaString = "###\n" for element in metaInformation: if type(metaInformation[element]) == set: metaInformation[element] = list(metaInformation[element]) metaString += "#@{0}:{1}\n".format(element, (metaInformation[element])) metaString += "###\n" # adding the string to model self.bngModel.metaString = metaString return metaString
[docs] def isSameNameDifferentCompartment(self, name): speciesList = [] for species in self.model.getListOfSpecies(): if species.getName() == name: speciesList.append(species.getCompartment()) return len(speciesList) == len(set(speciesList))
[docs] def getRawSpecies(self, species, parameters=[], logEntries=True): """ *species* is the element whose SBML information we will extract this method gets information directly from an SBML related to a particular species. It returns id,initialConcentration,(bool)isconstant and isboundary, and the compartment It also accounts for the fact that sometimes ppl use the same name for molecules with different identifiers """ identifier = species.getId() name = species.getName() if name == "": name = identifier if species.isSetInitialConcentration(): initialValue = species.getInitialConcentration() # if we have compartments and a set volume # we need to convert this to species counts # species_comp = species.compartment # compVol = self.compartmentDict[species_comp] # # TODO: Handle unit conversions here if units are given # initialValue *= float(compVol) # if not species.getHasOnlySubstanceUnits(): # initialValue *= 6.022e23 else: initialValue = species.getInitialAmount() # we need to ensure we have concentrations species_comp = species.compartment compVol = self.compartmentDict[species_comp] # TODO: Handle unit conversions here if units are given # if compVol > 0: # initialValue /= float(compVol) # TODO: If using moles of stuff # initialValue *= 6.022e23 isConstant = species.getConstant() isBoundary = species.getBoundaryCondition() # FIXME: this condition means that a variable/species can be changed # by rules and/or events. this means that we effectively need a variable # changed by a function that tracks this value, and all references # to this observable have to be changed to the referrencing variable. # http://sbml.org/Software/libSBML/docs/java-api/org/sbml/libsbml/Species.html if isBoundary and not isConstant: # isConstant = True if ( not species.isSetInitialConcentration() and not species.isSetInitialAmount() ): initialValue = 1 compartment = species.getCompartment() boundaryCondition = species.getBoundaryCondition() standardizedName = standardizeName(name) # if its a species that appends the compartment name remove it, it is not necessary in bionetgen if standardizedName.endswith("_{0}".format(compartment)): standardizedName = standardizedName.replace("_{0}".format(compartment), "") # if standardizedName in ['Source','Trash','Sink']: # standardizedName = '0' if standardizedName in parameters: standardizedName = "sp_{0}".format(standardizedName) # it cannot start with a number if standardizedName[:1].isdigit(): standardizedName = "s" + standardizedName # two species cannot have the same name. Ids are unique but less informative, however typically species can be differentiated # by compartment if logEntries and standardizedName != "0": if standardizedName in self.speciesMemory: if len(list(self.model.getListOfCompartments())) == 1: standardizedName += "_" + species.getId() else: # we can differentiate by compartment tag, no need to attach it to the name # however an actual check needs to be made to make sure that the modeler is actually # changing compartment information. If not, use id to differentiate. if not self.isSameNameDifferentCompartment(species.getName()): standardizedName += "_" + species.getId() self.speciesMemory.append(standardizedName) if boundaryCondition: if standardizedName not in self.boundaryConditionVariables: self.boundaryConditionVariables.append(standardizedName) self.speciesDictionary[identifier] = standardizedName returnID = identifier if self.useID else self.speciesDictionary[identifier] if species.isSetSubstanceUnits(): self.speciesUnits[returnID] = species.getSubstanceUnits() values = {} values["returnID"] = returnID if species.isSetInitialConcentration(): values["initialConcentration"] = initialValue values["initialAmount"] = -1 elif species.isSetInitialAmount(): values["initialAmount"] = initialValue values["initialConcentration"] = -1 else: values["initialAmount"] = -1 values["initialConcentration"] = -1 values["isConstant"] = isConstant values["isBoundary"] = isBoundary values["compartment"] = compartment values["name"] = name values["identifier"] = identifier if len(species.conversion_factor) > 0: # we have a conversion factor to account for values["conversionFactor"] = species.conversion_factor else: values["conversionFactor"] = None return values
[docs] def getPrunnedTree(self, math, remainderPatterns, artificialObservables={}): """ walks through a series of * nodes and removes the remainder reactant factors arg:remainderPatterns: argumetns to be removed from the tree it also changes references to time variables to the keyword 'Time' artificialObservables: species that are changed through an sbml assignment rule. their usage in bng requires special handling. """ swapDict = {libsbml.AST_NAME_TIME: "Time"} for node in [ x for x in (math.getLeftChild(), math.getRightChild()) if x != None ]: if node.getType() in list(swapDict.keys()): node.setName(swapDict[node.getType()]) if math.getCharacter() == "+" and math.getNumChildren() == 1: math = math.getLeftChild() while (math.getCharacter() == "*" or math.getCharacter() == "/") and len( remainderPatterns ) > 0: if libsbml.formulaToString(math.getLeftChild()) in remainderPatterns: remainderPatterns.remove(libsbml.formulaToString(math.getLeftChild())) if math.getCharacter() == "*": math = math.getRightChild() else: math.getLeftChild().setValue(1) elif libsbml.formulaToString(math.getRightChild()) in remainderPatterns: remainderPatterns.remove(libsbml.formulaToString(math.getRightChild())) math = math.getLeftChild() else: if (math.getLeftChild().getCharacter()) == "*": math.replaceChild( 0, self.getPrunnedTree(math.getLeftChild(), remainderPatterns) ) if (math.getRightChild().getCharacter()) == "*": math.replaceChild( math.getNumChildren() - 1, self.getPrunnedTree(math.getRightChild(), remainderPatterns), ) break return math
[docs] def getIsTreeNegative(self, math): """ walks through a series of * nodes and detects whether there's a negative factor fixme: we should actually test if the number of negative factors is odd right now we are relying on the modelers not being malicious when writing their rates laws. """ if math.getCharacter() == "*" or math.getCharacter() == "/": if math.getLeftChild().isUMinus(): math.setCharacter("+") # math.getLeftChild().setValue(long(libsbml.formulaToString(math.getLeftChild())[1:])) return True elif math.getLeftChild().getNumChildren() == 0 and libsbml.formulaToString( math.getLeftChild() ).startswith("-"): math.getLeftChild().setValue( int(libsbml.formulaToString(math.getLeftChild())[1:]) ) return True elif math.getRightChild().isUMinus(): math.setCharacter("+") # math.getRightChild().setValue(long(libsbml.formulaToString(math.getRightChild())[1:])) return True elif math.getRightChild().getNumChildren() == 0 and libsbml.formulaToString( math.getRightChild() ).startswith("-"): math.getRightChild().setValue( int(libsbml.formulaToString(math.getRightChild())[1:]) ) return True else: if (math.getLeftChild().getCharacter()) in ["*", "/", "-"]: if self.getIsTreeNegative(math.getLeftChild()): return True if (math.getRightChild().getCharacter()) in ["*", "/", "-"]: if self.getIsTreeNegative(math.getRightChild()): return True elif math.getCharacter() == "-" and math.getNumChildren() == 1: math.setCharacter("+") return True return False
[docs] def getUnitDefinitions(self): self.unitDictionary = {} for unitDefinition in self.model.getListOfUnitDefinitions(): unitList = [] for unit in unitDefinition.getListOfUnits(): correctedUnit = libsbml.Unit.convertToSI(unit) # unitList.append({'kind':unit.getKind(), 'scale':unit.getScale(),'multiplier':unit.getMultiplier(), 'exponent': unit.getExponent(), 'name':name}) for unit2 in correctedUnit.getListOfUnits(): name = ( unit.getName() if unit.getName() else unitDefinition.getName() ) unitList.append( { "kind": unit2.getKind(), "scale": unit2.getScale(), "multiplier": unit2.getMultiplier(), "exponent": unit2.getExponent(), "name": name, } ) self.unitDictionary[unitDefinition.getId()] = unitList return self.unitDictionary
[docs] def preProcessStoichiometry(self, reactants): """ checks for reactants with the same name in the reactant list. This is mainly to account for reactants that dont have the stoichiometry flag properly set and instead appear repeated """ uniqueReactantDict = defaultdict(int) for reactant in reactants: uniqueReactantDict[reactant[0]] += reactant[1] return [(x, uniqueReactantDict[x]) for x in uniqueReactantDict]
[docs] def removeFactorFromMath(self, math, reactants, products, artificialObservables): """ it also adds symmetry factors. this checks for symmetry in the species names artificialObservables: species names that are changed through assignment rules. their use requires special care when calculating a rate """ ifStack = Counter() remainderPatterns = [] highStoichoiMetryFactor = 1 processedReactants = self.preProcessStoichiometry(reactants) # ASS: I'm doing a hack, this is a flag to indicate # that a species appears on both sides of a reaction bothSides = False for x in processedReactants: # this is the symmtery factor for the rate constant highStoichoiMetryFactor *= factorial(x[1]) y = [i[1] for i in products if i[0] == x[0]] if len(y) > 0: bothSides = True y = y[0] if len(y) > 0 else 0 # TODO: check if this actually keeps the correct dynamics # this is basically there to address the case where theres more products # than reactants (synthesis) if x[1] > y: highStoichoiMetryFactor /= comb(int(x[1]), int(y), exact=True) # print("HSMF comb: {}".format(highStoichoiMetryFactor)) for counter in range(0, int(x[1])): remainderPatterns.append(x[0]) # for x in products: # highStoichoiMetryFactor /= math.factorial(x[1]) # remainderPatterns = [x[0] for x in reactants] # print remainderPatterns, artificialObservables.keys() math = self.getPrunnedTree(math, remainderPatterns) rateR = libsbml.formulaToString(math) for element in remainderPatterns: ifStack.update([element]) # This is here to denote if we are removing an extra factor # of one of the elements (on top of the one that might or might # not be there) so we have our element in the denominator # hence the special treatment. for element in ifStack: if ifStack[element] > 1: # ASS - removing if statement from functional rate definitions # rateR = '{1}/({0}^{2} + __epsilon__)'.format(element, rateR, ifStack[element]) rateR = "{1}/({0}^{2})".format(element, rateR, ifStack[element]) self.write_epsilon = True else: # ASS - removing if statement from functional rate definitions # rateR = '{1}/({0} + __epsilon__)'.format(element, rateR) rateR = "{1}/({0})".format(element, rateR) # self.write_epsilon = True numFactors = max(math.getNumChildren(), len(ifStack)) if pymath.isinf(highStoichoiMetryFactor): rateR = "{0} * 1e20".format(rateR) logMess( "ERROR:SIM204", 'Found usage of "inf" inside function {0}'.format(rateR) ) elif highStoichoiMetryFactor != 1 and bothSides: # ASS # there is something wrong here, this multiplies regular # rate constant by the highStoichiometry value and it's simply # incorrect. # Update: I think not multiplying is correct for most cases, # I think this changes when we have a species on both sides # of the reaction. Then, and only then, this parsing is relevant # I believe. # Update: hence the "bothSides" flag rateR = "{0} * {1}".format(rateR, int(highStoichoiMetryFactor)) # we are adding a factor to the rate so we need to account for it when # we are constructing the bngl equation (we dont want constrant expressions in there) numFactors = max(1, numFactors) # print("reactants and products: {} \n {} \n".format(reactants, products)) # print("HSMF: {}".format(highStoichoiMetryFactor)) # print("rateR: {}".format(rateR)) # print("numFactors: {}".format(numFactors)) return rateR, numFactors
[docs] def isAmount(self, reactantName): for species in self.model.getListOfSpecies(): if species.getName() == reactantName: if species.isSetInitialAmount(): return True return False
[docs] def calculate_factor(self, react, prod, expr, removed): """ Calculates factors from reactants and products? """ ifStack = Counter() remainderPatterns = [] highStoichoiMetryFactor = 1 processedReactants = self.preProcessStoichiometry(react) # ASS: I'm doing a hack, this is a flag to indicate # that a species appears on both sides of a reaction bothSides = False for x in processedReactants: # this is the symmtery factor for the rate constant highStoichoiMetryFactor *= factorial(x[1]) y = [i[1] for i in prod if i[0] == x[0]] if len(y) > 0: bothSides = True y = y[0] if len(y) > 0 else 0 if x[1] > y: highStoichoiMetryFactor /= comb(int(x[1]), int(y), exact=True) for counter in range(0, int(x[1])): remainderPatterns.append(x[0]) for rem in removed: if rem in remainderPatterns: remainderPatterns.remove(rem) for element in remainderPatterns: ifStack.update([element]) # Replace certain symbols s.t. libsbml.parseFormula gets it expr = expr.replace("**", "^") numFactors = max(libsbml.parseFormula(expr).getNumChildren(), len(ifStack)) if pymath.isinf(highStoichoiMetryFactor): logMess( "ERROR:SIM204", 'Found usage of "inf" inside function. This might be due to non-integer stoichiometry as well.', ) elif highStoichoiMetryFactor != 1 and bothSides: numFactors = max(1, numFactors) return numFactors
[docs] def find_all_symbols(self, math, reactionID): time_warn = False curr_keys = list(self.all_syms.keys()) if math is None: logMess( "ERROR:SIM211", f"Math for reaction with ID '{reactionID}' is not defined", ) raise TranslationException( f"ERROR:SIM211: Math for reaction with ID '{reactionID}' is not defined" ) l = math.getListOfNodes() replace_dict = {} for inode in range(l.getSize()): node = l.get(inode) # Sympy doesn't like "def" in our string name = node.getName() if name == "def": self.all_syms["def"] = sympy.symbols("__DEF__") replace_dict["def"] = "__DEF__" continue if name == "lambda": self.all_syms["lambda"] = sympy.symbols("__LAMBDA__") replace_dict["lambda"] = "__LAMBDA__" continue if name == "as": self.all_syms["as"] = sympy.symbols("__AS__") replace_dict["as"] = "__AS__" continue if name == "del": self.all_syms["del"] = sympy.symbols("__DEL__") replace_dict["del"] = "__DEL__" continue if name == "time": # Need to raise warning but do it # once per reaction if time_warn: pass else: time_warn = True if reactionID is not None: logMess( "WARNING:RATE004", "Kinetic law of reaction {} is using the time function. The time functionality in BNG might not work as expected and throw off the translation.".format( reactionID ), ) else: logMess( "WARNING:RATE004", "At least one reactions kinetic law is using the time function. The time functionality in BNG might not work as expected and throw off the translation.", ) if name is not None: if name not in curr_keys: self.all_syms[name] = sympy.symbols(name) # let's parse the formula and get non-numerical symbols form = libsbml.formulaToString(math) # If we need to replace anything # TODO: Replace all of these with regexp for it in replace_dict.items(): form = form.replace(it[0], it[1]) # Let's also pool this in used_symbols for sym in self.all_syms.keys(): if sym not in self.used_symbols: self.used_symbols.append(sym) # Sympy doesn't allow and/not/or to be used # outside what it deems to be acceptable # TODO: Replace all of these with regexp if "piecewise(" in form: form = form.replace("piecewise(", "sympyPiece(") replace_dict["piecewise"] = "sympyPiece" if "gt(" in form: form = form.replace("gt(", "sympyGT(") replace_dict["gt"] = "sympyGT" if "geq(" in form: form = form.replace("geq(", "sympyGEQ(") replace_dict["geq"] = "sympyGEQ" if "lt(" in form: form = form.replace("lt(", "sympyLT(") replace_dict["lt"] = "sympyLT" if "leq(" in form: form = form.replace("leq(", "sympyLEQ(") replace_dict["leq"] = "sympyLEQ" if "if(" in form: form = form.replace("if(", "sympyIF(") replace_dict["if"] = "sympyIF" if "and(" in form: form = form.replace("and(", "sympyAnd(") replace_dict["and"] = "sympyAnd" # TODO: "or(" catches stuff like "floor(" and other # potential functions. This needs to be extended # to more potential or statements (e.g. *or(, +or( etc # the same goes for other functions too but this is # particularly a problem for this one if " or(" in form: form = form.replace("or(", "sympyOr(") replace_dict["or"] = "sympyOr" if "not(" in form: form = form.replace("not(", "sympyNot(") replace_dict["not"] = "sympyNot" return form, replace_dict
[docs] def analyzeReactionRate( self, math, compartmentList, reversible, rReactant, rProduct, reactionID, parameterFunctions, rModifier=[], sbmlFunctions={}, split_rxn=False, ): """ This functions attempts to obtain the left and right hand sides of a rate reaction function given a MathML tree. It also removes compartments and chemical factors from the function Keyword arguments: math -- the MathML math object compartmentList -- a list of all the compartments in the system reversible -- a boolean indicated whether there's should be a backward rate rReactant -- a string list of the reactants. rProduct -- a string list of the products. sbmlFunctions -- a list of possible nested functiosn that we need to remove """ # We turn the AST object to string first mathstring = libsbml.formulaToString(math) # unroll sbml functions mathstring = unrollSBMLFunction(mathstring, sbmlFunctions) temp = libsbml.parseFormula(mathstring) # if we could parse back after all those modifications... if temp: math = temp # I'm removing the function defined and putting the code in here # directly react = rReactant prod = rProduct rateL = None rateR = None # Let's pull everything in the formula as symbols to use # with sympify form, replace_dict = self.find_all_symbols(math, reactionID) # let's pull all names all_names = [i[0] for i in react] + [i[0] for i in prod] # SymPy is wonderful, _clash1 avoids built-ins like E, I etc # FIXME:can we adjust the assignment rule stuff here? try: sym = sympy.sympify(form, locals=self.all_syms) except SympifyError as e: logMess( "ERROR:SYMP001", "Sympy couldn't parse a function, sorry but we can't handle this function definition.", ) raise e # TranslationException(reactionID) # Remove compartments if we use them. # if not self.noCompartment: compartments_to_remove = [sympy.symbols(comp) for comp in compartmentList] # TODO: This is not fully correct, we need to know what # compartment is on what side which is not currently # being provided to this function for comp in compartments_to_remove: if comp in sym.atoms(): # Further issue, I know that this should be # a multiplication but for BMD2 this is actually a # problem? In fact, it looks like this is the case # for regular mass action in SBML? # This doesn't look right and it is a current # hack? n, d = sym.as_numer_denom() if comp in n.atoms(): sym = sym / comp elif comp in d.atoms(): sym = sym * comp else: pass # If we are splitting, we don't need to do much if split_rxn: rate = str(sym).replace("**", "^") for it in replace_dict.items(): rate = rate.replace(it[1], it[0]) return rate, "", 1, 1, False, split_rxn # expand and take the terms out as left and right exp = sympy.expand(sym) # This shows if we can get X - Y ###### SPLIT RXN ####### # TODO: Figure out if something CAN be mass action # and if not, just skip the rest and use split_rxn ###### SPLIT RXN ####### if exp.is_Add: react_expr, prod_expr = self.gather_terms(exp) # l,r = exp.as_two_terms() # Let's also ensure that we have a + and - term # if str(l).startswith("-") or str(r).startswith("-"): # if str(l).startswith("-"): # fwd_expr = r # back_expr = l # else: # fwd_expr = l # back_expr = r if react_expr is None: # this means we only have negative # values in the kinetic laws. I guess # this is possible if one of the parameters # turns it into + overall rate. We'll turn this into # unidirectional after multiplying by -, also will # raise a warning logMess( "WARNING:RATE005", "All terms in reaction {} are negative. It's possible that one of the parameters make the overall rate positive but ensure that this is what you expect".format( reactionID ), ) react_expr = -1 * prod_expr prod_expr = None if prod_expr is not None: # Also get and parse the symbols react_bols = [x[0] for x in react] prod_bols = [x[0] for x in prod] react_symbols = sympy.symbols(react_bols) prod_symbols = sympy.symbols(prod_bols) # Now we can manipulate it # react_expr = fwd_expr removedL = [] for ibol, bol in enumerate(react_symbols): stoi = int(react[ibol][1]) # Now we can remove it react_expr = react_expr / (bol**stoi) removedL += [str(bol) for i in range(stoi)] # Check if we can get 0 in the denominator add_eps_react = False n, d = react_expr.as_numer_denom() for ibol, bol in enumerate(react_symbols): if bol in d.atoms(): d = d.subs(bol, 0) if d == 0: # logMess('WARNING:RATE001', 'Denominator of rate constant in reaction {} can be 0. We are adding a small value epsilon to avoid discontinuities which can cause small errors in the model.'.format(reactionID)) add_eps_react = True # let's instead split the rxn split_rxn = True rate = str(sym).replace("**", "^") for it in replace_dict.items(): rate = rate.replace(it[1], it[0]) return rate, "", 1, 1, False, split_rxn # prod_expr = back_expr removedR = [] for ibol, bol in enumerate(prod_symbols): stoi = int(prod[ibol][1]) # Now we can remove it prod_expr = prod_expr / (bol**stoi) removedR += [str(bol) for i in range(stoi)] # Check if we can get 0 in the denominator add_eps_prod = False n, d = prod_expr.as_numer_denom() for ibol, bol in enumerate(prod_symbols): if bol in d.atoms(): d = d.subs(bol, 0) if d == 0: # logMess('WARNING:RATE001', 'Denominator of rate constant in reaction {} can be 0. We are adding a small value epsilon to avoid discontinuities which can cause small errors in the model.'.format(reactionID)) # let's instead split the rxn split_rxn = True rate = str(sym).replace("**", "^") add_eps_prod = True for it in replace_dict.items(): rate = rate.replace(it[1], it[0]) return rate, "", 1, 1, False, split_rxn # prod_expr = prod_expr * -1 # TODO: We still need to figure out if we have # our reactant/products in our expressions and # if so set the nl/nr values accordingly # Reproducing current behavior + expansion re_proc = react_expr.nsimplify().evalf().simplify() pe_proc = prod_expr.nsimplify().evalf().simplify() # Adding epsilon if we have to if add_eps_react: # n,d = re_proc.as_numer_denom() # rateL = "(" + str(n) + ")/(" + str(d) + "+__epsilon__)" # self.write_epsilon = True # add_eps_prod = True # instead splitting the reaction split_rxn = True rate = str(sym).replace("**", "^") for it in replace_dict.items(): rate = rate.replace(it[1], it[0]) return rate, "", 1, 1, False, split_rxn else: rateL = str(re_proc) if add_eps_prod: # n,d = pe_proc.as_numer_denom() # rateR = "(" + str(n) + ")/(" + str(d) + "+__epsilon__)" # self.write_epsilon = True # add_eps_prod = True # instead splitting the reaction split_rxn = True rate = str(sym).replace("**", "^") for it in replace_dict.items(): rate = rate.replace(it[1], it[0]) return rate, "", 1, 1, False, split_rxn else: rateR = str(pe_proc) nl = self.calculate_factor(react, prod, rateL, removedL) nr = self.calculate_factor(prod, react, rateR, removedR) # nl, nr = 2, 2 # BNG power function is ^ and not ** rateL = rateL.replace("**", "^") rateR = rateR.replace("**", "^") reversible = True if rateL is None: # if not simply reversible, rely on the SBML spec if reversible: logMess( "WARNING:RATE002", "SBML claims reversibility of reaction {} but the kinetic law is not easily separated to forward and backward rate constants. Assuming irreversible reaction instead.".format( reactionID ), ) # Also get and parse the symbols react_bols = [x[0] for x in react] react_symbols = sympy.symbols(react_bols) # Now we can manipulate it react_expr = exp removedL = [] for ibol, bol in enumerate(react_symbols): stoi = int(react[ibol][1]) # Now we can remove it react_expr = react_expr / (bol**stoi) removedL += [str(bol) for i in range(stoi)] # Check if we can get the denominator to be 0 add_eps_react = False n, d = react_expr.as_numer_denom() for ibol, bol in enumerate(react_symbols): if bol in d.atoms(): d = d.subs(bol, 0) if d == 0: # logMess('WARNING:RATE001', 'Denominator of rate constant in reaction {} can be 0. We are adding a small value epsilon to avoid discontinuities which can cause small errors in the model.'.format(reactionID)) # add_eps_react = True # instead splitting the reaction split_rxn = True rate = str(sym).replace("**", "^") for it in replace_dict.items(): rate = rate.replace(it[1], it[0]) return rate, "", 1, 1, False, split_rxn re_proc = react_expr.nsimplify().evalf().simplify() if add_eps_react: # n,d = re_proc.as_numer_denom() # rateL = "(" + str(n) + ")/(" + str(d) + "+__epsilon__)" # self.write_epsilon = True # instead splitting the reaction split_rxn = True rate = str(sym).replace("**", "^") for it in replace_dict.items(): rate = rate.replace(it[1], it[0]) return rate, "", 1, 1, False, split_rxn else: rateL = str(re_proc) nl = self.calculate_factor(react, prod, rateL, removedL) rateL = rateL.replace("**", "^") # Make unidirectional rateR = "0" nr = -1 reversible = False # Check to ensure we don't have a negative # rate constant if rateL.startswith("-"): logMess( "WARNING:RATE003", "After analysis, forward rate constant of reaction {} starts with - which is likely an error.".format( reactionID ), ) MrateL = rateL MrateR = rateR Mnl, Mnr = nl, nr uRev = reversible for it in replace_dict.items(): MrateL = MrateL.replace(it[1], it[0]) MrateR = MrateR.replace(it[1], it[0]) return MrateL, MrateR, Mnl, Mnr, uRev, split_rxn
def __getRawRules( self, reaction, symmetryFactors, parameterFunctions, translator, sbmlfunctions ): zerospecies = ["emptyset", "trash", "sink", "source"] split_rxn = False if self.useID: reactant = [ ( reactant.getSpecies(), reactant.getStoichiometry(), reactant.getSpecies(), ) for reactant in reaction.getListOfReactants() if reactant.getSpecies().lower() not in zerospecies and reactant.getStoichiometry() not in [0, "0"] ] product = [ (product.getSpecies(), product.getStoichiometry(), product.getSpecies()) for product in reaction.getListOfProducts() if product.getSpecies().lower() not in zerospecies and product.getStoichiometry() not in [0, "0"] ] else: reactant = [ ( self.speciesDictionary[rElement.getSpecies()], rElement.getStoichiometry(), rElement.getSpecies(), ) for rElement in reaction.getListOfReactants() if self.speciesDictionary[rElement.getSpecies()].lower() not in zerospecies and rElement.getStoichiometry() not in [0, "0"] ] product = [ ( self.speciesDictionary[rProduct.getSpecies()], rProduct.getStoichiometry(), rProduct.getSpecies(), ) for rProduct in reaction.getListOfProducts() if self.speciesDictionary[rProduct.getSpecies()].lower() not in zerospecies and rProduct.getStoichiometry() not in [0, "0"] ] kineticLaw = reaction.getKineticLaw() reversible = reaction.getReversible() if kineticLaw is None: return { "reactants": reactant, "products": product, "parameters": [], "rates": ["0", "0"], "reversible": reversible, "reactionID": reaction.getId(), "numbers": [0, 0], "modifiers": None, "split_rxn": split_rxn, } rReactant = [] rProduct = [] # in case a given species was defined as the zero molecule don't include it in the rate correction method for x in reaction.getListOfReactants(): # check if we are splitting the reaction due to conversion factors if x.getSpecies() in self.bngModel.molecules: if self.bngModel.molecules[x.getSpecies()].conversionFactor is not None: split_rxn = True if ( x.getSpecies().lower() not in zerospecies and x.getStoichiometry() not in [ 0, "0", ] ): if not x.getConstant() and pymath.isnan(x.getStoichiometry()): logMess( "ERROR:SIM241", "BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format( reaction.getId() ), ) raise TranslationException( "ERROR:SIM241: BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format( reaction.getId() ) ) else: speciesName = ( self.speciesDictionary[x.getSpecies()] if x.getSpecies() in self.speciesDictionary else "" ) if ( speciesName in translator and str(translator[speciesName]) == "0" ): continue rReactant.append((x.getSpecies(), x.getStoichiometry())) for x in reaction.getListOfProducts(): # check if we are splitting the reaction due to conversion factors if x.getSpecies() in self.bngModel.molecules: if self.bngModel.molecules[x.getSpecies()].conversionFactor is not None: split_rxn = True if ( x.getSpecies().lower() not in zerospecies and x.getStoichiometry() not in [ 0, "0", ] ): if not x.getConstant() and pymath.isnan(x.getStoichiometry()): logMess( "ERROR:SIM207", "BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format( reaction.getId() ), ) raise TranslationException( "ERROR:SIM207: BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format( reaction.getId() ) ) else: speciesName = ( self.speciesDictionary[x.getSpecies()] if x.getSpecies() in self.speciesDictionary else "" ) if ( speciesName in translator and str(translator[speciesName]) == "0" ): continue rProduct.append((x.getSpecies(), x.getStoichiometry())) # rReactant = [(x.getSpecies(), x.getStoichiometry()) for x in reaction.getListOfReactants() if x.getSpecies() not in ['EmptySet']] # rProduct = [(x.getSpecies(), x.getStoichiometry()) for x in reaction.getListOfProducts() if x.getSpecies() not in ['EmptySet']] rModifiers = [ x.getSpecies() for x in reaction.getListOfModifiers() if x.getSpecies() != "EmptySet" ] parameters = [ (parameter.getId(), parameter.getValue(), parameter.getUnits()) for parameter in kineticLaw.getListOfParameters() ] rateL = rateR = nl = nr = None if True: # TODO: For some reason creating a deepcopy of this screws everything up, even # though its what we should be doing # update: apparently the solution was to use copy instead of deepcopy. This is because # the underlying swig code in c was causing conflicts when copied. make sure this actually works math = copy(kineticLaw.getMath()) math = math.deepCopy() # get a list of compartments so that we can remove them compartmentList = [] for compartment in self.model.getListOfCompartments(): compartmentList.append(compartment.getId()) # check for non-integer stoichiometry stoi = sum([r[1] for r in rReactant]) + sum([l[1] for l in rProduct]) if pymath.isnan(stoi): logMess( "ERROR:SIM208", f"Stoichiometry for reaction {reaction.getId()} is NaN, can't translate.", ) raise TranslationException( f"ERROR:SIM241 Stoichiometry for reaction {reaction.getId()} is NaN, can't translate." ) try: if stoi != int(stoi): non_integer_stoi = True else: non_integer_stoi = False except: logMess( "ERROR:SIM209", f"Stoichiometry for reaction {reaction.getId()} is {stoi}, can't translate.", ) raise TranslationException( f"ERROR:SIM208 Stoichiometry for reaction {reaction.getId()} is {stoi}, can't translate." ) # remove compartments from expression. also separate left hand and right hand side # check reaction molecularity and determine if it's # possible to have a mass action to begin with r_stoi = sum([r[1] for r in rReactant]) if r_stoi > 3 or non_integer_stoi: # modified analyze reaction rate call because we'll be splitting # the raection downstream split_rxn = True rateL, rateR, nl, nr, uReversible, split_rxn = self.analyzeReactionRate( math, compartmentList, reversible, rReactant, rProduct, reaction.getId(), parameterFunctions, rModifiers, sbmlfunctions, split_rxn, ) if rateR == "0": reversible = False # FIXME: make sure this actually works if symmetryFactors[0] > 1: rateL = "({0})*({1})".format(rateL, symmetryFactors[0]) if symmetryFactors[1] > 1: rateR = "({0})*({1})".format(rateR, symmetryFactors[1]) # we need to resolve observables BEFORE we do this for obs_key in self.obs_map: resL = re.search(r"(\W|^){0}(\W|$)".format(obs_key), rateL) if resL is not None: rateL = re.sub( r"(\W|^)({0})(\W|$)".format(obs_key), r"\g<1>{0}\g<3>".format(self.obs_map[obs_key]), rateL, ) resR = re.search(r"(\W|^){0}(\W|$)".format(obs_key), rateR) if resR is not None: rateR = re.sub( r"(\W|^)({0})(\W|$)".format(obs_key), r"\g<1>{0}\g<3>".format(self.obs_map[obs_key]), rateR, ) if not self.useID: srateL = self.convertToName(rateL) srateR = self.convertToName(rateR) else: srateL = rateL srateR = rateR if uReversible: pass # return compartments if the reaction is unimolecular # they were removed in the first palce because its easier to handle # around the equation in tree form when it has less terms """ if len(self.model.getListOfCompartments()) > 0: for compartment in (self.model.getListOfCompartments()): if compartment.getId() not in compartmentList: if len(rReactant) != 2: rateL = '{0} * {1}'.format(rateL,compartment.getSize()) if len(rProduct) != 2: rateR = '{0} * {1}'.format(rateR,compartment.getSize()) """ return { "reactants": reactant, "products": product, "parameters": parameters, "rates": [srateL, srateR], "reversible": uReversible, "reactionID": reaction.getId(), "numbers": [nl, nr], "modifiers": rModifiers, "split_rxn": split_rxn, }
[docs] def getReactionCenter(self, reactant, product, translator): rcomponent = Counter() pcomponent = Counter() for element in reactant: if element[0] in translator: for molecule in translator[element[0]].molecules: for component in molecule.components: molecule.sort() rcomponent.update( Counter( [ ( molecule.name, component.name, len(component.bonds) > 0, component.activeState, ) ] ) ) for element in product: if element[0] in translator: for molecule in translator[element[0]].molecules: molecule.sort() for component in molecule.components: molecule.sort() pcomponent.update( Counter( [ ( molecule.name, component.name, len(component.bonds) > 0, component.activeState, ) ] ) ) reactionCenter = [ (x[0], x[1]) for x in rcomponent for y in pcomponent if (x[0], x[1]) == (y[0], y[1]) and x != y and rcomponent[x] != pcomponent[y] ] rreactionCenter = [ (x[0], x[1]) for x in pcomponent for y in rcomponent if (x[0], x[1]) == (y[0], y[1]) and x != y and pcomponent[x] != rcomponent[y] ] return reactionCenter, rreactionCenter
[docs] def updateComponentCount(self, counterArray, reference, updateValue): for element in counterArray: if reference in counterArray[element]: counterArray[element][reference] += updateValue
[docs] def reduceComponentSymmetryFactors(self, reaction, translator, functions): """ create symmetry factors for reactions with components and species with identical names. This checks for symmetry in the components names then. """ # FIXME: This is entirely broken zerospecies = ["emptyset", "trash", "sink", "source"] if self.useID: reactant = [ (rElement.getSpecies(), rElement.getStoichiometry()) for rElement in reaction.getListOfReactants() if rElement.getSpecies() != "EmptySet" ] product = [ (product.getSpecies(), product.getStoichiometry()) for product in reaction.getListOfProducts() if product.getSpecies() != "EmptySet" ] else: reactant = [ ( self.speciesDictionary[rElement.getSpecies()], rElement.getStoichiometry(), ) for rElement in reaction.getListOfReactants() ] product = [ ( self.speciesDictionary[rProduct.getSpecies()], rProduct.getStoichiometry(), ) for rProduct in reaction.getListOfProducts() ] kineticLaw = reaction.getKineticLaw() reversible = reaction.getReversible() if kineticLaw is None: return 1, 1 rReactant = rProduct = [] for x in reaction.getListOfReactants(): if ( x.getSpecies().lower() not in zerospecies and x.getStoichiometry() not in [0, "0"] and pymath.isnan(x.getStoichiometry()) ): if not x.getConstant(): logMess( "ERROR:SIM241", "BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format( reaction.getId() ), ) return 1, 1 else: rReactant.append(x.getSpecies(), x.getStoichiometry()) for x in reaction.getListOfProducts(): if ( x.getSpecies().lower() not in zerospecies and x.getStoichiometry() not in [0, "0"] and pymath.isnan(x.getStoichiometry()) ): if not x.getConstant(): logMess( "ERROR:SIM241", "BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format( reaction.getId() ), ) return 1, 1 else: rProduct.append(x.getSpecies(), x.getStoichiometry()) # TODO: For some reason creating a deepcopy of this screws everything up, even # though its what we should be doing rcomponent = defaultdict(Counter) pcomponent = defaultdict(Counter) # get the total count of components in the reactants and products # e.g. components across diffent species freactionCenter, breactionCenter = self.getReactionCenter( reactant, product, translator ) for element in reactant: if element[0] in translator: for molecule in translator[element[0]].molecules: for component in molecule.components: molecule.sort() componentList = Counter([(molecule.signature(freactionCenter))]) for _ in range(0, int(element[1])): rcomponent[ ( molecule.name, component.name, len(component.bonds) > 0, component.activeState, ) ].update(componentList) for element in product: if element[0] in translator: for molecule in translator[element[0]].molecules: molecule.sort() for component in molecule.components: componentList = Counter([(molecule.signature(breactionCenter))]) for _ in range(0, int(element[1])): pcomponent[ ( molecule.name, component.name, len(component.bonds) > 0, component.activeState, ) ].update(componentList) """ only keep information for reaction centers """ reactionCenters = [ (x[0], x[1]) for x in rcomponent for y in pcomponent if (x[0], x[1]) == (y[0], y[1]) and x != y ] rcomponent = { x: rcomponent[x] for x in rcomponent if (x[0], x[1]) in reactionCenters } pcomponent = { x: pcomponent[x] for x in pcomponent if (x[0], x[1]) in reactionCenters } # is the number of components across products and reactants the same? # eg is there any DeleteMolecules action pcorrectionFactor = 1 rcorrectionFactor = 1 rStack = [] pStack = [] """ if a reaction can take place in several ways account for it in the reaction rate (this is specially important in dimer and trimer binding) pcomponent[element] < rcomponent[element] asks if an specific instance of a component decreases in number from a reactant to a product for example if there are 3 A(b)'s and one binds, we will have 2 A(b)'s in the product """ rcomponentTemp = deepcopy(rcomponent) pcomponentTemp = deepcopy(pcomponent) # calculate actual symmetry factors for key in rcomponent: if key in pcomponent: for element in rcomponent[key]: if rcomponent[key] == 1: continue # if theres a component on one side of the equation that # appears a different number of times on the other side of the equation if element in pcomponent[key]: if ( pcomponent[key][element] < rcomponent[key][element] and set([key[0].lower(), key[1].lower()]) not in rStack ): rcorrectionFactor *= comb( rcomponent[key][element], pcomponent[key][element], exact=1, ) rStack.append(set([key[0].lower(), key[1].lower()])) # once we choose a component for a previous action # this limits the options for subsequent actions # although substracting one from the target sites # may not be the right option. double check. self.updateComponentCount(pcomponent, element, -1) else: for element2 in pcomponent[key]: if ( pcomponent[key][element2] < rcomponent[key][element] and set([key[0].lower(), key[1].lower()]) not in rStack ): rcorrectionFactor *= comb( rcomponent[key][element], pcomponent[key][element2], exact=1, ) rStack.append(set([key[0].lower(), key[1].lower()])) self.updateComponentCount(pcomponent, element2, -1) rcomponent = rcomponentTemp pcomponent = pcomponentTemp if reversible: for key in pcomponent: if key in rcomponent: for element in pcomponent[key]: if pcomponent[key] == 1: continue if element in rcomponent[key]: if ( rcomponent[key][element] < pcomponent[key][element] and set([key[0].lower(), key[1].lower()]) not in pStack ): pcorrectionFactor *= comb( pcomponent[key][element], rcomponent[key][element], exact=1, ) pStack.append(set([key[0].lower(), key[1].lower()])) self.updateComponentCount(rcomponent, element, -1) else: for element2 in rcomponent[key]: if ( rcomponent[key][element2] < pcomponent[key][element] and set([key[0].lower(), key[1].lower()]) not in pStack ): pcorrectionFactor *= comb( pcomponent[key][element], rcomponent[key][element2], exact=1, ) pStack.append(set([key[0].lower(), key[1].lower()])) self.updateComponentCount(rcomponent, element2, -1) return rcorrectionFactor, pcorrectionFactor
[docs] def convertToName(self, rate): for element in sorted(self.speciesDictionary, key=len, reverse=True): while True: oldRate = rate if element in rate: rate = re.sub( r"(\W|^)({0})(\W|$)".format(element), r"\1{0}\3".format(self.speciesDictionary[element]), rate, ) if rate != oldRate: continue break # rate = rate.replace(element,self.speciesDictionary[element]) return rate
def __getRawCompartments(self, compartment): """ Private method used by the getCompartments method """ idid = compartment.getId() name = compartment.getName() size = compartment.getSize() # volume messes up the reactions # size = 1.0 dimensions = compartment.getSpatialDimensions() if dimensions in [0, 1]: logMess( "WARNING:SIM103", "{1}-D compartments are not supported. Changing for 2-D compartments for {0}. Please verify this does not affect simulation".format( name, dimensions ), ) dimensions = 2 if not compartment.getConstant(): # this compartment is not constant in size, not supported right now logMess( "ERROR:SIM203", f"Compartments with variable size is not supported, compartment: {idid}", ) return idid, dimensions, size, name def __getRawFunctions(self, function): math = function[1].getMath() name = function[1].getId() return name, libsbml.formulaToString(math)
[docs] def getSBMLFunctions(self): functions = {} for function in enumerate(self.model.getListOfFunctionDefinitions()): functionInfo = self.__getRawFunctions(function) if functionInfo[1] is None: logMess( "ERROR:SIM212", f"Math for function with ID '{functionInfo[0]}' is not defined", ) raise TranslationException( f"ERROR:SIM212: Math for function with ID '{functionInfo[0]}' is not defined" ) functions[functionInfo[0]] = writer.bnglFunction( functionInfo[1], functionInfo[0], [], reactionDict=self.reactionDictionary, ) return functions
[docs] def getCompartments(self): """ Returns an array of triples, where each triple is defined as (compartmentName,dimensions,size) """ compartments = [] # ASS - return empty if we are not using compartments units = None if self.noCompartment: return compartments unitDefinitions = self.getUnitDefinitions() if "volume" in unitDefinitions: units = "#volume units: {0}".format( "*".join([x["name"] for x in unitDefinitions["volume"]]) ) compartments.append(units) else: units = "#volume units: L" compartments.append(units) for _, compartment in enumerate(self.model.getListOfCompartments()): compartmentInfo = self.__getRawCompartments(compartment) comp_obj = self.bngModel.make_compartment() # ASS - removing "cell" as default compartment name = compartmentInfo[0] # making the obj before adding comp_obj.Id = name comp_obj.dim = compartmentInfo[1] comp_obj.size = compartmentInfo[2] comp_obj.unit = units if name != compartmentInfo[3] and compartmentInfo[3] != "": compartments.append( "%s %d %s #%s" % (name, compartmentInfo[1], compartmentInfo[2], compartmentInfo[3]) ) comp_obj.cmt = compartmentInfo[3] else: compartments.append( "%s %d %s" % (name, compartmentInfo[1], compartmentInfo[2]) ) self.bngModel.add_compartment(comp_obj) return compartments
[docs] def updateFunctionReference(self, reaction, updatedReferences): newRate = reaction[3] for reference in updatedReferences: newRate = re.sub( r"(\W|^)({0})(\W|$)".format(reference), r"\1{0}\3".format(updatedReferences[reference]), newRate, ) return newRate
[docs] def getSymmetryFactors(self, reaction): zerospecies = ["emptyset", "trash", "sink", "source"] if self.useID: reactant = [ (rElement.getSpecies(), rElement.getStoichiometry()) for rElement in reaction.getListOfReactants() if rElement.getSpecies() != "EmptySet" ] product = [ (product.getSpecies(), product.getStoichiometry()) for product in reaction.getListOfProducts() if product.getSpecies() != "EmptySet" ] else: reactant = [ ( self.speciesDictionary[rElement.getSpecies()], rElement.getStoichiometry(), ) for rElement in reaction.getListOfReactants() ] product = [ ( self.speciesDictionary[rProduct.getSpecies()], rProduct.getStoichiometry(), ) for rProduct in reaction.getListOfProducts() ] react_counts = {} for react in reactant: if react[0] in react_counts: react_counts[react[0]] += react[1] else: react_counts[react[0]] = react[1] if len(react_counts) == 0: lfact = 1 else: lfact = max(react_counts.values()) prod_counts = {} for prod in product: if prod[0] in prod_counts: prod_counts[prod[0]] += prod[1] else: prod_counts[prod[0]] = prod[1] if len(prod_counts) == 0: rfact = 1 else: rfact = max(prod_counts.values()) return lfact, rfact
[docs] def getReactions( self, translator={}, isCompartments=False, extraParameters={}, atomize=False, parameterFunctions={}, database=None, ): """ @returns: a triple containing the parameters,reactions,functions """ # @FIXME:this part of the code is there so that we only generate the functions list once through different # iterations of this call. This is because we cannot create a clone of the 'math' object for this # reaction and it is being permanently changed every call. It's ugly but it works. Change for something # better when we figure out how to clone the math object if not hasattr(self.getReactions, "functionFlag"): self.getReactions.__func__.functionFlag = False or (not atomize) reactions = [] reactionStructure = [] parameters = [] functions = [] # We want to keep track of the molecules/species we # actually used in the reactions functionTitle = "fRate" self.unitDefinitions = self.getUnitDefinitions() database.rawreactions = [] if len(self.model.getListOfReactions()) == 0: logMess( "WARNING:SIM104", "Model contains no natural reactions, all reactions are produced by SBML rules", ) # At the end of this loop, we want to go in and adjust any cases # where a value that can go to 0 is in the denominator for index, reaction in enumerate(self.model.getListOfReactions()): parameterDict = {} currParamConv = {} # symmetry factors for components with the same name # FIXME: This reduceComponentSymmetryFactors is completely broken # and will only give 1,1 right now # sl, sr = self.reduceComponentSymmetryFactors( # reaction, translator, functions # ) sl, sr = self.getSymmetryFactors(reaction) sbmlfunctions = self.getSBMLFunctions() try: rawRules = self.__getRawRules( reaction, [sl, sr], parameterFunctions, translator, sbmlfunctions ) database.rawreactions.append(rawRules) except TranslationException as e: if database and database.ignore: reactions.append( "#Reaction {0} is not correctly translated. check log for details".format( e.value ) ) continue else: raise TranslationException(e.value + " during reaction processing") # making the rule object rule_obj = self.bngModel.make_rule() rule_obj.rule_ind = index + 1 rule_obj.parse_raw(rawRules) rule_obj.symm_factors = [sl, sr] # Let's add our molecules for r in rawRules["reactants"]: if r[0] not in self.used_molecules: self.used_molecules.append(r[0]) for p in rawRules["products"]: if p[0] not in self.used_molecules: self.used_molecules.append(p[0]) if len(rawRules["parameters"]) > 0: for parameter in rawRules["parameters"]: """ if self.convertSubstanceUnits: #newParameter = self.convertToStandardUnits(parameter[1],self.unitDictionary) parameters.append('r%d_%s %f' % (index+1, parameter[0], parameter[1])) parameterDict[parameter[0]] = parameter[1] else: """ curr_param_name = "r%d_%s" % (index + 1, parameter[0]) curr_param = "r%d_%s %f" % (index + 1, parameter[0], parameter[1]) parameters.append(curr_param) parameterDict[parameter[0]] = parameter[1] currParamConv[parameter[0]] = curr_param_name # add these to parameters in obj param_obj = self.bngModel.make_parameter() param_obj.Id = curr_param_name param_obj.val = parameter[1] param_obj.rxn_ind = index + 1 self.bngModel.add_parameter(param_obj) compartmentList = [] compartmentList.extend( [ [self.__getRawCompartments(x)[0], self.__getRawCompartments(x)[2]] for x in self.model.getListOfCompartments() ] ) threshold = 0 if rule_obj.raw_num[0] > threshold or rule_obj.raw_rates[0] in translator: functionName = "%s%d()" % (functionTitle, index) else: # append reactionNumbers to parameterNames finalString = str(rule_obj.raw_rates[0]) for parameter in parameterDict: finalString = re.sub( r"(\W|^)({0})(\W|$)".format(parameter), r"\1{0}\3".format("r{0}_{1}".format(index + 1, parameter)), finalString, ) functionName = finalString if self.getReactions.functionFlag and "delay" in rule_obj.raw_rates[0]: logMess( "ERROR:SIM202", "BNG cannot handle delay functions in function %s" % functionName, ) fobj = self.bngModel.make_function() fobj.Id = functionName fobj.rule_ptr = rule_obj fobj.compartmentList = compartmentList if rule_obj.reversible: if ( rule_obj.raw_num[0] > threshold or rule_obj.raw_rates[0] in translator ): fobj.definition = rule_obj.raw_rates[0] if self.getReactions.functionFlag: # local parameter replacement flag if self.replaceLocParams: fstr = writer.bnglFunction( rule_obj.raw_rates[0], functionName, rawRules["reactants"], compartmentList, parameterDict, self.reactionDictionary, ) functions.append(fstr) fobj.local_dict = parameterDict self.bngModel.add_function(fobj) else: fstr = writer.bnglFunction( rule_obj.raw_rates[0], functionName, rawRules["reactants"], compartmentList, currParamConv, self.reactionDictionary, ) functions.append(fstr) fobj.local_dict = currParamConv self.bngModel.add_function(fobj) if ( rawRules["numbers"][1] > threshold or rule_obj.raw_rates[1] in translator ): functionName2 = "%s%dm()" % (functionTitle, index) fobj_2 = self.bngModel.make_function() fobj_2.Id = functionName2 fobj_2.rule_ptr = rule_obj fobj_2.definition = rule_obj.raw_rates[1] fobj_2.compartmentList = compartmentList if self.getReactions.functionFlag: # local parameter replacement flag if self.replaceLocParams: functions.append( writer.bnglFunction( rule_obj.raw_rates[1], functionName2, rule_obj.raw_prod, compartmentList, parameterDict, self.reactionDictionary, ) ) fobj_2.local_dict = parameterDict self.bngModel.add_function(fobj_2) else: functions.append( writer.bnglFunction( rule_obj.raw_rates[1], functionName2, rule_obj.raw_prod, compartmentList, currParamConv, self.reactionDictionary, ) ) fobj_2.local_dict = currParamConv self.bngModel.add_function(fobj_2) self.reactionDictionary[ rawRules["reactionID"] ] = "({0} - {1})".format(functionName, functionName2) finalRateStr = "{0},{1}".format(functionName, functionName2) rule_obj.rate_cts = (functionName, functionName2) else: finalString = str(rawRules["rates"][1]) for parameter in parameterDict: finalString = re.sub( r"(\W|^)({0})(\W|$)".format(parameter), r"\1{0}\3".format("r{0}_{1}".format(index + 1, parameter)), finalString, ) finalRateStr = "{0},{1}".format(functionName, finalString) rule_obj.rate_cts = (functionName, finalString) else: if ( rawRules["numbers"][0] > threshold or rawRules["rates"][0] in translator ): fobj.definition = rule_obj.raw_rates[0] if self.getReactions.functionFlag: # local parameter replacement flag if self.replaceLocParams: functions.append( writer.bnglFunction( rawRules["rates"][0], functionName, rawRules["reactants"], compartmentList, parameterDict, self.reactionDictionary, ) ) fobj.local_dict = parameterDict self.bngModel.add_function(fobj) else: functions.append( writer.bnglFunction( rawRules["rates"][0], functionName, rawRules["reactants"], compartmentList, currParamConv, self.reactionDictionary, ) ) fobj.local_dict = currParamConv self.bngModel.add_function(fobj) self.reactionDictionary[rawRules["reactionID"]] = "{0}".format( functionName ) finalRateStr = functionName rule_obj.rate_cts = (functionName,) reactants = [x for x in rawRules["reactants"]] products = [x for x in rawRules["products"]] modifierComment = ( "#Modifiers({0})".format(", ".join(rawRules["modifiers"])) if rawRules["modifiers"] else "" ) #### ADD RXN SEP HERE #### if rule_obj.raw_splt: # if we are splitting, we want to remove symmetry factor # from a fRate because it will spread over all reactants # and products if "fRate" in rule_obj.rate_cts[0]: # we have functional rate constant if int(rule_obj.symm_factors[0]) != 1: # we have a non-zero symmetry factor if rule_obj.rate_cts[0] not in self.bngModel.functions: logMess( "ERROR:SIM206", "Rate constant function needs adjusting but can't find function: {}" % functionName, ) defn = self.bngModel.functions[rule_obj.rate_cts[0]].definition self.bngModel.functions[ rule_obj.rate_cts[0] ].definition = f"({defn})/({rule_obj.symm_factors[0]})" if rule_obj.reversible: logMess( "ERROR:SIM205", "Splitting a reversible reaction, please check if correct, function: {}" % functionName, ) ctr = 0 # Now we write a single reaction for each # member with modified reaction rate constants # first RHS for reactant in reactants: r = reactant stoi = r[1] if int(stoi) != 1.0: nRateStr = "-1*{}*({})".format(stoi, finalRateStr) else: nRateStr = "-1*({})".format(finalRateStr) nr = (r[0], 1.0, r[2]) # adjust reaction name rxn_name = rawRules["reactionID"] + "_reactants_" + str(ctr) rxn_str = writer.bnglReaction( [], [nr], nRateStr, self.tags, translator, ( isCompartments or ( (len(reactants) == 0 or len(products) == 0) and self.getReactions.__func__.functionFlag ) ), rawRules["reversible"], reactionName=rxn_name, comment=modifierComment, ) reactions.append(rxn_str) # same thing for the model nrule_obj = self.bngModel.make_rule() nrule_obj.parse_raw(rawRules) nrule_obj.reversible = False nrule_obj.Id = rxn_name nrule_obj.rate_cts = (nRateStr,) nrule_obj.reactants = [] nrule_obj.products = [nr] self.bngModel.add_rule(nrule_obj) # tick the ctr ctr += 1 # then LHS ctr = 0 for product in products: p = product stoi = p[1] if int(stoi) != 1.0: nRateStr = "{}*{}".format(stoi, finalRateStr) else: nRateStr = "{}".format(finalRateStr) np = (p[0], 1.0, p[2]) # adjust reaction name rxn_name = rawRules["reactionID"] + "_products_" + str(ctr) rxn_str = writer.bnglReaction( [], [np], nRateStr, self.tags, translator, ( isCompartments or ( (len(reactants) == 0 or len(products) == 0) and self.getReactions.__func__.functionFlag ) ), rawRules["reversible"], reactionName=rxn_name, comment=modifierComment, ) reactions.append(rxn_str) # same thing for the model nrule_obj = self.bngModel.make_rule() nrule_obj.parse_raw(rawRules) nrule_obj.reversible = False nrule_obj.Id = rxn_name nrule_obj.rate_cts = (nRateStr,) nrule_obj.reactants = [] nrule_obj.products = [np] self.bngModel.add_rule(nrule_obj) # tick the ctr ctr += 1 #### END RXN SEP #### else: # add the rule rule_obj.reactants = reactants rule_obj.products = products self.bngModel.add_rule(rule_obj) rxn_str = writer.bnglReaction( reactants, products, finalRateStr, self.tags, translator, ( isCompartments or ( (len(reactants) == 0 or len(products) == 0) and self.getReactions.__func__.functionFlag ) ), rawRules["reversible"], reactionName=rawRules["reactionID"], comment=modifierComment, ) reactions.append(rxn_str) if atomize: self.getReactions.__func__.functionFlag = True self.bngModel.tags = self.tags return parameters, reactions, functions
[docs] def gather_terms(self, exp): pos, neg = [], [] l, r = exp.as_two_terms() resolve = [l, r] while len(resolve) > 0: elem = resolve.pop() if elem.is_Add: l, r = elem.as_two_terms() resolve += [l, r] else: # TODO: Do we have a better check? if str(elem).startswith("-"): neg.append(elem) else: pos.append(elem) # FIXME: Return None correctly l, r = None, None if len(pos) > 0: l = pos.pop(0) if len(pos) > 0: for e in pos: l += e if len(neg) > 0: r = -1 * neg.pop(0) if len(neg) > 0: for e in neg: r += -1 * e return l, r
def __getRawAssignmentRules(self, arule): """ Function to pull out rate constants (if going to be encoded as a reaction) and properties of an assignment rule to be parsed downstream. """ variable = arule.getVariable() var = sympy.symbols(variable) if variable not in self.all_syms.keys(): self.all_syms[variable] = var # Sympy stuff form, replace_dict = self.find_all_symbols(arule.getMath(), None) # We might need this for debugging complicated # piecewise function forms. try: sym = sympy.sympify(form, locals=self.all_syms) except SympifyError as e: logMess( "ERROR:SYMP001", "Sympy couldn't parse a function, sorry but we can't handle this function definition.", ) raise TranslationException(arule.getId()) # if it's an assignment, it's going to be # encoded as a unidirectional rxn if not arule.isAssignment(): # expand and take the terms out as left and right exp = sympy.expand(sym) rateL = None rateR = None # This shows if we can get X - Y if exp.is_Add: react_expr, prod_expr = self.gather_terms(exp) if react_expr is None: # TODO: LogMess this print("no forward reaction rate?") # Let's also ensure that we have a + and - term elif prod_expr is not None: # Remove mass action prod_expr = prod_expr / var # Check for epsilon add_eps_prod = False n, d = prod_expr.as_numer_denom() for bol in d.atoms(): d = d.subs(bol, 0) if d == 0: add_eps_prod = True # Reproducing current behavior + expansion re_proc = react_expr.nsimplify().evalf().simplify() pe_proc = prod_expr.nsimplify().evalf().simplify() # Adding epsilon if we have to rateL = str(re_proc) if add_eps_prod: rateL = str(sym) rateR = "0" else: rateR = str(pe_proc) rateL = rateL.replace("**", "^") rateR = rateR.replace("**", "^") # slap back in the replacement stuff for it in replace_dict.items(): rateL = rateL.replace(it[1], it[0]) rateR = rateR.replace(it[1], it[0]) if rateL is None: # if not simply reversible, rely on the SBML spec react_expr = exp re_proc = react_expr.nsimplify().evalf().simplify() rateL = str(re_proc) rateL = rateL.replace("**", "^") # Make unidirectional rateR = "0" # slap back in the replacement stuff for it in replace_dict.items(): rateL = rateL.replace(it[1], it[0]) else: rateL = str(sym) rateR = "0" # slap back in the replacement stuff for it in replace_dict.items(): rateL = rateL.replace(it[1], it[0]) if not self.useID: rateL = self.convertToName(rateL) rateR = self.convertToName(rateR) return variable, [rateL, rateR], arule.isAssignment(), arule.isRate()
[docs] def adjustInitialConditions( self, parameters, initialConditions, artificialObservables, observables, functions, ): """ assignment rules require further adjustment after parsed to their initial values. """ # We'll need this to write nicely prnter = StrPrinter({"full_prec": False}) # This gets the initial conditions we have to adjust initCondsToAdjust = [ x for x in artificialObservables.keys() if x.endswith("_ar") ] # return original values if nothing needs adjusted if len(initCondsToAdjust) == 0: return initialConditions logMess("DEBUG:ARUL001", "adjusting assignment rule initial conditions") # These are the formulas we need to calculate dynamically formulaToAdjustWith = {} for cond in initCondsToAdjust: form = artificialObservables[cond] name = form.split("=")[0].split()[0] for func in functions: splt = func.split("=") if len(splt) > 2: logMess( "WARNING:ARUL002", "There are multiple = signs in function {}".format(splt[0]), ) if name == splt[0].split()[0]: # This avoids issues with multiple = signs formulaToAdjustWith[name] = splt[-1] # Now we have the formulas to dynamically calculate, let's get a parameter dictionary to replace parameters down the line param_dict = dict([(x.split()[0], x.split()[1]) for x in parameters]) # we need to parse observables to map # names to species descriptions obs_map = {} for iobs, obs in enumerate(observables): splt = obs.split() obs_map[splt[2]] = splt[1] # Now that we have observables, we can map the names # to their initial values to be replaced later initValMap = {} initCondSplit = [] for initCond in initialConditions: splt = initCond.split() initCondSplit.append(splt) # I'm a bit vary of this, not sure if this is # the only way the $ might appear honestly # keep an eye out for bugs here if splt[0].startswith("$"): check_name = splt[0][1:] else: check_name = splt[0] # if the name is in the observable species defs if check_name in obs_map.keys(): # we slap that into our initial value map initValMap[obs_map[check_name]] = splt[1] # This is to allow us to just use .split to be able to # separate species/parameter stuff from math stuff sympy_forms = {} for item in formulaToAdjustWith.items(): name, form = item try: sympy_forms[name] = sympy.sympify(form, locals=self.all_syms) except: sympy_forms[name] = None sympy_subs = {} # We need to replace the values in sympy formulas for item in sympy_forms.items(): name, form = item if form is not None: # Looping over parameters and replacing for par_item in param_dict.items(): pn, pv = par_item pv = pv.replace("(", "").replace(")", "") try: ps, pvs = sympy.symbols(pn), sympy.Number(pv) except: logMess( "ERROR:SYMP003", "Sympy couldn't parse an initial condition: {} and {}, this might cause issues for the initial state of your model".format( pn, pv ), ) continue form = form.subs(ps, pvs) # Replacing species from initial conditions for spec_item in initValMap.items(): sn, sv = spec_item sv = sv.replace("(", "").replace(")", "") try: ss, svs = sympy.symbols(sn), sympy.Number(sv) except: logMess( "ERROR:SYMP003", "Sympy couldn't parse an initial condition: {} and {}, this might cause issues for the initial state of your model".format( sn, sv ), ) continue form = form.subs(ss, svs) # And now converting the rest to zeros so we have a value # that we can fully evaluate for atm in form.atoms(sympy.Symbol): form = form.subs(atm, 0) sympy_subs[name] = form.evalf() else: sympy_subs[name] = 0 # Evaluate them using Sympy adjustedInitialValues = [] for item in sympy_subs.items(): name, val = item try: val_str = prnter.doprint(val) if val_str == "nan": val = 0 adjustedInitialValues.append(val) except: adjustedInitialValues.append(0) # Re-write initial conditions adjustedInitialConditions = [] for iic, initCond in enumerate(initCondSplit): toAdjustCheck = initCond[-1].replace("#", "") + "_ar" if toAdjustCheck in initCondsToAdjust: val = adjustedInitialValues[initCondsToAdjust.index(toAdjustCheck)] initCondSplit[iic][1] = "{}".format(prnter.doprint(val)) adjustedInitialConditions.append(" ".join(initCond)) else: adjustedInitialConditions.append(" ".join(initCond)) return adjustedInitialConditions
[docs] def getAssignmentRules( self, zparams, parameters, molecules, observablesDict, translator ): """ this method obtains an SBML rate rules and assignment rules. They require special handling since rules are often both defined as rules and parameters initialized as 0, so they need to be removed from the parameters list """ # FIXME: This function removes compartment info and this leads to mis-replacement of variables downstream. e.g. Calc@ER and Calc@MIT both gets written as Calc and downstream the replacement is wrong. # FIXME: This function gets a list of observables which sometimes are turned into assignment rules but then are not updated in the observablesDict. E.g. X_comp1 gets in, X_ar is created and you can't have BOTH X_comp1 in a reaction AND X_ar adjusting X itself. You MUST pick one, if both are happening raise and error and exit out. For now I'll say if we have _ar then we replace the X_comp1 with X_ar and test. # Going to use this to match names and remove params # if need be param_map = dict([(x.split()[0], x) for x in parameters]) compartmentList = [] compartmentList.extend( [ [self.__getRawCompartments(x)[0], self.__getRawCompartments(x)[2]] for x in self.model.getListOfCompartments() ] ) arules = [] aParameters = {} zRules = zparams removeParameters = [] artificialReactions = [] artificialObservables = {} nonamecounter = 0 for arule in self.model.getListOfRules(): rawArule = self.__getRawAssignmentRules(arule) # print("looping over assignment rules") arule_obj = self.bngModel.make_arule() arule_obj.parse_raw(rawArule) arule_obj.compartmentList = compartmentList # rule has no name if arule_obj.Id == "": logMess( "ERROR:SIM215", "atomizer has found an sbml rule without a name. {0}".format( rawArule[1:] ), ) rawArule = list(rawArule) new_arule_name = "noname{0}".format(nonamecounter) rawArule[0] = new_arule_name nonamecounter += 1 arule_obj.Id = new_arule_name if arule_obj.isRate == True: # it is a rate rule if rawArule[0] in self.boundaryConditionVariables: logMess( "WARNING:SIM105", "rate rules ({0}) are not properly supported in BioNetGen simulator".format( rawArule[0] ), ) rateLaw1 = arule_obj.rates[0] rateLaw2 = arule_obj.rates[1] # TODO: Add to bngModel functions arate_name = "arRate{0}".format(rawArule[0]) func_str = writer.bnglFunction( rateLaw1, arate_name, [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) arules.append(func_str) if rateLaw2 != "0": # TODO: Add to bngModel functions armrate_name = "armRate{0}".format(rawArule[0]) func2_str = writer.bnglFunction( rateLaw2, armrate_name, [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) arules.append(func2_str) # ASS2019 - I'm not sure if this is the right place to fix the tags. Basically, up until this point, the artificial reactions don't have tags. This results in the 0 <-> A type reactions to lack a compartment, leading to a non-functional BNGL file. I think the better solution might be during rule (SBML rule, not BNGL rule) parsing and update the parser/SBML2BNGL tags instead. try: comp = self.tags[rawArule[0]] except KeyError: # ASS - We need to default to an existing compartment if we are going # to remove @cell as a default compartment if not self.noCompartment: self.tags[rawArule[0]] = "@" + compartmentList[0][0] # ASS - If self.useID is set, use the ID value, not the name if self.useID: self.used_molecules.append(rawArule[0]) if rateLaw2 == "0": rxn_str = writer.bnglReaction( [], [[rawArule[0], 1, rawArule[0]]], "{0}".format("arRate{0}".format(rawArule[0])), self.tags, translator, isCompartments=True, comment="#rateLaw", reversible=False, ) else: rxn_str = writer.bnglReaction( [], [[rawArule[0], 1, rawArule[0]]], "{0},{1}".format( "arRate{0}".format(rawArule[0]), "armRate{0}".format(rawArule[0]), ), self.tags, translator, isCompartments=True, comment="#rateLaw", ) artificialReactions.append(rxn_str) else: self.used_molecules.append(self.convertToName(rawArule[0]).strip()) if rateLaw2 == "0": rxn_str = writer.bnglReaction( [], [[self.convertToName(rawArule[0]).strip(), 1, rawArule[0]]], "{0}".format("arRate{0}".format(rawArule[0])), self.tags, translator, isCompartments=True, comment="#rateLaw", reversible=False, ) else: rxn_str = writer.bnglReaction( [], [[self.convertToName(rawArule[0]).strip(), 1, rawArule[0]]], "{0},{1}".format( "arRate{0}".format(rawArule[0]), "armRate{0}".format(rawArule[0]), ), self.tags, translator, isCompartments=True, comment="#rateLaw", ) artificialReactions.append(rxn_str) if rawArule[0] in zparams: removeParameters.append("{0} 0".format(rawArule[0])) zRules.remove(rawArule[0]) else: for element in parameters: # TODO: if for whatever reason a rate rule # was defined as a parameter that is not 0 # remove it. This might not be exact behavior if re.search("^{0}\s".format(rawArule[0]), element): logMess( "WARNING:SIM106", "Parameter {0} corresponds both as a non zero parameter \ and a rate rule, verify behavior".format( element ), ) removeParameters.append(element) # it is an assigment rule elif arule_obj.isAssignment: """ normal species observables references in functions keep the format <speciesName>_<compartment> in function references, and observables dict keeps track of that. however when a species is defined by an assignment function we wish to keep track of reference <speciesName> that points to a standard BNGL function """ # it was originially defined as a zero parameter, so delete it from the parameter list definition if rawArule[0] in zRules: # dont show assignment rules as parameters zRules.remove(rawArule[0]) matches = [ molecules[x] for x in molecules if molecules[x]["name"] == rawArule[0] ] if matches: if matches[0]["isBoundary"]: artificialObservables[ rawArule[0] + "_ar" ] = writer.bnglFunction( rawArule[1][0], rawArule[0] + "_ar()", [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) self.arule_map[rawArule[0]] = rawArule[0] + "_ar" if rawArule[0] in observablesDict: observablesDict[rawArule[0]] = rawArule[0] + "_ar" continue else: logMess( "ERROR:SIM201", "Variables that are both changed by an assignment rule and reactions are not \ supported in BioNetGen simulator. The variable will be split into two".format( rawArule[0] ), ) artificialObservables[ rawArule[0] + "_ar" ] = writer.bnglFunction( rawArule[1][0], rawArule[0] + "_ar()", [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) self.arule_map[rawArule[0]] = rawArule[0] + "_ar" if rawArule[0] in observablesDict: observablesDict[rawArule[0]] = rawArule[0] + "_ar" continue elif rawArule[0] in [observablesDict[x] for x in observablesDict]: artificialObservables[ rawArule[0] + "_ar" ] = writer.bnglFunction( rawArule[1][0], rawArule[0] + "_ar()", [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) self.arule_map[rawArule[0]] = rawArule[0] + "_ar" if rawArule[0] in observablesDict: observablesDict[rawArule[0]] = rawArule[0] + "_ar" continue elif rawArule[0] in molecules: if molecules[rawArule[0]]["isBoundary"]: # We should probably re-write this with the name since that's what's used other places name = molecules[rawArule[0]]["returnID"] artificialObservables[name + "_ar"] = writer.bnglFunction( rawArule[1][0], name + "_ar()", [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) self.arule_map[rawArule[0]] = name + "_ar" # TODO: Let's store what we know are assignment rules. We can maybe assume that, if something has an assignment rule, it can't in turn be in a reaction? If this is wrong, we can't model this anyway, so we should probably just make an assumption and let people know. self.only_assignment_dict[name] = name + "_ar" self.bngModel.add_arule(arule_obj) continue else: # if not boundary but is a species, Jose # is turning this into an assignment rule # with a different name (uses ID). # It looks as if the goal was to handle # both situations via renaming. # FIXME: This is very likely broken but # I'm not 100% sure how it breaks things. # TODO: Check, if we have this in observables we need to adjust the observablesDict because we are writing an assignment rule for this instead name = molecules[rawArule[0]]["returnID"] artificialObservables[name + "_ar"] = writer.bnglFunction( rawArule[1][0], name + "_ar()", [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) self.arule_map[rawArule[0]] = name + "_ar" self.only_assignment_dict[name] = name + "_ar" if name in observablesDict: observablesDict[name] = name + "_ar" self.bngModel.add_arule(arule_obj) continue else: # check if it is defined as an observable # FIXME: This doesn't check for parameter namespace # TODO: What is going on here? candidates = [ idx for idx, x in enumerate(observablesDict) if rawArule[0] == x ] assigObsFlag = False for idx in candidates: # if re.search('\s{0}\s'.format(rawArule[0]),observables[idx]): artificialObservables[ rawArule[0] + "_ar" ] = writer.bnglFunction( rawArule[1][0], rawArule[0] + "_ar()", [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) self.arule_map[rawArule[0]] = rawArule[0] + "_ar" assigObsFlag = True break if assigObsFlag: continue # if its not a param/species/observable # TODO: now, if we replace this with the returnID do we # overlap with it's other possible uses? # name = molecules[rawArule[0]]['returnID'] # self.only_assignment_dict[name] = name+"_ar" # artificialObservables[name+'_ar'] = writer.bnglFunction(rawArule[1][0],name+'()',[],compartments=compartmentList,reactionDict=self.reactionDictionary) # This doesn't actually check for clashes with # parameter namespace if rawArule[0] in param_map.keys(): removeParameters.append(param_map[rawArule[0]]) artificialObservables[rawArule[0] + "_ar"] = writer.bnglFunction( rawArule[1][0], rawArule[0] + "()", [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) self.arule_map[rawArule[0]] = rawArule[0] + "_ar" else: """ if for whatever reason you have a rule that is not assigment or rate and it is initialized as a non zero parameter, give it a new name """ if rawArule[0] not in zparams: ruleName = "ar" + rawArule[0] else: ruleName = rawArule[0] zRules.remove(rawArule[0]) arules.append( writer.bnglFunction( rawArule[1][0], ruleName, [], compartments=compartmentList, reactionDict=self.reactionDictionary, ) ) aParameters[rawArule[0]] = "ar" + rawArule[0] """ elif rawArule[2] == True: for parameter in parameters: if re.search('^{0}\s'.format(rawArule[0]),parameter): print '////',rawArule[0] """ # we can't decide any of this here, we need the # full model setup before deciding self.bngModel.add_arule(arule_obj) return ( aParameters, arules, zRules, artificialReactions, removeParameters, artificialObservables, )
[docs] def convertToStandardUnits(self, parameterValue, unitDefinition): for factor in unitDefinition: if factor["scale"] != 0: parameterValue *= 10 ** factor["scale"] if factor["multiplier"] != 1: parameterValue *= factor["multiplier"] if factor["exponent"] != 1: parameterValue **= factor["exponent"] return parameterValue
[docs] def convertToStandardUnitString(self, parameterValue, unitDefinition): for factor in unitDefinition: if factor["multiplier"] != 1: parameterValue = "({0} * {1})".format( parameterValue, factor["multiplier"] ) if factor["exponent"] != 1: parameterValue = "({0} ^ {1})".format( parameterValue, factor["exponent"] ) if factor["scale"] != 0: parameterValue = "({0} * 1e{1})".format(parameterValue, factor["scale"]) return parameterValue
def __getRawParameters(self, parameter): parameterSpecs = {} parameterSpecs["id"] = parameter.getId() parameterSpecs["value"] = parameter.getValue() parameterSpecs["name"] = parameter.getName() parameterSpecs["units"] = parameter.getUnits() return parameterSpecs
[docs] def getParameters(self): parameters = [] zparam = [] for parameter in self.model.getListOfParameters(): parameterSpecs = ( parameter.getId(), parameter.getValue(), parameter.getConstant(), parameter.getUnits(), ) # reserved keywords param_obj = self.bngModel.make_parameter() if parameterSpecs[0] == "e": # TODO: raise a warning parameterSpecs = ("__e__", parameterSpecs[1]) self.param_repl["e"] = "__e__" if parameterSpecs[1] == 0: zparam.append(parameterSpecs[0]) else: if parameter.getUnits() != "": parameters.append( "{0} {1} #units:{2}".format( parameterSpecs[0], parameterSpecs[1], parameter.getUnits() ) ) param_obj.units = parameter.getUnits() else: parameters.append( "{0} {1}".format(parameterSpecs[0], parameterSpecs[1]) ) if len(parameterSpecs) == 3: param_obj.Id = parameterSpecs[0] param_obj.val = parameterSpecs[1] param_obj.cts = parameterSpecs[2] else: param_obj.Id = parameterSpecs[0] param_obj.val = parameterSpecs[1] param_obj.cts = True self.bngModel.add_parameter(param_obj) self.bngModel.param_repl = self.param_repl return parameters, zparam
# ASS - trying to implement a more complete solution for # non compartmental models
[docs] def check_noCompartment(self, parameters=[]): self.compartmentDict = {} # self.compartmentDict[""] = 1 comps = [i for i in self.model.getListOfCompartments()] if len(comps) == 0: self.noCompartment = True self.bngModel.noCompartment = True return for compartment in self.model.getListOfCompartments(): self.compartmentDict[compartment.getId()] = compartment.getSize() self.noCompartment = False self.bngModel.noCompartment = False # Get all rawSpecies allRawSpecies = map( lambda x: self.getRawSpecies(x, parameters), self.model.getListOfSpecies() ) # Pull all used compartments in those speices and make a set allUsedCompartments = set(map(lambda x: x["compartment"], allRawSpecies)) # Basically, if we are using a single compartment of volume 1, # we are not actually using a compartment, not from the POV of # BNGL. To improve readability, we want to just write a regular # BNGL model instead of a cBNGL model. Especially true since # this is the case for most SBML models. if len(allUsedCompartments) == 1: # We are using only 1 compartment, check volume # FIXME: We will try removing the compartment # if only one is used # self.noCompartment = True # self.bngModel.noCompartment = True if self.compartmentDict[allUsedCompartments.pop()] == 1: # we have 1 compartment and it's volume is 1 # just don't use compartments. self.noCompartment = True self.bngModel.noCompartment = True elif len(allUsedCompartments) > 1: if not self.multi_comp_warn: logMess( "WARNING:COMP001", "Multiple compartments are used, please note that Atomizer does not automatically try to infer your compartment topology which is important for how rules fire in cBNGL. Make sure your comparment topology is set correctly after translation", ) self.multi_comp_warn = True self.speciesMemory = []
[docs] def getSpecies(self, translator={}, parameters=[]): """ in sbml parameters and species have their own namespace. not so in bionetgen, so we need to rename things if they share the same name """ def default_to_regular(d): if isinstance(d, defaultdict): d = {k: default_to_regular(v) for k, v in d.items()} return d # gotta reset the bngModel everytime # this function is called self.bngModel._reset() # find concentration units unitDefinitions = self.getUnitDefinitions() if "substance" in unitDefinitions: substance = "*".join([x["name"] for x in unitDefinitions["substance"]]) else: substance = "mol" if "volume" in unitDefinitions: volume = "/".join([x["name"] for x in unitDefinitions["volume"]]) else: volume = "L" concentrationUnits = "{0}/{1}".format(substance, volume) annotationMoleculesText = {} moleculesText = [] speciesText = [] observablesText = [] observablesDict = {} names = [] rawSpeciesName = list(translator.keys()) speciesTranslationDict = {} compartmentDict = {} compartmentDict[""] = 1 speciesAnnotationInfo = default_to_regular(self.getFullAnnotation()) annotationInfo = {"moleculeTypes": {}, "species": {}} for compartment in self.model.getListOfCompartments(): compartmentDict[compartment.getId()] = compartment.getSize() unitFlag = True for species in self.model.getListOfSpecies(): # making molecule and seed species objs for # the obj based model molec_obj = self.bngModel.make_molecule() spec_obj = self.bngModel.make_species() # rawSpecies = self.getRawSpecies(species, parameters) # letting the objs parse the rawSpecies molec_obj.parse_raw(rawSpecies) spec_obj.parse_raw(rawSpecies) if rawSpecies["compartment"] != "": # ASS - First change for "noCompartments" if self.noCompartment: rawSpecies["compartment"] = "" self.tags[rawSpecies["identifier"]] = "" else: self.tags[rawSpecies["identifier"]] = "@%s" % ( rawSpecies["compartment"] ) if rawSpecies["returnID"] in translator: if rawSpecies["returnID"] in rawSpeciesName: rawSpeciesName.remove(rawSpecies["returnID"]) if ( translator[rawSpecies["returnID"]].getSize() == 1 and translator[rawSpecies["returnID"]].molecules[0].name not in names and translator[rawSpecies["returnID"]].molecules[0].name not in rawSpeciesName ): names.append(translator[rawSpecies["returnID"]].molecules[0].name) annotationTemp = [] if rawSpecies["returnID"] in speciesAnnotationInfo: for annotation in speciesAnnotationInfo[rawSpecies["returnID"]]: parts = annotation.split("_") header = annotationHeader[parts[0]] qual = parts[1].lower() + "".join( [x.capitalize() for x in parts[2:]] ) entry = ", ".join( [ ":".join(x.split("/")[-2:]) for x in speciesAnnotationInfo[ rawSpecies["returnID"] ][annotation] ] ) annotationTemp.append( "#^ {0}:{1} {2}".format(header, qual, entry) ) # we'll add this to our model self.bngModel.add_molecule(molec_obj) mtext = translator[rawSpecies["returnID"]].str2() moleculesText.append(mtext) if rawSpecies["returnID"] in speciesAnnotationInfo: annotationInfo["moleculeTypes"][ translator[rawSpecies["returnID"]].str2() ] = annotationTemp del speciesAnnotationInfo[rawSpecies["returnID"]] # TODO: Not sure if there are more examples of this # but glucose in 380 has both a normal species AND # a boundary species separately # elif rawSpecies['isBoundary']: # self.bngModel.add_molecule(molec_obj) else: # we'll add this to our model self.bngModel.add_molecule(molec_obj) mtext = rawSpecies["returnID"] + "()" moleculesText.append(mtext) if rawSpecies["returnID"] in speciesAnnotationInfo: annotationInfo["moleculeTypes"][ rawSpecies["returnID"] ] = speciesAnnotationInfo[rawSpecies["returnID"]] del speciesAnnotationInfo[rawSpecies["returnID"]] # if rawSpecies['identifier'] == 'glx' and len(translator) > 0: temp = "$" if rawSpecies["isConstant"] != 0 else "" tmp = ( translator[str(rawSpecies["returnID"])] if rawSpecies["returnID"] in translator else rawSpecies["returnID"] + "()" ) # this determines the name to be written if ( rawSpecies["initialConcentration"] > 0 or rawSpecies["initialAmount"] > 0 ): tmp2 = temp if rawSpecies["identifier"] in self.tags: tmp2 = self.tags[rawSpecies["identifier"]] if rawSpecies["initialAmount"] > 0.0: # Removing the compartment section if we are not using it if self.noCompartment: speciesText.append( "{1}{2} {3} #{4} #{5}".format( tmp2, temp, str(tmp), rawSpecies["initialAmount"], rawSpecies["returnID"], rawSpecies["identifier"], ) ) else: speciesText.append( "{0}:{1}{2} {3} #{4} #{5}".format( tmp2, temp, str(tmp), rawSpecies["initialAmount"], rawSpecies["returnID"], rawSpecies["identifier"], ) ) elif rawSpecies["initialConcentration"] > 0.0: if self.isConversion: # convert to molecule counts if "substance" in unitDefinitions: newParameterStr = self.convertToStandardUnitString( rawSpecies["initialConcentration"], unitDefinitions["substance"], ) newParameter = self.convertToStandardUnits( rawSpecies["initialConcentration"], unitDefinitions["substance"], ) # conversion to moles else: newParameter = rawSpecies["initialConcentration"] newParameterStr = str(rawSpecies["initialConcentration"]) newParameter = ( newParameter * 6.022e23 ) # convertion to molecule counts # get compartment size if self.noCompartment: compartmentSize = 1.0 else: compartmentSize = self.model.getCompartment( rawSpecies["compartment"] ).getSize() newParameter = compartmentSize * newParameter if unitFlag: if self.noCompartment: speciesText.append( "{1}{2} {3} # {4}mol/L * 6.022e23/mol *{7}L #{5} #{6}".format( tmp2, temp, str(tmp), newParameter, newParameterStr, rawSpecies["returnID"], rawSpecies["identifier"], compartmentSize, concentrationUnits, ) ) else: speciesText.append( "{0}:{1}{2} {3} # {4}mol/L * 6.022e23/mol *{7}L #{5} #{6}".format( tmp2, temp, str(tmp), newParameter, newParameterStr, rawSpecies["returnID"], rawSpecies["identifier"], compartmentSize, concentrationUnits, ) ) unitFlag = False else: if self.noCompartment: speciesText.append( "{1}{2} {3} #original {4}{8} #{5} #{6}".format( tmp2, temp, str(tmp), newParameter, rawSpecies["initialConcentration"], rawSpecies["returnID"], rawSpecies["identifier"], compartmentSize, concentrationUnits, ) ) else: speciesText.append( "{0}:{1}{2} {3} #original {4}{8} #{5} #{6}".format( tmp2, temp, str(tmp), newParameter, rawSpecies["initialConcentration"], rawSpecies["returnID"], rawSpecies["identifier"], compartmentSize, concentrationUnits, ) ) else: if self.noCompartment: speciesText.append( "{1}{2} {3} #{4} #{5}".format( tmp2, temp, str(tmp), rawSpecies["initialConcentration"], rawSpecies["returnID"], rawSpecies["identifier"], ) ) else: speciesText.append( "{0}:{1}{2} {3} #{4} #{5}".format( tmp2, temp, str(tmp), rawSpecies["initialConcentration"], rawSpecies["returnID"], rawSpecies["identifier"], ) ) elif rawSpecies["isConstant"]: if self.noCompartment: speciesText.append( "{1}{2} {3} #{4} #{5}".format( tmp2, temp, str(tmp), 0, rawSpecies["returnID"], rawSpecies["identifier"], ) ) else: speciesText.append( "{0}:{1}{2} {3} #{4} #{5}".format( tmp2, temp, str(tmp), 0, rawSpecies["returnID"], rawSpecies["identifier"], ) ) self.bngModel.add_species(spec_obj) if rawSpecies["returnID"] == "e": modifiedName = "__e__" else: modifiedName = rawSpecies["returnID"] # user defined zero molecules are not included in the observable list if str(tmp) != "0": if ( rawSpecies["compartment"] != "" and len(list(self.model.getListOfCompartments())) > 1 ): self.obs_names.append(modifiedName) # self.obs_map[rawSpecies["identifier"]] = "{0}_{1}".format( # modifiedName, rawSpecies["compartment"] # ) # observablesText.append( # "Species {0}_{3} @{3}:{1} #{2}".format( # modifiedName, # tmp, # rawSpecies["name"], # rawSpecies["compartment"], # ) # ) # observablesDict[modifiedName] = "{0}_{1}".format( # modifiedName, rawSpecies["compartment"] # ) self.obs_map[rawSpecies["identifier"]] = "{0}".format( modifiedName, rawSpecies["compartment"] ) observablesText.append( "Species {0} @{3}:{1} #{2}".format( modifiedName, tmp, rawSpecies["name"], rawSpecies["compartment"], ) ) observablesDict[modifiedName] = "{0}".format( modifiedName, rawSpecies["compartment"] ) else: # ASS - Is this not supposed to be the version without compartments? self.obs_names.append(modifiedName) self.obs_map[rawSpecies["identifier"]] = modifiedName observablesText.append( "Species {0} {1} #{2}".format( modifiedName, tmp, rawSpecies["name"] ) ) observablesDict[modifiedName] = "{0}".format(modifiedName) speciesTranslationDict[rawSpecies["identifier"]] = tmp # add the observable in the model obs_obj = self.bngModel.make_observable() obs_obj.parse_raw(rawSpecies) obs_obj.Id = modifiedName self.bngModel.add_observable(obs_obj) # TODO: make sure this is replicated in bngModel sorted(rawSpeciesName, key=len) for species in rawSpeciesName: if ( translator[species].getSize() == 1 and translator[species].molecules[0].name not in names ): names.append(translator[species].molecules[0].name) mtext = translator[species].str2() moleculesText.append(mtext) molec_obj = self.bngModel.make_molecule() molec_obj.Id = species # TODO: Make sure we need str2 and not # just str if not str(translator[species]) in self.bngModel.molecules: molec_obj.name = str(translator[species]) else: molec_obj.name = translator[species].str2() self.bngModel.add_molecule(molec_obj) annotationInfo["species"] = speciesAnnotationInfo self.speciesMemory = [] return ( list(set(moleculesText)), speciesText, observablesText, speciesTranslationDict, observablesDict, annotationInfo, )
[docs] def getInitialAssignments( self, translator, param, zparam, molecules, initialConditions ): """ process the get initial assignments section. This can be used to initialize parameters or species, so we have to account for both checking both arrays """ param2 = param zparam2 = zparam initialConditions2 = initialConditions pparam = {} initCondMap = {} for element in param: pparam[element.split(" ")[0]] = (element.split(" ")[1], None) for element in zparam: pparam[element] = (0, None) for species in self.model.getListOfSpecies(): tmp = self.getRawSpecies(species) name = tmp["returnID"] constant = ( "$" if (species.getConstant() or species.getBoundaryCondition()) else "" ) if name in translator: if self.noCompartment: extendedStr = "{1}{0}".format(translator[name], constant) else: extendedStr = "@{0}:{2}{1}".format( species.getCompartment(), translator[name], constant ) else: std_name = standardizeName(tmp["name"]) + "()" # ASS - deal with no compartment case if self.noCompartment: extendedStr = "{1}{0}".format(std_name, constant) else: extendedStr = "@{0}:{2}{1}".format( tmp["compartment"], std_name, constant ) initConc = ( species.getInitialConcentration() if species.isSetInitialConcentration() else species.getInitialAmount() ) pparam[species.getId()] = (initConc, extendedStr) from copy import copy for initialAssignment in self.model.getListOfInitialAssignments(): symbol = initialAssignment.getSymbol() math = libsbml.formulaToString(initialAssignment.getMath()) for element in pparam: if math is None: logMess( "ERROR:SIM210", f"Initial assignment for '{element}' has no math defined", ) raise TranslationException( f"ERROR:SIM210: Initial assignment for '{element}' has no math defined" ) if element in math: math = re.sub( r"(\W|^)({0})(\W|$)".format(element), r"\1({0})\3".format(pparam[element][0]), math, ) # removing non bngl math elements for their equivalents math = writer.bnglFunction(math, "", []).split(" = ")[1] param2 = [x for x in param if "{0} ".format(symbol) not in x] zparam2 = [x for x in zparam if symbol != x] """ if (len(param2) != len(param)) or (len(zparam) != len(zparam2)): param2.append('{0} {1}'.format(symbol,math)) param = param2 zparam = zparam2 """ try: if pparam[symbol][1] == None: param2.append("{0} {1}".format(symbol, math)) param = param2 zparam = zparam2 else: initialConditions2 = [ x for x in initialConditions if "#{0}".format(symbol) not in x ] initialConditions2.append( "{0} {1} #{2}".format(pparam[symbol][1], math, symbol) ) initCondMap[ "{0} {1} #{2}".format(pparam[symbol][1], math, symbol) ] = symbol initialConditions = initialConditions2 except: logMess( "WARNING:RATE006", f"Can't find what {symbol} is referring to, some initial conditions might be incorrect.", ) continue return param, zparam, initialConditions, initCondMap
[docs] def getSpeciesAnnotation(self): if self.speciesAnnotation: return self.speciesAnnotation self.speciesAnnotation = defaultdict(list) for species in self.model.getListOfSpecies(): rawSpecies = self.getRawSpecies(species, logEntries=False) annotationXML = species.getAnnotation() lista = libsbml.CVTermList() libsbml.RDFAnnotationParser.parseRDFAnnotation(annotationXML, lista) if lista.getSize() == 0: self.speciesAnnotation[rawSpecies["returnID"]] = [] else: for idx in range(lista.getSize()): self.speciesAnnotation[rawSpecies["returnID"]].append( lista.get(idx).getResources() ) return self.speciesAnnotation
[docs] def getFullAnnotation(self): if self.speciesAnnotationDict: return self.speciesAnnotationDict self.speciesAnnotationDict = defaultdict(lambda: defaultdict(list)) for species in self.model.getListOfSpecies(): rawSpecies = self.getRawSpecies(species, logEntries=False) annotationXML = species.getAnnotation() lista = libsbml.CVTermList() libsbml.RDFAnnotationParser.parseRDFAnnotation(annotationXML, lista) if lista.getSize() == 0: continue else: for idx in range(lista.getSize()): for idx2 in range(0, lista.get(idx).getResources().getLength()): resource = lista.get(idx).getResources().getValue(idx2) qualifierType = lista.get(idx).getQualifierType() qualifierDescription = ( bioqual[lista.get(idx).getBiologicalQualifierType()] if qualifierType else modqual[lista.get(idx).getModelQualifierType()] ) self.speciesAnnotationDict[rawSpecies["returnID"]][ qualifierDescription ].append(resource) return self.speciesAnnotationDict
[docs] def getModelAnnotation(self): modelAnnotation = [] annotationXML = self.model.getAnnotation() lista = libsbml.CVTermList() libsbml.RDFAnnotationParser.parseRDFAnnotation(annotationXML, lista) if lista.getSize() == 0: modelAnnotations = [] else: tempDict = {} for element in [2, 3, 4, 5, 6]: if lista.get(element) == None: continue tempDict[lista.get(element).getBiologicalQualifierType()] = lista.get( element ) if 3 in tempDict: modelAnnotation = tempDict[3].getResources() elif 0 in tempDict and ( "GO" in tempDict[0].getResources().getValue(1) or "kegg" in tempDict[0].getResources().getValue(1) ): modelAnnotation = tempDict[0].getResources() elif 5 in tempDict: modelAnnotation = tempDict[5].getResources() else: if lista.get(3) != None and ( "GO" in lista.get(3).getResources().getValue(0) or "kegg" in lista.get(3).getResources().getValue(0) ): modelAnnotation = lista.get(3).getResources() elif lista.get(4) != None and ( "GO" in lista.get(4).getResources().getValue(0) or "kegg" in lista.get(4).getResources().getValue(0) ): modelAnnotation = lista.get(4).getResources() elif lista.get(5) != None and ( "GO" in lista.get(5).getResources().getValue(0) or "kegg" in lista.get(5).getResources().getValue(0) ): modelAnnotation = lista.get(5).getResources() else: if lista.get(3) != None and ( "reactome" in lista.get(3).getResources().getValue(0) ): modelAnnotation = lista.get(3).getResources() elif lista.get(4) != None and ( "reactome" in lista.get(4).getResources().getValue(0) ): modelAnnotation = lista.get(4).getResources() elif lista.get(5) != None and ( "reactome" in lista.get(5).getResources().getValue(0) ): modelAnnotation = lista.get(5).getResources() elif lista.get(2) != None: modelAnnotation = lista.get(2).getResources() return modelAnnotation
[docs] def getSpeciesInfo(self, name): return self.getRawSpecies(self.model.getSpecies(name))
[docs] def writeLog(self, translator): rawSpecies = [self.getRawSpecies(x) for x in self.model.getListOfSpecies()] log["species"].extend([x[0] for x in rawSpecies if x[0] not in translator]) logString = "" # species stuff if len(log["species"]) > 0: logString += "Species we couldn't recognize:\n" for element in log["species"]: logString += "\t%s\n" % element if len(log["reactions"]) > 0: logString += "Reactions we couldn't infer more about due to \ insufficient information:" for element in log["reactions"]: logString += "\t%s + %s -> %s\n" % ( element[0][0], element[0][1], element[1], ) return logString
[docs] def getStandardName(self, name): if name in self.speciesDictionary: return self.speciesDictionary[name] return name
[docs]def standardizeName(name): """ Remove stuff not used by bngl """ name2 = name sbml2BnglTranslationDict = { "^": "", "'": "", "*": "m", " ": "_", "#": "sh", ":": "_", "α": "a", "β": "b", "γ": "g", " ": "", "+": "pl", "/": "_", ":": "_", "-": "_", ".": "_", "?": "unkn", ",": "_", "(": "", ")": "", "[": "", "]": "", # "(": "__", # ")": "__", # "[": "__", # "]": "__", ">": "_", "<": "_", } for element in sbml2BnglTranslationDict: name = name.replace(element, sbml2BnglTranslationDict[element]) name = re.sub("[\W]", "", name) return name