add(genes.rotamers): 'with_original' flag to control whether to inclu…

…de inputfile chi angles in rotamer library

gaudi/genes/rotamers.py

@@ -4,17 +4,17 @@
 ##############
 # GaudiMM: Genetic Algorithms with Unrestricted
 # Descriptors for Intuitive Molecular Modeling
-# 
+#
 # https://github.com/insilichem/gaudi
 #
 # Copyright 2017 Jaime Rodriguez-Guerra, Jean-Didier Marechal
-# 
+#
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
-# 
+#
 #      http://www.apache.org/licenses/LICENSE-2.0
-# 
+#
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -79,6 +79,10 @@ class Rotamers(GeneProvider):
     library : {'Dunbrack', 'Dynameomics'}
         The rotamer library to use.
 
+    with_original: bool, defaults to True
+        Whether to include the original set of chi angles as part of the
+        rotamer library.
+
     Attributes
     ----------
     allele : list of float
@@ -94,11 +98,12 @@ class Rotamers(GeneProvider):
     # to not have any rotamers
     _residues_without_rotamers = set(('ALA', 'GLY'))
 
-    def __init__(self, residues=None, library='Dunbrack', **kwargs):
+    def __init__(self, residues=None, library='Dunbrack', with_original=True, **kwargs):
         GeneProvider.__init__(self, **kwargs)
         self._kwargs = kwargs
         self._residues = residues
         self.library = library
+        self.with_original = with_original
         self.allele = []
         # set caches
         try:
@@ -110,7 +115,7 @@ def __init__(self, residues=None, library='Dunbrack', **kwargs):
         except KeyError:
             self.rotamers = self._cache[self.name + '_rotamers'] = OrderedDict()
 
-        
+
     def __ready__(self):
         """
         Second stage of initialization.
@@ -125,7 +130,7 @@ def __ready__(self):
                 self.patch_residue(r)
                 self.residues[(molname, r.id.position)] = r
                 self.allele.append(random.random())
-    
+
     def __deepcopy__(self, memo):
         new = self.__class__(residues=self._residues, library=self.library, **self._kwargs)
         new.residues = self.residues
@@ -136,7 +141,7 @@ def express(self):
         for ((molname, pos), residue), i in zip(self.residues.items(), self.allele):
             if residue.type not in self._residues_without_rotamers:
                 try:
-                    rotamers = self.retrieve_rotamers(molname, pos, residue, 
+                    rotamers = self.retrieve_rotamers(molname, pos, residue,
                                                       library=self.library.title())
                 except NoResidueRotamersError:  # ALA, GLY...
                     logger.warn('%s/%s (%s) has no rotamers', molname, pos, residue.type)
@@ -166,10 +171,12 @@ def sort_by_rmsd(ref, query):
 
             chis = getRotamerParams(residue, lib=library)[2]
             rotamers_mols = getRotamers(residue, lib=library)[1]
+            reference = residue if self.with_original else rotamers_mols[0].residues[0]
             rotamers_and_chis = zip(rotamers_mols, [c.chis for c in chis])
-            rotamers_and_chis.sort(key=lambda rc: sort_by_rmsd(residue, rc[0]))
-            rotamers = (self.all_chis(residue),) + zip(*rotamers_and_chis)[1]
-            print(len(rotamers))
+            rotamers_and_chis.sort(key=lambda rc: sort_by_rmsd(reference, rc[0]))
+            rotamers = zip(*rotamers_and_chis)[1]
+            if self.with_original:
+                rotamers = (self.all_chis(residue),) + rotamers
             self.rotamers[(molecule, position)] = rotamers
             for rot in rotamers_mols:
                 rot.destroy()

tests/test_genes_rotamers.py

@@ -4,17 +4,17 @@
 ##############
 # GaudiMM: Genetic Algorithms with Unrestricted
 # Descriptors for Intuitive Molecular Modeling
-# 
+#
 # https://github.com/insilichem/gaudi
 #
 # Copyright 2017 Jaime Rodriguez-Guerra, Jean-Didier Marechal
-# 
+#
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
-# 
+#
 #      http://www.apache.org/licenses/LICENSE-2.0
-# 
+#
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -29,7 +29,7 @@
 
 def rotamers(individual, path, position, seed):
     individual.genes['Molecule'] = Molecule(parent=individual, path=datapath(path))
-    individual.genes['Rotamers'] = rotamers = Rotamers(parent=individual,
+    individual.genes['Rotamers'] = rotamers = Rotamers(parent=individual, with_original=True,
                                                        residues=[('Molecule', position)])
     individual.__ready__()
     individual.__expression_hooks__()
@@ -38,9 +38,9 @@ def rotamers(individual, path, position, seed):
 
 
 @pytest.mark.parametrize("path, position, seed, restype, original_chis, new_chis", [
-    ('4c3w_protein.mol2', 5, 0, 'ARG', [179.734, 178.061, 60.608, 90.076], 
+    ('4c3w_protein.mol2', 5, 0, 'ARG', [179.734, 178.061, 60.608, 90.076],
                                        [179.734, 178.061, 60.608, 90.076]),
-    ('4c3w_protein.mol2', 5, 0.015, 'ARG', [179.734, 178.061, 60.608, 90.076], 
+    ('4c3w_protein.mol2', 5, 0.015, 'ARG', [179.734, 178.061, 60.608, 90.076],
                                            [-174.7, -174.4, 69.2, 79.1]),
 ])
 def test_rotamers(individual, path, position, seed, restype, original_chis, new_chis):
@@ -49,7 +49,7 @@ def test_rotamers(individual, path, position, seed, restype, original_chis, new_
     # Cache initial coordinates for the residue and ALL the alpha carbons in the protein
     alpha_carbons_coords = [a.xformCoord() for a in residue.molecule.atoms if a.name == 'CA']
     residue_coords = [a.xformCoord() for a in residue.atoms]
-    
+
     with expressed(individual):
         assert residue.id.position == position
         assert residue.type == restype