Merge pull request #23 from lightdock/0.9.2

Merges branch 0.9.2 into master. Branch 0.9.2 adds support for flip flag, fixes reference counting errors and updates to the last version of ProDy.

README.md

@@ -4,7 +4,7 @@
 [![Supported versions](https://img.shields.io/pypi/pyversions/lightdock.svg)](https://pypi.org/project/lightdock)
 [![Build Status](https://travis-ci.com/lightdock/lightdock.svg?branch=master)](https://travis-ci.com/lightdock/lightdock)
 [![Code Coverage](https://codecov.io/gh/lightdock/lightdock/branch/master/graph/badge.svg)](https://codecov.io/gh/lightdock/lightdock)
-[![Commit activity](https://img.shields.io/github/commit-activity/m/lightdock/lightdock.svg)](https://github.com/lightdock/lightdock/commits/master)
+[![Downloads](https://static.pepy.tech/personalized-badge/lightdock?period=total&units=international_system&left_color=black&right_color=orange&left_text=Downloads)](https://pepy.tech/project/lightdock)
 
 <p align="center">
     <img src="https://lightdock.org/assets/images/lightdock_logo.png">
@@ -26,13 +26,17 @@ LightDock protocol and the updates to make use of residue restraints have been p
 Jorge Roel-Touris, Alexandre M.J.J. Bonvin, [Brian Jiménez-García](http://bjimenezgarcia.com)<br>
 *Bioinformatics*, btz642; doi: [https://doi.org/10.1093/bioinformatics/btz642](https://doi.org/10.1093/bioinformatics/btz642)
 
+**Integrative Modeling of Membrane-associated Protein Assemblies**<br>
+Jorge Roel-Touris, [Brian Jiménez-García](https://bjimenezgarcia.com) & Alexandre M.J.J. Bonvin<br>
+*Nat Commun* **11**, 6210 (2020); doi: [https://doi.org/10.1038/s41467-020-20076-5](https://doi.org/10.1038/s41467-020-20076-5)
 
 ## 3. Installation
 
 Lightdock software is compatible and it has been tested with the followings OS:
 
 * **macOS**: El Capitan, Sierra, High Sierra, Mojave, Catalina.
 * **GNU/Linux**: Ubuntu 16+, Debian Stretch+, Scientific Linux 6+, CentOS 6+.
+* **100% compatible with [Google Colab platform](https://colab.research.google.com/)**
 
 Microsoft Windows is not officially supported, despite many parts of the protocol might be able to run. Please use it at your own risk. If you wish to contribute testing and developing LightDock for Windows, please contact us.
 
@@ -42,13 +46,12 @@ LightDock has the following dependencies:
 * NumPy (<http://www.numpy.org/>)
 * Scipy (<http://www.scipy.org/>)
 * Cython (<http://cython.org/>)
-* BioPython (<http://biopython.org>)
 * ProDy (<http://prody.csb.pitt.edu/>)
 * Freesasa (<http://freesasa.github.io/>)
 
 Optional dependencies are:
 
-* MPI4py (if you plan to use MPI support which is experimental at the moment, <http://pythonhosted.org/mpi4py/>)
+* MPI4py (if you plan to use MPI support, <http://pythonhosted.org/mpi4py/>)
 
 ### 3.2. Install LightDock
 The fastest way to install LightDock is to use `pip`:
@@ -74,8 +77,6 @@ Please make sure dependencies are already installed (via pip, package manager, e
 * numpy>=1.17.1
 * scipy>=1.3.1
 * cython>=0.29.13
-* biopython>=1.74
-* pyparsing>=2.4.2
 * prody>=1.10.11
 * freesasa>=2.0.3
 
@@ -130,9 +131,9 @@ The complete documentation about how to run the LightDock protocol and several t
 
 LightDock is being actively developed and some issues may arise or you may need extra help to run LightDock. In those cases, there are two main ways to get help:
 
-1. Read the [FAQ](https://lightdock.org/tutorials/faq) in case your problem is known
+1. Read the [FAQ](https://lightdock.org/tutorials/faq) in case your problem was already reported
 2. Open a [new issue in this repository](https://github.com/lightdock/lightdock/issues/new)
-3. Or write an email to <lightdocking@gmail.com>
+3. Or write an email to <lightdocking@gmail.com> (we will do our best to answer your questions as soon as possible)
 
 ## 7. LICENSE
 

bin/lightdock3_setup.py

@@ -157,6 +157,7 @@
             args.transmembrane,
             args.write_starting_positions,
             args.swarm_radius,
+            args.flip,
         )
         if len(starting_points_files) != args.swarms:
             args.swarms = len(starting_points_files)

lightdock/prep/poses.py

@@ -20,6 +20,7 @@
 from lightdock.prep.geometry import create_bild_file
 from lightdock.structure.residue import Residue
 from lightdock.error.lightdock_errors import LightDockWarning, MembraneSetupError
+from scipy.spatial.transform import Rotation
 
 
 def get_random_point_within_sphere(number_generator, radius):
@@ -67,7 +68,16 @@ def quaternion_from_vectors(a, b):
     return Quaternion(real_part, w[0], w[1], w[2]).normalize()
 
 
-def get_quaternion_for_restraint(rec_residue, lig_residue, tx, ty, tz, rt, lt):
+def mirror_vector(v, axis, theta=np.pi):
+    axis = normalize_vector(axis)
+    rot = Rotation.from_rotvec(theta * axis)
+    new_v = rot.apply(v)
+    return new_v
+
+
+def get_quaternion_for_restraint(
+    rec_residue, lig_residue, tx, ty, tz, rt, lt, number_generator, flip=False
+):
     """Calculates the quaternion required for orienting the ligand towards the restraint"""
     r_ca = rec_residue.get_calpha()
     # Deal with possible DNA nucleotides
@@ -78,6 +88,7 @@ def get_quaternion_for_restraint(rec_residue, lig_residue, tx, ty, tz, rt, lt):
     if not l_ca:
         l_ca = lig_residue.get_atom("P")
 
+    # Center restraints pair at origin
     rx = r_ca.x + rt[0]
     ry = r_ca.y + rt[1]
     rz = r_ca.z + rt[2]
@@ -92,6 +103,17 @@ def get_quaternion_for_restraint(rec_residue, lig_residue, tx, ty, tz, rt, lt):
 
     q = quaternion_from_vectors(a, b)
 
+    # If flip mode is activated, toss a coin
+    if flip:
+        p_flip = number_generator()
+        if p_flip > 0.5:
+            # Mirror ligand restraint
+            new_lr = q.rotate([lx, ly, lz])
+
+            v_rot = mirror_vector(np.array(new_lr), np.array([tx, ty, tz]))
+
+            q = Quaternion(0, v_rot[0], v_rot[1], v_rot[2])
+
     return q
 
 
@@ -108,6 +130,7 @@ def populate_poses(
     receptor_restraints=None,
     ligand_restraints=None,
     ligand_diameter=1.0,
+    flip=False,
 ):
     """Creates new poses around a given center and a given radius"""
     new_poses = []
@@ -135,6 +158,8 @@ def populate_poses(
         distances.sort(key=lambda tup: tup[1])
         closest_residues = [x[0] for x in distances[:10]]
 
+    # Uncomment for fixed translations:
+    # x, y, z = get_random_point_within_sphere(number_generator, radius)
     for _ in range(to_generate):
         # First calculate a random translation within the swarm sphere
         x, y, z = get_random_point_within_sphere(number_generator, radius)
@@ -154,7 +179,15 @@ def populate_poses(
             ]
             # Calculate the quaternion which rotates the ligand to point to the given receptor restraint
             q = get_quaternion_for_restraint(
-                rec_residue, lig_residue, tx, ty, tz, rec_translation, lig_translation
+                rec_residue,
+                lig_residue,
+                tx,
+                ty,
+                tz,
+                rec_translation,
+                lig_translation,
+                number_generator,
+                flip,
             )
 
         # Only restraints in the ligand partner
@@ -178,8 +211,17 @@ def populate_poses(
                 number_generator.randint(0, len(ligand_restraints) - 1)
             ]
             q = get_quaternion_for_restraint(
-                rec_residue, lig_residue, tx, ty, tz, rec_translation, lig_translation
+                rec_residue,
+                lig_residue,
+                tx,
+                ty,
+                tz,
+                rec_translation,
+                lig_translation,
+                number_generator,
+                flip,
             )
+
         # No restraints at all
         else:
             q = Quaternion.random(number_generator)
@@ -374,6 +416,7 @@ def calculate_initial_poses(
     is_transmembrane=False,
     writing_starting_positions=False,
     swarm_radius=10.0,
+    flip=False,
 ):
     """Calculates the starting points for each of the glowworms using the center of swarms"""
 
@@ -424,7 +467,6 @@ def calculate_initial_poses(
     pdb_file_name = os.path.join(dest_folder, SWARM_CENTERS_FILE)
     create_pdb_from_points(pdb_file_name, swarm_centers)
 
-    _ = ligand.center_of_coordinates()
     positions_files = []
 
     for swarm_id, swarm_center in enumerate(swarm_centers):
@@ -441,6 +483,7 @@ def calculate_initial_poses(
             receptor_restraints,
             ligand_restraints,
             ligand_diameter,
+            flip,
         )
         if writing_starting_positions:
             # Save poses as pdb file

lightdock/prep/simulation.py

@@ -169,6 +169,7 @@ def calculate_starting_positions(
     is_transmembrane=False,
     write_starting_positions=False,
     swarm_radius=DEFAULT_SWARM_RADIUS,
+    flip=False,
 ):
     """Defines the starting positions of each glowworm in the simulation.
 
@@ -198,6 +199,7 @@ def calculate_starting_positions(
             is_transmembrane,
             write_starting_positions,
             swarm_radius,
+            flip,
         )
         log.info(f"Generated {len(starting_points_files)} positions files")
     else:

lightdock/scoring/sd/energy/c/sd.c

@@ -43,8 +43,6 @@ static PyObject * sd_calculate_energy(PyObject *self, PyObject *args) {
 
         interface_cutoff2 = interface_cutoff*interface_cutoff;
 
-        descr = PyArray_DescrFromType(NPY_DOUBLE);
-
         tmp0 = PyObject_GetAttrString(receptor_coordinates, "coordinates");
         tmp1 = PyObject_GetAttrString(ligand_coordinates, "coordinates");
 
@@ -53,9 +51,11 @@ static PyObject * sd_calculate_energy(PyObject *self, PyObject *args) {
 
         dims[1] = 3;
         dims[0] = rec_len;
+        descr = PyArray_DescrFromType(NPY_DOUBLE);
         PyArray_AsCArray((PyObject **)&tmp0, (void **)&rec_array, dims, 2, descr);
 
         dims[0] = lig_len;
+        descr = PyArray_DescrFromType(NPY_DOUBLE);
         PyArray_AsCArray((PyObject **)&tmp1, (void **)&lig_array, dims, 2, descr);
 
         total_elec = 0.0;
@@ -121,15 +121,10 @@ static PyObject * sd_calculate_energy(PyObject *self, PyObject *args) {
         }
 
         // Free structures
-        Py_DECREF(rec_c_charges);
-        Py_DECREF(lig_c_charges);
-        Py_DECREF(rec_c_vdw);
-        Py_DECREF(lig_c_vdw);
-        Py_DECREF(rec_c_vdw_radii);
-        Py_DECREF(lig_c_vdw_radii);
-        Py_DECREF(descr);
         PyArray_Free(tmp0, rec_array);
         PyArray_Free(tmp1, lig_array);
+        Py_DECREF(tmp0);
+        Py_DECREF(tmp1);
     }
 
     interface_receptor = realloc(interface_receptor, interface_len*sizeof(unsigned int));

lightdock/scoring/sipper/c/sipper.c

@@ -41,8 +41,6 @@ static PyObject * calculate_sipper(PyObject *self, PyObject *args) {
 
         interface_cutoff2 = interface_cutoff*interface_cutoff;
 
-        descr = PyArray_DescrFromType(NPY_DOUBLE);
-
         tmp0 = PyObject_GetAttrString(receptor_coordinates, "coordinates");
         tmp1 = PyObject_GetAttrString(ligand_coordinates, "coordinates");
 
@@ -51,9 +49,11 @@ static PyObject * calculate_sipper(PyObject *self, PyObject *args) {
 
         dims[1] = 3;
         dims[0] = rec_len;
+        descr = PyArray_DescrFromType(NPY_DOUBLE);
         PyArray_AsCArray((PyObject **)&tmp0, (void **)&rec_array, dims, 2, descr);
 
         dims[0] = lig_len;
+        descr = PyArray_DescrFromType(NPY_DOUBLE);
         PyArray_AsCArray((PyObject **)&tmp1, (void **)&lig_array, dims, 2, descr);
 
         // Get pointers to the Python array structures

lightdock/scoring/vdw/energy/c/cvdw.c

@@ -39,8 +39,6 @@ static PyObject * calculate_vdw(PyObject *self, PyObject *args) {
 
         interface_cutoff2 = interface_cutoff*interface_cutoff;
 
-        descr = PyArray_DescrFromType(NPY_DOUBLE);
-
         tmp0 = PyObject_GetAttrString(receptor_coordinates, "coordinates");
         tmp1 = PyObject_GetAttrString(ligand_coordinates, "coordinates");
 
@@ -49,9 +47,11 @@ static PyObject * calculate_vdw(PyObject *self, PyObject *args) {
 
         dims[1] = 3;
         dims[0] = rec_len;
+        descr = PyArray_DescrFromType(NPY_DOUBLE);
         PyArray_AsCArray((PyObject **)&tmp0, (void **)&rec_array, dims, 2, descr);
 
         dims[0] = lig_len;
+        descr = PyArray_DescrFromType(NPY_DOUBLE);
         PyArray_AsCArray((PyObject **)&tmp1, (void **)&lig_array, dims, 2, descr);
 
         // Get pointers to the Python array structures
@@ -99,6 +99,8 @@ static PyObject * calculate_vdw(PyObject *self, PyObject *args) {
         // Free structures
         PyArray_Free(tmp0, rec_array);
         PyArray_Free(tmp1, lig_array);
+        Py_DECREF(tmp0);
+        Py_DECREF(tmp1);
     }
 
     // Return a tuple with the following values for calculated energies:

lightdock/structure/complex.py

@@ -94,7 +94,8 @@ def get_nm_mask(self):
         """Calculates the mask on atoms to apply ANM"""
         mask = []
         for residue in self.residues:
-            if residue.is_standard() or residue.is_nucleic():
+            # is_protein should accept common modified cases supported by ProDy
+            if residue.is_protein() or residue.is_nucleic():
                 mask.extend([True] * len(residue.atoms))
             else:
                 mask.extend([False] * len(residue.atoms))

lightdock/structure/nm.py

@@ -21,7 +21,7 @@ def calculate_nmodes(pdb_file_name, n_modes, molecule):
     backbone_atoms = prody_molecule.select("name CA")
     if not backbone_atoms:
         # If previous step has failed, maybe we're dealing with DNA
-        backbone_atoms = prody_molecule.select("nucleic and name C4'")
+        backbone_atoms = prody_molecule.select("name C4'")
     if not backbone_atoms:
         raise NormalModesCalculationError(
             "Error selecting backbone atoms (protein or nucleic)"

lightdock/test/bin/golden_data/regression_dfire2_long/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 50,
     "ligand_pdb": "2UUY_lig.pdb",
     "membrane": false,

lightdock/test/bin/golden_data/regression_dfire2_restraints/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 25,
     "ligand_pdb": "2UUY_lig.pdb",
     "ligand_restraints": {

lightdock/test/bin/golden_data/regression_dfire2_short/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 25,
     "ligand_pdb": "2UUY_lig.pdb",
     "membrane": false,

lightdock/test/bin/golden_data/regression_dfire_long/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 20,
     "ligand_pdb": "1PPE_lig.pdb",
     "membrane": false,

lightdock/test/bin/golden_data/regression_dfire_short/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 25,
     "ligand_pdb": "2UUY_lig.pdb",
     "membrane": false,

lightdock/test/bin/golden_data/regression_dna_restraints/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 50,
     "ligand_pdb": "1DIZ_lig.pdb.H",
     "ligand_restraints": {

lightdock/test/bin/golden_data/regression_dna_short/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 25,
     "ligand_pdb": "1DIZ_lig.pdb.H",
     "membrane": false,

lightdock/test/bin/golden_data/regression_fastdfire_long/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 50,
     "ligand_pdb": "2UUY_lig.pdb",
     "membrane": false,

lightdock/test/bin/golden_data/regression_fastdfire_min/setup.json

@@ -2,6 +2,7 @@
     "anm_lig": 10,
     "anm_rec": 10,
     "anm_seed": 324324,
+    "flip": false,
     "glowworms": 200,
     "ligand_pdb": "2UUY_lig.pdb",
     "membrane": false,