Skip to content

Commit

Permalink
Add test to check that kb is not changed by model generator
Browse files Browse the repository at this point in the history
  • Loading branch information
@YinHoon
YinHoon committed Apr 17, 2019
1 parent a46b114 commit 76f9a04
Show file tree
Hide file tree
Showing 2 changed files with 65 additions and 34 deletions.
62 changes: 48 additions & 14 deletions tests/eukaryote/test_initialize_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,10 @@ def setUp(self):
'>chrX\nATGCGTCA\n'
'>chrM\nATGAARAARTTYCTCCTCACNCCNCTCTAATTT\n')

self.kb = wc_kb.KnowledgeBase()
cell = self.kb.cell = wc_kb.Cell()
self.kb = wc_kb.KnowledgeBase(id='test_kb', version='0.0.1')
cell = self.kb.cell = wc_kb.Cell(id='test_cell')

ref = wc_kb.core.Reference(id='ref', authors='John Smith', year=2018, comments='No comment')
ref = wc_kb.core.Reference(cell=cell, id='ref', authors='John Smith', year=2018, comments='No comment')
cell.parameters.create(id='cell_volume', value=10400., references=[ref])
cell.parameters.create(id='mean_doubling_time', value=20., units=unit_registry.parse_units('hour'))

Expand All @@ -65,43 +65,48 @@ def setUp(self):

chr1 = wc_kb.core.DnaSpeciesType(cell=cell, id='chr1', name='chromosome 1', ploidy=2,
sequence_path=self.sequence_path, circular=False, double_stranded=False)
gene1 = wc_kb.eukaryote_schema.GeneLocus(polymer=chr1, start=1, end=36)
gene1 = wc_kb.eukaryote_schema.GeneLocus(cell=cell, id='gene1', polymer=chr1, start=1, end=36)
exon1 = wc_kb.eukaryote_schema.GenericLocus(start=4, end=36)
transcript1 = wc_kb.eukaryote_schema.TranscriptSpeciesType(cell=cell, id='trans1',
name='transcript1', gene=gene1, exons=[exon1])
transcript1_spec = wc_kb.core.Species(species_type=transcript1, compartment=nucleus)
transcript1_conc = wc_kb.core.Concentration(cell=cell, species=transcript1_spec, value=0.02)
transcript1_conc.id = transcript1_conc.serialize()
cds1 = wc_kb.eukaryote_schema.GenericLocus(start=4, end=36)
prot1 = wc_kb.eukaryote_schema.ProteinSpeciesType(cell=cell, id='prot1', name='protein1',
transcript=transcript1, coding_regions=[cds1])
prot1_spec = wc_kb.core.Species(species_type=prot1, compartment=nucleus)
prot1_conc = wc_kb.core.Concentration(cell=cell, species=prot1_spec, value=0.03)
prot1_conc.id = prot1_conc.serialize()

chrX = wc_kb.core.DnaSpeciesType(cell=cell, id='chrX', name='chromosome X', ploidy=1,
sequence_path=self.sequence_path, circular=False, double_stranded=False)
gene2 = wc_kb.eukaryote_schema.GeneLocus(polymer=chrX, start=1, end=4)
gene2 = wc_kb.eukaryote_schema.GeneLocus(cell=cell, id='gene2', polymer=chrX, start=1, end=4)
exon2 = wc_kb.eukaryote_schema.GenericLocus(start=1, end=4)
transcript2 = wc_kb.eukaryote_schema.TranscriptSpeciesType(cell=cell, id='trans2',
name='transcript2', gene=gene2, exons=[exon2])
transcript2_spec = wc_kb.core.Species(species_type=transcript2, compartment=nucleus)
transcript2_conc = wc_kb.core.Concentration(cell=cell, species=transcript2_spec, value=0.01)
transcript2_conc = wc_kb.core.Concentration(cell=cell, species=transcript2_spec, value=0.01)
transcript2_conc.id = transcript2_conc.serialize()

chrM = wc_kb.core.DnaSpeciesType(cell=cell, id='chrM', name='mitochondrial chromosome', ploidy=150,
sequence_path=self.sequence_path, circular=False, double_stranded=False)
gene3 = wc_kb.eukaryote_schema.GeneLocus(polymer=chrM, start=1, end=33)
gene3 = wc_kb.eukaryote_schema.GeneLocus(cell=cell, id='gene3', polymer=chrM, start=1, end=33)
exon3 = wc_kb.eukaryote_schema.GenericLocus(start=1, end=30)
transcript3 = wc_kb.eukaryote_schema.TranscriptSpeciesType(cell=cell, id='trans3',
name='transcript3', gene=gene3, exons=[exon3])
transcript3_spec = wc_kb.core.Species(species_type=transcript3, compartment=mito)
transcript3_conc = wc_kb.core.Concentration(cell=cell, species=transcript3_spec, value=0.05)
transcript3_conc.id = transcript3_conc.serialize()
cds3 = wc_kb.eukaryote_schema.GenericLocus(start=1, end=30)
prot3 = wc_kb.eukaryote_schema.ProteinSpeciesType(cell=cell, id='prot3', name='protein3',
transcript=transcript3, coding_regions=[cds3])
prot3_spec = wc_kb.core.Species(species_type=prot3, compartment=mito)
prot3_conc = wc_kb.core.Concentration(cell=cell, species=prot3_spec, value=0.1)
prot3_conc.id = prot3_conc.serialize()

met1 = wc_kb.core.MetaboliteSpeciesType(cell=cell, id='met1', name='metabolite1')
met1.properties.append(wc_kb.core.SpeciesTypeProperty(property='structure', value=
met1_structure = wc_kb.core.SpeciesTypeProperty(property='structure', species_type=met1, value=
'InChI=1S'
'/C10H14N5O7P'
'/c11-8-5-9(13-2-12-8)15(3-14-5)10-7(17)6(16)4(22-10)1-21-23(18,19)20'
Expand All @@ -110,11 +115,12 @@ def setUp(self):
'/m1'
'/s1',
value_type=kb_onto['string']
))
)
met1_structure.id = met1_structure.gen_id()
met1_spec1 = wc_kb.core.Species(species_type=met1, compartment=nucleus)
met1_conc1 = wc_kb.core.Concentration(cell=cell, species=met1_spec1, value=0.5)
met1_conc1.id = met1_conc1.serialize()
met1_spec2 = wc_kb.core.Species(species_type=met1, compartment=extra)

species_type_coeff1 = wc_kb.core.SpeciesTypeCoefficient(species_type=prot1, coefficient=2)
species_type_coeff2 = wc_kb.core.SpeciesTypeCoefficient(species_type=met1, coefficient=3)
complex1 = wc_kb.core.ComplexSpeciesType(cell=cell, id='comp1', name='complex1',
Expand All @@ -126,8 +132,10 @@ def setUp(self):
subunits=[species_type_coeff3, species_type_coeff4])
complex2_spec1 = wc_kb.core.Species(species_type=complex2, compartment=nucleus)
complex2_conc1 = wc_kb.core.Concentration(cell=cell, species=complex2_spec1, value=0.)
complex2_conc1.id = complex2_conc1.serialize()
complex2_spec2 = wc_kb.core.Species(species_type=complex2, compartment=extra)
complex2_conc2 = wc_kb.core.Concentration(cell=cell, species=complex2_spec2, value=0.)
complex2_conc2.id = complex2_conc2.serialize()

expr1 = wc_kb.core.ObservableExpression(expression = '2.5 * prot1[n] + 1.3 * prot3[m]',
species = [prot1_spec, prot3_spec])
Expand Down Expand Up @@ -162,11 +170,13 @@ def setUp(self):
reaction=reaction1,
direction=wc_kb.core.RateLawDirection.forward,
expression=forward_exp)
forward_rate_law.id = forward_rate_law.gen_id()
backward_rate_law = wc_kb.core.RateLaw(
reaction=reaction1,
direction=wc_kb.core.RateLawDirection.backward,
expression=backward_exp)

expression=backward_exp)
backward_rate_law.id = backward_rate_law.gen_id()

def tearDown(self):
shutil.rmtree(self.tmp_dirname)

Expand Down Expand Up @@ -212,15 +222,20 @@ def test_gen_compartments_and_parameters(self):
self.assertEqual(model.parameters.get_one(id='mean_doubling_time').units, unit_registry.parse_units('s'))
self.assertEqual(model.parameters.get_one(id='mean_doubling_time').value, 20*3600)

def test_time_unit_conversion(self):

self.kb.cell.parameters.get_one(id='mean_doubling_time').units = unit_registry.parse_units('minute')
model = core.EukaryoteModelGenerator(self.kb,
component_generators=[initialize_model.InitializeModel],
options={'component': {'InitializeModel': self.set_options([])}}).run()
self.assertEqual(model.parameters.get_one(id='mean_doubling_time').units, unit_registry.parse_units('s'))
self.assertEqual(model.parameters.get_one(id='mean_doubling_time').value, 20*60)

self.kb.cell.parameters.get_one(id='mean_doubling_time').units = unit_registry.parse_units('s')
model = core.EukaryoteModelGenerator(self.kb,
component_generators=[initialize_model.InitializeModel],
options={'component': {'InitializeModel': self.set_options([])}}).run()
self.assertEqual(model.parameters.get_one(id='mean_doubling_time').units, unit_registry.parse_units('s'))
self.assertEqual(model.parameters.get_one(id='mean_doubling_time').value, 20)

def test_gen_dna(self):

model = core.EukaryoteModelGenerator(self.kb,
Expand Down Expand Up @@ -474,3 +489,22 @@ def test_gen_kb_rate_laws(self):
[model.observables.get_one(id='obs2')])
self.assertEqual(model.rate_laws.get_one(id='r1_kb-backward').expression.functions,
[model.functions.get_one(id='volume_n')])

def test_unchanged_kb(self):

kb2 = self.kb.copy()

model = core.EukaryoteModelGenerator(self.kb,
component_generators=[initialize_model.InitializeModel],
options={'component': {'InitializeModel': self.set_options([
'gen_dna',
'gen_transcripts',
'gen_protein',
'gen_metabolites',
'gen_complexes',
'gen_distribution_init_concentrations',
'gen_observables',
'gen_kb_reactions',
'gen_kb_rate_laws'])}}).run()

self.assertTrue(kb2.is_equal(self.kb))
37 changes: 17 additions & 20 deletions wc_model_gen/eukaryote/initialize_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,22 +133,7 @@ def gen_parameters(self):
Avogadro = model.parameters.create(id='Avogadro',
type = None,
value = scipy.constants.Avogadro,
units = unit_registry.parse_units('molecule mol^-1'))

# Standardize the units of doubling time
if kb.cell.parameters.get_one(id='mean_doubling_time'):
doubling_time_kb = kb.cell.parameters.get_one(id='mean_doubling_time')
else:
raise ValueError('The cell object does not have the parameter mean_doubling_time')

if not isinstance(doubling_time_kb.units, unit_registry.Unit):
ValueError('Unsupported units "{}"'.format(doubling_time_kb.units))

expr = unit_registry.parse_expression(str(doubling_time_kb.units))
scale = expr.to(unit_registry.parse_units('second'))
conversion_factor = scale.magnitude
doubling_time_kb.value *= conversion_factor
doubling_time_kb.units = unit_registry.parse_units('s')
units = unit_registry.parse_units('molecule mol^-1'))

# Create parameters from kb
for param in kb.cell.parameters:
Expand Down Expand Up @@ -181,7 +166,19 @@ def gen_parameters(self):
for identifier in param.identifiers:
identifier_model = wc_lang.Identifier(model=model,
namespace=identifier.namespace, id=identifier.id)
model_param.identifiers.append(identifier_model)
model_param.identifiers.append(identifier_model)

# Standardize the units of doubling time
if model.parameters.get_one(id='mean_doubling_time'):
model_doubling_time = model.parameters.get_one(id='mean_doubling_time')
else:
raise ValueError('The cell object does not have the parameter mean_doubling_time')

expr = unit_registry.parse_expression(str(model_doubling_time.units))
scale = expr.to(unit_registry.parse_units('second'))
conversion_factor = scale.magnitude
model_doubling_time.value *= conversion_factor
model_doubling_time.units = unit_registry.parse_units('s')

def gen_dna(self):
kb = self.knowledge_base
Expand Down Expand Up @@ -474,7 +471,7 @@ def gen_kb_rate_laws(self):
all_observables = {}
all_volumes = {}

kb_expression = kb_rate_law.expression.expression
model_expression = kb_rate_law.expression.expression

for observable in kb_rate_law.expression.observables:
all_observables[observable.id] = model.observables.get_one(id=observable.id)
Expand All @@ -492,10 +489,10 @@ def gen_kb_rate_laws(self):
if 'K_m' in param.id:
volume = model.compartments.get_one(id=param.id[param.id.rfind('_')+1:]).init_density.function_expressions[0].function
unit_adjusted_term = '{} * {} * {}'.format(param.id, Avogadro.id, volume.id)
kb_expression = kb_expression.replace(param.id, unit_adjusted_term)
model_expression = model_expression.replace(param.id, unit_adjusted_term)

rate_law_expression, error = wc_lang.RateLawExpression.deserialize(
kb_expression, {
model_expression, {
wc_lang.Parameter: all_parameters,
wc_lang.Species: all_species,
wc_lang.Observable: all_observables,
Expand Down

0 comments on commit 76f9a04

Please sign in to comment.