Commit 5ad01702 authored by Bharath Ramsundar's avatar Bharath Ramsundar
Browse files

Start on voxel testing

parent c89bc79f
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -27,3 +27,4 @@ from deepchem.feat.adjacency_fingerprints import AdjacencyFingerprint
from deepchem.feat.smiles_featurizers import SmilesToSeq, SmilesToImage
from deepchem.feat.materials_featurizers import ElementPropertyFingerprint, SineCoulombMatrix, StructureGraphFeaturizer
from deepchem.feat.contact_fingerprints import ContactCircularFingerprint
from deepchem.feat.contact_fingerprints import ContactCircularVoxelizer
+9 −0
Original line number Diff line number Diff line
@@ -237,6 +237,11 @@ class AtomicConvFeaturizer(ComplexFeaturizer):
  z-matrix for the molecule which is an array of shape
  (N_atoms,) that contains the atomic number of that atom.

  If `reduce_to_contacts` is set to `True`, only those atoms which are
  in a contact region of the molecular complex are considered. This
  can dramatically reduce the number of atoms that need to be
  featurized so this feature is set to True by default.

  Since the featurization computes these three quantities for
  each of the two molecules and the complex, a total of 9
  quantities are returned for each complex. Note that for
@@ -310,6 +315,10 @@ class AtomicConvFeaturizer(ComplexFeaturizer):
    if self.reduce_to_contacts:
      fragments = reduce_molecular_complex_to_contacts(fragments, self.cutoff)
    coords = [frag[0] for frag in fragments]
    ################################################
    print("[coord.shape for coord in coords]")
    print([coord.shape for coord in coords])
    ################################################
    mols = [frag[1] for frag in fragments]
    #system_mol = rdkit_util.merge_molecules(mols)
    system_mol = rdkit_util.merge_molecular_fragments(mols)
+48 −37
Original line number Diff line number Diff line
@@ -90,14 +90,6 @@ class ContactCircularFingerprint(ComplexFeaturizer):
    """
    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
    ----------
    molecular_complex: Object
@@ -171,20 +163,35 @@ class ContactCircularVoxelizer(ComplexFeaturizer):
    self.voxel_width = voxel_width
    self.voxels_per_edge = int(self.box_width / self.voxel_width)

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

    Parameters
    ----------
    mol: object
      Representation of the molecule
    protein: object
      Representation of the protein
    molecular_complex: Object
      A representation of a molecular complex, produced by
      `rdkit_util.load_complex`.
    """
    (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein
    distances = compute_pairwise_distances(prot_xyz, lig_xyz)
    return [

    # TODO(rbharath): This is a little tricky in the generalized
    # regime, but we need to find a way to compute the centroid. My
    # idea is that we can compute the centroid of the contact
    # atoms.and use this to recenter all the fragments.
    try:
      fragments = rdkit_util.load_complex(molecular_complex, add_hydrogens=False)

    except MoleculeLoadException:
      logger.warning("This molecule cannot be loaded by Rdkit. Returning None")
      return None
    pairwise_features = []
    # We compute pairwise contact fingerprints
    for (frag1, frag2) in itertools.combinations(fragments, 2):
      distances = compute_pairwise_distances(frag1[0], frag2[0])
      xyzs = [frag1[0], frag2[0]]
      #(lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein
      #distances = compute_pairwise_distances(prot_xyz, lig_xyz)
      pairwise_features.append(
          sum([
              voxelize(
                  convert_atom_to_voxel,
@@ -195,14 +202,18 @@ class ContactCircularVoxelizer(ComplexFeaturizer):
                  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,
              for xyz, ecfp_dict in zip(xyzs,
                                        featurize_contacts_ecfp(
                                            frag1,
                                            frag2,
                                            distances,
                                            cutoff=self.cutoff,
                                            ecfp_degree=self.radius))
          ])
    ]
      )
    ############################################
    print("[feat.shape for feat in pairwise_features]")
    print([feat.shape for feat in pairwise_features])
    ############################################
    # Features are of shape (voxels_per_edge, voxels_per_edge, voxels_per_edge, num_feat) so we should concatenate on the last axis.
    return np.concatenate(pairwise_features, axis=-1)
+27 −2
Original line number Diff line number Diff line
@@ -648,11 +648,36 @@ def get_contact_atom_indices(fragments, cutoff=4.5):
  #return atoms_to_keep


def compute_contact_centroid(molecular_complex, cutoff=4.5):
  """Computes the (x,y,z) centroid of the contact regions of this molecular complex.

  For a molecular complex, it's necessary for various featurizations
  that compute voxel grids to find a reasonable center for the
  voxelization. This function computes the centroid of all the contact
  atoms, defined as an atom that's within `cutoff` Angstroms of an
  atom from a different molecule.

  Parameters
  ----------
  molecular_complex: Object
    A representation of a molecular complex, produced by
    `rdkit_util.load_complex`.
  cutoff: float, optional
    The distance in Angstroms considered for computing contacts.
  """
  fragments = reduce_molecular_complex_to_contacts(molecular_complex, cutoff)
  coords = [frag[0] for frag in fragments]
  contact_coords = merge_molecules_xyz(coords) 
  centroid = np.mean(contact_coords, axis=0)
  return (centroid)

def compute_centroid(coordinates):
  """Compute the x,y,z centroid of provided coordinates
  """Compute the (x,y,z) centroid of provided coordinates

  Parameters
  ----------
  coordinates: np.ndarray
    Shape (N, 3), where N is number atoms.
    Shape `(N, 3)`, where `N` is the number of atoms.
  """
  centroid = np.mean(coordinates, axis=0)
  return (centroid)
+9 −2
Original line number Diff line number Diff line
@@ -4,14 +4,21 @@ import logging
from deepchem.molnet.load_function.pdbbind_datasets import get_pdbbind_molecular_complex_files

complex_files = get_pdbbind_molecular_complex_files(subset="core", version="v2015", interactions="protein-ligand", load_binding_pocket=False)
core_subset = complex_files[:2]
n_complexes = len(complex_files)
n_featurize = n_complexes
core_subset = complex_files[:n_featurize]

max_num_neighbors = 4
# Cutoff in angstroms
neighbor_cutoff = 4
featurizer = dc.feat.AtomicConvFeaturizer(
    frag_max_atoms=70,
    frag_max_atoms=[100, 50],
    max_num_neighbors=max_num_neighbors,
    neighbor_cutoff=neighbor_cutoff)
features, failures = featurizer.featurize_complexes(
    core_subset, parallelize=False)

print("features.shape")
print(features.shape)
print("%d failures" % len(failures))
assert features.shape == (n_featurize -len(failures), 2)
Loading