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

Preliminary PDBbind implementation

parent 929fd14b

deepchem/feat/__init__.py

+1 −0
Original line number Diff line number Diff line
@@ -16,3 +16,4 @@ from deepchem.feat.graph_features import ConvMolFeaturizer 
from deepchem.feat.fingerprints import CircularFingerprint

from deepchem.feat.basic import RDKitDescriptors

from deepchem.feat.coulomb_matrices import CoulombMatrixEig

from deepchem.feat.grid_featurizer import GridFeaturizer

deepchem/feat/grid_featurizer.py

+55 −24
Original line number Diff line number Diff line
@@ -6,20 +6,19 @@ __author__ = "Bharath Ramsundar and Evan Feinberg" 
__copyright__ = "Copyright 2016, Stanford University"

__license__ = "GPL"



from copy import deepcopy

import numpy as np

import os

import shutil

import time

from collections import deque

import tempfile

import hashlib

import sys

import numpy as np

from copy import deepcopy

import openbabel as ob

from collections import deque

from functools import partial

from deepchem.feat import ComplexFeaturizer

from deepchem.utils.save import log

import tempfile

import os

import shutil

import time





"""
@@ -30,7 +29,11 @@ def get_xyz_from_ob(ob_mol): 
  returns an m x 3 np array of 3d coords

  of given openbabel molecule

  """



  ######################################################### DEBUG

  #print("get_xyz_from_ob")

  #print("ob_mol.NumAtoms()")

  #print(ob_mol.NumAtoms())

  ######################################################### DEBUG

  xyz = np.zeros((ob_mol.NumAtoms(), 3))

  for i, atom in enumerate(ob.OBMolAtomIter(ob_mol)):

    xyz[i, 0] = atom.x()
@@ -38,6 +41,16 @@ def get_xyz_from_ob(ob_mol): 
    xyz[i, 2] = atom.z()

  return(xyz)



def get_ligand_filetype(ligand_filename):

  """Returns the filetype of ligand."""

  if ".mol2" in ligand_filename:

    return ".mol2"

  elif ".sdf" in ligand_filename:

    return "sdf"

  elif ".pdb" in ligand_filename:

    return ".pdb"

  else:

    raise ValueError("Unrecognized_filename")



def load_molecule(molecule_file, remove_hydrogens=True,

                  calc_charges=False):
@@ -875,15 +888,19 @@ class GridFeaturizer(ComplexFeaturizer): 
                        "S3", "S3+", "S2", "So2", "Sox" "Sac" "SO", "P3", 

                        "P", "P3+", "F", "Cl", "Br", "I"]



  def _featurize_complex(self, ligand_pdb_lines, protein_pdb_lines):

  def _featurize_complex(self, ligand_ext, ligand_lines, protein_pdb_lines):

    tempdir = tempfile.mkdtemp()



    ############################################################### DEBUG

    #print("ligand_lines")

    #print(ligand_lines)

    ############################################################### DEBUG

    ############################################################## TIMING

    time1 = time.time()

    ############################################################## TIMING

    ligand_pdb_file = os.path.join(tempdir, "ligand.pdb")

    with open(ligand_pdb_file, "w") as mol_f:

      mol_f.writelines(ligand_pdb_lines)

    ligand_file = os.path.join(tempdir, "ligand.%s" % ligand_ext)

    with open(ligand_file, "w") as mol_f:

      mol_f.writelines(ligand_lines)

    ############################################################## TIMING

    time2 = time.time()

    log("TIMING: Writing ligand took %0.3f s" % (time2-time1), self.verbose)
@@ -900,32 +917,36 @@ class GridFeaturizer(ComplexFeaturizer): 
    log("TIMING: Writing protein took %0.3f s" % (time2-time1), self.verbose)

    ############################################################## TIMING



    features_dict = self._transform(protein_pdb_file, ligand_pdb_file)

    features_dict = self._transform(protein_pdb_file, ligand_file)

    shutil.rmtree(tempdir)

    return features_dict.values()



  def featurize_complexes(self, mol_pdbs, protein_pdbs, log_every_n=1000):

  def featurize_complexes(self, mol_files, protein_pdbs, log_every_n=1000):

    """

    Calculate features for mol/protein complexes.



    Parameters

    ----------

    mol_pdbs: list

      List of PDBs for molecules. Each PDB should be a list of lines of the

      PDB file.

    mols: list

      List of PDB filenames for molecules.

    protein_pdbs: list

      List of PDBs for proteins. Each PDB should be a list of lines of the

      PDB file.

      List of PDB filenames for proteins.

    """

    features = []

    for i, (mol_pdb, protein_pdb) in enumerate(zip(mol_pdbs, protein_pdbs)):

    for i, (mol_file, protein_pdb) in enumerate(zip(mol_files, protein_pdbs)):

      if i % log_every_n == 0:

        log("Featurizing %d / %d" % (i, len(mol_pdbs)))

      features += self._featurize_complex(mol_pdb, protein_pdb)

        log("Featurizing %d / %d" % (i, len(mol_files)))

      ligand_ext = get_ligand_filetype(mol_file)

      with open(mol_file) as mol_f:

        mol_lines = mol_f.readlines()

      with open(protein_pdb) as protein_file:

        protein_pdb_lines = protein_file.readlines()

      features += self._featurize_complex(ligand_ext, mol_lines,

                                          protein_pdb_lines)

    features = np.asarray(features)

    return features



  def _transform(self, protein_pdb, ligand_pdb):

  def _transform(self, protein_pdb, ligand_file):

    """Computes featurization of protein/ligand complex.



    Takes as input files (strings) for pdb of the protein, pdb of the ligand,
@@ -947,6 +968,12 @@ class GridFeaturizer(ComplexFeaturizer): 
    if not self.ligand_only:

      protein_xyz, protein_ob = load_molecule(

          protein_pdb, calc_charges=False)

      ############################################################ DEBUG

      #print("protein_pdb")

      #print(protein_pdb)

      #print("protein_xyz")

      #print(protein_xyz)

      ############################################################ DEBUG

    ############################################################## TIMING

    time2 = time.time()

    log("TIMING: Loading protein coordinates took %0.3f s" % (time2-time1),
@@ -956,7 +983,7 @@ class GridFeaturizer(ComplexFeaturizer): 
    time1 = time.time()

    ############################################################## TIMING

    ligand_xyz, ligand_ob = load_molecule(

        ligand_pdb, calc_charges=False)

        ligand_file, calc_charges=False)

    ############################################################## TIMING

    time2 = time.time()

    log("TIMING: Loading ligand coordinates took %0.3f s" % (time2-time1),
@@ -999,6 +1026,10 @@ class GridFeaturizer(ComplexFeaturizer): 
          _featurize_binding_pocket_ecfp(

              protein_xyz, protein_ob, ligand_xyz, ligand_ob,

              pairwise_distances, cutoff=4.5, ecfp_degree=self.ecfp_degree))

      ################################################################ DEBUG

      #print("protein_ecfp_dict")

      #print(protein_ecfp_dict)

      ################################################################ DEBUG

      ############################################################## TIMING

      time2 = time.time()

      log("TIMING: ecfp voxel computataion took %0.3f s" % (time2-time1),

examples/pdbbind/get_pdbbind.sh

0 → 100755
+4 −0
Original line number Diff line number Diff line 
echo "Pulling pdbbind dataset from deepchem"

wget http://deepchem.io.s3-website-us-west-1.amazonaws.com/datasets/pdbbind_v2015.tar.gz

echo "Extracting pdbbind structures"

tar -zxvf pdbbind_v2015.tar.gz

examples/pdbbind/pdbbind_datasets.py

+45 −92
Original line number Diff line number Diff line
@@ -11,13 +11,7 @@ import numpy as np 
import pandas as pd

import shutil

from rdkit import Chem

from deepchem.utils.save import load_from_disk

from deepchem.data import DiskDataset

from deepchem.featurizers.fingerprints import CircularFingerprint

from deepchem.trans import BalancingTransformer

from deepchem.featurizers.nnscore import NNScoreComplexFeaturizer

from deepchem.featurizers.grid_featurizer import GridFeaturizer

from deepchem.featurizers.atomic_coordinates import NeighborListComplexAtomicCoordinates

import deepchem as dc



def load_pdbbind_labels(labels_file):

  """Loads pdbbind labels as dataframe"""
@@ -34,52 +28,33 @@ def load_pdbbind_labels(labels_file): 
               "ignore-this-field", "reference", "ligand name"))

  return contents_df



def compute_pdbbind_grid_feature(compound_featurizers, complex_featurizers,

                                 pdb_subdir, pdb_code):

def compute_pdbbind_features(grid_featurizer, pdb_subdir, pdb_code):

  """Compute features for a given complex"""

  protein_file = os.path.join(pdb_subdir, "%s_protein.pdb" % pdb_code)

  ligand_file = os.path.join(pdb_subdir, "%s_ligand.sdf" % pdb_code)

  rdkit_mol = next(Chem.SDMolSupplier(str(ligand_file)))



  all_features = []

  for complex_featurizer in complex_featurizers:

    features = complex_featurizer.featurize_complexes(

  ################################################################ DEBUG

  #print("protein_file")

  #print(protein_file)

  #print("ligand_file")

  #print(ligand_file)

  ################################################################ DEBUG

  features = grid_featurizer.featurize_complexes(

    [ligand_file], [protein_file])

    all_features.append(np.squeeze(features))

  

  for compound_featurizer in compound_featurizers:

    features = np.squeeze(compound_featurizer.featurize([rdkit_mol]))

    all_features.append(features)



  features = np.concatenate(all_features)

  features = np.squeeze(features)

  return features



def compute_pdbbind_coordinate_features(

    complex_featurizer, pdb_subdir, pdb_code):

  """Compute features for a given complex"""

  protein_file = os.path.join(pdb_subdir, "%s_protein.pdb" % pdb_code)

  ligand_file = os.path.join(pdb_subdir, "%s_ligand.sdf" % pdb_code)



  feature = complex_featurizer.featurize_complexes(

    [ligand_file], [protein_file])

  return feature



def load_core_pdbbind_coordinates(pdbbind_dir, base_dir, reload=True):

def load_core_pdbbind_grid(split="index", feat="grid"):

  """Load PDBBind datasets. Does not do train/test split"""

  # Set some global variables up top

  regen = False

  neighbor_cutoff = 4

  max_num_neighbors = 10



  # Create some directories for analysis

  # The base_dir holds the results of all analysis

  current_dir = os.path.dirname(os.path.realpath(__file__))

  pdbbind_dir = os.path.join(current_dir, "v2015")

  #Make directories to store the raw and featurized datasets.

  data_dir = os.path.join(base_dir, "dataset")



  # Load PDBBind dataset

  labels_file = os.path.join(pdbbind_dir, "INDEX_core_data.2013")

  pdb_subdirs = os.path.join(pdbbind_dir, "website-core-set")

  tasks = ["-logKd/Ki"]

  print("About to load contents.")

  contents_df = load_pdbbind_labels(labels_file)
@@ -87,62 +62,23 @@ def load_core_pdbbind_coordinates(pdbbind_dir, base_dir, reload=True): 
  y = np.array([float(val) for val in contents_df["-logKd/Ki"].values])



  # Define featurizers

  featurizer = NeighborListComplexAtomicCoordinates(

      max_num_neighbors, neighbor_cutoff)

  

  # Featurize Dataset

  features = []

  for ind, pdb_code in enumerate(ids):

    print("Processing %s" % str(pdb_code))

    pdb_subdir = os.path.join(pdb_subdirs, pdb_code)

    computed_feature = compute_pdbbind_coordinate_features(

        featurizer, pdb_subdir, pdb_code)

    features.append(computed_feature)

  X = np.array(features, dtype-object)

  w = np.ones_like(y)

   

  dataset = DiskDataset.from_numpy(data_dir, X, y, w, ids)

  transformers = []

  

  return tasks, dataset, transformers



def load_core_pdbbind_grid(pdbbind_dir, base_dir, reload=True):

  """Load PDBBind datasets. Does not do train/test split"""

  # Set some global variables up top

  regen = False



  # Create some directories for analysis

  # The base_dir holds the results of all analysis

  if not reload:

    if os.path.exists(base_dir):

      shutil.rmtree(base_dir)

  if not os.path.exists(base_dir):

    os.makedirs(base_dir)

  current_dir = os.path.dirname(os.path.realpath(__file__))

  #Make directories to store the raw and featurized datasets.

  data_dir = os.path.join(base_dir, "dataset")



  # Load PDBBind dataset

  labels_file = os.path.join(pdbbind_dir, "INDEX_core_data.2013")

  pdb_subdirs = os.path.join(pdbbind_dir, "website-core-set")

  tasks = ["-logKd/Ki"]

  print("About to load contents.")

  contents_df = load_pdbbind_labels(labels_file)

  ids = contents_df["PDB code"].values

  y = np.array([float(val) for val in contents_df["-logKd/Ki"].values])



  # Define featurizers

  grid_featurizer = GridFeaturizer(

  if feat == "grid":

    featurizer = dc.feat.GridFeaturizer(

        voxel_width=16.0, feature_types="voxel_combined",

        # TODO(rbharath, enf): Figure out why pi_stack is slow and cation_pi

        # causes segfaults.

        #voxel_feature_types=["ecfp", "splif", "hbond", "pi_stack", "cation_pi",

        #"salt_bridge"], ecfp_power=9, splif_power=9,

      voxel_feature_types=["ecfp", "splif", "hbond", 

      "salt_bridge"], ecfp_power=9, splif_power=9,

        voxel_feature_types=["ecfp", "splif", "hbond", "salt_bridge"],

        ecfp_power=9, splif_power=9,

        parallel=True, flatten=True)

  compound_featurizers = [CircularFingerprint(size=1024)]

  complex_featurizers = [grid_featurizer]

  elif feat == "coord":

    neighbor_cutoff = 4

    max_num_neighbors = 10

    featurizer = dc.feat.NeighborListComplexAtomicCoordinates(

        max_num_neighbors, neighbor_cutoff)

  else:

    raise ValueError("feat not defined.")

  

  # Featurize Dataset

  features = []
@@ -150,9 +86,13 @@ def load_core_pdbbind_grid(pdbbind_dir, base_dir, reload=True): 
  y_inds = []

  for ind, pdb_code in enumerate(ids):

    print("Processing %s" % str(pdb_code))

    pdb_subdir = os.path.join(pdb_subdirs, pdb_code)

    computed_feature = compute_pdbbind_grid_feature(

        compound_featurizers, complex_featurizers, pdb_subdir, pdb_code)

    pdb_subdir = os.path.join(pdbbind_dir, pdb_code)

    ######################################################## DEBUG

    #print("pdb_subdir, pdb_code")

    #print(pdb_subdir, pdb_code)

    ######################################################## DEBUG

    computed_feature = compute_pdbbind_features(

        featurizer, pdb_subdir, pdb_code)

    if feature_len is None:

      feature_len = len(computed_feature)

    if len(computed_feature) != feature_len:
@@ -160,11 +100,24 @@ def load_core_pdbbind_grid(pdbbind_dir, base_dir, reload=True): 
      continue

    y_inds.append(ind)

    features.append(computed_feature)

    ######################################################## DEBUG

    #print("np.count_nonzero(computed_feature)")

    #print(np.count_nonzero(computed_feature))

    #print("computed_feature")

    #print(computed_feature)

    ##assert 0 == 1

    ######################################################## DEBUG

  y = y[y_inds]

  X = np.vstack(features)

  w = np.ones_like(y)

   

  dataset = DiskDataset.from_numpy(data_dir, X, y, w, ids)

  dataset = dc.data.DiskDataset.from_numpy(X, y, w, ids)

  transformers = []



  return tasks, dataset, transformers
 
  splitters = {'index': dc.splits.IndexSplitter(),

               'random': dc.splits.RandomSplitter(),

               'scaffold': dc.splits.ScaffoldSplitter()}

  splitter = splitters[split]

  train, valid, test = splitter.train_valid_test_split(dataset)

  

  return tasks, (train, valid, test), transformers

examples/pdbbind/pdbbind_tf.py

+13 −60
Original line number Diff line number Diff line 
"""

Script that trains Sklearn models on PDBbind dataset.

"""



from __future__ import print_function

from __future__ import division

from __future__ import unicode_literals
@@ -10,79 +9,33 @@ __author__ = "Bharath Ramsundar" 
__copyright__ = "Copyright 2016, Stanford University"

__license__ = "GPL"



import os

import sys

import tempfile

import deepchem as dc

import numpy as np

import numpy.random



import sys

import shutil

from pdbbind_datasets import load_core_pdbbind_grid

from deepchem.featurizers.featurize import DataLoader

from deepchem.hyper import HyperparamOpt

from deepchem import metrics

from deepchem.metrics import Metric

from deepchem.models.tensorflow_models import TensorflowModel

from deepchem.models.tensorflow_models.fcnet import TensorflowMultiTaskRegressor

from deepchem.utils.evaluate import Evaluator

from deepchem.splits import RandomSplitter

from deepchem.featurizers.atomic_coordinates import AtomicCoordinates

from deepchem.data import DiskDataset



verbosity = "high"

base_dir = "/tmp/PDBBIND-ATOMICNET"

if os.path.exists(base_dir):

  shutil.rmtree(base_dir)

os.makedirs(base_dir)



feature_dir = os.path.join(base_dir, "feature")

train_dir = os.path.join(base_dir, "train")

valid_dir = os.path.join(base_dir, "valid")

test_dir = os.path.join(base_dir, "test")

model_dir = os.path.join(base_dir, "model")



# REPLACE WITH DOWNLOADED PDBBIND EXAMPLE

pdbbind_dir = "/tmp/deep-docking/datasets/pdbbind"

pdbbind_tasks, dataset, transformers = load_core_pdbbind_grid(

    pdbbind_dir, base_dir)

# For stable runs 

np.random.seed(123)



print("About to perform train/valid/test split.")

num_train = .8 * len(dataset)

X, y, w, ids = (dataset.X, dataset.y, dataset.w, dataset.ids)

pdbbind_tasks, pdbbind_datasets, transformers = load_core_pdbbind_grid()

train_dataset, valid_dataset, test_dataset = pdbbind_datasets 



X_train, X_valid = X[:num_train], X[num_train:]

y_train, y_valid = y[:num_train], y[num_train:]

w_train, w_valid = w[:num_train], w[num_train:]

ids_train, ids_valid = ids[:num_train], ids[num_train:]

metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)



train_dataset = DiskDataset.from_numpy(train_dir, X_train, y_train,

                                   w_train, ids_train, pdbbind_tasks)

valid_dataset = DiskDataset.from_numpy(valid_dir, X_valid, y_valid,

                                   w_valid, ids_valid, pdbbind_tasks)



classification_metric = Metric(metrics.pearson_r2_score, verbosity=verbosity,

                               mode="regression")



n_features = dataset.get_data_shape()[0]

tensorflow_model = TensorflowMultiTaskRegressor(

    len(pdbbind_tasks), n_features, model_dir, dropouts=[.25],

    learning_rate=0.0003, weight_init_stddevs=[.1],

    batch_size=64, verbosity=verbosity)

model = TensorflowModel(tensorflow_model, model_dir)

n_features = train_dataset.X.shape[1]

model = dc.models.TensorflowMultiTaskRegressor(

    len(pdbbind_tasks), n_features, dropouts=[.25], learning_rate=0.0003,

    weight_init_stddevs=[.1], batch_size=64)



# Fit trained model

model.fit(train_dataset, nb_epoch=20)

model.save()



train_evaluator = Evaluator(model, train_dataset, transformers, verbosity=verbosity)

train_scores = train_evaluator.compute_model_performance([classification_metric])

print("Evaluating model")

train_scores = model.evaluate(train_dataset, [metric], transformers)

valid_scores = model.evaluate(valid_dataset, [metric], transformers)



print("Train scores")

print(train_scores)



valid_evaluator = Evaluator(model, valid_dataset, transformers, verbosity=verbosity)

valid_scores = valid_evaluator.compute_model_performance([classification_metric])



print("Validation scores")

print(valid_scores)

CRA Git | Maintained and supported by SUSTech CRA and CCSE