1import sys
2import argparse
3
4import CDPL.Chem as Chem
5import CDPL.Pharm as Pharm
6import CDPL.Math as Math
7
8
9def parseArgs() -> argparse.Namespace:
10 parser = argparse.ArgumentParser(description='Aligns a set of input molecules onto a given reference pharmacophore.')
11
12 parser.add_argument('-r',
13 dest='ref_ph4_file',
14 required=True,
15 metavar='<file>',
16 help='Reference pharmacophore input file (*.pml, *.cdf)')
17 parser.add_argument('-i',
18 dest='in_file',
19 required=True,
20 metavar='<file>',
21 help='Molecule input file')
22 parser.add_argument('-o',
23 dest='out_file',
24 required=True,
25 metavar='<file>',
26 help='Aligned molecule output file')
27 parser.add_argument('-n',
28 dest='num_out_almnts',
29 required=False,
30 metavar='<integer>',
31 default=1,
32 help='Number of top-ranked alignment solutions to output per molecule (default: best alignment solution only)',
33 type=int)
34 parser.add_argument('-x',
35 dest='exhaustive',
36 required=False,
37 action='store_true',
38 default=False,
39 help='Perform an exhaustive alignment search (default: false)')
40 parser.add_argument('-d',
41 dest='min_pose_rmsd',
42 required=False,
43 metavar='<float>',
44 default=0.0,
45 help='Minimum required RMSD between two consecutively output molecule alignment poses (default: 0.0)',
46 type=float)
47 parser.add_argument('-q',
48 dest='quiet',
49 required=False,
50 action='store_true',
51 default=False,
52 help='Disable progress output (default: false)')
53 parser.add_argument('-p',
54 dest='pos_only',
55 required=False,
56 action='store_true',
57 default=False,
58 help='Ignore feature orientations, feature position matching only (default: false)')
59
60 return parser.parse_args()
61
62# reads and returns the specified alignment reference pharmacophore
63def readRefPharmacophore(filename: str) -> Pharm.Pharmacophore:
64 # create pharmacophore reader instance
65 reader = Pharm.PharmacophoreReader(filename)
66
67 # create an instance of the default implementation of the Pharm.Pharmacophore interface
68 ph4 = Pharm.BasicPharmacophore()
69
70 try:
71 if not reader.read(ph4): # read reference pharmacophore
72 sys.exit('Error: reading reference pharmacophore failed')
73
74 except Exception as e: # handle exception raised in case of severe read errors
75 sys.exit('Error: reading reference pharmacophore failed: ' + str(e))
76
77 return ph4
78
79# generates and returns the pharmacophore of the specified molecule
80def genPharmacophore(mol: Chem.Molecule) -> Pharm.Pharmacophore:
81 Pharm.prepareForPharmacophoreGeneration(mol) # call utility function preparing the molecule for pharmacophore generation
82
83 ph4_gen = Pharm.DefaultPharmacophoreGenerator() # create an instance of the pharmacophore generator default implementation
84 ph4 = Pharm.BasicPharmacophore() # create an instance of the default implementation of the Pharm.Pharmacophore interface
85
86 ph4_gen.generate(mol, ph4) # generate the pharmacophore
87
88 return ph4
89
90# remove feature orientation informations and set the feature geometry to Pharm.FeatureGeometry.SPHERE
91def clearFeatureOrientations(ph4: Pharm.BasicPharmacophore) -> None:
92 for ftr in ph4:
93 Pharm.clearOrientation(ftr)
94 Pharm.setGeometry(ftr, Pharm.FeatureGeometry.SPHERE)
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, strut)
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: %s' % (mol_id, str(e)))
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()