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

Adding in Vina pose generation

parent 62802c9c

deepchem/dock/pose_generation.py

+91 −3
Original line number Diff line number Diff line
@@ -9,9 +9,45 @@ __author__ = "Bharath Ramsundar" 
__copyright__ = "Copyright 2016, Stanford University"

__license__ = "GPL"



import numpy as np

import os

import pybel

import tempfile

from deepchem.feat import hydrogenate_and_compute_partial_charges

from subprocess import call



def write_conf(receptor_filename, ligand_filename, centroid, box_dims,

               conf_filename, exhaustiveness=None):

  """Writes Vina configuration file to disk."""

  with open(conf_filename, "wb") as f:

    f.write("receptor = %s\n" % receptor_filename)

    f.write("ligand = %s\n\n" % ligand_filename)



    f.write("center_x = %f\n" % centroid[0])

    f.write("center_y = %f\n" % centroid[1])

    f.write("center_z = %f\n\n" % centroid[2])



    f.write("size_x = %f\n" % box_dims[0])

    f.write("size_y = %f\n" % box_dims[1])

    f.write("size_z = %f\n\n" % box_dims[2])



    if exhaustiveness is not None:

      f.write("exhaustiveness = %d\n" % exhaustiveness)



def get_molecule_data(pybel_molecule):

  """Uses pybel to compute centroid and range of molecule (Angstroms)."""

  atom_positions = []

  for atom in pybel_molecule:

    atom_positions.append(atom.coords)

  num_atoms = len(atom_positions)

  protein_xyz = np.asarray(atom_positions)

  protein_centroid = np.mean(protein_xyz, axis=0)

  protein_max = np.max(protein_xyz, axis=0)

  protein_min = np.min(protein_xyz, axis=0)

  protein_range = protein_max - protein_min

  return protein_centroid, protein_range





class VinaPoseGenerator(object):



  def __init__(self):
@@ -20,6 +56,8 @@ class VinaPoseGenerator(object): 
    self.vina_dir = os.path.join(current_dir, "autodock_vina_1_1_2_linux_x86")

    if not os.path.exists(self.vina_dir):

      print("Vina not available. Downloading")

      # TODO(rbharath): May want to move this file to S3 so we can ensure it's

      # always available.

      wget_cmd = "wget http://vina.scripps.edu/download/autodock_vina_1_1_2_linux_x86.tgz"

      call(wget_cmd.split())

      print("Downloaded Vina. Extracting")
@@ -34,6 +72,56 @@ class VinaPoseGenerator(object): 
    self.vina_cmd = os.path.join(self.vina_dir, "bin/vina")

      



  def generate_poses(self, pdb_file, ligand_file):

    """Generates the docked complex and outputs the file."""

    pass
 
  def generate_poses(self, protein_file, ligand_file, out_dir=None):

    """Generates the docked complex and outputs files for docked complex."""

    if out_dir is None:

      out_dir = tempfile.mkdtemp()



    # Prepare receptor 

    receptor_name = os.path.basename(protein_file).split(".")[0]

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

    print("receptor_name")

    print(receptor_name)

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

    protein_hyd = os.path.join(out_dir, "%s.pdb" % receptor_name)

    protein_pdbqt = os.path.join(out_dir, "%s.pdbqt" % receptor_name)

    hydrogenate_and_compute_partial_charges(protein_file, "pdb",

                                            hyd_output=protein_hyd,

                                            pdbqt_output=protein_pdbqt,

                                            protein=True)

    # Get protein centroid and range

    receptor_pybel = next(pybel.readfile(str("pdb"), str(protein_hyd)))

    protein_centroid, protein_range = get_molecule_data(receptor_pybel)

    box_dims = protein_range + 5.0



    # Prepare receptor

    ligand_name = os.path.basename(ligand_file).split(".")[0]

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

    print("ligand_name")

    print(ligand_name)

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

    ligand_hyd = os.path.join(out_dir, "%s.pdb" % ligand_name)

    ligand_pdbqt = os.path.join(out_dir, "%s.pdbqt" % ligand_name)



    # TODO(rbharath): Generalize this so can support mol2 files as well.

    hydrogenate_and_compute_partial_charges(ligand_file, "sdf",

                                            hyd_output=ligand_hyd,

                                            pdbqt_output=ligand_pdbqt,

                                            protein=False)



    # Write Vina conf file

    conf_file = os.path.join(out_dir, "conf.txt")

    write_conf(protein_pdbqt, ligand_pdbqt, protein_centroid,

               box_dims, conf_file, exhaustiveness=1)



    # Define locations of log and output files

    log_file = os.path.join(out_dir, "%s_log.txt" % ligand_name)

    out_pdbqt = os.path.join(out_dir, "%s_docked.pdbqt" % ligand_name)

    # TODO(rbharath): Let user specify the number of poses required.

    print("About to call Vina")

    call("%s --config %s --log %s --out %s"

         % (self.vina_cmd, conf_file, log_file, out_pdbqt), shell=True)

    # TODO(rbharath): Convert the output pdbqt to a pdb file.



    # Return docked files 

    return protein_hyd, out_pdbqt

deepchem/dock/tests/test_pose_generation.py

+16 −1
Original line number Diff line number Diff line
@@ -21,7 +21,22 @@ class TestPoseGeneration(unittest.TestCase): 
  Does sanity checks on pose generation. 

  """



  def test_vina_initialization(self):

    """Test that VinaPoseGenerator can be initialized."""

    # Note this may download autodock Vina...

    vpg = dc.dock.VinaPoseGenerator()



  def test_vina_poses(self):

    """Test that VinaPoseGenerator Functions."""

    """Test that VinaPoseGenerator creates pose files."""

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

    protein_file = os.path.join(current_dir, "1jld_protein.pdb")

    ligand_file = os.path.join(current_dir, "1jld_ligand.sdf")

  

    # Note this may download autodock Vina...

    vpg = dc.dock.VinaPoseGenerator()

    protein_pose_file, ligand_pose_file = vpg.generate_poses(

        protein_file, ligand_file, out_dir="/tmp")



    # Check returned files exist

    assert os.path.exists(protein_pose_file)

    assert os.path.exists(ligand_pose_file)

deepchem/feat/__init__.py

+1 −0
Original line number Diff line number Diff line
@@ -17,3 +17,4 @@ 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

from deepchem.feat.nnscore_utils import hydrogenate_and_compute_partial_charges

deepchem/feat/nnscore_utils.py

+38 −15
Original line number Diff line number Diff line
@@ -51,7 +51,7 @@ def hydrogenate_and_compute_partial_charges(input_file, input_format, 
                                            hyd_output=None,

                                            pdbqt_output=None,

                                            protein=True,

                                            verbose=False):

                                            verbose=True):

  """Outputs a hydrogenated pdb and a pdbqt with partial charges.



  Takes an input file in specified format. Generates two outputs:
@@ -77,40 +77,63 @@ def hydrogenate_and_compute_partial_charges(input_file, input_format, 
  if verbose:

    print("Create pdb with hydrogens added")

  hyd_conversion = openbabel.OBConversion()

  hyd_conversion.SetInAndOutFormats(str(input_format), str("pdb"))

  hyd_conv = hyd_conversion.SetInAndOutFormats(str(input_format), str("pdb"))

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

  print("input_format")

  print(input_format)

  print("hyd_conv")

  print(hyd_conv)

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

  mol = openbabel.OBMol()

  hyd_conversion.ReadFile(mol, str(input_file))

  # AddHydrogens(not-polaronly, correctForPH, pH)

  mol.AddHydrogens(False, True, 7.4)

  hyd_conversion.WriteFile(mol, str(hyd_output))

  hyd_out = hyd_conversion.WriteFile(mol, str(hyd_output))

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

  print("hyd_out")

  print(hyd_out)

  print("os.path.exists(hyd_output)")

  print(os.path.exists(hyd_output))

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



  if verbose:

    print("Create a pdbqt file from the hydrogenated pdb above.")

  charge_conversion = openbabel.OBConversion()

  charge_conversion.SetInAndOutFormats(str("pdb"), str("pdbqt"))

  charge_conv = charge_conversion.SetInAndOutFormats(str("pdb"), str("pdbqt"))

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

  print("charge_conv")

  print(charge_conv)

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



  if protein and verbose:

    print("Make protein rigid.")

  if protein:

    charge_conversion.AddOption(str("c"), charge_conversion.OUTOPTIONS)

    print("Make protein rigid.")

    charge_conversion.AddOption(str("r"), charge_conversion.OUTOPTIONS)

  if verbose:

    charge_conversion.AddOption(str("c"), charge_conversion.OUTOPTIONS)

  print("Preserve hydrogens")

  charge_conversion.AddOption(str("h"), charge_conversion.OUTOPTIONS)

  if verbose:

  print("Preserve atom indices")

  charge_conversion.AddOption(str("p"), charge_conversion.OUTOPTIONS)

  if verbose:

  print("preserve atom indices.")

  charge_conversion.AddOption(str("n"), charge_conversion.OUTOPTIONS)



  if verbose:

  print("About to run obabel conversion.")

  mol = openbabel.OBMol()

  charge_conversion.ReadFile(mol, str(hyd_output))

  force_partial_charge_computation(mol)

  charge_conversion.WriteFile(mol, str(pdbqt_output))



  if protein:

    print("Removing ROOT/ENDROOT/TORSDOF")

    with open(pdbqt_output) as f:

      pdbqt_lines = f.readlines()

    filtered_lines = []

    for line in pdbqt_lines:

      if "ROOT" in line or "ENDROOT" in line or "TORSDOF" in line:

        continue

      filtered_lines.append(line)

    with open(pdbqt_output, "w") as f:

      f.writelines(filtered_lines)



class AromaticRing(object):

  """Holds information about an aromatic ring."""

  def __init__(self, center, indices, plane_coeff, radius):

CRA Git | Maintained and supported by SUSTech CRA and CCSE