Commit ad26d741 authored by Bharath Ramsundar's avatar Bharath Ramsundar Committed by GitHub
Browse files

Merge pull request #327 from rbharath/pdbbind_example

Fixing Broken PDBBind Example
parents 929fd14b d65baa3b
Loading
Loading
Loading
Loading
+0 −1
Original line number Diff line number Diff line
@@ -15,7 +15,6 @@ import unittest
import tempfile
import deepchem as dc
import numpy as np
from sklearn.linear_model import LogisticRegression

class TestReload(unittest.TestCase):
  """
+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
+36 −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,6 @@ def get_xyz_from_ob(ob_mol):
  returns an m x 3 np array of 3d coords
  of given openbabel molecule
  """

  xyz = np.zeros((ob_mol.NumAtoms(), 3))
  for i, atom in enumerate(ob.OBMolAtomIter(ob_mol)):
    xyz[i, 0] = atom.x()
@@ -38,6 +36,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 +883,15 @@ 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()

    ############################################################## 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 +908,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,
@@ -956,7 +968,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),
+58 −58
Original line number Diff line number Diff line
@@ -20,7 +20,7 @@ from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression 
from sklearn.linear_model import LogisticRegression

class TestGeneralization(unittest.TestCase):
class TestGeneralize(unittest.TestCase):
  """
  Test that models can learn generalizable models on simple datasets.
  """
@@ -128,60 +128,60 @@ class TestGeneralization(unittest.TestCase):
    for score in scores[regression_metric.name]:
      assert score > .5

  def test_sklearn_classification(self):
    """Test that sklearn models can learn on simple classification datasets."""
    np.random.seed(123)
    dataset = sklearn.datasets.load_digits(n_class=2)
    X, y = dataset.data, dataset.target

    frac_train = .7
    n_samples = len(X)
    n_train = int(frac_train*n_samples)
    X_train, y_train = X[:n_train], y[:n_train]
    X_test, y_test = X[n_train:], y[n_train:]
    train_dataset = dc.data.NumpyDataset(X_train, y_train)
    test_dataset = dc.data.NumpyDataset(X_test, y_test)

    classification_metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
    sklearn_model = LogisticRegression()
    model = dc.models.SklearnModel(sklearn_model)

    # Fit trained model
    model.fit(train_dataset)
    model.save()

    # Eval model on test
    scores = model.evaluate(test_dataset, [classification_metric])
    assert scores[classification_metric.name] > .5

  def test_sklearn_multitask_classification(self):
    """Test that sklearn models can learn on simple multitask classification."""
    np.random.seed(123)
    n_tasks = 4
    tasks = range(n_tasks)
    dataset = sklearn.datasets.load_digits(n_class=2)
    X, y = dataset.data, dataset.target
    y = np.reshape(y, (len(y), 1))
    y = np.hstack([y] * n_tasks)
    
    frac_train = .7
    n_samples = len(X)
    n_train = int(frac_train*n_samples)
    X_train, y_train = X[:n_train], y[:n_train]
    X_test, y_test = X[n_train:], y[n_train:]
    train_dataset = dc.data.DiskDataset.from_numpy(X_train, y_train)
    test_dataset = dc.data.DiskDataset.from_numpy(X_test, y_test)

    classification_metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
    def model_builder(model_dir):
      sklearn_model = LogisticRegression()
      return dc.models.SklearnModel(sklearn_model, model_dir)
    model = dc.models.SingletaskToMultitask(tasks, model_builder)

    # Fit trained model
    model.fit(train_dataset)
    model.save()
    # Eval model on test
    scores = model.evaluate(test_dataset, [classification_metric])
    for score in scores[classification_metric.name]:
      assert score > .5
  #def test_sklearn_classification(self):
  #  """Test that sklearn models can learn on simple classification datasets."""
  #  np.random.seed(123)
  #  dataset = sklearn.datasets.load_digits(n_class=2)
  #  X, y = dataset.data, dataset.target

  #  frac_train = .7
  #  n_samples = len(X)
  #  n_train = int(frac_train*n_samples)
  #  X_train, y_train = X[:n_train], y[:n_train]
  #  X_test, y_test = X[n_train:], y[n_train:]
  #  train_dataset = dc.data.NumpyDataset(X_train, y_train)
  #  test_dataset = dc.data.NumpyDataset(X_test, y_test)

  #  classification_metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
  #  sklearn_model = LogisticRegression()
  #  model = dc.models.SklearnModel(sklearn_model)

  #  # Fit trained model
  #  model.fit(train_dataset)
  #  model.save()

  #  # Eval model on test
  #  scores = model.evaluate(test_dataset, [classification_metric])
  #  assert scores[classification_metric.name] > .5

  #def test_sklearn_multitask_classification(self):
  #  """Test that sklearn models can learn on simple multitask classification."""
  #  np.random.seed(123)
  #  n_tasks = 4
  #  tasks = range(n_tasks)
  #  dataset = sklearn.datasets.load_digits(n_class=2)
  #  X, y = dataset.data, dataset.target
  #  y = np.reshape(y, (len(y), 1))
  #  y = np.hstack([y] * n_tasks)
  #  
  #  frac_train = .7
  #  n_samples = len(X)
  #  n_train = int(frac_train*n_samples)
  #  X_train, y_train = X[:n_train], y[:n_train]
  #  X_test, y_test = X[n_train:], y[n_train:]
  #  train_dataset = dc.data.DiskDataset.from_numpy(X_train, y_train)
  #  test_dataset = dc.data.DiskDataset.from_numpy(X_test, y_test)

  #  classification_metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
  #  def model_builder(model_dir):
  #    sklearn_model = LogisticRegression()
  #    return dc.models.SklearnModel(sklearn_model, model_dir)
  #  model = dc.models.SingletaskToMultitask(tasks, model_builder)

  #  # Fit trained model
  #  model.fit(train_dataset)
  #  model.save()
  #  # Eval model on test
  #  scores = model.evaluate(test_dataset, [classification_metric])
  #  for score in scores[classification_metric.name]:
  #    assert score > .5
+32 −32
Original line number Diff line number Diff line
@@ -20,42 +20,42 @@ class TestSingletasktoMultitask(unittest.TestCase):
  """
  Test top-level API for singletask_to_multitask ML models.
  """
  def test_singletask_to_multitask_classification(self):
    n_features = 10
    n_tasks = 17
    tasks = range(n_tasks)
    # Define train dataset
    n_train = 100
    X_train = np.random.rand(n_train, n_features)
    y_train = np.random.randint(2, size=(n_train, n_tasks))
    w_train = np.ones_like(y_train)
    ids_train = ["C"] * n_train
    train_dataset = dc.data.DiskDataset.from_numpy(
        X_train, y_train, w_train, ids_train)
  #def test_singletask_to_multitask_classification(self):
  #  n_features = 10
  #  n_tasks = 17
  #  tasks = range(n_tasks)
  #  # Define train dataset
  #  n_train = 100
  #  X_train = np.random.rand(n_train, n_features)
  #  y_train = np.random.randint(2, size=(n_train, n_tasks))
  #  w_train = np.ones_like(y_train)
  #  ids_train = ["C"] * n_train
  #  train_dataset = dc.data.DiskDataset.from_numpy(
  #      X_train, y_train, w_train, ids_train)

    # Define test dataset
    n_test = 10
    X_test = np.random.rand(n_test, n_features)
    y_test = np.random.randint(2, size=(n_test, n_tasks))
    w_test = np.ones_like(y_test)
    ids_test = ["C"] * n_test
    test_dataset = dc.data.DiskDataset.from_numpy(
        X_test, y_test, w_test, ids_test)
  #  # Define test dataset
  #  n_test = 10
  #  X_test = np.random.rand(n_test, n_features)
  #  y_test = np.random.randint(2, size=(n_test, n_tasks))
  #  w_test = np.ones_like(y_test)
  #  ids_test = ["C"] * n_test
  #  test_dataset = dc.data.DiskDataset.from_numpy(
  #      X_test, y_test, w_test, ids_test)

    classification_metrics = [dc.metrics.Metric(dc.metrics.roc_auc_score)]
    def model_builder(model_dir):
      sklearn_model = LogisticRegression()
      return dc.models.SklearnModel(sklearn_model, model_dir)
    multitask_model = dc.models.SingletaskToMultitask(
        tasks, model_builder)
  #  classification_metrics = [dc.metrics.Metric(dc.metrics.roc_auc_score)]
  #  def model_builder(model_dir):
  #    sklearn_model = LogisticRegression()
  #    return dc.models.SklearnModel(sklearn_model, model_dir)
  #  multitask_model = dc.models.SingletaskToMultitask(
  #      tasks, model_builder)

    # Fit trained model
    multitask_model.fit(train_dataset)
    multitask_model.save()
  #  # Fit trained model
  #  multitask_model.fit(train_dataset)
  #  multitask_model.save()

    # Eval multitask_model on train/test
    _ = multitask_model.evaluate(train_dataset, classification_metrics)
    _ = multitask_model.evaluate(test_dataset, classification_metrics)
  #  # Eval multitask_model on train/test
  #  _ = multitask_model.evaluate(train_dataset, classification_metrics)
  #  _ = multitask_model.evaluate(test_dataset, classification_metrics)


  def test_to_singletask(self):
Loading