diff --git a/build/lib/cometspy/__init__.py b/build/lib/cometspy/__init__.py new file mode 100644 index 0000000..f455f37 --- /dev/null +++ b/build/lib/cometspy/__init__.py @@ -0,0 +1,48 @@ +""" +*a python package for making and running COMETS simulations.* + +cometspy is a straight-forward python interface to the COMETS program. +Use for dynamic flux-balance analysis, with multiple models, heterogenous +spatio-temporal environments, and evolution. + +**To use this package, you must also have the actual COMETS program installed.** + +* For more information on COMETS, see https://runcomets.org + +* For COMETS development, see https://github.com/segrelab/comets + +* For development of this package, see http://github.com/segrelab/cometspy + +**cometspy workflow:** + +models -> layout, layout + params -> comets, comets.run() + +**cometspy hello.world** + +>>> import cobra.test +>>> import cometspy as c +>>> # make a model from a cobra model, open exchange reactions, and give a pop +>>> tb = cobra.test.create_test_model("textbook") +>>> m = c.model(tb) +>>> m.initial_pop = [0, 0, 1.e-4] +>>> m.open_exchanges() +>>> # make a layout with the model, and give it three nutrients +>>> l = c.layout([m]) +>>> l.set_specific_metabolite("glc__D_e", 0.01) +>>> l.set_specific_metabolite("nh4_e", 1000, static = True) +>>> l.set_specific_metabolite("pi_e", 1000, static = True) +>>> # make params and change one default +>>> p = c.params() +>>> p.set_param("maxCycles", 100) +>>> # make a simulation and run it! +>>> sim = c.comets(l, p) +>>> sim.run() +>>> print(sim.total_biomass) +>>> # saved data objects are pandas.DataFrames, and therefore can be plotted +>>> # sim.total_biomass.plot(x = "cycle") + +""" +from cometspy.comets import comets +from cometspy.model import model +from cometspy.layout import layout +from cometspy.params import params diff --git a/build/lib/cometspy/comets.py b/build/lib/cometspy/comets.py new file mode 100644 index 0000000..e3f4ea4 --- /dev/null +++ b/build/lib/cometspy/comets.py @@ -0,0 +1,771 @@ + +''' +The comets module runs COMETS simulations and stores output. + +Generally, a comets object is created just before calling run(). Afterwards, +any saved data (e.g. total_biomass) can be accessed from the object. + + + +''' + +import re +import subprocess as sp +import pandas as pd +import os +import glob +import numpy as np +import platform + +__author__ = "Djordje Bajic, Jean Vila, Jeremy Chacon, Ilija Dukovski" +__copyright__ = "Copyright 2024, The COMETS Consortium" +__credits__ = ["Djordje Bajic", "Jean Vila", "Jeremy Chacon", "Ilija Dukovski"] +__license__ = "MIT" +__version__ = "0.5.0" +__comets_compatibility__ = "2.12.0" # version of comets this was tested with (except signaling) +__comets_compatibility_signaling__ = "2.12.0" # version signaling was tested with + +__maintainer__ = "Djordje Bajic" +__email__ = "djordje.bajic@yale.edu" +__status__ = "Beta" + + +def _readlines_file(filename): + f = open(filename, 'r') + f_lines = f.readlines() + f.close() + return f_lines + + +class comets: + """ + the main simulation object to run COMETS + + The comets class stores a previously-made layout and params, and + interacts with COMETS to run the simulation. It also stores any + generated output, which could include total_biomass, biomass, + media, or fluxes, as set in the params object. std_out from the COMETS + simulation is saved in the attribute run_output and can be useful to + examine to debug errors. + + When creating a comets object, the optional relative_dir path is useful + when one is to run multiple simulations simultaneously, otherwise + temporary files may overwrite each other. + + Parameters + ---------- + + layout : layout + a cometspy.layout containing models and media information + parameters : params + a cometspy.params containing specified simulation parameters + relative_dir : str, optional + a directory to place temporary simulation files. + + Attributes + ---------- + + layout : cometspy.layout + the layout containing cometspy.model[s] and media information + params : cometspy.params + the params containing simulation and default biological parameters + working_dir : str + the directory at which to save temporary sim files + GUROBI_HOME : str + the directory where GUROBI exists on the system + COMETS_HOME : str + the directory where COMETS exists on the system + VERSION : str + the version of comets, read from the files at COMETS_HOME + classpath_pieces : dict + classpath separated into library name (key) and location (value) + JAVA_CLASSPATH : str + a generated (overwritable) string containing the java classpath + run_output : str + generated object containing text from COMETS sim's std_out + run_errors : str + generated object containing text from COMETS sim's std_err + total_biomass : pandas.DataFrame + generated object containing total biomass from simulation + biomass : pandas.DataFrame + generated object containing spatially-explicit biomass from sim + specific_media : pandas.DataFrame + generated object containing information on selected media from sim + media : pandas.Dataframe + generated object containing spatially-explicit media from sim + fluxes_by_species : dict{model_id : pandas.DataFrame} + generated object containing each species' spatial-explicit fluxes + genotypes : pandas.DataFrame + generated object containing genotypes if an evolution sim was run + + Examples + -------- + + >>> # assume layout and params have already been created + >>> # and cometspy was imported as c + >>> sim = c.comets(layout, params) + >>> sim.run() # this could take from seconds to hours + >>> print(sim.run_output) # to see the std_out + >>> sim.total_biomass.plot(x = "cycle") + >>> # assume params.all_params["writeBiomassLog"] == True, and that + >>> # params.all_params["BiomassLogRate"] = 1000, and that + >>> # params.all_params["maxCycles"] = 5000 + >>> im = sim.get_biomass_image(layout.models[0].id, 4000) + >>> from matplotlib import pyplot as plt + >>> plt.imshow(im / np.max(im)) + + + """ + + def __init__(self, layout, + parameters, relative_dir : str =''): + + # define instance variables + self.working_dir = os.getcwd() + '/' + relative_dir + try: + self.GUROBI_HOME = os.environ['GUROBI_COMETS_HOME'] + os.environ['GUROBI_HOME'] = self.GUROBI_HOME + except: + try: + self.GUROBI_HOME = os.environ['GUROBI_HOME'] + except: + try: + self.GUROBI_HOME = os.environ['COMETS_GUROBI_HOME'] + except: + self.GUROBI_HOME = '' + print("could not find environmental variable GUROBI_COMETS_HOME or GUROBI_HOME or COMETS_GUROBI_HOME") + print("COMETS will not work with GUROBI until this is solved. ") + print("Here is a solution:") + print(" 1. import os and set os.environ['GUROBI_HOME'] then try to make a comets object again") + print(" e.g. import os") + print(" os.environ['GUROBI_HOME'] = 'C:\\\\gurobi902\\\\win64'") + self.COMETS_HOME = os.environ['COMETS_HOME'] + self.VERSION = os.path.splitext(os.listdir(os.environ['COMETS_HOME'] + + '/bin')[0])[0] + + # set default classpaths, which users may change + self.__build_default_classpath_pieces() + self.__build_and_set_classpath() + self.__test_classpath_pieces() + + # check to see if user has the libraries where expected + + self.layout = layout + self.parameters = parameters + + # dealing with output files + self.parameters.set_param("useLogNameTimeStamp", False) + self.parameters.set_param("TotalBiomassLogName", + "totalbiomasslog" + '_' + hex(id(self))) + self.parameters.set_param("BiomassLogName", + "biomasslog" + '_' + hex(id(self))) + self.parameters.set_param("FluxLogName", + "fluxlog" + '_' + hex(id(self))) + self.parameters.set_param("MediaLogName", + "medialog" + '_' + hex(id(self))) + + def __build_default_classpath_pieces(self): + """ + sets up what it thinks the classpath should be + """ + self.classpath_pieces = {} + self.classpath_pieces['gurobi'] = (self.GUROBI_HOME + + '/lib/gurobi.jar') + + self.classpath_pieces['junit'] = glob.glob(self.COMETS_HOME + + '/lib/junit' + '/**/*junit*', + recursive=True)[0] + self.classpath_pieces['hamcrest'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/*hamcrest*', + recursive=True)[0] + + self.classpath_pieces['jogl_all'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/jogl-all.jar', + recursive=True)[0] + self.classpath_pieces['gluegen_rt'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/gluegen-rt.jar', + recursive=True)[0] + self.classpath_pieces['gluegen'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/gluegen.jar', + recursive=True)[0] + self.classpath_pieces['gluegen_rt_natives'] = glob.glob( + self.COMETS_HOME + '/lib' + '/**/gluegen-rt-natives-linux-amd64.jar', + recursive=True)[0] + self.classpath_pieces['jogl_all_natives'] = glob.glob( + self.COMETS_HOME + '/lib' + '/**/jogl-all-natives-linux-amd64.jar', + recursive=True)[0] + + self.classpath_pieces['jmatio'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/jmatio.jar', + recursive=True)[0] + self.classpath_pieces['jmat'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/jmatio.jar', + recursive=True)[0] + + self.classpath_pieces['concurrent'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/concurrent.jar', + recursive=True)[0] + self.classpath_pieces['colt'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/colt.jar', + recursive=True)[0] + + __lang3 = glob.glob(self.COMETS_HOME + + '/lib' + '/**/commons-lang3*jar', + recursive=True) + self.classpath_pieces['lang3'] = [i for i in __lang3 + if 'test' not in i + and 'sources' not in i][0] + + __math3 = glob.glob(self.COMETS_HOME + + '/lib' + '/**/commons-math3*jar', + recursive=True) + self.classpath_pieces['math3'] = [i for i in __math3 + if 'test' not in i + and 'sources' not in i + and 'tools' not in i + and 'javadoc' not in i][0] + + self.classpath_pieces['jdistlib'] = glob.glob(self.COMETS_HOME + + '/lib' + '/**/*jdistlib*', + recursive=True)[0] + self.classpath_pieces['bin'] = (self.COMETS_HOME + + '/bin/' + self.VERSION + '.jar') + + + def __build_and_set_classpath(self): + ''' builds the JAVA_CLASSPATH from the pieces currently in + self.classpath_pieces ''' + paths = list(self.classpath_pieces.values()) + if platform.system() == 'Windows': + classpath = ';'.join(paths) + classpath = '\"' + classpath + '\"' + self.JAVA_LIB = '\"' + self.GUROBI_HOME + '/lib;' + self.GUROBI_HOME+ '/bin;'+ self.COMETS_HOME+ '/lib/jogl/jogamp-all-platforms/lib'+ '\"' + + else: + classpath = ':'.join(paths) + self.JAVA_CLASSPATH = classpath + + def __test_classpath_pieces(self): + ''' checks to see if there is a file at each location in classpath + pieces. If not, warns the user that comets will not work without the + libraries. Tells the user to either edit those pieces (if in linux) + or just set the classpath directly''' + if platform.system() == 'Windows': + return # Windows uses the script, so classpath doesn't matter as long as env variable set + broken_pieces = self.__get_broken_classpath_pieces() + if len(broken_pieces) == 0: + pass # yay! class files are where we hoped + else: + print('Warning: java class libraries cannot be found') + print('These are the expected locations for dependencies:') + + print('Dependency \t\t\t expected path') + print('__________ \t\t\t _____________') + for key, value in broken_pieces.items(): + print('{}\t\t\t{}'.format(key, value)) + print('\n You have two options to fix this problem:') + print('1. set each class path correctly by doing:') + print(' comets.set_classpath(libraryname, path)') + print(' e.g. comets.set_classpath(\'hamcrest\', ' + + '\'/home/chaco001/comets/junit/hamcrest-core-1.3.' + + 'jar\')\n') + print(' note that versions dont always have to ' + + 'exactly match, but you\'re on your own if they ' + + 'don\'t\n') + print('2. fully define the classpath yourself by ' + + 'overwriting comets.JAVA_CLASSPATH') + print(' look at the current comets.JAVA_CLASSPATH ' + + 'to see how this should look.') + + def __get_broken_classpath_pieces(self): + ''' checks to see if there is a file at each location in classpath + pieces. Saves the pieces where there is no file and returns them as a + dictionary, where the key is the common name of the class library and + the value is the path ''' # + + broken_pieces = {} # + for key, value in self.classpath_pieces.items(): + if not os.path.isfile(value): # + broken_pieces[key] = value + return(broken_pieces) + + def set_classpath(self, libraryname : str, path : str): + """ + sets the full path to a required specified java library + + This can be used to set non-default classpaths. Note that currently, + it does not work for windows, because windows runs COMETS slightly + differently than Unix. + + Parameters + ---------- + + libraryname : str + the library for which the new path is being supplied. + path : str + the full path, including file name, to the library + + Examples + -------- + + >>> sim = c.comets(layout, params) + >>> sim.set_classpath("jmatio", "/opt/jmatio/jmatio.jar") + + """ + self.classpath_pieces[libraryname] = path + self.__build_and_set_classpath() + + def run(self, delete_files : bool = True): + """ + run a COMETS simulation + + This runs the COMETS simulation specified by the layout and params + objects supplied to this comets object. It creates the files needed + for COMETS into the current directory or the optional one specified + when this object was created. + + Once complete (or if an error has occurred), this tries to read + simulation logs including data, as well as the std_out in the + run_output attribute. + + If the optional delete_files is set to False, then temporary files and + data log files are not deleted. They are deleted by default. + + Parameters + ---------- + + delete_files : bool, optional + Whether to delete simulation and log files. The default is True. + + Examples + -------- + + >>> sim = c.comets(layout, params) + >>> sim.run(delete_files = True) + >>> print(sim.run_output) + >>> print(sim.total_biomass) + + """ + print('\nRunning COMETS simulation ...') + print('\nDebug Here ...') + + # If evolution is true, write the biomass but not the total biomass log + if self.parameters.all_params['evolution']: + self.parameters.all_params['writeTotalBiomassLog'] = False + self.parameters.all_params['writeBiomassLog'] = True + + to_append = '_' + hex(id(self)) + # write the files for comets in working_dir + c_global = self.working_dir + '.current_global' + to_append + c_package = self.working_dir + '.current_package' + to_append + c_script = self.working_dir + '.current_script' + to_append + + self.layout.write_necessary_files(self.working_dir, to_append) + + # self.layout.write_layout(self.working_dir + '.current_layout') + self.parameters.write_params(c_global, c_package) + + if os.path.isfile(c_script): + os.remove(c_script) + with open(c_script, 'a') as f: + f.write('load_comets_parameters ' + '.current_global' + to_append + '\n') + f.writelines('load_package_parameters ' + '.current_package' + to_append + '\n') + f.writelines('load_layout ' + '.current_layout' + to_append) + + if platform.system() == 'Windows': + self.cmd = ('\"' + self.COMETS_HOME + + '\\comets_scr' + '\" \"' + + c_script + + '\"') + else: + # simulate + self.cmd = ('java -classpath ' + self.JAVA_CLASSPATH + + # ' -Djava.library.path=' + self.D_JAVA_LIB_PATH + + ' edu.bu.segrelab.comets.Comets -loader' + + ' edu.bu.segrelab.comets.fba.FBACometsLoader' + + ' -script "' + c_script + '"') + + p = sp.Popen(self.cmd, + cwd = self.working_dir, + shell=True, stdout=sp.PIPE, stderr=sp.STDOUT) + + self.run_output, self.run_errors = p.communicate() + self.run_output = self.run_output.decode('ascii','ignore') + + if self.run_errors is not None: + self.run_errors = self.run_errors.decode('ascii','ignore') + else: + self.run_errors = "STDERR empty." + + # Raise RuntimeError if simulation had nonzero exit + self.__analyze_run_output() + + # '''----------- READ OUTPUT ---------------------------------------''' + # Read total biomass output + if self.parameters.all_params['writeTotalBiomassLog']: + tbmf = _readlines_file( + self.working_dir + self.parameters.all_params['TotalBiomassLogName']) + tbmf = [x.replace(",",".") for x in tbmf] # for systems that use commas as decimal place + self.total_biomass = pd.DataFrame([re.split(r'\t+', x.strip()) + for x in tbmf], + columns=['cycle'] + + self.layout.get_model_ids()) + self.total_biomass = self.total_biomass.astype('float') + self.total_biomass.cycle = self.total_biomass.cycle.astype('int') + if delete_files: + os.remove(self.working_dir + self.parameters.all_params['TotalBiomassLogName']) + + # Read flux + if self.parameters.all_params['writeFluxLog']: + + max_rows = 4 + max([len(m.reactions) for m in self.layout.models]) + + self.fluxes = pd.read_csv(self.working_dir + self.parameters.all_params['FluxLogName'], + delim_whitespace=True, + header=None, names=range(max_rows)) + # deal with commas-as-decimals + if any([isinstance(self.fluxes.iloc[0,i], str) for i in range(self.fluxes.shape[1])]): + self.fluxes = pd.read_csv(self.working_dir + self.parameters.all_params['FluxLogName'], + decimal = ",", delim_whitespace=True, + header=None, names=range(max_rows)) + if delete_files: + os.remove(self.working_dir + self.parameters.all_params['FluxLogName']) + self.__build_readable_flux_object() + + # Read media logs + if self.parameters.all_params['writeMediaLog']: + self.media = pd.read_csv(self.working_dir + self.parameters.all_params[ + 'MediaLogName'], delim_whitespace=True, names=('metabolite', + 'cycle', 'x', + 'y', + 'conc_mmol')) + # deal with commas-as-decimals + if isinstance(self.media.loc[0, "conc_mmol"], str): + self.media = pd.read_csv(self.working_dir + self.parameters.all_params[ + 'MediaLogName'], + decimal = ",", delim_whitespace=True, names=('metabolite', + 'cycle', 'x', + 'y', + 'conc_mmol')) + if delete_files: + os.remove(self.working_dir + self.parameters.all_params['MediaLogName']) + + # Read spatial biomass log + if self.parameters.all_params['writeBiomassLog']: + self.biomass = pd.read_csv(self.working_dir + self.parameters.all_params[ + 'BiomassLogName'], header=None, delimiter=r'\s+', names=['cycle', 'x', 'y', + 'species', 'biomass']) + # deal with commas-as-decimals + if isinstance(self.biomass.loc[0,"biomass"], str): + self.biomass = pd.read_csv(self.working_dir + self.parameters.all_params[ + 'BiomassLogName'], header=None, + decimal = ",", + delimiter=r'\s+', names=['cycle', 'x', 'y','species', 'biomass']) + + # cut off extension added by toolbox + self.biomass['species'] = [sp[:-4] if '.cmd' in sp else sp for sp in self.biomass.species] + + if delete_files: + os.remove(self.working_dir + self.parameters.all_params['BiomassLogName']) + + # Read evolution-related logs + if 'evolution' in list(self.parameters.all_params.keys()): + if self.parameters.all_params['evolution']: + genotypes_out_file = 'GENOTYPES_' + self.parameters.all_params[ + 'BiomassLogName'] + self.genotypes = pd.read_csv(genotypes_out_file, + header=None, delimiter=r'\s+', + names=['Ancestor', + 'Mutation', + 'Species']) + if delete_files: + os.remove(self.working_dir + genotypes_out_file) + + # Read specific media output + if self.parameters.all_params['writeSpecificMediaLog']: + spec_med_file = self.working_dir + self.parameters.all_params['SpecificMediaLogName'] + self.specific_media = pd.read_csv(spec_med_file, delimiter=r'\s+') + # deal with commas-as-decimals + if any([isinstance(self.specific_media.iloc[0,i], str) for i in range(3, self.specific_media.shape[1])]): + self.specific_media = pd.read_csv(spec_med_file, decimal = ",",delimiter=r'\s+') + + if delete_files: + os.remove(self.working_dir + self.parameters.all_params['SpecificMediaLogName']) + + # clean workspace + if delete_files: + self.layout.delete_model_files(self.working_dir) + os.remove(c_global) + os.remove(c_package) + os.remove(c_script) + os.remove(self.working_dir + '.current_layout' + to_append) + os.remove(self.working_dir + 'COMETS_manifest.txt') # todo: stop writing this in java + print('Done!') + + def __build_readable_flux_object(self): + """ comets.fluxes is an odd beast, where the column position has a + different meaning depending on what model the row is about. Therefore, + this function creates separate dataframes, stored in a dictionary with + model_id as a key, that are much more human-readable.""" + + self.fluxes_by_species = {} + for i in range(len(self.layout.models)): + model_num = i + 1 + + model_id = self.layout.models[model_num - 1].id + model_rxn_names = list(self.layout.models[ + model_num - 1].reactions.REACTION_NAMES) + model_rxn_len = len(model_rxn_names) + + sub_df = self.fluxes.loc[self.fluxes[3] == model_num] + + # this tosses extraneous columns and the model num column + sub_df = sub_df.drop(sub_df.columns[model_rxn_len+4: len(sub_df.columns)], + axis=1) + sub_df = sub_df.drop(sub_df.columns[3], axis=1) + sub_df.columns = ["cycle", "x", "y"] + model_rxn_names + self.fluxes_by_species[model_id] = sub_df + + def __analyze_run_output(self): + if "End of simulation" in self.run_output: + return + else: + print("Error: COMETS simulation did not complete\n") + print(" examine comets.run_output for the full java trace\n") + print(" if we detect a common reason, it will be stated in the RuntimeError at the bottom") + + if "Could not find or load main class edu.bu.segrelab.comets.Comets" in self.run_output: + message = "Could not find or load main class edu.bu.segrelab.comets.Comets\n" + message += "check if comets.version and comets.classpath_pieces['bin'] \n" + message += "point to an actual comets.jar file\n" + message += "this problem may be associated with a malformed\n" + message += "os.environ['COMETS_HOME'] environmental variable\n" + message += "that can be overwritten by, for example, \n" + message += ">>> import os\n" + message += ">>> os.environ['COMETS_HOME'] = '/home/comets/'" + raise RuntimeError(f"COMETS simulation did not complete:\n {message}") + + loc = self.run_output.find("NoClassDefFoundError") + if loc != -1: + error_string = self.run_output[(loc+22):(loc+100)] + missing_class = error_string.split("\n")[0] + if missing_class[0:6] == "gurobi": + message = "JAVA could not find gurobi.\n" + message += "try the following: \n" + message += ">>> import os\n" + message += ">>> os.environ['GUROBI_COMETS_HOME']\n" + message += "if there is nothing there try setting that variable\n" + message += "to the location of gurobi.jar, for example:\n" + message += ">>> os.environ['GUROBI_COMETS_HOME'] = '/opt/gurobi900/linux64'" + else: + message = f"JAVA could not find a needed class: {missing_class}\n" + message += "make sure it is in your java classpath\n" + message += "this can be changed with comets.set_classpath()\n" + message += "if in Unix. In Windows, it suggests that something changed\n" + message += "with the dependencies installed alongside COMETS" + raise RuntimeError(f"COMETS simulation did not complete:\n {message}") + + message = "undetected reason. examine comets.run_output for JAVA trace" + raise RuntimeError(f"COMETS simulation did not complete:\n {message}") + + + def get_metabolite_image(self, met : str, cycle : int) -> np.array: + """ + returns an image of metabolite concentrations at a given cycle + + This will only work if media was saved at the current cycle. This + requires the following parameters to have been set: + + params.set_param("writeMediaLog",True) + params.set_param("MediaLogRate", n) # n is in cycles + + Notes + ----- + + There may be a bug where cycles are + 1 what they should be. We will + fix this soon but for now be aware you may need to +1 your desired + cycle. The bug does not affect anything within the simulation. + + Parameters + ---------- + + met : str + the name of the metabolite + cycle : int + the cycle to get the metabolite data + + Returns + ------- + + A 2d numpy array which can be visualized like an image. + + Examples + -------- + + >>> sim.run() + >>> # assume sim.params.all_params["MediaLogRate"] = 100 and that + >>> # sim.params.all_params["maxCycles"] = 1000 + >>> im = sim.get_metabolite_image("ac_e", 500) + >>> from matplotlib import pyplot as plt # may need to be installed + >>> plt.imshow(im) + + """ + if not self.parameters.all_params['writeMediaLog']: + raise ValueError("media log was not recorded during simulation") + if met not in list(self.layout.media.metabolite): + raise NameError("met " + met + " is not in layout.media.metabolite") + if cycle not in list(np.unique(self.media['cycle'])): + raise ValueError('media was not saved at the desired cycle. try another.') + im = np.zeros((self.layout.grid[0], self.layout.grid[1])) + aux = self.media.loc[np.logical_and(self.media['cycle'] == cycle, + self.media['metabolite'] == met)] + for index, row in aux.iterrows(): + im[int(row['x']-1), int(row['y']-1)] = row['conc_mmol'] + return(im) + + def get_biomass_image(self, model_id : str, cycle : int) -> np.array: + """ + returns an image of biomass concentrations at a given cycle + + This will only work if biomass was saved at the current cycle. This + requires the following parameters to have been set: + >>> params.set_param("writeBiomassLog",True) + >>> params.set_param("BiomassLogRate", n) # n is in cycles + + + Parameters + ---------- + + model_id : str + the id of the model to get biomass data on + cycle : int + the cycle to get the biomass data + + Returns + ------- + + A 2d numpy array which can be visualized like an image. + + Examples + -------- + + >>> sim.run() + >>> # assume sim.params.all_params["BiomassLogRate"] = 100 and that + >>> # sim.params.all_params["maxCycles"] = 1000 + >>> im = sim.get_biomass_image("iJO1366", 500) # e.g. the ecoli model + >>> from matplotlib import pyplot as plt # may need to be installed + >>> plt.imshow(im) + + """ + if not self.parameters.all_params['writeBiomassLog']: + raise ValueError("biomass log was not recorded during simulation") + if model_id not in list(np.unique(self.biomass['species'])): + raise NameError("model " + model.id + " is not one of the model ids") + if cycle not in list(np.unique(self.biomass['cycle'])): + raise ValueError('biomass was not saved at the desired cycle. try another.') + im = np.zeros((self.layout.grid[0], self.layout.grid[1])) + aux = self.biomass.loc[np.logical_and(self.biomass['cycle'] == cycle, + self.biomass['species'] == model_id), :] + for index, row in aux.iterrows(): + im[int(row['x']-1), int(row['y']-1)] = row['biomass'] + return(im) + + def get_flux_image(self, model_id : str, + reaction_id : str, cycle : int) -> np.array: + """ + returns a 2d numpy array showing fluxes at a given cycle + + This will only work if flux was saved at the current cycle. This + requires the following parameters to have been set: + + params.set_param("writeFluxLog",True) + params.set_param("FluxLogRate", n) # n is in cycles + + Parameters + ---------- + + model_id : str + the id of the model about which to get fluxes + reaction_id : str + the id of the reaction about which to get fluxes + cycle : int + the cycle at which to get fluxes + + Returns + ------- + a 2d numpy array which can be visualized like an image + + Examples + -------- + + >>> sim.run() + >>> # assume sim.params.all_params["FluxLogRate"] = 100 and that + >>> # sim.params.all_params["maxCycles"] = 1000 + >>> # assume a model used was iJO1366 + >>> im = sim.get_flux_image("iJO1366", "EX_ac_e", 500) + >>> from matplotlib import pyplot as plt # may need to be installed + >>> plt.imshow(im) + + """ + if not self.parameters.all_params['writeFluxLog']: + raise ValueError("flux log was not recorded during simulation") + if model_id not in [m.id for m in self.layout.models]: + raise NameError("model " + model_id + " is not one of the model ids") + im = np.zeros((self.layout.grid[0], self.layout.grid[1])) + temp_fluxes = self.fluxes_by_species[model_id] + if cycle not in list(np.unique(temp_fluxes['cycle'])): + raise ValueError('flux was not saved at the desired cycle. try another.') + if reaction_id not in list(temp_fluxes.columns): + raise NameError("reaction_id " + reaction_id + + " is not a reaction in the desired model") + aux = temp_fluxes.loc[temp_fluxes['cycle'] == cycle, :] + for index, row in aux.iterrows(): + im[int(row['x']-1), int(row['y']-1)] = row[reaction_id] + return(im) + + def get_metabolite_time_series(self, upper_threshold : float = 1000.) -> pd.DataFrame: + """ + returns a pandas DataFrame containing extracellular metabolite time series + + Parameters + ---------- + + upper_threshold : float (optional) + metabolites ever above this are not returned + """ + total_media = self.media.groupby(by = ["metabolite", "cycle"]).agg(func = sum).reset_index().drop(columns = ["x", "y"]) + total_media = total_media.pivot(columns = "metabolite", values = "conc_mmol", index = ["cycle"]).reset_index().fillna(0.) + exceeded_threshold = [x for x in total_media.min().index[total_media.min() > upper_threshold] if x != "cycle"] + total_media = total_media.drop(columns = exceeded_threshold) + return(total_media) + + def get_species_exchange_fluxes(self, model_id : str, threshold : float = 0.) -> pd.DataFrame: + """ + returns a pandas DataFrame containing one model's exchange flux time series + + Parameters + ---------- + + model_id : str + id of the model + threshold : float (optional) + abs(flux) must exceed this to be returned + """ + fluxes = self.fluxes_by_species[model_id].copy() + fluxes = fluxes.groupby(by = "cycle").agg(func = sum).reset_index().drop(columns = ["x", "y"]) + not_exch = [x for x in fluxes.columns if "EX_" not in x and x != "cycle"] + fluxes = fluxes.drop(columns =not_exch) + didnt_exceed_threshold = [x for x in np.abs(fluxes).max().index[np.abs(fluxes).max() < threshold] if x != "cycle"] + fluxes = fluxes.drop(columns = didnt_exceed_threshold) + return(fluxes) + +# TODO: check for manual changes within layout that may not have triggered flags. See layout.write_layout for details +# TODO: fix read_comets_layout to always expect text addresses of comets model files +# TODO: remove comets manifest (preferably, dont write it) +# TODO: find quicker reading solution than the pd.read_csv stringIO hack +# TODO: solve weird rounding errors when reading from comets model +# TODO: give warning when unknown parameter is set +# TODO: write parameters in single file +# TODO: model biomass should be added in the layout "add_model" method, and not as a model class field +# TODO: adding models seems to remove media that has been previously set up +# TODO: write documentation properly for all functions and classes diff --git a/build/lib/cometspy/layout.py b/build/lib/cometspy/layout.py new file mode 100644 index 0000000..9cce0a0 --- /dev/null +++ b/build/lib/cometspy/layout.py @@ -0,0 +1,1500 @@ +import pandas as pd +import os +import numpy as np +import re +import math + +from cometspy.model import model + + +class layout: + ''' + object containing model and environmental definitions for a COMETS sim. + + The layout object is an essential object needed to start a comets + simulation, along with a params object. For a simulation to run, the layout + must contain model(s). Typically, the layout will be created after + creating models, and these models will be given in a list during layout + creation. Subsequently, methods like + set_specific_metabolite() will be used to create environmental media + definitions. Additionally, spatial information can be assigned, such as + the spatial dimensions in lattice units (layout.grid), spatial barriers, + etc. + + Parameters + ---------- + + input_obj : list(cometspy.model) or str, optional + either a list of cometspy models, or a path to a previously-written + COMETS layout to load. + + Examples + -------- + + >>> import cometspy as c + >>> import cobra.test + >>> e_cobra = cobra.test.create_test_model('textbook') + >>> e = c.model(e_cobra) + >>> e.open_exchanges() + >>> l = c.layout([e]) + >>> # use the media from the textbook in the COMETS media + >>> for key in e_cobra.medium.keys(): + >>> l.set_specific_metabolite(key[3:], 10.) + >>> l.display_current_media() + + Attributes + ---------- + models : list(cometspy.model) + a list of cometspy models. do not modify directly. use add_model() + grid : list of size two + the number of boxes in the x and y direction. default = [1, 1] + media : pandas.DataFrame + info on met initial values and diffusion. see set_specific_metabolite + refresh : list + info on met added per time. see set_specific_refresh + local_media : dict + info on spatial-explicit media. see set_specific_metabolite_at_location + local_refresh : dict + info on met added per time to loc. see set_specific_refresh_at_location + local_static : dict + info on constant met value at loc. see set_specific_static_at_location + default_diff_c : float + the baseline diffusion constant for all metabolites + barriers : list(tuple) + list of impenetrable locations. see add_barriers + region_map : numpy.array or None + integer array specifying regions. see set_region_map + region_parameters : dict + region-specific default diffusion params. see set_region_parameters + reactions : list + list of extracellular reactions. see add_external_reaction + periodic_media : list + list of media undergoing periodic cycles. see set_global_periodic_media + + + + + ''' + def __init__(self, input_obj : list = None): + + # define an empty layout that can be filled later + self.models = [] + self.grid = [1, 1] + self.media = pd.DataFrame(columns=['metabolite', + 'init_amount', + 'diff_c', + 'g_static', + 'g_static_val', + 'g_refresh']) + + # local_media is a dictionary with locations as keys, and as values, + # another dict with metabolite names as keys and amounts as values + # this information sets initial, location-specific media amounts. + self.local_media = {} + self.refresh = [] + self.local_refresh = {} + self.local_static = {} + self.initial_pop_type = "custom" # JMC not sure purpose of this + self.initial_pop = [] + self.all_exchanged_mets = [] + + self.default_diff_c = 5.0e-6 + self.default_g_static = 0 + self.default_g_static_val = 0 + self.default_g_refresh = 0 + + self.barriers = [] + + self.reactions = [] + self.periodic_media = [] + + self.region_map = None + self.region_parameters = {} + + self.models_pairs_frictions =[] + self.models_substrate_frictions =[] + + self.__local_media_flag = False + self.__diffusion_flag = False + self.__refresh_flag = False + self.__static_flag = False + self.__barrier_flag = False + self.__region_flag = False + self.__ext_rxns_flag = False + self.__periodic_media_flag = False + self.__models_frictions_flag = False + self.__models_pairs_frictions_flag = False + + if input_obj is None: + print('building empty layout model\nmodels will need to be added' + + ' with layout.add_model()') + elif isinstance(input_obj, str): + if not os.path.isfile(input_obj): + raise IOError(' when running comets.layout(), input_obj' + + ' is a string, and therefore should be a path' + + ' to a layout; however, no file could be found' + + ' at that path destination') + self.read_comets_layout(input_obj) + else: + if not isinstance(input_obj, list): + input_obj = [input_obj] # probably just one cobra model + self.models = input_obj + self.update_models() + + def set_region_parameters(self, region : int, + diffusion : list, + friction : float): + """ + set regions-specific diffusion and friction values + + COMETS can have different regions with different substrate + diffusivities and frictions. Here, you set those parameters. + + This does not affect a simulation unless a region map is also set, + using the layout.set_region_map() function. + + Parameters + ---------- + + region : int + the region (set by set_region_map) that uses the other parameters + diffusion : list(float) + a list containing diffusion constants (cm2/s) for each metabolite + friction : float + a single value for the friction observed in this region + + Examples + -------- + + >>> # this example assumes you have already created a cometspy model + >>> # called "e" and have imported cometspy as c and numpy as np + >>> # Here we are making a 2x2 world with two regions in which region + >>> # 1 metabolites have diffusion constant 1.e-6 cm2/s and in region + >>> # 2 metabolites have diffusion constant 1.e-7 cm2/s, and in both + >>> # cases mets have default friction constant 1. + >>> l = c.layout([e]) + >>> l.grid = [2,2] + >>> region_map = np.array([[1, 1], [2, 2]], dtype = int) + >>> l.set_region_map(region_map) + >>> n_mets = length(l.media) # how many metabolites there are + >>> l.set_region_parameters(1, [1.e-6] * n_mets, 1.) + >>> l.set_region_parameters(1, [1.e-7] * n_mets, 1.) + + + """ + if not self.__region_flag: + print("Warning: You are setting region parameters but a region" + + "map has not been set. Use layout.set_region_map() or these" + + "parameters will be unused") + self.region_parameters[region] = [diffusion, friction] + + def set_region_map(self, region_map : np.array): + """ + set a region_map dictating different regions in a spatial simulation + + COMETS can have different regions with different substrate diffusivities + and frictions. Here, you set the map defining the regions. Specifically, + you provide either: + 1) a numpy array whose shape == layout.grid, or + 2) a list of lists whose first length is grid[0] and second len is grid[1] + + Populating these objects should be integer values, beginning at 1 and + incrementing only, that define the different grid areas. These are + intimately connected to region_parameters, which are set with + layout.set_region_parameters() + + Parameters + ---------- + region_map : numpy.array(int) or list(list) + a 2d array of size layout.grid containing integers, or a list of + lists which can be coerced to an integer array using np.array() + + Examples + -------- + see set_region_parameters for an example + """ + if isinstance(region_map, list): + region_map = np.array(region_map) + if not tuple(self.grid) == region_map.shape: + raise ValueError("the shape of your region map must be the " + + "same as the grid size. specifically, \n" + + "tuple(layout.grid) == region_map.shape\n" + + "must be True after region_map = np.array(region_map)") + self.region_map = region_map + self.__region_flag = True + + def add_external_reaction(self, + rxnName : str, + metabolites : list, + stoichiometry : list, **kwargs): + """ + adds an extracellular reaction to the layout + + External reactions are reactions not tied to a growing (living) model. + They have rates defined in certain kwargs. There is no enyme + concentration; it is assumed fixed. + + Parameters + ---------- + rxnName : str + name of the external reaction + metabolites : list(str) + list of metabolite names (strs) + stoichiometry : list(float) + stoichiometry of the metabolites. must be same order as metabolites + kwargs : ['Kcat' or 'Km' or 'K'] + rate parameters. Either Kcat and Km, or K must be set. + + Examples + -------- + + >>> # set an external reaction that converts lactose into glucose + gal + >>> # assumes import cometspy as c and a model called m has been made + >>> rxn_name = 'lactase' + >>> metabolites = ['lcts_e', 'glc__D_e', 'gal_e'] + >>> stoichiometry = [-1., 1., 1.] + >>> K = 0.2 # mmol / hr + >>> l = c.layout([m]) # m is a cometspy.model that has exchanges for the mets above + >>> l.add_external_reaction(rxn_name, metabolites, stoichiometry, + K = K) + + """ + + ext_rxn = {'Name': rxnName, + 'metabolites': metabolites, + 'stoichiometry': stoichiometry} + + for key, value in kwargs.items(): + if key not in ['Kcat', 'Km', 'K']: + print('Warning: Parameter ' + key + ' i not recognized and ' + + 'will be ignored. Please set either Kcat and Km for' + + ' enzymatic reactions, or K for non catalyzed ones') + else: + ext_rxn[key] = value + + if 'Kcat' in ext_rxn and len([i for i in ext_rxn['stoichiometry'] + if i < 0]) > 1: + print('Warning: Enzymatic reactions are only allowed to have' + + 'one reactant') + + self.reactions.append(ext_rxn) + self.__ext_rxns_flag = True + + def set_global_periodic_media(self, + metabolite : str, function : str, + amplitude : float, period : float, + phase : float, offset : float): + """ + set cyclical changes in extracellular nutrient environment + + Parameters + ---------- + metabolite : str + name of the metabolite under a cycle + function : str, one of ['step', 'sin', 'cos', 'half_sin', 'half_cos'] + name of function that defines changing metabolite concentration + amplitude : float + amplitude of the function in mmol + period : float + duration of the function period in hr + phase : float + horizontal phase of the period (hr) + offset : float + vertical offset of the function (mmol) + + + """ + + if (metabolite not in self.media['metabolite'].values): + raise ValueError('the metabolite is not present in the media') + if (function not in ['step', 'sin', 'cos', 'half_sin', 'half_cos']): + raise ValueError(function + ': function unknown') + + self.periodic_media.append([self.media.index[self.media['metabolite'] + == metabolite][0], + function, amplitude, period, + phase, offset]) + self.__periodic_media_flag = True + + def read_comets_layout(self, input_obj : str): + """ + load a comets layout from a file + + If a COMETS layout file has been previously made, this function can + load that file into the current cometspy layout object. + + This is an unusual way of working with cometspy with rare uses. + + Parameters + ---------- + input_obj : str + the path to the existing COMETS layout file to be loaded + + """ + # .. load layout file + f_lines = [s for s in _read_file(input_obj).splitlines() if s] + filedata_string = os.linesep.join(f_lines) + end_blocks = [] + for i in range(0, len(f_lines)): + if '//' in f_lines[i]: + end_blocks.append(i) + + # '''----------- GRID ------------------------------------------''' + self.grid = [int(i) for i in f_lines[2].split()[1:]] + if len(self.grid) < 2: + print('\n Warning: Grid must contain only two values, but it \n' + + ' currently contains ' + len(self.grid) + + '\nCheck your layout file.') + + # '''----------- MODELS ----------------------------------------''' + ''' + Models can be specified in layout as either comets format models + or .xml format (sbml cobra compliant) + + ''' + # right now, assume all models in layouts are strings leading to + # comets model files + + # models need initial pop, so lets grab that first + + # '''----------- INITIAL POPULATION ----------------------------''' + lin_initpop = re.split('initial_pop', + filedata_string)[0].count('\n') + lin_initpop_end = next(x for x in end_blocks if x > lin_initpop) + + g_initpop = f_lines[lin_initpop].split()[1:] + + # TODO: I think we should deprecate these, it makes things difficult + # then, we could just generate these on-the-fly using the py toolbox, + # and have the initial_pop always appear to be 'custom' type to COMETS + # DB totally agree + + if (len(g_initpop) > 0 and g_initpop[0] in ['random', + 'random_rect', + 'filled', + 'filled_rect', + 'square']): + self.initial_pop_type = g_initpop[0] + self.initial_pop = [float(x) for x in g_initpop[1:]] + else: + self.initial_pop_type = 'custom' + + # .. local initial population values + lin_initpop += 1 + + # list of lists of lists. first level per-model, then per-location + temp_init_pop_for_models = [[] for x in + range(len(f_lines[0].split()[1:]))] + + for i in range(lin_initpop, lin_initpop_end): + ipop_spec = [float(x) for x in + f_lines[i].split()] + if len(ipop_spec)-2 != len(temp_init_pop_for_models): + print('\nWarning: Some initial population lines are corrupt.\n' + + 'Check your layout file') + if (ipop_spec[0] >= self.grid[0] or ipop_spec[1] >= self.grid[1]): + print('\nWarning: Some initial population values fall outside' + + '\nof the defined grid size. Check your layout file.') + else: + for j in range(len(ipop_spec)-2): + if ipop_spec[j+2] != 0.0: + if len(temp_init_pop_for_models[j]) == 0: + temp_init_pop_for_models[j] = [[ipop_spec[0], + ipop_spec[1], + ipop_spec[j+2]]] + else: + temp_init_pop_for_models[j].append([ipop_spec[0], + ipop_spec[1], + ipop_spec[j+2]]) + + models = f_lines[0].split()[1:] + if len(models) > 0: + for i, model_path in enumerate(models): + curr_model = model(model_path) + # TODO: get the initial pop information for each model, because the models own that info + curr_model.initial_pop = temp_init_pop_for_models[i] + self.add_model(curr_model) + self.update_models() + else: + print('Warning: No models in layout') + + # '''----------- MEDIA DESCRIPTION -----------------------------''' + lin_media = re.split('world_media', + filedata_string)[0].count('\n') + 1 + lin_media_end = next(x for x in end_blocks if x > lin_media) + + media_names = [] + media_conc = [] + for i in range(lin_media, lin_media_end): + metabolite = f_lines[i].split() + media_names.append(metabolite[0]) + media_conc.append(float(metabolite[1])) + + self.media['metabolite'] = media_names + self.media['init_amount'] = media_conc + + # '''----------- MEDIA DIFFUSION -------------------------------''' + self.__diffusion_flag = False + if 'DIFFUSION' in filedata_string: + self.__diffusion_flag = True + lin_diff = re.split('diffusion_constants', + filedata_string)[0].count('\n') + lin_diff_end = next(x for x in end_blocks if x > lin_diff) + + self.default_diff_c = float(re.findall(r'\S+', f_lines[lin_diff]. + strip())[1]) + + for i in range(lin_diff+1, lin_diff_end): + diff_spec = [float(x) for x in f_lines[i].split()] + if diff_spec[0] > len(self.media.metabolite)-1: + print('\n Warning: Corrupt line ' + str(i) + ' in ' + + 'diffusion values. \nLine not written. Check your ' + + 'layout file.') + else: + self.media.loc[int(diff_spec[0]), + 'diff_c'] = diff_spec[1] + + self.__local_media_flag = False + if 'MEDIA' in set(filedata_string.upper().strip().split()): + self.__local_media_flag = True + lin_media = [x for x in range(len(f_lines)) + if f_lines[x].strip().split()[0].upper() == + 'MEDIA'][0]+1 + lin_media_end = next(x for x in end_blocks if x > lin_media) + + for i in range(lin_media, lin_media_end): + media_spec = [float(x) for x in f_lines[i].split()] + if len(media_spec) != len(self.media.metabolite)+2: + print('\nWarning: Some local "media" lines are corrupt\n ' + + '(wrong number of entries). Check your layout file.') + elif (media_spec[0] >= self.grid[0] or + media_spec[1] >= self.grid[1]): + print('\nWarning: Some local "media" lines are corrupt\n' + + '(coordinates outside of the defined grid)\n' + + 'Check your layout file') + else: + loc = (int(media_spec[0]), int(media_spec[1])) + self.local_media[loc] = {} + media_spec = media_spec[2:] + for j in range(len(media_spec)): + if media_spec[j] != 0: + self.local_media[loc][ + self.all_exchanged_mets[j]] = media_spec[j] + + # '''----------- MEDIA REFRESH----------------------------------''' + # .. global refresh values + self.__refresh_flag = False + if 'REFRESH' in filedata_string.upper(): + self.__refresh_flag = True + lin_refr = re.split('REFRESH', + filedata_string.upper())[0].count('\n') + lin_refr_end = next(x for x in end_blocks if x > lin_refr) + + g_refresh = [float(x) for x in f_lines[lin_refr].split()[1:]] + + if len(g_refresh) != len(media_names): + print('\nWarning: Some local refresh lines are corrupt\n ' + + '(wrong number of entries). Check your layout file.') + + else: + self.media['g_refresh'] = g_refresh + + # .. local refresh values + lin_refr += 1 + + for i in range(lin_refr, lin_refr_end): + refr_spec = [float(x) for x in f_lines[i].split()] + if len(refr_spec) != len(self.media.metabolite)+2: + print('\nWarning: Some local "refresh" lines are corrupt\n ' + + '(wrong number of entries). Check your layout file.') + + elif (refr_spec[0] >= self.grid[0] or + refr_spec[1] >= self.grid[1]): + print('\nWarning: Some local "refresh" lines are corrupt\n' + + '(coordinates outside of the defined grid)\n' + + 'Check your layout file') + else: + loc = (int(refr_spec[0]), int(refr_spec[1])) + self.local_refresh[loc] = {} + refr_spec = refr_spec[2:] + for j in range(len(refr_spec)): + if refr_spec[j] != 0: + self.local_refresh[loc][ + self.all_exchanged_mets[j]] = refr_spec[j] + + # region-based information (substrate diffusivity,friction, layout) + self.__region_flag = False + + if 'SUBSTRATE_LAYOUT' in filedata_string.upper(): + lin_substrate = re.split('SUBSTRATE_LAYOUT', + filedata_string.upper())[0].count('\n') + lin_substrate_end = next(x for x in end_blocks if x > lin_substrate) + region_map_data = [] + for i in range(lin_substrate+1, lin_substrate_end): + region_map_data.append([int(x) for x in f_lines[i].split()]) + region_map_data = np.array(region_map_data, dtype=int) + if region_map_data.shape != tuple(self.grid): + print('\nWarning: Some substrate_layout lines are ' + + ' longer or shorter than the grid width, or there are more' + + ' lines than the grid length. Check your layout file.' + + 'Check your layout file.') + self.__region_flag = True + self.region_map = region_map_data + + if 'SUBSTRATE_DIFFUSIVITY' in filedata_string.upper(): + lin_substrate = re.split('SUBSTRATE_DIFFUSIVITY', + filedata_string.upper())[0].count('\n') + lin_substrate_end = next(x for x in end_blocks if x > lin_substrate) + self.region_parameters = {} + region = 1 + for i in range(lin_substrate+1, lin_substrate_end): + self.region_parameters[region] = [None, None] + self.region_parameters[region][0] = [float(x) + for x in f_lines[i].split()] + if len(self.region_parameters[ + region][0]) != len(self.media.metabolite): + print('\nWarning: Some substrate_diffusivity lines are ' + + ' longer or shorter than the number of metabolites\n' + + 'Check your layout file.') + region += 1 + + if 'SUBSTRATE_FRICTION' in filedata_string.upper(): + lin_substrate = re.split('SUBSTRATE_FRICTION', + filedata_string.upper())[0].count('\n') + lin_substrate_end = next(x for x in end_blocks if x > lin_substrate) + region = 1 + for i in range(lin_substrate+1, lin_substrate_end): + self.region_parameters[region][1] = float(f_lines[i].split()[0]) + region += 1 + + # '''----------- STATIC MEDIA ----------------------------------''' + # .. global static values + self.__static_flag = False + if 'STATIC' in filedata_string.upper(): + self.__static_flag = True + lin_static = re.split('STATIC', + filedata_string.upper())[0].count('\n') + lin_stat_end = next(x for x in end_blocks if x > lin_static) + + g_static = [float(x) for x in f_lines[lin_static].split()[1:]] + + if len(g_static) != 2*len(self.media.metabolite): + print('\nWarning: Wrong number of global static values .\n' + + 'Check your layout file') + + else: + self.media.loc[:, 'g_static'] = [int(x) + for x in g_static[0::2]] + self.media.loc[:, 'g_static_val'] = [float(x) for x in + g_static[1::2]] + + # .. local static values + lin_static += 1 + for i in range(lin_static, lin_stat_end): + stat_spec = [float(x) for x in f_lines[i].split()] + if len(stat_spec) != (2*len(self.media.metabolite))+2: + print('\nWarning: Wrong number of local static values at some\n' + + 'lines. Check your layout file.') + elif (stat_spec[0] >= self.grid[0] or + stat_spec[1] >= self.grid[1]): + print('\nWarning: Some local "static" lines have coordinates\n' + + 'that fall outside of the defined grid. Check your layout file') + else: + loc = (int(stat_spec[0]), int(stat_spec[1])) + self.local_static[loc] = {} + stat_spec = stat_spec[2:] + for j in range(int(len(stat_spec)/2)): + if stat_spec[j*2] != 0: + self.local_static[loc][ + self.all_exchanged_mets[j]] = stat_spec[j*2+1] + + def get_model_ids(self) -> list: + """ + returns a list of the ids of the models in this layout. + + """ + ids = [x.id for x in self.models] + return(ids) + + def write_necessary_files(self, working_dir : str, to_append = ""): + """ + writes the layout and the model files to file + + This method is called by comets.run(). Alternatively it may be called + directly if one wishes to inspect comets layout and model files. + + A path must be specified if called directly. + + Note: the layout file will be called ".current_layout" and therefore + be hidden from file browsers by default. + + Parameters + ---------- + working_dir : str, + The directory where the files will be written. + to_append : str, + String to append to written filenames + + """ + self.__check_if_initial_pops_in_range() + self.write_layout(working_dir, to_append) + self.write_model_files(working_dir) + + def write_model_files(self, working_dir=""): + '''writes each model file''' + for m in self.models: + m.write_comets_model(working_dir) + + def delete_model_files(self, working_dir): + """ + deletes model files in specified directory + """ + for m in self.models: + m.delete_comets_model(working_dir) + + def display_current_media(self): + """ + a utility function to show current non-zero media amounts + + """ + print(self.media[self.media['init_amount'] != 0.0]) + + def add_barriers(self, barriers : list): + """ + adds impenetrable barriers to the layout at specified locations + + Barriers can be used in COMETS to create impenetrable spatial + locations, for example to simulate rocks or to make a rounded + simulation. Neither biomass nor metabolites can diffuse through + barriers. + + In COMETS, locations are zero-ordered. + + Parameters + ---------- + + barriers : list(tuple) + A list containing tuples of integers of x,y locations. + + Examples + -------- + + >>> # make a diagonal line of barriers + >>> l = c.layout([model1, model2]) # assumes model1,2 are made + >>> l.grid = [5,5] + >>> l.add_barriers([(0,0),(1,1),(2,2),(3,3),(4,4)]) + + """ + # first see if they provided only one barrier not in a nested list, and if + # so, put it into a list + if len(barriers) == 2: + if isinstance(barriers[0], int): + barriers = [barriers] + # now check each barrier and make sure it has 2 ints that fit within the grid size + for b in barriers: + try: + if len(b) != 2 or b[0] >= self.grid[0] or b[1] >= self.grid[1]: + raise ValueError + self.barriers.append((int(b[0]), int(b[1]))) + except ValueError: + print('ERROR ADDING BARRIERS in add_barriers\n') + print("expecting barriers to be a list of tuples of coordinates which" + + " fit within the current grid") + print(" such as layout.grid = [5,5]") + print(" barriers = [(0,0),(1,1),(2,2),(4,4)]") + print(" layout.add_barriers(barriers)") + if len(self.barriers) > 0: + self.__barrier_flag = True + self.barriers = list(set(self.barriers)) + + def set_metabolite_diffusion(self, diffusion_constant : float): + """ + sets the default diffusion constant for all metabolites + + Specific diffusion constants can be set subsequent to setting this. + + Parameters + ---------- + diffusion_constant : float + the diffusion constant, in cm2 / s, that all metabolites use + + """ + try: + if not isinstance(diffusion_constant, float): + raise ValueError + except ValueError: + print("ERROR, diffusion_constant must be a float") + return + print("Warning: set_metabolite_diffusion overwrites all diffusion constants\nthis should be performed prior to setting specific metabolite diffusion") + self.default_diff_c = diffusion_constant + self.media.loc[:, "diff_c"] = diffusion_constant + self.__diffusion_flag = True + + def __check_if_diffusion_flag_should_be_set(self): + """ Internal function. If someone directly changes the diff_c in the + media dataframe, cometspy will be unaware and not set __diffusion_flag = True, + rendering that change moot. This function is run when doing write_layout() + to flip the flag if necessary """ + if len(np.unique(self.media.diff_c)) > 1: + self.__diffusion_flag = True + elif self.media.diff_c[0] != self.default_diff_c: + self.__diffusion_flag = True + + def set_specific_metabolite_diffusion(self, met : str, + diffusion_constant : float): + """ + sets the diffusion constant of a specific metabolite + + Parameters + ---------- + met : str + name of the metabolite + diffusion_constant : float + the diffusion constant, in cm2/s, of the named metabolite + + Examples + -------- + >>> l = c.layout([m]) + >>> l.set_specific_metabolite_diffusion("ac_e", 5.2e-6) + >>> l.set_specific_metabolite_diffusion("gal_e", 4.e-6) + + + """ + try: + if not isinstance(diffusion_constant, float): + raise ValueError + except ValueError: + print("ERROR, diffusion_constant must be a float") + return + try: + if len(self.media.loc[self.media.metabolite == met, "diff_c"]) != 1: + raise ValueError + except ValueError: + print("ERROR, met not in layout.media.metabolite list\ncannot change metabolite diffusion constant") + self.media.loc[self.media.metabolite == met, "diff_c"] = diffusion_constant + self.__diffusion_flag = True + + def set_specific_metabolite(self, met : str, amount : float, + static : bool = False): + """ + sets the initial (or constant) value for a metabolite across space + + This is the most common way to set a metabolite concentration. It + sets the initial value of the named metabolite to 'amount' in every + location. Optionally, if static = True, this concentration is fixed + during the simulation. + + Parameters + ---------- + met : str + name of the metabolite + amount : float + initial value for the metabolite in each box, in mmol + static : bool, optional + DEFAULT = False. If True, the 'amount' is fixed over time. + + Examples + -------- + + >>> l = c.layout([model]) + >>> l.set_specific_metabolite("glc__D_e", 0.005) + >>> l.set_specific_metabolite("o2_e", 15, static = True) + >>> # a common thing to do is populate this from the cobra model, as + >>> # shown here: + >>> import cobra.test + >>> import cometspy as c + >>> cobra_model = cobra.test.create_test_model("textbook") + >>> model = c.model(cobra_model) + >>> model.open_exchanges() + >>> l = c.layout([model]) + >>> for key, value in cobra_model.medium.items(): + >>> l.set_specific_metabolite(key[3:], value) + + + """ + if met in set(self.media['metabolite']): + self.media.loc[self.media['metabolite'] == met, + 'init_amount'] = amount + if static: + self.media.loc[self.media['metabolite'] == met, + 'g_static'] = 1 + self.media.loc[self.media['metabolite'] == met, + 'g_static_val'] = amount + else: + newrow = {'metabolite': met, + 'g_refresh': self.default_g_refresh, + 'g_static': 1 if static else self.default_g_static, + 'g_static_val': amount if static else self.default_g_static_val, + 'init_amount': amount, + 'diff_c': self.default_diff_c} + newrow = pd.DataFrame([newrow], columns=newrow.keys()) + self.media = pd.concat([self.media, + newrow], + axis=0, sort=False, ignore_index=True) + print('Warning: The added metabolite (' + met + ') is not ' + + 'able to be taken up by any of the current models') + + def set_specific_metabolite_at_location(self, met : str, + location : tuple, + amount : float): + """ + set an initial value of a metabolite at a specific spatial location + + In contrast to set_specific_metabolite, which sets initial values + universally over space, this method sets a specific metabolite + initial value at one location only. + + Parameters + ---------- + met : str + name of the metabolite + location : tuple + a tuple of integers specifying a specific location, e.g. (3,2) + amount : float + the initial value for the metabolite at the location, in mmol + + Examples + -------- + >>> l = c.layout([model]) + >>> l.grid = [10,1] + >>> l.set_specific_metabolite_at_location("glc__D_e", (9,0), 0.005) + + Notes + ----- + Remember that COMETS locations are zero-ordered. + + """ + if met not in self.all_exchanged_mets: + raise Exception('met is not in the list of exchangeable mets') + self.__local_media_flag = True + if location not in list(self.local_media.keys()): + self.local_media[location] = {} + self.local_media[location][met] = amount + + def set_specific_refresh(self, met : str, amount : float): + """ + set the amount to increment a metabolite continuously through time + + metabolites can be "refreshed" or "dripped in" during a COMETS + simulation. This sets the amount added each hour, which is divided + up by the number of time steps. It applies that refresh to all boxes + in space. + + + Parameters + ---------- + met : str + name of the metabolite refreshed + amount : float + the amount added to (or subtracted from) each box, in mmol / hr. + + + """ + if met in self.media['metabolite'].values: + self.media.loc[self.media['metabolite'] == met, + 'g_refresh'] = amount + self.__refresh_flag = True + else: + raise Exception("the specified metabolite " + met + + "is not in the medium; add it first!") + + def set_specific_refresh_at_location(self, met : str, + location : tuple, + amount : float): + """ + sets the amount to adjust a metabolite at a specific location in time + + Like set_specific_refresh, but for a specific lacation. This + Method is used to specify how much a metabolite is increased (or + decreased) each hour in a specific location. + + Parameters + ---------- + met : str + the name of the metabolite + location : tuple + the (x, y) location to adjust, with x and y of type int + amount : float + the amount added to (or subtracted from) the box, in mmol / hr + + + """ + if met not in self.all_exchanged_mets: + raise Exception('met is not in the list of exchangeable mets') + self.__refresh_flag = True + if location not in list(self.local_refresh.keys()): + self.local_refresh[location] = {} + self.local_refresh[location][met] = amount + + def set_specific_static(self, met : str, amount : float): + """ + sets a metabolite to a fixed value in every location + + Parameters + ---------- + met : str + name of the metabolite + amount : float + the amount, in mmol, that is in each box at each time step + + """ + if met in self.media['metabolite'].values: + self.media.loc[self.media['metabolite'] == met, + 'g_static'] = 1 + self.media.loc[self.media['metabolite'] == met, + 'g_static_val'] = amount + self.__static_flag = True + else: + print("the specified metabolite " + met + + "is not in the medium; add it first!") + + def set_specific_static_at_location(self, met : str, + location : tuple, + amount : float): + """ + sets a metabolite to a fixed value at a specific location + + Parameters + ---------- + met : str + name of the metabolite + location : tuple + the (x, y) location of the fixed metabolite amount. x, y are int + amount : float + the amount, in mmol, that is in each box at each time step + + + """ + if met not in self.all_exchanged_mets: + raise Exception('met is not in the list of exchangeable mets') + self.__static_flag = True + if location not in list(self.local_static.keys()): + self.local_static[location] = {} + self.local_static[location][met] = amount + + def set_models_substrate_frictions(self, + friction : list): + """ + sets the values of the friction paremeters between the organisms and substrate + + Parameters + ---------- + model : int + the order of the model + location : tuple + the (x, y) location of the fixed metabolite amount. x, y are int + amount : float + the amount, in mmol, that is in each box at each time step + + + """ + self.__models_frictions_flag = True + self.models_substrate_frictions=friction + + def set_models_pairs_frictions(self, + pairs_frictions : list): + """ + sets the values of the friction paremeters between the organisms and substrate + + Parameters + ---------- + model : int + the order of the model + location : tuple + the (x, y) location of the fixed metabolite amount. x, y are int + amount : float + the amount, in mmol, that is in each box at each time step + + + """ + self.__models_pairs_frictions_flag = True + self.models_pairs_frictions=pairs_frictions + + def add_typical_trace_metabolites(self, amount : float=1000.0, + static : bool =True): + """ + adds common BIGG-style metabolites to the environment. This only + works if your models use metabolites which are formatted like 'ca2_e' + and then should still only be used as a starting point. + + Parameters + ---------- + amount : float, optional + the amount, in mmol, of these metabolites. Default is 1000.0. + static : bool, optional + The default is True. If false, these metabolites are depletable. + + + """ + trace_metabolites = ['ca2_e', + 'cl_e', + 'cobalt2_e', + 'cu2_e', + 'fe2_e', + 'fe3_e', + 'h_e', + 'k_e', + 'h2o_e', + 'mg2_e', + 'mn2_e', + 'mobd_e', + 'na1_e', + 'ni2_e', + 'nh4_e', + 'o2_e', + 'pi_e', + 'so4_e', + 'zn2_e'] + + for met in trace_metabolites: + if met in set(self.media['metabolite']): + self.media.loc[self.media['metabolite'] == met, + 'init_amount'] = amount + if static: + self.media.loc[self.media['metabolite'] == met, + 'g_static'] = 1 + self.media.loc[self.media['metabolite'] == met, + 'g_static_val'] = amount + + else: + newrow = {'metabolite': met, + 'g_refresh': self.default_g_refresh, + 'g_static': 1 if static else self.default_g_static, + 'g_static_val': amount if static else self.default_g_static_val, + 'init_amount': amount, + 'diff_c': self.default_diff_c} + newrow = pd.DataFrame([newrow], columns=newrow.keys()) + self.media = pd.concat([self.media, + newrow], + axis=0, sort=False) + self.media = self.media.reset_index(drop=True) + + def write_layout(self, working_dir : str, to_append = ""): + """ + writes just the COMETS layout file to the supplied path + + Parameters + ---------- + working_dir : str + the path to the directory where .current_layout will be written + + + """ + # right now we only check if a user manually set a diff_c. Ideally, + # we should check for manual changes to everything. Alternatively, + # we should print all blocks no matter what. + self.__check_if_diffusion_flag_should_be_set() + outfile = working_dir + ".current_layout" + to_append + if os.path.isfile(outfile): + os.remove(outfile) + + lyt = open(outfile, 'a') + self.__write_models_and_world_grid_chunk(lyt) + self.__write_media_chunk(lyt) + self.__write_diffusion_chunk(lyt) + self.__write_local_media_chunk(lyt) + self.__write_refresh_chunk(lyt) + self.__write_static_chunk(lyt) + self.__write_barrier_chunk(lyt) + self.__write_regions_chunk(lyt) + self.__write_periodic_media_chunk(lyt) + self.__write_models_frictions_chunk(lyt) + self.__write_models_pairs_frictions_chunk(lyt) + lyt.write(r' //' + '\n') + + self.__write_initial_pop_chunk(lyt) + self.__write_ext_rxns_chunk(lyt) + lyt.close() + + def __write_models_and_world_grid_chunk(self, lyt): + """ writes the top 3 lines to the open lyt file""" + + model_file_line = "{}.cmd".format(".cmd ".join(self.get_model_ids())).split(" ") + model_file_line = "".join(["./" + _ + " " for _ in model_file_line]) + model_file_line = "model_file " + model_file_line + "\n" + lyt.write(model_file_line) + lyt.write(' model_world\n') + + lyt.write(' grid_size ' + + ' '.join([str(x) for x in self.grid]) + '\n') + + def __write_media_chunk(self, lyt): + """ used by write_layout to write the global media information to the + open lyt file """ + lyt.write(' world_media\n') + for i in range(0, len(self.media)): + lyt.write(' ' + self.media.metabolite[i] + + ' ' + str(self.media.init_amount[i]) + '\n') + lyt.write(r' //' + '\n') + + def __write_local_media_chunk(self, lyt): + """ used by write_layout to write the location-specific initial + metabolite data""" + if self.__local_media_flag: + lyt.write(' media\n') + locs = list(self.local_media.keys()) + for loc in locs: + # this chunk goes in order, not by name, so must get met number + # for each location, make a list with zeros for each met. Put + # non-zero numbers where the self.local_media tells us to + met_amounts_in_order = [0] * len(self.all_exchanged_mets) + for met in list(self.local_media[loc].keys()): + met_amounts_in_order[ + self.__get_met_number(met)] = self.local_media[loc][met] + lyt.write(' ') + lyt.write('{} {} '.format(loc[0], loc[1])) + lyt.write(' '.join(str(x) for x in met_amounts_in_order)) + lyt.write('\n') + lyt.write(' //\n') + + def __write_refresh_chunk(self, lyt): + if self.__refresh_flag: + lyt.write(' media_refresh ' + + ' '.join([str(x) for x in self.media. + g_refresh.tolist()]) + + '\n') + locs = list(self.local_refresh.keys()) + if len(locs) > 0: + for loc in locs: + met_amounts_in_order = [0] * len(self.all_exchanged_mets) + for met in list(self.local_refresh[loc].keys()): + met_amounts_in_order[self.__get_met_number(met)] = self.local_refresh[loc][met] + met_amounts_in_order.insert(0, loc[1]) + met_amounts_in_order.insert(0, loc[0]) + lyt.write(' ' + + ' '.join([str(x) for x in met_amounts_in_order]) + + '\n') + lyt.write(r' //' + '\n') + + def __write_static_chunk(self, lyt): + if self.__static_flag: + g_static_line = [None]*(len(self.media)*2) + g_static_line[::2] = self.media.g_static + g_static_line[1::2] = self.media.g_static_val + lyt.write(' static_media ' + + ' '.join([str(x) for x in g_static_line]) + '\n') + locs = list(self.local_static.keys()) + if len(locs) > 0: + for loc in locs: + # this is 2 * len because there is a pair of values for each met + # the first value is a flag--0 if not static, 1 if static + # the second value is the amount if it is static + met_amounts_in_order = [0] * 2 * len(self.all_exchanged_mets) + for met in list(self.local_static[loc].keys()): + met_amounts_in_order[self.__get_met_number(met) * 2] = 1 + met_amounts_in_order[ + self.__get_met_number(met) * 2 + 1] = self.local_static[ + loc][met] + met_amounts_in_order.insert(0, loc[1]) + met_amounts_in_order.insert(0, loc[0]) + lyt.write(' ' + + ' '.join([str(x) for x in + met_amounts_in_order]) + + '\n') + lyt.write(r' //' + '\n') + + def __write_diffusion_chunk(self, lyt): + """ used by write_layout to write the metab-specific + diffusion data to the open lyt file """ + + if self.__diffusion_flag: + lyt.write(' diffusion_constants ' + + str(self.default_diff_c) + + '\n') + for i in range(0, len(self.media)): + if not math.isnan(self.media.diff_c[i]): + lyt.write(' ' + str(i) + ' ' + + str(self.media.diff_c[i]) + '\n') + lyt.write(r' //' + '\n') + + def __write_barrier_chunk(self, lyt): + """ used by write_layout to write the barrier section to the open lyt file """ + if self.__barrier_flag: + lyt.write(' barrier\n') + for barrier in self.barriers: + lyt.write(' {} {}\n'.format(barrier[0], barrier[1])) + lyt.write(' //\n') + + def __write_ext_rxns_chunk(self, lyt): + """ used by write_layout to write the external reactions section + to the open lyt file + """ + reactants = [] + enzymes = [] + products = [] + + if self.__ext_rxns_flag: + for i, rxn in enumerate(self.reactions): + + current_reactants = [self.media.index[ + self.media['metabolite'] == + rxn['metabolites'][k]].tolist()[0]+1 + for k in range(len(rxn['metabolites'])) + if rxn['stoichiometry'][k] < 0] + + current_products = [self.media.index[ + self.media['metabolite'] == + rxn['metabolites'][k]].tolist()[0]+1 + for k in range(len(rxn['metabolites'])) + if rxn['stoichiometry'][k] > 0] + + current_react_stoich = [k for k in rxn['stoichiometry'] + if k < 0] + + current_prod_stoich = [k for k in rxn['stoichiometry'] + if k > 0] + + for ind, k in enumerate(current_reactants): + if ind == 0: + + cl = (' ' + str(i+1) # reaction + + ' ' + str(k) # metabolite + + ' ' + str(-current_react_stoich[ind]) # stoich + + ' ' + + str([rxn['K'] if 'K' in rxn else rxn['Km']][0]) + + '\n') + reactants.append(cl) + else: + cl = (' ' + str(i+1) + + ' ' + str(k) + + ' ' + str(-current_react_stoich[ind]) + + ' ' + '\n') + reactants.append(cl) + + for ind, k in enumerate(current_products): + cl = (' ' + str(i+1) + + ' ' + str(k) + + ' ' + str(current_prod_stoich[ind]) + + ' ' + '\n') + products.append(cl) + + if 'Kcat' in rxn: + cl = (' ' + str(i+1) + + ' ' + str(rxn['Kcat']) + + '\n') + enzymes.append(cl) + + # write the reaction lines + lyt.write('reactions\n') + lyt.write(' reactants\n') + for i in reactants: + lyt.write(i) + + lyt.write(' enzymes\n') + for i in enzymes: + lyt.write(i) + + lyt.write(' products\n') + for i in products: + lyt.write(i) + lyt.write('//\n') + + def __write_periodic_media_chunk(self, lyt): + """ used by write_layout to write the periodic media + """ + if self.__periodic_media_flag: + lyt.write(' periodic_media global\n') + for media in self.periodic_media: + lyt.write(' {} {} {} {} {} {}\n'. + format(media[0], media[1], media[2], + media[3], media[4], media[5])) + lyt.write(' //\n') + + def __write_regions_chunk(self, lyt): + """ used by write_layout to write the regions section to the open lyt file + specifically this section includes "substrate_diffusivity" "substrate_friction" + and "substrate_layout". + """ + if self.__region_flag: + keys = list(self.region_parameters.keys()) + keys.sort() + lyt.write(' substrate_diffusivity\n') + for key in keys: + diff = [str(x) for x in self.region_parameters[key][0]] + line = " " + " ".join(diff) + "\n" + lyt.write(line) + lyt.write(" //\n") + lyt.write(" substrate_friction\n") + for key in keys: + fric = self.region_parameters[key][1] + line = " " + str(fric) + "\n" + lyt.write(line) + lyt.write(" //\n") + lyt.write(" substrate_layout\n") + for i in range(self.region_map.shape[0]): + for j in range(self.region_map.shape[1]): + lyt.write(" ") + lyt.write(str(self.region_map[i, j])) + lyt.write("\n") + lyt.write(" //\n") + + def __write_initial_pop_chunk(self, lyt): + """ writes the initial pop to the open + lyt file and adds the closing //s """ + if (self.initial_pop_type == 'custom'): + lyt.write(' initial_pop\n') + for i in self.initial_pop: + lyt.write(' ' + str(int(i[0])) + ' ' + str(int(i[1])) + + ' ' + ' '.join([str(x) for x in i[2:]]) + + '\n') + else: + # TODO: test this part and fix, probably not functional currently + lyt.write(' initial_pop ' + self.initial_pop_type + + ' '.join([str(x) for x in self.initial_pop]) + + '\n') + lyt.write(r' //' + '\n') + lyt.write(r'//' + '\n') + + def __write_models_frictions_chunk(self, lyt): + """ writes the models frictions to the open + lyt file and adds the closing //s """ + if self.__models_frictions_flag: + print('HERE') + lyt.write(' modelsFriction\n') + for i in self.models_substrate_frictions: + print('HERE1') + lyt.write(' ' + str(self.models_substrate_frictions.index(i))+ ' ' + str(i) + '\n') + lyt.write(r' //' + '\n') + + def __write_models_pairs_frictions_chunk(self, lyt): + """ writes the intermodels pairs frictions to the open + lyt file and adds the closing //s """ + if self.__models_pairs_frictions_flag: + lyt.write(' interModelPairsFriction\n') + for i in range(len(self.models_pairs_frictions)): + for j in range(len(self.models_pairs_frictions)): + lyt.write(' ' + str(i) + ' ' +str(j) + ' ' + str(self.models_pairs_frictions[i][j]) + '\n') + lyt.write(r' //' + '\n') + + def update_models(self): + """ updates layout properties when models are added / removed + + This is usually just called internally. However, if you manually + change layout.models, then this method should be called to make other + necessary changes. + + """ + self.__build_initial_pop() + self.__build_exchanged_mets() + self.__add_new_mets_to_media() + + def __build_initial_pop(self): + # This counts how many models there are. then it goes through + # each model, and makes a new initial pop line of the right length + n_models = len(self.models) + initial_pop = [] + for i, m in enumerate(self.models): + if not isinstance(m.initial_pop[0], list): + m.initial_pop = [m.initial_pop] + for pop in m.initial_pop: + curr_line = [0] * (n_models + 2) + curr_line[0] = pop[0] + curr_line[1] = pop[1] + curr_line[i+2] = pop[2] + initial_pop.append(curr_line) + self.initial_pop = initial_pop + self.__resolve_initial_pop_line_dups() + + def __resolve_initial_pop_line_dups(self): + # sometimes each model has the same x,y, this resolves that + init_pop = self.initial_pop + init_pop_dict = {} + for row in init_pop: + if tuple(row[0:2]) in init_pop_dict.keys(): + init_pop_dict[tuple(row[0:2])] += np.array(row[2:]) + else: + init_pop_dict[tuple(row[0:2])] = np.array(row[2:]) + init_pop_fixed = [] + for key, value in init_pop_dict.items(): + loc = list(key) + loc.extend(list(value)) + init_pop_fixed.append(loc) + self.initial_pop = init_pop_fixed + + def __add_new_mets_to_media(self): + # usually run right after build_exchange mets, to add any new mets + # to the media data.frame + + for met in self.all_exchanged_mets: + if met not in self.media['metabolite'].values: + new_row = pd.DataFrame.from_dict({'metabolite': [met], + 'init_amount': [0], + 'diff_c': [self.default_diff_c], + 'g_static': [self.default_g_static], + 'g_static_val': [ + self.default_g_static_val], + 'g_refresh': [self.default_g_refresh]}) + self.media = pd.concat([self.media, new_row], + ignore_index=True, sort=True) + + def __build_exchanged_mets(self): + # goes through each model, grabs its exchange met names, and bundles + # them into a single list + all_exchanged_mets = [] + for m in self.models: + all_exchanged_mets.extend(m.get_exchange_metabolites()) + all_exchanged_mets = sorted(list(set(list(all_exchanged_mets)))) + self.all_exchanged_mets = all_exchanged_mets + + + def add_model(self, model): + """ + add a cometspy model to this layout + + This is the preferred way to add models to existing layout. + + Parameters + ---------- + + model : cometspy.model + a cometspy.model with an initial_pop + + Examples + -------- + + >>> import cobra.test + >>> import cometspy as c + >>> layout = c.layout() + >>> ecoli = cobra.test.create_test_model("ecoli") + >>> salmonella = cobra.test.create_test_model("salmonella") + >>> ecoli = c.model(ecoli) + >>> salmonella = c.model(salmonella) + >>> ecoli.open_exchanges() + >>> salmonella.open_exchanges() + >>> ecoli.initial_pop = [0, 0, 1.e-7] + >>> salmonella.initial_pop = [0, 0, 1.e-7] + >>> layout.add_model(ecoli) + >>> layout.add_model(salmonella) + + + """ + self.models.append(model) + self.update_models() + + def __check_if_initial_pops_in_range(self): + """ + checks if initial pops for each model are in the grid or raises error + + A common error is users setting model.initial_pop but not setting + layout.grid. Since they occur in different places, this method is + useful for a check before running. + + """ + for m in self.models: + for founder in m.initial_pop: + if (founder[0] < 0) or (founder[1] > (self.grid[1]-1)): + message = "the initial pop of a model is outside of layout.grid." + message += f" Either increase layout.grid or adjust {m.id}'s initial_pop" + raise ValueError(message) + + def __get_met_number(self, met): + """ returns the met number (of the external mets) given a name """ + met_number = [x for x in range(len(self.all_exchanged_mets)) if + self.all_exchanged_mets[x] == met][0] + return(met_number) + +def _read_file(filename): + f = open(filename, 'r') + f_lines = f.read() + f.close() + return f_lines diff --git a/build/lib/cometspy/model.py b/build/lib/cometspy/model.py new file mode 100644 index 0000000..9e5a41f --- /dev/null +++ b/build/lib/cometspy/model.py @@ -0,0 +1,1105 @@ +import pandas as pd +import os +import numpy as np +import re +import cobra +import io + + +class model: + """ a COMETS metabolic model to use in a layout + + A model contains information including the stoichiometric matrix, + reaction boundaries and Michaelis-Menten parameters, founder populations, + and various attributes related to motion in space. + + A model is usually initiated by supplying a cobrapy model, and to be used + in a COMETS simulation must minimally be given an initial population: + + Parameters + ---------- + + model : cobra.Model or str, optional + Either a cobrapy Model, a path to a .cmd COMETS model, a path to a + cobra sbml model which could be loaded with cobra.io.read_sbml_model() + or None. + + + Attributes + ---------- + + initial_pop : list + list of lists of spatial founder pops. e.g.[[x, y, grams], [x,y,grams]] + id : str + the name of the model. should be unique. + reactions : pandas.DataFrame + DataFrame containing reaction information including boundaries. + smat : pandas.DataFrame + DataFrame containing the stoichiometric matrix + metabolites : pandas.DataFrame + single-column DataFrame containing the names of the metabolites + signals : pandas.Dataframe + DataFrame on reactions whose bounds depend on metabolite concentrations + default_vmax : float + default uptake vmax in mmol / gDW /hr for Michaelis-Menten kinetics + default_km : float + default half-max metabolite concentration (M) for M-M kinetics + default_hill : float + default hill coefficient for hill-like uptake kinetics + default_bounds : list + list of two floats that are the default lower and upper reaction bounds + obj_style : str + one of MAXIMIZE_OBJECTIVE_FLUX (fba) or MAX_OBJECTIVE_MIN_TOTAL (pfba) + optimizer : str + one of "GUROBI" or "GLPK". not all functionality works with GLPK + + Examples + -------- + + >>> import cobra.test + >>> import cometspy as c + >>> ecoli = cobra.test.create_test_model("ecoli") + >>> model = c.model(ecoli) + >>> model.initial_pop = [0, 0, 1.e-12] # puts 1.e-12 gDW of biomass at 0,0 + >>> print(model.id) + iJO1366 + + """ + def __init__(self, model : cobra.Model = None, randomtag : bool = False): + self.initial_pop = [[0, 0, 0.0]] + + if randomtag: + self.id = '_' + hex(id(self)) + else: + self.id = "" + + self.reactions = pd.DataFrame(columns=['REACTION_NAMES', 'ID', + 'LB', 'UB', 'EXCH', + 'EXCH_IND', 'V_MAX', + 'KM', 'HILL']) + self.smat = pd.DataFrame(columns=['metabolite', + 'rxn', + 's_coef']) + self.metabolites = pd.DataFrame(columns=['METABOLITE_NAMES']) + self.signals = pd.DataFrame(columns=['REACTION_NUMBER', + 'EXCH_IND', + 'BOUND', + 'FUNCTION', + 'PARAMETERS', + 'REACTION_NAMES', 'EXCH'], + dtype=object) + # multitoxins is a special case of signaling + self.multitoxins = pd.DataFrame(columns = ['REACTION_NUMBER', + 'EXCH_INDS', + 'BOUND', + 'KMS', + 'HILLS', + 'VMAX', + 'REACTION_NAME', + 'EXCH_NAMES'], + dtype = object) + self.light = [] + + self.vmax_flag = False + self.km_flag = False + self.hill_flag = False + self.convection_flag = False + self.light_flag = False + + self.nonlinear_diffusion_flag = False + self.neutral_drift_flag = False + self.noise_variance_flag = False + self.default_vmax = 10 + self.default_km = 1 + self.default_hill = 1 + self.default_bounds = [0, 1000] + self.objective = None + self.optimizer = 'GUROBI' + self.obj_style = 'MAXIMIZE_OBJECTIVE_FLUX' + + if model is not None: + if isinstance(model, cobra.Model): + self.load_cobra_model(model, randomtag) + else: # assume it is a path + if model[-3:] == "cmd": + self.read_comets_model(model, randomtag) + else: + self.read_cobra_model(model, randomtag) + + def get_reaction_names(self) -> list: + """ returns a list of reaction names""" + return(list(self.reactions['REACTION_NAMES'])) + + def add_multitoxin(self, rxn_num, exch_ids, bound, + vmax, kms, hills): + """ + Adds a signaling relationship where >= 1 signal close a bound + + This incorporates a hill-function reduction of (usually) an upper bound + It is useful for incorporating multiple, additively-functioning signals. + + bound = vmax * MULT((km_i^h_i) / (km_i^h_i + M_i^h_i)) + where MULT is a multiplicative function, M_i is the molarity of the + metabolite in the given exch_id, and km_i and h_i are the km and h for + for that signal. + + Parameters + ---------- + rxn_num : int + number of the affected reaction + exch_ids : list[int] + exchange numbers of the signals + bound : str + "ub" or "lb" + vmax : float + maximum bound in absence of signals + kms : list[float] + kms for the signals + hills : list[float] + hill coefficients for the signals + + Returns + ------- + None. + + """ + rxn_name = self.reactions.loc[self.reactions.ID == rxn_num, + 'REACTION_NAMES'].item() + rxn_num = str(rxn_num) + exch_names = [list(self.get_exchange_metabolites())[exch_id-1] for + exch_id in exch_ids] + new_row = pd.DataFrame({'REACTION_NUMBER': rxn_num, + 'EXCH_INDS': 1, + 'BOUND': bound, + 'KMS': 1, + 'HILLS': 1, + 'VMAX' :vmax, + 'REACTION_NAME': rxn_name, + 'EXCH_NAMES': ""}, + index=[0], + dtype=object) + new_row.loc[0, 'KMS'] = kms + new_row.loc[0, 'HILLS'] = hills + new_row.loc[0, "EXCH_INDS"] = exch_ids + new_row.loc[0, 'EXCH_NAMES'] = exch_names + self.multitoxins = self.multitoxins.append(new_row, ignore_index=True) + + def add_signal(self, rxn_num, exch_id, bound, + function, parms): + """adds a signal to the reaction rxn_num that changes bounds + + Parameters + ---------- + rxn_num : int or 'death' + reaction affected (int) or whether the signal causes death 'death' + exch_id : int + the number of the exchange reaction. see model.metabolites + bound : str 'ub' or 'lb' + specifies whether the upper or lower bound is affected + function : str 'linear' or 'bounded_linear' or 'generalized_logistic' + specifies the function relating the metabolite conc. to the bound + parms : list (float) + a list of floats for the function + + """ + + if str(rxn_num).lower().strip() == 'death': + rxn_name = 'death' + rxn_num = 'death' + else: + rxn_name = self.reactions.loc[self.reactions.ID == rxn_num, + 'REACTION_NAMES'].item() + rxn_num = str(rxn_num) + # note: rxn_num matches the pandas DataFrame in this object. + # however, COMETS actually wants that rxn_num - 1, which + # is why it does that during the saving process, and why + # you will see it that way in saved model files. + exch_name = list(self.get_exchange_metabolites())[exch_id-1] + new_row = pd.DataFrame({'REACTION_NUMBER': rxn_num, + 'EXCH_IND': exch_id, + 'BOUND': bound, + 'FUNCTION': function, + 'PARAMETERS': 1, + 'REACTION_NAMES': rxn_name, + 'EXCH': exch_name}, + index=[0], + dtype=object) + new_row.loc[0, 'PARAMETERS'] = parms + self.signals = self.signals.append(new_row, ignore_index=True) + + def add_neutral_drift_parameter(self, neutralDriftSigma): + """ sets the neutralDriftSigma parameter """ + if not isinstance(neutralDriftSigma, float): + raise ValueError("neutralDriftSigma must be a float") + self.neutral_drift_flag = True + self.neutralDriftSigma = neutralDriftSigma + + def add_nonlinear_diffusion_parameters(self, + zero : float=1., + n : float =1., + exponent : float=1., + hilln : float=10., + hillk : float=0.9): + """ sets the model to use non-linear diffusion biomass spread. + + This also requires one set the biomassMotionStyle to + 'ConvNonlin Diffusion 2D' in the comets.params object (see Example) + + Parameters + ---------- + zero : float, optional + n : float, optional + exponent : float, optional + hilln : float, optional + hillk : float, optional + + Example + -------- + + >>> import cobra.test + >>> import cometspy as c + >>> model = c.model(cobra.test.create_test_model("ecoli")) + >>> model.add_nonlinear_diffusion_parameters(1., 1., 2., 5., 0.5) + >>> params = c.params() + >>> params.set_param("biomassMotionStyle", "ConvNonlin Diffusion 2D") + + """ + + for parm in [zero, n, + exponent, hilln, + hillk]: + if not isinstance(parm, float): + raise ValueError('all nonlinear diffusion terms must be float') + self.nonlinear_diffusion_flag = True + self.nonlinear_diffusion_parameters = {'convNonLinDiffZero': zero, + 'convNonlinDiffN': n, + 'convNonlinDiffExponent': exponent, + 'convNonlinDiffHillN': hilln, + 'convNonlinDiffHillK': hillk} + + def add_light(self, reaction : str, + abs_coefficient : float, + abs_base : float): + """Causes a reaction to function in response to light + + Parameters + ---------- + + reaction : str + the name of the reaction affected + abs_coefficient : float + the absorption coefficient + abs_base : float + absorption baseline + + """ + if (reaction not in self.reactions['REACTION_NAMES'].values): + raise ValueError('the reaction is not present in the model') + self.light.append([reaction, abs_coefficient, abs_base]) + self.light_flag = True + + def add_convection_parameters(self, packedDensity : float =1., + elasticModulus : float =1., + frictionConstant : float =1., + convDiffConstant : float =1.): + """ sets the parameters for biomass spread to use convection + + This also requires one set the biomassMotionStyle to + 'Convection 2D' in the comets.params object (see Example) + + Parameters + ---------- + + packedDensity : float, optional + elasticModulus : float, optional + frictionConstant : float, optional + convDiffConstant : float, optional + + Example + ------- + + import cobra.test + import cometspy as c + model = c.model(cobra.test.create_test_model("ecoli")) + model.add_convection_parameters(1., 0.5, 2., 10.e-5) + params = c.params() + params.set_param("biomassMotionStyle", "Convection 2D") + + """ + if not isinstance(packedDensity, float): + raise ValueError('packed_density must be a float') + if not isinstance(elasticModulus, float): + raise ValueError('elasticModulus must be a float') + if not isinstance(frictionConstant, float): + raise ValueError('frictionConstant must be a float') + if not isinstance(convDiffConstant, float): + raise ValueError('convDiffConstant must be a float') + self.convection_flag = True + self.convection_parameters = {'packedDensity': packedDensity, + 'elasticModulus': elasticModulus, + 'frictionConstant': frictionConstant, + 'convDiffConstant': convDiffConstant} + + def add_noise_variance_parameter(self, noiseVariance : float): + """ sets the noise variance parameter + + Parameters + ---------- + noiseVariance : float + + """ + if not isinstance(noiseVariance, float): + raise ValueError('noiseVariance must be a float') + self.noise_variance_flag = True + self.noise_variance = noiseVariance + + def ensure_sinks_are_not_exchanges(self, suffix : str = '_c'): + """ set exchange reactions ending in suffix (default = '_c') to sink + + Parameters + ---------- + suffix : str, optional + the suffix to look for in exchange reactions. default = '_c' + + """ + rxn_names = [rxn for rxn in self.reactions.REACTION_NAMES.to_list() if rxn[-2:] == suffix] + for rxn in rxn_names: + self.reactions.loc[self.reactions.REACTION_NAMES == rxn, 'EXCH'] = False + self.reactions.loc[self.reactions.REACTION_NAMES == rxn, 'EXCH_IND'] = 0 + + def open_exchanges(self, lower_bound : float = -1000., upper_bound : float = 1000.): + """ sets all exchange boundaries to the specifed values + + This method is useful when first making a COMETS model so that medium + definitions are not carried over from cobra into COMETS. + + Parameters + ---------- + + lower_bound : float, optional + default = -1000. in units of mmol / gDW / hr + upper_bound : float, optional + default = 1000. in units of mmol / gDW / hr + + """ + + self.reactions.loc[self.reactions.EXCH,'LB'] = lower_bound + self.reactions.loc[self.reactions.EXCH,'UB'] = upper_bound + + def get_exchange_metabolites(self) -> list: + """ returns a list of the names of the exchange metabolites """ + exchmets = pd.merge(self.reactions.loc[self.reactions['EXCH'], 'ID'], + self.smat, + left_on='ID', right_on='rxn', + how='inner')['metabolite'] + exchmets = self.metabolites.iloc[exchmets-1] + return(exchmets.METABOLITE_NAMES) + + def change_bounds(self, reaction : str, lower_bound : float, upper_bound : float): + """ changes the bounds of the specified reaction + + Parameters + ---------- + + reaction : str + the name of the reaction + lower_bound : float + the new value of the reaction's lower bound in mmol / gDW / hr + upper_bound : float + the new value of the reaction's upper bound in mmol / gDW / hr + + """ + if reaction not in self.reactions['REACTION_NAMES'].values: + print('reaction couldnt be found') + return + self.reactions.loc[self.reactions['REACTION_NAMES'] == reaction, + 'LB'] = lower_bound + self.reactions.loc[self.reactions['REACTION_NAMES'] == reaction, + 'UB'] = upper_bound + + def get_bounds(self, reaction : str) -> tuple: + """returns a tuple (lb, ub) for the specified reaction + + Parameters + ---------- + + reaction : str + the name of the reaction + + """ + if reaction not in self.reactions['REACTION_NAMES'].values: + print('reaction couldnt be found') + return + lb = float(self.reactions.loc[self.reactions[ + 'REACTION_NAMES'] == reaction, 'LB']) + ub = float(self.reactions.loc[self.reactions[ + 'REACTION_NAMES'] == reaction, 'UB']) + return((lb, ub)) + + def change_vmax(self, reaction : str, vmax : float): + """ changes the vmax of the specified reaction. + + Parameters + ---------- + + reaction : str + the name of the reaction + vmax : float + the new vmax value for the reaction, for Monod (Michaelis-Menten) + kinetics + + """ + if reaction not in self.reactions['REACTION_NAMES'].values: + print('reaction couldnt be found') + return + self.vmax_flag = True + self.reactions.loc[self.reactions[ + 'REACTION_NAMES'] == reaction, 'V_MAX'] = vmax + + def change_km(self, reaction, km): + """ changes the km of the specified reaction. + + Parameters + ---------- + + reaction : str + the name of the reaction + km : float + the new km value for the reaction, for Monod (Michaelis-Menten) + kinetics + + """ + if reaction not in self.reactions['REACTION_NAMES'].values: + print('reaction couldnt be found') + return + self.km_flag = True + self.reactions.loc[self.reactions[ + 'REACTION_NAMES'] == reaction, 'KM'] = km + + def change_hill(self, reaction, hill): + """ changes the hill coefficient of the specified reaction. + + Parameters + ---------- + + reaction : str + the name of the reaction + hill : float + the new hill coef. for the reaction, for Hill-like kinetics + + """ + if reaction not in self.reactions['REACTION_NAMES'].values: + print('reaction couldnt be found') + return + self.hill_flag = True + self.reactions.loc[self.reactions[ + 'REACTION_NAMES'] == reaction, 'HILL'] = hill + + def change_optimizer(self, optimizer): + """ changes the optimizer to a specified one. + + Parameters + ---------- + + optimizer : str + the name of the optimizer + + """ + self.optimizer = optimizer + + def read_cobra_model(self, path : str, randomtag : bool = False): + """ reads a cobra model from a file and loads it into this model + + This is an alternative way to initialize a COMETS model, by supplying + it with the path to a cobra model that can be read using + cobra.io.read_sbml_model. + + Parameters + ---------- + + path : str + path to a cobra model in sbml format + + Example + ------- + + >>> import cometspy as c + >>> path_to_file = "./my_new_model.xml" + >>> model = c.model() + >>> model.read_cobra_model(path_to_file) + + """ + curr_m = cobra.io.read_sbml_model(path) + self.load_cobra_model(curr_m, randomtag) + + def load_cobra_model(self, curr_m : cobra.Model, randomtag : bool = False): + """ creates the COMETS model from the supplied cobra model + + This is usually used internally when creating a COMETS model. + + Parameters + ---------- + + curr_m : cobra.Model + the cobra model object to be converted to the COMETS model object + + Example + ------- + + >>> import cometspy as c + >>> import cobra.test + >>> ecoli = cobra.test.create_test_model('ecoli') + >>> model = c.model() + >>> model.load_cobra_model(ecoli) + + """ + self.id = curr_m.id + if randomtag: + self.id = self.id + '_' + hex(id(self)) + + # reactions and their features + reaction_list = curr_m.reactions + self.reactions['REACTION_NAMES'] = [str(x).split(':')[0] for + x in reaction_list] + self.reactions['ID'] = [k for k in + range(1, len(reaction_list)+1)] + self.reactions['LB'] = [x.lower_bound for x in reaction_list] + self.reactions['UB'] = [x.upper_bound for x in reaction_list] + + self.reactions['EXCH'] = [True if (len(k.metabolites) == 1) & + (list(k.metabolites. + values())[0] == (-1)) & + ('DM_' not in k.id) + else False for k in reaction_list] + + exch = self.reactions.loc[self.reactions['EXCH'], 'ID'].tolist() + self.reactions['EXCH_IND'] = [exch.index(x)+1 + if x in exch else 0 + for x in self.reactions['ID']] + + self.reactions['V_MAX'] = [k.Vmax + if hasattr(k, 'Vmax') + else float('NaN') + for k in reaction_list] + + if not self.reactions.V_MAX.isnull().all(): + self.vmax_flag = True + + self.reactions['KM'] = [k.Km + if hasattr(k, 'Km') + else float('NaN') + for k in reaction_list] + + if not self.reactions.KM.isnull().all(): + self.km_flag = True + + self.reactions['HILL'] = [k.Hill + if hasattr(k, 'Hill') + else float('NaN') + for k in reaction_list] + + if not self.reactions.HILL.isnull().all(): + self.hill_flag = True + + if self.vmax_flag: + if hasattr(curr_m, 'default_vmax'): + self.default_vmax = curr_m.default_vmax + + if self.km_flag: + if hasattr(curr_m, 'default_km'): + self.default_km = curr_m.default_km + + if self.hill_flag: + if hasattr(curr_m, 'default_hill'): + self.default_hill = curr_m.default_hill + + # Metabolites + metabolite_list = curr_m.metabolites + self.metabolites['METABOLITE_NAMES'] = [str(x) for + x in metabolite_list] + + # S matrix + for index, row in self.reactions.iterrows(): + rxn = curr_m.reactions.get_by_id( + row['REACTION_NAMES']) + rxn_num = row['ID'] + rxn_mets = [1+list(self.metabolites[ + 'METABOLITE_NAMES']).index( + x.id) for x in rxn.metabolites] + met_s_coefs = list(rxn.metabolites.values()) + + cdf = pd.DataFrame({'metabolite': rxn_mets, + 'rxn': [rxn_num]*len(rxn_mets), + 's_coef': met_s_coefs}) + cdf = cdf.sort_values('metabolite') + self.smat = pd.concat([self.smat, cdf]) + + self.smat = self.smat.sort_values(by=['metabolite', 'rxn']) + + # The rest of stuff + if hasattr(curr_m, 'default_bounds'): + self.default_bounds = curr_m.default_bounds + + obj = {str(x).split(':')[0]:x.objective_coefficient + for x in reaction_list + if x.objective_coefficient != 0} + obj = {rx: -1 if coef < 0 else 1 for rx, coef in obj.items()} + + self.objective = [int(self.reactions[self.reactions. + REACTION_NAMES == rx]['ID']) * coef for rx, coef in obj.items()] + + if hasattr(curr_m, 'comets_optimizer'): + self.optimizer = curr_m.comets_optimizer + + if hasattr(curr_m, 'comets_obj_style'): + self.obj_style = curr_m.comets_obj_style + + def read_comets_model(self, path : str, randomtag : bool = False): + """ an alternative way to create a COMETS model by reading from a file + + This allows a user to load previously-saved COMETS models. The + contents populate this object's attributes. + + Parameters + ---------- + + path : str + the path to the COMETS model file + + Example + ------- + + >>> import cometspy as c + >>> path_to_model = "./iJO1366.cmd" # this must actually exist + >>> model = c.model() + >>> model.read_comets_model(path_to_model) + + """ + self.id = os.path.splitext(os.path.basename(path))[0] + self.id + if randomtag: + self.id = self.id + '_' + hex(id(self)) + + # in this way, its robust to empty lines: + m_f_lines = [s for s in _read_file(path).splitlines() if s] + m_filedata_string = os.linesep.join(m_f_lines) + ends = [] + for k in range(0, len(m_f_lines)): + if '//' in m_f_lines[k]: + ends.append(k) + + # '''----------- S MATRIX ------------------------------''' + lin_smat = re.split('SMATRIX', + m_filedata_string)[0].count('\n') + lin_smat_end = next(x for x in ends if x > lin_smat) + + self.smat = pd.read_csv(io.StringIO('\n'.join(m_f_lines[ + lin_smat:lin_smat_end])), + delimiter=r'\s+', + skipinitialspace=True) + self.smat.columns = ['metabolite', 'rxn', 's_coef'] + + # '''----------- REACTIONS AND BOUNDS-------------------''' + lin_rxns = re.split('REACTION_NAMES', + m_filedata_string)[0].count('\n') + lin_rxns_end = next(x for x in + ends if x > lin_rxns) + + rxn = pd.read_csv(io.StringIO('\n'.join(m_f_lines[ + lin_rxns:lin_rxns_end])), + delimiter=r'\s+', + skipinitialspace=True) + + rxn['ID'] = range(1, len(rxn)+1) + + lin_bnds = re.split('BOUNDS', + m_filedata_string)[0].count('\n') + lin_bnds_end = next(x for x in ends if x > lin_bnds) + + bnds = pd.read_csv(io.StringIO('\n'.join(m_f_lines[ + lin_bnds:lin_bnds_end])), + delimiter=r'\s+', + skipinitialspace=True) + + default_bounds = [float(bnds.columns[1]), + float(bnds.columns[2])] + + bnds.columns = ['ID', 'LB', 'UB'] + reactions = pd.merge(rxn, bnds, + left_on='ID', right_on='ID', + how='left') + reactions.LB.fillna(default_bounds[0], inplace=True) + reactions.UB.fillna(default_bounds[1], inplace=True) + + # '''----------- METABOLITES ---------------------------''' + lin_mets = re.split('METABOLITE_NAMES', + m_filedata_string)[0].count('\n') + lin_mets_end = next(x for x in ends if x > lin_mets) + + metabolites = pd.read_csv(io.StringIO('\n'.join(m_f_lines[ + lin_mets:lin_mets_end])), + delimiter=r'\s+', + skipinitialspace=True) + + # '''----------- EXCHANGE RXNS -------------------------''' + lin_exch = re.split('EXCHANGE_REACTIONS', + m_filedata_string)[0].count('\n')+1 + exch = [int(k) for k in re.findall(r'\S+', + m_f_lines[lin_exch]. + strip())] + + reactions['EXCH'] = [True if x in exch else False + for x in reactions['ID']] + reactions['EXCH_IND'] = [exch.index(x)+1 + if x in exch else 0 + for x in reactions['ID']] + + # '''----------- VMAX VALUES --------------------------''' + if 'VMAX_VALUES' in m_filedata_string: + self.vmax_flag = True + lin_vmax = re.split('VMAX_VALUES', + m_filedata_string)[0].count('\n') + lin_vmax_end = next(x for x in ends if x > lin_vmax) + + Vmax = pd.read_csv(io.StringIO('\n'.join(m_f_lines[ + lin_vmax:lin_vmax_end])), + delimiter=r'\s+', + skipinitialspace=True) + + Vmax.columns = ['EXCH_IND', 'V_MAX'] + + reactions = pd.merge(reactions, Vmax, + left_on='EXCH_IND', + right_on='EXCH_IND', + how='left') + self.default_vmax = float(m_f_lines[lin_vmax-1].split()[1]) + else: + reactions['V_MAX'] = np.NaN + + # '''----------- VMAX VALUES --------------------------''' + if 'KM_VALUES' in m_filedata_string: + self.km_flag = True + lin_km = re.split('KM_VALUES', + m_filedata_string)[0].count('\n') + lin_km_end = next(x for x in ends if x > lin_km) + + Km = pd.read_csv(io.StringIO('\n'.join(m_f_lines[ + lin_km:lin_km_end])), + delimiter=r'\s+', + skipinitialspace=True) + Km.columns = ['EXCH_IND', 'KM'] + + reactions = pd.merge(reactions, Km, + left_on='EXCH_IND', + right_on='EXCH_IND', + how='left') + self.default_km = float(m_f_lines[lin_km-1].split()[1]) + else: + reactions['KM'] = np.NaN + + # '''----------- VMAX VALUES --------------------------''' + if 'HILL_COEFFICIENTS' in m_filedata_string: + self.hill_flag = True + lin_hill = re.split('HILL_COEFFICIENTS', + m_filedata_string)[0].count('\n') + lin_hill_end = next(x for x in ends if x > lin_hill) + + Hill = pd.read_csv(io.StringIO('\n'.join(m_f_lines[ + lin_hill:lin_hill_end])), + delimiter=r'\s+', + skipinitialspace=True) + Hill.columns = ['EXCH_IND', 'HILL'] + + reactions = pd.merge(reactions, Hill, + left_on='EXCH_IND', + right_on='EXCH_IND', + how='left') + self.default_hill = float(m_f_lines[lin_hill-1].split()[1]) + else: + reactions['HILL'] = np.NaN + + # '''----------- OBJECTIVE -----------------------------''' + lin_obj = re.split('OBJECTIVE', + m_filedata_string)[0].count('\n')+1 + self.objective = int(m_f_lines[lin_obj].strip()) + + # '''----------- OBJECTIVE STYLE -----------------------''' + if 'OBJECTIVE_STYLE' in m_filedata_string: + lin_obj_st = re.split('OBJECTIVE_STYLE', + m_filedata_string)[0].count( + '\n')+1 + self.obj_style = m_f_lines[lin_obj_st].strip() + + # '''----------- OPTIMIZER -----------------------------''' + if 'OPTIMIZER' in m_filedata_string: + lin_opt = re.split('OPTIMIZER', + m_filedata_string)[0].count('\n') + self.optimizer = m_f_lines[lin_opt].split()[1] + + # '''--------------neutral drift------------------------''' + if "neutralDrift" in m_filedata_string: + lin_obj_st = re.split('neutralDrift', + m_filedata_string)[0].count( + '\n') + if "TRUE" == m_f_lines[lin_obj_st].strip().split()[1].upper(): + self.neutral_drift_flag = True + self.neutralDriftSigma = 0. + if "neutralDriftsigma" in m_filedata_string: + lin_opt = re.split('neutralDriftsigma', + m_filedata_string)[0].count('\n') + self.neutralDriftSigma = float(m_f_lines[lin_opt].split()[1]) + + # '''--------------convection---------------------------''' + for parm in ['packedDensity', 'elasticModulus', + 'frictionConstant', 'convDiffConstant']: + if parm in m_filedata_string: + lin_obj_st = re.split(parm, + m_filedata_string)[0].count( + '\n') + parm_value = float(m_f_lines[lin_obj_st].strip().split()[1]) + try: + self.convection_parameters[parm] = parm_value + except: # TODO change bare except statements to ifelse + self.convection_flag = True + self.convection_parameters = {'packedDensity': 1., + 'elasticModulus': 1., + 'frictionConstant': 1., + 'convDiffConstant': 1.} + self.convection_parameters[parm] = parm_value + + # '''--------------non-linear diffusion---------------------------''' + for parm in ['convNonLinDiffZero', 'convNonlinDiffN', 'convNonlinDiffExponent', + 'convNonlinDiffHillN', 'convNonlinDiffHillK']: + if parm in m_filedata_string: + lin_obj_st = re.split(parm, + m_filedata_string)[0].count( + '\n') + parm_value = float(m_f_lines[lin_obj_st].strip().split()[1]) + try: + self.nonlinear_diffusion_parameters[parm] = parm_value + except: # TODO change bare except statements to ifelse + self.nonlinear_diffusion_flag = True + self.nonlinear_diffusion_parameters = {'convNonLinDiffZero': 1., + 'convNonlinDiffN': 1., + 'convNonlinDiffExponent': 1., + 'convNonlinDiffHillN': 10., + 'convNonlinDiffHillK': .9} + self.nonlinear_diffusion_parameters[parm] = parm_value + + # '''-----------noise variance-----------------''' + if 'noiseVariance' in m_filedata_string: + lin_obj_st = re.split('noiseVariance', + m_filedata_string)[0].count( + '\n') + noiseVariance = float(m_f_lines[lin_obj_st].strip().split()[1]) + + self.noise_variance_flag = True + self.noise_variance = noiseVariance + # assign the dataframes we just built + self.reactions = reactions + self.metabolites = metabolites + + def delete_comets_model(self, working_dir = None): + """ deletes a file version of this model if it exists. + + Parameters + ---------- + + working_dir : str, optional + the directory where the comets model file exists + + """ + path_to_delete = "" + if working_dir is not None: + path_to_delete = working_dir + path_to_delete = path_to_delete + self.id + '.cmd' + os.remove(path_to_delete) + + def write_comets_model(self, working_dir : str = None): + """ writes the COMETS model object to a file + + This writes the current model object in COMETS format to file, either + in the current working directory or in the directory specified by + the optional parameter working_dir. It is mostly used by the comets + object to run simulations, though a user may with to write these to + examine directly or to save. + + Parameters + ---------- + working_dir : str, optional + a directory path to put COMETS files. for example "./data_files/" + + """ + path_to_write = "" + if working_dir is not None: + path_to_write = working_dir + path_to_write = path_to_write + self.id + '.cmd' + + # format variables for writing comets model + bnd = self.reactions.loc[(self.reactions['LB'] + != self.default_bounds[0]) | + (self.reactions['UB'] != + self.default_bounds[1]), + ['ID', 'LB', 'UB']].astype( + str).apply(lambda x: ' '.join(x), + axis=1) + bnd = ' ' + bnd.astype(str) + + rxn_n = ' ' + self.reactions['REACTION_NAMES'].astype(str) + + met_n = ' ' + self.metabolites.astype(str) + + smat = self.smat.astype(str).apply(lambda x: + ' '.join(x), axis=1) + smat = ' ' + smat.astype(str) + + exch_r = ' '.join([str(x) for x in + self.reactions.loc[self.reactions.EXCH, 'ID']]) + + # optional fields (vmax,km, hill) + if self.vmax_flag: + Vmax = self.reactions.loc[self.reactions['V_MAX'].notnull(), + ['EXCH_IND', 'V_MAX']] + Vmax = Vmax.astype(str).apply(lambda x: + ' '.join(x), axis=1) + Vmax = ' ' + Vmax.astype(str) + + if self.km_flag: + Km = self.reactions.loc[self.reactions['KM'].notnull(), + ['EXCH_IND', 'KM']] + Km = Km.astype(str).apply(lambda x: + ' '.join(x), axis=1) + Km = ' ' + Km.astype(str) + + if self.hill_flag: + Hill = self.reactions.loc[self.reactions['HILL'].notnull(), + ['EXCH_IND', 'HILL']] + Hill = Hill.astype(str).apply(lambda x: + ' '.join(x), axis=1) + Hill = ' ' + Hill.astype(str) + + if os.path.isfile(path_to_write): + os.remove(path_to_write) + + with open(path_to_write, 'a') as f: + + f.write('SMATRIX ' + str(len(self.metabolites)) + + ' ' + str(len(self.reactions)) + '\n') + smat.to_csv(f, mode='a', lineterminator = '\n', header=False, index=False) + f.write(r'//' + '\n') + + f.write('BOUNDS ' + + str(self.default_bounds[0]) + ' ' + + str(self.default_bounds[1]) + '\n') + bnd.to_csv(f, mode='a', lineterminator = '\n', header=False, index=False) + f.write(r'//' + '\n') + + f.write('OBJECTIVE\n' + + ' ' + ' '.join([str(rx) for rx in self.objective]) + '\n') + f.write(r'//' + '\n') + + f.write('METABOLITE_NAMES\n') + met_n.to_csv(f, mode='a', lineterminator = '\n', header=False, index=False) + f.write(r'//' + '\n') + + f.write('REACTION_NAMES\n') + rxn_n.to_csv(f, mode='a', lineterminator = '\n', header=False, index=False) + f.write(r'//' + '\n') + + f.write('EXCHANGE_REACTIONS\n') + f.write(' ' + exch_r + '\n') + f.write(r'//' + '\n') + + if self.vmax_flag: + f.write('VMAX_VALUES ' + + str(self.default_vmax) + '\n') + Vmax.to_csv(f, mode='a', lineterminator = '\n', header=False, index=False) + f.write(r'//' + '\n') + + if self.km_flag: + f.write('KM_VALUES ' + + str(self.default_km) + '\n') + Km.to_csv(f, mode='a', lineterminator = '\n', header=False, index=False) + f.write(r'//' + '\n') + + if self.hill_flag: + f.write('HILL_VALUES ' + + str(self.default_hill) + '\n') + Hill.to_csv(f, mode='a', lineterminator = '\n', header=False, index=False) + f.write(r'//' + '\n') + + if self.light_flag: + f.write('LIGHT\n') + for lrxn in self.light: + lrxn_ind = str(int(self.reactions.ID[ + self.reactions['REACTION_NAMES'] == lrxn[0]])) + f.write(' {} {} {}\n'.format(lrxn_ind, + lrxn[1], lrxn[2])) + f.write(r'//' + '\n') + + if self.signals.size > 0 or self.multitoxins.size > 0: + f.write('MET_REACTION_SIGNAL\n') + if self.signals.size > 0: + sub_signals = self.signals.drop(['REACTION_NAMES', 'EXCH'], + axis='columns') + col_names = list(self.signals.drop(['REACTION_NAMES', + 'EXCH', 'PARAMETERS'], + axis='columns').columns) + for idx in sub_signals.index: + row = sub_signals.drop(['PARAMETERS'], axis='columns').iloc[idx, :] + n_parms = len(sub_signals.PARAMETERS[idx]) + curr_col_names = col_names + [str(i) for i in range(n_parms)] + temp_df = pd.DataFrame(columns=curr_col_names) + if row.loc['REACTION_NUMBER'] != 'death': + row.loc['REACTION_NUMBER'] = str(int(row.loc['REACTION_NUMBER'])) + temp_df.loc[0, 'REACTION_NUMBER'] = row.loc['REACTION_NUMBER'] + temp_df.loc[0, 'EXCH_IND'] = row.loc['EXCH_IND'] + temp_df.loc[0, 'BOUND'] = row.loc['BOUND'] + temp_df.loc[0, 'FUNCTION'] = row.loc['FUNCTION'] + for i in range(n_parms): + temp_df.loc[0, str(i)] = sub_signals.PARAMETERS[idx][i] + temp_df.to_csv(f, mode='a', lineterminator = '\n', sep=' ', header=False, index=False) + if self.multitoxins.size > 0: + for idx in self.multitoxins.index: + rxn_num = self.multitoxins.loc[idx,"REACTION_NUMBER"] + bound = self.multitoxins.loc[idx,"BOUND"] + exch_inds = ','.join([str(ind) for ind in self.multitoxins.loc[idx,"EXCH_INDS"]]) + vmax = self.multitoxins.loc[idx, 'VMAX'] + kms = ','.join([str(ind) for ind in self.multitoxins.loc[idx,"KMS"]]) + hills = ','.join([str(ind) for ind in self.multitoxins.loc[idx,"HILLS"]]) + curr_line = f"multitoxin {rxn_num} {exch_inds} {bound} {vmax} {kms} {hills}" + f.write(curr_line + '\n') + f.write(r'//' + '\n') + + if self.convection_flag: + for key, value in self.convection_parameters.items(): + f.write(key + ' ' + str(value) + '\n') + f.write(r'//' + '\n') + + if self.nonlinear_diffusion_flag: + for key, value in self.nonlinear_diffusion_parameters.items(): + f.write(key + ' ' + str(value) + '\n') + f.write(r'//' + '\n') + + if self.noise_variance_flag: + f.write('noiseVariance' + ' ' + + str(self.noise_variance) + '\n') + f.write(r'//' + '\n') + + if self.neutral_drift_flag: + f.write("neutralDrift true\n//\n") + f.write("neutralDriftSigma " + str(self.neutralDriftSigma) + "\n//\n") + + f.write('OBJECTIVE_STYLE\n' + self.obj_style + '\n') + f.write(r'//' + '\n') + + f.write('OPTIMIZER ' + self.optimizer + '\n') + f.write(r'//' + '\n') + + +def _read_file(filename: str) -> str: + """ helper function to read non-rectangular files. + """ + f = open(filename, 'r') + f_lines = f.read() + f.close() + return f_lines diff --git a/build/lib/cometspy/params.py b/build/lib/cometspy/params.py new file mode 100644 index 0000000..4c06b8e --- /dev/null +++ b/build/lib/cometspy/params.py @@ -0,0 +1,399 @@ +import pandas as pd +import os + + +def __isfloat(value): + try: + float(value) + return True + except ValueError: + return False + + +class params: + ''' + object storing simulation-related and some default biological parameters + + The params object is an essential object needed to start a comets + simulation, along with a layout object containing models. The params + object stores simulation information such as timeStep and maxCycles, as + well as many default biological parameters such as defaultKm and + defaultDiffC (the default metabolite diffusion constant). + + All of the parameters are stored in a dictionary called 'all_params' and + can either be adjusted directly or using set_param(). + + When params are written to a COMETS object for COMETS to process, two + files are generated ('.current_global' and 'current_package') + + Parameters + ---------- + global_params : str, optional + optional path to an existing 'global_params' file to load + package_params : str, optional + optional path to an existing 'package_params' file to load + + Examples + -------- + + >>> # make a params object and change two params + >>> import cometspy as c + >>> params = c.params() + >>> params.set_param("timeStep", 0.01) # hours + >>> params.set_param("spaceWidth", 0.5) # cm + >>> # view all the params + >>> params.show_params() + >>> # access the params directly + >>> params.all_params # this is a dict + >>> # use a params object in a comets simulation + >>> import cobra.test + >>> model = c.model(cobra.test.create_test_model("textbook")) + >>> model.initial_pop = [0, 0, 1.e-7] + >>> model.open_exchanges() + >>> layout = c.layout([model]) + >>> # normally we'd now alter the media in the layout + >>> layout.add_typical_trace_metabolites() + >>> layout.set_specific_metabolite('lcts_e', 0.05) # mmol + >>> sim = c.comets(layout, params) + >>> sim.run() + + Attributes + ---------- + + all_params : dict + contains every possible param as keys, with their value as the value + params_units : dict + dictionary containing the units for each possible parameter + + ''' + def show_params(self): + """ + utility function to show all possible parameters and their units + """ + + pdparams = pd.DataFrame({'VALUE': pd.Series(self.all_params), + 'UNITS': pd.Series(self.param_units)}) + return pdparams + + def set_param(self, name : str, value): + """ + alters a specific parameter value + + See show_params() for a list of possible parameters to adjust. + + Parameters + ---------- + name : str + the name of a specific parameter to change + value : variable + the type of the value depends on the parameter. + + Examples + -------- + >>> p = cometspy.params() + >>> p.set_param("writeBiomassLog", True) + >>> p.set_param("maxCycles", 100) + >>> p.set_param("BiomassLogRate",5) + >>> p.set_param("defaultVmax", 10.24) + + """ + if name in self.all_params: + self.all_params[name] = value + else: + print('Parameter ' + name + ' does not exist') + + def get_param(self, name : str) -> object: + """ + returns the value associated with the given parameter name + + Parameters + ---------- + name : str + the name of the parameter whose value is desired + + Returns + ------- + object + the type of object depends on the parameter requested + + """ + if name in self.all_params: + return self.all_params[name] + else: + print('Parameter ' + name + ' does not exist') + + def __init__(self, global_params : str =None, + package_params : str =None): + self.all_params = {'writeSpecificMediaLog': False, + 'specificMediaLogRate': 1, + 'specificMedia': 'ac_e', + 'SpecificMediaLogName': 'specific_media', + 'BiomassLogName': 'biomass', + 'BiomassLogRate': 1, + 'biomassLogFormat': 'COMETS', + 'FluxLogName': 'flux_out', + 'FluxLogRate': 5, + 'fluxLogFormat': 'COMETS', + 'MediaLogName': 'media_out', + 'MediaLogRate': 5, + 'mediaLogFormat': 'COMETS', + 'TotalBiomassLogName': 'total_biomass_out', + 'maxCycles': 100, + 'saveslideshow': False, + 'totalBiomassLogRate': 1, + 'useLogNameTimeStamp': False, + 'writeBiomassLog': False, + 'writeFluxLog': False, + 'writeMediaLog': False, + 'writeTotalBiomassLog': True, + 'batchDilution': False, + 'dilFactor': 10, + 'dilTime': 2, + 'cellSize': 1e-13, + 'allowCellOverlap': True, + 'deathRate': 0, + 'defaultHill': 1, + 'defaultKm': 0.01, + 'defaultVmax': 10, + 'defaultAlpha': 1, + 'defaultW': 10, + 'defaultDiffConst': 1e-5, + 'exchangestyle': 'Monod Style', + 'flowDiffRate': 3e-9, + 'growthDiffRate': 0, + 'maxSpaceBiomass': 0.1, + 'minSpaceBiomass': 0.25e-10, + 'numDiffPerStep': 10, + 'numRunThreads': 1, + 'showCycleCount': True, + 'showCycleTime': False, + 'spaceWidth': 0.02, + 'timeStep': 0.1, + 'toroidalWorld': False, + 'simulateActivation': False, + 'activateRate': 0.001, + 'randomSeed': 0, + 'colorRelative': True, + 'slideshowColorRelative': True, + 'slideshowRate': 1, + 'slideshowLayer': 0, + 'slideshowExt': 'png', + 'biomassMotionStyle': 'Diffusion 2D(Crank-Nicolson)', + 'numExRxnSubsteps': 5, + 'costlyGenome': False, + 'geneFractionalCost': 1e-4, + 'evolution': False, + 'mutRate': 0, + 'addRate': 0, + 'metaboliteDilutionRate': 0.} + self.all_params = dict(sorted(self.all_params.items(), + key=lambda x: x[0])) + + self.param_units = {'writeSpecificMediaLog': "logical", + 'specificMediaLogRate': "cycles", + 'specificMedia': "comma separated string", + 'SpecificMediaLogName': '', + 'BiomassLogName': '', + 'BiomassLogRate': 'cycles', + 'biomassLogFormat': '', + 'FluxLogName': '', + 'FluxLogRate': 'cycles', + 'fluxLogFormat': '', + 'MediaLogName': '', + 'MediaLogRate': 'cycles', + 'mediaLogFormat': '', + 'TotalBiomassLogName': '', + 'maxCycles': 'cycles', + 'saveslideshow': 'logical', + 'totalBiomassLogRate': 'cycles', + 'useLogNameTimeStamp': 'logical', + 'writeBiomassLog': 'logical', + 'writeFluxLog': 'logical', + 'writeMediaLog': 'logical', + 'writeTotalBiomassLog': 'logical', + 'batchDilution': 'logical', + 'dilFactor': '', + 'dilTime': 'hr.', + 'cellSize': 'gr. cell dry weight', + 'allowCellOverlap': 'logical', + 'deathRate': 'percent/cycle', + 'defaultHill': 'unitless', + 'defaultKm': 'M (molar conc.)', + 'defaultVmax': 'mmol /gr. CDW /hr', + 'defaultAlpha': 'slope', + 'defaultW': 'M (molar conc.)', + 'defaultDiffConst': 'cm2/s', + 'exchangestyle': 'uptake type', + 'flowDiffRate': 'cm2/s', + 'growthDiffRate': 'cm2/s', + 'maxSpaceBiomass': 'gr. cell dry weight', + 'minSpaceBiomass': 'gr. cell dry weight', + 'numDiffPerStep': '', + 'numRunThreads': '', + 'showCycleCount': 'logical', + 'showCycleTime': 'logical', + 'spaceWidth': 'cm', + 'timeStep': 'hr.', + 'toroidalWorld': 'logical', + 'simulateActivation': 'logical', + 'activateRate': '', + 'randomSeed': 'numeric', + 'colorRelative': 'logical', + 'slideshowColorRelative': 'logical', + 'slideshowRate': 'cycles', + 'slideshowLayer': 'numeric', + 'slideshowExt': 'extension', + 'biomassMotionStyle': 'type of motion', + 'numExRxnSubsteps': 'numeric', + 'costlyGenome': 'logical', + 'geneFractionalCost': 'percent growth rate^(1/2) /reaction', + 'evolution': 'logical', + 'mutRate': 'deletions /reaction /generation', + 'addRate': 'additions /generation', + 'metaboliteDilutionRate': 'percent/cycle'} + self.param_units = dict(sorted(self.param_units.items(), + key=lambda x: x[0])) + + self.all_type = {'writeSpecificMediaLog': 'global', + 'specificMediaLogRate': 'global', + 'specificMedia': 'global', + 'SpecificMediaLogName': 'global', + 'BiomassLogName': 'global', + 'BiomassLogRate': 'global', + 'biomassLogFormat': 'global', + 'FluxLogName': 'global', + 'FluxLogRate': 'global', + 'fluxLogFormat': 'global', + 'MediaLogName': 'global', + 'MediaLogRate': 'global', + 'mediaLogFormat': 'global', + 'TotalBiomassLogName': 'global', + 'maxCycles': 'package', + 'saveslideshow': 'global', + 'totalBiomassLogRate': 'global', + 'useLogNameTimeStamp': 'global', + 'writeBiomassLog': 'global', + 'writeFluxLog': 'global', + 'writeMediaLog': 'global', + 'writeTotalBiomassLog': 'global', + 'batchDilution': 'global', + 'dilFactor': 'global', + 'dilTime': 'global', + 'cellSize': 'global', + 'allowCellOverlap': 'package', + 'deathRate': 'package', + 'defaultHill': 'package', + 'defaultKm': 'package', + 'defaultVmax': 'package', + 'defaultW': 'package', + 'defaultAlpha': 'package', + 'defaultDiffConst': 'package', + 'exchangestyle': 'package', + 'flowDiffRate': 'package', + 'growthDiffRate': 'package', + 'maxSpaceBiomass': 'package', + 'minSpaceBiomass': 'package', + 'numDiffPerStep': 'package', + 'numRunThreads': 'package', + 'showCycleCount': 'package', + 'showCycleTime': 'package', + 'spaceWidth': 'package', + 'timeStep': 'package', + 'toroidalWorld': 'package', + 'simulateActivation': 'global', + 'activateRate': 'global', + 'randomSeed': 'global', + 'colorRelative': 'global', + 'slideshowColorRelative': 'global', + 'slideshowRate': 'global', + 'slideshowLayer': 'global', + 'slideshowExt': 'global', + 'biomassMotionStyle': 'package', + 'numExRxnSubsteps': 'package', + 'costlyGenome': 'global', + 'geneFractionalCost': 'global', + 'evolution': 'package', + 'mutRate': 'package', + 'addRate': 'package', + 'metaboliteDilutionRate': 'package'} + self.all_type = dict(sorted(self.all_type.items(), + key=lambda x: x[0])) + + # .. parse parameters files to python type variables + if global_params is not None: + with open(global_params) as f: + for line in f: + if '=' in line: + k, v = line.split(' = ') + if v.strip() == 'true': + self.all_params[k.strip()] = True + elif v.strip() == 'false': + self.all_params[k.strip()] = False + elif v.strip().isdigit(): + self.all_params[k.strip()] = int(v.strip()) + elif __isfloat(v.strip()): + self.all_params[k.strip()] = float(v.strip()) + else: + self.all_params[k.strip()] = v.strip() + + if package_params is not None: + with open(package_params) as f: + for line in f: + if '=' in line: + k, v = line.split(' = ') + if v.strip() == 'true': + self.all_params[k.strip()] = True + elif v.strip() == 'false': + self.all_params[k.strip()] = False + elif v.strip().isdigit(): + self.all_params[k.strip()] = int(v.strip()) + elif __isfloat(v.strip()): + self.all_params[k.strip()] = float(v.strip()) + else: + self.all_params[k.strip()] = v.strip() + + # Additional processing. + # If evolution is true, we dont want to write the total biomass log + if self.all_params['evolution']: + self.all_params['writeTotalBiomassLog'] = False + self.all_params['writeBiomassLog'] = True + + def write_params(self, out_glb : str, out_pkg : str): + """ + writes params data to the specified files + + Usually this is only used internally, though a user could use it to + pre-generate params files. + + Parameters + ---------- + out_glb : str + the path and name of the 'global' parameters file to be written + out_pkg : str + the path and name of the 'package' parameters file to be written + + + """ + if os.path.isfile(out_glb): + os.remove(out_glb) + + if os.path.isfile(out_pkg): + os.remove(out_pkg) + + # convert booleans to java format before writing + towrite_params = {} + for k, v in self.all_params.items(): + if v is True: + towrite_params[k] = 'true' + elif v is False: + towrite_params[k] = 'false' + else: + towrite_params[k] = str(v) + + with open(out_glb, 'a') as glb, open(out_pkg, 'a') as pkg: + for k, v in towrite_params.items(): + if self.all_type[k] == 'global': + glb.writelines(k + ' = ' + v + '\n') + else: + pkg.writelines(k + ' = ' + v + '\n') + diff --git a/build/lib/cometspy/utils.py b/build/lib/cometspy/utils.py new file mode 100644 index 0000000..baf59c0 --- /dev/null +++ b/build/lib/cometspy/utils.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" + +The utils module contains helper functions generating spatial patterns. +""" + +from .layout import layout +from .params import params +import random + +def pick_random_locations(n : int, + xrange : tuple, yrange : tuple, + forbidden_locs : set = set()) -> list: + """ + returns a list of n x,y tuples corresponding to locations in the range given + + Parameters + ---------- + + n : int + number of locations desired + xrange : tuple + the x-range (min, max) of the x range possible + yrange : tuple + the y-range (min, max) of the y range possible + forbidden_locs : set, optional + A list of tuples that cannot be chosen. + + Returns + ------- + + list + a list of (x,y) values. + + Examples + -------- + + >>> from cometspy.utils import pick_random_locations + >>> locs = pick_random_locations(3, (0, 10), (0,10)) + >>> locs + + """ + pickable_locs = [(x,y) for x in range(xrange[0], xrange[1]) for y in range(yrange[0], yrange[1])] + pickable_locs = list(set(pickable_locs).difference(set(forbidden_locs))) + if len(pickable_locs) < n: + print("there are fewer available locations than n, returning all available locs") + locs = pickable_locs + else: + locs = random.sample(pickable_locs, n) + + return(locs) + +def grow_rocks(n : int, + xrange : tuple, + yrange : tuple, + mean_size: int) -> list: + """ + grows simple simulated rocks by adding random adjacent points from seeds + + n number of seed points are generated first with pick_random_locations. + Then, mean_size * n - n additional points are added to these seed locations. + For each new point, one random location out of the set of all possible + unoccupied locations next to occupied locations are chosen. Only + lattice points directly to the NSEW are considered. This process is + repeated until all new points are assigned. This usually results in + 'rocks' of different shapes and sizes. + + This function can be very slow (> 1 min) when n * mean_size > 2000 + + Parameters + ---------- + n : int + number of seed points to generate rocks from + xrange : tuple + x range possible, e.g. (0, 5) + yrange : tuple + y range possible, e.g. (0, 10) + mean_size : int + average size in lattice units of a generated rock + + Returns + ------- + list + list of all points generated including seed points + + """ + grow_points = n * mean_size - n + locs = pick_random_locations(n, xrange, yrange) + if (n*mean_size) > ((xrange[1]-xrange[0]) * (yrange[1]-yrange[0])): + print("more rocks asked for than space possible, try again") + return + for i in range(grow_points): + adjacents = _find_unoccupied_adjacent(locs, xrange, yrange) + locs.append(random.sample(adjacents, 1)[0]) + return(locs) + +def _find_unoccupied_adjacent(locs, xrange, yrange): + """ + returns the unoccupied adjacent locations to the given locations + + this is not the set, so unoccupied locations adjacent to > 1 occupied + location will have more spots + """ + locs = [(loc[0], loc[1]) for loc in locs] + adjacent = [] + for loc in locs: + if ((loc[0]-1, loc[1]) not in locs) and ((loc[0]-1) >= xrange[0]): + adjacent.append((loc[0]-1, loc[1])) + if ((loc[0]+1, loc[1]) not in locs) and ((loc[0]+1) <= xrange[1]): + adjacent.append((loc[0]+1, loc[1])) + if ((loc[0], loc[1]-1) not in locs) and ((loc[1]-1) >= yrange[0]): + adjacent.append((loc[0], loc[1]-1)) + if ((loc[0], loc[1]+1) not in locs) and ((loc[1]+1) <= yrange[1]): + adjacent.append((loc[0], loc[1]+1)) + #adjacent = list(set(adjacent)) + return(adjacent) + +def chemostat(models : list, + reservoir_media : dict, + dilution_rate : float) -> tuple: + """ + helper function to let a user skip some steps when generating a chemostat + + This sets relevant simulation parameters (e.g. deathRate, + metaboliteDilutionRate) and layout values (e.g. refresh media) based upon + a "reservoir" definition and a dilution rate. + + It generates a layout that has the reservoir_media as the initial values, + as well as set it to drip in / out based upon the dilution rate. + + The returned layout and params can be further modified before supplying + to a comets object if desired. + + Parameters + ---------- + + models : list(cometspy.model) + list of cometspy.model(s) with initial_pop set to use in the sim + reservoir_media : dict + media definition with metabolite names as keys as mmol amounts as values + dilution_rate : float + the dilution rate of the chemostat, in 1/hr + + Returns + ------- + + tuple (layout, params) + a cometspy.layout object and a cometspy.params object + + Examples + -------- + + >>> import cobra.test + >>> import cometspy as c + >>> from cometspy.utils import chemostat + >>> # make a model from a cobra model, open exchange reactions, and give a pop + >>> tb = cobra.test.create_test_model("textbook") + >>> m = c.model(tb) + >>> m.initial_pop = [0, 0, 1.e-4] + >>> m.open_exchanges() + >>> reservoir = {'glc__D_e' : 0.01, 'nh4_e' : 1000., 'pi_e' : 1000.} + >>> layout, params = chemostat([m], reservoir, 0.1) + >>> params.set_param("maxCycles", 100) + >>> sim = c.comets(layout, params) + >>> sim.run() + >>> print(sim.total_biomass) + + """ + mylayout = layout(models) + for key, value in reservoir_media.items(): + mylayout.set_specific_metabolite(key, value) + mylayout.set_specific_refresh(key, value * dilution_rate) + + parameters = params() + parameters.all_params['metaboliteDilutionRate'] = dilution_rate + parameters.all_params['deathRate'] = dilution_rate + return((mylayout, parameters)) diff --git a/cometspy.egg-info/PKG-INFO b/cometspy.egg-info/PKG-INFO index 6cd7f02..2eec238 100644 --- a/cometspy.egg-info/PKG-INFO +++ b/cometspy.egg-info/PKG-INFO @@ -1,13 +1,14 @@ Metadata-Version: 2.1 Name: cometspy -Version: 0.4.20 +Version: 0.5.0 Summary: The Python interface to COMETS Home-page: https://github.com/segrelab/cometspy -Download-URL: https://github.com/segrelab/cometspy/archive/v0.4.16.tar.gz +Download-URL: https://github.com/segrelab/cometspy/archive/v0.5.0.tar.gz Author: The COMETSPy Core Team Author-email: djordje.bajic@yale.edu License: GPL Keywords: metabolism,dynamic,flux,balance,analysis,spatial,evolution +Platform: UNKNOWN Classifier: Development Status :: 4 - Beta Classifier: Intended Audience :: Science/Research Classifier: License :: OSI Approved :: MIT License @@ -15,6 +16,6 @@ Classifier: Programming Language :: Python :: 3.6 Classifier: Programming Language :: Python :: 3.7 Classifier: Programming Language :: Python :: 3.8 License-File: LICENSE -Requires-Dist: numpy -Requires-Dist: cobra -Requires-Dist: pandas>=1.0.0 + +UNKNOWN + diff --git a/cometspy/__pycache__/comets.cpython-39.pyc b/cometspy/__pycache__/comets.cpython-39.pyc index 0580fac..7adfcb1 100644 Binary files a/cometspy/__pycache__/comets.cpython-39.pyc and b/cometspy/__pycache__/comets.cpython-39.pyc differ diff --git a/cometspy/comets.py b/cometspy/comets.py index e3f4ea4..1eda65e 100644 --- a/cometspy/comets.py +++ b/cometspy/comets.py @@ -21,9 +21,9 @@ __copyright__ = "Copyright 2024, The COMETS Consortium" __credits__ = ["Djordje Bajic", "Jean Vila", "Jeremy Chacon", "Ilija Dukovski"] __license__ = "MIT" -__version__ = "0.5.0" -__comets_compatibility__ = "2.12.0" # version of comets this was tested with (except signaling) -__comets_compatibility_signaling__ = "2.12.0" # version signaling was tested with +__version__ = "0.5.1" +__comets_compatibility__ = "2.12.2" # version of comets this was tested with (except signaling) +__comets_compatibility_signaling__ = "2.12.2" # version signaling was tested with __maintainer__ = "Djordje Bajic" __email__ = "djordje.bajic@yale.edu" diff --git a/cometspy/model.py b/cometspy/model.py index 9e5a41f..185b23f 100644 --- a/cometspy/model.py +++ b/cometspy/model.py @@ -117,6 +117,10 @@ def __init__(self, model : cobra.Model = None, randomtag : bool = False): self.optimizer = 'GUROBI' self.obj_style = 'MAXIMIZE_OBJECTIVE_FLUX' + self.multispecies_convmodel_flag = False + + + if model is not None: if isinstance(model, cobra.Model): self.load_cobra_model(model, randomtag) @@ -339,6 +343,21 @@ def add_convection_parameters(self, packedDensity : float =1., 'elasticModulus': elasticModulus, 'frictionConstant': frictionConstant, 'convDiffConstant': convDiffConstant} + + def add_multspecies_convmodel_parameters(self, pressureKappa : float, pressureExponent : float, packBiomass : float, maxPressure : float): + if not isinstance(pressureKappa, float): + raise ValueError('pressureKappa must be a float') + if not isinstance(pressureExponent, float): + raise ValueError('pressureExponent must be a float') + if not isinstance(packBiomass, float): + raise ValueError('packBiomass must be a float') + if not isinstance(maxPressure, float): + raise ValueError('maxPressure must be a float') + self.multispecies_convmodel_flag = True + self.multimodel_parameters = {'pressureKappa': pressureKappa, + 'pressureExponent': pressureExponent, + 'packBiomass': packBiomass, + 'maxPressure': maxPressure} def add_noise_variance_parameter(self, noiseVariance : float): """ sets the noise variance parameter @@ -1074,6 +1093,11 @@ def write_comets_model(self, working_dir : str = None): for key, value in self.convection_parameters.items(): f.write(key + ' ' + str(value) + '\n') f.write(r'//' + '\n') + + if self.multispecies_convmodel_flag: + for key, value in self.multimodel_parameters.items(): + f.write(key + ' ' + str(value) + '\n') + f.write(r'//' + '\n') if self.nonlinear_diffusion_flag: for key, value in self.nonlinear_diffusion_parameters.items(): diff --git a/dist/cometspy-0.4.11.tar.gz b/dist/cometspy-0.4.11.tar.gz deleted file mode 100644 index 6fe6beb..0000000 Binary files a/dist/cometspy-0.4.11.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.12.tar.gz b/dist/cometspy-0.4.12.tar.gz deleted file mode 100644 index f303669..0000000 Binary files a/dist/cometspy-0.4.12.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.13.tar.gz b/dist/cometspy-0.4.13.tar.gz deleted file mode 100644 index e964d7c..0000000 Binary files a/dist/cometspy-0.4.13.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.14.tar.gz b/dist/cometspy-0.4.14.tar.gz deleted file mode 100644 index a436ac5..0000000 Binary files a/dist/cometspy-0.4.14.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.15.tar.gz b/dist/cometspy-0.4.15.tar.gz deleted file mode 100644 index 2ea254f..0000000 Binary files a/dist/cometspy-0.4.15.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.16.tar.gz b/dist/cometspy-0.4.16.tar.gz deleted file mode 100644 index 7e07394..0000000 Binary files a/dist/cometspy-0.4.16.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.17.tar.gz b/dist/cometspy-0.4.17.tar.gz deleted file mode 100644 index 5584d06..0000000 Binary files a/dist/cometspy-0.4.17.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.18.tar.gz b/dist/cometspy-0.4.18.tar.gz deleted file mode 100644 index eb7f83b..0000000 Binary files a/dist/cometspy-0.4.18.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.19-py3-none-any.whl b/dist/cometspy-0.4.19-py3-none-any.whl deleted file mode 100644 index b1e2a55..0000000 Binary files a/dist/cometspy-0.4.19-py3-none-any.whl and /dev/null differ diff --git a/dist/cometspy-0.4.19.tar.gz b/dist/cometspy-0.4.19.tar.gz deleted file mode 100644 index 4a76632..0000000 Binary files a/dist/cometspy-0.4.19.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.20-py3-none-any.whl b/dist/cometspy-0.4.20-py3-none-any.whl deleted file mode 100644 index 7e2dbe1..0000000 Binary files a/dist/cometspy-0.4.20-py3-none-any.whl and /dev/null differ diff --git a/dist/cometspy-0.4.20.tar.gz b/dist/cometspy-0.4.20.tar.gz deleted file mode 100644 index c25daeb..0000000 Binary files a/dist/cometspy-0.4.20.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.6.tar.gz b/dist/cometspy-0.4.6.tar.gz deleted file mode 100644 index 4fec1a2..0000000 Binary files a/dist/cometspy-0.4.6.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.4.7.tar.gz b/dist/cometspy-0.4.7.tar.gz deleted file mode 100644 index 94cb1d2..0000000 Binary files a/dist/cometspy-0.4.7.tar.gz and /dev/null differ diff --git a/dist/cometspy-0.5.0-py3-none-any.whl b/dist/cometspy-0.5.0-py3-none-any.whl new file mode 100644 index 0000000..29fa6ea Binary files /dev/null and b/dist/cometspy-0.5.0-py3-none-any.whl differ diff --git a/dist/cometspy-0.5.0.tar.gz b/dist/cometspy-0.5.0.tar.gz new file mode 100644 index 0000000..6f92e07 Binary files /dev/null and b/dist/cometspy-0.5.0.tar.gz differ diff --git a/setup.py b/setup.py index 8f9a4d6..e6c20bb 100644 --- a/setup.py +++ b/setup.py @@ -4,13 +4,13 @@ setup( name='cometspy', packages=['cometspy'], - version='0.4.20', + version='0.5.0', license='GPL', description='The Python interface to COMETS', author='The COMETSPy Core Team', author_email='djordje.bajic@yale.edu', url='https://github.com/segrelab/cometspy', - download_url='https://github.com/segrelab/cometspy/archive/v0.4.16.tar.gz', # New releases here!! + download_url='https://github.com/segrelab/cometspy/archive/v0.5.0.tar.gz', # New releases here!! keywords=['metabolism', 'dynamic', 'flux', 'balance', 'analysis', 'spatial', 'evolution'], install_requires=[ # I get to this in a second