automatic detection of external compartment

scripts/smetana

@@ -59,7 +59,6 @@ if __name__ == '__main__':
     parser.add_argument('-v', '--verbose', action='store_true', dest='verbose', help="Switch to verbose mode")
     parser.add_argument('-z', '--zeros', action='store_true', dest='zeros', help="Include entries with zero score.")
     parser.add_argument('--solver', help="Change default solver (current options: 'gurobi', 'cplex').")
-    parser.add_argument('--ext', help="Extracellular compartment identifier in the models (default 'C_e').")
     parser.add_argument('--molweight', action='store_true', help="Use molecular weight minimization (recomended).")
     parser.add_argument('--lp', action='store_true', help=argparse.SUPPRESS)
     parser.add_argument('--exclude', help="List of compounds to exclude from calculations (e.g.: inorganic compounds).")
@@ -87,8 +86,6 @@ if __name__ == '__main__':
     if args.solver:
         set_default_solver(args.solver)
 
-    ext_comp_id = args.ext if args.ext else 'C_e'
-
     main(
         models=args.models,
         communities=args.communities,
@@ -98,7 +95,6 @@ if __name__ == '__main__':
         media=args.media,
         mediadb=args.mediadb,
         zeros=args.zeros,
-        ext_comp_id=ext_comp_id,
         verbose=args.verbose,
         min_mol_weight=args.molweight,
         use_lp=args.lp,

setup.py

@@ -41,6 +41,6 @@
     packages=find_packages(include=['smetana']),
     setup_requires=['setuptools_scm'],
     url='https://github.com/cdanielmachado/smetana',
-    version='1.1.0',
+    version='1.2.0',
     zip_safe=False,
 )

smetana/__init__.py

@@ -4,5 +4,5 @@
 
 __author__ = "Daniel Machado"
 __email__ = '[email protected]'
-__version__ = '1.1.0'
+__version__ = '1.2.0'
 

smetana/interface.py

@@ -313,7 +313,7 @@ def export_results(mode, output, data, debug_data, zeros):
 
 
 def main(models, communities=None, mode=None, output=None, flavor=None, media=None, mediadb=None, zeros=False,
-         ext_comp_id=None, verbose=False, min_mol_weight=False, use_lp=False, exclude=None, debug=False,
+         verbose=False, min_mol_weight=False, use_lp=False, exclude=None, debug=False,
          other=None, n=1, p=1, ignore_coupling=False):
 
     other_models = other if mode == "biotic" else None
@@ -331,7 +331,7 @@ def main(models, communities=None, mode=None, output=None, flavor=None, media=No
             print("Loading community: " + comm_id)
 
         comm_models = [model_cache.get_model(org_id, reset_id=True) for org_id in organisms]
-        community = Community(comm_id, comm_models, copy_models=False, extracellular_compartment_id=ext_comp_id)
+        community = Community(comm_id, comm_models, copy_models=False)
 
         for medium in media:
 

smetana/legacy.py

@@ -28,12 +28,13 @@ def __init__(self, original_reaction=None, organism_reaction=None, original_meta
         self.community_exchange_reaction = community_exchange_reaction
 
     def __repr__(self):
-        return "<orig_m: {}, org_m: {}, ex_m: {}, orig_r: {}, org_r: {}, exch_r: {}>".format(self.original_metabolite,
-                                                                                             self.organism_metabolite,
-                                                                                             self.extracellular_metabolite,
-                                                                                             self.original_reaction,
-                                                                                             self.organism_reaction,
-                                                                                             self.community_exchange_reaction)
+        repr_str =  "<orig_m: {}, org_m: {}, ex_m: {}, orig_r: {}, org_r: {}, exch_r: {}>"
+        return repr_str.format(self.original_metabolite,
+                               self.organism_metabolite,
+                               self.extracellular_metabolite,
+                               self.original_reaction,
+                               self.organism_reaction,
+                               self.community_exchange_reaction)
 
 
 class Community(object):
@@ -44,7 +45,7 @@ class Community(object):
     into a single multi-species model (CBModel) which is compatible with most types of constraint-based methods.
     """
 
-    def __init__(self, community_id, models=None, copy_models=True, extracellular_compartment_id="C_e",
+    def __init__(self, community_id, models=None, copy_models=True,
                  merge_extracellular_compartments=False, create_biomass=True, interacting=True,
                  exchanged_metabolites_blacklist=set()):
         """
@@ -66,7 +67,6 @@ def __init__(self, community_id, models=None, copy_models=True, extracellular_co
 
         self.id = community_id
         self._organisms = AttrOrderedDict()
-        self._extracellular_compartment = extracellular_compartment_id  # TODO: maybe merge and compartment id arguments should be merged?
         self._merge_extracellular_compartments = merge_extracellular_compartments
         self._create_biomass = create_biomass
         self._merged_model = None
@@ -268,12 +268,6 @@ def _copy_object(obj, org_id, compartment=None):
 
             return new_obj
 
-        models_missing_extracelullar_compartment = [m.id for m in self._organisms.values()
-                                                    if self._extracellular_compartment not in m.compartments]
-        if models_missing_extracelullar_compartment:
-            raise RuntimeError("Extracellular compartment '{}' missing from models: '{}'".format(
-                self._extracellular_compartment, "', '".join(models_missing_extracelullar_compartment)))
-
         models_missing_biomass = [m.id for m in self._organisms.values() if not m.biomass_reaction]
         if models_missing_biomass:
             raise RuntimeError("Biomass reaction not found in models: {}".format("', '".join(models_missing_biomass)))
@@ -299,14 +293,14 @@ def _copy_object(obj, org_id, compartment=None):
                 merged_model.add_compartment(export_pool_compartment)
 
             for c_id, comp in model.compartments.items():
-                if c_id != self._extracellular_compartment or not self._merge_extracellular_compartments:
+                if not comp.external or not self._merge_extracellular_compartments:
                     new_comp = _copy_object(comp, org_id)
                     merged_model.add_compartment(new_comp)
                 elif c_id not in merged_model.compartments:
                     merged_model.add_compartment(deepcopy(comp))
 
             for m_id, met in model.metabolites.items():
-                if met.compartment != self._extracellular_compartment or not self._merge_extracellular_compartments:
+                if not model.compartments[met.compartment].external or not self._merge_extracellular_compartments:
                     new_met = _copy_object(met, org_id, _id_pattern(met.compartment, org_id))
                     merged_model.add_metabolite(new_met)
                 elif m_id not in merged_model.metabolites:
@@ -354,15 +348,15 @@ def _copy_object(obj, org_id, compartment=None):
 
                     for m_id, coeff in rxn.stoichiometry.items():
                         m_blacklisted = m_id in self._exchanged_metabolites_blacklist
-                        if (model.metabolites[m_id].compartment != self._extracellular_compartment
+                        if (not model.compartments[model.metabolites[m_id].compartment].external
                                 or not self._merge_extracellular_compartments):
                             del new_rxn.stoichiometry[m_id]
                             new_id = _id_pattern(m_id, org_id)
                             new_rxn.stoichiometry[new_id] = coeff
 
                         if is_exchange:
                             new_rxn.reaction_type = ReactionType.OTHER
-                            if (model.metabolites[m_id].compartment == self._extracellular_compartment
+                            if (model.compartments[model.metabolites[m_id].compartment].external
                                     and not self._merge_extracellular_compartments):
                                 # TODO: if m_id in self._exchanged_metabolites_blacklist:
                                 pool_id = _id_pattern(m_id, "pool")
@@ -404,8 +398,8 @@ def _copy_object(obj, org_id, compartment=None):
 
                         # Add biomass metabolite to biomass equation
                         m_id = _id_pattern('Biomass', org_id)
-                        name = _name_pattern('Framed biomass', org_id)
-                        comp = 'pool' if not self._merge_extracellular_compartments else self._extracellular_compartment
+                        name = _name_pattern('Community biomass', org_id)
+                        comp = 'pool'
                         biomass_met = Metabolite(m_id, name, comp)
                         merged_model.add_metabolite(biomass_met)
                         new_rxn.stoichiometry[m_id] = 1
@@ -435,7 +429,7 @@ def _copy_object(obj, org_id, compartment=None):
 
                         m_id = _id_pattern('Biomass', org_id)
                         name = _name_pattern('Biomass', org_id)
-                        comp = 'pool' if not self._merge_extracellular_compartments else self._extracellular_compartment
+                        comp = 'pool'
                         biomass_met = Metabolite(m_id, name, comp)
                         merged_model.add_metabolite(biomass_met)
                         new_rxn.stoichiometry[m_id] = 1
@@ -485,7 +479,6 @@ def copy(self, merge_extracellular_compartments=None, copy_models=None, interact
 
         copy_community = Community(self.id, models=list(self._organisms.values()),
                                    copy_models=copy_models, create_biomass=create_biomass,
-                                   extracellular_compartment_id=self._extracellular_compartment,
                                    merge_extracellular_compartments=merge_extracellular_compartments,
                                    interacting=interacting,
                                    exchanged_metabolites_blacklist=exchanged_metabolites_blacklist)

tests/output/test_detailed.tsv

@@ -1,12 +1,12 @@
 community	medium	receiver	donor	compound	scs	mus	mps	smetana
-all	M9	ec_glc_ko	ec_nh4_ko	M_ac_e	1.0	0.11	1	0.11
-all	M9	ec_glc_ko	ec_nh4_ko	M_acald_e	1.0	0.28	1	0.28
-all	M9	ec_glc_ko	ec_nh4_ko	M_akg_e	1.0	0.25	1	0.25
-all	M9	ec_glc_ko	ec_nh4_ko	M_co2_e	1.0	0.02	1	0.02
-all	M9	ec_glc_ko	ec_nh4_ko	M_etoh_e	1.0	0.31	1	0.31
-all	M9	ec_glc_ko	ec_nh4_ko	M_lac__D_e	1.0	0.37	1	0.37
-all	M9	ec_glc_ko	ec_nh4_ko	M_pyr_e	1.0	0.24	1	0.24
-all	M9	ec_glc_ko	ec_nh4_ko	M_succ_e	1.0	0.09	1	0.09
+all	M9	ec_glc_ko	ec_nh4_ko	M_ac_e	1.0	0.11627906976744186	1	0.11627906976744186
+all	M9	ec_glc_ko	ec_nh4_ko	M_acald_e	1.0	0.2558139534883721	1	0.2558139534883721
+all	M9	ec_glc_ko	ec_nh4_ko	M_akg_e	1.0	0.23255813953488372	1	0.23255813953488372
+all	M9	ec_glc_ko	ec_nh4_ko	M_etoh_e	1.0	0.23255813953488372	1	0.23255813953488372
+all	M9	ec_glc_ko	ec_nh4_ko	M_lac__D_e	1.0	0.23255813953488372	1	0.23255813953488372
+all	M9	ec_glc_ko	ec_nh4_ko	M_pyr_e	1.0	0.2558139534883721	1	0.2558139534883721
+all	M9	ec_glc_ko	ec_nh4_ko	M_succ_e	1.0	0.046511627906976744	1	0.046511627906976744
+all	M9	ec_nh4_ko	ec_glc_ko	M_acald_e	1.0	0.07692307692307693	1	0.07692307692307693
 all	M9	ec_nh4_ko	ec_glc_ko	M_glu__L_e	1.0	1.0	1	1.0
-all	M9	ec_nh4_ko	ec_glc_ko	M_lac__D_e	1.0	0.08333333333333333	1	0.08333333333333333
-all	M9	ec_nh4_ko	ec_glc_ko	M_pyr_e	1.0	0.08333333333333333	1	0.08333333333333333
+all	M9	ec_nh4_ko	ec_glc_ko	M_lac__D_e	1.0	0.07692307692307693	1	0.07692307692307693
+all	M9	ec_nh4_ko	ec_glc_ko	M_pyr_e	1.0	0.07692307692307693	1	0.07692307692307693

tests/test_smetana.py

@@ -9,7 +9,7 @@
 class TestGlobal(unittest.TestCase):
     def test_global(self):
         main(["tests/data/ec_*_ko.xml"], mode="global", output="tests/output/test", media="M9,LB",
-             mediadb="tests/data/media_db.tsv", ext_comp_id="e", exclude="tests/data/inorganic.txt")
+             mediadb="tests/data/media_db.tsv", exclude="tests/data/inorganic.txt")
         df = pd.read_csv("tests/output/test_global.tsv")
         self.assertEqual(df.shape[0], 2)
 
@@ -18,7 +18,7 @@ class TestDetailed(unittest.TestCase):
 
     def test_detailed(self):
         main(["tests/data/ec_*.xml"], mode="detailed", output="tests/output/test", media="M9,LB",
-             mediadb="tests/data/media_db.tsv", ext_comp_id="e", exclude="tests/data/inorganic.txt")
+             mediadb="tests/data/media_db.tsv", exclude="tests/data/inorganic.txt")
         df = pd.read_csv("tests/output/test_detailed.tsv")
         self.assertGreater(df.shape[0], 5)
         self.assertLess(df.shape[0], 15)