♻️ create complex featurizers directory

from deepchem.feat.graph_features import ConvMolFeaturizer

from deepchem.feat.graph_features import WeaveFeaturizer

from deepchem.feat.rdkit_grid_featurizer import RdkitGridFeaturizer

from deepchem.feat.binding_pocket_features import BindingPocketFeaturizer

from deepchem.feat.atomic_coordinates import AtomicCoordinates

from deepchem.feat.atomic_coordinates import NeighborListComplexAtomicCoordinates

# molecule featurizers

from deepchem.feat.molecule_featurizers import MolGraphConvFeaturizer

from deepchem.feat.molecule_featurizers import AtomicCoordinates

from deepchem.feat.molecule_featurizers import CircularFingerprint

from deepchem.feat.molecule_featurizers import CoulombMatrix

from deepchem.feat.molecule_featurizers import CoulombMatrixEig

from deepchem.feat.molecule_featurizers import MordredDescriptors

from deepchem.feat.molecule_featurizers import Mol2VecFingerprint

from deepchem.feat.molecule_featurizers import MolGraphConvFeaturizer

from deepchem.feat.molecule_featurizers import OneHotFeaturizer

from deepchem.feat.molecule_featurizers import RawFeaturizer

from deepchem.feat.molecule_featurizers import RDKitDescriptors

from deepchem.feat.molecule_featurizers import SmilesToImage

from deepchem.feat.molecule_featurizers import SmilesToSeq, create_char_to_idx

# complex featurizers

from deepchem.feat.complex_featurizers import RdkitGridFeaturizer

from deepchem.feat.complex_featurizers import NeighborListAtomicCoordinates

from deepchem.feat.complex_featurizers import NeighborListComplexAtomicCoordinates

from deepchem.feat.complex_featurizers import ComplexNeighborListFragmentAtomicCoordinates

# material featurizers

from deepchem.feat.material_featurizers import ElementPropertyFingerprint

from deepchem.feat.material_featurizers import SineCoulombMatrix

"""

Featurizers for complex.

"""

# flake8: noqa

from deepchem.feat.complex_featurizers.rdkit_grid_featurizer import RdkitGridFeaturizer

from deepchem.feat.complex_featurizers.complex_atomic_coordinates import NeighborListAtomicCoordinates

from deepchem.feat.complex_featurizers.complex_atomic_coordinates import NeighborListComplexAtomicCoordinates

from deepchem.feat.complex_featurizers.complex_atomic_coordinates import ComplexNeighborListFragmentAtomicCoordinates

Atomic coordinate featurizer.

"""

import logging

import numpy as np

from deepchem.feat import Featurizer

from deepchem.feat import ComplexFeaturizer

from deepchem.feat.base_classes import Featurizer, ComplexFeaturizer

from deepchem.feat.molecule_featurizers import AtomicCoordinates

from deepchem.utils.data_utils import pad_array

from deepchem.utils.rdkit_utils import MoleculeLoadException, get_xyz_from_mol, \

  load_molecule, merge_molecules_xyz, merge_molecules

class AtomicCoordinates(Featurizer):

  """

  Nx3 matrix of Cartesian coordinates [Angstrom]

  """

  name = ['atomic_coordinates']

  def _featurize(self, mol):

    """

    Calculate atomic coodinates.

    Parameters

    ----------

    mol : RDKit Mol

          Molecule.

    """

    N = mol.GetNumAtoms()

    coords = np.zeros((N, 3))

    # RDKit stores atomic coordinates in Angstrom. Atomic unit of length is the

    # bohr (1 bohr = 0.529177 Angstrom). Converting units makes gradient calculation

    # consistent with most QM software packages.

    coords_in_bohr = [

        mol.GetConformer(0).GetAtomPosition(i).__idiv__(0.52917721092)

        for i in range(N)

    ]

    for atom in range(N):

      coords[atom, 0] = coords_in_bohr[atom].x

      coords[atom, 1] = coords_in_bohr[atom].y

      coords[atom, 2] = coords_in_bohr[atom].z

    coords = [coords]

    return coords

def compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors,

                          periodic_box_size):

  """Computes a neighbor list from atom coordinates."""

  return neighbor_list

def get_coords(mol):

  """

  Gets coordinates in Angstrom for RDKit mol.

  """

  N = mol.GetNumAtoms()

  coords = np.zeros((N, 3))

  coords_raw = [mol.GetConformer(0).GetAtomPosition(i) for i in range(N)]

  for atom in range(N):

    coords[atom, 0] = coords_raw[atom].x

    coords[atom, 1] = coords_raw[atom].y

    coords[atom, 2] = coords_raw[atom].z

  return coords

class NeighborListAtomicCoordinates(Featurizer):

  """

  Adjacency List of neighbors in 3-space

    self.periodic_box_size = periodic_box_size

    # Type of data created by this featurizer

    self.dtype = object

    self.coordinates_featurizer = AtomicCoordinates()

    self.bohr_coords_featurizer = AtomicCoordinates(use_bohr=True)

    self.coords_featurizer = AtomicCoordinates(use_bohr=False)

  def _featurize(self, mol):

    """

        To be featurized.

    """

    # TODO(rbharath): Should this return a list?

    bohr_coords = self.coordinates_featurizer._featurize(mol)[0]

    coords = get_coords(mol)

    bohr_coords = self.bohr_coords_featurizer._featurize(mol)

    coords = self.coords_featurizer._featurize(mol)

    neighbor_list = compute_neighbor_list(coords, self.neighbor_cutoff,

                                          self.max_num_neighbors,

                                          self.periodic_box_size)

    self.neighbor_cutoff = neighbor_cutoff

    # Type of data created by this featurizer

    self.dtype = object

    self.coordinates_featurizer = AtomicCoordinates()

  def _featurize(self, mol_pdb_file, protein_pdb_file):

    """

# 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 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)

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

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

  return (R)

logger = logging.getLogger(__name__)

def rotate_molecules(mol_coordinates_list):

  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)

import numpy as np

import deepchem as dc

from deepchem.feat.base_classes import MolecularFeaturizer

from deepchem.feat.atomic_coordinates import ComplexNeighborListFragmentAtomicCoordinates

from deepchem.feat.complex_featurizers import ComplexNeighborListFragmentAtomicCoordinates

from deepchem.feat.mol_graphs import ConvMol, WeaveMol

from deepchem.data import DiskDataset

import logging