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

Cleaning up rdkit grid featurizer

parent 1a2050ba

deepchem/feat/grid_featurizers.py

+72 −4
Original line number Diff line number Diff line
@@ -30,7 +30,7 @@ def compute_charge_dictionary(molecule): 
class ChargeVoxelizer(ComplexFeaturizer):

  """Localize partial charges of atoms in macromolecular complexes.



  Given a macromolecular compelx made up of multiple

  Given a macromolecular complex 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
@@ -88,7 +88,7 @@ class ChargeVoxelizer(ComplexFeaturizer): 
class SaltBridgeVoxelizer(ComplexFeaturizer):

  """Localize salt bridges between atoms in macromolecular complexes.



  Given a macromolecular compelx made up of multiple

  Given a macromolecular complex made up of multiple

  constitutent molecules, compute salt bridges between atoms in

  the macromolecular complex. For each atom, localize this salt

  bridge in the voxel in which it originated to create a local
@@ -152,7 +152,7 @@ class SaltBridgeVoxelizer(ComplexFeaturizer): 
class CationPiVoxelizer(ComplexFeaturizer):

  """Localize cation-Pi interactions between atoms in macromolecular complexes.



  Given a macromolecular compelx made up of multiple

  Given a macromolecular complex made up of multiple

  constitutent molecules, compute cation-Pi between atoms in

  the macromolecular complex. For each atom, localize this salt

  bridge in the voxel in which it originated to create a local
@@ -226,7 +226,7 @@ class CationPiVoxelizer(ComplexFeaturizer): 
class PiStackVoxelizer(ComplexFeaturizer):

  """Localize Pi stacking interactions between atoms in macromolecular complexes.



  Given a macromolecular compelx made up of multiple

  Given a macromolecular complex made up of multiple

  constitutent molecules, compute pi-stacking interactions

  between atoms in the macromolecular complex. For each atom,

  localize this salt bridge in the voxel in which it originated
@@ -323,3 +323,71 @@ class PiStackVoxelizer(ComplexFeaturizer): 
        feature_dict=ligand_pi_t,

        nb_channel=1)

    return [pi_parallel_tensor, pi_t_tensor]

    

class HydrogenBondVoxelizer(ComplexFeaturizer):

  """Localize hydrogen bonds between atoms in macromolecular complexes.



  Given a macromolecular complex made up of multiple

  constitutent molecules, compute hydrogen bonds between atoms

  in the macromolecular complex. For each atom, localize this

  hydrogen bond in the voxel in which it originated to create a

  local hydrogen bond array. Note that if atoms in two

  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

  the number of hydrogen bonds at each voxel.

  """

  def __init__(self, 

               distance_cutoff=5.0,

               angle_cutoff=40.0,

               box_width=16.0,

               voxel_width=1.0):

    """

    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.

    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.distance_cutoff = distance_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

    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,

            None, (prot_xyz, lig_xyz),

            feature_list=hbond_list,

            nb_channel=1) for hbond_list in compute_hydrogen_bonds(

                prot_xyz, prot_rdk, lig_xyz, lig_rdk,

                distances, self.distance_cutoff,

                self.angle_cutoff)

    ]

deepchem/feat/rdkit_grid_featurizer.py

+11 −40
Original line number Diff line number Diff line
@@ -21,6 +21,8 @@ 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 HydrogenBondVoxelizer

from deepchem.utils.rdkit_util import compute_hydrogen_bonds

from deepchem.utils.rdkit_util import load_molecule

from deepchem.utils.rdkit_util import compute_centroid

from deepchem.utils.rdkit_util import subtract_centroid
@@ -56,45 +58,6 @@ VOXEL_FEATURES = [ 
]





def compute_hbonds_in_range(protein, protein_xyz, ligand, ligand_xyz,

                            pairwise_distances, hbond_dist_bin,

                            hbond_angle_cutoff):

  """

  Find all pairs of (protein_index_i, ligand_index_j) that hydrogen bond given

  a distance bin and an angle cutoff.

  """



  contacts = np.nonzero((pairwise_distances > hbond_dist_bin[0]) &

                        (pairwise_distances < hbond_dist_bin[1]))

  contacts = zip(contacts[0], contacts[1])

  hydrogen_bond_contacts = []

  for contact in contacts:

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

                        hbond_angle_cutoff):

      hydrogen_bond_contacts.append(contact)

  return hydrogen_bond_contacts





def compute_hydrogen_bonds(protein_xyz, protein, ligand_xyz, ligand,

                           pairwise_distances, hbond_dist_bins,

                           hbond_angle_cutoffs):

  """Computes hydrogen bonds between proteins and ligands.



  Returns a list of sublists. Each sublist is a series of tuples of

  (protein_index_i, ligand_index_j) that represent a hydrogen bond. Each sublist

  represents a different type of hydrogen bond.

  """



  hbond_contacts = []

  for i, hbond_dist_bin in enumerate(hbond_dist_bins):

    hbond_angle_cutoff = hbond_angle_cutoffs[i]

    hbond_contacts.append(

        compute_hbonds_in_range(protein, protein_xyz, ligand, ligand_xyz,

                                pairwise_distances, hbond_dist_bin,

                                hbond_angle_cutoff))

  return (hbond_contacts)





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).
@@ -160,6 +123,8 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 
    pi_stack_angle_cutoff: 30.0

    cation_pi_dist_cutoff: 6.5

    cation_pi_angle_cutoff: 30.0

    hbond_dist_cutoff: 4.0

    hbond_angle_cutoff: 40.0

    """



    # list of features that require sanitized molecules
@@ -183,6 +148,8 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 
        'pi_stack_angle_cutoff': 30.0,

        'cation_pi_dist_cutoff': 6.5,

        'cation_pi_angle_cutoff': 30.0,

        'hbond_dist_cutoff': 4.0,

        'hbond_angle_cutoff': 40.0

    }



    # update with cutoffs specified by the user
@@ -206,7 +173,6 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 
    ignored_features = []

    if self.sanitize is False:

      ignored_features += require_sanitized

    #ignored_features += not_implemented



    # Intantiate Featurizers

    self.ecfp_featurizer = CircularFingerprint(
@@ -247,6 +213,11 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 
        angle_cutoff=self.cutoffs['pi_stack_angle_cutoff'],

        box_width=self.box_width,

        voxel_width=self.voxel_width)

    self.hbond_voxelizer = HydrogenBondVoxelizer(

        distance_cutoff=self.cutoffs['hbond_dist_cutoff'],

        angle_cutoff=self.cutoffs['hbond_angle_cutoff'],

        box_width=self.box_width,

        voxel_width=self.voxel_width)



    # parse provided feature types

    for feature_type in feature_types:

deepchem/feat/tests/test_rdkit_grid_features.py

+2 −0
Original line number Diff line number Diff line
@@ -337,6 +337,8 @@ class TestRdkitGridFeaturizer(unittest.TestCase): 
        'pi_stack_angle_cutoff': 15.0,

        'cation_pi_dist_cutoff': 5.5,

        'cation_pi_angle_cutoff': 20.0,

        'hbond_dist_cutoff': 4.0,

        'hbond_angle_cutoff': 40.0

    }

    rgf_featurizer = rgf.RdkitGridFeaturizer(**custom_cutoffs)

    self.assertEqual(rgf_featurizer.cutoffs, custom_cutoffs)

deepchem/splits/splitters.py

+10 −9
Original line number Diff line number Diff line
@@ -375,18 +375,19 @@ class RandomGroupSplitter(Splitter): 
class RandomStratifiedSplitter(Splitter):

  """RandomStratified Splitter class.



  For sparse multitask datasets, a standard split offers no guarantees

  that the splits will have any activate compounds. This class guarantees

  that each task will have a proportional split of the activates in a

  split. TO do this, a ragged split is performed with different numbers

  of compounds taken from each task. Thus, the length of the split arrays

  may exceed the split of the original array. That said, no datapoint is

  copied to more than one split, so correctness is still ensured.

  For sparse multitask datasets, a standard split offers no

  guarantees that the splits will have any activate compounds.

  This class guarantees that each task will have a proportional

  split of the activates in a split. TO do this, a ragged split

  is performed with different numbers of compounds taken from

  each task. Thus, the length of the split arrays may exceed the

  split of the original array. That said, no datapoint is copied

  to more than one split, so correctness is still ensured.



  Note that this splitter is only valid for boolean label data.



  TODO(rbharath): This splitter should be refactored to match style of

  other splitter classes.

  TODO(rbharath): This splitter should be refactored to match

  style of other splitter classes.

  """



  def __generate_required_hits(self, w, frac_split):

deepchem/utils/rdkit_util.py

+45 −4
Original line number Diff line number Diff line
@@ -20,6 +20,8 @@ from deepchem.utils.pdbqt_utils import convert_protein_to_pdbqt 
from deepchem.utils.geometry_utils import angle_between

from deepchem.utils.geometry_utils import is_angle_within_cutoff

from deepchem.utils.geometry_utils import generate_random_rotation_matrix

from rdkit import Chem

from rdkit.Chem.rdchem import AtomValenceException



logger = logging.getLogger(__name__)
@@ -403,6 +405,7 @@ def merge_molecules(molecules): 
      combined = rdmolops.CombineMols(combined, nextmol)

    return combined





def is_hydrogen_bond(protein_xyz,

                     protein,

                     ligand_xyz,
@@ -443,12 +446,50 @@ def is_hydrogen_bond(protein_xyz, 


    for hydrogen_xyz in hydrogens:

      hydrogen_to_ligand = ligand_atom_xyz - hydrogen_xyz

      hydrogen_to_protein = protein_atom_xyz - hydrogen_xyzh

      if math.abs(180 - angle_between(hydrogen_to_protein, hydrogen_to_ligand) * 180.0 / np.pi) <= hbond_angle_cutoff:

      hydrogen_to_protein = protein_atom_xyz - hydrogen_xyz

      if np.abs(180 - angle_between(hydrogen_to_protein, hydrogen_to_ligand) * 180.0 / np.pi) <= hbond_angle_cutoff:

        return True

  return False





def compute_hbonds_in_range(protein, protein_xyz, ligand, ligand_xyz,

                            pairwise_distances, hbond_dist_bin,

                            hbond_angle_cutoff):

  """

  Find all pairs of (protein_index_i, ligand_index_j) that hydrogen bond given



  a distance bin and an angle cutoff.

  """



  contacts = np.nonzero((pairwise_distances > hbond_dist_bin[0]) &

                        (pairwise_distances < hbond_dist_bin[1]))

  contacts = zip(contacts[0], contacts[1])

  hydrogen_bond_contacts = []

  for contact in contacts:

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

                        hbond_angle_cutoff):

      hydrogen_bond_contacts.append(contact)

  return hydrogen_bond_contacts



def compute_hydrogen_bonds(protein_xyz, protein, ligand_xyz, ligand,

                           pairwise_distances, hbond_dist_bins,

                           hbond_angle_cutoffs):

  """Computes hydrogen bonds between proteins and ligands.



  Returns a list of sublists. Each sublist is a series of tuples

  of (protein_index_i, ligand_index_j) that represent a hydrogen

  bond. Each sublist represents a different type of hydrogen

  bond.

  """



  hbond_contacts = []

  for i, hbond_dist_bin in enumerate(hbond_dist_bins):

    hbond_angle_cutoff = hbond_angle_cutoffs[i]

    hbond_contacts.append(

        compute_hbonds_in_range(protein, protein_xyz, ligand, ligand_xyz,

                                pairwise_distances, hbond_dist_bin,

                                hbond_angle_cutoff))

  return (hbond_contacts)



def compute_salt_bridges(first,

                         second,

CRA Git | Maintained and supported by SUSTech CRA and CCSE