Testing new protocol for restraints

Author: brianjimenez

bin/post/lgd_filter_membrane.py

@@ -0,0 +1,139 @@
+#!/usr/bin/env python
+
+"""Filter LightDock final swarm results depending on the compatibility with the membrane"""
+
+from __future__ import print_function
+import sys
+import os
+import argparse
+import shutil
+import re
+from prody.measure.contacts import Contacts
+from prody import parsePDB, confProDy
+from lightdock.util.logger import LoggingManager
+from lightdock.util.analysis import read_ranking_file
+from lightdock.pdbutil.PDBIO import parse_complex_from_file
+from lightdock.structure.complex import Complex
+
+
+# Disable ProDy output
+confProDy(verbosity='info')
+filtered_folder = 'filtered'
+
+log = LoggingManager.get_logger('lgd_filter_membrane')
+
+
+def get_structures(ranking, base_path='.'):
+    structures = []
+    for rank in ranking:
+        swarm_id = rank.id_cluster
+        glowworm_id = rank.id_glowworm
+        structures.append(os.path.join(base_path, 
+                                       'swarm_{}'.format(swarm_id), 
+                                       'lightdock_{}.pdb'.format(glowworm_id)))
+    return structures
+
+
+def get_restraints(restraints_file):
+    restraints_receptor = set()
+    restraints_ligand = set()
+    with open(restraints_file) as handle:
+        for line in handle:
+            line = line.rstrip(os.linesep)
+            if line:
+                if line.startswith('R'):
+                    restraints_receptor.add(line.split()[-1])
+                if line.startswith('L'):
+                    restraints_ligand.add(line.split()[-1])
+    return restraints_receptor, restraints_ligand
+
+
+def calculate_membrane_height(parsed_receptor_file, restraints):
+    atoms, residues, chains = parse_complex_from_file(parsed_receptor_file)
+    receptor = Complex(chains, atoms)
+    z_coord = []
+    for restraint in restraints:
+        chain_id, residue_name, residue_number = restraint.split(".")
+        residue = receptor.get_residue(chain_id, residue_name, residue_number)
+        ca = residue.get_calpha()
+        z_coord.append(ca.z)
+    return min(z_coord)
+
+
+def parse_command_line():
+    """Parses command line arguments"""
+    parser = argparse.ArgumentParser(prog='lgd_filter_restraints')
+
+    parser.add_argument("ranking_file", help="Path of ranking to be used", metavar="ranking_file")
+    parser.add_argument("restraints_file", help="File including restraints", metavar="restraints_file")
+    parser.add_argument("parsed_receptor_file", help="Receptor PDB parsed by LightDock", metavar="parsed_receptor_file")
+    parser.add_argument("receptor_chains", help="Chains on the receptor partner", metavar="receptor_chains")
+    parser.add_argument("ligand_chains", help="Chains on the receptor partner", metavar="ligand_chains")
+    parser.add_argument("--cutoff", "-cutoff", "-c", help="Interaction cutoff",
+                            dest="cutoff", type=float, default=1.0)
+
+    return parser.parse_args()
+
+
+if __name__ == '__main__':
+ 
+    # Parse command line
+    args = parse_command_line()
+
+    log.info("Cutoff for membrane is {:3.1f}A".format(args.cutoff))
+
+    # Get ranking
+    ranking = read_ranking_file(args.ranking_file)
+
+    # Get all the PDB structures in a given directory
+    base_path = os.path.abspath(os.path.dirname(args.ranking_file))
+    structures = get_structures(ranking, base_path)
+
+    restraints_receptor, restraints_ligand = get_restraints(args.restraints_file)
+
+    membrane_height_z = calculate_membrane_height(args.parsed_receptor_file, restraints_receptor)
+
+    if os.path.exists(filtered_folder):
+        raise SystemExit("Folder {} already exists".format(filtered_folder))
+    else:
+        os.makedirs(filtered_folder)
+
+    filter_passed = {}
+    percentages = {}
+    for pdb_file in structures:
+        try:
+            swarm_id = int(re.findall(r'swarm_\d+', pdb_file)[0].split('_')[-1])
+            glowworm_id = int(re.findall(r'lightdock_\d+', pdb_file)[0].split('_')[-1])
+
+            # Read molecule and split by receptor and ligand
+            molecule = parsePDB(pdb_file)
+            ca_ligand = molecule.select('protein and chain {} and calpha'.format(args.ligand_chains))
+
+            # Contacts on ligand side
+            out = 0
+            for ca in ca_ligand:
+                coord = ca.getCoords()
+                if coord[-1] >= membrane_height_z:
+                    out += 1
+            perc = out / float(len(ca_ligand))
+            
+            if perc >= args.cutoff:
+            percentages[(swarm_id, glowworm_id)] = perc
+                shutil.copyfile(pdb_file, os.path.join(filtered_folder, 'swarm_{}_{}.pdb'.format(swarm_id, glowworm_id)))
+                try:
+                    filter_passed[swarm_id].append(glowworm_id)
+                except:
+                    filter_passed[swarm_id] = [glowworm_id]
+            print("{:40s}  {:5.3f}".format(pdb_file, perc))
+
+        except Exception, e:
+            log.error('Filtering has failed for structure {}. Please see error:'.format(pdb_file))
+            log.error(str(e))
+
+
+    filtered_ranking = os.path.join(filtered_folder, 'rank_filtered.list')
+    with open(filtered_ranking, 'w') as handle:
+        for rank in ranking:
+            if rank.id_cluster in filter_passed and rank.id_glowworm in filter_passed[rank.id_cluster]:
+                handle.write('swarm_{}_{}.pdb   {:5.3f}  {:5.3f}'.format(rank.id_cluster, 
+                    rank.id_glowworm, rank.scoring, percentages[(rank.id_cluster, rank.id_glowworm)]) + os.linesep)

bin/simulation/lightdock_setup.py

@@ -89,7 +89,7 @@
                                                              receptor_restraints, ligand_restraints, 
                                                              rec_translation, lig_translation,
                                                              args.ftdock_file, args.use_anm, args.anm_seed,
-                                                             args.anm_rec, args.anm_lig)
+                                                             args.anm_rec, args.anm_lig, args.membrane)
         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)

lightdock/prep/poses.py

@@ -2,6 +2,7 @@
 
 import os
 import operator
+import numpy as np
 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
@@ -24,16 +25,40 @@ def get_random_point_within_sphere(number_generator, radius):
         if x**2 + y**2 + z**2 <= r2:
             return x, y, z
 
+def get_quaternion_for_restraint(rec_residue, lig_residue, cx, cy, cz):
+    """Calculates the quaternion required for orienting the ligand towards the restraint"""
+    r_ca = rec_residue.get_calpha()
+    l_ca = lig_residue.get_calpha()
+    a = np.array([r_ca.x, r_ca.y, r_ca.z])  
+    b = np.array([l_ca.x-cx, l_ca.y-cy, l_ca.z-cz])
+    c = np.cross(a, b)
+    d = np.dot(a, b)
 
-def populate_poses(to_generate, center, radius, number_generator, rng_nm=None, rec_nm=0, lig_nm=0):
+    s = np.sqrt( (1+d)*2 )
+    invs = 1. / s
+    x = c[0] * invs
+    y = c[1] * invs
+    z = c[2] * invs
+    w = s * 0.5
+
+    return Quaternion(w=w, x=x, y=y, z=z).normalize()
+
+
+def populate_poses(to_generate, center, radius, number_generator, rng_nm=None, rec_nm=0, lig_nm=0,
+                    receptor_restraints=None, ligand_restraints=None):
     """Creates new poses around a given center and a given radius"""
     new_poses = []
     for _ in xrange(to_generate):
         x, y, z = get_random_point_within_sphere(number_generator, radius)
         tx = center[0] + x
         ty = center[1] + y
         tz = center[2] + z
-        q = Quaternion.random(number_generator)
+        if receptor_restraints and ligand_restraints:
+            rec_residue = receptor_restraints[number_generator.randint(0, len(receptor_restraints)-1)]
+            lig_residue = ligand_restraints[number_generator.randint(0, len(ligand_restraints)-1)]
+            q = get_quaternion_for_restraint(rec_residue, lig_residue, center[0], center[1], center[2])
+        else:
+            q = Quaternion.random(number_generator)
         op_vector = [tx, ty, tz, q.w, q.x, q.y, q.z]
         if rng_nm:
             if rec_nm > 0:
@@ -53,7 +78,8 @@ 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):
+def apply_restraints(swarm_centers, receptor_restraints, distance_cutoff, translation, 
+                     is_membrane=False, swarms_per_restraint=10):
 
     closer_swarms = []
     for i, residue in enumerate(receptor_restraints):
@@ -72,14 +98,25 @@ def apply_restraints(swarm_centers, receptor_restraints, distance_cutoff, transl
             if swarms_considered == swarms_per_restraint:
                 break
     closer_swarms = list(set(closer_swarms))
-    return closer_swarms
+    if is_membrane:
+        # Requieres the receptor to be aligned in the Z axis with the membrane
+        min_z = min([residue.get_calpha().z for residue in receptor_restraints]) + translation[2]
+        compatible = []
+        for swarm_id, center in enumerate(swarm_centers):
+            if swarm_id in closer_swarms:
+                if center[2] >= min_z:
+                    compatible.append(swarm_id)
+        return compatible
+    else:
+        return closer_swarms
 
 
 def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
                             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
+                            dest_folder, ftdock_file='', nm_mode=False, nm_seed=0, rec_nm=0, lig_nm=0,
+                            is_membrane=False):
+    """Calculates the starting points for each of the glowworms using the center of swarms
     and FTDock poses.
     """
     # Random number generator for poses
@@ -91,18 +128,19 @@ def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
     else:
         rng_nm = None
 
-    # Calculate cluster centers
-    cluster_centers, receptor_diameter, ligand_diameter = calculate_surface_points(receptor, 
+    # Calculate swarm centers
+    swarm_centers, receptor_diameter, ligand_diameter = calculate_surface_points(receptor, 
                                                                                    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]
+        filtered_swarms = apply_restraints(swarm_centers, receptor_restraints, ligand_diameter / 2., 
+                                           rec_translation, is_membrane=is_membrane)
+        swarm_centers = [swarm_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)
+    create_pdb_from_points(pdb_file_name, swarm_centers)
 
     ligand_center = ligand.center_of_coordinates()
     radius = 10.    # ligand_diameter / 2.
@@ -114,7 +152,7 @@ def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
             raise NotImplementedError('Using FTDock poses with NM is not supported')
 
         poses = FTDockCoordinatesParser.get_list_of_poses(ftdock_file, ligand_center)
-        clusters = classify_ftdock_poses(poses, cluster_centers, radius)
+        clusters = classify_ftdock_poses(poses, swarm_centers, radius)
 
         for cluster_id, ftdock_poses in clusters.iteritems():
             # Translate FTDock poses into lightdock poses
@@ -131,7 +169,7 @@ def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
             # Populate new poses if needed
             if len(poses) < num_glowworms:
                 needed = num_glowworms - len(poses)
-                poses.extend(populate_poses(needed, cluster_centers[cluster_id], radius, rng))
+                poses.extend(populate_poses(needed, swarm_centers[cluster_id], radius, rng))
 
             # Save poses as pdb file
             pdb_file_name = os.path.join(dest_folder, '%s_%s.pdb' % (DEFAULT_PDB_STARTING_PREFIX, cluster_id))
@@ -144,8 +182,9 @@ def calculate_initial_poses(receptor, ligand, num_clusters, num_glowworms,
             positions_files.append(pos_file_name)
             create_bild_file(bild_file_name, poses)
     else:
-        for cluster_id, cluster_center in enumerate(cluster_centers):
-            poses = populate_poses(num_glowworms, cluster_center, radius, rng, rng_nm, rec_nm, lig_nm)
+        for cluster_id, cluster_center in enumerate(swarm_centers):
+            poses = populate_poses(num_glowworms, cluster_center, radius, rng, rng_nm, rec_nm, lig_nm,
+                                    receptor_restraints, ligand_restraints)
             # Save poses as pdb file
             pdb_file_name = os.path.join(dest_folder, '%s_%s.pdb' % (DEFAULT_PDB_STARTING_PREFIX, cluster_id))
             create_pdb_from_points(pdb_file_name, [[pose[0], pose[1], pose[2]] for pose in poses[:num_glowworms]])

lightdock/prep/simulation.py

@@ -94,7 +94,8 @@ def check_starting_file(file_name, glowworms, use_anm, anm_rec, anm_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):
+    use_anm=False, anm_seed=0, anm_rec=DEFAULT_NMODES_REC, anm_lig=DEFAULT_NMODES_LIG,
+    is_membrane=False):
     """Defines the starting positions of each glowworm in the simulation.
 
     If the init_folder already exists, uses the starting positions from this folder.
@@ -110,7 +111,8 @@ def calculate_starting_positions(receptor, ligand, swarms, glowworms, starting_p
                                                         rec_translation, lig_translation,
                                                         init_folder,
                                                         ftdock_file, use_anm,
-                                                        anm_seed, anm_rec, anm_lig)
+                                                        anm_seed, anm_rec, anm_lig,
+                                                        is_membrane)
         log.info("Generated %d positions files" % len(starting_points_files))
     else:
         if receptor_restraints:

lightdock/scoring/fastdfire/driver.py

@@ -140,9 +140,11 @@ def __call__(self, receptor, receptor_coordinates, ligand, ligand_coordinates):
         energy, interface_receptor, interface_ligand = calculate_dfire(receptor, ligand, 
                                                                        self.potential.dfire_energy, 
                                                                        receptor_coordinates, ligand_coordinates)
-        perc_receptor_restraints = ScoringFunction.restraints_satisfied(receptor.restraints, set(interface_receptor))
-        perc_ligand_restraints = ScoringFunction.restraints_satisfied(ligand.restraints, set(interface_ligand))
-        return energy + perc_receptor_restraints * energy + perc_ligand_restraints * energy
+        # Code to consider contacts in the interface
+        #perc_receptor_restraints = ScoringFunction.restraints_satisfied(receptor.restraints, set(interface_receptor))
+        #perc_ligand_restraints = ScoringFunction.restraints_satisfied(ligand.restraints, set(interface_ligand))
+        #return energy + perc_receptor_restraints * energy + perc_ligand_restraints * energy
+        return energy
 
 # Needed to dynamically load the scoring functions from command line
 DefinedScoringFunction = DFIRE

lightdock/util/parser.py

@@ -101,6 +101,9 @@ def __init__(self):
         parser.add_argument("-rst", "--rst", help="Restraints file", 
                             dest="restraints", type=CommandLineParser.valid_file,
                             metavar="restraints", default=None)
+        # Membrane setup
+        parser.add_argument("-membrane", "--membrane", help="Enables the extra filter for membrane restraints", 
+                            dest="membrane", action='store_true', default=False)
 
         self.args = parser.parse_args()
 

lightdock/version.py

@@ -1,3 +1,3 @@
 """Framework version"""
 
-CURRENT_VERSION = "0.5.3"
+CURRENT_VERSION = "0.5.4"