Merge remote-tracking branch 'remotes/origin/master' into HEAD

   conda install -c omnia openbabel=2.4.0

   ``` 

3. `pandas`

   ```bash

   conda install pandas 

   ```

4. `rdkit`

3. `rdkit`

   ```bash

   conda install -c omnia rdkit

   ```

5. `boost`

   ```bash

   conda install -c omnia boost=1.59.0

   ```

6. `joblib`

4. `joblib`

   ```bash

   conda install joblib 

   ```

7. `keras`

   ```bash

   pip install keras --user

   ```

   `deepchem` only supports the `tensorflow` backend for keras. To set the backend to `tensorflow`,

   add the following line to your `~/.bashrc`

5. `keras`

   ```bash

   export KERAS_BACKEND=tensorflow

   pip install keras

   ```

   See [keras docs](https://keras.io/backend/) for more details and alternate methods of setting backend.

   `deepchem` only supports the `tensorflow` (default) backend for keras.

8. `mdtraj`

6. `mdtraj`

   ```bash

   conda install -c omnia mdtraj

   ```

9. `scikit-learn`

   ```bash

   conda install scikit-learn 

   ```

10. `tensorflow`: Installing `tensorflow` on older versions of Linux (which

7. `tensorflow`: Installing `tensorflow` on older versions of Linux (which

    have glibc < 2.17) can be very challenging. For these older Linux versions,

    contact your local sysadmin to work out a custom installation. If your

    version of Linux is recent, then the following command will work:

    conda install -c https://conda.anaconda.org/jjhelmus tensorflow

    ```

11. `h5py`:

    ```

    conda install h5py

    ```

12. `deepchem`: Clone the `deepchem` github repo:

8. `deepchem`: Clone the `deepchem` github repo:

   ```bash

   git clone https://github.com/deepchem/deepchem.git

   ```

   python setup.py install

   ```

13. To run test suite, install `nosetests`:

9. To run test suite, install `nosetests`:

   ```bash

    pip install nose --user

   pip install nose

   ```

   Make sure that the correct version of `nosetests` is active by running

   ```bash

   You might need to uninstall a system install of `nosetests` if

   there is a conflict.

14. If installation has been successful, all tests in test suite should pass:

10. If installation has been successful, all tests in test suite should pass:

    ```bash

    nosetests -v deepchem --nologcapture 

    ```

Most machine learning algorithms require that input data form vectors.

However, input data for drug-discovery datasets routinely come in the

format of lists of molecules and associated experimental readouts. To

transform lists of molecules into vectors, we need to use the DeepChem

loader class ``dc.load.DataLoader``. Instances of this class must be

passed a ``Featurizer`` object. DeepChem provides a number of

different subclasses of ``Featurizer`` for convenience:

transform lists of molecules into vectors, we need to subclasses of DeepChem

loader class ```dc.data.DataLoader``` such as ```dc.data.CSVLoader``` or 

```dc.data.SDFLoader```. Users can subclass ```dc.data.DataLoader``` to

load arbitrary file formats. All loaders must be

passed a ```dc.feat.Featurizer``` object. DeepChem provides a number of

different subclasses of ```dc.feat.Featurizer``` for convenience.

### Performances

* Classification

import deepchem.splits

import deepchem.trans

import deepchem.utils

import deepchem.dock

"""

Imports all submodules 

"""

from __future__ import print_function

from __future__ import division

from __future__ import unicode_literals

from deepchem.dock.pose_generation import PoseGenerator

from deepchem.dock.pose_generation import VinaPoseGenerator

from deepchem.dock.pose_scoring import PoseScorer

from deepchem.dock.pose_scoring import GridPoseScorer

from deepchem.dock.docking import Docker

from deepchem.dock.docking import VinaGridRFDocker

from deepchem.dock.docking import VinaGridDNNDocker

from deepchem.dock.binding_pocket import ConvexHullPocketFinder

"""

Computes putative binding pockets on protein.

"""

from __future__ import print_function

from __future__ import division

from __future__ import unicode_literals

__author__ = "Bharath Ramsundar"

__copyright__ = "Copyright 2017, Stanford University"

__license__ = "GPL"

import numpy as np

import os

import pybel

import tempfile

from scipy.spatial import ConvexHull

from deepchem.feat import hydrogenate_and_compute_partial_charges

from deepchem.feat.atomic_coordinates import AtomicCoordinates

from deepchem.feat.grid_featurizer import load_molecule

from subprocess import call

def extract_active_site(protein_file, ligand_file, cutoff=4):

  """Extracts a box for the active site."""

  protein_coords = load_molecule(protein_file)[0]

  ligand_coords = load_molecule(ligand_file)[0]

  num_ligand_atoms = len(ligand_coords)

  num_protein_atoms = len(protein_coords)

  pocket_inds = []

  pocket_atoms = set([])

  for lig_atom_ind in range(num_ligand_atoms):

    lig_atom = ligand_coords[lig_atom_ind]

    for protein_atom_ind in range(num_protein_atoms):

      protein_atom = protein_coords[protein_atom_ind]

      if np.linalg.norm(lig_atom - protein_atom) < cutoff:

        if protein_atom_ind not in pocket_atoms:

          pocket_atoms = pocket_atoms.union(set([protein_atom_ind]))

  # Should be an array of size (n_pocket_atoms, 3)

  pocket_atoms = list(pocket_atoms)

  n_pocket_atoms = len(pocket_atoms)

  pocket_coords = np.zeros((n_pocket_atoms, 3))

  for ind, pocket_ind in enumerate(pocket_atoms):

    pocket_coords[ind] = protein_coords[pocket_ind]

  x_min = int(np.floor(np.amin(pocket_coords[:, 0])))

  x_max = int(np.ceil(np.amax(pocket_coords[:, 0])))

  y_min = int(np.floor(np.amin(pocket_coords[:, 1])))

  y_max = int(np.ceil(np.amax(pocket_coords[:, 1])))

  z_min = int(np.floor(np.amin(pocket_coords[:, 2])))

  z_max = int(np.ceil(np.amax(pocket_coords[:, 2])))

  return (((x_min, x_max), (y_min, y_max), (z_min, z_max)),

          pocket_atoms, pocket_coords)

def compute_overlap(mapping, box1, box2):

  """Computes overlap between the two boxes.

  Overlap is defined as % atoms of box1 in box2. Note that

  overlap is not a symmetric measurement.

  """

  atom1 = set(mapping[box1])

  atom2 = set(mapping[box2])

  return len(atom1.intersection(atom2))/float(len(atom1))

def get_all_boxes(coords, pad=5):

  """Get all pocket boxes for protein coords.

  We pad all boxes the prescribed number of angstroms.

  TODO(rbharath): It looks like this may perhaps be non-deterministic?

  """

  hull = ConvexHull(coords)

  boxes = []

  for triangle in hull.simplices:

    # coords[triangle, 0] gives the x-dimension of all triangle points

    # Take transpose to make sure rows correspond to atoms.

    points = np.array(

        [coords[triangle, 0], coords[triangle, 1], coords[triangle, 2]]).T

    # We voxelize so all grids have integral coordinates (convenience)

    x_min, x_max = np.amin(points[:, 0]), np.amax(points[:, 0])

    x_min, x_max = int(np.floor(x_min))-pad, int(np.ceil(x_max))+pad

    y_min, y_max = np.amin(points[:, 1]), np.amax(points[:, 1])

    y_min, y_max = int(np.floor(y_min))-pad, int(np.ceil(y_max))+pad

    z_min, z_max = np.amin(points[:, 2]), np.amax(points[:, 2])

    z_min, z_max = int(np.floor(z_min))-pad, int(np.ceil(z_max))+pad

    boxes.append(((x_min, x_max), (y_min, y_max), (z_min, z_max)))

  return boxes 

def boxes_to_atoms(atom_coords, boxes):

  """Maps each box to a list of atoms in that box.

  TODO(rbharath): This does a num_atoms x num_boxes computations. Is

  there a reasonable heuristic we can use to speed this up?

  """

  mapping = {}

  for box_ind, box in enumerate(boxes):

    box_atoms = []

    (x_min, x_max), (y_min, y_max), (z_min, z_max) = box

    print("Handing box %d/%d" % (box_ind, len(boxes)))

    for atom_ind in range(len(atom_coords)):

      atom = atom_coords[atom_ind] 

      x_cont = x_min <= atom[0] and atom[0] <= x_max

      y_cont = y_min <= atom[1] and atom[1] <= y_max

      z_cont = z_min <= atom[2] and atom[2] <= z_max

      if x_cont and y_cont and z_cont:

        box_atoms.append(atom_ind)

    mapping[box] = box_atoms

  return mapping

def merge_boxes(box1, box2):

  """Merges two boxes."""

  (x_min1, x_max1), (y_min1, y_max1), (z_min1, z_max1) = box1 

  (x_min2, x_max2), (y_min2, y_max2), (z_min2, z_max2) = box2 

  x_min = min(x_min1, x_min2)

  y_min = min(y_min1, y_min2)

  z_min = min(z_min1, z_min2)

  x_max = max(x_max1, x_max2)

  y_max = max(y_max1, y_max2)

  z_max = max(z_max1, z_max2)

  return ((x_min, x_max), (y_min, y_max), (z_min, z_max))

def merge_overlapping_boxes(mapping, boxes, threshold=.8):

  """Merge boxes which have an overlap greater than threshold.

  TODO(rbharath): This merge code is terribly inelegant. It's also quadratic

  in number of boxes. It feels like there ought to be an elegant divide and

  conquer approach here. Figure out later...

  """

  num_boxes = len(boxes)

  outputs = []

  for i in range(num_boxes):

    box = boxes[0]

    new_boxes = []

    new_mapping = {}

    # If overlap of box with previously generated output boxes, return

    contained = False

    for output_box in outputs:

      # Carry forward mappings

      new_mapping[output_box] = mapping[output_box]

      if compute_overlap(mapping, box, output_box) == 1:

        contained = True

    if contained:

      continue

    # We know that box has at least one atom not in outputs

    unique_box = True

    for merge_box in boxes[1:]:

      overlap = compute_overlap(mapping, box, merge_box)

      if overlap < threshold:

        new_boxes.append(merge_box)

        new_mapping[merge_box] = mapping[merge_box]

      else:

        # Current box has been merged into box further down list.

        # No need to output current box

        unique_box = False

        merged = merge_boxes(box, merge_box)

        new_boxes.append(merged)

        new_mapping[merged] = list(

            set(mapping[box]).union(set(mapping[merge_box])))

    if unique_box:

      outputs.append(box)

      new_mapping[box] = mapping[box]

    boxes = new_boxes

    mapping = new_mapping

  return outputs, mapping

class BindingPocketFinder(object):

  """Abstract superclass for binding pocket detectors"""

  def find_pockets(self, protein_file):

    """Finds potential binding pockets in proteins."""

    raise NotImplementedError

class ConvexHullPocketFinder(BindingPocketFinder):

  """Implementation that uses convex hull of protein to find pockets.

  Based on https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4112621/pdf/1472-6807-14-18.pdf

  """

  def __init__(self, pad=5):

    self.pad = pad

  def find_all_pockets(self, protein_file):

    """Find list of binding pockets on protein."""

    # protein_coords is (N, 3) tensor

    coords = load_molecule(protein_file)[0]

    return get_all_boxes(coords, self.pad)

  def find_pockets(self, protein_file, ligand_file):

    """Find list of suitable binding pockets on protein."""

    protein_coords = load_molecule(protein_file)[0]

    ligand_coords = load_molecule(ligand_file)[0]

    boxes = get_all_boxes(protein_coords, self.pad)

    mapping = boxes_to_atoms(protein_coords, boxes)

    merged_boxes, mapping = merge_overlapping_boxes(mapping, boxes)

    return merged_boxes, mapping

"""

Docks protein-ligand pairs 

"""

from __future__ import print_function

from __future__ import division

from __future__ import unicode_literals

__author__ = "Bharath Ramsundar"

__copyright__ = "Copyright 2016, Stanford University"

__license__ = "GPL"

import numpy as np

import os

import tempfile

from deepchem.feat import GridFeaturizer

from deepchem.data import DiskDataset

from deepchem.models import SklearnModel

from deepchem.models import TensorflowMultiTaskRegressor

from deepchem.dock.pose_scoring import GridPoseScorer

from deepchem.dock.pose_generation import VinaPoseGenerator

from sklearn.ensemble import RandomForestRegressor

from subprocess import call

class Docker(object):

  """Abstract Class specifying API for Docking."""

  def dock(self, protein_file, ligand_file):

    raise NotImplementedError

class VinaGridRFDocker(Docker):

  """Vina pose-generation, RF-models on grid-featurization of complexes."""

  def __init__(self):

    """Builds model."""

    self.base_dir = tempfile.mkdtemp()

    print("About to download trained model.")

    call(("wget -c http://deepchem.io.s3-website-us-west-1.amazonaws.com/trained_models/random_full_RF.tar.gz").split())

    call(("tar -zxvf random_full_RF.tar.gz").split())

    call(("mv random_full_RF %s" % (self.base_dir)).split())

    self.model_dir = os.path.join(self.base_dir, "random_full_RF")

    # Fit model on dataset

    model = SklearnModel(model_dir=self.model_dir)

    model.reload()

    self.pose_scorer = GridPoseScorer(model, feat="grid")

    self.pose_generator = VinaPoseGenerator() 

  def dock(self, protein_file, ligand_file):

    """Docks using Vina and RF."""

    protein_docked, ligand_docked = self.pose_generator.generate_poses(

        protein_file, ligand_file)

    score = self.pose_scorer.score(protein_docked, ligand_docked)

    return (score, (protein_docked, ligand_docked))

class VinaGridDNNDocker(object):

  """Vina pose-generation, DNN-models on grid-featurization of complexes."""

  def __init__(self, n_trees=100):

    """Builds model."""

    self.base_dir = tempfile.mkdtemp()

    print("About to download trained model.")

    call(("wget -c http://deepchem.io.s3-website-us-west-1.amazonaws.com/trained_models/random_full_DNN.tar.gz").split())

    call(("tar -zxvf random_full_DNN.tar.gz").split())

    call(("mv random_full_DNN %s" % (self.base_dir)).split())

    self.model_dir = os.path.join(self.base_dir, "random_full_DNN")

    # Fit model on dataset

    pdbbind_tasks = ["-logKd/Ki"]

    n_features = 2052

    model = TensorflowMultiTaskRegressor(

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

        learning_rate=0.0003, weight_init_stddevs=[.1], batch_size=64)

    model.reload()

    self.pose_scorer = GridPoseScorer(model, feat="grid")

    self.pose_generator = VinaPoseGenerator() 

  def dock(self, protein_file, ligand_file):

    """Docks using Vina and DNNs."""

    protein_docked, ligand_docked = self.pose_generator.generate_poses(

        protein_file, ligand_file)

    score = self.pose_scorer.score(protein_docked, ligand_docked)

    return (score, (protein_docked, ligand_docked))