Commit 7b1f6500 authored by marta-sd's avatar marta-sd
Browse files

new implementation of pi interactions

parent 637f10b1

deepchem/feat/rdkit_grid_featurizer.py

+196 −76
Original line number Diff line number Diff line
@@ -12,6 +12,7 @@ from warnings import warn 
import time

import tempfile

import hashlib

from collections import Counter

from rdkit import Chem

from rdkit.Chem import AllChem

from deepchem.utils.rdkit_util import load_molecule
@@ -162,6 +163,9 @@ def angle_between(vector_i, vector_j): 
  0.000000

  >>> print("%0.06f" % angle_between((1, 0, 0), (-1, 0, 0)))

  3.141593



  Note that this function always returns the smaller of the two angles between

  the vectors (value between 0 and pi).

  """

  vector_i_u = unit_vector(vector_i)

  vector_j_u = unit_vector(vector_j)
@@ -373,71 +377,60 @@ def featurize_splif(protein_xyz, protein, ligand_xyz, ligand, contact_bins, 
  return (splif_dicts)





def compute_ring_center(mol, ring):

  ring_xyz = np.zeros((len(ring._path), 3))

  for i, atom_idx in enumerate(ring._path):

    atom = mol.GetAtom(int(atom_idx))

    ring_xyz[i, :] = [atom.x(), atom.y(), atom.z()]

def compute_ring_center(mol, ring_indices):

  conformer = mol.GetConformer()

  ring_xyz = np.zeros((len(ring_indices), 3))

  for i, atom_idx in enumerate(ring_indices):

    atom_position = conformer.GetAtomPosition(atom_idx)

    ring_xyz[i] = np.array(atom_position)

  ring_centroid = compute_centroid(ring_xyz)

  return ring_centroid





def compute_ring_normal(mol, ring):

def compute_ring_normal(mol, ring_indices):

  conformer = mol.GetConformer()

  points = np.zeros((3, 3))

  for i, atom_idx in enumerate(ring._path):

    if i == 3: break

    atom = mol.GetAtom(int(atom_idx))

    points[i, :] = [atom.x(), atom.y(), atom.z()]

  for i, atom_idx in enumerate(ring_indices[:3]):

    atom_position = conformer.GetAtomPosition(atom_idx)

    points[i] = np.array(atom_position)



  v1 = points[1, :] - points[0, :]

  v2 = points[2, :] - points[0, :]

  v1 = points[1] - points[0]

  v2 = points[2] - points[0]

  normal = np.cross(v1, v2)

  return normal





def is_pi_parallel(protein_ring_center, protein_ring_normal, ligand_ring_center,

                   ligand_ring_normal):

  dist = np.linalg.norm(protein_ring_center - ligand_ring_center)

  angle = angle_between(protein_ring_normal, ligand_ring_normal) * 180 / np.pi

  if ((np.abs(angle) < 30.0) or

      (np.abs(angle) > 150.0 and np.abs(angle) < 210.0) or

      (np.abs(angle) > 330.0 and np.abs(angle) < 360.0)):

    if dist < 8.0:

def is_pi_parallel(ring1_center,

                   ring1_normal,

                   ring2_center,

                   ring2_normal,

                   dist_cutoff=8.0,

                   angle_cutoff=30.0):

  dist = np.linalg.norm(ring1_center - ring2_center)

  angle = angle_between(ring1_normal, ring2_normal) * 180 / np.pi

  if ((angle < angle_cutoff or angle > 180.0 - angle_cutoff) and

      dist < dist_cutoff):

    return True



  else:

    return False





def is_pi_t(protein_ring_center, protein_ring_normal, ligand_ring_center,

            ligand_ring_normal):

  dist = np.linalg.norm(protein_ring_center - ligand_ring_center)

  angle = angle_between(protein_ring_normal, ligand_ring_normal) * 180 / np.pi

  if ((np.abs(angle) > 60.0 and np.abs(angle) < 120.0) or

      (np.abs(angle) > 240.0 and np.abs(angle) < 300.0)):

    if dist < 5.5:

def is_pi_t(ring1_center,

            ring1_normal,

            ring2_center,

            ring2_normal,

            dist_cutoff=5.5,

            angle_cutoff=30.0):

  dist = np.linalg.norm(ring1_center - ring2_center)

  angle = angle_between(ring1_normal, ring2_normal) * 180 / np.pi

  if ((90.0 - angle_cutoff < angle < 90.0 + angle_cutoff) and

      dist < dist_cutoff):

    return True



  return False





def update_feature_dict(feature_dict, idxs=None, indices=None):

  if idxs is not None:

    indices = []

    for idx in idxs:

      indices.append(idx - 1)



  for index in indices:

    if index not in feature_dict.keys():

      feature_dict[index] = 1

  else:

      feature_dict[index] += 1



  return feature_dict

    return False





def compute_pi_stack(protein_xyz,

                     protein,

                     ligand_xyz,

def compute_pi_stack(protein,

                     ligand,

                     pairwise_distances=None,

                     dist_cutoff=4.4,
@@ -458,35 +451,118 @@ def compute_pi_stack(protein_xyz, 
            for each atom in ligand and in protein:

              add to list of atom indices

  """

  raise NotImplementedError("This function is not implemented yet")





def is_cation_pi(cation_position, ring_center, ring_normal):

  protein_pi_parallel = Counter()

  protein_pi_t = Counter()

  ligand_pi_parallel = Counter()

  ligand_pi_t = Counter()



  protein_aromatic_rings = []

  ligand_aromatic_rings = []

  for mol, ring_list in ((protein, protein_aromatic_rings),

                         (ligand, ligand_aromatic_rings)):

    aromatic_atoms = set(atom.GetIdx() for atom in mol.GetAromaticAtoms())

    for ring in Chem.GetSymmSSSR(mol):

      if set(ring).issubset(aromatic_atoms):

        ring_center = compute_ring_center(mol, ring)

        ring_normal = compute_ring_normal(mol, ring)

        ring_list.append((ring, ring_center, ring_normal))



  counted_pairs_parallel = set()

  counted_pairs_t = set()

  for prot_ring, prot_ring_center, prot_ring_normal in protein_aromatic_rings:

    for lig_ring, lig_ring_center, lig_ring_normal in ligand_aromatic_rings:

      if is_pi_parallel(

          prot_ring_center,

          prot_ring_normal,

          lig_ring_center,

          lig_ring_normal,

          angle_cutoff=angle_cutoff,

          dist_cutoff=dist_cutoff):

        prot_to_update = set()

        lig_to_update = set()

        for i in prot_ring:

          for j in lig_ring:

            if (i, j) not in counted_pairs_parallel:

              prot_to_update.add(i)

              lig_to_update.add(j)

              counted_pairs_parallel.add((i, j))



        protein_pi_parallel.update(prot_to_update)

        ligand_pi_parallel.update(lig_to_update)



      if is_pi_t(

          prot_ring_center,

          prot_ring_normal,

          lig_ring_center,

          lig_ring_normal,

          angle_cutoff=angle_cutoff,

          dist_cutoff=dist_cutoff):

        prot_to_update = set()

        lig_to_update = set()

        for i in prot_ring:

          for j in lig_ring:

            if (i, j) not in counted_pairs_t:

              prot_to_update.add(i)

              lig_to_update.add(j)

              counted_pairs_t.add((i, j))



        protein_pi_t.update(prot_to_update)

        ligand_pi_t.update(lig_to_update)



  return (protein_pi_t, protein_pi_parallel, ligand_pi_t, ligand_pi_parallel)





def is_cation_pi(cation_position,

                 ring_center,

                 ring_normal,

                 dist_cutoff=6.5,

                 angle_cutoff=30.0):

  cation_to_ring_vec = cation_position - ring_center

  dist = np.linalg.norm(cation_to_ring_vec)

  angle = angle_between(cation_to_ring_vec, ring_normal) * 180. / np.pi

  if dist < 6.5:

    if ((np.abs(angle) < 30.0) or

        (np.abs(angle) > 150.0 and np.abs(angle) < 210.0) or

        (np.abs(angle) > 330.0 and np.abs(angle) < 360.0)):

  if ((angle < angle_cutoff or angle > 180.0 - angle_cutoff) and

      (dist < dist_cutoff)):

    return True



  else:

    return False





def compute_cation_pi(protein, ligand, protein_cation_pi, ligand_cation_pi):

  raise NotImplementedError("This function is not implemented yet")

def compute_cation_pi(mol1, mol2, charge_tolerance=0.01):

  """Finds aromatic rings in mo1 interacting with cations in mol2"""

  mol1_pi = Counter()

  mol2_cation = Counter()

  conformer = mol2.GetConformer()



  aromatic_atoms = set(atom.GetIdx() for atom in mol1.GetAromaticAtoms())

  rings = [list(r) for r in Chem.GetSymmSSSR(mol1)]



  for ring in rings:

    if set(ring).issubset(aromatic_atoms):

      ring_center = compute_ring_center(mol1, ring)

      ring_normal = compute_ring_normal(mol1, ring)



def compute_binding_pocket_cation_pi(protein_xyz, protein, ligand_xyz, ligand):

  protein_cation_pi = {}

  ligand_cation_pi = {}

      for atom in mol2.GetAtoms():

        if atom.GetFormalCharge() > 1.0 - charge_tolerance:

          cation_position = np.array(conformer.GetAtomPosition(atom.GetIdx()))

          if is_cation_pi(cation_position, ring_center, ring_normal):

            mol1_pi.update(ring)

            mol2_cation.update([atom.GetIndex()])

  return mol1_pi, mol2_cation



  (protein_cation_pi, ligand_cation_pi) = compute_cation_pi(

      protein, ligand, protein_cation_pi, ligand_cation_pi)

  (ligand_cation_pi, protein_cation_pi) = compute_cation_pi(

      ligand, protein, ligand_cation_pi, protein_cation_pi)

  return (protein_cation_pi, ligand_cation_pi)



def compute_binding_pocket_cation_pi(protein, ligand):

  protein_pi, ligand_cation = compute_cation_pi(protein, ligand)

  ligand_pi, protein_cation = compute_cation_pi(ligand, protein)



  protein_cation_pi = Counter()

  protein_cation_pi.update(protein_pi)

  protein_cation_pi.update(protein_cation)



  ligand_cation_pi = Counter()

  ligand_cation_pi.update(ligand_pi)

  ligand_cation_pi.update(ligand_cation)



  return protein_cation_pi, ligand_cation_pi





def get_partial_charge(atom):
@@ -500,9 +576,8 @@ def get_partial_charge(atom): 




def get_formal_charge(atom):

  warn(

      'get_formal_charge function is deprecated, use get_partial_charge instead',

      DeprecationWarning)

  warn('get_formal_charge function is deprecated and will be removed'

       ' in version 1.4, use get_partial_charge instead', DeprecationWarning)

  return get_partial_charge(atom)
@@ -667,7 +742,8 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 
    ]

    for arg in deprecated_args:

      if arg in kwargs:

        warn('%s argument was removed and it is ignored' % arg,

        warn('%s argument was removed and it is ignored,'

             ' using it will result in error in version 1.4' % arg,

             DeprecationWarning)



    self.verbose = verbose
@@ -748,6 +824,36 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 
                self.hbond_angle_cutoffs)]

    }



    def voxelize_pi_stack(prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances):

      protein_pi_t, protein_pi_parallel, ligand_pi_t, ligand_pi_parallel = (

          compute_pi_stack(prot_rdk, lig_rdk, distances))

      pi_parallel_tensor = self._voxelize(

          convert_atom_to_voxel,

          None,

          prot_xyz,

          feature_dict=protein_pi_parallel,

          nb_channel=1)

      pi_parallel_tensor += self._voxelize(

          convert_atom_to_voxel,

          None,

          lig_xyz,

          feature_dict=ligand_pi_parallel,

          nb_channel=1)



      pi_t_tensor = self._voxelize(

          convert_atom_to_voxel,

          None,

          prot_xyz,

          feature_dict=protein_pi_t,

          nb_channel=1)

      pi_t_tensor += self._voxelize(

          convert_atom_to_voxel,

          None,

          lig_xyz,

          feature_dict=ligand_pi_t,

          nb_channel=1)

      return [pi_parallel_tensor, pi_t_tensor]



    self.VOXEL_FEATURES = {

        'ecfp': lambda prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances:

            [sum([self._voxelize(
@@ -839,7 +945,21 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 
                distances,

                self.hbond_dist_bins,

                self.hbond_angle_cutoffs)

            ]

            ],

        'pi_stack': voxelize_pi_stack,



        'cation_pi': lambda prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances:

            [sum([self._voxelize(

                convert_atom_to_voxel,

                None,

                xyz,

                feature_dict=cation_pi_dict,

                nb_channel=1

            ) for xyz, cation_pi_dict in zip(

                (prot_xyz, lig_xyz), compute_binding_pocket_cation_pi(

                    prot_rdk,

                    lig_rdk,

                ))])],

    }



    if feature_types is None:

deepchem/feat/tests/3ws9_ligand.sdf

0 → 100644
+117 −0
Original line number Diff line number Diff line 
3ws9_ligand



Created by X-TOOL on Sat Nov 28 16:04:28 2015

 46 50  0  0  0  0  0  0  0  0999 V2000

   17.9520    7.9520   55.8870  C 0  0  0  2  0  3

   19.0490    8.5940   56.4540  C 0  0  0  2  0  3

   17.8800    6.5720   55.8140  C 0  0  0  2  0  3

   24.2190    5.8580   55.5190  C 0  0  0  2  0  3

   22.8760    6.1220   55.6870  C 0  0  0  2  0  3

   20.0970    7.8390   56.9560  C 0  0  0  2  0  3

   18.9340    5.8230   56.3200  C 0  0  0  2  0  3

   24.4600    4.6370   57.5640  C 0  0  0  2  0  3

   17.0690    1.9920   61.1550  C 0  0  0  2  0  3

   24.9970    5.1470   56.4080  C 0  0  0  1  0  3

   23.1230    4.8880   57.7500  C 0  0  0  1  0  3

   22.3640    5.6030   56.8470  C 0  0  0  1  0  3

   20.0370    6.4530   56.8870  C 0  0  0  1  0  3

   17.5790    2.4330   62.3580  C 0  0  0  1  0  3

   18.9920    2.8850   60.7480  C 0  0  0  1  0  3

   21.1470    4.9560   58.5550  C 0  0  0  1  0  3

   17.1390    2.4280   63.6860  C 0  0  0  2  0  3

   17.9060    2.9730   64.6400  C 0  0  0  2  0  3

   19.1960    3.5540   64.1810  C 0  0  0  2  0  3

   26.4440    4.9360   56.1160  C 0  0  0  4  0  4

   20.2170    3.4070   60.1380  C 0  0  0  3  0  4

   20.0030    4.7400   59.4610  C 0  0  0  3  0  4

   17.9710    2.2830   60.1590  N 0  0  0  1  0  2

   22.3580    4.4990   58.8010  N 0  0  0  1  0  2

   19.6510    3.5860   62.9680  N 0  0  0  1  0  2

   21.0990    5.6500   57.3710  N 0  0  0  1  0  3

   18.8000    3.0020   62.0940  N 0  0  0  1  0  3

   17.1354    8.5465   55.4936  H 0  0  0  1  0  1

   19.0834    9.6764   56.5026  H 0  0  0  1  0  1

   17.0185    6.0857   55.3709  H 0  0  0  1  0  1

   24.6977    6.2368   54.6232  H 0  0  0  1  0  1

   22.2847    6.6838   54.9729  H 0  0  0  1  0  1

   20.9576    8.3266   57.3994  H 0  0  0  1  0  1

   18.8971    4.7406   56.2729  H 0  0  0  1  0  1

   25.0500    4.0757   58.2795  H 0  0  0  1  0  1

   16.1099    1.4951   61.0167  H 0  0  0  1  0  1

   16.1785    1.9837   63.9429  H 0  0  0  1  0  1

   17.6038    2.9895   65.6861  H 0  0  0  1  0  1

   19.8268    4.0012   64.9478  H 0  0  0  1  0  1

   26.6971    5.4095   55.1559  H 0  0  0  1  0  1

   26.6532    3.8576   56.0593  H 0  0  0  1  0  1

   27.0491    5.3853   56.9172  H 0  0  0  1  0  1

   20.5740    2.6844   59.3894  H 0  0  0  1  0  1

   20.9781    3.5281   60.9229  H 0  0  0  1  0  1

   19.0648    4.7251   58.8869  H 0  0  0  1  0  1

   19.9604    5.5421   60.2125  H 0  0  0  1  0  1

  2  1  4  0  0  1

  3  1  4  0  0  1

  6  2  4  0  0  1

  7  3  4  0  0  1

  5  4  4  0  0  1

  4 10  4  0  0  1

 12  5  4  0  0  1

 13  6  4  0  0  1

 13  7  4  0  0  1

  8 10  4  0  0  1

 11  8  4  0  0  1

 14  9  2  0  0  1

 23  9  1  0  0  1

 10 20  1  0  0  2

 12 11  4  0  0  1

 24 11  1  0  0  1

 26 12  1  0  0  1

 26 13  1  0  0  2

 14 17  1  0  0  1

 27 14  1  0  0  1

 21 15  1  0  0  2

 15 23  2  0  0  1

 15 27  1  0  0  1

 22 16  1  0  0  2

 16 24  2  0  0  1

 16 26  1  0  0  1

 17 18  2  0  0  1

 18 19  1  0  0  1

 19 25  2  0  0  1

 21 22  1  0  0  2

 25 27  1  0  0  1

  1 28  1  0  0  2

  2 29  1  0  0  2

  3 30  1  0  0  2

  4 31  1  0  0  2

  5 32  1  0  0  2

  6 33  1  0  0  2

  7 34  1  0  0  2

  8 35  1  0  0  2

  9 36  1  0  0  2

 17 37  1  0  0  2

 18 38  1  0  0  2

 19 39  1  0  0  2

 20 40  1  0  0  2

 20 41  1  0  0  2

 20 42  1  0  0  2

 21 43  1  0  0  2

 21 44  1  0  0  2

 22 45  1  0  0  2

 22 46  1  0  0  2

M  END

> <MOLECULAR_FORMULA>

C22H19N5



> <MOLECULAR_WEIGHT>

353.3



> <NUM_HB_ATOMS>

5  



> <NUM_ROTOR>

3  



> <XLOGP2>

4.58 



$$$$

deepchem/feat/tests/3ws9_protein_fixer_rdkit.pdb

0 → 100644
+4456 −0

File added.

Preview size limit exceeded, changes collapsed.

deepchem/feat/tests/test_rdkit_grid_features.py

+208 −69

File changed.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE