Support for different number of ANM modes

Author: brianjimenez

bin/post/lgd_generate_conformations.py

@@ -13,6 +13,7 @@
 from lightdock.structure.complex import Complex
 from lightdock.mathutil.cython.quaternion import Quaternion
 from lightdock.structure.nm import read_nmodes
+from lightdock.prep.simulation import get_setup_from_file
 
 
 log = LoggingManager.get_logger('generate_conformations')
@@ -59,7 +60,7 @@ def get_lightdock_structures(input_file):
     return file_names
 
 
-def parse_output_file(lightdock_output):
+def parse_output_file(lightdock_output, num_anm_rec, num_anm_lig):
     translations = []
     rotations = []
     receptor_ids = []
@@ -79,9 +80,9 @@ def parse_output_file(lightdock_output):
             coord = line[1:last].split(',')
             translations.append([float(coord[0]), float(coord[1]), float(coord[2])])
             rotations.append(Quaternion(float(coord[3]), float(coord[4]), float(coord[5]), float(coord[6])))
-            if len(coord) == (7 + DEFAULT_NMODES_REC + DEFAULT_NMODES_LIG):
-                rec_extents.append(np.array([float(x) for x in coord[7:7+DEFAULT_NMODES_REC]]))
-                lig_extents.append(np.array([float(x) for x in coord[-DEFAULT_NMODES_LIG:]]))
+            if len(coord) > 7:
+                rec_extents.append(np.array([float(x) for x in coord[7:7+num_anm_rec]]))
+                lig_extents.append(np.array([float(x) for x in coord[-num_anm_lig:]]))
             raw_data = line[last+1:].split()
             receptor_id = int(raw_data[0])
             ligand_id = int(raw_data[1])
@@ -105,9 +106,21 @@ def parse_output_file(lightdock_output):
                         type=valid_file, metavar="lightdock_output")
     # Number of glowworms
     parser.add_argument("glowworms", help="number of glowworms", type=valid_integer_number)
+    # Optional, setup file
+    parser.add_argument("--setup", "-setup", "-s", help="Simulation setup file",
+                            dest="setup_file", metavar="setup_file", type=valid_file, default=None)
 
     args = parser.parse_args()
 
+    # Load setup configuration if provided
+    setup = get_setup_from_file(args.setup_file) if args.setup_file else None
+
+    num_anm_rec = DEFAULT_NMODES_REC
+    num_anm_lig = DEFAULT_NMODES_LIG
+    if setup and setup['use_anm']:
+        num_anm_rec = setup['anm_rec']
+        num_anm_lig = setup['anm_lig']
+
     # Receptor
     structures = []
     for structure in get_lightdock_structures(args.receptor_structures):
@@ -128,7 +141,7 @@ def parse_output_file(lightdock_output):
 
     # Output file
     translations, rotations, receptor_ids, ligand_ids, \
-        rec_extents, lig_extents = parse_output_file(args.lightdock_output)
+        rec_extents, lig_extents = parse_output_file(args.lightdock_output, num_anm_rec, num_anm_lig)
 
     found_conformations = len(translations)
     num_conformations = args.glowworms
@@ -157,22 +170,30 @@ def parse_output_file(lightdock_output):
         # Use normal modes if provided:
         if len(rec_extents):
             try:
-                for nm in range(DEFAULT_NMODES_REC):
+                for nm in range(num_anm_rec):
                     receptor_pose.coordinates += nmodes_rec[nm] * rec_extents[i][nm]
-            except ValueError, e:
+            except ValueError:
                 log.error("Problem found on calculating ANM for receptor:")
                 log.error("Number of atom coordinates is: %s" % str(receptor_pose.coordinates.shape))
                 log.error("Number of ANM is: %s" % str(nmodes_rec.shape))
                 raise SystemExit
+            except IndexError:
+                log.error("Problem found on calculating ANM for receptor:")
+                log.error("If you have used anm_rec different than default, please use --setup")
+                raise SystemExit
         if len(lig_extents):
             try:
-                for nm in range(DEFAULT_NMODES_LIG):
+                for nm in range(num_anm_lig):
                     ligand_pose.coordinates += nmodes_lig[nm] * lig_extents[i][nm]
-            except ValueError, e:
+            except ValueError:
                 log.error("Problem found on calculating ANM for ligand:")
                 log.error("Number of atom coordinates is: %s" % str(receptor_pose.coordinates.shape))
                 log.error("Number of ANM is: %s" % str(nmodes_rec.shape))
                 raise SystemExit
+            except IndexError:
+                log.error("Problem found on calculating ANM for ligand:")
+                log.error("If you have used anm_lig different than default, please use --setup")
+                raise SystemExit
 
         # We rotate first, ligand it's at initial position
         ligand_pose.rotate(rotations[i])

bin/simulation/docking_mpi.py

@@ -38,6 +38,8 @@ def set_gso(number_of_glowworms, adapters, scoring_functions, initial_positions,
     else:
         gso_parameters = GSOParameters()
     builder = LightdockGSOBuilder()
+    if not use_anm:
+        anm_rec = anm_lig = 0
     gso = builder.create_from_file(number_of_glowworms, random_number_generator, gso_parameters,
                                    adapters, scoring_functions, bounding_box, initial_positions,
                                    step_translation, step_rotation, nmodes_step, local_minimization)
@@ -112,8 +114,16 @@ def run_simulation(parser):
         ligand.move_to_origin()
 
         if args.use_anm:
-            receptor.n_modes = read_nmodes("%s%s" % (DEFAULT_REC_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
-            ligand.n_modes = read_nmodes("%s%s" % (DEFAULT_LIG_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
+            try:
+                receptor.n_modes = read_nmodes("%s%s" % (DEFAULT_REC_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
+            except:
+                log.warning("No ANM found for receptor molecule")
+                receptor.n_modes = None
+            try:
+                ligand.n_modes = read_nmodes("%s%s" % (DEFAULT_LIG_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
+            except:
+                log.warning("No ANM found for ligand molecule")
+                ligand.n_modes = None
 
         starting_points_files = load_starting_positions(args.swarms, args.glowworms, args.use_anm,
                                                         args.anm_rec, args.anm_lig)

bin/simulation/docking_multiprocessing.py

@@ -37,9 +37,12 @@ def set_gso(number_of_glowworms, adapters, scoring_functions, initial_positions,
     else:
         gso_parameters = GSOParameters()
     builder = LightdockGSOBuilder()
+    if not use_anm:
+        anm_rec = anm_lig = 0
     gso = builder.create_from_file(number_of_glowworms, random_number_generator, gso_parameters,
                                    adapters, scoring_functions, bounding_box, initial_positions,
-                                   step_translation, step_rotation, nmodes_step, local_minimization)
+                                   step_translation, step_rotation, nmodes_step, local_minimization,
+                                   anm_rec, anm_lig)
     return gso
 
 
@@ -121,8 +124,16 @@ def run_simulation(parser):
         ligand.move_to_origin()
 
         if args.use_anm:
-            receptor.n_modes = read_nmodes("%s%s" % (DEFAULT_REC_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
-            ligand.n_modes = read_nmodes("%s%s" % (DEFAULT_LIG_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
+            try:
+                receptor.n_modes = read_nmodes("%s%s" % (DEFAULT_REC_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
+            except:
+                log.warning("No ANM found for receptor molecule")
+                receptor.n_modes = None
+            try:
+                ligand.n_modes = read_nmodes("%s%s" % (DEFAULT_LIG_NM_FILE, NUMPY_FILE_SAVE_EXTENSION) )
+            except:
+                log.warning("No ANM found for ligand molecule")
+                ligand.n_modes = None
 
         starting_points_files = load_starting_positions(args.swarms, args.glowworms, args.use_anm,
                                                         args.anm_rec, args.anm_lig)

bin/simulation/lightdock_setup.py

@@ -61,8 +61,12 @@
 
         # Calculate and save ANM if required
         if args.use_anm:
-            calculate_anm(receptor, args.anm_rec, DEFAULT_REC_NM_FILE)
-            calculate_anm(ligand, args.anm_lig, DEFAULT_LIG_NM_FILE)
+            if args.anm_rec > 0:
+                log.info("Calculating ANM for receptor molecule...")
+                calculate_anm(receptor, args.anm_rec, DEFAULT_REC_NM_FILE)
+            if args.anm_lig > 0:
+                log.info("Calculating ANM for ligand molecule...")
+                calculate_anm(ligand, args.anm_lig, DEFAULT_LIG_NM_FILE)
 
         # Parse restraints if any:
         receptor_restraints = ligand_restraints = None

lightdock/gso/algorithm.py

@@ -107,13 +107,14 @@ class LightdockGSOBuilder(object):
     """Creates a GSO simulation object for the LightDock framework"""
     def create_from_file(self, number_of_glowworms, random_number_generator, gso_parameters,
                          adapters, scoring_functions, bounding_box, initial_population_file,
-                         step_translation, step_rotation, step_nmodes, local_minimization):
+                         step_translation, step_rotation, step_nmodes, local_minimization,
+                         anm_rec, anm_lig):
         """Creates a new GSO instance of the algorithm reading the initial position of the glowworms
         agents from initial_population_file and using the scoring function adapter.
         """
         self._initializer = LightdockFromFileInitializer(adapters, scoring_functions, number_of_glowworms,
                                                          gso_parameters, bounding_box.dimension,
                                                          initial_population_file, step_translation, step_rotation,
-                                                         random_number_generator, step_nmodes)
+                                                         random_number_generator, step_nmodes, anm_rec, anm_lig)
         return GSO(self._initializer.generate_glowworms(), gso_parameters, random_number_generator,
                    local_minimization=local_minimization)

lightdock/gso/initializer.py

@@ -88,7 +88,7 @@ class LightdockFromFileInitializer(Initializer):
     """This initializer takes into account the complex provided by the adapter"""
     def __init__(self, adapters, scoring_functions, number_of_glowworms, gso_parameters,
                  dimensions, initial_population_file, step_translation, step_rotation,
-                 random_number_generator, step_nmodes):
+                 random_number_generator, step_nmodes, anm_rec, anm_lig):
         super(LightdockFromFileInitializer, self).__init__(scoring_functions, number_of_glowworms, gso_parameters)
         self.dimensions = dimensions
         self.initial_population_file = initial_population_file
@@ -98,6 +98,8 @@ def __init__(self, adapters, scoring_functions, number_of_glowworms, gso_paramet
         self.step_nmodes = step_nmodes
         # Patch to not mess with old simulations
         self.random_number_generator = MTGenerator(random_number_generator.seed)
+        self.anm_rec = anm_rec
+        self.anm_lig = anm_lig
 
     def generate_landscape_positions(self):
         """Generates a list of landscape positions that have been read
@@ -122,6 +124,5 @@ def generate_landscape_positions(self):
                                                              adapter.receptor_model, adapter.ligand_model,
                                                              receptor_index, ligand_index,
                                                              self.step_translation, self.step_rotation,
-                                                             self.step_nmodes))
-
+                                                             self.step_nmodes, self.anm_rec, self.anm_lig))
         return positions

lightdock/gso/searchspace/landscape.py

@@ -97,24 +97,22 @@ class DockingLandscapePosition(LandscapePosition):
     """
     def __init__(self, scoring_function, coordinates, receptor, ligand, receptor_id=0, ligand_id=0,
                  step_translation=DEFAULT_TRANSLATION_STEP, step_rotation=DEFAULT_ROTATION_STEP,
-                 step_nmodes=DEFAULT_NMODES_STEP):
+                 step_nmodes=0, num_rec_nmodes=0, num_lig_nmodes=0):
         self.objective_function = scoring_function
         self.translation = np.array(coordinates[:3])
         self.rotation = Quaternion(coordinates[3], coordinates[4], coordinates[5], coordinates[6])
-        # Copy ANM information if required
-        if len(coordinates) > 7 and len(coordinates) == (7 + DEFAULT_NMODES_REC + DEFAULT_NMODES_LIG):
-            self.rec_extent = np.array(coordinates[7:7+DEFAULT_NMODES_REC])
-            self.lig_extent = np.array(coordinates[-DEFAULT_NMODES_LIG:])
-        else:
-            self.rec_extent = np.array([])
-            self.lig_extent = np.array([])
         self.receptor = receptor
         self.ligand = ligand
         self.receptor_id = receptor_id
         self.ligand_id = ligand_id
         self.step_translation = step_translation
         self.step_rotation = step_rotation
         self.step_nmodes = step_nmodes
+        self.num_rec_nmodes = num_rec_nmodes
+        self.num_lig_nmodes = num_lig_nmodes
+        # Copy ANM information if required
+        self.rec_extent = np.array(coordinates[7:7+self.num_rec_nmodes]) if self.num_rec_nmodes > 0 else np.array([])
+        self.lig_extent = np.array(coordinates[-self.num_lig_nmodes:]) if self.num_lig_nmodes > 0 else np.array([])
         # This part is important, each position needs to retain its own pose coordinates
         self.receptor_pose = self.receptor.coordinates[self.receptor_id].clone()
         self.ligand_pose = self.ligand.coordinates[self.ligand_id].clone()
@@ -129,7 +127,8 @@ def clone(self):
         return DockingLandscapePosition(self.objective_function, coordinates,
                                         self.receptor, self.ligand,
                                         self.receptor_id, self.ligand_id,
-                                        self.step_translation, self.step_rotation, self.step_nmodes)
+                                        self.step_translation, self.step_rotation, self.step_nmodes,
+                                        self.num_rec_nmodes, self.num_lig_nmodes)
 
     def evaluate_objective_function(self, receptor_structure_id=None, ligand_structure_id=None):
         """Evaluates the objective function at the given coordinates"""
@@ -147,13 +146,11 @@ def evaluate_objective_function(self, receptor_structure_id=None, ligand_structu
         self.ligand_reference_points = self.ligand.reference_points.clone()
 
         # Use normal modes if provided:
-        num_rec_nmodes = len(self.rec_extent)
-        if num_rec_nmodes:
-            for i in range(num_rec_nmodes):
+        if self.num_rec_nmodes > 0:
+            for i in range(self.num_rec_nmodes):
                 self.receptor_pose.coordinates += self.receptor.n_modes[i] * self.rec_extent[i]
-        num_lig_nmodes = len(self.lig_extent)
-        if num_lig_nmodes:
-            for i in range(num_lig_nmodes):
+        if self.num_lig_nmodes > 0:
+            for i in range(self.num_lig_nmodes):
                 self.ligand_pose.coordinates += self.ligand.n_modes[i] * self.lig_extent[i]
 
         # We rotate first, ligand it's at initial position
@@ -198,14 +195,14 @@ def move(self, other):
             # Rotation (Quaternion SLERP)
             self.rotation = self.rotation.slerp(other.rotation, self.step_rotation)
             # NModes
-            if len(self.rec_extent):
+            if self.num_rec_nmodes > 0:
                 delta_x = other.rec_extent - self.rec_extent
                 n = np.linalg.norm(delta_x)
                 # Only move if required
                 if not np.allclose([0.0], [n]):
                     delta_x *= (self.step_nmodes / n)
                     self.rec_extent += delta_x
-            if len(self.lig_extent):
+            if self.num_lig_nmodes > 0:
                 delta_x = other.lig_extent - self.lig_extent
                 n = np.linalg.norm(delta_x)
                 # Only move if required
@@ -238,8 +235,8 @@ def update_landscape_position(self, optimized_vector):
         """Updates the current pose"""
         self.translation = optimized_vector[:3]
         self.rotation = Quaternion(optimized_vector[3], optimized_vector[4], optimized_vector[5], optimized_vector[6])
-        self.rec_extent = optimized_vector[7:7+DEFAULT_NMODES_REC]
-        self.lig_extent = optimized_vector[-DEFAULT_NMODES_LIG:]
+        self.rec_extent = optimized_vector[7:7+self.num_rec_nmodes] if self.num_rec_nmodes > 0 else np.array([])
+        self.lig_extent = optimized_vector[-self.num_lig_nmodes:] if self.num_lig_nmodes > 0 else np.array([])
 
     def minimize(self):
         """Returns the new scoring after minimizing this landscape position using a local non-grandient

lightdock/test/gso/test_initializer.py

@@ -109,7 +109,7 @@ def test_create_swarm(self):
         initializer = LightdockFromFileInitializer([self.adapter], [self.scoring_function],
                                                    number_of_glowworms, gso_parameters, 
                                                    7, self.golden_data_path+'initial_positions_1PPE.txt', 
-                                                   0.5, 0.5, random_number_generator, 0.5)
+                                                   0.5, 0.5, random_number_generator, 0.5, 10, 10)
         swarm = initializer.generate_glowworms()
 
         assert number_of_glowworms == swarm.get_size()
@@ -123,7 +123,7 @@ def test_generate_landscape_positions_without_coordinates(self):
         initializer = LightdockFromFileInitializer(self.adapter, self.scoring_function, 
                                                    number_of_glowworms, gso_parameters, 
                                                    7, self.golden_data_path+'initial_positions_empty.txt', 
-                                                   0.5, 0.5, random_number_generator, 0.5)
+                                                   0.5, 0.5, random_number_generator, 0.5, 10, 10)
         swarm = initializer.generate_glowworms()
 
         assert swarm.get_size() > 0
@@ -137,7 +137,7 @@ def test_generate_landscape_positions_num_glowworms_different(self):
         initializer = LightdockFromFileInitializer(self.adapter, self.scoring_function, 
                                                    number_of_glowworms, gso_parameters, 
                                                    7, self.golden_data_path+'initial_positions_1PPE.txt', 
-                                                   0.5, 0.5, random_number_generator, 0.5)
+                                                   0.5, 0.5, random_number_generator, 0.5, 10, 10)
         swarm = initializer.generate_glowworms()
 
         assert swarm.get_size() > 0

lightdock/test/gso/test_lightdockbuilder.py

@@ -53,7 +53,7 @@ def test_LightDockGSOBuilder_using_FromFileInitializer(self):
         gso = builder.create_from_file(number_of_glowworms, self.random_number_generator, 
                                        self.gso_parameters, [adapter], [scoring_function], self.bounding_box,
                                        self.golden_data_path+'initial_positions_1PPE.txt',
-                                       0.5, 0.5, 0.5, local_minimization=False)
+                                       0.5, 0.5, 0.5, False, 10, 10)
 
         assert 5 == gso.swarm.get_size()