models.py

# -*- coding: utf-8 *-*
import subprocess
import xml.dom.minidom as mdom
import collectiontools.reactions as reactions
sbmlStam = """
<sbml xmlns="http://www.sbml.org/sbml/level2" level="2" version="1" xmlns:html="http://www.w3.org/1999/xhtml">
<model id="%(id)s" name="%(name)s">
<listOfUnitDefinitions>
    <unitDefinition id="mmol_per_gDW_per_hr">
        <listOfUnits>
            <unit kind="mole" scale="-3"/>
            <unit kind="gram" exponent="-1"/>
            <unit kind="second" multiplier=".00027777" exponent="-1"/>
        </listOfUnits>
    </unitDefinition>
</listOfUnitDefinitions>
<listOfCompartments>
<compartment id="Extracellular"/>
<compartment id="Cytosol" outside="Extracellular"/>
</listOfCompartments>
<listOfSpecies>
</listOfSpecies>
<listOfReactions>
</listOfReactions>
</model></sbml>
"""
kineticLawStam = """
<kineticLaw>
    <math xmlns="http://www.w3.org/1998/Math/MathML">
            <ci> FLUX_VALUE </ci>
    </math>
    <listOfParameters>
        <parameter id="LOWER_BOUND" value="%(lowerbound)s" units="mmol_per_gDW_per_hr"/>
        <parameter id="UPPER_BOUND" value="%(upperbound)s" units="mmol_per_gDW_per_hr"/>
        <parameter id="OBJECTIVE_COEFFICIENT" value="0.0"/>
        <parameter id="FLUX_VALUE" value="0.0" units="mmol_per_gDW_per_hr"/>
    </listOfParameters>
</kineticLaw>
"""
def readSurreyModel(surreyFile):
    surreyFile = open(surreyFile,"r")
    model = list()
    for line in surreyFile.readlines():
        newReaction = _parseReactionLine(line)
        model.append(newReaction)
    surreyFile.close()
    return model
def readSbml(sbmlFile):
    outputName = sbmlFile[:sbmlFile.rfind(".")]
    outputName += ".surrey"
    args = ["python", "parseSbml.py", sbmlFile, outputName]
    subprocess.call(args)
    model = readSurreyModel(outputName)
    return model
def writeSurreyModel(reactions,targetfile):
    with open(targetfile,"w") as f:
        for r in reactions:
            f.write(str(r) + "\n")
def writeMinimalSBML(reactions,targetFile,modelName):
    with open(targetFile,"w") as tFile:
        xmlDoc = _assembleDocument(reactions,modelName)
        xmlString = xmlDoc.toprettyxml()
        tFile.write(xmlString)
        
def _assembleDocument(reactions,modelName):
    allComps = getAllInvolvedCompounds(reactions)
    xmlDoc = _sbmlDocStam(modelName)
    for comp in allComps:
        _compNode(comp,xmlDoc)
    for reaction in reactions:
        _reactNode(reaction,xmlDoc)
    print(xmlDoc)
    return xmlDoc
def _sbmlDocStam(modelName):
    stam = sbmlStam % {"name":modelName,"id":modelName}
    xmlDoc = mdom.parseString(stam)
    return xmlDoc
def _reactNode(reaction,xmlDoc):
    node = xmlDoc.createElement("reaction")
    reactList = xmlDoc.getElementsByTagName("listOfReactions")[0]
    _setAttribute(node,"name",reaction.name)
    _setAttribute(node,"id",reaction.name)
    rev = "True"
    if reaction.lowerbound == 0:
        rev = "False"
    _setAttribute(node,"reversible",rev)
    _addReactants(node,reaction,xmlDoc)
    _decorateReactNode(node,xmlDoc,reaction)
    reactList.appendChild(node)
def _decorateReactNode(reactNode,xmlDoc,reaction):
    notes = xmlDoc.createElement("notes")
    reactNode.appendChild(notes)
    kineticLaw = kineticLawStam % {"lowerbound":reaction.lowerbound,"upperbound":reaction.upperbound}
    kineticLaw = mdom.parseString(kineticLaw)
    reactNode.appendChild(kineticLaw.getElementsByTagName("kineticLaw")[0])
def _addReactants(node,reaction,xmlDoc):
    substrates = reaction.getSubstrates()
    products = reaction.getProducts()
    subsList = xmlDoc.createElement("listOfReactants")
    prodList = xmlDoc.createElement("listOfProducts")
    for sub in substrates:
        stoichiometry = str(reaction._substrates[sub])
        sRef = xmlDoc.createElement("speciesReference")
        _setAttribute(sRef,"species",sub)
        _setAttribute(sRef,"stoichiometry",stoichiometry)
        subsList.appendChild(sRef)
    node.appendChild(subsList)
    for prod in products:
        stoichiometry = str(reaction._products[prod])
        sRef = xmlDoc.createElement("speciesReference")
        _setAttribute(sRef,"species",prod)
        _setAttribute(sRef,"stoichiometry",stoichiometry)
        prodList.appendChild(sRef)
    node.appendChild(prodList)
def _compNode(compName,xmlDoc):
    node = xmlDoc.createElement("species")
    _setAttribute(node,"name",compName)
    _setAttribute(node,"id",compName)     
    #_setAttribute(node,"compartment",compartment)
    boundraryCond = isExternal(compName)
    _setAttribute(node,"boundaryCondition",str(boundraryCond))
    _setAttribute(node,"charge","0")
    speciesList = xmlDoc.getElementsByTagName("listOfSpecies")[0]
    speciesList.appendChild(node)
def _setAttribute(node,attribute,value):
    #node.attributes.setNamedItem(attribute)
    #node.attributes.getNamedItem(attribute).nodeValue = value
    node.setAttribute(attribute,value)
def writeListOfExternals(model,targetFile):
    with open(targetFile,"w") as f:
        f.write("\\")
        allCompounds = set()
        for r in model:
            allCompounds = allCompounds.union(r.getCompounds())
        for c in allCompounds:
            if isExternal(c):
                f.write(c+" ")
def producingReaction(reactions,metabolite):
    for r in reactions:
        if r.isProduct(metabolite):
            return r
    raise ValueError("No reaction of origin for:"+metabolite )
def consumingReaction(reactions,metabolite):
    for r in reactions:
        if r.isSubstrate(metabolite):
            return r
    raise ValueError("No reaction consuming:"+metabolite )
def getReactionByName(reactions,name):
    for r in reactions:
        if r.name == name:
            return r
    raise Exception("could not find reaction:"+name)
def isExternal(compound):
    if compound[-4:] == '_ext':
        return True
    else:
        return False
def _parseReactionLine(line):
    terms = line.split('\t')
    reaction = reactions.Reaction(terms[0])
    loBound = terms[2]
    loBound = loBound.strip()
    upBound = terms[3]
    upBound = upBound.strip()
    reaction.upperbound = upBound
    reaction.lowerbound = loBound
    equationLine = terms[1].split('=')
    substrates = equationLine[0].split('+')
    products = equationLine[1].split('+')
    for s in substrates:
        s = s.strip()
        s = s.split()
        if len(s) > 1:
            reaction.addSubstrate(s[1], s[0])
        else:
            reaction.addSubstrate(s[0], 1.0)
    for s in products:
        s = s.strip()
        s = s.split()
        if len(s) > 1:
            reaction.addProduct(s[1], s[0])
        else:
            reaction.addProduct(s[0], 1.0)
    return reaction
def getAllExternals(model):
    allComps = getAllInvolvedCompounds(model)
    externals = list()
    for comp in allComps:
        if isExternal(comp):
            externals.append(comp)
    return externals
def getTransportReaction(model,metabolite):
    for r in model:
        if r.isInvolved(metabolite):
            return r
def getAllInvolvedCompounds(reactions):
    allComps = set()
    for r in reactions:
        allComps = allComps.union(r.getCompounds())
    return allComps

def test_writeMinimalSBML():
    testPath = r"sbmlmodel.sbml"
    testModel = readSbml(testPath)
    writeMinimalSBML(testModel,"testMinimalSbml.sbml","testModel")
if __name__ == "__main__":
    test_writeMinimalSBML()