Revision as of 22:54, 6 May 2013 by Oliv Eidam (talk | contribs) (Created page with " search a substructure in DOCK poses and calculate atom distances= Sometimes it is interesting to search for a substructure in your DOCK output poses and calculate...")
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search search a substructure in DOCK poses and calculate atom distances

Sometimes it is interesting to search for a substructure in your DOCK output poses and calculate distances to a reference point. This is what (courtesy from Nir) can do for you. For example, lots of opioid ligands exhibit a piperidine ring. The charged nitrogen of the piperidine ring forms a salt bridge with the conserved Asp[3.32] in the delta, mu and kappa opioid receptor crystal structures, and you may be interested to interrogate how many of your docked opioid ligands place the piperidine nitrogen within salt bridge distance of the conserved Aspartate.


Edit (see code below):
1. Define a substructure for your ligand(s) by providing smiles or smarts for your substructure (variable ss);
2. Define an atom within that substructure to calculate distance from a reference point (variable attack_site); the first atom in your smiles or smarts is 0, the second is 1, the third is 2, etc.

Then run providing coordinates and coordinates of the reference point: topdock.pdb -36.457 9.996 6.888 > filter.log

This outputs a line for each ligand that is within a the defined vicinity threshold t, which is 15.5 Angstrom in the example below. The last number in the output (filter.log) is the distance (4.59A below).

Example output in filter.log:

1       cmpd24  -55.98  -45.59  -29.67  24.45   -5.17 4.5911097889 code

#!/usr/bin/env python
import os
import sys
import math

from openeye.oechem import *
from openeye.oeomega import *

#usage ./ topdock.pdb -36.457   9.996   6.888

def calc_dist(p1,p2):
    return math.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2+(p1[2]-p2[2])**2)

#main program
#input file is the first argument
infile = sys.argv[1]
#x y z coords of target cys are the next 3 args 
cov_xyz = [0,0,0]
cov_xyz[0] = float(sys.argv[2])
cov_xyz[1] = float(sys.argv[3])
cov_xyz[2] = float(sys.argv[4])

ifs = oemolistream(infile)
mol = OEGraphMol()
lig_xyz = OEDoubleArray(3)

#define the vicinity threashold
t = 15.5
#define the reactive warhead and identify the attack site

#piperidine smiles; distance to nitrogen measured
ss = OESubSearch("CN1CCCCC1")
attack_site = 1

#amide smarts, distance to oxygen measured
#ss = OESubSearch("[NX3H1][CX3](=[OX1])[#6]")
#attack_site = 2

while OEReadMolecule(ifs, mol):
    for count,match in enumerate(ss.Match(mol,"true")):
        for ma in match.GetAtoms():
            #starting count from zero and target carbon is 5
            if (ma.pattern.GetIdx()==attack_site): 
                #print mol.GetTitle(),calc_dist(lig_xyz,cov_xyz)
                if d<t:
                    print mol.GetTitle(),d