Source code for cbmpy.CBWrite

"""
CBMPy: CBWrite module
=====================
PySCeS Constraint Based Modelling (http://cbmpy.sourceforge.net)
Copyright (C) 2009-2015 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBWrite.py 427 2016-04-05 09:28:55Z bgoli $)

"""

# preparing for Python 3 port
from __future__ import division, print_function
from __future__ import absolute_import
#from __future__ import unicode_literals

if 'cDir' in vars():
    cDir=vars()['cDir']
else:
    cDir=None
import os, time, numpy, zipfile
# this is a hack that needs to be streamlined a bit
try:
    import cStringIO as csio
except ImportError:
    import io as csio

from . import CBTools, CBXML, CBDataStruct


from .CBConfig import __CBCONFIG__ as __CBCONFIG__
__DEBUG__ = __CBCONFIG__['DEBUG']
__version__ = __CBCONFIG__['VERSION']

HAVE_SYMPY = False
try:
    import sympy
    if int(sympy.__version__.split('.')[1]) >= 7 and int(sympy.__version__.split('.')[2]) >= 5:
        HAVE_SYMPY = True
    elif int(sympy.__version__.split('.')[0]) >= 1:
        HAVE_SYMPY = True
    else:
        del sympy
        print('\nWARNING: SymPy version 0.7.5 or newer is required for symbolic matrix support.')
except ImportError:
    HAVE_SYMPY = False
    print('\nERROR: SymPy import error (required for symbolic matrix support only).')

_HAVE_XLWT_ = False
try:
    import xlwt
    _HAVE_XLWT_ = True
except ImportError:
    print('\nINFO: No xlwt module available, Excel spreadsheet creation disabled')

[docs]def writeSBML3FBC(fba, fname, directory=None, gpr_from_annot=False,\ add_groups=False, add_cbmpy_annot=True, add_cobra_annot=False,\ xoptions={'fbc_version': 1, 'validate' : False, 'compress_bounds' : True}): """ Takes an FBA model object and writes it to file as SBML L3 FBC: - *fba* an fba model object - *fname* the model will be written as XML to *fname* - *directory* [default=None] if defined it is prepended to fname - *gpr_from_annot* [default=True] if enabled will attempt to add the gene protein associations from the annotations if no gene protein association objects exist - *add_cbmpy_annot* [default=True] add CBMPy KeyValueData annotation. Replaces <notes> - *add_cobra_annot* [default=True] add COBRA <notes> annotation - *xoptions* extended options - *validate* [default=False] validate the output SBML file - *compress_bounds* [default=False] try compress output flux bound parameters """ sbml_level_version = (3,1) autofix=True return_fbc=False return CBXML.sbml_writeSBML3FBC(fba, fname, directory, sbml_level_version, autofix, return_fbc,\ gpr_from_annot, add_groups, add_cbmpy_annot, add_cobra_annot, xoptions)
[docs]def writeSBML3FBCV2(fba, fname, directory=None, gpr_from_annot=False, add_groups=False, add_cbmpy_annot=True, add_cobra_annot=False,\ validate=False, compress_bounds=True): """ Takes an FBA model object and writes it to file as SBML L3 FBCv2 : - *fba* an fba model object - *fname* the model will be written as XML to *fname* - *directory* [default=None] if defined it is prepended to fname - *gpr_from_annot* [default=False] if enabled will attempt to add the gene protein associations from the annotations - *add_groups* [default=False] add SBML3 groups (if supported by libSBML) - *add_cbmpy_annot* [default=True] add CBMPy KeyValueData annotation. Replaces <notes> - *add_cobra_annot* [default=False] add COBRA <notes> annotation - *validate* [default=False] validate the output SBML file - *compress_bounds* [default=True] try compress output flux bound parameters """ xoptions = {'fbc_version': 2, 'validate' : validate, 'compress_bounds' : compress_bounds} sbml_level_version=(3,1) autofix=True, return_fbc=False #if fbc_version == 2: #add_cobra_annot = False return CBXML.sbml_writeSBML3FBC(fba, fname, directory, sbml_level_version, autofix, return_fbc,\ gpr_from_annot, add_groups, add_cbmpy_annot, add_cobra_annot, xoptions)
[docs]def writeCOBRASBML(fba, fname, directory=None): """ Takes an FBA model object and writes it to file as a COBRA compatible : - *fba* an fba model object - *fname* the model will be written as XML to *fname* - *directory* [default=None] if defined it is prepended to fname """ return CBXML.sbml_writeCOBRASBML(fba, fname, directory)
[docs]def writeSBML2FBA(fba, fname, directory=None, sbml_level_version=None): """ Takes an FBA model object and writes it to file as SBML L2 with FBA annotations. Note if you want to write BiGG/FAME style annotations then you must use *sbml_level_version=(2,1)* - *fba* an fba model object - *fname* the model will be written as XML to *fname* - *sbml_level_version* [default=None] a tuple containing the SBML level and version e.g. (2,1) This is a utility wrapper for the function `CBXML.sbml_writeSBML2FBA` """ CBXML.sbml_writeSBML2FBA(fba, fname, directory, sbml_level_version)
[docs]def writeSensitivitiesToCSV(sensitivities, fname): """ Write out a sensitivity report using the objective sensitivities and bound sensitivity dictionaries created by e.g. cplx_getSensitivities(). - *sensitivity* tuple containing - *obj_sens* dictionary of objective coefficient sensitivities (per flux) - *rhs_sens* dictionary of constraint rhs sensitivities (per constraint) - *bound_sens* dictionary of bound sensitivities (per flux) - *fname* output filename e.g. fname.csv """ obj_sens = sensitivities[0] rhs_sens = sensitivities[1] bound_sens = sensitivities[2] F = file(fname+'_flux_sensitivity.csv', 'w') head = "Flux,Reduced cost,OCS low,OC value,OCS high,LB low,LB high,UB low,UB high" F.write(head+'\n') for j in obj_sens: rc = obj_sens[j][0] lcs = obj_sens[j][1] ocv = obj_sens[j][2] ucs = obj_sens[j][3] if j in bound_sens: lbs = bound_sens[j][0] lb = bound_sens[j][1] ub = bound_sens[j][2] ubs = bound_sens[j][3] else: lbs = 0 lb = 0 ub = 0 ubs = 0 F.write('%s,%s,%s,%s,%s,%s,%s,%s,%s\n' % (j,rc,lcs,ocv,ucs,lbs,lb,ub,ubs)) for j in bound_sens: if j not in obj_sens: rc = 0 lcs = 0 ocv = 0 ucs = 0 lbs = bound_sens[j][0] lb = bound_sens[j][1] ub = bound_sens[j][2] ubs = bound_sens[j][3] F.write('%s,%s,%s,%s,%s,%s,%s,%s,%s\n' % (j,rc,lcs,ocv,ucs,lbs,lb,ub,ubs)) F.flush() F.close() F = file(fname+'_constraint_sensitivity.csv', 'w') F.write('Constraint,RHS low,RHS,RHS high\n') for c in rhs_sens: F.write('%s,%s,%s,%s\n' % (c, rhs_sens[c][0], rhs_sens[c][1], rhs_sens[c][2])) F.flush() F.close()
[docs]def WriteModelRaw(fba, work_dir=cDir): """ INFO: this method will be deprecated please update your scripts to use \"writeModelRaw()\" """ print('\nINFO: this method will be deprecated please update your scripts to use \"writeModelRaw()\"\n') time.sleep(1) writeModelRaw(fba, work_dir)
[docs]def writeModelRaw(fba, work_dir=cDir): """ Writes a fba (actually just dumps it) to a text file. - *fba* an instantiated FBAmodel instance - *work_dir* directory designated for output """ if work_dir == None: work_dir = os.getcwd() FF = file(os.path.join(work_dir,'WriteModelRawOutput.txt'), 'w') FF.write('Species information:\n\n') for s in fba.species: FF.write('%s: value=%f, is_boundary=%s, name=%s\n' % (s.getPid(), s.value, int(s.is_boundary), s.name)) FF.write('\nReaction information:\n\n') for r in fba.reactions: FF.write('%s: reversible=%s, name=%s\n' % (r.getPid(), int(r.reversible), r.name)) FF.write('\t%s\n' % r.getSpeciesIds()) FF.write('\nConstraint information:\n\n') for c in fba.flux_bounds: FF.write('%s %s %f\n' % (c.reaction, c.operation, c.value)) FF.write('\nObjective information:\n\n') FF.write('Active Objective: %s (%s)\n' % (fba.objectives[fba.activeObjIdx].getPid(), fba.objectives[fba.activeObjIdx].operation)) for o in fba.objectives: FF.write('%s: %s, %s\n' % (o.getPid(), o.operation, o.getFluxObjectiveReactions())) if hasattr(fba.N, 'shape'): FF.write('\nStoichiometric information:\n\n') FF.write('N-matrix dimensions = (%s,%s)\n' % fba.N.shape) FF.close() print('WriteModelRaw has written a file to {}'.format(os.path.join(work_dir,'WriteModelRawOutput.txt')))
[docs]def BuildLPFluxBounds(fba, use_rational=False): """ Build and return a csio that contains the flux bounds in LP format """ if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational IO') FFS = csio.StringIO() c2s = {} for c in fba.flux_bounds: if __DEBUG__: print('%s: %s' % (c.getPid(), c.operation)) minv = None maxv = None oper = None R = c.reaction if c.operation in ['less','lessEqual']: maxv = c.value oper = c.operation if oper == 'less': oper = '<=' if __DEBUG__: print('LP Bounds operator must be <= not <') else: oper = '<=' if c.operation in ['greater','greaterEqual']: minv = c.value oper = c.operation if oper == 'greater': oper = '<=' if __DEBUG__: print('LP Bounds operator must be <= not <') else: oper = '<=' if c.operation == 'equals': maxv = c.value oper = '=' if maxv != None: if use_rational: c2s.update({R+'r' : '%s %s' % (oper, sympy.Rational(maxv))}) else: if maxv == float('inf'): c2s.update({R+'r' : '%s +%s' % (oper, maxv)}) else: c2s.update({R+'r' : '%s %s' % (oper, maxv)}) elif minv != None: if use_rational: c2s.update({R+'l' : '%s %s' % (sympy.Rational(minv), oper)}) else: c2s.update({R+'l' : '%s %s' % (minv, oper)}) c2sk = list(c2s) c2sk.sort() for r in (c2sk): rht = None lht = None R = r[:-1] if R+'r' in c2s: rht = c2s[R+'r'] c2s.pop(R+'r') if R+'l' in c2s: lht = c2s[R+'l'] c2s.pop(R+'l') if lht != None and rht != None: if __DEBUG__: print('{} {} {}'.format(lht, R, rht)) ## FFS.write('%s: %s %s %s\n' % (R, lht, R, rht)) FFS.write('%s %s %s\n' % (lht, R, rht)) # THIS MUST BE SO for GLPK elif lht == None and rht == None: if __DEBUG__: print('Skipping: ({}, {}, {})'.format(lht, R, rht)) elif lht != None: if __DEBUG__: print('{} {}'.format(lht, R)) #FFS.write('%s: %s %s\n' % (R, lht, R)) FFS.write('%s %s\n' % (lht, R)) elif rht != None: if __DEBUG__: print('{} {}'.format(R, rht)) #FFS.write('%s: %s %s\n' % (R, R, rht)) FFS.write('%s %s\n' % (R, rht)) else: print('CONFUSION: ({}, {}, {})'.format(lht, R, rht)) return FFS
[docs]def BuildLPConstraints(fba, use_rational=False): """ Build and return a csio that contains constraint constructed from the StoichiometeryLP object - *fba* an fba model object which has a stoichiometry - *use_rational* write rational number output [default=False] """ if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational output') rebuild_stoich = False if not hasattr(fba,'N') or fba.N == None: rebuild_stoich = True else: if len([s for s in fba.species if not s.is_boundary]) != fba.N.array.shape[0]: rebuild_stoich = True elif len(fba.reactions) != fba.N.array.shape[1]: rebuild_stoich = True if rebuild_stoich: print('\nWarning FBA object has inconsistant stoichiometric matrix, rebuilding it now.') CBTools.addStoichToFBAModel(fba) constr = {} for r in range(fba.N.array.shape[0]): rowName = fba.N.row[r] RCon = [] for c in range(fba.N.array.shape[1]): colName = fba.N.col[c] colCoef = fba.N.array[r,c] if colCoef != 0.0: RCon.append((colCoef, colName)) constr.update({rowName : RCon}) FFS = csio.StringIO() constrsk = list(constr) constrsk.sort() ## for r in constrsk: for r in range(len(fba.N.row)): if len(constr[fba.N.row[r]]) > 0: FFS.write(' %s: ' % fba.N.row[r]) ## FFS.write(' ') for col in constr[fba.N.row[r]]: if use_rational: if col[0] > 0.0: FFS.write('+%s %s ' % (sympy.Rational('%s' % col[0]), col[1])) else: FFS.write('%s %s ' % (sympy.Rational('%s' % col[0]), col[1])) else: if col[0] > 0.0: if col[0] == 1.0: FFS.write('+ %s ' % (col[1])) else: FFS.write('+%.20f %s ' % (col[0], col[1])) else: if col[0] == -1.0: FFS.write('- %s ' % (col[1])) else: FFS.write('%.20f %s ' % (col[0], col[1])) operator = fba.N.operators[r].upper() if operator == 'E' or operator == '=': operator = '=' elif operator == 'G' or operator == '>' or operator == '>=': operator = '>=' elif operator == 'L' or operator == '<' or operator == '<=': operator = '<=' else: raise RuntimeError('\nINVALID operator: %s' % operator) FFS.write('%s ' % operator) FFS.write('%s\n' % fba.N.RHS[r]) return FFS
[docs]def BuildLPUserConstraints(fba, use_rational=False): """ Build and return a csio that contains constraint constructed from the StoichiometeryLP object - *fba* an fba model object which has a stoichiometry - *use_rational* write rational number output [default=False] """ if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational output') rebuild_stoich = False assert fba.user_constraints != None, "\nNo user constraints to build" if not hasattr(fba,'CM') or fba.CM == None: rebuild_stoich = True if rebuild_stoich: print('\nWarning FBA object has inconsistant user constraint matrix, rebuilding it now.') CBTools.addStoichToFBAModel(fba) constr = {} for r in range(fba.CM.array.shape[0]): rowName = fba.CM.row[r] RCon = [] for c in range(fba.CM.array.shape[1]): colName = fba.CM.col[c] colCoef = fba.CM.array[r,c] if colCoef != 0.0: RCon.append((colCoef, colName)) constr.update({rowName : RCon}) FFS = csio.StringIO() constrsk = list(constr) constrsk.sort() ## for r in constrsk: for r in range(len(fba.CM.row)): if len(constr[fba.CM.row[r]]) > 0: FFS.write(' %s: ' % fba.CM.row[r]) ## FFS.write(' ') for col in constr[fba.CM.row[r]]: if use_rational: if col[0] > 0.0: FFS.write('+%s %s ' % (sympy.Rational('%s' % col[0]), col[1])) else: FFS.write('%s %s ' % (sympy.Rational('%s' % col[0]), col[1])) else: if col[0] > 0.0: if col[0] == 1.0: FFS.write('+ %s ' % (col[1])) else: FFS.write('+%.20f %s ' % (col[0], col[1])) else: if col[0] == -1.0: FFS.write('- %s ' % (col[1])) else: FFS.write('%.20f %s ' % (col[0], col[1])) operator = fba.CM.operators[r] if operator == 'E' or operator == '=': operator = '=' elif operator == 'G' or operator == '>' or operator == '>=': operator = '>=' elif operator == 'L' or operator == '<' or operator == '<=': operator = '<=' else: raise RuntimeError('\nINVALID operator: %s' % operator) FFS.write('%s ' % operator) FFS.write('%s\n' % fba.CM.RHS[r]) return FFS
[docs]def BuildLPConstraintsRelaxed(fba): """ Build and return a csio that contains the constaints in LP format Relaxed refers to dS/dt >= 0 """ raise DeprecationWarning("\nThis method is deprecated") if not hasattr(fba,'N') or fba.N == None: print('\nWarning FBA object has no stoichiometric matrix constructing it now.') CBTools.addStoichToFBAModel(fba) time.sleep(1) constr = {} for r in range(fba.N.array.shape[0]): rowName = fba.N.row[r] RCon = [] for c in range(fba.N.array.shape[1]): colName = fba.N.col[c] colCoef = fba.N.array[r,c] if colCoef != 0.0: RCon.append((colCoef, colName)) constr.update({rowName : RCon}) FFS = csio.StringIO() constrsk = list(constr) constrsk.sort() for r in fba.N.row: if len(constr[r]) > 0: FFS.write('%s: ' % r) ## FFS.write(' ') for col in constr[r]: if col[0] > 0.0: FFS.write('+%.20f %s ' % (col[0], col[1])) else: FFS.write('%.20f %s ' % (col[0], col[1])) FFS.write('>= 0\n') return FFS
[docs]def BuildLPConstraintsStrict(fba, use_rational=False): """ Build and return a csio that contains the constaints in LP format Strict refers to dS/dt = 0 """ ## print "Consider using the new BuildLPConstraints() method" if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational IO') if not hasattr(fba,'N') or fba.N == None: print('\nWarning FBA object has no stoichiometric matrix constructing it now.') CBTools.addStoichToFBAModel(fba) time.sleep(1) constr = {} for r in range(fba.N.array.shape[0]): rowName = fba.N.row[r] RCon = [] for c in range(fba.N.array.shape[1]): colName = fba.N.col[c] colCoef = fba.N.array[r,c] if colCoef != 0.0: RCon.append((colCoef, colName)) constr.update({rowName : RCon}) FFS = csio.StringIO() constrsk = list(constr) constrsk.sort() ## for r in constrsk: for r in fba.N.row: if len(constr[r]) > 0: FFS.write(' %s: ' % r) ## FFS.write(' ') for col in constr[r]: if use_rational: if col[0] > 0.0: FFS.write('+%s %s ' % (sympy.Rational('%s' % col[0]), col[1])) else: FFS.write('%s %s ' % (sympy.Rational('%s' % col[0]), col[1])) else: if col[0] > 0.0: if col[0] == 1.0: FFS.write('+ %s ' % (col[1])) else: FFS.write('+%.20f %s ' % (col[0], col[1])) else: if col[0] == -1.0: FFS.write('- %s ' % (col[1])) else: FFS.write('%.20f %s ' % (col[0], col[1])) if len(constr[r]) != 0: FFS.write('= 0\n') else: FFS.write('\n') return FFS
[docs]def BuildLPConstraintsMath(fba, use_rational=False): """ Build and return a csio that contains the constaints in LP format Strict refers to dS/dt => 0 and dS/dt <= 0 """ raise DeprecationWarning("\nThis method is deprecated") if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational IO') if not hasattr(fba,'N') or fba.N == None: print('\nWarning FBA object has no stoichiometric matrix constructing it now.') CBTools.addStoichToFBAModel(fba) time.sleep(1) constr = {} for r in range(fba.N.array.shape[0]): rowName = fba.N.row[r] RCon = [] for c in range(fba.N.array.shape[1]): colName = fba.N.col[c] colCoef = fba.N.array[r,c] if colCoef != 0.0: RCon.append((colCoef, colName)) constr.update({rowName : RCon}) FFS = csio.StringIO() constrsk = list(constr) constrsk.sort() for r in fba.N.row: if len(constr[r]) > 0: FFS.write('%sn1: ' % r) ## FFS.write(' ') for col in constr[r]: if col[0] > 0.0: FFS.write('+%s %s ' % (col[0], col[1])) else: FFS.write('%s %s ' % (col[0], col[1])) # This is a fudge FFS.write('>= 0\n') for r in fba.N.row: if len(constr[r]) > 0: FFS.write('%sn2: ' % r) for col in constr[r]: if -col[0] > 0.0: FFS.write('+%s %s ' % (-col[0], col[1])) else: FFS.write('%s %s ' % (-col[0], col[1])) # This is a fudge FFS.write('>= 0\n') return FFS
[docs]def WriteModelLPOld(fba, work_dir=None, multisymb=' ', lpt=True, constraint_mode='strict', use_rational=False, format='%s'): """ INFO: this method will be deprecated please update your scripts to use \"writeModelLPOld()\" """ print('\nINFO: this method will be deprecated please update your scripts to use \"writeModelLPOld()\"\n') time.sleep(1) writeModelLPOld(fba, work_dir, multisymb, lpt, constraint_mode, use_rational, format)
[docs]def writeModelLPOld(fba, work_dir=None, multisymb=' ', lpt=True, constraint_mode='strict', use_rational=False, format='%s'): """ Writes a fba as an LP/LPT - *fba* an instantiated FBAmodel instance - *work_dir* directory designated for output - *multisymb* the multiplication symbol (default: <space>) - *lpt* the file format (default: True for lpt) or False for lp """ print("\nTHIS FUNCTION IS DEPRECATED\n") if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational IO') FNAME = None if work_dir == None: FnameTmp = fba.getPid() else: FnameTmp = os.path.join(work_dir, fba.getPid()) if use_rational: FnameTmp = FnameTmp+'.rat' if not lpt: FNAME = FnameTmp+'.lp' FF = file(FNAME, 'w') FF.write('Problem\n %s\n\n' % FnameTmp) else: FNAME = FnameTmp+'.lp' FF = file(FNAME, 'w') FF.write('\\\\ %s \n\n' % FnameTmp) objO = fba.objectives[fba.activeObjIdx].operation.lower() objO = objO[0].upper() + objO[1:] FF.write('%s\n' % objO) objStr = '%s_objf: ' % fba.objectives[fba.activeObjIdx].getPid() for fObj in fba.objectives[fba.activeObjIdx].fluxObjectives: sign = None nc = 0.0 try: nc = float(fObj.coefficient) except ValueError: print('Suspected rational number ({}) detected in fluxObjective {}'.format(fObj.coefficient, fObj.getPid())) if nc >= 0.0: sign = '+' else: sign = '-' # TODO: if use_rational is not used simply try and evaluate the coefficient string with # sympy.Rational.evalf() and use this as the value for nc ... also remove use_rational case if use_rational: objStr += ' %s%s%s' % (sympy.Rational(fObj.coefficient), multisymb, fObj.reaction) ## FF.write('%s: %s%s%s\n' % (fobj0.reaction, sympy.Rational(fobj0.coefficient), multisymb, fobj0.reaction)) else: objStr += ' %s %s%s%s' % (sign, abs(nc), multisymb, fObj.reaction) ## FF.write('%s: %s%s%s\n' % (fobj0.reaction, fobj0.coefficient, multisymb, fobj0.reaction)) FF.write('%s\n' % objStr) if constraint_mode == 'math': CONST = BuildLPConstraintsMath(fba) elif constraint_mode == 'relaxed': CONST = BuildLPConstraintsRelaxed(fba) else: CONST = BuildLPConstraintsStrict(fba, use_rational) CONST.seek(0) BOUNDS = BuildLPFluxBounds(fba, use_rational) BOUNDS.seek(0) if __DEBUG__: print(CONST.read()); CONST.seek(0) if __DEBUG__: print(BOUNDS.read()); BOUNDS.seek(0) FF.write('\nSubject To\n') FF.write(CONST.read()) FF.write('\nBounds\n') FF.write(BOUNDS.read()) FF.write('\nEND\n') FF.close() print('writeModelLP has written a file to {}'.format(NAME)) return FNAME
[docs]def WriteModelLP(fba, work_dir=None, fname=None, multisymb=' ', format='%s', use_rational=False, constraint_mode=None, quiet=False): """ INFO: this method will be deprecated please update your scripts to use \"writeModelLP()\" """ print('\nINFO: this method will be deprecated please update your scripts to use \"writeModelLP()\"\n') time.sleep(1) writeModelLP(fba, work_dir, fname, multisymb, format, use_rational, constraint_mode, quiet)
[docs]def writeModelLP(fba, work_dir=None, fname=None, multisymb=' ', format='%s', use_rational=False, constraint_mode=None, quiet=False): """ Writes an FBA object as an LP in CPLEX LP format - *fba* an instantiated FBAmodel instance - *work_dir* directory designated for output - *fname* the file name [default=fba.getPid()] - *multisymb* the multiplication symbol (default: <space>) - *format* the number format of the output - *use_rational* output rational numbers [default=False] - *quiet* [default=False] supress information messages """ if constraint_mode != None: print("\nConstraint_mode has been deprecated") time.sleep(5) if not _HAVE_SYMPY_ and use_rational: use_rational = False print('\nWarning switching to floating point arithmetic: install Sympy for rational IO') time.sleep(2) FNAME = None if fname == None: fname = fba.getPid() if work_dir != None: FnameTmp = os.path.join(work_dir, fname) else: FnameTmp = fname FNAME = FnameTmp+'.lp' FF = file(FNAME, 'w') FF.write('\\\\ %s \n\n' % FnameTmp) if len(fba.objectives) > 0: ## print fba.objectives if fba.objectives[fba.activeObjIdx].operation == None: print('\nWARNING: Objective function \"{}\" has no \"operation\" defined assuming \"maximize\"'.format(fba.objectives[fba.activeObjIdx].getPid())) fba.objectives[fba.activeObjIdx].operation = 'maximize' time.sleep(2) objO = fba.objectives[fba.activeObjIdx].operation.lower() objO = objO[0].upper() + objO[1:] FF.write('%s\n' % objO) objStr = '%s_objf: ' % fba.objectives[fba.activeObjIdx].getPid() for fObj in fba.objectives[fba.activeObjIdx].fluxObjectives: sign = None nc = 0.0 try: nc = float(fObj.coefficient) except ValueError: if _HAVE_SYMPY_: nc = sympy.Rational(fObj.coefficient).evalf() else: raise ValuError( 'Invalid coefficient (%s) detected in fluxObjective %s' % (fObj.coefficient, fObj.getPid())) if nc >= 0.0: sign = '+' else: sign = '-' if use_rational: objStr += ' %s %s%s%s' % (sign, sympy.Rational(abs(nc)), multisymb, fObj.reaction) else: objStr += ' %s %s%s%s' % (sign, abs(nc), multisymb, fObj.reaction) else: objStr = '\n\\\\ No objectives defined\n\n' FF.write('%s\n' % objStr) CONST = BuildLPConstraints(fba, use_rational) CONST.seek(0) if fba.user_constraints != None and len(fba.user_constraints) > 0: UCONST = BuildLPUserConstraints(fba, use_rational=False) UCONST.seek(0) BOUNDS = BuildLPFluxBounds(fba, use_rational) BOUNDS.seek(0) if __DEBUG__: print(CONST.read()); CONST.seek(0) if __DEBUG__: print(BOUNDS.read()); BOUNDS.seek(0) FF.write('\nSubject To\n') FF.write(CONST.read()) if fba.user_constraints != None and len(fba.user_constraints) > 0: FF.write('\\\\UserConstraints\n') FF.write(UCONST.read()) FF.write('\nBounds\n') FF.write(BOUNDS.read()) FF.write('\nEND\n') FF.close() if not quiet: print('writeModelLP has written a file to {}'.format(FNAME)) return FNAME
[docs]def BuildHformatFluxBounds(fba, infinity_replace=None, use_rational=False): """ Build and return a csio that contains the flux bounds in H format - *fba* a PySCeS-CBM FBA object - *infinity_replace* [default=None] if defined this is the abs(value) of +-<infinity> """ LBs = {} UBs = {} for c in fba.flux_bounds: ## print '%s: %s' % (c.getPid(), c.operation) minv = None maxv = None oper = None R = c.reaction if c.operation in ['less','lessEqual']: if infinity_replace != None and numpy.isposinf([c.value])[0]: maxv = infinity_replace elif infinity_replace != None and numpy.isneginf([c.value])[0]: maxv = -infinity_replace else: maxv = c.value ## print 'maxv', maxv oper = c.operation if oper == 'less': oper = '<=' ## print 'LP Bounds operator must be <= not <' else: oper = '<=' if c.operation in ['greater','greaterEqual']: if infinity_replace != None and numpy.isposinf([c.value])[0]: minv = infinity_replace elif infinity_replace != None and numpy.isneginf([c.value])[0]: minv = -infinity_replace else: minv = c.value ## print 'minv', minv oper = c.operation if oper == 'greater': oper = '<=' ## print 'LP Bounds operator must be <= not <' else: oper = '<=' if c.operation == 'equals': if infinity_replace != None and numpy.isposinf([c.value])[0]: maxv = infinity_replace elif infinity_replace != None and numpy.isneginf([c.value])[0]: maxv = -infinity_replace else: maxv = c.value oper = '=' ## print 'maxv2', maxv if maxv != None: if not use_rational: UBs.update({R : float(maxv)}) else: UBs.update({R : maxv}) elif minv != None: if not use_rational: LBs.update({R : float(minv)}) else: LBs.update({R : minv}) if __DEBUG__: print(' ') print(LBs) print(UBs) BsRHS = [] if not use_rational: LBm = numpy.zeros((len(LBs), fba.N.shape[1])) UBm = numpy.zeros((len(UBs), fba.N.shape[1])) else: if _HAVE_SYMPY_: LBm = sympy.zeros(len(LBs), fba.N.shape[1]) UBm = sympy.zeros(len(UBs), fba.N.shape[1]) print('\nUsing rational bounds') else: raise RuntimeError('\nError: SymPy required for rational operations') LBskeys = list(LBs) for lb in range(len(LBskeys)): if not use_rational: LBm[lb, fba.N.col.index(LBskeys[lb])] = 1.0 BsRHS.append(LBs[LBskeys[lb]]) else: LBm[lb, fba.N.col.index(LBskeys[lb])] = sympy.Rational(1.0).limit_denominator(__CBCONFIG__['SYMPY_DENOM_LIMIT']) print(LBs[LBskeys[lb]]) BsRHS.append(sympy.Rational(LBs[LBskeys[lb]]).limit_denominator(__CBCONFIG__['SYMPY_DENOM_LIMIT'])) if __DEBUG__: print(fba.N.col) print(LBm) print(UBm) print(BsRHS) UBskeys = list(UBs) for ub in range(len(UBskeys)): if not use_rational: UBm[ub, fba.N.col.index(UBskeys[ub])] = -1.0 BsRHS.append(-UBs[UBskeys[ub]]) else: UBm[ub, fba.N.col.index(UBskeys[ub])] = sympy.Rational(-1.0).limit_denominator(__CBCONFIG__['SYMPY_DENOM_LIMIT']) BsRHS.append(-sympy.Rational(UBs[UBskeys[ub]]).limit_denominator(__CBCONFIG__['SYMPY_DENOM_LIMIT'])) del LBskeys, UBskeys if __DEBUG__: print(fba.N.col) print(LBm) print(UBm) print(BsRHS) print(' ') return numpy.vstack([LBm, UBm]), BsRHS
[docs]def WriteModelHFormatFBA(fba, work_dir=None, use_rational=False, fullLP=True, format='%s', infinity_replace=None): """ INFO: this method will be deprecated please update your scripts to use \"writeModelHFormatFBA2()\" """ print('\nINFO: this method will be deprecated please update your scripts to use \"writeModelHFormatFBA2()\"\n') time.sleep(1) writeModelHFormatFBA(fba, work_dir, use_rational, fullLP, format, infinity_replace)
[docs]def writeModelHFormatFBA(fba, work_dir=None, use_rational=False, fullLP=True, format='%s', infinity_replace=None): """ Write an FBA-LP in polynomial H-Format file. This version has been replaced by `writeModelHFormatFBA2()` but is kept for backwards compatability. - *fba* a PySCeS-CBM FBA object - *Work_dir* [default=None] the output directory - *use_rational* [default=false] use rational numbers in output (requires sympy) - *fullLP* [default=True] include the default objective function as a maximization target - *format* [default='%s'] the number format string - *infinity_replace* [default=None] if defined this is the abs(value) of +-<infinity> """ print('\nwriteModelHFormatFBA is deprecated please use: writeModelHFormatFBA2\n') if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational IO') M = fba LHS = M.N.array.copy() RHS = [0.0 for e in range(M.N.shape[0])] if __DEBUG__: print(LHS) print(RHS) LHS = numpy.vstack([LHS, -M.N.array.copy()]) RHS += [0.0 for e in range(M.N.shape[0])] if __DEBUG__: print(LHS) print(RHS) BsLHS, BsRHS = BuildHformatFluxBounds(M, infinity_replace=infinity_replace) if __DEBUG__: print(BsLHS) print(BsRHS) LHS = numpy.vstack([LHS, BsLHS]) RHS += BsRHS del BsLHS, BsRHS if __DEBUG__: print(LHS) print(RHS) if not use_rational: name = M.getPid().replace('.xml', '') + '.ine' else: name = M.getPid().replace('.xml', '') + '_r.ine' RHS = numpy.array(RHS,'d') RHS.shape = (len(RHS), 1) if __DEBUG__: print(RHS) ## LP = numpy.hstack([LHS, RHS]) OBJ_FUNC = numpy.zeros(LHS.shape[1]+1) for j in range(LHS.shape[1]): # first objective function, first flux objective if __DEBUG__: print(M.objectives[0].fluxObjectives[0].reaction, M.N.col[j]) if M.objectives[0].fluxObjectives[0].reaction == M.N.col[j]: OBJ_FUNC[j] = float(M.objectives[0].fluxObjectives[0].coefficient) if __DEBUG__: print(OBJ_FUNC) # for Ax >= B Hformat wants -B A >= 0 LP = numpy.hstack([-RHS, LHS]) OBJ_FUNC = numpy.hstack([-OBJ_FUNC[-1], OBJ_FUNC[:-1]]) if __DEBUG__: print(OBJ_FUNC) print(LP) del LHS, RHS if work_dir == None: Fname = name else: assert os.path.exists(work_dir), '\nJanee ...' Fname = os.path.join(work_dir, name) F = file(Fname, 'w') F.write('* %s\n\n' % name) F.write('H-representation\n\nbegin\n') NUM_TYPE = 'real' if use_rational: NUM_TYPE = 'rational' F.write('%s %s %s\n' % (LP.shape[0], LP.shape[1], NUM_TYPE)) strW = format+' ' for r in range(LP.shape[0]): for c in range(LP.shape[1]): if not use_rational: if LP[r,c] == 0.0 or LP[r,c] == -0.0: LP[r,c] = 0.0 F.write(strW % LP[r,c]) else: F.write('%s ' % sympy.Rational(format % LP[r,c])) F.write('\n') if fullLP: F.write('end\nlponly\n') F.write('maximize\n') for o in OBJ_FUNC: if not use_rational: F.write(strW % o) else: F.write('%s ' % sympy.Rational(format % o)) else: F.write('end\n') F.write('\n') F.close() F = file(Fname.replace('.ine','')+'.columns.txt', 'w') for j in range(M.N.array.shape[1]): F.write('%s,%s\n' % (j, M.N.col[j])) F.write('\n') F.close() return Fname
[docs]def WriteModelHFormatFBA2(fba, fname=None, work_dir=None, use_rational=False, fullLP=True, format='%s', infinity_replace=None): """ INFO: this method will be deprecated please update your scripts to use \"writeModelHFormatFBA2()\" """ print('\nINFO: this method will be deprecated please update your scripts to use \"writeModelHFormatFBA2()\"\n') time.sleep(1) writeModelHFormatFBA2(fba, fname, work_dir, use_rational, fullLP, format, infinity_replace)
[docs]def writeModelHFormatFBA2(fba, fname=None, work_dir=None, use_rational=False, fullLP=True, format='%s', infinity_replace=None): """ Write an FBA-LP in polynomial H-Format file. This is an improved version of `WriteModelHFormatFBA()` which it replaces. Note that if a SymPy matrix is used as input then use_rational is automatically enabled. - *fba* a PySCeS-CBM FBA object - *fname* [default=None] the output filename, fba.getPid() if not defined - *Work_dir* [default=None] the output directory - *use_rational* [default=false] use rational numbers in output (requires sympy) - *fullLP* [default=True] include the default objective function as a maximization target - *format* [default='%s'] the number format string - *infinity_replace* [default=None] if defined this is the abs(value) of +-<infinity> """ use_rational_old = use_rational if _HAVE_SYMPY_ and fba.N.__array_type__ == sympy.MutableDenseMatrix: use_rational = True print('INFO: using rational matrix') M = fba LHS = M.N.array.copy() #RHS = [0.0 for e in range(M.N.shape[0])] RHS = M.N.RHS.tolist() if __DEBUG__: print(LHS) print(RHS) LHS = numpy.vstack([LHS, -M.N.array.copy()]) tmp = -M.N.RHS RHS += tmp.tolist() del tmp if __DEBUG__: print(LHS) print(RHS) BsLHS, BsRHS = BuildHformatFluxBounds(M, infinity_replace=infinity_replace, use_rational=use_rational) if use_rational: BsLHS = sympy.Matrix(BsLHS) den_lim = __CBCONFIG__['SYMPY_DENOM_LIMIT'] BsRHS = [sympy.Rational(i).limit_denominator(den_lim) for i in BsRHS] for r in range(BsLHS.shape[0]): for c in range(BsLHS.shape[1]): BsLHS[r,c] = sympy.Rational(BsLHS[r,c]).limit_denominator(den_lim) if __DEBUG__: print(BsLHS) print(BsRHS) LHS = numpy.vstack([LHS, BsLHS]) RHS += BsRHS del BsLHS, BsRHS if __DEBUG__: print(LHS) print(RHS) RHS = numpy.array(RHS) RHS.shape = (len(RHS), 1) if __DEBUG__: print(RHS) ## LP = numpy.hstack([LHS, RHS]) OBJ_FUNC = numpy.zeros(LHS.shape[1]+1) objIdx = M.activeObjIdx for j in range(LHS.shape[1]): for fo in range(len(M.objectives[objIdx].getFluxObjectiveReactions())): if M.objectives[objIdx].fluxObjectives[fo].reaction == M.N.col[j]: print(M.objectives[objIdx].fluxObjectives[fo].reaction, M.N.col[j]) OBJ_FUNC[j] = float(M.objectives[objIdx].fluxObjectives[fo].coefficient) ## print OBJ_FUNC # for Ax >= B Hformat wants -B A >= 0 LP = numpy.hstack([-RHS, LHS]) OBJ_FUNC = numpy.hstack([-OBJ_FUNC[-1], OBJ_FUNC[:-1]]) if __DEBUG__: print(OBJ_FUNC) print(LP) del LHS, RHS if work_dir != None: assert os.path.exists(work_dir), '\nJanee ...' fname = os.path.join(work_dir, fname) if fname == None: fname = M.getPid().replace('.xml', '') if not use_rational: fname += '.ine' else: fname += '_r.ine' F = file(fname, 'w') F.write('* %s\n' % os.path.split(fname)[-1]) F.write('H-representation\nbegin\n') NUM_TYPE = 'real' if use_rational: NUM_TYPE = 'rational' F.write('%s %s %s\n' % (LP.shape[0], LP.shape[1], NUM_TYPE)) strW = format+' ' for r in range(LP.shape[0]): for c in range(LP.shape[1]): if use_rational or use_rational_old: F.write('%s ' % LP[r,c]) else: if LP[r,c] == 0.0 or LP[r,c] == -0.0: LP[r,c] = 0.0 F.write(strW % LP[r,c]) F.write('\n') if fullLP: F.write('end\nlponly\n') F.write('maximize\n') # check if Hformat has a minimize kw for o in OBJ_FUNC: #if not use_rational: #F.write(strW % o) #else: #F.write('%s ' % sympy.Rational(format % o)) if use_rational or use_rational_old: F.write('%s ' % sympy.Rational(format % o)) else: if o == 0.0 or o == -0.0: LP[r,c] = 0.0 F.write(strW % o) # then we can use use this ## F.write('%s\n' % M.objectives[M.activeObjIdx].operation) ## if M.activeObjIdx].operation == 'maximize': ## for o in OBJ_FUNC: ## if not use_rational: ## F.write(strW % o) ## else: ## F.write('%s ' % sympy.Rational(format % o)) ## else: ## for o in OBJ_FUNC: ## o = -o ## if not use_rational: ## F.write(strW % o) ## else: ## F.write('%s ' % sympy.Rational(format % o)) else: F.write('end\n') F.write('\n') F.close() F = file(fname.replace('.ine','')+'.columns.txt', 'w') for j in range(M.N.array.shape[1]): F.write('%s,%s\n' % (j, M.N.col[j])) F.write('\n') F.close() return fname
[docs]def writeStoichiometricMatrix(fba, fname=None, work_dir=None, use_rational=False, fullLP=True, format='%s', infinity_replace=None): """ Write an FBA-LP in polynomial H-Format file. This is an improved version of `WriteModelHFormatFBA()` which it replaces but is kept for backwards compatability. - *fba* a PySCeS-CBM FBA object - *fname* [default=None] the output filename, fba.getPid() if not defined - *Work_dir* [default=None] the output directory - *use_rational* [default=false] use rational numbers in output (requires sympy) - *fullLP* [default=True] include the default objective function as a maximization target - *format* [default='%s'] the number format string - *infinity_replace* [default=None] if defined this is the abs(value) of +-<infinity> """ if not _HAVE_SYMPY_ and use_rational: use_rational = False print('Warning: install Sympy for rational IO') M = fba LHS = M.N.array.copy() RHS = [0.0 for e in range(M.N.shape[0])] if __DEBUG__: print(LHS) print(RHS) #LHS = numpy.vstack([LHS, -M.N.array.copy()]) RHS += [0.0 for e in range(M.N.shape[0])] if __DEBUG__: print(LHS) print(RHS) #BsLHS, BsRHS = BuildHformatFluxBounds(M, infinity_replace=infinity_replace) if __DEBUG__: print(BsLHS) print(BsRHS) #LHS = numpy.vstack([LHS, BsLHS]) #RHS += BsRHS #del BsLHS, BsRHS if __DEBUG__: print(LHS) print(RHS) #RHS = numpy.array(RHS,'d') #RHS.shape = (len(RHS), 1) if __DEBUG__: print(RHS) ## LP = numpy.hstack([LHS, RHS]) #OBJ_FUNC = numpy.zeros(LHS.shape[1]+1) objIdx = M.activeObjIdx #for j in range(LHS.shape[1]): #for fo in range(len(M.objectives[objIdx].getFluxObjectiveReactions())): #if M.objectives[objIdx].fluxObjectives[fo].reaction == M.N.col[j]: #print(M.objectives[objIdx].fluxObjectives[fo].reaction, M.N.col[j]) #OBJ_FUNC[j] = float(M.objectives[objIdx].fluxObjectives[fo].coefficient) ### print OBJ_FUNC # for Ax >= B Hformat wants -B A >= 0 #LP = numpy.hstack([-RHS, LHS]) #OBJ_FUNC = numpy.hstack([-OBJ_FUNC[-1], OBJ_FUNC[:-1]]) LP = LHS if __DEBUG__: print(OBJ_FUNC) print(LP) del LHS, RHS if work_dir != None: assert os.path.exists(work_dir), '\nJanee ...' fname = os.path.join(work_dir, fname) if fname == None: fname = M.getPid().replace('.xml', '') if not use_rational: fname += '.ine' else: fname += '_r.ine' F = file(fname, 'w') #F.write('* %s\n' % os.path.split(fname)[-1]) #F.write('H-representation\nbegin\n') #NUM_TYPE = 'real' #if use_rational: #NUM_TYPE = 'rational' #F.write('%s %s %s\n' % (LP.shape[0], LP.shape[1], NUM_TYPE)) strW = format+' ' for r in range(LP.shape[0]): for c in range(LP.shape[1]): if not use_rational: if LP[r,c] == 0.0 or LP[r,c] == -0.0: LP[r,c] = 0.0 F.write(strW % LP[r,c]) else: ## print LP[r,c] F.write('%s ' % sympy.Rational(format % LP[r,c])) F.write('\n') #if fullLP: #F.write('end\nlponly\n') #F.write('maximize\n') # check if Hformat has a minimize kw #for o in OBJ_FUNC: #if not use_rational: #F.write(strW % o) #else: #F.write('%s ' % sympy.Rational(format % o)) # then we can use use this ## F.write('%s\n' % M.objectives[M.activeObjIdx].operation) ## if M.activeObjIdx].operation == 'maximize': ## for o in OBJ_FUNC: ## if not use_rational: ## F.write(strW % o) ## else: ## F.write('%s ' % sympy.Rational(format % o)) ## else: ## for o in OBJ_FUNC: ## o = -o ## if not use_rational: ## F.write(strW % o) ## else: ## F.write('%s ' % sympy.Rational(format % o)) #else: #F.write('end\n') #F.write('\n') F.close() F = file(fname.replace('.ine','')+'.columns.txt', 'w') for j in range(M.N.array.shape[1]): F.write('%s,%s\n' % (j, M.N.col[j])) F.write('\n') F.close() return fname
def writeListToLP(fname, obj=None, const=None, bnds=None, work_dir=None, objtype='maximize'): if work_dir == None: work_dir = os.getcwd() F = file(os.path.join(work_dir, fname+'.lp'), 'w') F.write("\\\\ %s\n" % fname) objtype = objtype.lower() if objtype == 'max': objtype = 'maximize' if objtype == 'min': objtype = 'minimize' if objtype in ['maximise', 'minimise']: objtype = objtype.replace('se','ze') assert objtype in ['maximize', 'minimize'], "\nobjtype must be ['maximize', 'minimize'] not %s" % objtype if obj != None: if objtype == 'maximize': F.write('\nMaximize\n') elif objtype == 'minimize': F.write('\nMinimize\n') for o in obj: F.write(' %s\n' % o) if const != None: F.write('\nSubject to\n') for c in const: F.write('%s\n' % c) if bnds != None: F.write('\nBounds\n') for b in bnds: F.write('%s\n' % b) F.write('\nEND\n') F.close() print('LP written to: {}.lp'.format(os.path.join(work_dir, fname))) return os.path.join(work_dir, fname+'.lp') def writeMinDistanceLP(fname, fbas, work_dir=None, ignoreDistance=[], with_protein_cost=False, constraint_mode='strict', moma=False): if work_dir == None: work_dir = os.getcwd() fC = [] objFname = '' for l in fbas: fC.append(len(l.reactions)) objFname += l.prefix fC = numpy.array(fC) if __DEBUG__: print(fC) print((fC == fC[0])) print(numpy.alltrue((fC == fC[0]))) if moma: assert len(fbas) == 2, '\nMOMA only defined for two inputs' else: assert numpy.alltrue((fC == fC[0])), '\nModels must have the same number of fluxes\n!' conL = [] # model flux_bounds initial_cnstr = [] if moma: initial_cnstr.append(BuildLPConstraintsStrict(fbas[0])) mt = fbas[1].clone() mt.setPrefix('mt_', 'all') initial_cnstr.append(BuildLPConstraintsStrict(mt)) del mt else: for f in fbas: if constraint_mode == 'math': initial_cnstr.append(BuildLPConstraintsMath(f)) elif constraint_mode == 'relaxed': initial_cnstr.append(BuildLPConstraintsRelaxed(f)) else: initial_cnstr.append(BuildLPConstraintsStrict(f)) for ib in initial_cnstr: ib.seek(0) for l in ib: conL.append(l.strip()) conL.append(' ') del initial_cnstr bndL = [] # model bounds initial_bnds = [] for f in fbas: initial_bnds.append(BuildLPFluxBounds(f)) for ib in initial_bnds: ib.seek(0) for l in ib: bndL.append(l.strip()) bndL.append(' ') del initial_bnds artVar = [] artVarX = [] #ignoreDistance = [] Combi = CBTools.ComboGen() Cnumber = 2 Cdata = '' unique_combinations = None for x in range(len(fbas)): Cdata += '%s' % x if __DEBUG__: print(Cdata) Combi.uniqueCombinations(Cdata, Cnumber, temp=[]) Combi.numberifyComb2Int() unique_combinations = Combi.combo_int if __DEBUG__: print('Data ({}):\n{}\n'.format(Cnumber, Cdata)) print('UniqueCombinStr:\n{}'.format(Combi.combo)) print('UniqueCombinations:\n{}'.format(unique_combinations)) zbase = 0 combcount = 1 for uq in unique_combinations: if __DEBUG__: print(uq) MD0 = fbas[uq[0]] MD0 = fbas[uq[0]] if moma: MD1 = fbas[uq[1]].clone() MD1.setPrefix('mt_', 'all') else: MD1 = fbas[uq[1]] RiD1 = MD0.getReactionIds() RiD2 = MD1.getReactionIds() for s in range(len(MD0.reactions)): if __DEBUG__: print(RiD1[s], RiD2[s]) if RiD1[s] not in ignoreDistance: if moma: ridx2 = RiD2.index(MD1.prefix+RiD1[s]) else: ridx2 = s ## av = 'z%s' % (zbase+s+1) av = 'zvar%s%s' % (combcount,RiD1[s].replace(MD0.prefix,'')) c1 = '%sa: %s - %s - %s <= 0.0' % (av, RiD1[s], RiD2[ridx2], av) ## c1 = '%s - %s - %s <= 0.0' % (RiD1[s], RiD2[s], av) c2 = '%sb: %s - %s + %s >= 0.0' % (av, RiD1[s], RiD2[ridx2], av) ## c2 = '%s - %s + %s >= 0.0' % (RiD1[s], RiD2[s], av) # add the protein cost if with_protein_cost: av = '%s %s' % (MD0.reactions[s].annotation['CBM_PEPTIDE_COST'], av) artVar.append(av) conL.append(c1) conL.append(c2) combcount += 1 zbase += len(MD0.reactions) objS = '%smulti: ' % objFname vcntr = 0 for v in artVar: objS += '%s + ' % v vcntr += 1 if vcntr >= 500: objS += '\n' vcntr = 0 objS = objS[:-3] if len(artVarX) >= 1: objS += ' \\* Ignored: ' for o in artVarX: objS += '%s ' % o objS += '*\\\n' objL = [objS] objFcnstr = [' '] assert len(f.objectives[f.activeObjIdx].getFluxObjectiveReactions()) == 1, "\nOnly single fluxObjectives dealt with at this time" if moma: f = fbas[0] objFcnstr.append('C_%s: %s >= %f' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\ f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], f.objectives[f.activeObjIdx].value)) else: for f in fbas: objFcnstr.append('C_%s: %s >= %f' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\ f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], f.objectives[f.activeObjIdx].value)) conL = conL + objFcnstr if __DEBUG__: print(objL) for c in conL: print(c) print(' ') for b in bndL: print(b) F = file(os.path.join(work_dir, fname+'.lp'), 'w') header = '\\\\ MultiInputMinimization: ' for f in fbas: header += '%s, ' % f.getPid() F.write('%s\n' % header[:-2]) F.write('\nMinimize\n') for o in objL: F.write('%s \n' % o) F.write('\nSubject to\n') for c in conL: F.write('%s \n' % c) F.write('\nBounds\n') for b in bndL: F.write(' %s \n' % b) F.write('END\n\n') F.close() print('LP written to: {}.lp'.format(os.path.join(work_dir, fname))) return os.path.join(work_dir, fname+'.lp') def writeMinDistanceLP_absL1(fname, fbas, work_dir=None, ignoreDistance=[], bigM=500, with_protein_cost=False, constraint_mode='strict', moma=False): if work_dir == None: work_dir = os.getcwd() fC = [] objFname = '' for l in fbas: fC.append(len(l.reactions)) objFname += l.prefix fC = numpy.array(fC) if __DEBUG__: print(fC) print((fC == fC[0])) print(numpy.alltrue((fC == fC[0]))) if moma: assert len(fbas) == 2, '\nMOMA only defined for two inputs' else: assert numpy.alltrue((fC == fC[0])), '\nModels must have the same number of fluxes\n!' conL = [] # model flux_bounds initial_cnstr = [] if moma: initial_cnstr.append(BuildLPConstraintsStrict(fbas[0])) mt = fbas[1].clone() mt.setPrefix('mt_', 'all') initial_cnstr.append(BuildLPConstraintsStrict(mt)) del mt else: for f in fbas: if constraint_mode == 'math': initial_cnstr.append(BuildLPConstraintsMath(f)) elif constraint_mode == 'relaxed': initial_cnstr.append(BuildLPConstraintsRelaxed(f)) else: initial_cnstr.append(BuildLPConstraintsStrict(f)) for ib in initial_cnstr: ib.seek(0) for l in ib: conL.append(l.strip()) conL.append(' ') del initial_cnstr bndL = [] # model bounds initial_bnds = [] bigMS = [] for f in fbas: initial_bnds.append(BuildLPFluxBounds(f)) ## bigMS.append(max([abs(float(v.value)) for v in f.reactions])) for ib in initial_bnds: ib.seek(0) for l in ib: bndL.append(l.strip()) bndL.append(' ') del initial_bnds Combi = CBTools.ComboGen() Cnumber = 2 Cdata = '' unique_combinations = None for x in range(len(fbas)): Cdata += '%s' % x if __DEBUG__: print(Cdata) Combi.uniqueCombinations(Cdata, Cnumber, temp=[]) Combi.numberifyComb2Int() unique_combinations = Combi.combo_int if __DEBUG__: print('Data ({}):\n{}\n'.format(Cnumber, Cdata)) print('UniqueCombinStr:\n{}'.format(Combi.combo)) print('UniqueCombinations:\n{}'.format(unique_combinations)) print(bigMS) print('\nbigM = ', bigM, '\n') ## ILPMETHOD = 'SK' # steven ILPMETHOD = 'GK' # gunnar zbase = 0 combcount = 1 artVar = [] artVarX = [] #ignoreDistance = [] boolVars = [] for uq in unique_combinations: if __DEBUG__: print(uq) MD0 = fbas[uq[0]] if moma: MD1 = fbas[uq[1]].clone() MD1.setPrefix('mt_', 'all') else: MD1 = fbas[uq[1]] RiD1 = MD0.getReactionIds() RiD2 = MD1.getReactionIds() for s in range(len(MD0.reactions)): if __DEBUG__: print(RiD1[s], RiD2[s]) if RiD1[s] not in ignoreDistance: if moma: ridx2 = RiD2.index(MD1.prefix+RiD1[s]) else: ridx2 = s Var1 = RiD1[s] bVar1 = 'xvar_%s' % Var1 absVar1 = 'absL_%s' % Var1 Var2 = RiD2[ridx2] bVar2 = 'xvar_%s' % Var2 absVar2 = 'absL_%s' % Var2 c0a = '\n' c0a += '%s - %s >= 0\n' % (absVar1, Var1) c0a += '%s + %s >= 0\n' % (absVar1, Var1) if ILPMETHOD == 'GK': # gunnar c0a += '%s + %s - %s %s <= 0\n' % (absVar1, Var1, bigM, bVar1) c0a += '%s - %s + %s %s <= %s\n' % (absVar1, Var1, bigM, bVar1, bigM) elif ILPMETHOD == 'SK': # steven c0a += '%s - %s - %s %s <= 0\n' % (absVar1, Var1, bigM, bVar1) c0a += '%s + %s + %s %s <= %s\n' % (absVar1, Var1, bigM, bVar1, bigM) c0b = '\n' c0b += '%s - %s >= 0\n' % (absVar2, Var2) c0b += '%s + %s >= 0\n' % (absVar2, Var2) if ILPMETHOD == 'GK': # gunnar c0b += '%s + %s - %s %s <= 0\n' % (absVar2, Var2, bigM, bVar2) c0b += '%s - %s + %s %s <= %s\n' % (absVar2, Var2, bigM, bVar2, bigM) elif ILPMETHOD == 'SK': # steven c0b += '%s - %s - %s %s <= 0\n' % (absVar2, Var2, bigM, bVar2) c0b += '%s + %s + %s %s <= %s\n' % (absVar2, Var2, bigM, bVar2, bigM) c0 = c0a + c0b av = 'zvar%s%s' % (combcount,Var1.replace(MD0.prefix,'')) if bVar1 not in boolVars: boolVars.append(bVar1) if bVar2 not in boolVars: boolVars.append(bVar2) c1 = '%sa: %s - %s - %s <= 0.0' % (av, absVar1, absVar2, av) c2 = '%sb: %s - %s + %s >= 0.0' % (av, absVar1, absVar2, av) # add the protein cost if with_protein_cost: av = '%s %s' % (MD0.reactions[s].annotation['CBM_PEPTIDE_COST'], av) artVar.append(av) conL.append(c0) conL.append(c1) conL.append(c2) combcount += 1 zbase += len(MD0.reactions) objS = '%smulti: ' % objFname vcntr = 0 for v in artVar: objS += '%s + ' % v vcntr += 1 if vcntr >= 500: objS += '\n' vcntr = 0 objS = objS[:-3] if len(artVarX) >= 1: objS += ' \\* Ignored: ' for o in artVarX: objS += '%s ' % o objS += '*\\\n' objL = [objS] objFcnstr = [' '] assert len(f.objectives[f.activeObjIdx].getFluxObjectiveReactions()) == 1, "\nOnly single fluxObjectives dealt with at this time" if moma: #OFvalue = f.objectives[f.activeObjIdx].value f = fbas[0] objFcnstr.append('C_%s: %s >= %f' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\ f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], f.objectives[f.activeObjIdx].value)) else: for f in fbas: objFcnstr.append('C_%s: %s >= %f' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\ f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], f.objectives[f.activeObjIdx].value)) conL = conL + objFcnstr if __DEBUG__: print(objL) for c in conL: print(c) print(' ') for b in bndL: print(b) F = file(os.path.join(work_dir, fname+'.lp'), 'w') header = '\\\\ MultiInputMinimization: ' for f in fbas: header += '%s, ' % f.getPid() F.write('%s\n' % header[:-2]) F.write('\nMinimize\n') for o in objL: F.write('%s \n' % o) F.write('\nSubject to\n') for c in conL: F.write('%s \n' % c) F.write('\nBounds\n') for b in bndL: F.write(' %s \n' % b) if len(boolVars) > 0: F.write('Binary\n') for b in boolVars: F.write(' %s \n' % b) F.write('\nEND\n') F.close() print('LP written to: {}.lp'.format(os.path.join(work_dir, fname))) return os.path.join(work_dir, fname+'.lp')
[docs]def writeMinDistanceLPwithCost(fname, fbas, work_dir=None, ignoreDistance=[], constraint_mode='strict'): """ For backwards compatability only """ print("\n\n**********\nDeprecation warning!\nPlease use writeMinDistanceLP(with_protein_cost=True) instead of writeMinDistanceLPwithCost()\n\n**********\n") writeMinDistanceLP(fname, fbas, work_dir=work_dir, ignoreDistance=ignoreDistance, with_protein_cost=True, constraint_mode=constraint_mode)
# CAN GO SOON JUST PUTTING IN FOR SVN SYNCH """ def writeMinDistanceLPwithCost(fname, fbas, work_dir=None, ignoreDistance=[], constraint_mode='strict'): if work_dir == None: work_dir = os.getcwd() fC = [] objFname = '' for l in fbas: fC.append(len(l.reactions)) objFname += l.prefix fC = numpy.array(fC) if __DEBUG__: print fC print (fC == fC[0]) print numpy.alltrue((fC == fC[0])) assert numpy.alltrue((fC == fC[0])), '\nModels must have the same number of fluxes\n!' conL = [] # model flux_bounds initial_cnstr = [] for f in fbas: if constraint_mode == 'math': initial_cnstr.append(BuildLPConstraintsMath(f)) elif constraint_mode == 'relaxed': initial_cnstr.append(BuildLPConstraintsRelaxed(f)) else: initial_cnstr.append(BuildLPConstraintsStrict(f)) for ib in initial_cnstr: ib.seek(0) for l in ib: conL.append(l.strip()) conL.append(' ') del initial_cnstr bndL = [] # model bounds initial_bnds = [] for f in fbas: initial_bnds.append(BuildLPFluxBounds(f)) for ib in initial_bnds: ib.seek(0) for l in ib: bndL.append(l.strip()) bndL.append(' ') del initial_bnds artVar = [] artVarX = [] ignoreDistance = [] Combi = CBTools.ComboGen() Cnumber = 2 Cdata = '' unique_combinations = None for x in range(len(fbas)): Cdata += '%s' % x if __DEBUG__: print Cdata Combi.uniqueCombinations(Cdata, Cnumber, temp=[]) Combi.numberifyComb2Int() unique_combinations = Combi.combo_int if __DEBUG__: print 'Data (%s):\n%s\n' % (Cnumber, Cdata) print 'UniqueCombinStr:\n%s' % Combi.combo print 'UniqueCombinations:\n%s' % unique_combinations zbase = 0 combcount = 1 for uq in unique_combinations: if __DEBUG__: print uq RiD1 = fbas[uq[0]].getReactionIds() RiD2 = fbas[uq[1]].getReactionIds() for s in range(len(fbas[uq[0]].reactions)): if __DEBUG__: print RiD1[s], RiD2[s] if RiD1[s] not in ignoreDistance: ## av = 'z%s' % (zbase+s+1) ## av = 'zvar%s' % (RiD1[s].replace(fbas[uq[0]].prefix,'')) av = 'zvar%s%s' % (combcount,RiD1[s].replace(fbas[uq[0]].prefix,'')) c1 = '%sa: %s - %s - %s <= 0.0' % (av, RiD1[s], RiD2[s], av) ## c1 = '%s - %s - %s <= 0.0' % (RiD1[s], RiD2[s], av) c2 = '%sb: %s - %s + %s >= 0.0' % (av, RiD1[s], RiD2[s], av) ## c2 = '%s - %s + %s >= 0.0' % (RiD1[s], RiD2[s], av) conL.append(c1) conL.append(c2) ## print RiD1[s], RiD2[s], av ## print fbas[uq[0]].reactions[s].getPid(), fbas[uq[1]].reactions[s].getPid() ## print fbas[uq[0]].reactions[s].annotation['CBM_PEPTIDE_COST'], fbas[uq[1]].reactions[s].annotation['CBM_PEPTIDE_COST'] av = '%s %s' % (fbas[uq[0]].reactions[s].annotation['CBM_PEPTIDE_COST'], av) ## print RiD1[s], RiD2[s], av artVar.append(av) else: # THIS MAY BE A BUG ## av = 'z%s' % (zbase+s+1) ## av = 'zvar%s' % (RiD1[s].replace(fbas[uq[0]].prefix,'')) av = 'zvar%s%s' % (combcount,RiD1[s].replace(fbas[uq[0]].prefix,'')) artVarX.append(av) c1 = '\\* %sa: %s - %s - %s <= 0.0 *\\' % (av, RiD1[s], RiD2[s], av) ## c1 = '\\* %s - %s - %s <= 0.0 *\\' % (RiD1[s], RiD2[s], av) c2 = '\\* %sb: %s - %s + %s >= 0.0 *\\' % (av, RiD1[s], RiD2[s], av) ## c2 = '\\* %s - %s + %s >= 0.0 *\\' % (RiD1[s], RiD2[s], av) conL.append(c1) conL.append(c2) combcount += 1 zbase += len(fbas[uq[0]].reactions) objS = '%smulti: ' % objFname vcntr = 0 for v in artVar: objS += '%s + ' % v vcntr += 1 if vcntr >= 500: objS += '\n' vcntr = 0 objS = objS[:-3] if len(artVarX) >= 1: objS += ' \\* Ignored: ' for o in artVarX: objS += '%s ' % o objS += '*\\\n' objL = [objS] objFcnstr = [' '] for f in fbas: objFcnstr.append('C_%s: %s >= %f' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\ f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], f.objectives[f.activeObjIdx].value)) conL = conL + objFcnstr if __DEBUG__: print objL for c in conL: print c print ' ' for b in bndL: print b F = file(os.path.join(work_dir, fname+'.lp'), 'w') header = '\\\\ MultiInputMinimization: ' for f in fbas: header += '%s, ' % f.getPid() F.write('%s\n' % header[:-2]) F.write('\nMinimize\n') for o in objL: F.write('%s \n' % o) F.write('\nSubject to\n') for c in conL: F.write('%s \n' % c) F.write('\nBounds\n') for b in bndL: F.write(' %s \n' % b) F.write('END\n\n') F.close() print 'LP written to: %s.lp' % os.path.join(work_dir, fname) return os.path.join(work_dir, fname+'.lp') """
[docs]def writeOptimalSolution(fba, fname, Dir=None, separator=',', only_exchange=False): """ This function writes the optimal solution to file - *fba* an instance of an PySCeSCBM model - *fname* the output filename - *Dir* [default=None] use current directory if not None - *separator* [default=','] the column separator - *only_exchange* [default=False] only output fluxes labelled as exchange reactions """ if Dir != None: assert os.path.exists(Dir), '\nPath does not exist' fname = os.path.join(Dir, fname) if separator == ',': fname_r = fname + '_solution.csv' else: fname_r = fname + '_solution.txt' objName = '' if len(fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()) > 1: for J in fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions(): objName += '%s_' % J objName = objName[:-1] else: objName = fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()[0] try: F = file(fname_r, 'w') except IOError: print('\nCSV file \"{}\" is locked by an external application (probably Excel) please close file and try again (or use a different filename).'.format(fname_r)) return cntr = 0 F.write('%s%s%s%s%s%s\n' % ('ObjectiveFunction',separator,objName,separator,separator,separator)) F.write('\"%s\"%s%s%s%s%s%s%s%s%s\"%s\"%s\"%s\"\n' % ('Reaction',separator,'Value',separator,'LowerBound',separator,'UpperBound',separator,'Reduced cost',separator,'Name',separator,'Gene association')) for r in fba.reactions: GO = False if not only_exchange: GO = True elif only_exchange and r.is_exchange: GO = True if GO: if 'GENE ASSOCIATION' in r.annotation: gene = r.annotation['GENE ASSOCIATION'] else: gene = 'none' bnds = fba.getReactionBounds(r.getPid()) Lbnd = -numpy.inf Ubnd = numpy.inf if bnds != None: if bnds[1] != None: Lbnd = bnds[1] if bnds[2] != None: Ubnd = bnds[2] F.write('\"%s\"%s%s%s%s%s%s%s%s%s\"%s\"%s\"%s\"\n' % (r.getPid(),separator,r.value,separator,Lbnd,separator,Ubnd,separator,r.reduced_cost,separator,r.name,separator,gene)) F.flush() F.close() print('Reactions exported to {}'.format(fname_r))
[docs]def writeModelInfoToFile(fba, fname, Dir=None, separator=',', only_exchange=False, met_type='all'): """ This function writes a CBModel to file - *fba* an instance of an PySCeSCBM model - *fname* the output filename - *Dir* [default=None] use directory if not None - *separator* [default=','] the column separator - *only_exchange* [default=False] only output fluxes labelled as exchange reactions - *type* [default='all'] only output certain type of species: 'all','boundary' or 'variable' """ writeReactionInfoToFile(fba, fname, Dir=Dir, separator=separator, only_exchange=only_exchange) writeSpeciesInfoToFile(fba, fname, Dir=Dir, separator=separator, met_type=met_type)
[docs]def writeReactionInfoToFile(fba, fname, Dir=None, separator=',', only_exchange=False): """ This function writes a CBModel to file - *fba* an instance of an PySCeSCBM model - *fname* the output filename - *Dir* [default=None] use directory if not None - *separator* [default=','] the column separator - *only_exchange* [default=False] only output fluxes labelled as exchange reactions """ if Dir != None: assert os.path.exists(Dir), '\nPath does not exist' fname = os.path.join(Dir, fname) if separator == ',': fname_r = fname + '.rxns.csv' else: fname_r = fname + '.rxns.txt' objName = '' if len(fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()) > 1: for J in fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions(): objName += '%s_' % J objName = objName[:-1] else: objName = fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()[0] try: F = file(fname_r, 'w') except IOError: print('\nOutput file \"{}\" is locked by an external application (probably Excel) please close file and try again (or use a different filename).'.format(fname_r)) return cntr = 0 F.write('%s%s%s%s%s%s\n' % ('ObjectiveFunction',separator,objName,separator,separator,separator)) F.write('\"%s\"%s%s%s%s%s\"%s\"%s\"%s\"%s\"%s\"\n' % ('Reaction',separator,'LowerBound',separator,'UpperBound',separator,'Name',separator,'Equation',separator,'Gene association')) for r in fba.reactions: if not only_exchange: GO = True elif only_exchange and r.is_exchange: GO = True if GO: if 'GENE ASSOCIATION' in r.annotation: gene = r.annotation['GENE ASSOCIATION'] else: gene = 'none' bnds = fba.getReactionBounds(r.getPid()) Lbnd = -numpy.inf Ubnd = numpy.inf if bnds != None: if bnds[1] != None: Lbnd = bnds[1] if bnds[2] != None: Ubnd = bnds[2] if r.reversible: equation = ' %s ' % __CBCONFIG__['REVERSIBLE_SYMBOL'] else: equation = ' %s ' % __CBCONFIG__['IRREVERSIBLE_SYMBOL'] subs = '' prods = '' for rr in r.reagents: if rr.coefficient > 0.0: if abs(rr.coefficient) == 1.0: prods += ' + %s' % (rr.species_ref) else: prods += ' + %s %s' % (abs(rr.coefficient), rr.species_ref) else: if abs(rr.coefficient) == 1.0: subs += ' + %s' % (rr.species_ref) else: subs += ' + %s %s' % (abs(rr.coefficient), rr.species_ref) subs = subs[3:] prods = prods[3:] equation = subs + equation + prods F.write('\"%s\"%s%s%s%s%s\"%s\"%s\"%s\"%s\"%s\"\n' % (r.getPid(),separator,Lbnd,separator,Ubnd,separator,r.name,separator,equation,separator,gene)) GO = False F.flush() F.close() print('Reactions exported to {}'.format(fname_r))
[docs]def writeSpeciesInfoToFile(fba, fname, Dir=None, separator=',', met_type='all'): """ This function writes a CBModel to file - *fba* an instance of an PySCeSCBM model - *fname* the output filename - *Dir* [default=None] use directory if not None - *separator* [default=','] the column separator - *met_type* [default='all'] only output certain type of species: 'all','boundary' or 'variable' """ if Dir != None: assert os.path.exists(Dir), '\nPath does not exist' fname = os.path.join(Dir, fname) if separator == ',': fname_s = fname + '.spec.csv' fname_sr = fname + '.s2r.csv' else: fname_s = fname + '.spec.txt' fname_sr = fname + '.s2r.txt' if met_type not in ['all','boundary','variable']: met_type = 'all' try: F2 = file(fname_s, 'w') F3 = file(fname_sr, 'w') except IOError: print('\nOutput file \"{}\" is locked by an external application (probably Excel) please close file and try again (or use a different filename).'.format(fname_s)) return maxReagentOf = 0 F2.write('\"%s\"%s%s%s%s%s\"%s\"%s\"%s\"%s\"%s\"\n' % ('Species',separator,'Name',separator,'Compartment',separator,'Fixed',separator,'ChemicalFormula',separator,'Charge')) for s in fba.species: if len(s.isReagentOf()) > maxReagentOf: maxReagentOf = len(s.reagent_of) GOS = False if met_type == 'all': GOS = True elif met_type == 'variable' and not s.is_boundary: GOS = True elif met_type == 'boundary' and s.is_boundary: GOS = True if GOS: F2.write('\"%s\"%s%s%s%s%s\"%s\"%s\"%s\"%s\"%s\"\n' % (s.getPid(),separator,s.getName(),separator,s.compartment,separator,s.is_boundary,separator,s.chemFormula,separator,s.charge)) for s in fba.species: srl = len(s.reagent_of) rgOut = '\"%s\"%s' % (s.getPid(), separator) GOS = False if met_type == 'all': GOS = True elif met_type == 'variable' and not s.is_boundary: GOS = True elif met_type == 'boundary' and s.is_boundary: GOS = True if GOS: for rg in range(maxReagentOf): if rg < srl: rgOut += '%s%s' % (s.reagent_of[rg], separator) else: rgOut += '%s' % separator F3.write('%s\n' % rgOut[:-1]) F2.flush() F2.close() F3.flush() F3.close() print('Species exported to {}'.format(fname_s)) print('Reagent map exported to {}'.format(fname_sr))
[docs]def printFBASolution(fba, include_all=False): """ Prints the FBA optimal solution to the screen. - *fba* an FBA model object - *include_all* include all variables """ OFflux = fba.objectives[fba.activeObjIdx].fluxObjectives[0].reaction OFvalue = fba.objectives[fba.activeObjIdx].value OFSense = fba.objectives[fba.activeObjIdx].operation print('\n\n**********\nModel: {}\n\n'.format(fba.getPid())) print('{} objective: {}\nOptimal value: {}\n\n'.format(OFSense,OFflux,OFvalue)) if include_all: for J in fba.reactions: print('{}: {}'.format(J.getPid(), J.value)) print('**********\n')
[docs]def exportModel(fba, fname=None, fmt='lp', work_dir=None, use_rational='both'): """ Export the FBA model in different formats: - *fba* the FBA model - *fname* [default=None] the exported filename if None then `fba.getPid()` is used - *fmt* [default='lp'] the export format can be one of: 'lp' (CPLEX), 'hformat' (Polyhedra), 'all' (both) - *use_rational* [default='both'] if *all* or *hformat* is specified should hformat files be written using rational math or not. The default *both* is the legacy behaviour and writes both. Note that 'hformat' ignores 'fname' and only uses fba.getPid() this is a legacy behaviour """ if work_dir == None: work_dir = os.getcwd() if fmt == 'all' or fmt == 'lp': writeModelLP(fba, work_dir=work_dir, fname=fname) if fmt == 'all' or fmt == 'hformat': if fname == None: fname = fba.getPid().replace('.xml', '') if use_rational == 'both': writeModelHFormatFBA2(fba, fname=fname, work_dir=work_dir) writeModelHFormatFBA2(fba, fname=fname, work_dir=work_dir, use_rational=True) elif use_rational == True: writeModelHFormatFBA2(fba, fname=fname, work_dir=work_dir, use_rational=True) elif use_rational == False: writeModelHFormatFBA2(fba, fname=fname, work_dir=work_dir)
[docs]def writeProteinCostToCSV(fba, fname): """ Writes the protein costs 'CBM_PEPTIDE_COST' annotation toa csv file. - *fba* an instantiated FBA object - *fname* the exported file name """ F = file(fname+'.costs.csv','w') F.write('rid,minL,maxL,avgL,cost\n') for R in fba.reactions: rid = R.getPid() pcost = '' avg_l = '' minL = '' maxL = '' if 'CBM_PEPTIDE_COST' in R.annotation: pcost = R.annotation['CBM_PEPTIDE_COST'] if 'CBM_PEPTIDE_LENGTH_MAX' in R.annotation: maxL = R.annotation['CBM_PEPTIDE_LENGTH_MAX'] if R.annotation['CBM_PEPTIDE_LENGTH_MAX'] == None: avg_l = 1 else: if 'CBM_AVG_PEPTIDE_LENGTH' in R.annotation: avg_l = R.annotation['CBM_AVG_PEPTIDE_LENGTH'] else: avg_l = 0 if 'CBM_PEPTIDE_LENGTH_MIN' in R.annotation: minL = R.annotation['CBM_PEPTIDE_LENGTH_MIN'] F.write('%s,%s,%s,%s,%s\n' % (rid,minL,maxL,avg_l,pcost)) F.flush() F.close() print('Protein costs written to file: {}'.format(fname+'.csv'))
[docs]def WriteFVAtoCSV(id, fva, names, Dir=None, fbaObj=None): """ INFO: this method will be deprecated please update your scripts to use \"writeFVAtoCSV()\" """ print('\nINFO: this method will be deprecated please update your scripts to use \"writeFVAtoCSV()\"\n') time.sleep(1) writeFVAtoCSV(fva, names, id, Dir, fbaObj)
#def writeFVAtoCSV(id, fva, names, Dir=None, fbaObj=None): #""" #INFO: this method will be deprecated please update your scripts to use \"writeFVAtoCSV\" #""" #print('\nINFO: this method will be deprecated please update your scripts to use \"writeFVAtoCSV\"\n') #time.sleep(1) #writeFVAtoCSV(id, fva, names, Dir, fbaObj)
[docs]def writeFVAtoCSV(fvadata, names, fname, Dir=None, fbaObj=None): """ Takes the resuls of a FluxVariabilityAnalysis method and writes it to a nice csv file. Note this method replaces the glpk/cplx_WriteFVAtoCSV methods. - *fvadata* FluxVariabilityAnalysis() OUTPUT_ARRAY - *names* FluxVariabilityAnalysis() OUTPUT_NAMES - *fname* filename_base for the CSV output - *Dir* [default=None] if set the output directory for the csv files - *fbaObj* [default=None] if supplied adds extra model information into the output tables """ if Dir != None: Dir = os.path.join(Dir, fname+'.fva.csv') else: Dir = fname+'.fva.csv' F = file(Dir, 'w') if fbaObj == None: F.write('name,optval,min,max,diff,red cost,minstat,maxstat\n') else: F.write('name,optval,min,max,diff,red cost,minstat,maxstat,"equation","subsystem","gene association","confidence level"\n') for Jidx in range(len(names)): if names[Jidx] != None: name = names[Jidx] optval = fvadata[Jidx][0] rc = fvadata[Jidx][1] min = fvadata[Jidx][2] max = fvadata[Jidx][3] diff = fvadata[Jidx][4] minstat = fvadata[Jidx][5] maxstat = fvadata[Jidx][6] if fbaObj == None: F.write('%s,%s,%s,%s,%s,%s,%s,%s\n' % (name, optval, min, max, diff, rc, minstat, maxstat)) else: xInf = [] Ro = fbaObj.getReaction(name) for k in ['Equation','SUBSYSTEM','GENE ASSOCIATION','Confidence Level']: if k in Ro.annotation: xInf.append(Ro.annotation[k]) else: xInf.append('') F.write('%s,%s,%s,%s,%s,%s,%s,%s,"%s","%s","%s",%s\n' % (name, optval, min, max, diff, rc, minstat, maxstat, xInf[0], xInf[1], xInf[2], xInf[3])) F.flush() F.close() print('FVA results written to: {}'.format(Dir))
[docs]def WriteFVAdata(fva, names, fname, work_dir=None, roundec=None, scale_min=False, appendfile=False, info=None): """ INFO: this method will be deprecated please update your scripts to use \"writeFVAdata()\" """ print('\nINFO: this method will be deprecated please update your scripts to use \"writeFVAdata()\"\n') time.sleep(1) writeFVAdata(fva, names, fname, work_dir, roundec, scale_min, appendfile, info)
[docs]def writeFVAdata(fvadata, names, fname, work_dir=None, roundec=None, scale_min=False, appendfile=False, info=None): """ Takes the resuls of a FluxVariabilityAnalysis method and writes it to a nice csv file. Note this method replaces the glpk/cplx_WriteFVAtoCSV methods. Data is output as a csv file with columns: FluxName, FVA_MIN, FVA_MAX, OPT_VAL, SPAN - *fvadata* FluxVariabilityAnalysis() FVA OUTPUT_ARRAY - *names* FluxVariabilityAnalysis() FVA OUTPUT_NAMES - *fname* filename_base for the CSV output - *work_dir* [default=None] if set the output directory for the csv files - *roundec* [default=None] an integer indicating at which decimal to round off output. Default is no rounding. - *scale_min* [default=False] normalise each flux such that that FVA_MIN = 0.0 - *appendfile* [default=False] instead of opening a new file try and append the data - *info* [default=None] a string added to the results as an extra column, useful with `appendfile` """ if work_dir != None: work_dir = os.path.join(work_dir, fname+'.fvadata.csv') else: work_dir = fname+'.fvadata.csv' if not appendfile: F = file(work_dir, 'w') else: F = file(work_dir, 'a') if info == None: F.write('%s,%s,%s,%s,%s\n' % ('Jid', 'min', 'max', 'optval', 'span')) else: F.write('%s,%s,%s,%s,%s,%s\n' % ('Jid', 'min', 'max', 'optval', 'span', 'info')) for Jidx in range(len(names)): if names[Jidx] != None: name = names[Jidx] max = fvadata[Jidx][3] min = fvadata[Jidx][2] optval = fvadata[Jidx][0] if roundec != None: max = round(max, roundec) min = round(min, roundec) optval = round(optval, roundec) if scale_min: if min > 0.0: max = max-min optval = optval-min min = 0.0 elif min <= 0.0: max = max + abs(min) optval = optval + abs(min) min = 0.0 if info == None: if roundec == None: F.write('%s,%s,%s,%s,%s\n' % (name, min, max, optval, abs(fvadata[Jidx][3] - fvadata[Jidx][2]))) else: F.write('%s,%s,%s,%s,%s\n' % (name, min, max, optval, round(abs(fvadata[Jidx][3] - fvadata[Jidx][2]), roundec))) else: if roundec == None: F.write('%s,%s,%s,%s,%s,%s\n' % (name, min, max, optval, abs(fvadata[Jidx][3] - fvadata[Jidx][2]), info)) else: F.write('%s,%s,%s,%s,%s,%s\n' % (name, min, max, optval, round(abs(fvadata[Jidx][3] - fvadata[Jidx][2]), roundec), info)) F.flush() F.close() print('FVAdata results written to: {}'.format(work_dir))
[docs]def writeSolutions(fname, sols=[], sep=',', extra_output=None, fba=None): """ Write 2 or more solutions where a solution is a dictionary of flux:value pairs: - *fname* the export filename - *sols* a list of dictionaries containing flux:value pairs (e.g. output by cmod.getReactionValues()) - *sep* [default=','] the column separator - *extra_output* [default=None] add detailed information to output e.g. reaction names by giving a CBModel object as an argument to *extra_output*. - *fba* an fba model that canbe used for extra_output """ assert len(sols) >= 2, "\nThere must be two or more solutions to work with" reac_ids = set([]) for s in sols: reac_ids = reac_ids.union(set(s)) reac_ids = list(reac_ids) reac_ids.sort() reac_names = [] ## reac_bnds = [] ## reac_eqns = [] if extra_output != None and fba != None: for r in reac_ids: reac_names.append(fba.getReaction(r).getName()) ## reac_bnds.append(fba.getReactionBounds(r.getPid())) F = file(fname+'.csv', 'w') for r in range(len(reac_ids)): row = '%s%s' % (reac_ids[r], sep) s_str = '' for s in range(len(sols)): if reac_ids[r] in sols[s]: s_str += '%s%s' % (sols[s][reac_ids[r]], sep) else: s_str += '%s%s' % ('\"none\"', sep) row += s_str if extra_output: row += '\"%s\"\n' % reac_names[r] else: row = row[:-1]+'\n' F.write(row) F.flush() F.close() print('\nSolutions written to: \"%s\"\n' % fname)
[docs]def generateBGID(num, prefix): """ Create a BGID generator, which is <prefix><num> where perfix is two letters num is padded to 6 figures - *num* the starting number - *prefix* the two letter prefix """ assert prefix[0].isalpha() and prefix[1].isalpha(), '\nInvalid prefix' while True: num += 1 yield '{}{:0>6}'.format(prefix, num)
[docs]def convertFloatToExcel(num, roundoff): """ Converts a float to Excel compatible "number" - *num* a number - *roundoff* the number of roundoff digits for round() """ ## keep this in sync with convertExcelToFloat if num is None: rval = '' elif numpy.isnan(num): rval = 'nan' elif numpy.isposinf(num): rval = 'inf' elif numpy.isneginf(num): rval = '-inf' else: rval = round(num, roundoff) return rval
[docs]def convertExcelToFloat(num): """ Converts an Excel "number" to a float - *num* a number """ ## keep this in sync with convertFloatToExcel if num == '': rval = None elif num is 'nan': rval = numpy.nan elif num == 'inf': rval = numpy.inf elif num == '-inf': rval = -numpy.inf else: rval = float(num) return rval
[docs]def writeModelToExcel97(fba, filename, roundoff=6): """ Exports the model as an Excel 97 spreadsheet - *fba* a CBMPy model instance - *filename* the filename of the workbook - *roundoff* [default=6] the number of digits to round off to """ if not _HAVE_XLWT_: print('\nERROR: Cannot create Excel file, XLWT package not available (http://pypi.python.org/pypi/xlwt)') return try: F = file('{}.xls'.format(filename), 'wb') F.close() except: print('\nERROR: cannot open file "{}"! Please close workbook before writing!\n'.format('{}.xls'.format(filename))) return fontB = xlwt.Font() fontB.bold = True fontI = xlwt.Font() fontI.italic = True fontBI = xlwt.Font() fontBI.bold = True fontBI.italic = True styleBold = xlwt.XFStyle() styleBold.font = fontB styleBoldC = xlwt.XFStyle() styleBoldC.font = fontB styleBoldC.alignment.horz = xlwt.Alignment.HORZ_CENTER styleItalic = xlwt.XFStyle() styleItalic.font = fontI styleItalicC = xlwt.XFStyle() styleItalicC.font = fontI styleItalicC.alignment.horz = xlwt.Alignment.HORZ_CENTER styleBoldItalic = xlwt.XFStyle() styleBoldItalic.font = fontBI styleBoldItalicC = xlwt.XFStyle() styleBoldItalicC.font = fontBI styleBoldItalicC.alignment.horz = xlwt.Alignment.HORZ_CENTER styleHyper = xlwt.easyxf('font: underline single') styleHyper.alignment.horz = xlwt.Alignment.HORZ_CENTER wb = xlwt.Workbook(encoding='utf-8') wsInf = wb.add_sheet('info') wsSol = wb.add_sheet('solution') wsRe = wb.add_sheet('reactions') wsMet = wb.add_sheet('metabolites') wsSt = wb.add_sheet('network_react') wsStR = wb.add_sheet('network_metab') wsMiriam = wb.add_sheet('miriam') wsComp = wb.add_sheet('compartments') wsGrp = wb.add_sheet('groups') ridx = 0 wsInf.write(ridx, 0, 'id', styleBold) wsInf.write(ridx, 1, fba.getPid()) ridx += 1 wsInf.write(ridx, 0, 'name', styleBold) wsInf.write(ridx, 1, fba.getName()) ridx += 1 wsInf.write(ridx, 0, 'objective', styleBold) try: wsInf.write(ridx, 1, fba.getActiveObjective().getOperation()) except: wsInf.write(ridx, 1, 'undefined') cidx = 2 try: for t_ in fba.getActiveObjective().getFluxObjectiveData(): try: wsInf.write(ridx, cidx, t_[0]) except: wsInf.write(ridx, cidx, 'undefined') try: wsInf.write(ridx, cidx+1, t_[1]) except: wsInf.write(ridx, cidx+1, 'undefined') cidx += 2 except: pass ridx += 1 wsInf.write(ridx, 0, 'obj. value', styleBold) try: wsInf.write(ridx, 1, fba.getActiveObjective().getValue()) except: wsInf.write(ridx, 1, 'undefined') ridx += 1 wsInf.write(ridx, 0, 'modelcreated', styleBold) if fba.DATE_CREATED != None: for t_ in range(len(fba.DATE_CREATED)): wsInf.write(ridx, t_+1, fba.DATE_CREATED[t_]) ridx += 1 wsInf.write(ridx, 0, 'lastmodified', styleBold) if fba.DATE_MODIFIED != None: for t_ in range(len(fba.DATE_MODIFIED)): wsInf.write(ridx, t_+1, fba.DATE_MODIFIED[t_]) ridx += 1 cTime = list(time.gmtime()) wsInf.write(ridx, 0, 'exported', styleBold) for t_ in range(len(cTime)): if t_<= 5: wsInf.write(ridx, t_+1, cTime[t_]) ridx += 1 ccridx = 0 for cr_ in fba.getModelCreators(): cidx = 0 ccridx += 1 wsInf.write(ridx, cidx, 'creator{}'.format(ccridx), styleBold) for ne_ in cr_: cidx += 1 wsInf.write(ridx, cidx, ne_) ridx += 1 ridx += 1 wsInf.write(ridx, 0, xlwt.Formula('HYPERLINK("{}";"{}")'.format('http://cbmpy.sourceforge.net', 'Exported by CBMPy {}'.format(__version__))), styleHyper) ridx += 1 #Scol = ['metabolite', 'value', 'shadow price', '', 'reaction', 'flux', 'reduced cost', 'FVA min', 'FVA max'] if fba.SCALED_REDUCED_COSTS: rcstr = 'scaled reduced cost' else: rcstr = 'reduced cost' Scol = ['reaction', 'flux', 'lower', 'upper', rcstr, 'FVA min', 'FVA max', 'FVA span', 'exchange', 'info', 'stoichiometry'] for s_ in range(len(Scol)): wsSol.write(0, s_, Scol[s_], styleBoldC) # metabolites Mlist = [] Mcols = ['id','name','charge','chemformula','compartment','fixed'] #Mcols = ['id','name','charge','chemformula','compartment','fixed','bgid'] MUcols = [] MiriQual = ["is","isEncodedBy","encodes","hasPart","hasProperty","hasTaxon","hasVersion","isDescribedBy",\ "isHomologTo","isPartOf","isPropertyOf","isVersionOf","occursIn"] wsMiriam.write(0, 0, 'id', styleBoldC) for q_ in range(len(MiriQual)): wsMiriam.write(0, q_+1, MiriQual[q_], styleBold) # get data mcntr = 1 for s_ in range(len(fba.species)): Mdi = {} Mdi['id'] = fba.species[s_].getPid() Mdi['name'] = fba.species[s_].getName() Mdi['charge'] = fba.species[s_].getCharge() Mdi['chemformula'] = fba.species[s_].getChemFormula() Mdi['compartment'] = fba.species[s_].compartment Mdi['fixed'] = fba.species[s_].is_boundary #Mdi['bgid'] = bgGen.next() udata = {} for k_ in fba.species[s_].getAnnotations(): if k_ not in MUcols: MUcols.append(k_) udata[k_] = fba.species[s_].annotation[k_] Mdi['data'] = udata Mdi['reagentof'] = fba.species[s_].isReagentOf() Mlist.append(Mdi) ## TODO: removed for now will appear in its own sheet at some point #wsSol.write(s_+1, 0, Mdi['id']) #Sval = fba.species[s_].getValue() #if Sval == None or numpy.isnan(Sval) or numpy.isinf(Sval): #Sval = '' #wsSol.write(s_+1, 1, Sval) #wsSol.write(s_+1, 2, fba.species[s_].shadow_price) annot = fba.species[s_].getMIRIAMannotations() if annot != None: for k_ in annot: for m_ in annot[k_]: wsMiriam.write(mcntr, 0, Mdi['id']) #if 'identifiers.org' in m_: #wsMiriam.write(mcntr, MiriQual.index(k_)+1, xlwt.Formula('HYPERLINK("{}";"{}")'.format(m_.replace('identifiers.org', 'info.identifiers.org'),\ #m_.rsplit('/',1)[1])), styleHyper) #else: #wsMiriam.write(mcntr, MiriQual.index(k_)+1, m_) wsMiriam.write(mcntr, MiriQual.index(k_)+1, m_) mcntr += 1 # reactions Rlist = [] #Rcols = ['id','name','reversible','lowerbound','upperbound','compartment','bgid'] Rcols = ['id','name','reversible','lowerbound','upperbound','compartment'] RUcols = [] try: AOIDS = fba.getActiveObjective().getFluxObjectiveReactions() except: AOIDS = [] rids = fba.getReactionIds() rids.sort() for r_ in range(len(fba.reactions)): Rdi = {} REAC = fba.getReaction(rids[r_]) Rdi['id'] = REAC.getPid() Rdi['name'] = REAC.getName() Rdi['reversible'] = REAC.reversible Rdi['compartment'] = REAC.compartment #Rdi['bgid'] = bgGen.next() Rdi['lowerbound'] = fba.getReactionLowerBound(Rdi['id']) Rdi['upperbound'] = fba.getReactionUpperBound(Rdi['id']) udata = {} for k_ in REAC.getAnnotations(): if k_ not in RUcols: RUcols.append(k_) udata[k_] = REAC.annotation[k_] Rdi['data'] = udata subs = [] prods = [] for rr_ in REAC.reagents: if rr_.getCoefficient() > 0.0: prods.append((abs(rr_.getCoefficient()), rr_.getSpecies())) else: subs.append((abs(rr_.getCoefficient()), rr_.getSpecies())) Rdi['substrates'] = subs Rdi['products'] = prods Rlist.append(Rdi) cstart = 0 bnds = list(fba.getReactionBounds(Rdi['id'])) # deal with equality if bnds[3] != None: bnds[1] = bnds[3] bnds[2] = bnds[3] if Rdi['id'] not in AOIDS: wsSol.write(r_+1, cstart, Rdi['id']) else: wsSol.write(r_+1, cstart, Rdi['id'], styleBold) rval = REAC.getValue() rval = convertFloatToExcel(rval, roundoff) wsSol.write(r_+1, cstart+1, rval) del rval wsSol.write(r_+1, cstart+2, convertFloatToExcel(bnds[1], roundoff)) wsSol.write(r_+1, cstart+3, convertFloatToExcel(bnds[2], roundoff)) wsSol.write(r_+1, cstart+4, convertFloatToExcel(REAC.reduced_cost, roundoff)) if REAC.fva_min != None: wsSol.write(r_+1, cstart+5, convertFloatToExcel(REAC.fva_min, roundoff)) if REAC.fva_max != None: wsSol.write(r_+1, cstart+6, convertFloatToExcel(REAC.fva_max, roundoff)) if REAC.fva_min != None and REAC.fva_max != None: wsSol.write(r_+1, cstart+7, convertFloatToExcel(REAC.fva_max-REAC.fva_min, roundoff)) if REAC.is_exchange: wsSol.write(r_+1, cstart+8, 'yes') #else: #wsSol.write(r_+1, cstart+8, 'no') wsSol.write(r_+1, cstart+9, xlwt.Formula('HYPERLINK("#reactions!A{}";"{}")'.format(r_+1, 'info')), styleHyper) wsSol.write(r_+1, cstart+10, xlwt.Formula('HYPERLINK("#network_react!A{}";"{}")'.format(((r_+1)*3)-2, 'stoich')), styleHyper) annot = REAC.getMIRIAMannotations() if annot != None: for k_ in annot: for m_ in annot[k_]: wsMiriam.write(mcntr, 0, Rdi['id']) #if 'identifiers.org' in m_: #wsMiriam.write(mcntr, MiriQual.index(k_)+1, xlwt.Formula('HYPERLINK("{}";"{}")'.format(m_.replace('identifiers.org', 'info.identifiers.org'),\ #m_.rsplit('/',1)[1])), styleHyper) #else: #wsMiriam.write(mcntr, MiriQual.index(k_)+1, m_) wsMiriam.write(mcntr, MiriQual.index(k_)+1, m_) mcntr += 1 compHead = ['id', 'name', 'size', 'dimensions', '# species', '# reactions'] compHeadU = [] for c_ in range(len(compHead)): wsComp.write(0, c_, compHead[c_], styleBoldC) ridx = 1 for c_ in fba.compartments: wsComp.write(ridx, 0, c_.getPid()) wsComp.write(ridx, 1, c_.getName()) wsComp.write(ridx, 2, convertFloatToExcel(c_.getSize(), roundoff)) wsComp.write(ridx, 3, c_.getDimensions()) wsComp.write(ridx, 4, len(c_.containsSpecies())) cRe = c_.containsReactions() if len(cRe) != 0: wsComp.write(ridx, 5, len(cRe)) annot = c_.getAnnotations() for k_ in annot: if k_ not in compHeadU: compHeadU.append(k_) for k_ in annot: wsComp.write(ridx, 6+compHeadU.index(k_), annot[k_]) annot = c_.getMIRIAMannotations() if annot != None: for k_ in annot: for m_ in annot[k_]: wsMiriam.write(mcntr, 0, c_.getPid()) wsMiriam.write(mcntr, MiriQual.index(k_)+1, m_) mcntr += 1 ridx += 1 for ud_ in range(0, len(compHeadU)): wsComp.write(0, 6+ud_, compHeadU[ud_], styleItalicC) ridx = 1 del ridx, c_ # write data # write metabolites for m_ in range(len(Mcols)): wsMet.write(0, m_, Mcols[m_], styleBold) for m_ in range(len(MUcols)): wsMet.write(0, len(Mcols)+m_, MUcols[m_], styleItalic) stridx = 0 for s_ in range(len(Mlist)): s = s_+1 wsMet.write(s, Mcols.index('id'), Mlist[s_]['id']) wsMet.write(s, Mcols.index('name'), Mlist[s_]['name']) wsMet.write(s, Mcols.index('charge'), Mlist[s_]['charge']) wsMet.write(s, Mcols.index('chemformula'), Mlist[s_]['chemformula']) wsMet.write(s, Mcols.index('compartment'), Mlist[s_]['compartment']) wsMet.write(s, Mcols.index('fixed'), Mlist[s_]['fixed']) #wsMet.write(s, Mcols.index('bgid'), Mlist[s_]['bgid']) for ud_ in Mlist[s_]['data']: if len(Mlist[s_]['data'][ud_]) < 30000: wsMet.write(s, len(Mcols)+MUcols.index(ud_), Mlist[s_]['data'][ud_]) else: wsMet.write(s, len(Mcols)+MUcols.index(ud_),'Data too long (more than 30000 characters)', styleBold) wsStR.write(stridx, 0, Mlist[s_]['id'], styleBold) cidx = 0 stridx += 1 for ro_ in range(len(Mlist[s_]['reagentof'])): wsStR.write(stridx, cidx, Mlist[s_]['reagentof'][ro_]) cidx += 1 if cidx >= 200: cidx = 0 stridx += 1 stridx += 1 # write reactions for r_ in range(len(Rcols)): wsRe.write(0, r_, Rcols[r_], styleBold) for r_ in range(len(RUcols)): try: wsRe.write(0, len(Rcols)+r_, RUcols[r_], styleItalic) except Exception as why: print(why) print('\nERROR: Exceeded Excel columns') ridx = 0 for r_ in range(len(Rlist)): r = r_+1 wsRe.write(r, Rcols.index('id'), Rlist[r_]['id']) wsRe.write(r, Rcols.index('name'), Rlist[r_]['name']) wsRe.write(r, Rcols.index('reversible'), Rlist[r_]['reversible']) wsRe.write(r, Rcols.index('lowerbound'), convertFloatToExcel(Rlist[r_]['lowerbound'], roundoff)) wsRe.write(r, Rcols.index('upperbound'), convertFloatToExcel(Rlist[r_]['upperbound'], roundoff)) wsRe.write(r, Rcols.index('compartment'), Rlist[r_]['compartment']) #wsRe.write(r, Rcols.index('bgid'), Rlist[r_]['bgid']) for ud_ in Rlist[r_]['data']: if Rlist[r_]['data'][ud_] == None: wsRe.write(r, len(Rcols)+RUcols.index(ud_), '') elif type(Rlist[r_]['data'][ud_]) == str and len(Rlist[r_]['data'][ud_]) >= 30000: wsRe.write(r, len(Rcols)+RUcols.index(ud_), 'Data too long (more than 30000 characters)', styleBold) else: try: wsRe.write(r, len(Rcols)+RUcols.index(ud_), Rlist[r_]['data'][ud_]) except: print('Exceeded Excel columns') wsSt.write(ridx, 0, Rlist[r_]['id'], styleBold) ridx += 1 wsSt.write(ridx, 0, 'substrates', styleBold) cidx = 1 for sr_ in range(len(Rlist[r_]['substrates'])): wsSt.write(ridx, cidx+sr_, Rlist[r_]['substrates'][sr_][0]) wsSt.write(ridx, cidx+sr_+1, Rlist[r_]['substrates'][sr_][1]) cidx += 1 ridx += 1 wsSt.write(ridx, 0, 'products', styleBold) cidx = 1 for pr_ in range(len(Rlist[r_]['products'])): wsSt.write(ridx, cidx+pr_, Rlist[r_]['products'][pr_][0]) wsSt.write(ridx, cidx+pr_+1, Rlist[r_]['products'][pr_][1]) cidx += 1 ridx += 1 gids = fba.getGroupIds() for g_ in range(len(gids)): mbrs = fba.getGroup(gids[g_]).getMemberIDs() wsGrp.write(0, g_, '{} ({})'.format(gids[g_], fba.getGroup(gids[g_]).getName()), styleBold) for m_ in range(len(mbrs)): wsGrp.write(m_+1, g_, str(mbrs[m_])) wb.save('{}.xls'.format(filename))
[docs]def writeModelToCOMBINEarchive(mod, fname=None, directory=None, sbmlname=None, withExcel=True, vc_given='CBMPy', vc_family='Software', vc_email='None', vc_org='cbmpy.sourceforge.net', add_cbmpy_annot=True, add_cobra_annot=True): """ Write a model in SBML and Excel format to a COMBINE archive using the following information: - *mod* a model object - *fname* the output base filename, archive will be <fname>.zip - *directory* [default=None] created the combine archive 'directory' - *sbmlname* [default='None'] If *sbmlname* is defined then SBML file is <sbmlname>.xml otherwise sbml will be <fname>.xml. - *withExcel* [default=True] include a human readable Excel spreadsheet version of the model - *vc_given* [default='CBMPy'] first name - *vc_family* [default='Software'] family name - *vc_email* [default='None'] email - *vc_org* [default='None'] organisation - *add_cbmpy_annot* [default=True] add CBMPy KeyValueData annotation. Replaces <notes> - *add_cobra_annot* [default=True] add COBRA <notes> annotation """ scTime = time.strftime('%Y-%m-%dT%H:%M:%S') + '%i:00' % (time.timezone/60/60) #self.writeSedXML(sedx=True) #sedxname = '%s.sed.omex' % (self.id) #sf = os.path.join(self.sedpath, sedxname) #self.__sedarchive__ = sf if directory != None: zfpath = os.path.join(directory, fname+'.sbex.zip') else: zfpath = fname+'.sbex.zip' zf = zipfile.ZipFile(zfpath, mode='w', compression=zipfile.ZIP_DEFLATED) if sbmlname != None: if sbmlname.endswith('.xml'): sbmlf = sbmlname xlf = sbmlname[:-4] else: sbmlf = sbmlname+'.xml' xlf = sbmlname else: sbmlf = fname+'.xml' xlf = fname ptmp = os.path.join(os.getcwd(), 'sedxtmp') if not os.path.exists(ptmp): os.makedirs(ptmp) assert os.path.exists(ptmp), "Could not create temporary archive directory: {}".format(ptmp) MFstr = '' MDstr = '' MFstr += '<omexManifest xmlns="http://identifiers.org/combine.specifications/omex-manifest">\n' MFstr += ' <content location="." format="http://identifiers.org/combine.specifications/omex"/>\n' #MFstr += ' <content location="./%s" format="http://identifiers.org/combine.specifications/sedml"/>\n' % os.path.split(self.__sedxml__)[-1] MFstr += ' <content location="./metadata.rdf" format="http://identifiers.org/combine.specifications/omex-metadata"/>\n' # SBML writeSBML3FBC(mod, sbmlf, ptmp, add_cbmpy_annot=add_cbmpy_annot, add_cobra_annot=add_cobra_annot) zf.write(os.path.join(ptmp, sbmlf), arcname=sbmlf) MFstr += ' <content location="./{}" format="http://identifiers.org/combine.specifications/sbml.level-3.version-1"/>\n'.format(sbmlf) # Excel if withExcel and _HAVE_XLWT_: writeModelToExcel97(mod, os.path.join(ptmp, xlf)) xlf += '.xls' zf.write(os.path.join(ptmp, xlf), arcname=xlf) MFstr += ' <content location="./{}" format="http://mediatypes.appspot.com/application/vnd.ms-excel"/>'.format(xlf) MF = file(os.path.join(ptmp, 'manifest.xml'), 'w') MF.write('<?xml version="1.0" encoding="utf-8"?>\n{}\n</omexManifest>\n'.format(MFstr)) MF.close() MD = file(os.path.join(ptmp, 'metadata.rdf'), 'w') MD.write('<?xml version="1.0" encoding="UTF-8"?>\n') MD.write('<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"\n') MD.write(' xmlns:dcterms="http://purl.org/dc/terms/"\n') MD.write(' xmlns:vCard="http://www.w3.org/2006/vcard/ns#"\n') MD.write(' xmlns:bqmodel="http://biomodels.net/models-qualifiers">\n') MD.write(' <rdf:Description rdf:about=".">\n') #MDstr += ' <dcterms:description>\n %s\n </dcterms:description>\n' % self.omex_description MDstr += ' <dcterms:creator>\n' MDstr += ' <rdf:Bag>\n' MDstr += ' <rdf:li rdf:parseType="Resource">\n' MDstr += ' <vCard:hasName rdf:parseType="Resource">\n' MDstr += ' <vCard:family-name>{}</vCard:family-name>\n'.format(vc_family) MDstr += ' <vCard:given-name>{}</vCard:given-name>\n'.format(vc_given) MDstr += ' </vCard:hasName>\n' MDstr += ' <vCard:hasEmail rdf:resource="{}" />\n'.format(vc_email) MDstr += ' <vCard:organization-name>\n' MDstr += ' {}\n'.format(vc_org) MDstr += ' </vCard:organization-name>\n' MDstr += ' </rdf:li>\n' MDstr += ' </rdf:Bag>\n' MDstr += ' </dcterms:creator>\n' MDstr += ' <dcterms:created rdf:parseType="Resource">\n' MDstr += ' <dcterms:W3CDTF>{}</dcterms:W3CDTF>\n'.format(scTime) MDstr += ' </dcterms:created>\n' MDstr += ' <dcterms:modified rdf:parseType="Resource">\n' MDstr += ' <dcterms:W3CDTF>{}</dcterms:W3CDTF>\n'.format(scTime) MDstr += ' </dcterms:modified>\n' MD.write('{}'.format(MDstr)) MD.write(' </rdf:Description>\n') MD.write('</rdf:RDF> \n') MD.close() zf.write(os.path.join(ptmp, 'manifest.xml'), arcname='manifest.xml') zf.write(os.path.join(ptmp, 'metadata.rdf'), arcname='metadata.rdf') zf.close() for f_ in os.listdir(ptmp): os.remove(os.path.join(ptmp, f_)) os.removedirs(ptmp) print('COMBINE archive created: {}'.format(fname+'.zip'))