⏪ revert changes in rdkit_grid_featurizer

# flake8: noqa

import logging

import time

import hashlib

from collections import Counter

from deepchem.utils.rdkit_utils import MoleculeLoadException, load_molecule

import numpy as np

from scipy.spatial.distance import cdist

from copy import deepcopy

from deepchem.feat.base_classes import ComplexFeaturizer

from deepchem.utils.rdkit_utils import MoleculeLoadException, load_molecule

from deepchem.utils.geometry_utils import angle_between

from deepchem.utils.geometry_utils import compute_centroid, subtract_centroid

from deepchem.utils.geometry_utils import generate_random_rotation_matrix

from deepchem.utils.geometry_utils import compute_pairwise_distances

from deepchem.utils.hash_utils import hash_ecfp, hash_ecfp_pair

from deepchem.utils.voxel_utils import convert_atom_to_voxel

from deepchem.utils.voxel_utils import convert_atom_pair_to_voxel

from deepchem.feat import ComplexFeaturizer

logger = logging.getLogger(__name__)

def compute_centroid(coordinates):

  """Compute the x,y,z centroid of provided coordinates

  coordinates: np.ndarray

    Shape (N, 3), where N is number atoms.

  """

  centroid = np.mean(coordinates, axis=0)

  return (centroid)

def generate_random__unit_vector():

  """Generate a random unit vector on the 3-sphere.

  citation:

  http://mathworld.wolfram.com/SpherePointPicking.html

  a. Choose random theta \element [0, 2*pi]

  b. Choose random z \element [-1, 1]

  c. Compute output vector u: (x,y,z) = (sqrt(1-z^2)*cos(theta), sqrt(1-z^2)*sin(theta),z)

  """

  theta = np.random.uniform(low=0.0, high=2 * np.pi)

  z = np.random.uniform(low=-1.0, high=1.0)

  u = np.array(

      [np.sqrt(1 - z**2) * np.cos(theta),

       np.sqrt(1 - z**2) * np.sin(theta), z])

  return (u)

def generate_random_rotation_matrix():

  """Generate a random rotation matrix in 3D.

  1. Generate a random unit vector u, randomly sampled from the unit

     3-sphere (see function generate_random__unit_vector() for details)

  2. Generate a second random unit vector v

    a. If absolute value of u \dot v > 0.99, repeat.

       (This is important for numerical stability. Intuition: we want them to

       be as linearly independent as possible or else the orthogonalized

       version of v will be much shorter in magnitude compared to u. I assume

       in Stack they took this from Gram-Schmidt orthogonalization?)

    b. v" = v - (u \dot v)*u, i.e. subtract out the component of v that's in

       u's direction

    c. normalize v" (this isn"t in Stack but I assume it must be done)

  3. find w = u \cross v"

  4. u, v", and w will form the columns of a rotation matrix, R. The

     intuition is that u, v" and w are, respectively, what the standard basis

     vectors e1, e2, and e3 will be mapped to under the transformation.

  """

  u = generate_random__unit_vector()

  v = generate_random__unit_vector()

  while np.abs(np.dot(u, v)) >= 0.99:

    v = generate_random__unit_vector()

  vp = v - (np.dot(u, v) * u)

  vp /= np.linalg.norm(vp)

  w = np.cross(u, vp)

  R = np.column_stack((u, vp, w))

  return (R)

def rotate_molecules(mol_coordinates_list):

  """Rotates provided molecular coordinates.

  return (rotated_coordinates_list)

def compute_pairwise_distances(protein_xyz, ligand_xyz):

  """Takes an input m x 3 and n x 3 np arrays of 3D coords of protein and ligand,

  respectively, and outputs an m x n np array of pairwise distances in Angstroms

  between protein and ligand atoms. entry (i,j) is dist between the i"th protein

  atom and the j"th ligand atom.

  """

  pairwise_distances = cdist(protein_xyz, ligand_xyz, metric='euclidean')

  return (pairwise_distances)

"""following two functions adapted from:

http://stackoverflow.com/questions/2827393/angles-between-two-n-dimensional-vectors-in-python

"""

def unit_vector(vector):

  """ Returns the unit vector of the vector.  """

  return vector / np.linalg.norm(vector)

def angle_between(vector_i, vector_j):

  """Returns the angle in radians between vectors "vector_i" and "vector_j"::

  >>> print("%0.06f" % angle_between((1, 0, 0), (0, 1, 0)))

  1.570796

  >>> print("%0.06f" % angle_between((1, 0, 0), (1, 0, 0)))

  0.000000

  >>> print("%0.06f" % angle_between((1, 0, 0), (-1, 0, 0)))

  3.141593

  Note that this function always returns the smaller of the two angles between

  the vectors (value between 0 and pi).

  """

  vector_i_u = unit_vector(vector_i)

  vector_j_u = unit_vector(vector_j)

  angle = np.arccos(np.dot(vector_i_u, vector_j_u))

  if np.isnan(angle):

    if np.allclose(vector_i_u, vector_j_u):

      return 0.0

    else:

      return np.pi

  return angle

def hash_sybyl(sybyl, sybyl_types):

  return (sybyl_types.index(sybyl))

def hash_ecfp(ecfp, power):

  """

  Returns an int of size 2^power representing that

  ECFP fragment. Input must be a string.

  """

  ecfp = ecfp.encode('utf-8')

  md5 = hashlib.md5()

  md5.update(ecfp)

  digest = md5.hexdigest()

  ecfp_hash = int(digest, 16) % (2**power)

  return (ecfp_hash)

def hash_ecfp_pair(ecfp_pair, power):

  """Returns an int of size 2^power representing that ECFP pair. Input must be

  a tuple of strings.

  """

  ecfp = "%s,%s" % (ecfp_pair[0], ecfp_pair[1])

  ecfp = ecfp.encode('utf-8')

  md5 = hashlib.md5()

  md5.update(ecfp)

  digest = md5.hexdigest()

  ecfp_hash = int(digest, 16) % (2**power)

  return (ecfp_hash)

def compute_all_ecfp(mol, indices=None, degree=2):

  """Obtain molecular fragment for all atoms emanating outward to given degree.

  For each fragment, compute SMILES string (for now) and hash to an int.

  cutoff: float

    Cutoff distance for contact consideration.

  """

  features_dict = {}

  if pairwise_distances is None:

    pairwise_distances = compute_pairwise_distances(protein_xyz, ligand_xyz)

  contacts = np.nonzero((pairwise_distances < cutoff))

  return salt_bridge_contacts

def is_angle_within_cutoff(vector_i, vector_j, hbond_angle_cutoff):

  angle = angle_between(vector_i, vector_j) * 180. / np.pi

  return (angle > (180 - hbond_angle_cutoff) and

          angle < (180. + hbond_angle_cutoff))

def is_hydrogen_bond(protein_xyz, protein, ligand_xyz, ligand, contact,

                     hbond_angle_cutoff):

  """

  return (hbond_contacts)

def convert_atom_to_voxel(molecule_xyz,

                          atom_index,

                          box_width,

                          voxel_width,

                          verbose=False):

  """Converts atom coordinates to an i,j,k grid index.

  Parameters

  ----------

  molecule_xyz: np.ndarray

    Array with coordinates of all atoms in the molecule, shape (N, 3)

  atom_index: int

    Index of an atom

  box_width: float

    Size of a box

  voxel_width: float

    Size of a voxel

  verbose: bool

    Print warnings when atom is outside of a box

  """

  indices = np.floor(

      (molecule_xyz[atom_index] + box_width / 2.0) / voxel_width).astype(int)

  if ((indices < 0) | (indices >= box_width / voxel_width)).any():

    if verbose:

      logger.warning('Coordinates are outside of the box (atom id = %s,'

                     ' coords xyz = %s, coords in box = %s' %

                     (atom_index, molecule_xyz[atom_index], indices))

  return ([indices])

def convert_atom_pair_to_voxel(molecule_xyz_tuple, atom_index_pair, box_width,

                               voxel_width):

  """Converts a pair of atoms to a list of i,j,k tuples."""

  indices_list = []

  indices_list.append(

      convert_atom_to_voxel(molecule_xyz_tuple[0], atom_index_pair[0],

                            box_width, voxel_width)[0])

  indices_list.append(

      convert_atom_to_voxel(molecule_xyz_tuple[1], atom_index_pair[1],

                            box_width, voxel_width)[0])

  return (indices_list)

def compute_charge_dictionary(molecule):

  """Create a dictionary with partial charges for each atom in the molecule.

  return charge_dictionary

def subtract_centroid(xyz, centroid):

  """Subtracts centroid from each coordinate.

  Subtracts the centroid, a numpy array of dim 3, from all coordinates of all

  atoms in the molecule

  """

  xyz -= np.transpose(centroid)

  return (xyz)

class RdkitGridFeaturizer(ComplexFeaturizer):

  """Featurizes protein-ligand complex using flat features or a 3D grid (in which

  each voxel is described with a vector of features).

      Filename for protein pdb file. 

    """

    try:

      # TIMING

      ############################################################## TIMING

      time1 = time.time()

      # TIMING

      ############################################################## TIMING

      protein_xyz, protein_rdk = load_molecule(

          protein_pdb_file, calc_charges=True, sanitize=self.sanitize)

      # TIMING

      ############################################################## TIMING

      time2 = time.time()

      logger.info(

          "TIMING: Loading protein coordinates took %0.3f s" % (time2 - time1),

          self.verbose)

      # TIMING

      # TIMING

      ############################################################## TIMING

      ############################################################## TIMING

      time1 = time.time()

      # TIMING

      ############################################################## TIMING

      ligand_xyz, ligand_rdk = load_molecule(

          mol_pdb_file, calc_charges=True, sanitize=self.sanitize)

      # TIMING

      ############################################################## TIMING

      time2 = time.time()

      logger.info(

          "TIMING: Loading ligand coordinates took %0.3f s" % (time2 - time1),

          self.verbose)

      # TIMING

      ############################################################## TIMING

    except MoleculeLoadException:

      logger.warning("Some molecules cannot be loaded by Rdkit. Skipping")

      return None

    # TIMING

    ############################################################## TIMING

    time1 = time.time()

    # TIMING

    ############################################################## TIMING

    centroid = compute_centroid(ligand_xyz)

    ligand_xyz = subtract_centroid(ligand_xyz, centroid)

    protein_xyz = subtract_centroid(protein_xyz, centroid)

    # TIMING

    ############################################################## TIMING

    time2 = time.time()

    logger.info("TIMING: Centroid processing took %0.3f s" % (time2 - time1),

                self.verbose)

    # TIMING

    ############################################################## TIMING

    pairwise_distances = compute_pairwise_distances(protein_xyz, ligand_xyz)

        self.assertIsInstance(mol_rdk, Mol)

        self.assertEqual(mol_xyz.shape, (num_atoms, 3))

  def test_generate_random_unit_vector(self):

    for _ in range(100):

      u = rgf.generate_random__unit_vector()

      # 3D vector with unit length

      self.assertEqual(u.shape, (3,))

      self.assertAlmostEqual(np.linalg.norm(u), 1.0)

  def test_generate_random_rotation_matrix(self):

    # very basic test, we check if rotations actually work in test_rotate_molecules

    for _ in range(100):

    distance = rgf.compute_pairwise_distances(coords1, coords2)

    self.assertTrue((distance == [[1, 2, 3], [0, 1, 2]]).all())

  def test_unit_vector(self):

    for _ in range(10):

      vector = np.random.rand(3)

      norm_vector = rgf.unit_vector(vector)

      self.assertAlmostEqual(np.linalg.norm(norm_vector), 1.0)

  def test_angle_between(self):

    for _ in range(10):

      v1 = np.random.rand(3,)

      for box_width in (10, 20, 40):

        for voxel_width in (0.5, 1, 2):

          voxel = rgf.convert_atom_to_voxel(xyz, idx, box_width, voxel_width)

          self.assertIsInstance(voxel, np.ndarray)

          self.assertEqual(voxel.shape, (3,))

          self.assertIs(voxel.dtype, np.dtype('int'))

          self.assertIsInstance(voxel, list)

          self.assertEqual(len(voxel), 1)

          self.assertIsInstance(voxel[0], np.ndarray)

          self.assertEqual(voxel[0].shape, (3,))

          self.assertIs(voxel[0].dtype, np.dtype('int'))

          # indices are positive

          self.assertTrue((voxel >= 0).all())

          self.assertTrue((voxel[0] >= 0).all())

          # coordinates were properly translated and scaled

          self.assertTrue(

              (voxel < (box_width + coords_range) / 2.0 / voxel_width).all())

              (voxel[0] < (box_width + coords_range) / 2.0 / voxel_width).all())

          self.assertTrue(

              np.allclose(voxel,

              np.allclose(voxel[0],

                          np.floor((xyz[idx] + box_width / 2.0) / voxel_width)))

    # for coordinates outside of the box function should properly transform them

      # but it is not implemented in 2.7 and buggy in 3.5 (issue 29620)

      voxel = rgf.convert_atom_to_voxel(*args)

      self.assertTrue(

          np.allclose(voxel, np.floor((args[0] + args[2] / 2.0) / args[3])))

          np.allclose(voxel[0], np.floor((args[0] + args[2] / 2.0) / args[3])))

  def test_convert_atom_pair_to_voxel(self):

    # 20 points with coords between -5 and 5, centered at 0

  def test_featurize_binding_pocket_ecfp(self):

    prot_xyz, prot_rdk = rgf.load_molecule(self.protein_file)

    lig_xyz, lig_rdk = rgf.load_molecule(self.ligand_file)

    distance = rgf.compute_pairwise_distances(prot_xyz, lig_xyz)

    distance = rgf.compute_pairwise_distances(

        protein_xyz=prot_xyz, ligand_xyz=lig_xyz)

    # check if results are the same if we provide precomputed distances

    prot_dict, lig_dict = rgf.featurize_binding_pocket_ecfp(

    lig_xyz, lig_rdk = rgf.load_molecule(self.ligand_file)

    prot_num_atoms = prot_rdk.GetNumAtoms()

    lig_num_atoms = lig_rdk.GetNumAtoms()

    distance = rgf.compute_pairwise_distances(prot_xyz, lig_xyz)

    distance = rgf.compute_pairwise_distances(

        protein_xyz=prot_xyz, ligand_xyz=lig_xyz)

    for bins in ((0, 2), (2, 3)):

      splif_dict = rgf.compute_splif_features_in_range(

  def test_featurize_splif(self):

    prot_xyz, prot_rdk = rgf.load_molecule(self.protein_file)

    lig_xyz, lig_rdk = rgf.load_molecule(self.ligand_file)

    distance = rgf.compute_pairwise_distances(prot_xyz, lig_xyz)

    distance = rgf.compute_pairwise_distances(

        protein_xyz=prot_xyz, ligand_xyz=lig_xyz)

    bins = [(1, 2), (2, 3)]