Commit 2f144e31 authored by Bharath Ramsundar's avatar Bharath Ramsundar
Browse files

Trying to hbond featurizers working

parent f6f80625
Loading
Loading
Loading
Loading
+13 −1
Original line number Diff line number Diff line
@@ -54,7 +54,15 @@ def featurize_binding_pocket_ecfp(protein_xyz,
class ContactCircularFingerprint(ComplexFeaturizer):
  """Compute (Morgan) fingerprints near contact points of macromolecular complexes.

  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.
  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.

  For a macromolecular complex, returns a vector of shape
  `(2*size,)`
  """

  def __init__(self, 
@@ -120,6 +128,10 @@ class ContactCircularVoxelizer(ComplexFeaturizer):
  fragments which are sub-molecules centered at atoms in the
  contact region. Localize these ECFP fingeprints at the voxel
  in which they originated.

  Featurizes a macromolecular complex into a tensor of shape
  `(voxels_per_edge, voxels_per_edge, voxels_per_edge, size)`
  where `voxels_per_edge = int(box_width/voxel_width)`.
  """

  def __init__(self, 
+101 −22
Original line number Diff line number Diff line
@@ -4,6 +4,8 @@ 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.hash_utils import hash_ecfp_pair
from deepchem.utils.hash_utils import vectorize
from deepchem.utils.voxel_utils import voxelize 
from deepchem.utils.rdkit_util import compute_salt_bridges
from deepchem.utils.rdkit_util import compute_binding_pocket_cation_pi
@@ -11,9 +13,14 @@ 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
from deepchem.utils.rdkit_util import compute_pi_stack
from deepchem.utils.rdkit_util import compute_hydrogen_bonds

logger = logging.getLogger(__name__)


HBOND_DIST_BINS = [(2.2, 2.5), (2.5, 3.2), (3.2, 4.0)]
HBOND_ANGLE_CUTOFFS = [5, 50, 90]

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

@@ -37,9 +44,9 @@ class ChargeVoxelizer(ComplexFeaturizer):
  a local charge array. Sum contributions 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.
  Creates a tensor output of shape `(voxels_per_edge,
  voxels_per_edge, voxels_per_edge, 1)` for each macromolecular
  complex that computes the effective charge at each voxel.
  """
  def __init__(self, 
               box_width=16.0,
@@ -96,8 +103,8 @@ class SaltBridgeVoxelizer(ComplexFeaturizer):
  interact in a salt-bridge, the interaction is double counted
  in both voxels.

  Creates a tensor output of shape (box_width, box_width,
  box_width, 1) for each macromolecular complex that computes
  Creates a tensor output of shape `(voxels_per_edge,
  voxels_per_edge, voxels_per_edge, 1)` for each macromolecular
  the number of salt bridges at each voxel.
  """
  def __init__(self, 
@@ -158,8 +165,8 @@ class CationPiVoxelizer(ComplexFeaturizer):
  bridge in the voxel in which it originated to create a local
  cation-pi array.

  Creates a tensor output of shape (box_width, box_width,
  box_width, 1) for each macromolecular complex that computes
  Creates a tensor output of shape `(voxels_per_edge,
  voxels_per_edge, voxels_per_edge, 1)` for each macromolecular
  the number of cation-pi interactions at each voxel.
  """
  def __init__(self, 
@@ -232,11 +239,11 @@ class PiStackVoxelizer(ComplexFeaturizer):
  localize this salt bridge in the voxel in which it originated
  to create a local pi-stacking array.

  Creates a tensor output of shape (box_width, box_width,
  box_width, 2) for each macromolecular complex. Each voxel has
  2 fields, with the first tracking the number of pi-pi
  parallel interactions, and the second tracking the number of
  pi-T interactions.
  Creates a tensor output of shape `(voxels_per_edge,
  voxels_per_edge, voxels_per_edge, 2)` for each macromolecular
  Each voxel has 2 fields, with the first tracking the number
  of pi-pi parallel interactions, and the second tracking the
  number of pi-T interactions.
  """
  def __init__(self, 
               distance_cutoff=4.4,
@@ -324,6 +331,75 @@ class PiStackVoxelizer(ComplexFeaturizer):
        nb_channel=1)
    return [pi_parallel_tensor, pi_t_tensor]

class HydrogenBondCounter(ComplexFeaturizer):
  """Counts hydrogen bonds between atoms in macromolecular complexes.

  Given a macromolecular complex made up of multiple
  constitutent molecules, count the number hydrogen bonds
  between atoms in the macromolecular complex.

  Creates a scalar output of shape `(1,)` for each
  macromolecular that computes the total number of hydrogen
  bonds.
  """
  def __init__(self, 
               distance_bins=None,
               angle_cutoffs=None):
    """
    Parameters
    ----------
    distance_bins: list[tuple] 
      List of hydgrogen bond distance bins. If not specified is
      set to default
      `[(2.2, 2.5), (2.5, 3.2), (3.2, 4.0)]`.
    angle_cutoffs: list[float]
      List of hydrogen bond angle cutoffs. Max allowed
      deviation from the ideal (180 deg) angle between
      hydrogen-atom1, hydrogen-atom2 vectors.If not specified
      is set to default `[5, 50, 90]`
    """
    if distance_bins is None:
      self.distance_bins = HBOND_DIST_BINS
    else:
      self.distance_bins = distance_bins
    if angle_cutoffs is None:
      self.angle_cutoffs = HBOND_ANGLE_CUTOFFS
    else:
      self.angle_cutoffs = angle_cutoffs
    ######################################################
    print("constructor")
    print("self.distance_bins")
    print(self.distance_bins)
    print("self.angle_cutoffs")
    print(self.angle_cutoffs)
    ######################################################

  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)
    ######################################################
    print("self.distance_bins")
    print(self.distance_bins)
    print("self.angle_cutoffs")
    print(self.angle_cutoffs)
    ######################################################
    return [
        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.distance_bins, self.angle_cutoffs)
    ]

class HydrogenBondVoxelizer(ComplexFeaturizer):
  """Localize hydrogen bonds between atoms in macromolecular complexes.

@@ -335,8 +411,8 @@ class HydrogenBondVoxelizer(ComplexFeaturizer):
  different voxels interact in a hydrogen bond, the interaction
  is double counted in both voxels.

  Creates a tensor output of shape (box_width, box_width,
  box_width, 1) for each macromolecular complex that computes
  Creates a tensor output of shape `(voxels_per_edge,
  voxels_per_edge, voxels_per_edge, 1)` for each macromolecular
  the number of hydrogen bonds at each voxel.
  """
  def __init__(self, 
@@ -347,13 +423,15 @@ class HydrogenBondVoxelizer(ComplexFeaturizer):
    """
    Parameters
    ----------
    distance_cutoff: float, optional (default 4.0)
      The distance in angstroms within which atoms must be to
      be considered for a hydrogen bond interaction between
      them.
    angle_cutoff: float, optional (default 40.0)
      Angle cutoff. Max allowed deviation from the ideal (180
      deg) angle between hydrogen-atom1, hydrogen-atom2 vectors.
    distance_bins: list[tuple] 
      List of hydgrogen bond distance bins. If not specified is
      set to default
      `[(2.2, 2.5), (2.5, 3.2), (3.2, 4.0)]`.
    angle_cutoffs: list[float]
      List of hydrogen bond angle cutoffs. Max allowed
      deviation from the ideal (180 deg) angle between
      hydrogen-atom1, hydrogen-atom2 vectors.If not specified
      is set to default `[5, 50, 90]`
    box_width: float, optional (default 16.0)
      Size of a box in which voxel features are calculated. Box
      is centered on a ligand centroid.
@@ -361,6 +439,7 @@ class HydrogenBondVoxelizer(ComplexFeaturizer):
      Size of a 3D voxel in a grid.
    """
    self.distance_cutoff = distance_cutoff
    self.angle_cutoff = angle_cutoff
    self.box_width = box_width
    self.voxel_width = voxel_width
    self.voxels_per_edge = int(self.box_width / self.voxel_width)
+6 −7
Original line number Diff line number Diff line
@@ -21,6 +21,7 @@ from deepchem.feat.grid_featurizers import ChargeVoxelizer
from deepchem.feat.grid_featurizers import SaltBridgeVoxelizer
from deepchem.feat.grid_featurizers import CationPiVoxelizer
from deepchem.feat.grid_featurizers import PiStackVoxelizer
from deepchem.feat.grid_featurizers import HydrogenBondCounter
from deepchem.feat.grid_featurizers import HydrogenBondVoxelizer
from deepchem.utils.rdkit_util import compute_hydrogen_bonds
from deepchem.utils.rdkit_util import load_molecule
@@ -213,6 +214,9 @@ class RdkitGridFeaturizer(ComplexFeaturizer):
        angle_cutoff=self.cutoffs['pi_stack_angle_cutoff'],
        box_width=self.box_width,
        voxel_width=self.voxel_width)
    self.hbond_counter = HydrogenBondCounter(
        distance_bins=self.cutoffs['hbond_dist_cutoff'],
        angle_cutoffs=self.cutoffs['hbond_angle_cutoff'])
    self.hbond_voxelizer = HydrogenBondVoxelizer(
        distance_cutoff=self.cutoffs['hbond_dist_cutoff'],
        angle_cutoff=self.cutoffs['hbond_angle_cutoff'],
@@ -272,12 +276,8 @@ class RdkitGridFeaturizer(ComplexFeaturizer):
      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=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'])
      ]
      return self.hbond_counter._featurize_complex((lig_xyz, lig_rdk),
                                                   (prot_xyz, prot_rdk))
    if feature_name == 'ecfp':
      return self.contact_voxelizer._featurize_complex((lig_xyz, lig_rdk),
                                                       (prot_xyz, prot_rdk))
@@ -352,7 +352,6 @@ 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)
+11 −0
Original line number Diff line number Diff line
@@ -79,6 +79,13 @@ class SplifFingerprint(ComplexFeaturizer):
  each ligand atom, it identifies close protein atoms. These
  atom pairs are expanded to 2D circular fragments and a
  fingerprint for the union is turned on in the bit vector.

  This is conceptually pretty similar to
  `ContactCircularFingerprint` but computes ECFP fragments only
  for direct contacts instead of the entire contact region.

  For a macromolecular complex, returns a vector of shape
  `(2*size,)`
  """

  def __init__(self, 
@@ -145,6 +152,10 @@ class SplifVoxelizer(ComplexFeaturizer):
  space, by assigning features to the voxel in which they
  originated. This technique may be useful for downstream
  learning methods such as convolutional networks.

  Featurizes a macromolecular complex into a tensor of shape
  `(voxels_per_edge, voxels_per_edge, voxels_per_edge, size)`
  where `voxels_per_edge = int(box_width/voxel_width)`.
  """

  def __init__(self, 
+1 −1
Original line number Diff line number Diff line
@@ -12,7 +12,7 @@ import tensorflow as tf
tf.random.set_seed(123)
import deepchem as dc

# Load Tox21 dataset
# Load QM7B dataset
tasks, datasets, transformers = dc.molnet.load_qm7b_from_mat()
train_dataset, valid_dataset, test_dataset = datasets