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