1import sys
2import os
3import argparse
4
5import CDPL.Chem as Chem
6import CDPL.Pharm as Pharm
7import CDPL.Math as Math
8
9
10def parseArgs() -> argparse.Namespace:
11 parser = argparse.ArgumentParser(description='Aligns a set of input molecules onto a given reference pharmacophore.')
12
13 parser.add_argument('-r',
14 dest='ref_ph4_file',
15 required=True,
16 metavar='<file>',
17 help='Reference pharmacophore input file (*.pml, *.cdf)')
18 parser.add_argument('-i',
19 dest='in_file',
20 required=True,
21 metavar='<file>',
22 help='Molecule input file')
23 parser.add_argument('-o',
24 dest='out_file',
25 required=True,
26 metavar='<file>',
27 help='Aligned molecule output file (*.pml, *.cdf)')
28 parser.add_argument('-n',
29 dest='num_out_almnts',
30 required=False,
31 metavar='<integer>',
32 default=1,
33 help='Number of top-ranked alignment solutions to output per molecule (default: best alignment solution only)',
34 type=int)
35 parser.add_argument('-x',
36 dest='exhaustive',
37 required=False,
38 action='store_true',
39 default=False,
40 help='Perform an exhaustive alignment search (default: false)')
41 parser.add_argument('-d',
42 dest='min_pose_rmsd',
43 required=False,
44 metavar='<float>',
45 default=0.0,
46 help='Minimum required RMSD between two consecutively output molecule alignment poses (default: 0.0)',
47 type=float)
48 parser.add_argument('-q',
49 dest='quiet',
50 required=False,
51 action='store_true',
52 default=False,
53 help='Disable progress output (default: false)')
54 parser.add_argument('-p',
55 dest='pos_only',
56 required=False,
57 action='store_true',
58 default=False,
59 help='Ignore feature orientations, feature position matching only (default: false)')
60
61 return parser.parse_args()
62
63# reads and returns the specified alignment reference pharmacophore
64def readRefPharmacophore(filename: str) -> Pharm.Pharmacophore:
65 # create pharmacophore reader instance
66 reader = Pharm.PharmacophoreReader(filename)
67
68 # create an instance of the default implementation of the Pharm.Pharmacophore interface
69 ph4 = Pharm.BasicPharmacophore()
70
71 try:
72 if not reader.read(ph4): # read reference pharmacophore
73 sys.exit('Error: reading reference pharmacophore failed')
74
75 except Exception as e: # handle exception raised in case of severe read errors
76 sys.exit('Error: reading reference pharmacophore failed: ' + str(e))
77
78 return ph4
79
80# generates and returns the pharmacophore of the specified molecule
81def genPharmacophore(mol: Chem.Molecule) -> Pharm.Pharmacophore:
82 Pharm.prepareForPharmacophoreGeneration(mol) # call utility function preparing the molecule for pharmacophore generation
83
84 ph4_gen = Pharm.DefaultPharmacophoreGenerator() # create an instance of the pharmacophore generator default implementation
85 ph4 = Pharm.BasicPharmacophore() # create an instance of the default implementation of the Pharm.Pharmacophore interface
86
87 ph4_gen.generate(mol, ph4) # generate the pharmacophore
88
89 return ph4
90
91# remove feature orientation informations and set the feature geometry to Pharm.FeatureGeometry.SPHERE
92def clearFeatureOrientations(ph4: Pharm.BasicPharmacophore) -> None:
93 for ftr in ph4:
94 Pharm.clearOrientation(ftr)
95 Pharm.setGeometry(ftr, Pharm.FeatureGeometry.SPHERE)
96
97def main() -> None:
98 args = parseArgs()
99
100 # read the reference pharmacophore
101 ref_ph4 = readRefPharmacophore(args.ref_ph4_file)
102
103 # create reader for input molecules (format specified by file extension)
104 mol_reader = Chem.MoleculeReader(args.in_file)
105
106 Chem.setMultiConfImportParameter(mol_reader, False) # treat conformers as individual molecules
107
108 # create writer for aligned molecules (format specified by file extension)
109 mol_writer = Chem.MolecularGraphWriter(args.out_file)
110
111 # create an instance of the default implementation of the Chem.Molecule interface
112 mol = Chem.BasicMolecule()
113
114 # create instance of class implementing the pharmacophore alignment algorithm
115 almnt = Pharm.PharmacophoreAlignment(True) # True = aligned features have to be within the tolerance spheres of the ref. features
116
117 if args.pos_only: # clear feature orientation information
118 clearFeatureOrientations(ref_ph4)
119
120 almnt.addFeatures(ref_ph4, True) # set reference features (True = first set = reference)
121 almnt.performExhaustiveSearch(args.exhaustive) # set minimum number of top. mapped feature pairs
122
123 # create pharmacophore fit score calculator instance
124 almnt_score = Pharm.PharmacophoreFitScore()
125
126 # read and process molecules one after the other until the end of input has been reached
127 try:
128 i = 1
129
130 while mol_reader.read(mol):
131 # compose a simple molecule identifier
132 mol_id = Chem.getName(mol).strip()
133
134 if mol_id == '':
135 mol_id = '#' + str(i) # fallback if name is empty
136 else:
137 mol_id = '\'%s\' (#%s)' % (mol_id, str(i))
138
139 if not args.quiet:
140 print('- Aligning molecule %s...' % mol_id)
141
142 try:
143 mol_ph4 = genPharmacophore(mol) # generate input molecule pharmacophore
144
145 if args.pos_only: # clear feature orientation information
146 clearFeatureOrientations(mol_ph4)
147
148 almnt.clearEntities(False) # clear features of previously aligned pharmacophore
149 almnt.addFeatures(mol_ph4, False) # specify features of the pharmacophore to align
150
151 almnt_solutions = [] # stores the found alignment solutions
152
153 while almnt.nextAlignment(): # iterate over all alignment solutions that can be found
154 score = almnt_score(ref_ph4, mol_ph4, almnt.getTransform()) # calculate alignment score
155 xform = Math.Matrix4D(almnt.getTransform()) # make a copy of the alignment transformation (mol. ph4 -> ref. ph4)
156
157 almnt_solutions.append((score, xform))
158
159 if not args.quiet:
160 print(' -> Found %s alignment solutions' % str(len(almnt_solutions)))
161
162 saved_coords = Math.Vector3DArray() # create data structure for storing 3D coordinates
163
164 Chem.get3DCoordinates(mol, saved_coords) # save the original atom coordinates
165
166 struct_data = None
167
168 if Chem.hasStructureData(mol): # get existing structure data if available
169 struct_data = Chem.getStructureData(mol)
170 else: # otherwise create and set new structure data
171 struct_data = Chem.StringDataBlock()
172
173 Chem.setStructureData(mol, strut)
174
175 # add alignment score entry to struct. data
176 struct_data.addEntry('<PharmFitScore>', '')
177
178 output_cnt = 0
179 last_pose = None
180
181 # order solutions by desc. alignment score
182 almnt_solutions = sorted(almnt_solutions, key=lambda entry: entry[0], reverse=True)
183
184 # output molecule alignment poses until the max. number of best output solutions has been reached
185 for solution in almnt_solutions:
186 if output_cnt == args.num_out_almnts:
187 break
188
189 curr_pose = Math.Vector3DArray(saved_coords)
190
191 Math.transform(curr_pose, solution[1]) # transform atom coordinates
192
193 # check whether the current pose is 'different enough' from
194 # the last pose to qualify for output
195 if args.min_pose_rmsd > 0.0 and last_pose and Math.calcRMSD(last_pose, curr_pose) < args.min_pose_rmsd:
196 continue
197
198 # apply the transformed atom coordinates
199 Chem.set3DCoordinates(mol, curr_pose)
200
201 # store alignment score in the struct. data entry
202 struct_data[len(struct_data) - 1].setData(format(solution[0], '.4f'))
203
204 try:
205 if not mol_writer.write(mol): # output the alignment pose of the molecule
206 sys.exit('Error: writing alignment pose of molecule %s failed: %s' % (mol_id, str(e)))
207
208 except Exception as e: # handle exception raised in case of severe write errors
209 sys.exit('Error: writing alignment pose of molecule %s failed: %s' % (mol_id, str(e)))
210
211 last_pose = curr_pose
212 output_cnt += 1
213
214 if not args.quiet:
215 print(' -> %s alignment poses saved' % str(output_cnt))
216
217 except Exception as e:
218 sys.exit('Error: pharmacophore alignment of molecule %s failed: %s' % (mol_id, str(e)))
219
220 i += 1
221
222 except Exception as e: # handle exception raised in case of severe read errors
223 sys.exit('Error: reading input molecule failed: ' + str(e))
224
225 mol_writer.close()
226 sys.exit(0)
227
228if __name__ == '__main__':
229 main()