""" ODDT's internal docking/scoring engines """
import numpy as np
import math
from oddt.spatial import distance, dihedral, rotate
[docs]def get_children(molecule, mother, restricted):
atoms = np.zeros(len(molecule.atoms), dtype=bool)
atoms[mother] = True
d = 1 # Pass first
prev = 0
while d > 0:
atoms[[n.idx0
for i in np.nonzero(atoms)[0] if i != restricted
for n in molecule.atoms[i].neighbors if n.idx0 != restricted]] = True
d = atoms.sum() - prev
prev = atoms.sum()
return atoms
[docs]def get_close_neighbors(molecule, a_idx, num_bonds=1):
atoms = np.zeros(len(molecule.atoms), dtype=bool)
atoms[a_idx] = True
for i in range(num_bonds):
atoms[[n.idx0
for i in np.nonzero(atoms)[0]
for n in molecule.atoms[i].neighbors]] = True
return atoms
[docs]def change_dihedral(coords, a1, a2, a3, a4, target_angle, rot_mask):
sin = math.sin(target_angle)
cos = math.cos(target_angle)
t = 1 - cos
v0 = coords[a2] - coords[a3]
v = (v0) / np.linalg.norm(v0)
rot_matrix = np.array([[t * v[0] * v[0] + cos,
t * v[0] * v[1] + sin * v[2],
t * v[0] * v[2] - sin * v[1]],
[t * v[0] * v[1] - sin * v[2],
t * v[1] * v[1] + cos,
t * v[1] * v[2] + sin * v[0]],
[t * v[0] * v[2] + sin * v[1],
t * v[1] * v[2] - sin * v[0],
t * v[2] * v[2] + cos]])
centroid = coords[a3]
coords = coords - centroid
# Old and slower version
# coords[rot_mask] = (coords[rot_mask,np.newaxis,:] * rot_matrix).sum(axis=2)
coords[rot_mask] = np.einsum("ij,jk->ik", coords[rot_mask], rot_matrix)
return coords + centroid
[docs]def num_rotors_pdbqt(lig):
i = 0
for atom in lig.atoms:
if atom.atomicnum == 1:
continue
num_local_rot = sum(int(b.isrotor) for b in atom.bonds)
if num_local_rot == 0:
pass
elif num_local_rot == 1:
i += 0.5
elif num_local_rot == 2:
i += 1.0
elif num_local_rot >= 3:
i += 0.5
return i
[docs]class vina_docking(object):
def __init__(self, rec, lig=None, box=None, box_size=1., weights=None):
self.box_size = box_size # TODO: Unify box
if rec:
self.set_protein(rec)
if lig:
self.set_ligand(lig)
self.set_box(box)
# constants
self.weights = weights or np.array((-0.0356, -0.00516, 0.840, -0.0351, -0.587, 0.0585))
self.mask_inter = {}
self.mask_intra = {}
[docs] def set_box(self, box):
if box is not None:
self.box = np.array(box)
# delete unused atoms
r = self.rec_dict['coords']
# X, Y, Z within box and cutoff
mask = (self.box[0][0] - 8 <= r[:, 0]) & (r[:, 0] <= self.box[1][0] + 8)
mask *= (self.box[0][1] - 8 <= r[:, 1]) & (r[:, 1] <= self.box[1][1] + 8)
mask *= (self.box[0][2] - 8 <= r[:, 2]) & (r[:, 2] <= self.box[1][2] + 8)
self.rec_dict = self.rec_dict#[mask]
else:
self.box = box
[docs] def set_protein(self, rec):
if rec is None:
self.rec_dict = None
self.mask_inter = {}
else:
self.rec_dict = rec.atom_dict[rec.atom_dict['atomicnum'] != 1].copy()
self.rec_dict = self.correct_radius(self.rec_dict)
self.mask_inter = {}
[docs] def set_ligand(self, lig):
lig_hvy_mask = (lig.atom_dict['atomicnum'] != 1)
self.lig_dict = lig.atom_dict[lig_hvy_mask].copy()
self.lig_dict = self.correct_radius(self.lig_dict)
self.num_rotors = num_rotors_pdbqt(lig)
self.mask_inter = {}
self.mask_intra = {}
# Find distant members (min 3 consecutive bonds)
mask = np.vstack([~get_close_neighbors(lig, i, num_bonds=3)
for i in range(len(lig.atoms))])
mask = mask[lig_hvy_mask[np.newaxis, :] * lig_hvy_mask[:, np.newaxis]]
self.lig_distant_members = mask.reshape(lig_hvy_mask.sum(), lig_hvy_mask.sum())
# prepare rotors dictionary
self.rotors = []
for b in lig.bonds:
if b.isrotor:
a2 = int(b.atoms[0].idx0)
a3 = int(b.atoms[1].idx0)
for n in b.atoms[0].neighbors:
if a3 != int(n.idx0) and n.atomicnum != 1:
a1 = int(n.idx0)
break
for n in b.atoms[1].neighbors:
if a2 != int(n.idx0) and n.atomicnum != 1:
a4 = int(n.idx0)
break
rot_mask = get_children(lig, a3, a2)[lig_hvy_mask]
# translate atom indicies to lig_dict indicies (heavy only)
a1 = np.argwhere(self.lig_dict['id'] == a1).flatten()[0]
a2 = np.argwhere(self.lig_dict['id'] == a2).flatten()[0]
a3 = np.argwhere(self.lig_dict['id'] == a3).flatten()[0]
a4 = np.argwhere(self.lig_dict['id'] == a4).flatten()[0]
# rotate smaller part of the molecule
if rot_mask.sum() > len(rot_mask):
rot_mask = -rot_mask
a4, a3, a2, a1 = a1, a2, a3, a4
self.rotors.append({'atoms': (a1, a2, a3, a4), 'mask': rot_mask})
# Setup cached ligand coords
self.lig = vina_ligand(self.lig_dict['coords'].copy(), len(self.rotors), self, self.box_size)
[docs] def set_coords(self, coords):
self.lig_dict['coords'] = coords
[docs] def score(self, coords=None):
return (self.score_inter(coords) * self.weights[:5]).sum() / (1 + self.weights[5] * self.num_rotors)
# inter = (self.score_inter(coords) * self.weights[:5]).sum()
# total = (self.score_total(coords) * self.weights[:5]).sum()
# return total/(1+self.weights[5]*self.num_rotors)
[docs] def weighted_total(self, coords=None):
return (self.score_total(coords) * self.weights[:5]).sum()
[docs] def score_total(self, coords=None):
return self.score_inter(coords) + self.score_intra(coords)
[docs] def weighted_inter(self, coords=None):
return (self.score_inter(coords) * self.weights[:5]).sum()
[docs] def weighted_intra(self, coords=None):
return (self.score_intra(coords) * self.weights[:5]).sum()
[docs] def score_inter(self, coords=None):
if coords is None:
coords = self.lig_dict['coords']
# Inter-molecular
r = distance(self.rec_dict['coords'], coords)
d = (r - self.rec_dict['radius'][:, np.newaxis] - self.lig_dict['radius'][np.newaxis, :])
mask = r < 8
inter = []
# Gauss 1
inter.append(np.exp(-(d[mask] / 0.5)**2).sum())
# Gauss 2
inter.append(np.exp(-((d[mask] - 3.) / 2.)**2).sum())
# Repiulsion
inter.append((d[(d < 0) & mask]**2).sum())
# Hydrophobic
if 'hyd' not in self.mask_inter:
self.mask_inter['hyd'] = ((self.rec_dict['ishydrophobe'] | self.rec_dict['ishalogen'])[:, np.newaxis] *
(self.lig_dict['ishydrophobe'] | self.lig_dict['ishalogen'])[np.newaxis, :])
mask_hyd = mask & self.mask_inter['hyd']
d_hyd = d[mask_hyd]
inter.append((d_hyd <= 0.5).sum() + (1.5 - d_hyd[(0.5 < d_hyd) & (d_hyd < 1.5)]).sum())
# H-Bonding
if 'da' not in self.mask_inter:
self.mask_inter['da'] = ((self.rec_dict['isdonor'] | self.rec_dict['ismetal'])[:, np.newaxis] *
self.lig_dict['isacceptor'][np.newaxis, :])
if 'ad' not in self.mask_inter:
self.mask_inter['ad'] = (self.rec_dict['isacceptor'][:, np.newaxis] *
(self.lig_dict['isdonor'] | self.lig_dict['ismetal'])[np.newaxis, :])
d_h = d[mask & (self.mask_inter['da'] | self.mask_inter['ad'])]
inter.append((d_h <= -0.7).sum() + (d_h[(-0.7 < d_h) & (d_h < 0)] / -0.7).sum())
return np.array(inter)
[docs] def score_intra(self, coords=None):
if coords is None:
coords = self.lig_dict['coords']
# Intra-molceular
r = distance(coords, coords)
d = (r - self.lig_dict['radius'][:, np.newaxis] - self.lig_dict['radius'][np.newaxis, :])
mask = self.lig_distant_members & (r < 8)
intra = []
# Gauss 1
intra.append(np.exp(-(d[mask] / 0.5)**2).sum())
# Gauss 2
intra.append(np.exp(-((d[mask] - 3.) / 2.)**2).sum())
# Repiulsion
intra.append((d[(d < 0) & mask]**2).sum())
# Hydrophobic
if 'hyd' not in self.mask_intra:
self.mask_intra['hyd'] = ((self.lig_dict['ishydrophobe'] | self.lig_dict['ishalogen'])[:, np.newaxis] *
(self.lig_dict['ishydrophobe'] | self.lig_dict['ishalogen'])[np.newaxis, :])
mask_hyd = mask & self.mask_intra['hyd']
d_hyd = d[mask_hyd]
intra.append((d_hyd <= 0.5).sum() + (1.5 - d_hyd[(0.5 < d_hyd) & (d_hyd < 1.5)]).sum())
# H-Bonding
if 'da' not in self.mask_intra:
self.mask_intra['da'] = ((self.lig_dict['isdonor'] | self.lig_dict['ismetal'])[..., np.newaxis] *
self.lig_dict['isacceptor'][np.newaxis, ...])
if 'ad' not in self.mask_intra:
self.mask_intra['ad'] = (self.lig_dict['isacceptor'][..., np.newaxis] *
(self.lig_dict['isdonor'] | self.lig_dict['ismetal'])[np.newaxis, ...])
d_h = d[mask & (self.mask_intra['da'] | self.mask_intra['ad'])]
intra.append((d_h <= -0.7).sum() + (d_h[(-0.7 < d_h) & (d_h < 0)] / -0.7).sum())
return np.array(intra)
[docs] def correct_radius(self, atom_dict):
vina_r = {6: 1.9,
7: 1.8,
8: 1.7,
9: 1.5,
15: 2.1,
16: 2.0,
17: 1.8,
35: 2.0,
53: 2.2,
}
for a, r in vina_r.items():
atom_dict['radius'][atom_dict['atomicnum'] == a] = r
# metals - 1.2 A
atom_dict['radius'][atom_dict['ismetal']] = 1.2
return atom_dict
[docs]class vina_ligand(object):
def __init__(self, c0, num_rotors, engine, box_size=1):
self.c0 = c0.copy()
self.x0 = np.zeros(6 + num_rotors)
self.c1 = c0.copy()
self.x1 = np.zeros_like(self.x0)
self.engine = engine
self.box_size = box_size
[docs] def mutate(self, x2, force=False):
delta_x0 = x2 - self.x0
delta_x1 = x2 - self.x1
if not force and (delta_x1 != 0).sum() <= 3:
return self._inc_mutate(delta_x1, self.c1)
elif not force and (delta_x0 != 0).sum() <= 3:
return self._inc_mutate(delta_x0, self.c0)
else:
return self._full_mutate(x2)
def _full_mutate(self, x):
c = self.c0.copy()
trans_vec = x[:3]
rot_vec = x[3:6]
rotors_vec = x[6:]
c = rotate(c, *rot_vec) + trans_vec * self.box_size
for i in range(len(rotors_vec)):
a = self.engine.rotors[i]['atoms']
mask = self.engine.rotors[i]['mask']
c = change_dihedral(c, a[0], a[1], a[2], a[3], rotors_vec[i], mask)
self.c1 = c.copy()
self.x1 = x.copy()
return c
def _inc_mutate(self, x, c):
c = c.copy()
trans_vec = x[:3]
rot_vec = x[3:6]
rotors_vec = x[6:]
if (rot_vec != 0).any():
c = rotate(c, *rot_vec)
if (trans_vec != 0).any():
c += trans_vec * self.box_size
for i in np.where(rotors_vec != 0)[0]:
a = self.engine.rotors[i]['atoms']
mask = self.engine.rotors[i]['mask']
c = change_dihedral(c, a[0], a[1], a[2], a[3], rotors_vec[i], mask)
return c