Commit 225b0934 authored by Bharath Ramsundar's avatar Bharath Ramsundar
Browse files

Starting rdkit grid featurizer cleanup

parent 8bfa5990
Loading
Loading
Loading
Loading
+2 −0
Original line number Diff line number Diff line
@@ -5,6 +5,7 @@ from deepchem.feat.base_classes import Featurizer
from deepchem.feat.base_classes import MolecularFeaturizer
from deepchem.feat.base_classes import StructureFeaturizer
from deepchem.feat.base_classes import CompositionFeaturizer
from deepchem.feat.base_classes import ReactionFeaturizer
from deepchem.feat.base_classes import ComplexFeaturizer
from deepchem.feat.base_classes import UserDefinedFeaturizer
from deepchem.feat.graph_features import ConvMolFeaturizer
@@ -21,6 +22,7 @@ from deepchem.feat.raw_featurizer import RawFeaturizer
from deepchem.feat.raw_featurizer import RawReactionFeaturizer
from deepchem.feat.atomic_coordinates import AtomicCoordinates
from deepchem.feat.atomic_coordinates import NeighborListComplexAtomicCoordinates
from deepchem.feat.atomic_coordinates import AtomicConvFeaturizer
from deepchem.feat.adjacency_fingerprints import AdjacencyFingerprint
from deepchem.feat.smiles_featurizers import SmilesToSeq, SmilesToImage
from deepchem.feat.materials_featurizers import ElementPropertyFingerprint, SineCoulombMatrix, StructureGraphFeaturizer
+46 −16
Original line number Diff line number Diff line
@@ -3,7 +3,8 @@ Atomic coordinate featurizer.
"""
import logging
import numpy as np
from deepchem.feat import Featurizer
import logging
from deepchem.feat import MolecularFeaturizer
from deepchem.feat import ComplexFeaturizer
from deepchem.utils import rdkit_util, pad_array
from deepchem.utils.rdkit_util import MoleculeLoadException
@@ -11,7 +12,7 @@ from deepchem.utils.rdkit_util import MoleculeLoadException
logger = logging.getLogger(__name__)


class AtomicCoordinates(Featurizer):
class AtomicCoordinates(MolecularFeaturizer):
  """
  Nx3 matrix of Cartesian coordinates [Angstrom]
  """
@@ -26,7 +27,6 @@ class AtomicCoordinates(Featurizer):
    mol : RDKit Mol
          Molecule.
    """

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

@@ -47,7 +47,7 @@ class AtomicCoordinates(Featurizer):
    return coords


def compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors,
def _compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors,
                           periodic_box_size):
  """Computes a neighbor list from atom coordinates."""
  N = coords.shape[0]
@@ -92,7 +92,7 @@ def get_coords(mol):
  return coords


class NeighborListAtomicCoordinates(Featurizer):
class NeighborListAtomicCoordinates(MolecularFeaturizer):
  """
  Adjacency List of neighbors in 3-space

@@ -138,7 +138,7 @@ class NeighborListAtomicCoordinates(Featurizer):
    # TODO(rbharath): Should this return a list?
    bohr_coords = self.coordinates_featurizer._featurize(mol)[0]
    coords = get_coords(mol)
    neighbor_list = compute_neighbor_list(coords, self.neighbor_cutoff,
    neighbor_list = _compute_neighbor_list(coords, self.neighbor_cutoff,
                                           self.max_num_neighbors,
                                           self.periodic_box_size)
    return (bohr_coords, neighbor_list)
@@ -178,13 +178,13 @@ class NeighborListComplexAtomicCoordinates(ComplexFeaturizer):
    protein_coords, protein_mol = rdkit_util.load_molecule(protein_pdb_file)
    system_coords = rdkit_util.merge_molecules_xyz([mol_coords, protein_coords])

    system_neighbor_list = compute_neighbor_list(
    system_neighbor_list = _compute_neighbor_list(
        system_coords, self.neighbor_cutoff, self.max_num_neighbors, None)

    return (system_coords, system_neighbor_list)


class ComplexNeighborListFragmentAtomicCoordinates(ComplexFeaturizer):
class AtomicConvFeaturizer(ComplexFeaturizer):
  """This class computes the featurization that corresponds to AtomicConvModel.

  This class computes featurizations needed for AtomicConvModel.
@@ -206,22 +206,35 @@ class ComplexNeighborListFragmentAtomicCoordinates(ComplexFeaturizer):
  """

  def __init__(self,
               frag1_num_atoms,
               frag2_num_atoms,
               frag_num_atoms,
               complex_num_atoms,
               max_num_neighbors,
               neighbor_cutoff,
               strip_hydrogens=True):
    self.frag1_num_atoms = frag1_num_atoms
    self.frag2_num_atoms = frag2_num_atoms
    self.complex_num_atoms = complex_num_atoms
    """Initialize an AtomicConvFeaturizer object.

    Parameters
    ----------
    frag_num_atoms: list[int]
      List of the number of atoms in each fragment.
    max_num_neighbors: int
      The maximum number of neighbors allowed
    neighbor_cutoff: float
      The distance in angstroms after which neighbors are cutoff.
    strip_hydrogens: bool, optional
      If true, remove hydrogens before featurizing.
    """
    # TODO(rbharath): extend to more fragments
    if len(frag_num_atoms) != 2:
      raise ValueError("Currently only supports two fragments")
    self.complex_num_atoms = sum(frag_num_atoms)
    self.max_num_neighbors = max_num_neighbors
    self.neighbor_cutoff = neighbor_cutoff
    self.strip_hydrogens = strip_hydrogens
    self.neighborlist_featurizer = NeighborListComplexAtomicCoordinates(
        self.max_num_neighbors, self.neighbor_cutoff)

  def _featurize(self, mol_pdb_file, protein_pdb_file):
    frag_coords = []
    frag_mols = []
    try:
      frag1_coords, frag1_mol = rdkit_util.load_molecule(
          mol_pdb_file, is_protein=False, sanitize=True, add_hydrogens=False)
@@ -265,7 +278,7 @@ class ComplexNeighborListFragmentAtomicCoordinates(ComplexFeaturizer):

  def featurize_mol(self, coords, mol, max_num_atoms):
    logging.info("Featurizing molecule of size: %d", len(mol.GetAtoms()))
    neighbor_list = compute_neighbor_list(coords, self.neighbor_cutoff,
    neighbor_list = _compute_neighbor_list(coords, self.neighbor_cutoff,
                                           self.max_num_neighbors, None)
    z = self.get_Z_matrix(mol, max_num_atoms)
    z = pad_array(z, max_num_atoms)
@@ -304,3 +317,20 @@ class ComplexNeighborListFragmentAtomicCoordinates(ComplexFeaturizer):
    mol = MoleculeShim(atomic_numbers)
    coords = coords[indexes_to_keep]
    return coords, mol


############################# Deprecation warning for old name of AtomicConvFeaturizer ###############################

DEPRECATION = "{} is deprecated and has been renamed to {} and will be removed in DeepChem 3.0."


class ComplexNeighborListFragmentAtomicCoordinates(AtomicConvFeaturizer):

  def __init__(self, *args, **kwargs):

    warnings.warn(
        DEPRECATION.format("ComplexNeighborListFragmentAtomicCoordinates",
                           "AtomicConvFeaturizer"), FutureWarning)

    super(ComplexNeighborListFragmentAtomicCoordinates, self).__init__(
        *args, **kwargs)
+0 −133
Original line number Diff line number Diff line
@@ -732,136 +732,3 @@ class WeaveFeaturizer(MolecularFeaturizer):
        graph_distance=self.graph_distance)

    return WeaveMol(nodes, pairs)


class AtomicConvFeaturizer(ComplexNeighborListFragmentAtomicCoordinates):
  """This class computes the Atomic Convolution features"""

  # TODO (VIGS25): Complete the description

  name = ['atomic_conv']

  def __init__(self,
               labels,
               neighbor_cutoff,
               frag1_num_atoms=70,
               frag2_num_atoms=634,
               complex_num_atoms=701,
               max_num_neighbors=12,
               batch_size=24,
               atom_types=[
                   6, 7., 8., 9., 11., 12., 15., 16., 17., 20., 25., 30., 35.,
                   53., -1.
               ],
               radial=[[
                   1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0,
                   7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0
               ], [0.0, 4.0, 8.0], [0.4]],
               layer_sizes=[32, 32, 16],
               strip_hydrogens=True,
               learning_rate=0.001,
               epochs=10):
    """
    Parameters

    labels: numpy.ndarray
      Labels which we want to predict using the model
    neighbor_cutoff: int
      TODO (VIGS25): Add description
    frag1_num_atoms: int
      Number of atoms in first fragment
    frag2_num_atoms: int
      Number of atoms in second fragment
    complex_num_atoms: int
      TODO (VIGS25) : Add description
    max_num_neighbors: int
      Maximum number of neighbors possible for an atom
    batch_size: int
      Batch size used for training and evaluation
    atom_types: list
      List of atoms recognized by model. Atoms are indicated by their
      nuclear numbers.
    radial: list
      TODO (VIGS25): Add description
    layer_sizes: list
      List of layer sizes for the AtomicConvolutional Network
    strip_hydrogens: bool
      Whether to remove hydrogens while computing neighbor features
    learning_rate: float
      Learning rate for training the model
    epochs: int
      Number of epochs to train the model for
    """

    self.atomic_conv_model = dc.models.atomic_conv.AtomicConvModel(
        frag1_num_atoms=frag1_num_atoms,
        frag2_num_atoms=frag2_num_atoms,
        complex_num_atoms=complex_num_atoms,
        max_num_neighbors=max_num_neighbors,
        batch_size=batch_size,
        atom_types=atom_types,
        radial=radial,
        layer_sizes=layer_sizes,
        learning_rate=learning_rate)

    super(AtomicConvFeaturizer, self).__init__(
        frag1_num_atoms=frag1_num_atoms,
        frag2_num_atoms=frag2_num_atoms,
        complex_num_atoms=complex_num_atoms,
        max_num_neighbors=max_num_neighbors,
        neighbor_cutoff=neighbor_cutoff,
        strip_hydrogens=strip_hydrogens)

    self.epochs = epochs
    self.labels = labels

  def featurize(self, mol_files, protein_files):
    features = []
    failures = []
    for i, (mol_file, protein_pdb) in enumerate(zip(mol_files, protein_files)):
      logging.info("Featurizing %d / %d" % (i, len(mol_files)))
      new_features = self._featurize(mol_file, protein_pdb)
      # Handle loading failures which return None
      if new_features is not None:
        features.append(new_features)
      else:
        failures.append(ind)

    features = np.asarray(features)
    labels = np.delete(self.labels, failures)
    dataset = DiskDataset.from_numpy(features, labels)

    # Fit atomic conv model
    self.atomic_conv_model.fit(dataset, nb_epoch=self.epochs)

    # Add the Atomic Convolution layers to fetches
    layers_to_fetch = [
        self.atomic_conv_model._frag1_conv, self.atomic_conv_model._frag2_conv,
        self.atomic_conv_model._complex_conv
    ]

    # Extract the atomic convolution features
    atomic_conv_features = list()
    batch_generator = self.atomic_conv_model.default_generator(
        dataset=dataset, epochs=1)

    for X, y, w in batch_generator:
      frag1_conv, frag2_conv, complex_conv = self.atomic_conv_model.predict_on_generator(
          [(X, y, w)], outputs=layers_to_fetch)
      concatenated = np.concatenate(
          [frag1_conv, frag2_conv, complex_conv], axis=1)
      atomic_conv_features.append(concatenated)

    batch_size = self.atomic_conv_model.batch_size

    if len(features) % batch_size != 0:
      num_batches = (len(features) // batch_size) + 1
      num_to_skip = num_batches * batch_size - len(features)
    else:
      num_to_skip = 0

    atomic_conv_features = np.asarray(atomic_conv_features)
    atomic_conv_features = atomic_conv_features[-num_to_skip:]
    atomic_conv_features = np.squeeze(atomic_conv_features)

    return atomic_conv_features, failures
+55 −269

File changed.

Preview size limit exceeded, changes collapsed.

+36 −38
Original line number Diff line number Diff line
@@ -10,7 +10,7 @@ from deepchem.feat.atomic_coordinates import get_coords
from deepchem.feat.atomic_coordinates import AtomicCoordinates
from deepchem.feat.atomic_coordinates import NeighborListAtomicCoordinates
from deepchem.feat.atomic_coordinates import NeighborListComplexAtomicCoordinates
from deepchem.feat.atomic_coordinates import ComplexNeighborListFragmentAtomicCoordinates
from deepchem.feat.atomic_coordinates import AtomicConvFeaturizer 

logger = logging.getLogger(__name__)

@@ -158,41 +158,39 @@ class TestAtomicCoordinates(unittest.TestCase):
    for atom in range(N):
      assert len(system_neighbor_list[atom]) <= max_num_neighbors

  def test_full_complex_featurization(self):
    """Unit test for AtomicConvFeaturizer."""
    dir_path = os.path.dirname(os.path.realpath(__file__))
    ligand_file = os.path.join(dir_path, "data/3zso_ligand_hyd.pdb")
    protein_file = os.path.join(dir_path, "data/3zso_protein.pdb")
    # Pulled from PDB files. For larger datasets with more PDBs, would use
    # max num atoms instead of exact.
    frag1_num_atoms = 44  # for ligand atoms
    frag2_num_atoms = 2336  # for protein atoms
    complex_num_atoms = 2380  # in total
    max_num_neighbors = 4
    # Cutoff in angstroms
    neighbor_cutoff = 4
    complex_featurizer = AtomicConvFeaturizer(
        frag1_num_atoms, frag2_num_atoms, complex_num_atoms, max_num_neighbors,
        neighbor_cutoff)
    (frag1_coords, frag1_neighbor_list, frag1_z, frag2_coords,
     frag2_neighbor_list, frag2_z, complex_coords,
     complex_neighbor_list, complex_z) = complex_featurizer._featurize_complex(
         ligand_file, protein_file)

# TODO(rbharath): This test will be uncommented in the next PR up on the docket.
#  def test_full_complex_featurization(self):
#    """Unit test for ComplexNeighborListFragmentAtomicCoordinates."""
#    dir_path = os.path.dirname(os.path.realpath(__file__))
#    ligand_file = os.path.join(dir_path, "data/3zso_ligand_hyd.pdb")
#    protein_file = os.path.join(dir_path, "data/3zso_protein.pdb")
#    # Pulled from PDB files. For larger datasets with more PDBs, would use
#    # max num atoms instead of exact.
#    frag1_num_atoms = 44  # for ligand atoms
#    frag2_num_atoms = 2336  # for protein atoms
#    complex_num_atoms = 2380  # in total
#    max_num_neighbors = 4
#    # Cutoff in angstroms
#    neighbor_cutoff = 4
#    complex_featurizer = ComplexNeighborListFragmentAtomicCoordinates(
#        frag1_num_atoms, frag2_num_atoms, complex_num_atoms, max_num_neighbors,
#        neighbor_cutoff)
#    (frag1_coords, frag1_neighbor_list, frag1_z, frag2_coords,
#     frag2_neighbor_list, frag2_z, complex_coords,
#     complex_neighbor_list, complex_z) = complex_featurizer._featurize_complex(
#         ligand_file, protein_file)
#
#    assert frag1_coords.shape == (frag1_num_atoms, 3)
#    self.assertEqual(
#        sorted(list(frag1_neighbor_list.keys())), list(range(frag1_num_atoms)))
#    self.assertEqual(frag1_z.shape, (frag1_num_atoms,))
#
#    self.assertEqual(frag2_coords.shape, (frag2_num_atoms, 3))
#    self.assertEqual(
#        sorted(list(frag2_neighbor_list.keys())), list(range(frag2_num_atoms)))
#    self.assertEqual(frag2_z.shape, (frag2_num_atoms,))
#
#    self.assertEqual(complex_coords.shape, (complex_num_atoms, 3))
#    self.assertEqual(
#        sorted(list(complex_neighbor_list.keys())),
#        list(range(complex_num_atoms)))
#    self.assertEqual(complex_z.shape, (complex_num_atoms,))
    self.assertEqual(frag1_coords.shape, (frag1_num_atoms, 3))
    self.assertEqual(
        sorted(list(frag1_neighbor_list.keys())), list(range(frag1_num_atoms)))
    self.assertEqual(frag1_z.shape, (frag1_num_atoms,))

    self.assertEqual(frag2_coords.shape, (frag2_num_atoms, 3))
    self.assertEqual(
        sorted(list(frag2_neighbor_list.keys())), list(range(frag2_num_atoms)))
    self.assertEqual(frag2_z.shape, (frag2_num_atoms,))

    self.assertEqual(complex_coords.shape, (complex_num_atoms, 3))
    self.assertEqual(
        sorted(list(complex_neighbor_list.keys())),
        list(range(complex_num_atoms)))
    self.assertEqual(complex_z.shape, (complex_num_atoms,))
Loading