Merge pull request #3 from brianjimenez/restraints

Restraints

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,35 @@
             calculate_anm(receptor, args.anm_rec, DEFAULT_REC_NM_FILE)
             calculate_anm(ligand, args.anm_lig, DEFAULT_LIG_NM_FILE)
 
+        # Parse restraints if any:
+        receptor_restraints = ligand_restraints = None
+        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)

docs/README.md

@@ -59,7 +59,7 @@ If you execute `lightdock_setup` several options will appear:
 ```bash
 usage: lightdock_setup [-h] [--seed_points STARTING_POINTS_SEED]
                        [-ft ftdock_file] [--noxt] [-anm] [--seed_anm ANM_SEED]
-                       [-anm_rec ANM_REC] [-anm_lig ANM_LIG]
+                       [-anm_rec ANM_REC] [-anm_lig ANM_LIG] [-rst restraints]
                        receptor_pdb_file ligand_pdb_file swarms glowworms
 lightdock_setup: error: too few arguments
 
@@ -85,6 +85,7 @@ Below there is a description of the rest of accepted paramenters of `lightdock_s
 - **--seed_anm** *ANM_SEED*: An integer can be specified as the seed used in the random number generator of ANM normal modes extents. 
 - **--anm_rec** *ANM_REC*: The number of non-trivial normal modes calculated for the recepetor in the ANM mode.
 - **--anm_lig** *ANM_LIG*: The number of non-trivial normal modes calculated for the ligand in the ANM mode.
+- **-rst** *restraints_file*: If `restraints_file` is provided, distance restraints on the receptor surface will be calculated to provided residues. See `2.2. Restraints` section for more information.
 
 ### 2.1. Results of the setup
 
@@ -101,7 +102,22 @@ After the execution of `lightdock_setup` script, several files and directories w
 - `*.xyz.npy`: Two files with this extension, one for the receptor and one for the ligand, which contains information about the minimum ellipsoid containing each of the structures in NumPy format.
 - `lightdock_rec.nm.npy`and `lightdock_lig.nm.npy`: If ANM is enabled, two files are created containing the normal modes information calculated by the ProDy library.
 
-### 2.2. Tips and tricks
+### 2.2. Distance restraints (Experimental)
+
+In version `0.5.1`, distance restraints on the receptor surface have been implemented. This new feature works in the following way:
+
+From a `restraints_file` file containing the following information:
+
+```
+R A.SER.150 
+R A.TYR.151
+```
+
+Where `R` means for `Receptor` (at the moment only residue restraints on the receptor are accepted) and then the rest of the line is a residue identifier in the format `Chain.Residue_Name.Residue_Number` in the original PDB file numeration.
+
+For each of these residues, only the closest swarms are considered during `lightdock_setup`. At the moment, only the **10 closest swarms** for each residue are considered. If some of the swarms overlap, then only one copy of that swarm is used.
+
+### 2.3. Tips and tricks
 
 - As a general rule of thumb, the receptor structure is the bigger molecule and the ligand the smaller. With bigger and smaller, the metric used is the longest diameter of the minimum ellipsoid containing the molecule. A script called `lgd_calculate_diameter.py` can be found in `$LIGHTDOCK_HOME/bin/support` path in order to calculate an approximation of that metric.
 

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,48 @@ 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)
+        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()
 

lightdock/version.py

@@ -1,3 +1,3 @@
 """Framework version"""
 
-CURRENT_VERSION = "0.5.0"
+CURRENT_VERSION = "0.5.1"