Commit cafb2dca authored by Bharath Ramsundar's avatar Bharath Ramsundar
Browse files

Continuing factoring

parent d6c0f61c
Loading
Loading
Loading
Loading
+78 −0
Original line number Diff line number Diff line
@@ -6,6 +6,8 @@ import logging
from deepchem.utils.hash_utils import hash_ecfp
from deepchem.feat import ComplexFeaturizer
from deepchem.utils.hash_utils import vectorize
from deepchem.utils.voxel_utils import voxelize 
from deepchem.utils.voxel_utils import convert_atom_to_voxel
from deepchem.utils.rdkit_util import compute_all_ecfp
from deepchem.utils.rdkit_util import compute_pairwise_distances

@@ -107,3 +109,79 @@ class ContactCircularFingerprint(ComplexFeaturizer):
            cutoff=self.cutoff,
            ecfp_degree=self.radius)
    ]

class ContactCircularVoxelizer(ComplexFeaturizer):
  """Computes ECFP fingerprints on a voxel grid.

  Given a macromolecular complex made up of multiple
  constituent molecules, first compute the contact points where
  atoms from different molecules come close to one another. For
  atoms within "contact regions," compute radial "ECFP"
  fragments which are sub-molecules centered at atoms in the
  contact region. Localize these ECFP fingeprints at the voxel
  in which they originated.
  """

  def __init__(self, 
               cutoff=4.5,
               radius=2,
               size=8,
               box_width=16.0,
               voxel_width=1.0):
    """
    Parameters
    ----------
    cutoff: float (default 4.5)
      Distance cutoff in angstroms for molecules in complex.
    radius : int, optional (default 2)
        Fingerprint radius.
    size : int, optional (default 8)
      Length of generated bit vector.
    box_width: float, optional (default 16.0)
      Size of a box in which voxel features are calculated. Box
      is centered on a ligand centroid.
    voxel_width: float, optional (default 1.0)
      Size of a 3D voxel in a grid.
    """
    self.cutoff = cutoff
    self.radius = radius
    self.size = size
    self.box_width = box_width
    self.voxel_width = voxel_width
    self.voxels_per_edge = int(self.box_width / self.voxel_width)

  def _featurize_complex(self, mol, protein):
    """
    Compute featurization for a single mol/protein complex

    Parameters
    ----------
    mol: object
      Representation of the molecule
    protein: object
      Representation of the protein
    """
    (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein
    distances = compute_pairwise_distances(prot_xyz, lig_xyz)
    return [
        sum([
            voxelize(
                convert_atom_to_voxel,
                self.voxels_per_edge,
                self.box_width,
                self.voxel_width,
                hash_ecfp,
                xyz,
                feature_dict=ecfp_dict,
                nb_channel=self.size)
            for xyz, ecfp_dict in zip((prot_xyz, lig_xyz),
                                      featurize_binding_pocket_ecfp(
                                          prot_xyz,
                                          prot_rdk,
                                          lig_xyz,
                                          lig_rdk,
                                          distances,
                                          cutoff=self.cutoff,
                                          ecfp_degree=self.radius))
        ])
    ]
+81 −0
Original line number Diff line number Diff line
"""
Compute various spatial fingerprints for macromolecular complexes.
"""
import logging
from deepchem.utils.rdkit_util import get_partial_charge
from deepchem.feat import ComplexFeaturizer
from deepchem.utils.voxel_utils import voxelize 

logger = logging.getLogger(__name__)

def compute_charge_dictionary(molecule):
  """Create a dictionary with partial charges for each atom in the molecule.

  This function assumes that the charges for the molecule are
  already computed (it can be done with
  rdkit_util.compute_charges(molecule))
  """

  charge_dictionary = {}
  for i, atom in enumerate(molecule.GetAtoms()):
    charge_dictionary[i] = get_partial_charge(atom)
  return charge_dictionary

class ChargeVoxelizer(ComplexFeaturizer):
  """Localize partial charges of atoms in macromolecular complexes.

  Given a macromolecular compelx made up of multiple
  constitutent molecules, compute the partial (Gasteiger
  charge) on each molecule. For each atom, localize this
  partial charge in the voxel in which it originated to create
  a local charge array. Sum contributes to get an effective
  charge at each voxel.

  Creates a tensor output of shape (box_width, box_width,
  box_width, 1) for each macromolecular complex that computes
  the effective charge at each voxel.
  """
  def __init__(self, 
               box_width=16.0,
               voxel_width=1.0):
    """
    Parameters
    ----------
    box_width: float, optional (default 16.0)
      Size of a box in which voxel features are calculated. Box
      is centered on a ligand centroid.
    voxel_width: float, optional (default 1.0)
      Size of a 3D voxel in a grid.
    """
    self.cutoff = cutoff
    self.box_width = box_width
    self.voxel_width = voxel_width
    self.voxels_per_edge = int(self.box_width / self.voxel_width)

  def _featurize_complex(self, mol, protein):
    """
    Compute featurization for a single mol/protein complex

    Parameters
    ----------
    mol: object
      Representation of the molecule
    protein: object
      Representation of the protein
    """
    (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein
    return [
        sum([
            voxelize(
                convert_atom_to_voxel,
                self.voxels_per_edge,
                self.box_width,
                self.voxel_width,
                None,
                xyz,
                feature_dict=compute_charge_dictionary(mol),
                nb_channel=1,
                dtype="np.float16")
            for xyz, mol in ((prot_xyz, prot_rdk), (lig_xyz, lig_rdk))
        ])
    ]
+21 −70
Original line number Diff line number Diff line
@@ -13,15 +13,14 @@ from warnings import warn
from collections import Counter
from deepchem.feat.fingerprints import CircularFingerprint
from deepchem.feat.contact_fingerprints import ContactCircularFingerprint
from deepchem.feat.contact_fingerprints import featurize_binding_pocket_ecfp
from deepchem.feat.splif_fingerprint import featurize_splif
from deepchem.feat.contact_fingerprints import ContactCircularVoxelizer
from deepchem.feat.splif_fingerprint import SplifFingerprint
from deepchem.feat.splif_fingerprint import SplifVoxelizer
from deepchem.utils.rdkit_util import load_molecule
from deepchem.utils.rdkit_util import compute_centroid
from deepchem.utils.rdkit_util import subtract_centroid
from deepchem.utils.rdkit_util import compute_ring_center
from deepchem.utils.rdkit_util import rotate_molecules
from deepchem.utils.rdkit_util import get_partial_charge
from deepchem.utils.rdkit_util import compute_pairwise_distances
from deepchem.utils.rdkit_util import is_salt_bridge
from deepchem.utils.rdkit_util import compute_salt_bridges
@@ -76,6 +75,7 @@ def featurize_binding_pocket_sybyl(protein_xyz,
                                   cutoff=7.0):
  """Computes Sybyl dicts for ligand and binding pocket of the protein.

  TODO(rbharath): Consider this comment on rdkit forums https://github.com/rdkit/rdkit/issues/1590
  Parameters
  ----------
  protein_xyz: np.ndarray
@@ -142,30 +142,6 @@ def compute_hydrogen_bonds(protein_xyz, protein, ligand_xyz, ligand,
  return (hbond_contacts)


def compute_charge_dictionary(molecule):
  """Create a dictionary with partial charges for each atom in the molecule.

  This function assumes that the charges for the molecule are already
  computed (it can be done with rdkit_util.compute_charges(molecule))
  """

  charge_dictionary = {}
  for i, atom in enumerate(molecule.GetAtoms()):
    charge_dictionary[i] = get_partial_charge(atom)
  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).
@@ -292,10 +268,22 @@ class RdkitGridFeaturizer(ComplexFeaturizer):
        cutoff=self.cutoffs['ecfp_cutoff'],
        radius=self.ecfp_degree,
        size=2**self.ecfp_power)
    self.contact_voxelizer = ContactCircularVoxelizer(
        cutoff=self.cutoffs['ecfp_cutoff'],
        radius=self.ecfp_degree,
        size=2**self.ecfp_power,
        box_width=self.box_width,
        voxel_width=self.voxel_width)
    self.splif_featurizer = SplifFingerprint(
        contact_bins=self.cutoffs['splif_contact_bins'],
        radius=self.ecfp_degree,
        size=2**self.splif_power)
    self.splif_voxelizer = SplifVoxelizer(
        contact_bins=self.cutoffs['splif_contact_bins'],
        radius=self.ecfp_degree,
        size=2**self.splif_power,
        box_width=self.box_width,
        voxel_width=self.voxel_width)

    # parse provided feature types
    for feature_type in feature_types:
@@ -351,59 +339,21 @@ class RdkitGridFeaturizer(ComplexFeaturizer):
      return self.contact_featurizer._featurize_complex((lig_xyz, lig_rdk),
                                                        (prot_xyz, prot_rdk))
    if feature_name == 'splif_hashed':
      #return [
      #    vectorize(
      #        hash_ecfp_pair,
      #        feature_dict=splif_dict,
      #        size=2**self.splif_power) for splif_dict in featurize_splif(
      #            prot_xyz, prot_rdk, lig_xyz, lig_rdk, self.cutoffs[
      #                'splif_contact_bins'], distances, self.ecfp_degree)
      #]
      return self.splif_featurizer._featurize_complex((lig_xyz, lig_rdk),
                                                      (prot_xyz, prot_rdk))
    if feature_name == 'hbond_count':
      return [
          vectorize(hash_ecfp_pair, feature_list=hbond_list, size=2**0)
          vectorize(hash_ecfp_pair, feature_list=hbond_list, size=1)
          for hbond_list in compute_hydrogen_bonds(
              prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances, self.cutoffs[
                  'hbond_dist_bins'], self.cutoffs['hbond_angle_cutoffs'])
      ]
    if feature_name == 'ecfp':
      return [
          sum([
              voxelize(
                  convert_atom_to_voxel,
                  self.voxels_per_edge,
                  self.box_width,
                  self.voxel_width,
                  hash_ecfp,
                  xyz,
                  feature_dict=ecfp_dict,
                  nb_channel=2**self.ecfp_power)
              for xyz, ecfp_dict in zip((prot_xyz, lig_xyz),
                                        featurize_binding_pocket_ecfp(
                                            prot_xyz,
                                            prot_rdk,
                                            lig_xyz,
                                            lig_rdk,
                                            distances,
                                            cutoff=self.cutoffs['ecfp_cutoff'],
                                            ecfp_degree=self.ecfp_degree))
          ])
      ]
      return self.contact_voxelizer._featurize_complex((lig_xyz, lig_rdk),
                                                       (prot_xyz, prot_rdk))
    if feature_name == 'splif':
      return [
          voxelize(
              convert_atom_pair_to_voxel,
              self.voxels_per_edge,
              self.box_width,
              self.voxel_width,
              hash_ecfp_pair, (prot_xyz, lig_xyz),
              feature_dict=splif_dict,
              nb_channel=2**self.splif_power) for splif_dict in featurize_splif(
                  prot_xyz, prot_rdk, lig_xyz, lig_rdk, self.cutoffs[
                      'splif_contact_bins'], distances, self.ecfp_degree)
      ]
      return self.splif_voxelizer._featurize_complex((lig_xyz, lig_rdk),
                                                     (prot_xyz, prot_rdk))
    if feature_name == 'sybyl':
      return [
          voxelize(
@@ -535,6 +485,7 @@ class RdkitGridFeaturizer(ComplexFeaturizer):
      return None

    time1 = time.time()
    # TODO(rbharath): Wait a minute, if we're subtract different centroids for each, have we ruined the alignment between protein and ligand. Figure this out before merging in.
    centroid = compute_centroid(ligand_xyz)
    ligand_xyz = subtract_centroid(ligand_xyz, centroid)
    protein_xyz = subtract_centroid(protein_xyz, centroid)
+85 −0
Original line number Diff line number Diff line
@@ -7,6 +7,9 @@ from deepchem.utils.hash_utils import hash_ecfp_pair
from deepchem.utils.rdkit_util import compute_all_ecfp
from deepchem.feat import ComplexFeaturizer
from deepchem.utils.hash_utils import vectorize
from deepchem.utils.voxel_utils import voxelize 
from deepchem.utils.voxel_utils import convert_atom_to_voxel
from deepchem.utils.voxel_utils import convert_atom_pair_to_voxel
from deepchem.utils.rdkit_util import compute_pairwise_distances

logger = logging.getLogger(__name__)
@@ -126,3 +129,85 @@ class SplifFingerprint(ComplexFeaturizer):
            size=self.size) for splif_dict in featurize_splif(
                prot_xyz, prot_rdk, lig_xyz, lig_rdk, self.contact_bins, distances, self.radius)
    ]

class SplifVoxelizer(ComplexFeaturizer):
  """Computes SPLIF voxel grid for a macromolecular complex.

  SPLIF fingerprints are based on a technique introduced in the
  following paper. 

  Da, C., and D. Kireev. "Structural protein–ligand interaction
  fingerprints (SPLIF) for structure-based virtual screening:
  method and benchmark study." Journal of chemical information
  and modeling 54.9 (2014): 2555-2561.

  The SPLIF voxelizer localizes local SPLIF descriptors in
  space, by assigning features to the voxel in which they
  originated. This technique may be useful for downstream
  learning methods such as convolutional networks.
  """

  def __init__(self, 
               contact_bins=None,
               radius=2,
               size=8,
               box_width=16.0,
               voxel_width=1.0):
    """
    Parameters
    ----------
    contact_bins: list[tuple] 
      List of contact bins. If not specified is set to default
      `[(0, 2.0), (2.0, 3.0), (3.0, 4.5)]`.
    radius : int, optional (default 2)
        Fingerprint radius used for circular fingerprints.
    size: int, optional (default 8)
      Length of generated bit vector.
    box_width: float, optional (default 16.0)
      Size of a box in which voxel features are calculated. Box
      is centered on a ligand centroid.
    voxel_width: float, optional (default 1.0)
      Size of a 3D voxel in a grid.
    """
    if contact_bins is None:
      self.contact_bins = SPLIF_CONTACT_BINS
    else:
      self.contact_bins = contact_bins
    self.size = size
    self.radius = radius
    self.box_width = box_width
    self.voxel_width = voxel_width
    self.voxels_per_edge = int(self.box_width / self.voxel_width)

  def _featurize_complex(self, mol, protein):
    """
    Compute featurization for a single mol/protein complex

    TODO(rbharath): This is very not ergonomic. I'd much prefer
    returning an vector instead of a list of two vectors. In
    addition, there's a question of efficiency.
    RdkitGridFeaturizer caches rotated versions etc internally.
    To make things work out of box, we are accepting that
    kludgey input. This needs to be cleaned up before full
    merge.

    Parameters
    ----------
    mol: object
      Representation of the molecule
    protein: object
      Representation of the protein
    """
    (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein
    distances = compute_pairwise_distances(prot_xyz, lig_xyz)
    return [
        voxelize(
            convert_atom_pair_to_voxel,
            self.voxels_per_edge,
            self.box_width,
            self.voxel_width,
            hash_ecfp_pair, (prot_xyz, lig_xyz),
            feature_dict=splif_dict,
            nb_channel=self.size) for splif_dict in featurize_splif(
                prot_xyz, prot_rdk, lig_xyz, lig_rdk, self.contact_bins, distances, self.radius)
    ]
+0 −0

Empty file added.