Filter.py:

Filter.py: search a substructure in DOCK poses and calculate atom distances

Opioid ligands exhibit a piperdine substructure and the charged nitrogen forms a salt bridge with Asp[3.32].

Sometimes it is interesting to search for a substructure in your DOCK output poses and calculate distances to a reference point. This is what filter.py (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 nociceptin 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 filter.py

Edit filter.py (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.
3. Find the x,y,z coordinates of the reference point, e.g.: -36.457 9.996 6.888

Then run filter.py providing pdb coordinates and coordinates of the reference point:

filter.py 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

filter.py code

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

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

#usage ./filter.py 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): 
                mol.GetCoords(ma.target,lig_xyz)
                #print mol.GetTitle(),calc_dist(lig_xyz,cov_xyz)
                d=calc_dist(lig_xyz,cov_xyz)
                if d<t:
                    print mol.GetTitle(),d
ifs.close()

Acknowledgments

Thanks to Nir for sharing this script and to OpenEye of course.