Mercurial > repos > bgruening > openbabel_subsearch

comparison distance_finder.py @ 15:9adf3fae2771 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/chemicaltoolbox/openbabel commit d9c51279c061a1da948a2582d5b502ca7573adbf
author bgruening
date Thu, 15 Aug 2024 11:04:41 +0000
parents bd678d7db2ae
children
comparison
equal deleted inserted replaced
14:c19608db6b15 15:9adf3fae2771
17 17
18 from openbabel import pybel 18 from openbabel import pybel
19 19
20 20
21 def log(*args, **kwargs): 21 def log(*args, **kwargs):
22 """Log output to STDERR 22 """Log output to STDERR"""
23 """ 23 print(*args, file=sys.stderr, **kwargs)
24 print(*args, file=sys.stderr, ** kwargs)
25 24
26 25
27 def execute(ligands_sdf, points_file, outfile): 26 def execute(ligands_sdf, points_file, outfile):
28 """ 27 """
29 :param ligands_sdf: A SDF with the 3D molecules to test 28 :param ligands_sdf: A SDF with the 3D molecules to test
33 """ 32 """
34 33
35 points = [] 34 points = []
36 35
37 # read the points 36 # read the points
38 with open(points_file, 'r') as f: 37 with open(points_file, "r") as f:
39 for line in f.readlines(): 38 for line in f.readlines():
40 line.strip() 39 line.strip()
41 if line: 40 if line:
42 p = line.split() 41 p = line.split()
43 if len(p) == 3: 42 if len(p) == 3:
44 points.append((float(p[0]), float(p[1]), float(p[2]))) 43 points.append((float(p[0]), float(p[1]), float(p[2])))
45 log("Read points", p) 44 log("Read points", p)
46 continue 45 continue
47 log("Failed to read line:", line) 46 log("Failed to read line:", line)
48 log('Found', len(points), 'atom points') 47 log("Found", len(points), "atom points")
49 48
50 sdf_writer = pybel.Outputfile("sdf", outfile, overwrite=True) 49 sdf_writer = pybel.Outputfile("sdf", outfile, overwrite=True)
51 50
52 count = 0 51 count = 0
53 for mol in pybel.readfile("sdf", ligands_sdf): 52 for mol in pybel.readfile("sdf", ligands_sdf):
54 count += 1 53 count += 1
55 if count % 50000 == 0: 54 if count % 50000 == 0:
56 log('Processed', count) 55 log("Processed", count)
57 56
58 try: 57 try:
59 # print("Processing mol", mol.title) 58 # print("Processing mol", mol.title)
60 clone = pybel.Molecule(mol) 59 clone = pybel.Molecule(mol)
61 clone.removeh() 60 clone.removeh()
68 for point in points: 67 for point in points:
69 p += 1 68 p += 1
70 distances = [] 69 distances = []
71 for i in coords: 70 for i in coords:
72 # calculates distance based on cartesian coordinates 71 # calculates distance based on cartesian coordinates
73 distance = math.sqrt((point[0] - i[0])**2 + (point[1] - i[1])**2 + (point[2] - i[2])**2) 72 distance = math.sqrt(
73 (point[0] - i[0]) ** 2
74 + (point[1] - i[1]) ** 2
75 + (point[2] - i[2]) ** 2
76 )
74 distances.append(distance) 77 distances.append(distance)
75 # log("distance:", distance) 78 # log("distance:", distance)
76 min_distance = min(distances) 79 min_distance = min(distances)
77 # log('Min:', min_distance) 80 # log('Min:', min_distance)
78 # log(count, p, min_distance) 81 # log(count, p, min_distance)
79 82
80 mol.data['distance' + str(p)] = min_distance 83 mol.data["distance" + str(p)] = min_distance
81 84
82 sdf_writer.write(mol) 85 sdf_writer.write(mol)
83 86
84 except Exception as e: 87 except Exception as e:
85 log('Failed to handle molecule: ' + str(e)) 88 log("Failed to handle molecule: " + str(e))
86 continue 89 continue
87 90
88 sdf_writer.close() 91 sdf_writer.close()
89 log('Wrote', count, 'molecules') 92 log("Wrote", count, "molecules")
90 93
91 94
92 def main(): 95 def main():
93 global work_dir 96 global work_dir
94 97
95 parser = argparse.ArgumentParser(description='XChem distances - measure distances to particular points') 98 parser = argparse.ArgumentParser(
96 parser.add_argument('-i', '--input', help="SDF containing the 3D molecules to score)") 99 description="XChem distances - measure distances to particular points"
97 parser.add_argument('-p', '--points', help="PDB format file with atoms") 100 )
98 parser.add_argument('-o', '--outfile', default='output.sdf', help="File name for results") 101 parser.add_argument(
102 "-i", "--input", help="SDF containing the 3D molecules to score)"
103 )
104 parser.add_argument("-p", "--points", help="PDB format file with atoms")
105 parser.add_argument(
106 "-o", "--outfile", default="output.sdf", help="File name for results"
107 )
99 108
100 args = parser.parse_args() 109 args = parser.parse_args()
101 log("XChem distances args: ", args) 110 log("XChem distances args: ", args)
102 111
103 execute(args.input, args.points, args.outfile) 112 execute(args.input, args.points, args.outfile)