re-arrange directories

to ensure package is pip-installable and executable

pymcadre/__init__.py

@@ -0,0 +1,2 @@
+__author__ = "Nantia Leonidou"
+__description__ = " "

pymcadre/check/__init__.py

@@ -0,0 +1,2 @@
+__author__ = "Nantia Leonidou"
+__description__ = " "

pymcadre/check/check_model_function.py

@@ -0,0 +1,115 @@
+__author__ = "Nantia Leonidou"
+__description__ = """ Check production of required metabolites  """
+
+from time import process_time
+from check.find_ex_reactions import *
+from check.set_metabolite_bounds import *
+from check.find_required_rxns import *
+from check.check_rxn_flux import *
+import os
+
+
+def check_model_function(model, *args):
+    """ Check production of required metabolites
+        Inputs: - a cobra model
+                - optional inputs: - 'required_mets', met_lst
+                                   - 'biomass', biomass_rxn
+                                   - 'media', media_def
+
+        Output:  - required_met_status
+                 - time """
+
+    ################################################################################################
+    #                               Parse input parameters                                         #
+    ################################################################################################
+    if len(args) != 0:
+        if len(args) % 2 == 0:  # if division between len(args) and 2 gives rest 0
+
+            options = ["required_mets", "biomass", "media"]
+
+            for i in range(0, len(args), 2):  # i goes from 1 to len(args) incrementally in steps of 2
+                arg_name = args[i]  # define arg_name, can be 'required_mets' or 'biomass' or 'media'
+                arg_val = args[i + 1]  # define arg_val, can be met_lst or biomass_rxn or media_def
+
+                # read input .txt file
+                with open(arg_val, 'r') as file:
+                    # list of all input metabolites
+                    met_lst = [line.rstrip('\n') for line in file]
+
+                # set 1 if one of the options is given as arguments, else set 0
+                option = [1 if arg_name == opt else 0 for opt in options]
+
+                # find indices of nonzero elements in the lst, called options
+                # returns list with indices depending on whether an input argument from *args matches one of the strings in list 'options'
+                # Example: if 'media' is given as an argument in the check_model_function function,
+                #          then the find_idx_option = [2], because the 'media' string is located in the index 2 in 'options' list
+                find_idx_option = [i for i, e in enumerate(option) if e != 0]
+
+                # if indices of nonzero elements are found
+                # biomass_rxn = []
+                # media_def = []
+                if len(option) != 0:
+                    for elem in find_idx_option:
+                        if elem == 0:
+                            arg_val = met_lst
+                        elif elem == 1:
+                            raise Exception("%s option not yet implemented." % arg_name)
+                            # arg_val = biomass_rxn
+                        elif elem == 2:
+                            raise Exception("%s option not yet implemented." % arg_name)
+                            # arg_val = media_def
+                # else-case
+                else:
+                    raise Exception("Unknown option %s." % arg_name)
+        else:
+            raise Exception("Incorrect number of input arguments to function.", os.path.basename(__file__))
+    ################################################################################################
+
+    # Start the stopwatch / counter
+    # returns float value of time in seconds
+    t1_start = process_time()
+
+    # Identify exchange reactions in the model
+    ex_rxns = find_ex_rxns(model)
+
+    # Turn off uptake of organic metabolites
+    if "media_def" in globals():  # check if variable named media_def is defined
+        set_media_ex_bounds(model, ex_rxns)  # not implemented in this version
+    else:
+        set_organic_met_bounds(model, ex_rxns)
+
+    # store reaction names
+    names = []  # rxn name
+    for r in model.reactions:
+        names.append(r.name)
+
+    ###########
+    # Allow uptake of glucose and CO2, only if they are in the model
+    ###########
+    if 'D-Glucose exchange' in names:  # Reaction id: EX_glc__D_e
+        idx = names.index('D-Glucose exchange')
+        model.reactions[idx].lower_bound = -5.0
+
+    if 'CO2 exchange' in names: # Reaction id: EX_co2_e
+        idx = names.index('CO2 exchange')
+        model.reactions[idx].lower_bound = -1000.0
+
+    # Add demand reactions for required metabolites
+    if "met_lst" in globals():
+        model, required_rxns = find_required_rxns(model, met_lst)
+    else:
+        required_rxns = []
+
+    # if "biomass_rxn" in globals():
+    #     required_rxns = np.concatenate(required_rxns, biomass_rxn) # not implemented in this version
+
+    inactive_required = check_rxn_flux(model, required_rxns)
+    required_met_status = int(not len(inactive_required))  # return negation of number corresponding to length of inactive_required list
+
+    # Stop the stopwatch / counter
+    t1_stop = process_time()
+    # compute elapsed time
+    time = t1_stop - t1_start
+
+    print('check_model_function done ....')
+    return required_met_status, time

pymcadre/check/check_rxn_flux.py

@@ -0,0 +1,83 @@
+__author__ = "Nantia Leonidou"
+__description__ = "This function uses the heuristic speed-up proposed by Jerby et al. in the MBA paper for performing a pseudo-FVA calculation. "
+
+# from cobra import *
+# from check.find_required_rxns import *
+import numpy as np
+
+
+def check_rxn_flux(model, required_rxns):
+    """
+        Input:
+            - model: cobra.io.core.model.Model
+            - required_rxns: list of required reactions to check their flux
+
+        Output:
+            - inactive_required
+    """
+    rxn_lst = required_rxns
+    inactive_required = []
+
+    while len(rxn_lst):  # while lst is not empty
+        # set number of reactions equal to number of input reactions
+        num_rxn_lst = len(rxn_lst)
+
+        # modify objective function
+        model.objective = rxn_lst
+        # optimize using glpk solver
+        fba_solution = model.optimize(objective_sense="maximize")
+        # get flux values of all reactions as an array
+        opt_max = fba_solution.fluxes.values
+
+        # If no solution was achieved while trying to maximize all reactions, skip the next
+        # step of checking individual reactions
+        if len(opt_max) == 0:
+            inactive_required.append(1)
+            break  # break the whole while-loop
+
+        # check whether reactions from model are included in input required reactions
+        # put 1 if reaction from model is in the list of required reactions
+        is_member = [1 if react.id in required_rxns else 0 for react in model.reactions]
+
+        # get fluxes for these specific reactions
+        required_flux = [opt_max[idx] for idx in range(len(is_member)) if is_member[idx]]
+
+        # absolute values of all fluxes
+        abs_fluxes = [abs(ele) for ele in required_flux]
+        # filter out those with a specific value --> those reactions will be kept
+        filtered = [1 if abs_fluxes[idx] >= 1e-6 else 0 for idx in range(len(abs_fluxes))]
+        # get related reaction IDs
+        active_required = [required_rxns[idx] for idx in range(len(filtered)) if filtered[idx]]
+
+        # update rxn_lst to be the difference between reactions' list and required active reactions
+        rxn_lst = (list(set(rxn_lst) - set(active_required)))
+        # get how many reactions were removed from required_rxns list --> updated length of rxns_lst
+        num_removed = num_rxn_lst - len(rxn_lst)
+
+        if not num_removed:  # if this if-case is fulfilled, throws an out-of-bounds error, similar happens to MATLAB code as well
+            # create a random series of numbers from 0 to len(rxn_lst)
+            rand_ind = np.random.permutation(len(rxn_lst))
+
+            # extract reaction from rxn_lst based on the position taken from the first number in the randomly generated series
+            i = rxn_lst[rand_ind[0]]
+
+            # change objective function
+            model.objective = i
+            # Maximize reaction i
+            fba_solution = model.optimize(objective_sense="maximize")
+            # get flux values of all reactions
+            opt_max = fba_solution.fluxes
+            if len(opt_max) == 0:
+                inactive_required.append(i)
+                break  # break the whole while-loop
+
+            abs_opt = [abs(ele) for ele in opt_max]
+            for n in abs_opt:
+                if n < 1e-6:
+                    inactive_required.append(i)
+                    break
+
+            rxn_lst.remove(i)
+
+    return inactive_required
+

pymcadre/check/find_ex_reactions.py

@@ -0,0 +1,26 @@
+__author__ = "Nantia Leonidou"
+__description__ = " Identification of all exchange reactions in a given Cobra model"
+
+
+from cobra import *
+
+
+def find_ex_rxns(model):
+
+    """ This function identifies reactions that exchange metabolites into and out of the extra- and intracellular space
+        (sinks, demands and exchanges)
+
+        Input:
+            - model: cobra.io.core.model.Model
+
+        Ouput:
+            - ex_rxns_ids: list of all identified reactions' IDs """
+
+    ex_rxns_ids = []
+
+    for react in model.reactions:
+
+        if len(react.metabolites) == 1:
+            ex_rxns_ids.append(react.id)
+
+    return ex_rxns_ids

pymcadre/check/find_organic_ex_rxns.py

@@ -0,0 +1,45 @@
+__author__ = "Nantia Leonidou"
+__description__ = """ Detection of all organic reactions given a cobra model and a list of exchange reactions 
+                   (sink, demand and exchange reactions) """
+
+
+# from cobra import *
+# from cobra.io.sbml import *
+# from check.find_ex_reactions import *
+
+def find_organic_ex_rxns(model, ex_rxns):
+    """ Detection of organic reactions included in all exchange reactions
+
+        Input:
+            - model: cobra.io.core.model.Model
+            - ex_rxns: list of reaction IDs
+
+        Output:
+            - organic_ex_rxns: list of organic reactions' IDs"""
+
+    # Organic metabolites contain carbon (C) and hydrogen (H) Thus, check molecular formulae
+    organic_metabolites = []
+    for met in model.metabolites:
+        if "C" in met.formula and "H" in met.formula:
+            organic_metabolites.append(met.id)
+        else:
+            # add metabolites, that are organic but don't contain only "C" and "H" in their chemical formulae OR include protein compounds
+            if ("C" in met.formula and "R" in met.formula) or ("X" in met.formula and "H" in met.formula) \
+                    or ("C" in met.formula and "X" in met.formula) or ("H" in met.formula and "R" in met.formula):
+                organic_metabolites.append(met.id)
+
+    # Find all organic reactions (i.e reactions evolving organic metabolites)
+    all_organic_reactions = []
+    for met in organic_metabolites:
+        # all organic reactions associated with current organic metabolite
+        org_rxns = model.metabolites.get_by_id(met).reactions
+        for react in org_rxns:
+            # put all of them in a list, after excluding reactions that may come up multiple times
+            if react.id not in all_organic_reactions:
+                all_organic_reactions.append(react.id)
+
+    # extract all exchange reactions that involve organic metabolites
+    # organic_ex_rxns = [react for react in ex_rxns if react in all_organic_reactions]
+    organic_ex_rxns = list(set(all_organic_reactions) & set(ex_rxns))
+
+    return organic_ex_rxns

pymcadre/check/find_required_rxns.py

@@ -0,0 +1,70 @@
+__author__ = "Nantia Leonidou"
+__description__ = " Detect required reactions and add them to the model"
+
+from cobra import *
+
+
+def add_metabolites(met_id_lst, coeffs_lst, compartment_lst):
+    """ Create a dictionary with metabolites and their coefficients
+        Input: lst with metabolites' IDs, list with repsective coefficients
+               and lst with respective compartment letters
+        Output: dictionary with metabolites"""
+
+    if len(met_id_lst) != len(coeffs_lst) or len(met_id_lst) != len(compartment_lst) or len(coeffs_lst) != len(
+            coeffs_lst):
+        print("ERROR: Input lists should have the same size")
+    else:
+        met_dict = {}
+        for i in range(len(met_id_lst)):
+            met_dict[Metabolite(met_id_lst[i], compartment=compartment_lst[i])] = coeffs_lst[i]
+        return met_dict
+
+
+def add_single_reaction(model, id, lb, up, met_dict):
+    """ Add reactions to a model.
+        Input: model, reaction id, lower/upper bound and dictionary with the metabolites to add
+        Output: model with a newly added reaction and the respective reaction."""
+
+    react = Reaction()
+    react.id = id
+    react.lower_bound, react.upper_bound = lb, up
+    # add metabolites by hand in the model, First checks whether the metabolites exist in the model
+    react.add_metabolites(met_dict)
+
+    return model, react
+
+
+def find_required_rxns(model, metabolites):
+    """ Input: model and list of metabollites in a .txt format (see example metabolites_humanModel.txt)
+        Output: model and list of only required reactions' id """
+
+    # define manually _coa metabolites
+    coa_mets = ['accoa_m', 'succoa_m', 'pmtcoa_c']
+
+    # read input file including diverse metabolites
+    with open(metabolites, 'r') as file:
+        # list of all input metabolites
+        met_lst = [line.rstrip('\n') for line in file]
+
+    # difference between coa_mets and input list of metabolites --> list of strings
+    diff = (list(set(met_lst) - set(coa_mets)))
+    print(diff)
+
+    required_rxns = []
+    # add metabolites from diff as demand reactions, if metabolite is in the model
+    model_mets = [m.id for m in model.metabolites]
+    for met in diff:
+        r = model.add_boundary(model.metabolites.get_by_id(met),type="demand")  # add_boundary if metabolite is already in the model
+        required_rxns.append(r.id)
+
+
+    react1 = add_single_reaction(model, "DM_accoa_m",  0.0,1000.0, met_dict={Metabolite('accoa_m', compartment='m'): -1, Metabolite('coa_m', compartment='m'): 1})[1]
+    required_rxns.append(react1.id)
+    react2 = add_single_reaction(model, "DM_succoa_m", 0.0,1000.0, met_dict={Metabolite('succoa_m', compartment='m'): -1, Metabolite('coa_m', compartment='m'): 1})[1]
+    required_rxns.append(react2.id)
+    react3 = add_single_reaction(model, "DM_pmtcoa_c", 0.0,1000, met_dict={Metabolite('pmtcoa_c', compartment='c'): -1, Metabolite('coa_c', compartment='c'): 1})[1]
+    required_rxns.append(react3.id)
+    model.add_reactions([react1,react2,react3])
+
+    return model, required_rxns
+

pymcadre/check/set_metabolite_bounds.py

@@ -0,0 +1,40 @@
+__author__ = "Nantia Leonidou"
+__description__ = " Modify metabolites' bounds "
+
+from check.find_organic_ex_rxns import *
+
+
+def set_organic_met_bounds(model, ex_rxns):
+    """ Set lower bound of organic reactions to 0.0, aiming the turn-off of the uptake
+
+        Input: a cobra model and a list of reaction IDs
+
+        Output: a cobra model with the updated bounds  """
+
+    # Identify all exchange reactions that include organic metabolites
+    organic_ex_rxns = find_organic_ex_rxns(model, ex_rxns)
+
+    # Reset all lower bounds for organic reactions to 0 to disable uptake
+    for react in organic_ex_rxns:
+        model.reactions.get_by_id(react).lower_bound = 0.0
+
+    #############
+    # check whether bounds were set to 0.0
+    # for r in organic_ex_rxns:
+    #     print(model.reactions.get_by_id(r).bounds)
+    ############
+
+    return model
+
+
+##############################################################################
+# not implemented in this version --> see MATLAB code
+def set_media_ex_bounds(model, ex_rxns):
+    # for rxns in ex_rxns:
+    #
+    #     # set all lower and upper bounds to 0
+    #     model.reactions.get_by_id(rxns).lower_bound = 0.0
+    #     model.reactions.get_by_id(rxns).upper_bound = 0.0
+    print("Not implemented yet.")
+    return model
+###############################################################################