Filter.py: 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 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 (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
#!/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-p2)**2+(p1-p2)**2+(p1-p2)**2) #main program #input file is the first argument infile = sys.argv #x y z coords of target cys are the next 3 args cov_xyz = [0,0,0] cov_xyz = float(sys.argv) cov_xyz = float(sys.argv) cov_xyz = float(sys.argv) 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()
Thanks to Nir for sharing this script and to OpenEye of course.