#!/usr/bin/python """ ########################################### cbraPytool A Python tool for constraint-based metabolic model reconstruction analysis (pre alpha verion 0.0) ########################################### It can: - Read SBML format model into CbModel format in python - Manual constraint-based model modification - Flux balance analysis - optimize the CbModel by -- standard FBA (maximizing the flux for biomass formation, i.e growth rate) -- MOMA (minimisation of metabolic adjustment) It depends on other python packages: - numpy - libsbml: to read model in SBML format - solver packages: -- lpSolve (for standard FBA) -- openopt (required for MOMA only) To run the example of simulating iND750 model: cbraPytool [-e] [sbmlFile] """ import sys from time import time from lp_solve import * from lp_maker import * from libsbml import * sbmlR = SBMLReader() from numpy import * from scikits.openopt import LP, MILP, NLP, QP import csv from copy import deepcopy # #-------------------------------------------------------------------------------------------------- # Solution for constraint-based model optimization problems class solution: def __init__(self, obj=None, x=None, dual=None, feasible=None, msg=None, gap=None, solver=None, elapsed=None, solOrigin=None): self.obj = obj; self.x = x; self.dual = dual self.feasible = feasible self.msg = msg self.gap = gap self.solver = solver self.elapsed = elapsed self.solOrigin = solOrigin def copy(self): return solution(obj=self.obj, x=self.x.copy(), dual=self.dual.copy(), feasible=self.feasible, msg=self.msg, gap=self.gap, solver=self.solver, elapsed=self.elapsed, solOrigin=self.solOrigin) #______________________________________________________________________________ # P-Graphs are used to describe network structure besides the stoichiometric matrix # required for the constraint-based models class PGraph: """A PGraph, or a Process Graph is a kind of bipartite graph, whose vertices can be divided into two disjoint sets M (material or metabolite set) and O (operation or reaction set) such that every edge connects a vertex in M to one in O; that is, M and O are independent sets. The constructor call is something like: pg = PGraph({'M':{('A',0):0,('B',0):0,('C',0):0}, 'O':{'R1':{'In':{('A',0):1,('B',0):1},'Out':{('C',0):1}}},directed=True) pg = PGraph({'M':{('A',0):0,('B',0):0,('C',0):0}, 'O':{'R1':{'In':{('A',0):1,('B',0):1},'Out':{('C',0):1}, 'Mod':{('E1',0)}, 'Wt':1}},directed=True) this makes a graph with 3 nodes, A, B, and C, with an edge of length 1 from A to B, and an edge of length 2 from A to C. You can also do: g = Graph({'A': {'B': 1, 'C': 2}, directed=False) This makes an undirected graph, so inverse links are also added. The graph stays undirected; if you add more links with g.connect('B', 'C', 3), then inverse link is also added. You can use g.nodes() to get a list of nodes, g.get('A') to get a dict of links out of A, and g.get('A', 'B') to get the length of the link from A to B. 'Lengths' can actually be any object at all, and nodes can be any hashable object.""" def __init__(self, dict=None, Mok=False, Ook=False, directed=True, oleadtom = None): self.dict = dict or {} self.directed = directed self.Mok = Mok self.Ook = Ook if directed: for k in self.dict['O']: self.dict['O'][k].setdefault('Rev', False) else: for k in self.dict['O']: self.dict['O'][k]['Rev'] = True self.oleadtom = oleadtom or {} def copy(self): return PGraph({'M':self.dict['M'].copy(), 'O': self.dict['O'].copy()}, Mok=self.Mok, Ook=self.Ook, directed=self.directed, oleadtom = self.oleadtom.copy()) def checkM(self): M = self.dict['M'] O = self.dict['O'] keys = set() keysM = M.keys() for k in O: keys.update(O[k]['In'].keys()) keys.update(O[k]['Out'].keys()) Mc = {} for k in keys: if k in keysM: Mc[k] = M[k] else: Mc[k] = 1 self.dict['M'] = Mc self.Mok = True def make_undirected(self): "Make a digraph into an undirected graph by adding symmetric edges." for k in self.dict['O']: self.dict['O'][k]['Rev'] = True def reverse1(self, A): """Update an operating unit A (key in the O-type dictionary) in a P-Graph. A should be a string. The reverse units is named '-'+A """ " Check existence of A in the dictionary" if A in self.dict['O']: Ar = {} Ar['In'] = self.dict['O'][A]['Out'] Ar['Out'] = self.dict['O'][A]['In'] exist = True i = 0 while exist: ArKey = A+'-' if ArKey in self.dict['O']: i = i+1 ArKey = ArKey + str(i) else: exist = False self.dict['O'][ArKey] = Ar return ArKey else: return None def reverse(self, A=None): """Update a set of operating units A (set of keys in the O-type dictionary) in a P-Graph. Keys should be of string type. """ if A == None: A = self.dict['O'].keys() Ar = [] for k in A: Ar.append(self.reverse1(k)) return Ar def cost(self, a=None): """Return an operation cost or sum of all the operation costs""" a = a or None if a != None: if type(a) is set or type(a) is list or type(a) is tuple: c = 0 for a1 in a: c = c + self.dict['O'][a1].setdefault('Cost',1) return c else: return self.dict['O'][a].setdefault('Cost', 1) c = 0 for a in self.dict['O']: c = c + self.dict['O'][a].setdefault('Cost', 1) return c "------ start of pgraph methods" def show(self, A): """ Print out a list of operations(s) in the graph A should be a dictionary or a list/set of keys. """ print '== Reaction list ==' for o in A: print '%s: ' %o print ' %s' %(self.dict['O'][o]['In']) if self.dict['O'][o]['Rev']: direction = '<->' else: direction = '->' print '%s %s' %(direction, self.dict['O'][o]['Out']) print '' def nodes(self): """Return a list of nodes in the p-graph, including both M-type and O-type nodes.""" Nds = self.dict['M'] Nds.update(self.dict['O']) return Nds def inputs(self, A=None): """Return a list of materials appear as inputs in sets of keys for a list of operating units""" keysIn = set() O = self.dict['O'] if A == None: A = O.keys() for key in A: keysIn.update(O[key]['In'].keys()) if O[key].has_key('Rev') and O[key]['Rev']: keysIn.update(O[key]['Out'].keys()) Inputs = {} for key in keysIn: Inputs[key] = self.dict['M'][key] return Inputs def outputs(self, A=None): "Return a list of materials appear as outputs in operating units" keysOut = set() O = self.dict['O'] if A == None: A = O.keys() for key in A: keysOut.update(O[key]['Out'].keys()) if O[key].has_key('Rev') and O[key]['Rev']: keysOut.update(O[key]['In'].keys()) Outputs = {} for key in keysOut: Outputs[key] = self.dict['M'][key] return Outputs def structmap(self, nodes, direction='b', sufficient=False): """ Return a a list of """ nodeO = {} nodeM = {} """ Separate two types of nodes """ if type(nodes) is dict: for k in nodes: if type(nodes[k]) is dict: nodeO[k] = nodes[k] else: nodeM[k] = nodes[k] else: for k in nodes: if k in self.dict['O'].keys(): nodeO[k] = self.dict['O'][k] else: nodeM[k] = self.dict['M'][k] nodesmap = {} "Dealing with structure mapping for operating unit nodes" keys = set() for key in nodeO: node = nodeO[key] if direction=='b': keys.update(node['In'].keys()) keys.update(node['Out'].keys()) elif direction=='-': keys.update(node['In'].keys()) if node.has_key('Rev') and node['Rev']: keys.update(node['Out'].keys()) else: keys.update(node['Out'].keys()) if node.has_key('Rev') and node['Rev']: keys.update(node['In'].keys()) for key in keys: nodesmap[key] = self.dict['M'][key] "Dealing with structure mapping for material nodes" keys = set(nodeM.keys()) O = self.dict['O'] for key in O: if direction=='b': O1keys = O[key]['In'].keys()+O[key]['Out'].keys() elif direction=='-': O1keys = O[key]['Out'].keys() if O[key].has_key('Rev') and O[key]['Rev']: O1keys = O1keys + O[key]['In'].keys() else: O1keys = O[key]['In'].keys() O1keys1 = O1keys; O1keys2 = [] if O[key].has_key('Rev') and O[key]['Rev']: O1keys = O1keys + O[key]['Out'].keys() if sufficient: O1keys2 = O[key]['Out'].keys() if direction == '+' and sufficient: if set(O1keys1).issubset(keys): nodesmap[key] = O[key] elif O[key].has_key('Rev') and O[key]['Rev']: if set(O1keys2).issubset(keys): nodesmap[key] = O[key] elif len(set(O1keys).intersection(keys)) > 0: nodesmap[key] = O[key] return nodesmap def links(self, vertex, direction='b'): """ Return the links (a node dictionary) that are connected to the given vertex (a key in the pgraph M type or O type dictionary)""" O = self.dict['O'] links = {} if vertex in O.keys(): if direction == 'b': links.update(O[vertex]['In']) links.update(O[vertex]['Out']) elif direction == '-': links.update(O[vertex]['In']) if O[vertex].has_key('Rev') and O[vertex]['Rev']: links.update(O[vertex]['Out']) elif direction == '+': links.update(O[vertex]['Out']) if O[vertex].has_key('Rev') and O[vertex]['Rev']: links.update(O[vertex]['In']) else: print "Direction can only be one of 'b','-' or '+'" links = None elif vertex in self.dict['M'].keys(): for key in O: if direction == 'b': if (vertex in O[key]['In'].keys()) or (vertex in O[key]['Out'].keys()): links[key] = O[key] elif direction == '-': if vertex in O[key]['Out'].keys(): links[key] = O[key] elif O[key].has_key('Rev') and O[key]['Rev']: if vertex in O[key]['In'].keys(): links[key] = O[key] elif direction == '+': if vertex in O[key]['In'].keys(): links[key] = O[key] elif O[key].has_key('Rev') and O[key]['Rev']: if vertex in O[key]['Out'].keys(): links[key] = O[key] else: print "Direction can only be one of 'b','-' or '+'" links = None else: links = None return links def dlinks(self, vertex, direction='b'): """ Return the number of links that are connected to the given vertex (a key in the pgraph M type or O type dictionary)""" O = self.dict['O'] d = 0 if vertex in O.keys(): if direction == 'b': d = len(O[vertex]['In']) + len(O[vertex]['Out']) elif direction == '-': d = len(O[vertex]['In']) if O[vertex].has_key('Rev') and O[vertex]['Rev']: d = d + len(O[vertex]['Out']) else: d = len(O[vertex]['Out']) if O[vertex].has_key('Rev') and O[vertex]['Rev']: d = d + len(O[vertex]['In']) elif vertex in self.dict['M'].keys(): for key in O: if direction == 'b': if vertex in O[key]['In'].keys(): d = d + 1 elif vertex in O[key]['Out'].keys(): d = d + 1 elif direction == '-': if vertex in O[key]['Out'].keys(): d = d + 1 elif O[key].has_key('Rev') and O[key]['Rev']: if vertex in O[key]['In'].keys(): d = d + 1 elif direction == '+': if vertex in O[key]['In'].keys(): d = d + 1 elif O[key].has_key('Rev') and O[key]['Rev']: if vertex in O[key]['Out'].keys(): d = d + 1 else: d = None return d def forward_reachable(self, start, goal=None): goal = goal or set() Okcur = set() Scur = set(start) while True: Oact = self.structmap(Scur, '+', sufficient=True) Mact = self.outputs(Oact) Scur = Scur.union(set(Mact)) Okact = set(Oact.keys()) if Okact == Okcur: break else: Ocur = Oact Okcur = Okact Mcur = {} for m in Scur: Mcur[m] = self.dict['M'].setdefault(m, 1) g_reachable = PGraph({'M':Mcur, 'O':Ocur}) if goal is None: return g_reachable else: goal_unreachable = set(goal) - set(Scur) return g_reachable, goal_unreachable # # ------------------------------------------------------------------------------------------------" # Constraint-based optimization for metabolic model reconstruction analysis class CbModel(PGraph): "Constraint based optimization problem for flux balance analysis." """ Note: objC - in object of cbmod 'Objc' - key in dict['O'] """ def __init__(self, dict={'M':{},'O':{},'S':{},'C':{}}, Mok=False, Ook=False, directed=True, oleadtom = None, Aeq=None, beq=None, A=None, b=None, lb=None, ub=None, intVars=None, lstM=None, lstR=None, lstEnz=None, lstMBound=None, objC=None, objOpt='max', flux=None, ninf = 1000.0, option = None): PGraph.__init__(self, dict=dict, Mok=Mok, Ook=Ook, directed=directed, oleadtom = oleadtom) self.Aeq = array(Aeq); self.beq = array(beq); self.A=array(A); self.b = array(b); self.lb = array(lb); self.ub = array(ub); self.intVars = intVars or []; self.lstM = array(lstM); self.lstR = array(lstR); self.lstEnz = array(lstEnz); self.lstMBound = array(lstMBound); self.objC = array(objC); self.flux = array(flux); self.objOpt = objOpt self.ninf = ninf or 1000.0 self.option = option or {} def copy(self): return CbModel({'M':self.dict['M'].copy(), 'O': self.dict['O'].copy(),'C':self.dict['C'].copy(), 'S': self.dict['S'].copy()}, self.Mok, self.Ook, self.directed, self.oleadtom.copy(), self.Aeq.copy(), self.beq.copy(), self.A.copy(), self.b.copy(), self.lb.copy(), self.ub.copy(), self.intVars, self.lstM.copy(), self.lstR.copy(), self.lstEnz.copy(), self.lstMBound.copy(), self.objC.copy(), self.objOpt, self.flux.copy(), self.ninf, self.option.copy()) def readSbmlCbModel(self, sbmlFile="Sc_iND750_GlcMM.xml", lubPreset=True, ninf=1000.0, objOpt='max'): [O, specs, comps] = self.readSbmlModel(sbmlFile=sbmlFile, lubPreset=lubPreset, ninf=ninf) self.getCbModel(O=O, S=specs, C=comps, lubPreset=lubPreset) return self def readSbmlModel(self, sbmlFile="Sc_iND750_GlcMM.xml", lubPreset=False, ninf=1000.0): doc = sbmlR.readSBML(sbmlFile) print doc.getNumErrors() print doc.getVersion() model = doc.getModel() # Read stochiometry matrix [O, specs, comps] = getSbmlModel(model, useKeggid = False, useCompart=False) print comps return O, specs, comps def getCbModel(self, O=None, S=None, C=None, lstR=None, lstMBound=None, lubPreset=True, ninf=None, objOpt=None): # lstR - a list of reaction IDs O = O or self.dict['O'] specs = S; if specs == None and 'S' in self.dict: specs = self.dict['S'] comps = C; if comps == None and 'C' in self.dict: comps = self.dict['C'] if lstR == None: lstR = self.lstR lstR = O.keys() ninf = ninf or self.ninf if objOpt == None: objOpt = self.objOpt lstEnz = [O[r1]['Enz'] for r1 in lstR] lstEnz = array(lstEnz) self.dict = {'M':{},'O':O,'S':specs,'C':comps} self.checkM() # Get dict['M'] which are metabolites involved in dict['O'] only nM = len(self.dict['M']) nR = len(self.dict['O']) lstM = self.dict['M'].keys() ind = arange(nM) STmat = zeros((nM, nR)) # Stoichiometric meatrix for all reactions in O for i1 in range(nR): r1 = O[lstR[i1]] for s1 in r1['In']: idx = ind.compress(array(lstM)==s1) STmat[idx, i1] = -r1['In'][s1] for s1 in r1['Out']: idx = ind.compress(array(lstM)==s1) STmat[idx, i1] = +r1['Out'][s1] lstM = array(lstM); lstR = array(lstR); objC = zeros(nR) for i1 in arange(nR): if 'Objc' in O[lstR[i1]]: objC[i1] = O[lstR[i1]]['Objc'] or 0 if sum(objC != 0) > 0: print ' - sum(objC != 0): ', sum(objC != 0) keyObj = lstR[objC!= 0][0] objM = O[keyObj]['Out'].keys() objMin = O[keyObj]['In'].keys() else: keyObj = '' objM = [] objMin = [] idxMobj = [arange(nM).compress(lstM==i)[0] for i in objM] idxMobjin = [arange(nM).compress(lstM==i)[0] for i in objMin] # get boundary condition: read all species whose BoundaryCondition is True if lstMBound is None: #lstBound = [specs[m1]['Bound'] for m1 in lstM] lstBound =[] for m1 in lstM: try: bnd = specs[m1]['Bound'] except: bnd = False lstBound.append(bnd) else: lstBound = [b1 in lstMBound for b1 in lstM] # non existing metabolite can't be bound lstInner = [not lstBound[i1] for i1 in range(nM)] idxBound = list(ind.compress(lstBound)) #idxInner = list(set(range(nM)) - set(idxBound) - set(idxMobj)) idxInner = list(ind.compress(lstInner)) lstMBound = lstM[idxBound] Aeq = STmat[idxInner,:] beq = zeros(len(Aeq)) lstM = lstM[idxInner] # initial setting for contraints for bound lb = -ninf * ones(nR); # unbounded ub = ninf * ones(nR); # unbounded flx = zeros(nR) for i1 in range(nR): r = O[lstR[i1]] if not r['Rev']: lb[i1] = 0 if lubPreset: if 'LB' in r: lb[i1] = r['LB'] if 'UB' in r: ub[i1] = r['UB'] if 'Flux' in r: if type(r['Flux']) in [float, int]: flx[i1] = r['Flux'] self.Aeq = Aeq; self.beq=beq; self.lb=lb; self.ub=ub; self.lstM = lstM; self.lstR=lstR; self.lstEnz=lstEnz; self.lstMBound=lstMBound self.objC=objC; self.flux=flx; self.ninf = max(abs(ub)) return self def getMBound(self): # get boundary condition: read all species whose BoundaryCondition is True if 'S' in self.dict: specs = self.dict['S'] lstBound = [specs[m1]['Bound'] or False for m1 in self.lstM] else: return array([]) idxBound = list(arange(len(self.lstM)).compress(lstBound)) lstMBound = lstM[idxBound] return lstMBound def resetBound(self, lubPreset=True, controlUptake=True, ninf = None): """ If lubPreset is True, the lower and upper bound will be set to as defined in self.dict['O'] If controlUptakself.e, all uptake rate will be set to zeros """ # initial setting for contraints for bound ninf = ninf or self.ninf nR = len(self.lstR) lb = -ninf * ones(nR); # unbounded ub = ninf * ones(nR); # unbounded notR=0 for i1 in range(nR): r = self.dict['O'][self.lstR[i1]] if lubPreset: if 'LB' in r: lb[i1] = r['LB'] if 'UB' in r: ub[i1] = r['UB'] if not r['Rev']: lb[i1] = 0 if controlUptake: for b1 in self.lstMBound: rs = self.links(b1, '+') if rs is None: continue for k1 in rs: if not rs[k1]['Rev']: r1 = k1; useUB = True; ub[self.lstR == r1] = 0 continue if b1 not in rs[k1]['In']: r1 = k1; useUB = False; lb[self.lstR == r1] = 0 #break else: r1 = k1; useUB = True ub[self.lstR == r1] = 0 if False: if useUB: ub[self.lstR == r1] = 0 else: lb[self.lstR == r1] = 0 self.lb = lb; self.ub = ub self.ninf = max(abs(ub)) return self def setUptakeBound(self, uptakeBound, printBound=False): """ update bounds for exchange reactions based on predefined uptake rate uptakeBound: a dictionary object storing the maximum uptake rate info lstMBound: a list array of boundary metabolites """ lb1 = self.lb.copy(); ub1 = self.ub.copy() lstR = self.lstR; for b1 in self.lstMBound: rs = self.links(b1, '+') if rs is None: continue if len(rs) > 1: print 'Bound', b1, ' has multiple reactions linked!' return(self) for k1 in rs: if not rs[k1]['Rev']: r1 = k1; useUB = True; break if b1 not in rs[k1]['In']: r1 = k1; useUB = False; else: r1 = k1; useUB = True if b1 in uptakeBound: bd1 = abs(uptakeBound[b1]) else: bd1 = 0 if useUB: ub1[lstR == r1] = bd1 else: lb1[lstR == r1] = -bd1 if printBound: print r1 + ': ' + str(self.dict['O'][r1]['Rev'])+str(self.lb[lstR==r1]) + str(self.ub[lstR==r1]) print '-----> ' + str(lb1[lstR == r1]) + str(ub1[lstR == r1]) self.lb = lb1; self.ub = ub1; return lb1, ub1 def getMidx(self, Ms): "Get the indices for metabolites in lstM" idxMobj = [arange(len(self.lstM)).compress(self.lstM==i)[0] for i in Ms] return idxMs def getRidx(self, Rs): "Get the indices for reactions in lstR" idxRs = [arange(len(self.lstR)).compress(self.lstR==r1)[0] for r1 in Rs] return idxRs def setObjR(self, objC, objOpt = 'max', updDict=True): """ objC: coefficent for objective function, or a set of indices/or reaction IDs whose reaction coefficient in objective function equal 1 """ nR = len(self.lstR) if type(objC) in [list, set]: objList = objC objC = zeros(nR) for i in objList: if i in range(nR): objC[i] = 1 elif type(i) is str: objC[self.lstR == i] = 1 objC = array(objC) if updDict: for i in range(nR): self.dict['O'][self.lstR[i]]['Objc'] = objC[i] self.objC = objC return self def getObjR(self): """Get the reaction IDs and relevant metabolites (reactants and products) with non-zero coefficient in objective function""" O = self.dict['O']; lstR = self.lstR; lstM = self.lstM; nR = len(lstR); nM = len(lstM); objC = self.objC if objC is None: objC = zeros(nR) for i1 in arange(nR): objC[i1] = O[lstR[i1]]['Objc'] or 0 keyObj = '' objM = [] objMin = [] if sum(objC != 0) > 0: print 'sum(objC != 0): ', sum(objC != 0) keyObj = lstR[objC!= 0] for keyObj1 in keyObj: objM = objM + O[keyObj1]['Out'].keys() objMin = objMin + O[keyObj1]['In'].keys() objR = keyObj objM = list(set(objM)); objMin = list(set(objMin)) return objR, objM, objMin def getMset(self, filter = None, ignoreCase=True): """ filter: dictionary containing attributes and corresponding values to filter out metabolite sets e.g. {'Comp': 'External', 'Bound': True, 'Annot': 'C05861', 'Contains':'CHEBI'} """ filter = filter or {} M = self.dict['M'].copy() S = self.dict['S'] filtered = [] for m1 in M: if m1 not in S: continue filterIn = False for f1 in filter: if f1 == 'Name': if 'Name' not in S[m1]: continue else: if ignoreCase: nam0 = kcpd.handle_synonyms(S[m1]['Name']) nam1 = kcpd.handle_synonyms(filter[f1]) filterIn = (nam1 == nam0) else: filterIn = (filter[f1] == S[m1]['Name']) if f1 == 'Annot': if 'Annot' not in S[m1]: continue else: if ignoreCase: filterIn = (filter[f1].lower() in S[m1]['Annot'].lower()) else: filterIn = (filter[f1] in S[m1]['Annot']) elif f1 == 'Comp': if 'Comp' not in S[m1]: continue else: filterIn = (filter[f1]==S[m1][f1]) if not filterIn: try: filterIn = (filter[f1] == self.dict['C'][S[m1][f1]]) except: filterIn = False elif f1 == 'Contains': if ignoreCase: filterIn = (str(filter[f1]).lower() in str(S[m1]).lower()) else: filterIn = (str(filter[f1]) in str(S[m1])) else: filterIn = (filter[f1] == S[m1][f1]) if filterIn: filtered.append(m1) return filtered def getRexc(self, by='Name'): """ Check the graph dict of the model and return a list of reaction ids who are for exchange. 1. Exchange reaction: Exchange reactions only have one non-zero (+1/-1) element in the corresponding column of the stoichiometric matrix. Uptake reactions are exchange reactions are exchange reactions with negative lower bounds. Check subsystem in dict['O']. """ O = self.dict['O'] lstRexc = [] if 'lstR' in self.__dict__.keys(): lstR = self.lstR else: lstR = O.keys() for r1 in lstR: key1 = O[r1]['Subs']; key1 = key1 or ''; key1 = key1.lower() key2 = O[r1]['Name']; key2 = key2 or ''; key2 = key2.lower() #print key1, 'exc' in key1, ' --- ', key2, 'exc' in key2 if ('exc' in key1) or ('exc' in key2): lstRexc.append(r1) return lstRexc def addBound(self, lstMBoundAdd, reversible=True, lubPreset=False): """ Add boundary metabolite and corresponding exchange reactions to relevant metabolite. Update pgraph dict and the constraint based model """ lstMBound = self.lstMBound.copy() lbmin = abs(min(self.lb)) ubmax = abs(max(self.ub)) for m1 in lstMBoundAdd: m1b = m1 + '_b' while m1b in self.dict['S']: print '[%s] Existing boundary metabolite ID : %s!!!', (m1, m1b) m1b = m1b+'_b' s1b = self.dict['S'][m1].copy() s1b['Bound'] = True self.dict['S'].update({m1b:s1b}) rbid = '_R_EX'+'_'+m1 while rbid in self.dict['O']: print '[%s] Existing boundary reaction ID : %s!!!', (m1, rbid) rbid = rbid+'_b' mn = self.dict['S'][m1]['Name'] if reversible: lb1 = -lbmin; ub1 = ubmax else: lb1 = 0; ub1 = ubmax rb = {'Name': '_'.join(['R', mn, 'exchange']), 'In':{m1: 1.0}, 'Out':{m1b: 1.0}, 'Enz':'', 'Rev': reversible, 'Objc':0.0, 'Flux':0.0, 'Subs':'_exchange','LB':lb1, 'UB':ub1} self.dict['O'].update({rbid:rb}) self.checkM() self.getCbModel(lubPreset=lubPreset) return self def getNumOfIsoEnz(self, KO, lstEnz=None): return None def getRwoEnz(self, KO, lstEnz=None, returnIndex=True, verbose=True): if lstEnz is None: lstEnz = self.lstEnz lstR = self.lstR nR = len(lstR) idxRmvd = set() if type(KO) is str: if ':' not in KO: KO = [KO] else: KO = KO.split(':') for KO1 in KO: lstRmvd1 = [(KO1.strip() != '' and self.disabledReactionEnz(enzModel=enz1, KO=KO1)) for enz1 in lstEnz.tolist()] idxRmvd1 = arange(nR).compress(lstRmvd1) idxRmvd = idxRmvd.union(set(idxRmvd1)) idxRmvd = list(idxRmvd) if verbose: print KO, ' :', len(idxRmvd), ' reactions disabled in model!' for i1 in idxRmvd: print self.dict['O'][lstR[i1]] if returnIndex: return idxRmvd else: return lstR[idxRmvd] def disabledReactionEnz(self, enzModel, KO=''): """check whether a certain KO (knockout) would disable the reaction given the annotated list of genes (enzModel) for that reaction input: enzModel- strings of enzymes that catalze the reaction, isoenzymes divided by '|', protein complexes enzymes divided by ':' KO - string of list of deleted enzymes divided by ':' output: True if the deletions of enzymes will diallowed the associated reactions; False otherwise. """ disabled = False if enzModel == '': return disabled if type(enzModel) is str: enzs = enzModel.split('|') else: enzs = enzModel if type(KO) is str: if ':' not in KO: KO = [KO] else: KO = KO.split(':') disabled = True for enz1 in enzs: if type(enz1) is str: enz1s = enz1.split(':') else: enz1s = set(enz1) koenz = set(KO).intersection(set(enz1s)) if len(koenz)==0: disabled = False break return disabled def setBinVars(self, binVars = None): binVars = binVars or [] intVars = list(set(self.intVars).union(set(binVars))) lb = self.lb.copy(); ub = self.ub.copy(); for b1 in binVars: lb[b1] = 0; ub[b1] = 1; self.lb = lb; self.ub = ub; self.intVars=intVars; return self def optCbModel(self, objC=None, Aeq=None, beq=None, A=None, b=None, lb=None, ub=None, objOpt=None, intVars=None, solver=None, getSolOrigin=True, option=None, cleanMemory=True): """ Constraint based model optimization """ if Aeq is None: Aeq = self.Aeq; if beq is None: beq = self.beq; if A is None: A = self.A; if b is None: b = self.b; if lb is None: lb = self.lb; if ub is None: ub = self.ub; if intVars is None: intVars = self.intVars if objC is None: objC = self.objC; sos = {} if 'sos' in self.__dict__.keys(): if sos is None: sos = self.sos objOpt = objOpt or self.objOpt solver = solver or 'lpsolve_orig' # or lpSolve option = option or {} method = 'LP'; if len(intVars) > 0: method = 'MILP'; print 'Method: ', method, ', ', objOpt, '; Solver: ', solver # For solvers not applied within openopt framework if solver == 'lpsolve_orig': if objOpt == 'min': signC = -1; else: signC = 1 intVars = list(intVars) if A.shape == (): va = Aeq.tolist(); vb = beq.tolist(); ve = zeros(len(va)); else: va = vstack((Aeq, A)).tolist(); vb = hstack((beq, b)).tolist() ve = zeros(len(va)); ve[len(Aeq):]=-1; ve=ve.tolist() vf = (signC*objC).tolist() vub = ub.tolist() vlb = lb.tolist() #print 'len of vf', len(vf), 'len va[0]', len(va), 'len va[0][0]', len(va[0]) "[v,x,duals] = lp_solve(f=vf,a=va,b=vb,e=ve,vlb=vlb,vub=vub,xint=intVars)" p = lp_maker(f=vf, a=va, b=vb, e=ve, vlb=vlb, vub=vub, xint=intVars) lpsolve('set_presolve', p, 1) # Presolve the rows: more stable solution for s1 in sos: lpsolve('add_sos', p, sos[s1]) for arg1 in option: exec("lpsolve('set_"+arg1+"', p, " + str(option[arg1]) + ")") timeStart = time() lpsolve('solve', p) msg = lpsolve('get_statustext', p, lpsolve('get_status', p)) gap = lpsolve('get_mip_gap', p, False) ff = lpsolve('get_objective', p) if ff != [] and ('INFEASIBLE' not in msg): ff = ff*signC presolve = lpsolve('get_presolve', p) if presolve == 0: xf = lpsolve('get_variables', p)[0] xf = array(xf) else: nx = len(objC) nrow = lpsolve('get_Norig_rows', p) xf = zeros(nx) for k in range(nx): xf[k] = lpsolve('get_var_primalresult', p, k+nrow+1) duals = xf feasible = True else: ff = nan; xf = nan*ones(len(objC)) duals = [] feasible = False solOrigin = p if getSolOrigin: solOrigin = solOrigin else: solOrigin = None elapsed = time() - timeStart; if cleanMemory: lpsolve('delete_lp', p) # Free the handle and its associated memory print 'Time Elapsed = ', elapsed print 'objFunValue: ', ff, '(feasible', feasible, ', ', msg, ')' print ' ' solCbM = solution(obj = ff, x = xf, dual = duals, feasible = feasible, msg = msg, gap = gap, solver = solver, elapsed = elapsed, solOrigin = solOrigin) return solCbM signC = 1; if method == 'LP': if objOpt == 'max': signC = -1; p = LP(f=signC*objC, A=None, Aeq=Aeq, b=None, beq=beq, lb=lb, ub=ub) elif method == 'MILP' and solver != 'lpsolve_orig': if objOpt == 'max': signC = -1; p = MILP(f=signC*objC, A=A, Aeq=Aeq, b=b, beq=beq, lb=lb, ub=ub, intVars=intVars) for arg1 in option: #print 'p.'+arg1 + '=' + str(option[arg1]) exec('p.'+arg1 + '=' + str(option[arg1])) sol = p.solve(solver) if not hasattr(sol, 'duals'): sol.duals = None solOrigin = None if getSolOrigin: solOrigin = sol solCbM = solution(obj = sol.ff*signC, x = sol.xf, dual = sol.duals, feasible = sol.isFeasible, msg = sol.msg, gap = sol.rf, solver = solver, elapsed = sol.elapsed['solver_time'], solOrigin = solOrigin) return solCbM def optTwoCbModels(self, mod1, flxWT=None, objC=None, Aeq=None, beq=None, A=None, b=None, lb=None, ub=None, objOpt=None, intVars=None, solver=None, getSolOrigin=True, option=None, cleanMemory=True): """ Constraint based model optimization optimize two CbModels Simultaneously solve two flux balance problems and minimize the difference between the two solutions self - the reference model for FBA analysis mod1 - the perturbed model from the reference model self flxWT - list flux values for the reference model """ solver = solver or 'ralg' # solver for QP problem only if objC is None: objC = self.objC if flxWT == None: solWT = self.optCbModel(objC, Aeq, beq, A, b, lb, ub, objOpt, intVars, None, getSolOrigin, option, cleanMemory) flxWT = solWT.x else: solWT = solution(obj = dot(flxWT, objC), x = flxWT, solver = '', elapsed = 0) if solver=='cvxopt_qp': #p1 = QP(H=diag(ones(mod1.Aeq.shape[1])), f=-2*flxWT, Aeq=mod1.Aeq, beq=mod1.beq, lb=mod1.lb, ub=mod1.ub) A=(ones(mod1.Aeq.shape[1])); print 'Rank A', rank(A) print 'Rank b', rank(array([1e12])) print '-2*flxWT', (-2*flxWT).shape p1 = QP(H=diag(ones(mod1.Aeq.shape[1])), f=-2*flxWT, A=(ones(mod1.Aeq.shape[1])), b=array([1e12]), Aeq=mod1.Aeq, beq=mod1.beq, lb=mod1.lb, ub=mod1.ub) sol = p1.solve('cvxopt_qp', iprint = 0) #sol = p1.solve('nlp:ralg') elif solver=='ralg': p1 = NLP(lambda x: ((flxWT-x)**2).sum(), x0=flxWT, Aeq=mod1.Aeq, b=None, beq=mod1.beq, lb=mod1.lb, ub=mod1.ub, iprint=0, maxIter = 1e5, contol=1e-5) p1.df = lambda x: 2*(x - flxWT) sol = p1.solve('ralg') solOrigin = None if getSolOrigin: solOrigin = sol if sol.isFeasible: obj = dot(sol.xf, objC) xf = sol.xf else: obj = nan; xf = nan*ones(len(objC)) solCbM = solution(obj = obj, x = xf, feasible = sol.isFeasible, msg = sol.msg, gap = sol.rf, solver = solver, elapsed = sol.elapsed['solver_time'], solOrigin = solOrigin) solCbM.xWT = flxWT return solCbM, solWT def deleteGenes(self, genes = ''): "A soft deletion of genes by fixing the flux to zero" lbd = self.lb.copy(); ubd = self.ub.copy() idxRmvd = self.getRwoEnz(genes) lbd[idxRmvd] = 0 ubd[idxRmvd] = 0 self.lb = lbd; self.ub = ubd; return self def removeReactions(self, lstRrm, updDict=True): " Have a hard removal of reactions from the model" " Temporarily not update dict['S']" if type(lstRrm) is str: lstRrm = [lstRrm] lbd = self.lb.copy(); ubd = self.ub.copy(); lstR = self.lstR.copy(); objC = self.objC.copy() Aeq = self.Aeq.copy(); O = self.dict['O']; S = self.dict['S'] for krm in lstRrm: lbd = lbd[lstR!=krm]; ubd = ubd[lstR!=krm] Aeqd = Aeq[:,lstR!=krm]; objCd = objC[lstR!=krm] lstRd = lstR[lstR!=krm]; if updDict: r1 = O.pop(krm) if updDict: self.checkM() self.lb = lbd; self.ub = ubd; self.Aeq = Aeqd; self.objC = objCd; self.lstR = lstRd; return self def addReactions(self, Ob, updDict=True): "Add or update reactions described in dict Ob" O = self.dict['O'].copy() O.update(Ob) self.checkM() def simGeneDeletion(self, lstGene=[''], lstEnz=None, objC=None, Aeq=None, beq=None, lb=None, ub=None, objOpt=None, fileOut=None, method='fba', solver=None, flxWT=None, option=None): """Simulate the FBA model with gene deletions lstGene - list of deletion genes (multiple deletion were connected with : lstEnz - lst of enzymes (annotated genes) f - function to optimize a linear combination of fluxes Aeq - Stoichiometric matrix for mass balance beq - the equality lb - lower bound ub - upper bound """ nR = len(self.dict['O']) if lstEnz is None: lstEnz = self.lstEnz if objC is None: objC = self.objC; lstR = self.lstR objOpt = objOpt or self.objOpt if Aeq is None: Aeq = self.Aeq; if beq is None: beq = self.beq if lb is None: lb = self.lb; if ub is None: ub = self.ub option = option or {} if fileOut!=None: fileOut = fileOut.rstrip('.csv')+'.csv' fo = open(fileOut, 'wb') csvwriter = csv.writer(fo) csvwriter.writerow(['deletant', 'f' + objOpt]) solWT = None i = 0 lstObj = [] lstSol = [] for KO in lstGene: i = i + 1 print '====================================' print 'Analyzing mutant ', i, '... ' idxRmvd = self.getRwoEnz(KO) for i1 in idxRmvd: print self.dict['O'][lstR[i1]] lbd = lb.copy(); ubd = ub.copy(); lbd[idxRmvd] = 0 ubd[idxRmvd] = 0 if method == 'fba': #solgd = self.optCbModel(solver=solver, lb=lbd, ub=ubd) solgd = self.optCbModel(lb=lbd, ub=ubd) elif method == 'moma': cbmd = self.copy() cbmd.lb = lbd; cbmd.ub = ubd; [solgd, solWT] = self.optTwoCbModels(cbmd, flxWT=flxWT, solver=solver, getSolOrigin=True, option=option) flxWT = solWT.x lstObj.append(solgd.obj) lstSol.append(solgd) if fileOut!=None: csvwriter.writerow([KO, str(solgd.obj)]) if fileOut!=None: fo.close() return lstObj, lstSol, solWT def addConstraints(self, constrR=None, constrM=None): """ Update the constraints specified in dict cstr Key to constrain over reaction R is either an integer corresponding to the reaction index in lstR or the reaction ID in string. Note: the constraints will be updated only if the new one is stricter than the old one. E.g. To add contraints over reaction 'R_1' with lower bound 0.05, and reaction number 10 (in lstR) with new lower and upper bound of 0 and 10. constrR = {'R_1': {'lb': 0.05}, 10:{'lb':0, 'ub': 10}} """ nR = len(self.lstR) lb1 = self.lb.copy(); ub1 = self.ub.copy() if constrR is not None: for i1 in constrR: c1 = constrR[i1] print i1 if type(i1) is not int: i1 = arange(nR).compress(self.lstR == i1)[0] if 'lb' in c1: lb1[i1] = max(lb1[i1], c1['lb']) if 'ub' in c1: ub1[i1] = min(ub1[i1], c1['ub']) self.lb = lb1; self.ub = ub1; return self def checkSolution_uptake(self, x, uptakeBound=None): """Check the uptake rate for nutrients given the solution flux vector x""" if uptakeBound is None: uptakeBound = self.lstMBound if type(uptakeBound) is dict: lstB = uptakeBound.keys() else: lstB = uptakeBound rB = [self.links(b1).keys() for b1 in lstB] nB = len(rB) rB1 = flxB = lbB = ubB = [] for n1 in range(nB): rs1 = rB[n1] for r1 in rs1: idx= self.lstR==r1 #print lstB[n1], ':', r1, ':', x[idx][0], (self.lb[idx][0], self.ub[idx][0]) print"%s : %s : %.5f ( %.5f, %.5f)" %(lstB[n1], r1, x[idx][0], self.lb[idx][0], self.ub[idx][0]) #print"%s [%s]: %s : %.5f ( %.5f, %.5f)" %(lstB[n1], self.dict['S'][lstB[n1]]['Name'], r1, x[idx][0], self.lb[idx][0], self.ub[idx][0]) rB1.append(r1); flxB.append(x[idx][0]); lbB.append(self.lb[idx][0]); ubB.append(self.ub[idx][0]) flxB = array(flxB); lbB = array(lbB); ubB = array(ubB); return rB, flxB, lbB, ubB def fva(self, objC=None, flxRef=None, lstv=None, fileOut=None, fixedObjC=False, constrR=None): """Simulate the FBA model with gene deletions f - function to optimize a linear combination of fluxes flxRef - reference flux solutions lstv - list of fluxes to be tested, default is to test all possible fluxes fileOut - name for file to store the results """ if objC is None: objC = self.objC lstR = self.lstR lb = self.lb.copy(); ub = self.ub.copy(); Aeq = self.Aeq; beq = self.beq if fileOut!=None: fileOut = fileOut.rstrip('.csv')+'.csv' fo = open(fileOut, 'wb') csvwriter = csv.writer(fo) csvwriter.writerow(['Reaction', 'lb', 'ub', 'flx_ref', 'flx_min', 'flx_max', 'flx_min_obj', 'flx_max_obj']) nR = Aeq.shape[1] i = 0 lst_sol_min = [] lst_sol_max = [] cbmod = self.copy() if constrR is not None: cbmod.addConstraints(constrR = constrR) if flxRef == None: # get reference (wildtype) flux distribution first! solRef = cbmod.optCbModel() flxRef = solRef.x fopt = dot(objC, flxRef) idxopt = arange(nR)[objC!=0] if lstv == None or len(lstv) == 0: lstv = array(list(set(range(nR)) - set(idxopt))) if lstR==None: lstR = vstack(([idxopt], lstv)) if fileOut!=None: csvwriter.writerow([lstR[idxopt], 0, 'inf', fopt]) lstvMin = []; lstvMax = []; lstObjMin = []; lstObjMax = []; lbv = lb.copy(); ubv = ub.copy(); if fixedObjC: lbv[idxopt] = fopt; ubv[idxopt] = fopt; cbmod.lb = lbv; cbmod.ub = ubv; print lstv for i in lstv: print 'Processing flux' + ':' + str(i) # Get minimum flux for v fv = zeros(nR); fv[i] = 1 solv = cbmod.optCbModel(objC=fv, objOpt='min') lst_sol_min.append(solv) fmin = solv.obj lstvMin.append(fmin) fobjMin = dot(objC, solv.x) lstObjMin.append(fobjMin) # Get maximum flux for v solv = cbmod.optCbModel(objC=fv, objOpt='max') lst_sol_max.append(solv) fmax = solv.obj lstvMax.append(fmax) fobjMax = dot(objC, solv.x) lstObjMax.append(fobjMax) if fileOut!=None: #csvwriter.writerow([lstR[i], lb[i], ub[i], flxRef[i], str(fmin), str(fmax)]) csvwriter.writerow([lstR[i], lb[i], ub[i], flxRef[i], fmin, fmax, fobjMin, fobjMax]) if fileOut!=None: fo.close() return lstvMin, lstvMax, lstObjMin, lstObjMax, lst_sol_min, lst_sol_max, lstv # # ------------------------------------------------------------------------------------------------- # utility used in xml reading #______________________________________________________________________________ # Queues: Stack, FIFOQueue, PriorityQueue class Queue: """Queue is an abstract class/interface. There are three types: Stack(): A Last In First Out Queue. FIFOQueue(): A First In First Out Queue. PriorityQueue(lt): Queue where items are sorted by lt, (default <). Each type supports the following methods and functions: q.append(item) -- add an item to the queue q.extend(items) -- equivalent to: for item in items: q.append(item) q.pop() -- return the top item from the queue len(q) -- number of items in q (also q.__len()) Note that isinstance(Stack(), Queue) is false, because we implement stacks as lists. If Python ever gets interfaces, Queue will be an interface.""" def __init__(self): abstract def extend(self, items): for item in items: self.append(item) def Stack(): """Return an empty list, suitable as a Last-In-First-Out Queue.""" return [] class FIFOQueue(Queue): """A First-In-First-Out Queue.""" def __init__(self): self.A = []; self.start = 0 def append(self, item): self.A.append(item) def __len__(self): return len(self.A) - self.start def extend(self, items): self.A.extend(items) def pop(self): e = self.A[self.start] self.start += 1 if self.start > 5 and self.start > len(self.A)/2: self.A = self.A[self.start:] self.start = 0 return e # #----------------------------------------------------------------------------------------------- # Modules for reading SBML file into CbModel # def readSBMLCompartments(modsbml): comps = {} ss = modsbml.getListOfCompartments() for s1 in ss: sid = s1.getId() nm1 = s1.getName() comps[sid] = nm1 return comps def readXMLNode(xmlnode = None): rdfLinks = '' if xmlnode is not None: fringe=FIFOQueue() fringe.append(xmlnode) while fringe: node1 = fringe.pop() if node1.isEnd(): rdfLinks = '|'.join([rdfLinks, node1.getAttributes().getValue(0)]) elif node1.isText(): #print node1.toString() rdfLinks = '|'.join([rdfLinks, node1.toString()]) else: nChild = node1.getNumChildren() children = [] for i in range(0, nChild): if node1.getChild(i) is not None: children.append(node1.getChild(i)) if len(children)>0 : fringe.extend(children) rdfLinks = rdfLinks.strip('|') return rdfLinks def readSBMLSpecies(modsbml): specs = {} ss = modsbml.getListOfSpecies() for s1 in ss: sid = s1.getId() nm1 = s1.getName() #sbo1 = s1.getSBOTerm() cmpt1 = s1.getCompartment() bnd1 = s1.getBoundaryCondition() # Get Notes note1 = readXMLNode(s1.getNotes()) # Get Annotation ann1 = readXMLNode(s1.getAnnotation()) specs[sid] = {'Name': nm1, 'Comp':cmpt1, 'Annot':ann1, 'Note':note1, 'Bound': bnd1} return specs def getSbmlModel(modsbml, useKeggid = False, useCompart=True, delComplex=True): specs = readSBMLSpecies(modsbml) comps = readSBMLCompartments(modsbml) lstRid = []; lstRnm = []; lstIn = []; lstOut = []; lstEC = []; lstLB = []; lstUB=[]; lstObjC = []; O={} rs = modsbml.getListOfReactions() for r1 in rs: rId = r1.getId(); rName = r1.getName(); rRev = r1.getReversible() isComplex = False outD = {} outSet = []; outs = r1.getListOfProducts() for s in outs: sid = s.getSpecies() if sid[0:2]=='Cx': isComplex =True if delComplex: break if useKeggid and 'kegg' in specs[sid]: skegg = specs[sid]['kegg'] if skegg == '': sid = modsbml.getSpecies(sid).getName() else: sid = skegg else: if useKeggid and sid[0:2]!='Cx': sid = modsbml.getSpecies(sid).getName() if useCompart: compart = modsbml.getSpecies(s.getSpecies()).getCompartment() compart = modsbml.getCompartment(compart).getName() sid = (compart, sid) outD[sid] = s.getStoichiometry() outSet.append(sid) if isComplex: if delComplex: continue rName = rName + ', Complex' inD = {}; inSet = [] ins = r1.getListOfReactants() for s in ins: sid = s.getSpecies() if useKeggid and 'kegg' in specs[sid]: skegg = specs[sid]['kegg'] if skegg == '': sid = modsbml.getSpecies(sid).getName() else: sid = skegg else: if useKeggid: sid = modsbml.getSpecies(sid).getName() if useCompart: compart = modsbml.getSpecies(s.getSpecies()).getCompartment() compart = modsbml.getCompartment(compart).getName() sid = (compart, sid) inD[sid] = s.getStoichiometry() inSet.append(sid) nmdfr = r1.getNumModifiers() rEnz = [] if nmdfr>0: for nm in range(0,nmdfr): rEnz.append(r1.getModifier(nm).getSpecies()) rEnz = '|'.join(rEnz) else: rEnz = '' notes = r1.getNotesString() if 'GENE_ASSOCIATION: ' in notes: rEnz = notes.split('GENE_ASSOCIATION: ')[1].split('')[0].replace('(','').replace(')','').replace(' and ', ':'). replace(' or ', '|').replace(' ','').strip() # Get annotations (references to EC code) ann1 = r1.getAnnotation() dbLinks = readXMLNode(ann1) EC = ''; # get EC or PubMed if dbLinks != '': lstDB = dbLinks.split('|') for s in lstDB: if s.find('ec-code')>=0: if s.find('#')>=0: EC = s.split('#')[1] # for consensus model (Aug2008) elif s.find('ec-code:')>=0: EC = s.split('ec-code:')[1] # for iIN800 model (Nov 2008) if EC == '': notes = r1.getNotesString() if 'PROTEIN_CLASS: ' in notes: EC = notes.split('PROTEIN_CLASS: ')[1].split('')[0].replace('(','').replace(')','').replace(' and ', ':').replace(' or ', '|').replace(' ','').strip() objc1 = ''; lb1 = ''; ub1 = ''; flx1 = ''; if r1.getKineticLaw() != None: pars = r1.getKineticLaw().getListOfParameters() for par in pars: if 'LOWER_BOUND' in [par.getId(), par.getName()]: lb1 = par.getValue() if 'UPPER_BOUND' in [par.getId(), par.getName()]: ub1 = par.getValue() if 'OBJECTIVE_COEFFICIENT' in [par.getId(), par.getName()]: objc1 = par.getValue() if 'FLUX_VALUE' in [par.getId(), par.getName()]: flx1 = par.getValue() # Get annotations (references to EC code) note1 = readXMLNode(r1.getNotes()) notes = note1.split('|'); ssys1 = '' for note1 in notes: if 'SUBSYSTEM:' in note1: ssys1 = note1.split('SUBSYSTEM: ')[1] lstIn.append(set(inSet)) lstOut.append(set(outSet)) lstRid.append(r1.getId()) lstRnm.append(r1.getName()) lstEC.append(EC) lstLB.append(lb1) lstUB.append(ub1) lstObjC.append(objc1) O[rId] = {'Name': rName, 'In':inD, 'Out':outD, 'Rev':rRev, 'Enz':rEnz, 'EC':EC, 'LB':lb1, 'UB':ub1, 'Objc':objc1, 'Flux':flx1, 'Subs':ssys1} # Get EC no if len(specs)>0: nEq = 0; nAdd = 0 for rId in O: r1 = O[rId] if r1['Enz'] != '' and r1['Enz'] in specs and 'EC' in specs[r1['Enz']]: ec1 = specs[r1['Enz']]['EC'] else: continue if ec1 == r1['EC']: if ec1!='': nEq = nEq+1 else: if r1['EC'] =='': r1['EC'] = ec1 nAdd = nAdd + 1 # borrow EC from Enz property for rId in O: r1 = O[rId] if r1['Enz']!='' and r1['Enz'] in specs: e1 = specs[r1['Enz']]['Name'] O[rId]['Enz']=e1 return O, specs, comps def get_reaction_enz(modsbml, useKeggid = False, useCompart=True, delComplex=True): specs = readSBMLSpecies(modsbml) comps = readSBMLCompartments(modsbml) lstRid = []; lstRnm = []; lstIn = []; lstOut = []; lstEC = []; lstLB = []; lstUB=[]; lstObjC = []; O={} rs = modsbml.getListOfReactions() for r1 in rs: rId = r1.getId(); rName = r1.getName(); rRev = r1.getReversible() isComplex = False outD = {} outSet = []; outs = r1.getListOfProducts() for s in outs: sid = s.getSpecies() if sid[0:2]=='Cx': isComplex =True if delComplex: break if useKeggid and 'kegg' in specs[sid]: skegg = specs[sid]['kegg'] if skegg == '': sid = modsbml.getSpecies(sid).getName() else: sid = skegg else: if useKeggid and sid[0:2]!='Cx': sid = modsbml.getSpecies(sid).getName() if useCompart: compart = modsbml.getSpecies(s.getSpecies()).getCompartment() compart = modsbml.getCompartment(compart).getName() sid = (compart, sid) outD[sid] = s.getStoichiometry() outSet.append(sid) if isComplex: if delComplex: continue rName = rName + ', Complex' inD = {}; inSet = [] ins = r1.getListOfReactants() for s in ins: sid = s.getSpecies() if useKeggid and 'kegg' in specs[sid]: skegg = specs[sid]['kegg'] if skegg == '': sid = modsbml.getSpecies(sid).getName() else: sid = skegg else: if useKeggid: sid = modsbml.getSpecies(sid).getName() if useCompart: compart = modsbml.getSpecies(s.getSpecies()).getCompartment() compart = modsbml.getCompartment(compart).getName() sid = (compart, sid) inD[sid] = s.getStoichiometry() inSet.append(sid) nmdfr = r1.getNumModifiers() rEnz = [] if nmdfr>0: for nm in range(0,nmdfr): rEnz.append(r1.getModifier(nm).getSpecies()) rEnz = '|'.join(rEnz) else: rEnz = '' notes = r1.getNotesString() if 'GENE_ASSOCIATION: ' in notes: rEnz = notes.split('GENE_ASSOCIATION: ')[1].split('')[0].replace('(','').replace(')','').replace(' and ', ':'). replace(' or ', '|').replace(' ','').strip() # Get annotations (references to EC code) ann1 = r1.getAnnotation() dbLinks = readXMLNode(ann1) EC = ''; # get EC or PubMed if dbLinks != '': lstDB = dbLinks.split('|') for s in lstDB: if s.find('ec-code')>=0: if s.find('#')>=0: EC = s.split('#')[1] # for consensus model (Aug2008) elif s.find('ec-code:')>=0: EC = s.split('ec-code:')[1] # for iIN800 model (Nov 2008) if EC == '': notes = r1.getNotesString() if 'PROTEIN_CLASS: ' in notes: EC = notes.split('PROTEIN_CLASS: ')[1].split('')[0].replace('(','').replace(')','').replace(' and ', ':').replace(' or ', '|').replace(' ','').strip() objc1 = ''; lb1 = ''; ub1 = ''; flx1 = ''; if r1.getKineticLaw() != None: pars = r1.getKineticLaw().getListOfParameters() for par in pars: if 'LOWER_BOUND' in [par.getId(), par.getName()]: lb1 = par.getValue() if 'UPPER_BOUND' in [par.getId(), par.getName()]: ub1 = par.getValue() if 'OBJECTIVE_COEFFICIENT' in [par.getId(), par.getName()]: objc1 = par.getValue() if 'FLUX_VALUE' in [par.getId(), par.getName()]: flx1 = par.getValue() # Get annotations (references to EC code) note1 = readXMLNode(r1.getNotes()) notes = note1.split('|'); ssys1 = '' for note1 in notes: if 'SUBSYSTEM:' in note1: ssys1 = note1.split('SUBSYSTEM: ')[1] lstIn.append(set(inSet)) lstOut.append(set(outSet)) lstRid.append(r1.getId()) lstRnm.append(r1.getName()) lstEC.append(EC) lstLB.append(lb1) lstUB.append(ub1) lstObjC.append(objc1) O[rId] = {'Name': rName, 'In':inD, 'Out':outD, 'Rev':rRev, 'Enz':rEnz, 'EC':EC, 'LB':lb1, 'UB':ub1, 'Objc':objc1, 'Flux':flx1, 'Subs':ssys1} lstR = [] lstRt = [] lstRt1 = [] for i in range(0,len(lstOut)): strin = str(lstIn[i]).replace('set', 'reactants').replace('[','').replace(']','') strout = str(lstOut[i]).replace('set', 'products').replace('[','').replace(']','') lstR.append(strin+'<->'+strout) lstRt.append((lstRid[i], strin+'<->'+strout)) lstRt1.append((lstRid[i], lstRnm[i], strin+'<->'+strout, lstEC[i])) dictRC = dict(lstRt) import operator # Sort according to the reactant and products lstRt_sorted = sorted(lstRt, key=operator.itemgetter(1)) lstR_sorted=map(operator.itemgetter(1), lstRt_sorted) lstRt1_sorted = sorted(lstRt1, key=operator.itemgetter(1)) lstRt_sorted=lstRt1_sorted # Get EC no if len(specs)>0: nEq = 0; nAdd = 0 for rId in O: r1 = O[rId] if r1['Enz'] != '' and r1['Enz'] in specs and 'EC' in specs[r1['Enz']]: ec1 = specs[r1['Enz']]['EC'] else: continue if ec1 == r1['EC']: if ec1!='': nEq = nEq+1 else: if r1['EC'] =='': r1['EC'] = ec1 nAdd = nAdd + 1 for rId in O: r1 = O[rId] if r1['Enz']!='' and r1['Enz'] in specs: #e1 = modsbml.getSpecies(r1['Enz']).getName() e1 = specs[r1['Enz']]['Name'] print e1 O[rId]['Enz']=e1 return O, specs, comps, lstRt_sorted, lstRt1_sorted # # --------------------------------------------------------------------------------------------- # Usage example: def example1_ind750fba(sbmlFile=None): """ Usage example1 Step 1. importing iND750 model in xml file Step 2. simulating wild-type and a deletion mutant's growth using standard FBA """ print "=== Example1: import and simulate a metabolic model ===" if sbmlFile == None: sbmlFile = "Sc_iND750_GlcMM.xml" if not os.path.isfile(sbmlFile): print "SBML file for model not found: %s" %(sbmlFile) return print "Step 1. Reading SBML file " cbm0 = CbModel() cbm0.readSbmlCbModel(sbmlFile=sbmlFile) # Define minimal media by setting maximum uptake rate for exchange metabolites (nutrients) # Metabolite IDs are specific for iND750 model uptakeBound_mm = {'M_o2_b': 6.3, 'M_nh4_b': 100, # Ammonium (nitrogen source) 'M_so4_b': 100, # sulfate 'M_pi_b': 0.89, # phospate 'M_glc_D_b': 22.6, # D-glucose 'M_his_L_b': 0.082, # histidin 'M_leu_L_b': 0.4, # Leu 'M_met_L_b': 0.0468, # Met 'M_k_b': 4.44, # Potassium "M_na1_b": 0.75, # Sodium "M_btn_b": 0.00000142, #Biotin "M_chol_b": 0.000092, "M_inost_b": 0.00193, "M_pnto_R_b": 0.0002, "M_ribflv_b": 0.00063, "M_ura_b": 0.4 } # Set uptake bound for the model cbm0.setUptakeBound(uptakeBound_mm) print "Step 2.1. Running standard fba simulation for wild-type " sol_wt = cbm0.optCbModel() print "Step 2.2. Checking uptake rate for exchange metabolites " rex = cbm0.checkSolution_uptake(sol_wt.x, uptakeBound_mm) print "Step 2.3 Simulating deletion mutant's growth by standard fba, " print " MOMA could be used by setting method = 'moma' instead of default 'fba' " [lstGr1, lstSol1, null] = cbm0.simGeneDeletion(lstGene=['YGR061C'], method = 'fba', fileOut='test1.csv') def main(): args = sys.argv[1:] if not args: sys.stderr.write("Usage:"+ __doc__+"\n") return if args[0][0:2] == '-e': if len(args)>1: sbmlFile = args[1] else: sbmlFile = None example1_ind750fba(sbmlFile) else: sys.stderr.write("Usage:"+ __doc__+"\n") return if __name__ == '__main__': main()