Skip to content

Commit

Permalink
Merge pull request #199 from jonls/gapfill-compartments
Browse files Browse the repository at this point in the history 
Use compartment information when gap-filling
  • Loading branch information
@jonls
jonls committed Mar 1, 2017
2 parents 2e75025 + 4b4e788 commit ee2e9a6
Show file tree
Hide file tree
Showing 5 changed files with 181 additions and 100 deletions.
6 changes: 5 additions & 1 deletion psamm/datasource/context.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,11 @@ def __init__(self, arg):
self._filepath = arg
else:
self._filepath = arg.filepath
self._basepath = os.path.dirname(self._filepath)

if self._filepath is not None:
self._basepath = os.path.dirname(self._filepath)
else:
self._basepath = None
else:
self._filepath = None
self._basepath = None
Expand Down
46 changes: 37 additions & 9 deletions psamm/fastgapfill.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,16 +52,44 @@ def create_extended_model(model, db_penalty=None, ex_penalty=None,

# Create metabolic model
model_extended = model.create_metabolic_model()
model_compartments = set(model_extended.compartments)
extra_compartments = model.extracellular_compartment

# Add exchange and transport reactions to database
logger.info('Adding database, exchange and transport reactions')
db_added = model_extended.add_all_database_reactions(model_compartments)
extra_compartment = model.extracellular_compartment
compartments, boundaries = model.parse_compartments()

compartment_ids = set(c.id for c in compartments)

# Add database reactions to extended model
if len(compartment_ids) > 0:
logger.info(
'Using all database reactions in compartments: {}...'.format(
', '.join('{}'.format(c) for c in compartment_ids)))
db_added = model_extended.add_all_database_reactions(compartment_ids)
else:
logger.warning(
'No compartments specified in the model; database reactions will'
' not be used! Add compartment specification to model to include'
' database reactions for those compartments.')
db_added = set()

# Add exchange reactions to extended model
logger.info(
'Using artificial exchange reactions for compartment: {}...'.format(
extra_compartment))
ex_added = model_extended.add_all_exchange_reactions(
extra_compartments, allow_duplicates=True)
tp_added = model_extended.add_all_transport_reactions(
extra_compartments, allow_duplicates=True)
extra_compartment, allow_duplicates=True)

# Add transport reactions to extended model
if len(boundaries) > 0:
logger.info(
'Using artificial transport reactions for the compartment'
' boundaries: {}...'.format(
'; '.join('{}<->{}'.format(c1, c2) for c1, c2 in boundaries)))
tp_added = model_extended.add_all_transport_reactions(
boundaries, allow_duplicates=True)
else:
logger.warning(
'No compartment boundaries specified in the model;'
' artificial transport reactions will not be used!')
tp_added = set()

# Add penalty weights on reactions
weights = {}
Expand Down
75 changes: 50 additions & 25 deletions psamm/metabolicmodel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -282,12 +282,15 @@ def add_all_exchange_reactions(self, compartment, allow_duplicates=False):
all_reactions[rx] = rxnid

added = set()
for compound in sorted(self.compounds):
initial_compounds = set(self.compounds)
for model_compound in initial_compounds:
compound = model_compound.in_compartment(compartment)
rxnid_ex = ('rxnex', compound)
if rxnid_ex in added:
continue

if not self._database.has_reaction(rxnid_ex):
reaction_ex = Reaction(Direction.Both, {
compound.in_compartment(compartment): -1
})
reaction_ex = Reaction(Direction.Both, {compound: -1})
if reaction_ex not in all_reactions:
self._database.set_reaction(rxnid_ex, reaction_ex)
else:
Expand All @@ -299,8 +302,20 @@ def add_all_exchange_reactions(self, compartment, allow_duplicates=False):

return added

def add_all_transport_reactions(self, compartment, allow_duplicates=False):
"""Add all transport reactions to database and to model"""
def add_all_transport_reactions(self, boundaries, allow_duplicates=False):
"""Add all transport reactions to database and to model.
Add transport reactions for all boundaries. Boundaries are defined
by pairs (2-tuples) of compartment IDs. Transport reactions are
added for all compounds in the model, not just for compounds in the
two boundary compartments.
Args:
boundaries: Set of compartment boundary pairs.
Returns:
Set of IDs of reactions that were added.
"""

all_reactions = {}
if not allow_duplicates:
Expand All @@ -310,26 +325,36 @@ def add_all_transport_reactions(self, compartment, allow_duplicates=False):
rx = self._database.get_reaction(rxnid)
all_reactions[rx] = rxnid

added = set()
for compound in sorted(self.compounds):
if compound.compartment == compartment:
# A transport reaction with exchange would not be valid
continue
boundary_pairs = set()
for source, dest in boundaries:
if source != dest:
boundary_pairs.add(tuple(sorted((source, dest))))

rxnid_tp = ('rxntp', compound)
if not self._database.has_reaction(rxnid_tp):
reaction_tp = Reaction(Direction.Both, {
compound.in_compartment(compartment): -1,
compound: 1
})
if reaction_tp not in all_reactions:
self._database.set_reaction(rxnid_tp, reaction_tp)
else:
rxnid_tp = all_reactions[reaction_tp]

if rxnid_tp not in self._reaction_set:
added.add(rxnid_tp)
self.add_reaction(rxnid_tp)
added = set()
initial_compounds = set(self.compounds)
for compound in initial_compounds:
for c1, c2 in boundary_pairs:
compound1 = compound.in_compartment(c1)
compound2 = compound.in_compartment(c2)
pair = compound1, compound2

rxnid_tp = ('rxntp',) + pair
if rxnid_tp in added:
continue

if not self._database.has_reaction(rxnid_tp):
reaction_tp = Reaction(Direction.Both, {
compound1: -1,
compound2: 1
})
if reaction_tp not in all_reactions:
self._database.set_reaction(rxnid_tp, reaction_tp)
else:
rxnid_tp = all_reactions[reaction_tp]

if rxnid_tp not in self._reaction_set:
added.add(rxnid_tp)
self.add_reaction(rxnid_tp)

return added

Expand Down
104 changes: 62 additions & 42 deletions psamm/tests/test_fastgapfill.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,40 +32,56 @@

class TestCreateExtendedModel(unittest.TestCase):
def setUp(self):
self._model_dir = tempfile.mkdtemp()
with open(os.path.join(self._model_dir, 'model.yaml'), 'w') as f:
f.write('\n'.join([
'---',
'reactions:',
' - id: rxn_1',
' equation: A[e] <=> B[c]',
' - id: rxn_2',
' equation: A[e] => C[c]',
' - id: rxn_3',
' equation: A[e] => D[e]',
' - id: rxn_4',
' equation: C[c] => D[e]',
'compounds:',
' - id: A',
' - id: B',
' - id: C',
' - id: D',
'media:',
' - compartment: e',
' compounds:',
' - id: A',
' reaction: rxn_5',
' - id: D',
' reaction: rxn_6',
'model:',
' - reactions:',
' - rxn_1',
' - rxn_2',
]))
self._model = NativeModel.load_model_from_path(self._model_dir)

def tearDown(self):
shutil.rmtree(self._model_dir)
self._model = NativeModel({
'compartments': [
{
'id': 'c',
'adjacent_to': 'e'
}, {
'id': 'e'
}
],
'reactions': [
{
'id': 'rxn_1',
'equation': 'A[e] <=> B[c]'
}, {
'id': 'rxn_2',
'equation': 'A[e] => C[c]'
}, {
'id': 'rxn_3',
'equation': 'A[e] => D[e]'
}, {
'id': 'rxn_4',
'equation': 'C[c] => D[e]'
}
],
'compounds': [
{'id': 'A'},
{'id': 'B'},
{'id': 'C'},
{'id': 'D'}
],
'media': [{
'compartment': 'e',
'compounds': [
{
'id': 'A',
'reaction': 'rxn_5'
},
{
'id': 'D',
'reaction': 'rxn_6'
}
]
}],
'model': [{
'reactions': [
'rxn_1',
'rxn_2'
]
}]
})

def test_create_model_extended(self):
expected_reactions = set([
Expand All @@ -75,22 +91,26 @@ def test_create_model_extended(self):
'rxn_4',
'rxn_5',
'rxn_6',
('rxntp', Compound('B', 'c')),
('rxntp', Compound('C', 'c')),
('rxntp', Compound('A', 'c'), Compound('A', 'e')),
('rxntp', Compound('B', 'c'), Compound('B', 'e')),
('rxntp', Compound('C', 'c'), Compound('C', 'e')),
('rxntp', Compound('D', 'c'), Compound('D', 'e')),
('rxnex', Compound('A', 'e')),
('rxnex', Compound('B', 'c')),
('rxnex', Compound('C', 'c')),
('rxnex', Compound('B', 'e')),
('rxnex', Compound('C', 'e')),
('rxnex', Compound('D', 'e'))
])

expected_weights = {
'rxn_3': 5.6,
'rxn_4': 1.0,
('rxntp', Compound('B', 'c')): 3.0,
('rxntp', Compound('C', 'c')): 3.0,
('rxntp', Compound('A', 'c'), Compound('A', 'e')): 3.0,
('rxntp', Compound('B', 'c'), Compound('B', 'e')): 3.0,
('rxntp', Compound('C', 'c'), Compound('C', 'e')): 3.0,
('rxntp', Compound('D', 'c'), Compound('D', 'e')): 3.0,
('rxnex', Compound('A', 'e')): 2.0,
('rxnex', Compound('B', 'c')): 2.0,
('rxnex', Compound('C', 'c')): 2.0,
('rxnex', Compound('B', 'e')): 2.0,
('rxnex', Compound('C', 'e')): 2.0,
('rxnex', Compound('D', 'e')): 2.0
}
penalties = {'rxn_3': 5.6}
Expand Down
50 changes: 27 additions & 23 deletions psamm/tests/test_metabolicmodel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,12 +47,12 @@ def test_load_model_with_reaction_subset(self):
class TestMetabolicModel(unittest.TestCase):
def setUp(self):
self.database = DictDatabase()
self.database.set_reaction('rxn_1', parse_reaction('=> (2) |A|'))
self.database.set_reaction('rxn_2', parse_reaction('|A| <=> |B|'))
self.database.set_reaction('rxn_3', parse_reaction('|A| => |D[e]|'))
self.database.set_reaction('rxn_4', parse_reaction('|A| => |C|'))
self.database.set_reaction('rxn_5', parse_reaction('|C| => |D[e]|'))
self.database.set_reaction('rxn_6', parse_reaction('|D[e]| =>'))
self.database.set_reaction('rxn_1', parse_reaction('=> (2) A[c]'))
self.database.set_reaction('rxn_2', parse_reaction('A[c] <=> B[c]'))
self.database.set_reaction('rxn_3', parse_reaction('A[c] => D[e]'))
self.database.set_reaction('rxn_4', parse_reaction('A[c] => C[c]'))
self.database.set_reaction('rxn_5', parse_reaction('C[c] => D[e]'))
self.database.set_reaction('rxn_6', parse_reaction('D[e] =>'))
self.model = MetabolicModel.load_model(
self.database, self.database.reactions)

Expand All @@ -66,14 +66,14 @@ def test_reaction_set(self):

def test_compound_set(self):
self.assertEqual(set(self.model.compounds),
{Compound('A'), Compound('B'),
Compound('C'), Compound('D', 'e')})
{Compound('A', 'c'), Compound('B', 'c'),
Compound('C', 'c'), Compound('D', 'e')})

def test_compartments(self):
self.assertEqual(set(self.model.compartments), {None, 'e'})
self.assertEqual(set(self.model.compartments), {'c', 'e'})

def test_add_reaction_new(self):
self.database.set_reaction('rxn_7', parse_reaction('|D[e]| => |E[e]|'))
self.database.set_reaction('rxn_7', parse_reaction('D[e] => E[e]'))
self.model.add_reaction('rxn_7')
self.assertIn('rxn_7', set(self.model.reactions))
self.assertIn(Compound('E', 'e'), set(self.model.compounds))
Expand All @@ -83,9 +83,9 @@ def test_add_reaction_existing(self):
self.assertEqual(
set(self.model.reactions),
{'rxn_1', 'rxn_2', 'rxn_3', 'rxn_4', 'rxn_5', 'rxn_6'})
self.assertEqual(
set(self.model.compounds),
{Compound('A'), Compound('B'), Compound('C'), Compound('D', 'e')})
self.assertEqual(set(self.model.compounds), {
Compound('A', 'c'), Compound('B', 'c'), Compound('C', 'c'),
Compound('D', 'e')})

def test_add_reaction_invalid(self):
with self.assertRaises(Exception):
Expand All @@ -96,9 +96,8 @@ def test_remove_reaction_existing(self):
self.assertEqual(
set(self.model.reactions),
{'rxn_1', 'rxn_3', 'rxn_4', 'rxn_5', 'rxn_6'})
self.assertEqual(
set(self.model.compounds),
{Compound('A'), Compound('C'), Compound('D', 'e')})
self.assertEqual(set(self.model.compounds), {
Compound('A', 'c'), Compound('C', 'c'), Compound('D', 'e')})

def test_is_reversible_on_reversible(self):
self.assertTrue(self.model.is_reversible('rxn_2'))
Expand All @@ -121,20 +120,24 @@ def test_is_exchange_on_empty(self):
self.assertFalse(self.model.is_exchange('rxn_7'))

def test_add_all_database_reactions(self):
self.database.set_reaction('rxn_7', parse_reaction('|D| => |E|'))
added = self.model.add_all_database_reactions({None, 'e'})
self.assertEqual(added, { 'rxn_7' })
# Should get added
self.database.set_reaction('rxn_7', parse_reaction('D[c] => E[c]'))
# Not added because of compartment
self.database.set_reaction('rxn_8', parse_reaction('D[c] => E[p]'))
added = self.model.add_all_database_reactions({'c', 'e'})
self.assertEqual(added, {'rxn_7'})
self.assertEqual(set(self.model.reactions), {
'rxn_1', 'rxn_2', 'rxn_3', 'rxn_4', 'rxn_5', 'rxn_6', 'rxn_7'
})

def test_add_all_database_reactions_none(self):
added = self.model.add_all_database_reactions({None, 'e'})
added = self.model.add_all_database_reactions({'c', 'e'})
self.assertEqual(added, set())
self.assertEqual(set(self.model.reactions), { 'rxn_1', 'rxn_2', 'rxn_3', 'rxn_4', 'rxn_5', 'rxn_6' })
self.assertEqual(set(self.model.reactions), {
'rxn_1', 'rxn_2', 'rxn_3', 'rxn_4', 'rxn_5', 'rxn_6'})

def test_add_all_transport_reactions(self):
added = self.model.add_all_transport_reactions('e')
added = self.model.add_all_transport_reactions({('e', 'c')})
for reaction in added:
compartments = tuple(c.compartment for c, _ in
self.model.get_reaction_values(reaction))
Expand Down Expand Up @@ -215,7 +218,8 @@ def test_limits_set_both_flux_bounds_delete_upper(self):
self.assertEqual(self.model.limits['rxn_2'].bounds, (-10, 1000))

def test_limits_iter(self):
self.assertEqual(set(iter(self.model.limits)), { 'rxn_1', 'rxn_2', 'rxn_3', 'rxn_4', 'rxn_5', 'rxn_6' })
self.assertEqual(set(iter(self.model.limits)), {
'rxn_1', 'rxn_2', 'rxn_3', 'rxn_4', 'rxn_5', 'rxn_6'})

def test_limits_len(self):
self.assertEqual(len(self.model.limits), 6)
Expand Down

0 comments on commit ee2e9a6

Please sign in to comment.