
* add recipe for improved-rdock
* style: fix format
* style: fix format again
* Fix location of the directory
* Fix copyright
* Fix according to the reviewer's comments.
* Revert "Fix according to the reviewer's comments."
This reverts commit 4877877daf
.
* Fix according to the reviewer's comments.
* style: fix format
* Update var/spack/repos/builtin/packages/improved-rdock/package.py
Co-authored-by: Tamara Dahlgren <35777542+tldahlgren@users.noreply.github.com>
Co-authored-by: Yuichi Otsuka <otsukay@riken.jp>
Co-authored-by: Tamara Dahlgren <35777542+tldahlgren@users.noreply.github.com>
277 lines
12 KiB
Diff
277 lines
12 KiB
Diff
diff -u -r -N a/bin/sdrmsd b/bin/sdrmsd
|
|
--- a/bin/sdrmsd 2020-10-15 13:34:27.000000000 +0900
|
|
+++ b/bin/sdrmsd 2020-10-19 09:21:43.000000000 +0900
|
|
@@ -11,7 +11,7 @@
|
|
# Date: 08-11-2013
|
|
|
|
import math
|
|
-import pybel
|
|
+from openbabel import pybel
|
|
import numpy as npy
|
|
import optparse
|
|
|
|
@@ -103,24 +103,24 @@
|
|
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,
|
|
+ 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("--threshold","-t",dest="threshold", action="store", nargs=1,
|
|
+ parser.add_option("--threshold","-t",dest="threshold", action="store", nargs=1,
|
|
help="Discard poses with RMSD < THRESHOLD with respect previous poses which where not rejected based on same principle. A Population SDField will be added to output SD with the population number.", type=float)
|
|
- parser.add_option("-o","--out", dest="outfilename", metavar="FILE", default=False,
|
|
+ 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.")
|
|
+ (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
|
|
+ return options, args
|
|
|
|
def updateCoords(obmol, newcoords):
|
|
"Update OBMol coordinates. newcoords is a numpy array"
|
|
@@ -133,8 +133,8 @@
|
|
for correct RMSD comparison. Only the lowest RMSD will be returned.
|
|
|
|
Returns:
|
|
- If fit=False: bestRMSD (float) RMSD of the best matching mapping.
|
|
- If fit=True: (bestRMSD, molecCoordinates) (float, npy.array) RMSD of best match and its molecule fitted coordinates.
|
|
+ If fit=False: bestRMSD (float) RMSD of the best matching mapping.
|
|
+ If fit=True: (bestRMSD, molecCoordinates) (float, npy.array) RMSD of best match and its molecule fitted coordinates.
|
|
"""
|
|
mappings = pybel.ob.vvpairUIntUInt()
|
|
bitvec = pybel.ob.OBBitVec()
|
|
@@ -148,18 +148,18 @@
|
|
posecoords = npy.array([atom.coords for atom in molec])[mappose]
|
|
resultrmsd = 999999999999
|
|
for mapping in mappings:
|
|
- automorph_coords = [None] * len(targetcoords)
|
|
- for x, y in mapping:
|
|
- automorph_coords[lookup.index(x)] = targetcoords[lookup.index(y)]
|
|
- 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
|
|
+ automorph_coords = [None] * len(targetcoords)
|
|
+ for x, y in mapping:
|
|
+ automorph_coords[lookup.index(x)] = targetcoords[lookup.index(y)]
|
|
+ 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
|
|
|
|
if fit:
|
|
return (resultrmsd, fitted_pose)
|
|
@@ -167,16 +167,16 @@
|
|
return resultrmsd
|
|
|
|
def saveMolecWithRMSD(outsdf, molec, rmsd, population=False):
|
|
- newData = pybel.ob.OBPairData()
|
|
+ newData = pybel.ob.OBPairData()
|
|
newData.SetAttribute("RMSD")
|
|
newData.SetValue('%.3f'%rmsd)
|
|
|
|
if population:
|
|
- popData = pybel.ob.OBPairData()
|
|
- popData.SetAttribute("Population")
|
|
- popData.SetValue('%i'%population)
|
|
- molec.OBMol.CloneData(popData)
|
|
-
|
|
+ popData = pybel.ob.OBPairData()
|
|
+ popData.SetAttribute("Population")
|
|
+ popData.SetValue('%i'%population)
|
|
+ molec.OBMol.CloneData(popData)
|
|
+
|
|
molec.OBMol.CloneData(newData) # Add new data
|
|
outsdf.write(molec)
|
|
|
|
@@ -184,13 +184,13 @@
|
|
import sys, os
|
|
|
|
(opts, args) = parseArguments()
|
|
-
|
|
+
|
|
xtal = args[0]
|
|
poses = args[1]
|
|
|
|
if not os.path.exists(xtal) or not os.path.exists(poses):
|
|
- sys.exit("Input files not found. Please check the path given is correct.")
|
|
-
|
|
+ sys.exit("Input files not found. Please check the path given is correct.")
|
|
+
|
|
fit = opts.fit
|
|
outfname = opts.outfilename
|
|
threshold = opts.threshold
|
|
@@ -202,71 +202,71 @@
|
|
|
|
#If outfname is defined, prepare an output SDF sink to write molecules
|
|
if outfname:
|
|
- outsdf = pybel.Outputfile('sdf', outfname, overwrite=True)
|
|
+ outsdf = pybel.Outputfile('sdf', outfname, overwrite=True)
|
|
|
|
# Find the RMSD between the crystal pose and each docked pose
|
|
dockedposes = pybel.readfile("sdf", poses)
|
|
- if fit: print "POSE\tRMSD_FIT"
|
|
- else: print "POSE\tRMSD_NOFIT"
|
|
+ if fit: print ("POSE\tRMSD_FIT")
|
|
+ else: print ("POSE\tRMSD_NOFIT")
|
|
skipped = []
|
|
- moleclist = {} # Save all poses with their dockid
|
|
- population = {} # Poses to be written
|
|
+ moleclist = {} # Save all poses with their dockid
|
|
+ population = {} # Poses to be written
|
|
outlist = {}
|
|
for docki, dockedpose in enumerate(dockedposes):
|
|
dockedpose.removeh()
|
|
- natoms = len(dockedpose.atoms)
|
|
- if natoms != crystalnumatoms:
|
|
- skipped.append(docki+1)
|
|
- continue
|
|
- if fit:
|
|
- resultrmsd, fitted_result = getAutomorphRMSD(crystal, dockedpose, fit=True)
|
|
- updateCoords(dockedpose, fitted_result)
|
|
- else:
|
|
- resultrmsd = getAutomorphRMSD(crystal, dockedpose, fit=False)
|
|
-
|
|
- if threshold:
|
|
- # Calculate RMSD between all previous poses
|
|
- # Discard if rmsd < FILTER threshold
|
|
- if moleclist:
|
|
- match = None
|
|
- bestmatchrmsd = 999999
|
|
- for did,prevmol in moleclist.iteritems():
|
|
- tmprmsd = getAutomorphRMSD(prevmol, dockedpose)
|
|
- if tmprmsd < threshold:
|
|
- if tmprmsd < bestmatchrmsd:
|
|
- bestmatchrmsd = tmprmsd
|
|
- match = did
|
|
-
|
|
- if match != None:
|
|
- # Do not write this one
|
|
- # sum one up to the matching previous molecule id
|
|
- print >> sys.stderr, "Pose %i matches pose %i with %.3f RMSD"%(docki+1, match+1, bestmatchrmsd)
|
|
- population[match] += 1
|
|
- else:
|
|
- # There's no match. Print info for this one and write to outsdf if needed
|
|
- # Save this one!
|
|
- if outfname: outlist[docki] = (dockedpose, resultrmsd)
|
|
- print "%d\t%.2f"%((docki+1),resultrmsd)
|
|
- moleclist[docki] = dockedpose
|
|
- population[docki] = 1
|
|
- else:
|
|
- # First molecule in list. Append for sure
|
|
- moleclist[docki] = dockedpose
|
|
- population[docki] = 1
|
|
- if outfname: outlist[docki] = (dockedpose, resultrmsd)
|
|
- else:
|
|
- # Just write best rmsd found and the molecule to outsdf if demanded
|
|
- if outfname: saveMolecWithRMSD(outsdf, dockedpose, resultrmsd)
|
|
- print "%d\t%.2f"%((docki+1),resultrmsd)
|
|
+ natoms = len(dockedpose.atoms)
|
|
+ if natoms != crystalnumatoms:
|
|
+ skipped.append(docki+1)
|
|
+ continue
|
|
+ if fit:
|
|
+ resultrmsd, fitted_result = getAutomorphRMSD(crystal, dockedpose, fit=True)
|
|
+ updateCoords(dockedpose, fitted_result)
|
|
+ else:
|
|
+ resultrmsd = getAutomorphRMSD(crystal, dockedpose, fit=False)
|
|
+
|
|
+ if threshold:
|
|
+ # Calculate RMSD between all previous poses
|
|
+ # Discard if rmsd < FILTER threshold
|
|
+ if moleclist:
|
|
+ match = None
|
|
+ bestmatchrmsd = 999999
|
|
+ for did,prevmol in moleclist.iteritems():
|
|
+ tmprmsd = getAutomorphRMSD(prevmol, dockedpose)
|
|
+ if tmprmsd < threshold:
|
|
+ if tmprmsd < bestmatchrmsd:
|
|
+ bestmatchrmsd = tmprmsd
|
|
+ match = did
|
|
+
|
|
+ if match != None:
|
|
+ # Do not write this one
|
|
+ # sum one up to the matching previous molecule id
|
|
+ print >> sys.stderr, "Pose %i matches pose %i with %.3f RMSD"%(docki+1, match+1, bestmatchrmsd)
|
|
+ population[match] += 1
|
|
+ else:
|
|
+ # There's no match. Print info for this one and write to outsdf if needed
|
|
+ # Save this one!
|
|
+ if outfname: outlist[docki] = (dockedpose, resultrmsd)
|
|
+ print ("%d\t%.2f"%((docki+1),resultrmsd))
|
|
+ moleclist[docki] = dockedpose
|
|
+ population[docki] = 1
|
|
+ else:
|
|
+ # First molecule in list. Append for sure
|
|
+ moleclist[docki] = dockedpose
|
|
+ population[docki] = 1
|
|
+ if outfname: outlist[docki] = (dockedpose, resultrmsd)
|
|
+ else:
|
|
+ # Just write best rmsd found and the molecule to outsdf if demanded
|
|
+ if outfname: saveMolecWithRMSD(outsdf, dockedpose, resultrmsd)
|
|
+ print ("%d\t%.2f"%((docki+1),resultrmsd))
|
|
|
|
if outlist:
|
|
- # Threshold applied and outlist need to be written
|
|
- for docki in sorted(outlist.iterkeys()):
|
|
- molrmsd = outlist[docki]
|
|
- # Get number of matchs in thresholding operation
|
|
- pop = population.get(docki)
|
|
- if not pop: pop = 1
|
|
- # Save molecule
|
|
- saveMolecWithRMSD(outsdf, molrmsd[0], molrmsd[1], pop)
|
|
-
|
|
- if skipped: print >> sys.stderr, "SKIPPED input molecules due to number of atom missmatch: %s"%skipped
|
|
+ # Threshold applied and outlist need to be written
|
|
+ for docki in sorted(outlist.iterkeys()):
|
|
+ molrmsd = outlist[docki]
|
|
+ # Get number of matchs in thresholding operation
|
|
+ pop = population.get(docki)
|
|
+ if not pop: pop = 1
|
|
+ # Save molecule
|
|
+ saveMolecWithRMSD(outsdf, molrmsd[0], molrmsd[1], pop)
|
|
+
|
|
+ if skipped: print("SKIPPED input molecules due to number of atom missmatch: %s"%skipped, file=sys.stderr)
|
|
diff -u -r -N a/build/test/RBT_HOME/check_test.py b/build/test/RBT_HOME/check_test.py
|
|
--- a/build/test/RBT_HOME/check_test.py 2020-10-14 11:48:36.000000000 +0900
|
|
+++ b/build/test/RBT_HOME/check_test.py 2020-10-19 09:23:31.000000000 +0900
|
|
@@ -21,6 +21,6 @@
|
|
error = 1
|
|
|
|
if error == 1:
|
|
- print "The test failed, please check the compilation is OK and no errors were raised."
|
|
+ print ("The test failed, please check the compilation is OK and no errors were raised.")
|
|
else:
|
|
- print "The test succeeded! The results agree with the reference ones.\nHave fun using rDock!!"
|
|
+ print ("The test succeeded! The results agree with the reference ones.\nHave fun using rDock!!")
|