Merge pull request #6 from lightdock/blocking_residues

Blocking residues

bin/lightdock3_setup.py

@@ -77,11 +77,13 @@
             # Calculate number of restraints in order to check them
             num_rec_active = len(restraints['receptor']['active'])
             num_rec_passive = len(restraints['receptor']['passive'])
+            num_rec_blocked = len(restraints['receptor']['blocked'])
             num_lig_active = len(restraints['ligand']['active'])
             num_lig_passive = len(restraints['ligand']['passive'])
 
             # Complain if not a single restraint has been defined, but restraints are enabled
-            if not num_rec_active and not num_rec_passive and not num_lig_active and not num_lig_passive:
+            if not num_rec_active and not num_rec_passive and not num_rec_blocked \
+                and not num_lig_active and not num_lig_passive:
                 raise LightDockError("Restraints file specified, but not a single restraint found")
 
             # Check if restraints correspond with real residues
@@ -93,12 +95,6 @@
             log.info("Number of receptor restraints is: %d (active), %d (passive)" % (num_rec_active, num_rec_passive))
             log.info("Number of ligand restraints is: %d (active), %d (passive)" % (num_lig_active, num_lig_passive))
 
-        rec_restraints = None
-        try:
-            rec_restraints = receptor_restraints['active'] + receptor_restraints['passive']
-        except:
-            pass
-
         lig_restraints = None
         try:
             lig_restraints = ligand_restraints['active'] + ligand_restraints['passive']
@@ -109,7 +105,7 @@
         starting_points_files = calculate_starting_positions(receptor, ligand, 
                                                              args.swarms, args.glowworms, 
                                                              args.starting_points_seed,
-                                                             rec_restraints, lig_restraints,
+                                                             receptor_restraints, lig_restraints,
                                                              rec_translation, lig_translation,
                                                              args.ftdock_file, args.use_anm, args.anm_seed,
                                                              args.anm_rec, args.anm_lig, args.membrane)

lightdock/prep/poses.py

@@ -97,6 +97,13 @@ def populate_poses(to_generate, center, radius, number_generator, rec_translatio
     """Creates new poses around a given center and a given radius"""
     new_poses = []
 
+    # Flatten if necessary receptor_restraints
+    if receptor_restraints:
+        try:
+            receptor_restraints = receptor_restraints['active'] + receptor_restraints['passive']
+        except TypeError:
+            pass
+
     # Calculate closest residue restraints
     closest_residues = []
     if receptor_restraints:
@@ -171,9 +178,10 @@ def create_file_from_poses(pos_file_name, poses):
     positions_file.close()
 
 
-def apply_restraints(swarm_centers, receptor_restraints, distance_cutoff, translation, 
-                     swarms_per_restraint=10):
-    """Filter out swarm centers which are not close to the given restraints"""
+def apply_restraints(swarm_centers, receptor_restraints, blocking_restraints, 
+                     distance_cutoff, translation, swarms_per_restraint=10):
+    """Filter out swarm centers which are not close to the given restraints or too close 
+    to blocking residues"""
     closer_swarms = []
     for i, residue in enumerate(receptor_restraints):
         distances = {}
@@ -200,7 +208,38 @@ def apply_restraints(swarm_centers, receptor_restraints, distance_cutoff, transl
 
     # Final filtered list of swarms
     new_swarm_centers = [swarm_centers[i] for i in closer_swarm_ids]
-    return new_swarm_centers
+
+    if not blocking_restraints:
+        return new_swarm_centers
+    else:
+        # We need to distinguish between two cases:
+        if len(closer_swarm_ids) > 0:
+            # We have a list of closer swarms to passive and active restraints, we need 
+            # to filter out from this list of closer swarms the ones closed to blocking
+            # restraints
+            centers_list = new_swarm_centers
+        else:
+            # From the total list of swarms, we will need to filter out the ones closed to
+            # the blocking residues. This is basically the same as the code for calculating
+            # closer_swarms
+            centers_list = swarm_centers
+
+        for i, residue in enumerate(blocking_restraints):
+            # We will use CA in case of protein, P in case of DNA
+            ca = residue.get_calpha()
+            if not ca:
+                ca = residue.get_atom('P')
+            to_remove = []
+            for center_id, center in enumerate(centers_list):
+                d = cdistance(ca.x + translation[0], ca.y + translation[1], ca.z + translation[2],
+                              center[0], center[1], center[2])
+                # Using ligand radius minus 5 Angstroms (10A is the swarm radius)
+                if d <= distance_cutoff - 5.:
+                    to_remove.append(center_id)
+            to_remove = list(set(to_remove))
+            centers_list = [centers_list[c] for c in range(len(centers_list)) if not c in to_remove]
+
+        return centers_list
 
 
 def estimate_membrane(z_coordinates, cutoff=10.0):
@@ -267,7 +306,8 @@ def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
                                                                                  is_membrane=is_membrane)
     # Filter swarms far from the restraints
     if receptor_restraints:
-        swarm_centers = apply_restraints(swarm_centers, receptor_restraints, 
+        regular_restraints = receptor_restraints['active'] + receptor_restraints['passive']
+        swarm_centers = apply_restraints(swarm_centers, regular_restraints, receptor_restraints['blocked'], 
                                          ligand_diameter / 2., rec_translation)
 
     # Filter out swarms which are not compatible with the explicit membrane

lightdock/test/prep/test_poses.py

@@ -6,12 +6,15 @@
 import numpy as np
 from nose.tools import assert_almost_equal
 from lightdock.prep.poses import normalize_vector, quaternion_from_vectors, get_quaternion_for_restraint, \
-        get_random_point_within_sphere, estimate_membrane, upper_layer, populate_poses, calculate_initial_poses
+        get_random_point_within_sphere, estimate_membrane, upper_layer, populate_poses, \
+        calculate_initial_poses, apply_restraints
 from lightdock.mathutil.cython.quaternion import Quaternion
 from lightdock.structure.residue import Residue
 from lightdock.mathutil.lrandom import MTGenerator, NMExtentGenerator
 from lightdock.pdbutil.PDBIO import parse_complex_from_file
 from lightdock.structure.complex import Complex
+from lightdock.structure.residue import Residue
+from lightdock.structure.atom import Atom
 from lightdock.prep.simulation import parse_restraints_file, get_restraints
 
 
@@ -355,3 +358,111 @@ def test_calculate_initial_poses(self):
 
         assert filecmp.cmp(positions_files[0], os.path.join(self.golden_data_path, '3p0g', 'init', 'initial_positions_0.dat'))
         assert filecmp.cmp(positions_files[1], os.path.join(self.golden_data_path, '3p0g', 'init', 'initial_positions_1.dat'))
+
+    def test_apply_restraints(self):
+        swarm_centers = [[13.861446721229493, -28.730453433364488, 0.024305878079852336],
+                         [-29.427540743780476, -2.960299745528308, -5.168493906227671],
+                         [24.244162925019417, 2.9322386764307304, -15.44286271885186],
+                         [2.0463687050893857, -13.94537658944639, -24.099372234598455],
+                         [-10.833120242061193, 13.517315520348813, -20.69794891072167],
+                         [-10.934383879949547, 24.302392364266915, 4.508648686944662],
+                         [-13.742376304140548, -18.137915011405447, 12.54861921589357],
+                         [-6.5003433948735, 7.824777098389764, 23.942006739050495],
+                         [17.345254847698744, 18.209853942307866, 6.384609403593564],
+                         [13.837506572006518, -7.152838167487336, 18.165186863034638]]
+        distance_cutoff = 19.522956365520056
+        translation = [-21.630461302211305, -5.510320024570013, -20.27116646191647]
+        atoms = [Atom(atom_number=1, atom_name='N', x=25.344, y=-10.048, z=26.435),
+                 Atom(atom_number=2, atom_name='CA', x=26.439, y=-10.798, z=26.994),
+                 Atom(atom_number=3, atom_name='C', x=26.646, y=-12.135, z=26.328),
+                 Atom(atom_number=4, atom_name='O', x=27.790, y=-12.531, z=26.028),
+                 Atom(atom_number=5, atom_name='CB', x=26.230, y=-10.933, z=28.503),
+                 Atom(atom_number=6, atom_name='CG', x=26.448, y=-9.540, z=29.191),
+                 Atom(atom_number=7, atom_name='CD', x=26.283, y=-9.676, z=30.689),
+                 Atom(atom_number=8, atom_name='CE', x=26.609, y=-8.392, z=31.454),
+                 Atom(atom_number=9, atom_name='NZ', x=26.210, y=-8.498, z=32.884)]
+        receptor_restraints = [Residue('LYS', 239, atoms=atoms)]
+
+        expected_new_centers = [[13.861446721229493, -28.730453433364488, 0.024305878079852336],
+                                [13.837506572006518, -7.152838167487336, 18.165186863034638]]
+
+        new_centers = apply_restraints(swarm_centers, receptor_restraints, [], distance_cutoff, translation)
+
+        assert np.allclose(expected_new_centers, new_centers)
+
+    def test_apply_restraints_with_blocking(self):
+        swarm_centers = [[13.861446721229493, -28.730453433364488, 0.024305878079852336],
+                         [-29.427540743780476, -2.960299745528308, -5.168493906227671],
+                         [24.244162925019417, 2.9322386764307304, -15.44286271885186],
+                         [2.0463687050893857, -13.94537658944639, -24.099372234598455],
+                         [-10.833120242061193, 13.517315520348813, -20.69794891072167],
+                         [-10.934383879949547, 24.302392364266915, 4.508648686944662],
+                         [-13.742376304140548, -18.137915011405447, 12.54861921589357],
+                         [-6.5003433948735, 7.824777098389764, 23.942006739050495],
+                         [17.345254847698744, 18.209853942307866, 6.384609403593564],
+                         [13.837506572006518, -7.152838167487336, 18.165186863034638]]
+        distance_cutoff = 19.522956365520056
+        translation = [-21.630461302211305, -5.510320024570013, -20.27116646191647]
+        atoms = [Atom(atom_number=1, atom_name='N', x=25.344, y=-10.048, z=26.435),
+                 Atom(atom_number=2, atom_name='CA', x=26.439, y=-10.798, z=26.994),
+                 Atom(atom_number=3, atom_name='C', x=26.646, y=-12.135, z=26.328),
+                 Atom(atom_number=4, atom_name='O', x=27.790, y=-12.531, z=26.028),
+                 Atom(atom_number=5, atom_name='CB', x=26.230, y=-10.933, z=28.503),
+                 Atom(atom_number=6, atom_name='CG', x=26.448, y=-9.540, z=29.191),
+                 Atom(atom_number=7, atom_name='CD', x=26.283, y=-9.676, z=30.689),
+                 Atom(atom_number=8, atom_name='CE', x=26.609, y=-8.392, z=31.454),
+                 Atom(atom_number=9, atom_name='NZ', x=26.210, y=-8.498, z=32.884)]
+        receptor_restraints = [Residue('LYS', 239, atoms=atoms)]
+
+        atoms = [Atom(atom_number=1, atom_name='N', x=30.795, y=-14.563, z=19.577),
+                 Atom(atom_number=2, atom_name='CA', x=30.969, y=-13.948, z=18.287),
+                 Atom(atom_number=3, atom_name='C', x=32.148, y=-12.988, z=18.291),
+                 Atom(atom_number=4, atom_name='O', x=32.601, y=-12.684, z=19.413),
+                 Atom(atom_number=5, atom_name='CB', x=29.742, y=-13.163, z=17.861),
+                 Atom(atom_number=6, atom_name='CG', x=28.586, y=-14.066, z=17.510),
+                 Atom(atom_number=7, atom_name='OD1', x=28.470, y=-14.548, z=16.359),
+                 Atom(atom_number=8, atom_name='ND2', x=27.742, y=-14.322, z=18.456)]
+        blocking_restraints = [Residue('ASN', 245, atoms=atoms)]
+        expected_new_centers = [[13.837506572006518, -7.152838167487336, 18.165186863034638]]
+
+        new_centers = apply_restraints(swarm_centers, receptor_restraints, blocking_restraints, 
+                                       distance_cutoff, translation)
+
+        assert np.allclose(expected_new_centers, new_centers)
+
+    def test_apply_restraints_only_blocking(self):
+        swarm_centers = [[13.861446721229493, -28.730453433364488, 0.024305878079852336],
+                         [-29.427540743780476, -2.960299745528308, -5.168493906227671],
+                         [24.244162925019417, 2.9322386764307304, -15.44286271885186],
+                         [2.0463687050893857, -13.94537658944639, -24.099372234598455],
+                         [-10.833120242061193, 13.517315520348813, -20.69794891072167],
+                         [-10.934383879949547, 24.302392364266915, 4.508648686944662],
+                         [-13.742376304140548, -18.137915011405447, 12.54861921589357],
+                         [-6.5003433948735, 7.824777098389764, 23.942006739050495],
+                         [17.345254847698744, 18.209853942307866, 6.384609403593564],
+                         [13.837506572006518, -7.152838167487336, 18.165186863034638]]
+        distance_cutoff = 19.522956365520056
+        translation = [-21.630461302211305, -5.510320024570013, -20.27116646191647]
+
+        atoms = [Atom(atom_number=1, atom_name='N', x=30.795, y=-14.563, z=19.577),
+                 Atom(atom_number=2, atom_name='CA', x=30.969, y=-13.948, z=18.287),
+                 Atom(atom_number=3, atom_name='C', x=32.148, y=-12.988, z=18.291),
+                 Atom(atom_number=4, atom_name='O', x=32.601, y=-12.684, z=19.413),
+                 Atom(atom_number=5, atom_name='CB', x=29.742, y=-13.163, z=17.861),
+                 Atom(atom_number=6, atom_name='CG', x=28.586, y=-14.066, z=17.510),
+                 Atom(atom_number=7, atom_name='OD1', x=28.470, y=-14.548, z=16.359),
+                 Atom(atom_number=8, atom_name='ND2', x=27.742, y=-14.322, z=18.456)]
+        blocking_restraints = [Residue('ASN', 245, atoms=atoms)]
+        expected_new_centers = [[-29.427540743780476, -2.960299745528308, -5.168493906227671],
+                                [24.244162925019417, 2.9322386764307304, -15.44286271885186],
+                                [2.0463687050893857, -13.94537658944639, -24.099372234598455],
+                                [-10.833120242061193, 13.517315520348813, -20.69794891072167],
+                                [-10.934383879949547, 24.302392364266915, 4.508648686944662],
+                                [-13.742376304140548, -18.137915011405447, 12.54861921589357],
+                                [-6.5003433948735, 7.824777098389764, 23.942006739050495],
+                                [17.345254847698744, 18.209853942307866, 6.384609403593564],
+                                [13.837506572006518, -7.152838167487336, 18.165186863034638]]
+
+        new_centers = apply_restraints(swarm_centers, [], blocking_restraints, distance_cutoff, translation)
+
+        assert np.allclose(expected_new_centers, new_centers)