Commit af94caf8 authored by sruizcarmona's avatar sruizcarmona

new sdtether and sdrmsd dani versions and prevented rbrmsd from compilation

parent 01e4d49a
#! /usr/bin/env python
#
# Calculate SMART RMSD with or without molecular superposition (FIT or NOFIT)
# Script distributed under GNU LGPL 3.0 along rDock software.
# Script distributed under GNU LGPL 3.0 along RDOCK software.
#
# This algorithm takes into account molecular automorphism. That is, it identifies
# molecules which are the same but might have atom orders changed and still be able to
# match the pairs and correctly calculate the RMSD.
#
# Author: Daniel Alvarez-Garcia
# Date: 08-11-2013
import math
import pybel
import numpy as npy
import optparse
def superpose3D(ref, target, weights=None,refmask=None,targetmask=None,returnRotMat=False):
"""superpose3D performs 3d superposition using a weighted Kabsch algorithm : http://dx.doi.org/10.1107%2FS0567739476001873 & doi: 10.1529/biophysj.105.066654
......@@ -97,26 +102,44 @@ def mapToCrystal(xtal, pose):
exit=mapper.MapUnique(pose.OBMol,mappingpose)
return mappingpose[0]
def parseArguments():
optparse.OptionParser.format_epilog = lambda self, formatter: self.epilog
epilog = """Args:
reference.sdf SDF file with the reference molecule.
input.sdf SDF file with the molecules to be compared to reference.\n"""
parser = optparse.OptionParser("usage: %prog [options] reference.sdf input.sdf", epilog=epilog)
parser.add_option("-f", "--fit",dest="fit", action="store_true", default=False,
help="Superpose molecules before RMSD calculation")
parser.add_option("-o","--out", dest="outfilename", metavar="FILE", default=False,
help="If declared, write an output SDF file with the input molecules with a new sdfield <RMSD>. If molecule was fitted, the fitted molecule coordinates will be saved.")
(options, args) = parser.parse_args()
#Check we have two arguments
if len(args) < 2:
parser.error("Incorrect number of arguments. Use -h or --help options to print help.")
return options, args
def updateCoords(obmol, newcoords):
"Update OBMol coordinates. newcoords is a numpy array"
for i,atom in enumerate(obmol):
atom.OBAtom.SetVector(*newcoords[i])
if __name__ == "__main__":
import sys
if len(sys.argv) < 3:
sys.exit("USAGE: smartRMSD.py reference.sd poses.sd [fit]")
(opts, args) = parseArguments()
xtal = sys.argv[1]
poses = sys.argv[2]
fit = False
try:
fit=sys.argv[3]
if fit == 'fit': fit=True
except:
fit = False
xtal = args[0]
poses = args[1]
fit = opts.fit
outfname = opts.outfilename
# Read crystal pose
crystal = next(pybel.readfile("sdf", xtal))
crystal.removeh()
crystalnumatoms = len(crystal.atoms)
# Find automorphisms involving only non-H atoms
mappings = pybel.ob.vvpairUIntUInt()
......@@ -126,13 +149,22 @@ if __name__ == "__main__":
lookup.append(i)
success = pybel.ob.FindAutomorphisms(crystal.OBMol, mappings)
#If outfname is defined, prepare an output SDF sink to write molecules
if outfname:
outsdf = pybel.Outputfile('sdf', outfname, overwrite=True)
# Find the RMSD between the crystal pose and each docked pose
xtalcoords = [atom.coords for atom in crystal]
dockedposes = pybel.readfile("sdf", poses)
if fit: print "POSE\tRMSD_FIT"
else: print "POSE\tRMSD_NOFIT"
skipped = []
for i, dockedpose in enumerate(dockedposes):
dockedpose.removeh()
natoms = len(dockedpose.atoms)
if natoms != crystalnumatoms:
skipped.append(i)
continue
mappose = mapToCrystal(crystal, dockedpose)
mappose = npy.array(mappose)
mappose = mappose[npy.argsort(mappose[:,0])][:,1]
......@@ -145,9 +177,28 @@ if __name__ == "__main__":
mapping_rmsd = rmsd(posecoords, automorph_coords)
if mapping_rmsd < resultrmsd:
resultrmsd = mapping_rmsd
fitted_result = False
if fit:
fitted_pose, fitted_rmsd = superpose3D(npy.array(automorph_coords), npy.array(posecoords))
if fitted_rmsd < resultrmsd:
resultrmsd = fitted_rmsd
fitted_result = fitted_pose
# Write molecule to ousdf if demanded
if outfname:
# Add RMSD field and save molecule
newData = pybel.ob.OBPairData()
newData.SetAttribute("RMSD")
newData.SetValue('%.3f'%resultrmsd)
dockedpose.OBMol.CloneData(newData) # Add new data
# Update position if fit option
if npy.any(fitted_result): updateCoords(dockedpose, fitted_result)
outsdf.write(dockedpose)
print "%d\t%.2f"%((i+1),resultrmsd)
if skipped: print >> sys.stderr, "SKIPPED input molecules due to number of atom missmatch: %d"%skipped
#! /usr/bin/env python
#
# Substitute for rbtether of rDock. Will align input molecules to a reference fragment defined by a smarts string,
# it will add a TETHERED ATOM property field to the output SDF that is correctly understood by rDock
# rDock will restrain the matching atom positions to the reference molecule coordinates.
# Substitute for rbtether of RDock. Will align input molecules to a reference fragment defined by a smarts string,
# it will add a TETHERED ATOM property field to the output SDF that is correctly understood by RDOCK
# RDock will restrain the matching atom positions to the reference molecule coordinates.
#
# Initially implemented with a conformational search algorithm to better match target coordinates.
# But had problems with OBabel FF generating non-sense conformers. So in this version the conformer search is commented out.
......@@ -10,7 +10,7 @@
# dimished or even vanish if the SMARTS string is defined for a rigid region (like a ring).
# I'm still trying to incorporate somehow this conformational search.
#
# Script distributed under GNU LGPL 3.0 along rDock software.
# Script distributed under GNU LGPL 3.0 along RDOCK software.
#
# Author: Daniel Alvarez-Garcia
# Date: 08-11-2013
......@@ -130,7 +130,7 @@ if __name__ == "__main__":
import sys
if len(sys.argv) != 5:
sys.exit("USAGE: tetherSDF.py reference.sdf input.sdf output.sdf 'SMARTS'")
sys.exit("USAGE: %s reference.sdf input.sdf output.sdf 'SMARTS'"%sys.argv[0])
refsdf = sys.argv[1]
molsdf = sys.argv[2]
......
......@@ -27,14 +27,23 @@ UNITTEST_LINKLIBS = -lRbt -ldl -lcppunit
# XB Build targets except the ones requiring daylight toolkits
all: exe unit_test
#XB 2013 version removed rbrms, new sdrmsd script makes it all
# Standalone rDock executables
exe: $(DESTDIR)/rbdock \
$(DESTDIR)/rbcavity \
$(DESTDIR)/rbrms \
$(DESTDIR)/rbmoegrid \
$(DESTDIR)/rblist \
$(DESTDIR)/rbcalcgrid
#old version
## Standalone rDock executables
#exe: $(DESTDIR)/rbdock \
# $(DESTDIR)/rbcavity \
# $(DESTDIR)/rbrms \
# $(DESTDIR)/rbmoegrid \
# $(DESTDIR)/rblist \
# $(DESTDIR)/rbcalcgrid
# These exes require the Daylight SMARTS/SMILES toolkits
dt_exe: $(DESTDIR)/rbtether \
$(DESTDIR)/smart_rms
......@@ -60,8 +69,9 @@ $(DESTDIR)/rbdock: $(SRCDIR)/rbdock.cxx $(DEPLIBS)
$(DESTDIR)/rbcalcgrid: $(SRCDIR)/rbcalcgrid.cxx $(DEPLIBS)
$(CXX) $(CXXFLAGS) $(SRCDIR)/rbcalcgrid.cxx $(LINKLIBS) -o [email protected]
$(DESTDIR)/rbrms: $(SRCDIR)/rbrms.cxx $(DEPLIBS)
$(CXX) $(CXXFLAGS) $(SRCDIR)/rbrms.cxx $(LINKLIBS) -o [email protected]
# XB comment following lines as rbrms is deprecated. sdrmsd makes it all
#$(DESTDIR)/rbrms: $(SRCDIR)/rbrms.cxx $(DEPLIBS)
# $(CXX) $(CXXFLAGS) $(SRCDIR)/rbrms.cxx $(LINKLIBS) -o [email protected]
# These exes depend on the Daylight SMARTS/SMILES toolkit libraries
# XB comment following lines as daylight is not needed in 2013
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment