ENH: implement mitochondrial expression system

etfl/core/genes.py

@@ -14,6 +14,7 @@
 from cobra import Gene
 from Bio.Seq import Seq
 from Bio.Alphabet import DNAAlphabet, RNAAlphabet, ProteinAlphabet
+from etfl.optim.constraints import SynthesisConstraint
 
 
 def make_sequence(sequence, seq_type):
@@ -31,7 +32,7 @@ def make_sequence(sequence, seq_type):
     else:
         raise TypeError('The type of the sequence argument should be either '
                         'string or a Bio.Seq')
-
+    
     return typed_sequence
 
 class ExpressedGene(Gene):
@@ -77,7 +78,27 @@ def transcribed_by(self):
     @transcribed_by.setter
     def transcribed_by(self,value):
         # TODO: Make this a setter that rewrites the adequate constraints
-        raise NotImplementedError()
+        if value != self._transcribed_by:
+            if self.model is None:
+                # Easy
+                self._transcribed_by = value
+            else:
+                # We need to make the model change the corresponding cstr
+                mod_id = self.id + '_transcription'
+                Cstr = self.model.get_constraints_of_type(SynthesisConstraint)
+                cstr_exist = 1
+                try:
+                    Cstr.get_by_id(mod_id)
+                except KeyError:
+                    cstr_exist = 0
+                if cstr_exist:
+                    # trancription constraint already exists
+                    raise NotImplementedError()
+                else:
+                    self._transcribed_by = value
+        else:
+            # Nothing to do here :)
+            pass
 
     @property
     def translated_by(self):
@@ -86,7 +107,27 @@ def translated_by(self):
     @translated_by.setter
     def translated_by(self,value):
         # TODO: Make this a setter that rewrites the adequate constraints
-        raise NotImplementedError()
+        if value != self._translated_by:
+            if self.model is None:
+                # Easy
+                self._translated_by = value
+            else:
+                # We need to make the model change the corresponding cstr
+                mod_id = self.id + '_translation'
+                Cstr = self.model.get_constraints_of_type(SynthesisConstraint)
+                cstr_exist = 1
+                try:
+                    Cstr.get_by_id(mod_id)
+                except KeyError:
+                    cstr_exist = 0
+                if cstr_exist:
+                    # trancription constraint already exists
+                    raise NotImplementedError()
+                else:
+                    self._translated_by = value
+        else:
+            # Nothing to do here :)
+            pass
 
 
     @property

etfl/core/memodel.py

@@ -259,24 +259,54 @@ def add_nucleotide_sequences(self, sequences):
         :param sequences:
         :return:
         """
-
+        
         for gene_id, seq in sequences.items():
             if gene_id in self.genes:
                 new = replace_by_me_gene(self, gene_id, seq)
-
+                
             else:
                 self.logger.warning('Model has no gene {}, Adding it'.format(gene_id))
                 new = ExpressedGene(id= gene_id, name = gene_id, sequence=seq)
                 self.add_genes([new])
 
             self._make_peptide_from_gene(gene_id)
+
+
+    def add_transcription_by(self, transcription_dict):
+
+        for gene_id, transcripted_by in transcription_dict.items():
+            # transcripted_by is a list of rnap(s)
+            try:
+                self.genes.get_by_id(gene_id).transcribed_by = transcripted_by
+            except KeyError:
+                continue
+
+    def add_translation_by(self, translation_dict):
+
+        for gene_id, translated_by in translation_dict.items():
+            # transcripted_by is a list of rnap(s)
+            try:
+                self.genes.get_by_id(gene_id).translated_by = translated_by
+            except KeyError:
+                continue
 
     def _make_peptide_from_gene(self, gene_id):
         free_pep = Peptide(id=gene_id,
                            name='Peptide, {}'.format(gene_id),
                            gene_id=gene_id)
         free_pep._model = self
         self.peptides += [free_pep]
+
+    def add_peptide_sequences(self, aa_sequences):
+
+        for pep_id, seq in aa_sequences.items():
+            if pep_id in self.peptides:
+                self.peptides.get_by_id(pep_id).peptide = seq
+            else:
+                self.logger.warning('Model has no peptide {}'.format(pep_id))
+                continue
+
+
 
     def add_dummies(self, nt_ratios, mrna_kdeg, mrna_length, aa_ratios,
                     enzyme_kdeg, peptide_length):
@@ -1405,36 +1435,41 @@ def _populate_rnap(self):
 
         self.regenerate_variables()
 
-        usage = self._get_rnap_total_capacity()
-
-        # usage = (sum_RMs + self.RNAPf[this_rnap.id].unscaled - this_rnap.concentration) / \
-        #         this_rnap.scaling_factor
-
-        # Create the capacity constraint
-        self.add_constraint(kind=TotalCapacity,
-                            hook=self,
-                            id_ = 'rnap',
-                            expr=usage,
-                            lb = 0,
-                            ub = 0,
-                            )
+        for rnap_id in self.rnap.keys():
+            usage = self._get_rnap_total_capacity(rnap_id = rnap_id)
+
+            # usage = (sum_RMs + self.RNAPf[this_rnap.id].unscaled - this_rnap.concentration) / \
+            #         this_rnap.scaling_factor
+
+            # Create the capacity constraint
+            self.add_constraint(kind=TotalCapacity,
+                                hook=self,
+                                id_ = rnap_id,
+                                expr=usage,
+                                lb = 0,
+                                ub = 0,
+                                )
 
         # update variable and constraints attributes
         self.regenerate_constraints()
         self.regenerate_variables()
 
-    def _get_rnap_total_capacity(self):
+    def _get_rnap_total_capacity(self, rnap_id):
         all_rnap_usage = self.get_variables_of_type(RNAPUsage)
-        sum_RMs = symbol_sum([x.unscaled for x in all_rnap_usage])
-        free_rnap = [self.get_variables_of_type(FreeEnzyme).get_by_id(x)
-                     for x in self.rnap]
+        # only for genes trascribed by this rnap
+        sum_RMs = symbol_sum([x.unscaled for x in all_rnap_usage \
+                             if x.hook.transcribed_by is None \
+                             or rnap_id in x.hook.transcribed_by])
+        free_rnap = [self.get_variables_of_type(FreeEnzyme).get_by_id(rnap_id)]
+
         # The total RNAP capacity constraint looks like
         # ΣRMi + Σ(free RNAPj) = Σ(Total RNAPj)
         usage = sum_RMs \
                 + sum([x.unscaled for x in free_rnap]) \
-                - sum([x.concentration for x in self.rnap.values()])
-        usage /= min([x.scaling_factor for x in self.rnap.values()])
+                - sum([self.rnap[rnap_id].concentration])
+        usage /= min([self.rnap[rnap_id].scaling_factor])
         return usage
+
 
     def apply_rnap_catalytic_constraint(self, reaction, queue):
         """
@@ -1615,41 +1650,46 @@ def _populate_ribosomes(self):
         # 3 -> Add ribosomal capacity constraint
         self.regenerate_variables()
 
-        free_ribosome = [self.get_variables_of_type(FreeEnzyme).get_by_id(x)
-                     for x in self.ribosome]
-        # CATCH : This is summing ~1500+ variable objects, and for a reason
-        # sympy does not like it. Let's cut it in smaller chunks and sum
-        # afterwards
-        # sum_RPs = sum(self.get_variables_of_type(RibosomeUsage))
-        all_ribosome_usage = self.get_variables_of_type(RibosomeUsage)
-
-        # sum_RPs = chunk_sum(all_ribosome_usage)
-        sum_RPs = symbol_sum([x.unscaled for x in all_ribosome_usage])
-
-        # The total RNAP capacity constraint looks like
-        # ΣRMi + Σ(free RNAPj) = Σ(Total RNAPj)
-        ribo_usage = sum_RPs \
-                + sum([x.unscaled for x in free_ribosome]) \
-                - sum([x.concentration for x in self.ribosome.values()])
-        ribo_usage /= min([x.scaling_factor for x in self.ribosome.values()])
-        # For nondimensionalization
-        # ribo_usage = (sum_RPs + self.Rf.unscaled - self.ribosome.concentration) \
-        #              / self.ribosome.scaling_factor
-
-        # Create the capacity constraint
-        self.add_constraint(kind=TotalCapacity,
-                            hook=self,
-                            id_='rib',
-                            expr=ribo_usage,
-                            lb = 0,
-                            ub = 0,
-                            )
+
+        for rib_id in self.ribosome.keys():
+            # CATCH : This is summing ~1500+ variable objects, and for a reason
+            # sympy does not like it. Let's cut it in smaller chunks and sum
+            # afterwards
+            # sum_RPs = sum(self.get_variables_of_type(RibosomeUsage))
+            ribo_usage = self._get_rib_total_capacity(rib_id = rib_id)
+            # For nondimensionalization
+            # ribo_usage = (sum_RPs + self.Rf.unscaled - self.ribosome.concentration) \
+            #              / self.ribosome.scaling_factor
+
+            # Create the capacity constraint
+            self.add_constraint(kind=TotalCapacity,
+                                hook=self,
+                                id_=rib_id,
+                                expr=ribo_usage,
+                                lb = 0,
+                                ub = 0,
+                                )
 
         # update variable and constraints attributes
         self.regenerate_constraints()
         self.regenerate_variables()
 
+    def _get_rib_total_capacity(self, rib_id):
+        all_ribosome_usage = self.get_variables_of_type(RibosomeUsage)
+        # only for genes traslated by this ribosome
+        sum_RPs = symbol_sum([x.unscaled for x in all_ribosome_usage \
+                              if x.hook.translated_by is None \
+                              or rib_id in x.hook.translated_by])
+        free_ribosome = [self.get_variables_of_type(FreeEnzyme).get_by_id(rib_id)]
 
+        # The total RNAP capacity constraint looks like
+        # ΣRMi + Σ(free RNAPj) = Σ(Total RNAPj)
+        usage = sum_RPs \
+                + sum([x.unscaled for x in free_ribosome]) \
+                - sum([self.ribosome[rib_id].concentration])
+        usage /= min([self.ribosome[rib_id].scaling_factor])
+        return usage
+
     def apply_ribosomal_catalytic_constraint(self, reaction):
         """
         Given a translation reaction, apply the constraint that links it with

etfl/io/dict.py

@@ -12,6 +12,9 @@
 
 from tqdm import tqdm
 
+from Bio.Seq import Seq
+from Bio.Alphabet import DNAAlphabet, RNAAlphabet, ProteinAlphabet
+
 import cobra.io.dict as cbd
 from cobra.exceptions import SolverNotFound
 from optlang.util import expr_to_json, parse_expr
@@ -51,10 +54,18 @@ def expressed_gene_to_dict(gene):
     obj['rna'] = str(gene.rna)
     obj['peptide'] = str(gene.peptide)
     obj['copy_number'] = str(gene.copy_number)
-    obj['transcribed_by'] = [x.id for x in gene.transcribed_by] \
-        if gene.transcribed_by is not None else None
-    obj['translated_by'] = [x.id for x in gene.translated_by]\
-        if gene.translated_by is not None else None
+    try:
+        obj['transcribed_by'] = [x.id for x in gene.transcribed_by] \
+            if gene.transcribed_by is not None else None
+    except AttributeError:
+        obj['transcribed_by'] = [x for x in gene.transcribed_by] \
+            if gene.transcribed_by is not None else None
+    try:
+        obj['translated_by'] = [x.id for x in gene.translated_by]\
+            if gene.translated_by is not None else None
+    except AttributeError:
+        obj['translated_by'] = [x for x in gene.translated_by] \
+            if gene.translated_by is not None else None            
 
     return obj
 
@@ -537,7 +548,7 @@ def init_me_model_from_dict(new, obj):
     new.growth_reaction = obj['growth_reaction']
 
     # Populate enzymes
-    # new.coupling_dict = rebuild_coupling_dict(new, obj['coupling_dict'])
+    # new.coupling_dict = rebuild_coupling_dict(, obj['coupling_dict'])
     new.add_enzymes([enzyme_from_dict(x) for x in obj['enzymes']], prep=False)
 
     # Make RNAP
@@ -892,8 +903,8 @@ def find_genes_from_dict(new, obj):
             g = replace_by_me_gene(new, gene_dict['id'], str(sequence))
             if key == 'expressed_genes':
                 # Newer models
-                g._rna            = gene_dict['rna']
-                g._peptide        = gene_dict['peptide']
+                g._rna            = Seq(gene_dict['rna'], alphabet=RNAAlphabet())
+                g._peptide        = Seq(gene_dict['peptide'], alphabet=ProteinAlphabet())
                 g._copy_number    = gene_dict['copy_number']
                 g._transcribed_by = [new.enzymes.get_by_id(e)
                                      for e in gene_dict['transcribed_by']] \