Apply restraints on receptor

Author: brianjimenez

bin/simulation/lightdock_setup.py

@@ -13,7 +13,8 @@
 from lightdock.util.parser import SetupCommandLineParser
 from lightdock.prep.simulation import read_input_structure, save_lightdock_structure, \
                                       calculate_starting_positions, prepare_results_environment, \
-                                      create_setup_file, calculate_anm
+                                      create_setup_file, calculate_anm, parse_restraints_file, \
+                                      get_restraints
 from lightdock.constants import DEFAULT_LIGHTDOCK_PREFIX, DEFAULT_ELLIPSOID_DATA_EXTENSION, \
                                 DEFAULT_NMODES_REC, DEFAULT_REC_NM_FILE, DEFAULT_NMODES_LIG, DEFAULT_LIG_NM_FILE
 from lightdock.mathutil.ellipsoid import MinimumVolumeEllipsoid
@@ -35,8 +36,8 @@
         ligand = read_input_structure(args.ligand_pdb, args.noxt)
 
         # Move structures to origin
-        receptor.move_to_origin()
-        ligand.move_to_origin()
+        rec_translation = receptor.move_to_origin()
+        lig_translation = ligand.move_to_origin()
 
         # Calculate reference points for receptor
         log.info("Calculating reference points for receptor %s..." % args.receptor_pdb)
@@ -63,12 +64,34 @@
             calculate_anm(receptor, args.anm_rec, DEFAULT_REC_NM_FILE)
             calculate_anm(ligand, args.anm_lig, DEFAULT_LIG_NM_FILE)
 
+        # Parse restraints if any:
+        if args.restraints:
+            log.info("Reading restraints from %s" % args.restraints)
+            restraints = parse_restraints_file(args.restraints)
+            # Check if restraints have been defined for the ligand, but not to the receptor
+            if len(restraints['ligand']) and not len(restrains['receptor']):
+                raise LightDockError("Restraints defined for ligand, but not receptor. Try switching structures.")
+
+            if not len(restraints['receptor']) and not len(restraints['ligand']):
+                raise LightDockError("Restraints file specified, but not a single restraint found")
+
+            # Check if restraints correspond with real residues
+            receptor_restraints = get_restraints(receptor, restraints['receptor'])
+            args.receptor_restraints = restraints['receptor']
+            ligand_restraints = get_restraints(ligand, restraints['ligand'])
+            args.ligand_restraints = restraints['ligand']
+
         # Calculate surface points (swarm centers) over receptor structure
         starting_points_files = calculate_starting_positions(receptor, ligand, 
                                                              args.swarms, args.glowworms, 
-                                                             args.starting_points_seed, 
+                                                             args.starting_points_seed,
+                                                             receptor_restraints, ligand_restraints, 
+                                                             rec_translation, lig_translation,
                                                              args.ftdock_file, args.use_anm, args.anm_seed,
                                                              args.anm_rec, args.anm_lig)
+        if len(starting_points_files) != args.swarms:
+            args.swarms = len(starting_points_files)
+            log.info('Number of swarms is %d after applying restraints' % args.swarms)
 
         # Create simulation folders
         prepare_results_environment(args.swarms)

lightdock/prep/poses.py

@@ -1,11 +1,13 @@
 """Module to prepare initial poses for docking"""
 
 import os
+import operator
 from lightdock.pdbutil.PDBIO import create_pdb_from_points
 from lightdock.prep.starting_points import calculate_surface_points
 from lightdock.prep.ftdock import FTDockCoordinatesParser, classify_ftdock_poses
 from lightdock.mathutil.lrandom import MTGenerator, NormalGenerator
 from lightdock.mathutil.cython.quaternion import Quaternion
+from lightdock.mathutil.cython.cutil import distance as cdistance
 from lightdock.constants import CLUSTERS_CENTERS_FILE,\
     DEFAULT_PDB_STARTING_PREFIX, DEFAULT_STARTING_PREFIX, DEFAULT_BILD_STARTING_PREFIX, DEFAULT_EXTENT_MU, \
     DEFAULT_EXTENT_SIGMA
@@ -51,8 +53,32 @@ 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):
+    
+    closer_swarms = []
+    for i, residue in enumerate(receptor_restraints):
+        distances = {}
+        ca = residue.get_calpha()
+        for swarm_id, center in enumerate(swarm_centers):
+            distances[swarm_id] = cdistance(ca.x + translation[0], ca.y + translation[1], ca.z + translation[2],
+                                            center[0], center[1], center[2])
+        sorted_distances = sorted(distances.items(), key=operator.itemgetter(1))
+        swarms_considered = 0
+        for swarm in sorted_distances:
+            swarm_id, distance = swarm[0], swarm[1]
+            if distance <= distance_cutoff:
+                closer_swarms.append(swarm_id)
+                swarms_considered += 1
+            if swarms_considered == swarms_per_restraint:
+                break
+    closer_swarms = list(set(closer_swarms))
+    return closer_swarms
+
+
 def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
-                            seed, dest_folder, ftdock_file='', nm_mode=False, nm_seed=0, rec_nm=0, lig_nm=0):
+                            seed, receptor_restraints, ligand_restraints, 
+                            rec_translation, lig_translation,
+                            dest_folder, ftdock_file='', nm_mode=False, nm_seed=0, rec_nm=0, lig_nm=0):
     """Calculates the starting points for each of the glowworms using the center of clusters
     and FTDock poses.
     """
@@ -70,6 +96,11 @@ def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
                                                                                    ligand, 
                                                                                    num_clusters,
                                                                                    distance_step=1.0)
+    # Filter cluster centers far from the restraints
+    if receptor_restraints:
+        filtered_swarms = apply_restraints(cluster_centers, receptor_restraints, ligand_diameter / 2., rec_translation)
+        cluster_centers = [cluster_centers[i] for i in filtered_swarms]
+
     pdb_file_name = os.path.join(dest_folder, CLUSTERS_CENTERS_FILE)
     create_pdb_from_points(pdb_file_name, cluster_centers)
 

lightdock/prep/simulation.py

@@ -92,8 +92,9 @@ def check_starting_file(file_name, glowworms, use_anm, anm_rec, anm_lig):
         return num_glowworms == glowworms
 
 
-def calculate_starting_positions(receptor, ligand, swarms, glowworms, starting_points_seed, 
-    ftdock_file=None, use_anm=False, anm_seed=0, anm_rec=DEFAULT_NMODES_REC, anm_lig=DEFAULT_NMODES_LIG):
+def calculate_starting_positions(receptor, ligand, swarms, glowworms, starting_points_seed,
+    receptor_restraints, ligand_restraints, rec_translation, lig_translation, ftdock_file=None, 
+    use_anm=False, anm_seed=0, anm_rec=DEFAULT_NMODES_REC, anm_lig=DEFAULT_NMODES_LIG):
     """Defines the starting positions of each glowworm in the simulation.
 
     If the init_folder already exists, uses the starting positions from this folder.
@@ -104,12 +105,19 @@ def calculate_starting_positions(receptor, ligand, swarms, glowworms, starting_p
         os.mkdir(init_folder)
         starting_points_files = calculate_initial_poses(receptor, ligand,
                                                         swarms, glowworms,
-                                                        starting_points_seed, init_folder,
+                                                        starting_points_seed, 
+                                                        receptor_restraints, ligand_restraints,
+                                                        rec_translation, lig_translation,
+                                                        init_folder,
                                                         ftdock_file, use_anm,
                                                         anm_seed, anm_rec, anm_lig)
         log.info("Generated %d positions files" % len(starting_points_files))
     else:
-        log.warning("Folder %s already exists, skipping calculation" % init_folder)
+        if receptor_restraints:
+            log.warning("Folder %s already exists and restraints apply. Check for consistency, possible error!" % init_folder)
+        else:
+            log.warning("Folder %s already exists, skipping calculation" % init_folder)
+
         pattern = os.path.join(DEFAULT_POSITIONS_FOLDER, "%s*.dat" % DEFAULT_STARTING_PREFIX)
         starting_points_files = glob.glob(pattern)
         if len(starting_points_files) != swarms:
@@ -197,3 +205,49 @@ def get_default_box(use_anm, anm_rec, anm_lig):
         boundaries.extend([Boundary(MIN_EXTENT, MAX_EXTENT) for _ in xrange(anm_lig)])
 
     return BoundingBox(boundaries)
+
+
+def parse_restraints_file(restraints_file_name):
+    with open(restraints_file_name) as input_restraints:
+        raw_restraints = [line.rstrip(os.linesep) for line in input_restraints.readlines()]
+        restraints = {'receptor': [], 'ligand': []}
+        for restraint in raw_restraints:
+            if restraint and restraint[0] in ['R', 'L']:
+                try:
+                    fields = restraint[2:].split('.')
+                    # Only consider first character if many
+                    chain_id = fields[0][0].upper()
+                    # Only considering 3 chars if many
+                    residue = fields[1][0:3].upper()
+                    # Check for integer
+                    residue_number = int(fields[2])
+                    parsed_restraint = "%s.%s.%d" % (chain_id, residue, residue_number)
+                    if restraint[0] == 'R':
+                        restraints['receptor'].append(parsed_restraint)
+                    else:
+                        restraints['ligand'].append(parsed_restraint)
+                except:
+                    log.warning('Ignoring restraints %s' % restraint)
+
+            # Make unique the restraints:
+            restraints['receptor'] = list(set(restraints['receptor']))
+            restraints['ligand'] = list(set(restraints['ligand']))
+
+        return restraints
+
+
+def get_restraints(structure, restraints):
+    """Check for each restraint in the format Chain.ResidueName.ResidueNumber in 
+    restraints if they exist in the given structure.
+    """
+    residues = []
+    for restraint in restraints:
+        chain_id, residue_name, residue_number = restraint.split('.')
+        residue = structure.get_residue(chain_id, residue_name, residue_number)
+        residue
+        if not residue:
+            raise LightDockError("Restraint %s not found in structure" % restraint)
+        else:
+            residues.append(residue)
+    return residues
+

lightdock/structure/complex.py

@@ -108,6 +108,15 @@ def move_to_origin(self):
         """Moves the structure to the origin of coordinates"""
         translation = [-1*c for c in self.center_of_coordinates()]
         self.translate(translation)
+        return translation
+
+    def get_residue(self, chain_id, residue_name, residue_number):
+        for chain in self.chains:
+            if chain_id == chain.cid:
+                for residue in chain.residues:
+                    if residue.name == residue_name and int(residue.number) == int(residue_number):
+                        return residue
+        return None
 
     def __getitem__(self, item):
         return self.atom_coordinates[item]

lightdock/structure/residue.py

@@ -85,6 +85,10 @@ def get_atom(self, atom_name):
                 return atom
         return None
 
+    def get_calpha(self):
+        """Get the Calpha atom"""
+        return self.get_atom('CA')
+
     def get_non_hydrogen_atoms(self):
         return [atom for atom in self.atoms if not atom.is_hydrogen()]
 

lightdock/util/parser.py

@@ -97,6 +97,10 @@ def __init__(self):
         parser.add_argument("-anm_lig", "--anm_lig", help="Number of ANM modes for ligand", 
                             type=SetupCommandLineParser.valid_integer_number,
                             dest="anm_lig", default=DEFAULT_NMODES_LIG)
+        # Restraints file
+        parser.add_argument("-rst", "--rst", help="Restraints file", 
+                            dest="restraints", type=CommandLineParser.valid_file,
+                            metavar="restraints", default=None)
 
         self.args = parser.parse_args()