1.1.3.3. Database Cleaning
The script clean_mol_db.py reads molecules from an input file, performs (optional) preprocessing, and then writes only those molecules that fulfill particular user-defined criteria to an output file.
Synopsis
python clean_mol_db.py [-h] -i <file> -o <file> [-d <file>] [-s] [-c] [-x <element list>] [-a <element list>] [-m <element count list>] [-M <element count list>] [-v <0|1|2|3>]
Mandatory options
- -i <file>
Input molecule file
- -o <file>
Output molecule file
Other options
- -h, --help
Show help message and exit
- -d <file>
Discarded molecule output file
- -s
Keep only the largest molecule component (default: false)
- -S
Remove atom and bond stereochemistry information (default: false)
- -r
Replace radical electrons by bonds to implicit hydrogens (default: false)
- -c
Minimize the number of charged atoms (default: false) by protonation/deprotonation and charge equalization
- -C <true|false>
Canonicalize output molecule (default: true)
- -x <element list>
List of excluded chem. elements (default: no elements are excluded)
- -a <element list>
List of allowed chem. elements (default: all elements are allowed)
- -m <element count list>
Minimum chem. element specific atom counts (default: no count limits)
- -M <element count list>
Maximum chem. element specific atom counts (default: no count limits)
- -v <0|1|2|3>
Verbosity level (default: 1; 0 -> no console output, 1 -> print summary, 2 -> verbose, 3 -> extra verbose)
The options -a and -x both require a list of chemical elements as argument. Chemical element lists are specified in the form <S>,…,<S> where <S> is the symbol of a chemical element or generic atom type. Supported generic types are:
Symbol |
Meaning |
|---|---|
M |
any metal |
MH |
any metal or hydrogen |
A |
any element except hydrogen |
AH |
any element |
* |
any element (equivalent to AH) |
X |
any halogen |
XH |
any halogen or hydrogen |
Q |
any element except hydrogen and carbon |
QH |
any element except carbon |
The options -m and -M both require a list of chemical element counts as argument.
Chemical element counts are specified in the form <S>:<N>,…,<S>:<N> where <S> is
the symbol of a chemical element or generic atom type (see above) and <N> is
the corresponding minimum or maximum count. If the count part is omitted and only <S>
gets specified then the count is assumed to be 1.
Example usage
python clean_mol_db.py -i <path/to/molecule/input/file> -o <path/to/molecule/output/file> -a C,H,N,O,S,P,F,Cl,Br,I -m C,A:3 -M F:9 -c -s -r
When executed as shown, the script will perform the following operations on each read input molecule (in the order listed):
Radical electrons replaced by bonds to implicit hydrogens (if any)
Reduction of the number of charged atoms (if any and if possible)
Removal of all but the largest molecular graph component (only if multi-comp. molecule)
Check whether the chem. element of each atom of the working molecule (= result of prev. steps) is either C, H, N, O, S, P, F, Cl, Br, or I.
Check whether the atom list of the working molecule contains at least one carbon and three heavy atoms
Check whether the atom list of the working molecule contains not more than 9 fluorine atoms
The first check that fails leads to a rejection of the molecule. Working molecules that pass all checks will be written to the specified output file.
Code
1import sys
2import argparse
3
4import CDPL.Chem as Chem
5import CDPL.MolProp as MolProp
6
7
8# performs a (optional) standardization of the argument molecule and then checks whether
9# it fulfills certain user-defined criteria for forwarding/rejection
10def processMolecule(mol: Chem.Molecule, args: argparse.Namespace) -> tuple:
11 derad = False
12
13 if args.deradicalize:
14 # will store the 'deradicalized' version of the argument mol
15 res_mol = Chem.BasicMolecule(mol)
16
17 # calculates the number of electrons for each atom that are not engaged in a bond although
18 # they should be according to the valence rules
19 Chem.calcUnpairedElectronCounts(res_mol, False)
20
21 # clear radical information
22 for atom in res_mol.atoms:
23 if Chem.getRadicalType(atom) != Chem.RadicalType.NOT_RADICALIC or Chem.getUnpairedElectronCount(atom) > 0:
24 Chem.clearRadicalType(atom)
25 Chem.clearUnpairedElectronCount(atom)
26 Chem.clearImplicitHydrogenCount(atom)
27 derad = True
28
29 if derad:
30 mol = res_mol
31
32 chgs_mod = False
33
34 if args.min_charges:
35 # will store the 'uncharged' version of the argument mol
36 res_mol = Chem.BasicMolecule()
37
38 # minimize the number of charged atoms using the corresponding functionality implemented
39 # in class Chem.ProtonationStateStandardizer
40 chgs_mod = Chem.ProtonationStateStandardizer().standardize(mol, res_mol, Chem.ProtonationStateStandardizer.MIN_CHARGED_ATOM_COUNT)
41
42 # if changes were made then from now on use the 'uncharged' molecule for further checks
43 if chgs_mod:
44 mol = res_mol
45
46 comps_strpd = False
47
48 if args.strip_comps:
49 # perceive components (if not done already)
50 comps = Chem.perceiveComponents(mol, False)
51
52 # find largest component (regarding heavy atom count) but only if multi-comp. molecule
53 if comps.size > 1: # check if multi-comp. molecule
54 lgst_comp = None
55 lgst_comp_hvy_atom_count = 0
56
57 for comp in Chem.getComponents(mol): # for each component
58 hvy_atom_count = MolProp.getHeavyAtomCount(comp) # calc. non-hydrogen atom count
59
60 if hvy_atom_count > lgst_comp_hvy_atom_count: # check if the largest comp. so far has been found
61 lgst_comp_hvy_atom_count = hvy_atom_count # if so, store for later use
62 lgst_comp = comp
63
64 # if the input molecule has structure data then pass them on (for SDF output)
65 if Chem.hasStructureData(mol):
66 Chem.setStructureData(lgst_comp, Chem.getStructureData(mol))
67
68 # for further checks use only the identified largest component
69 comps_strpd = True
70 mol = lgst_comp
71
72 # calc. implicit hydrogen counts (if not done already)
73 Chem.calcImplicitHydrogenCounts(mol, False)
74
75 # create instance of class MolProp.ElementHistogram for storing the per-element atom counts
76 # of the molecule (or its largest comp.)
77 elem_histo = MolProp.ElementHistogram()
78
79 # get per-element atom counts
80 MolProp.generateElementHistogram(mol, elem_histo, False)
81
82 # check if the found chem. elements are all in the set of allowed elements (if specified)
83 if not checkAllowedElements(elem_histo, args.allowed_elements):
84 return (None, 'prohibited element detected')
85
86 # check if none of the found chem. elements is in the set of excluded elements (if specified)
87 if not checkExcludedElements(elem_histo, args.excluded_elements):
88 return (None, 'prohibited element detected')
89
90 # check if all specified minium chem. element atom counts are reached
91 if not checkMinAtomCounts(elem_histo, args.min_atom_counts):
92 return (None, 'minimum atom counts not reached')
93
94 # check if none of the specified maximum chem. element atom counts gets exceeded
95 if not checkMaxAtomCounts(elem_histo, args.max_atom_counts):
96 return (None, 'maximum atom count exceeded')
97
98 if chgs_mod:
99 log_msg = 'reduced charged atom count'
100 else:
101 log_msg = ''
102
103 if derad:
104 if log_msg:
105 log_msg += ', '
106
107 log_msg += 'removed radical information'
108
109 if comps_strpd:
110 if log_msg:
111 log_msg += ', '
112
113 log_msg += 'removed components'
114
115 if args.strip_stereo:
116 for atom in mol.atoms: # clear atom stereo descriptors
117 Chem.clearStereoDescriptor(atom)
118
119 for bond in mol.bonds: # clear bond stereo descriptors
120 Chem.clearStereoDescriptor(bond)
121
122 if args.canonicalize:
123 Chem.calcBasicProperties(mol, False) # calc. basic properties required for canonicalization
124 Chem.calcCanonicalNumbering(mol, True) # calc. canon. numbering of atoms
125 Chem.canonicalize(mol) # order atoms/bonds according to canon. numbering
126
127 return (mol, log_msg)
128
129# checks if all found chem. elements are in the set of allowed elements (if specified)
130def checkAllowedElements(elem_histo: MolProp.ElementHistogram, elem_list: list) -> bool:
131 if not elem_list:
132 return True
133
134 for atom_type in elem_histo.getKeys():
135 allowed = False
136
137 for all_atom_type in elem_list:
138 if Chem.atomTypesMatch(all_atom_type, atom_type):
139 allowed = True
140 break
141
142 if not allowed:
143 return False
144
145 return True
146
147# checks if none of the found chem. elements is in the set of excluded elements (if specified)
148def checkExcludedElements(elem_histo: MolProp.ElementHistogram, elem_list: list) -> bool:
149 if not elem_list:
150 return True
151
152 for atom_type in elem_histo.getKeys():
153 for x_atom_type in elem_list:
154 if Chem.atomTypesMatch(x_atom_type, atom_type):
155 return False
156
157 return True
158
159# return the total number of atoms matching the specified atom type (element or generic type)
160def getNumAtomsOfType(elem_histo: MolProp.ElementHistogram, qry_atom_type: int) -> int:
161 tot_count = 0
162
163 for atom_type, count in elem_histo.getEntries():
164 # check if the elem. histogram entry matches the specified atom type (element or generic type)
165 # if so, add the stored count the total count
166 if Chem.atomTypesMatch(qry_atom_type, atom_type):
167 tot_count += count
168
169 return tot_count
170
171# checks if all the specified minium counts of atoms matching a particular element or
172# generic type are reached
173def checkMinAtomCounts(elem_histo: MolProp.ElementHistogram, atom_counts: dict) -> bool:
174 if not atom_counts:
175 return True
176
177 for atom_type, min_count in atom_counts.items():
178 if getNumAtomsOfType(elem_histo, atom_type) < min_count:
179 return False
180
181 return True
182
183# checks if none of the specified maximum counts of atoms matching a particular element or
184# generic type gets exceeded
185def checkMaxAtomCounts(elem_histo: MolProp.ElementHistogram, atom_counts: dict) -> bool:
186 if not atom_counts:
187 return True
188
189 for atom_type, max_count in atom_counts.items():
190 if getNumAtomsOfType(elem_histo, atom_type) > max_count:
191 return False
192
193 return True
194
195# converts a comma separated list of chem. element/generic type symbols into a list of the
196# corresponding numeric atom types defined in Chem.AtomType
197def parseElementList(elem_list: str) -> list:
198 if not elem_list:
199 return []
200
201 atom_types = []
202
203 for elem in elem_list.split(','):
204 elem = elem.strip()
205
206 if not elem: # zero length?
207 continue
208
209 atom_type = Chem.AtomDictionary.getType(elem)
210
211 if atom_type == Chem.AtomType.UNKNOWN:
212 sys.exit('Error: unknown chemical element/generic type \'%s\'' % elem)
213
214 atom_types.append(atom_type)
215
216 return atom_types
217
218# converts a comma separated list of chem. element (or generic type) symbol/count pairs (sep. by colon) into a dictionary
219# mapping the corresponding numeric atom types (defined in Chem.AtomType) to an integer number
220def parseElementCountList(elem_count_list: str) -> dict:
221 if not elem_count_list:
222 return {}
223
224 atom_type_counts = {}
225
226 for elem_count in elem_count_list.split(','):
227 elem_count = elem_count.strip()
228
229 if not elem_count: # zero length?
230 continue
231
232 elem_count = elem_count.split(':')
233 elem_count[0] = elem_count[0].strip()
234 atom_type = Chem.AtomDictionary.getType(elem_count[0])
235
236 if atom_type == Chem.AtomType.UNKNOWN:
237 sys.exit('Error: unknown chemical element/generic type \'%s\'' % elem_count[0])
238
239 count = 1
240
241 if len(elem_count) > 1: # has count spec.?
242 count = int(elem_count[1])
243
244 atom_type_counts[atom_type] = count
245
246 return atom_type_counts
247
248def parseArgs() -> argparse.Namespace:
249 def strtobool(value: str) -> bool:
250 value = value.lower()
251
252 if value in ('y', 'yes', 'on', '1', 'true', 't'):
253 return True
254
255 return False
256
257 parser = argparse.ArgumentParser(description='Strips compounds that fulfill particular user-defined criteria from a molecule database.')
258
259 parser.add_argument('-i',
260 dest='in_file',
261 required=True,
262 metavar='<file>',
263 help='Input molecule file')
264 parser.add_argument('-o',
265 dest='out_file',
266 required=True,
267 metavar='<file>',
268 help='Output molecule file')
269 parser.add_argument('-d',
270 dest='disc_file',
271 required=False,
272 metavar='<file>',
273 help='Discarded molecule output file')
274 parser.add_argument('-s',
275 dest='strip_comps',
276 required=False,
277 action='store_true',
278 default=False,
279 help='Keep only the largest molecule component (default: false)')
280 parser.add_argument('-S',
281 dest='strip_stereo',
282 required=False,
283 action='store_true',
284 default=False,
285 help='Remove atom and bond stereochemistry information (default: false)')
286 parser.add_argument('-c',
287 dest='min_charges',
288 required=False,
289 action='store_true',
290 default=False,
291 help='Minimize number of charged atoms (default: false)')
292 parser.add_argument('-r',
293 dest='deradicalize',
294 required=False,
295 action='store_true',
296 default=False,
297 help='Replace radical electrons by bonds to implicit hydrogens (default: false)')
298 parser.add_argument('-C',
299 dest='canonicalize',
300 required=False,
301 metavar='<true|false>',
302 type=lambda x:bool(strtobool(x)),
303 default=True,
304 help='Canonicalize output molecule (default: true)')
305 parser.add_argument('-x',
306 dest='excluded_elements',
307 required=False,
308 metavar="<element list>",
309 help='List of excluded chem. elements (default: no elements are excluded)')
310 parser.add_argument('-a',
311 dest='allowed_elements',
312 required=False,
313 metavar="<element list>",
314 help='List of allowed chem. elements (default: all elements are allowed)')
315 parser.add_argument('-m',
316 dest='min_atom_counts',
317 required=False,
318 metavar="<element count list>",
319 help='Minimum chem. element specific atom counts (default: no count limits)')
320 parser.add_argument('-M',
321 dest='max_atom_counts',
322 required=False,
323 metavar="<element count list>",
324 help='Maximum chem. element specific atom counts (default: no count limits)')
325 parser.add_argument('-v',
326 dest='verb_level',
327 required=False,
328 metavar='<0|1|2|3>',
329 choices=range(0, 4),
330 default=1,
331 help='Verbosity level (default: 1; 0 -> no console output, 1 -> print summary, 2 -> verbose, 3 -> extra verbose)',
332 type=int)
333
334 return parser.parse_args()
335
336def main() -> None:
337 args = parseArgs() # process command line arguments
338
339 # convert specified lists of chem. element/gen. type symbols into a corresponding list of
340 # numeric atom types (defined in Chem.Atomtype)
341 args.allowed_elements = parseElementList(args.allowed_elements)
342 args.excluded_elements = parseElementList(args.excluded_elements)
343
344 # convert specified comma separated lists of chem. element (or generic type) symbol/count pairs (sep. by colon) into
345 # corresponding dictionaries mapping numeric atom types (defined in Chem.AtomType) to integers
346 args.min_atom_counts = parseElementCountList(args.min_atom_counts)
347 args.max_atom_counts = parseElementCountList(args.max_atom_counts)
348
349 # create reader for input molecules (format specified by file extension)
350 reader = Chem.MoleculeReader(args.in_file)
351
352 # create writer for molecules passing the checks (format specified by file extension)
353 writer = Chem.MolecularGraphWriter(args.out_file)
354
355 # if canonicalization is enabled and the output format is SMILES then generate canonical SMILES directly
356 if args.canonicalize and (writer.getDataFormat() == Chem.DataFormat.SMILES or \
357 writer.getDataFormat() == Chem.DataFormat.SMILES_GZ or \
358 writer.getDataFormat() == Chem.DataFormat.SMILES_BZ2):
359 args.canonicalize = False
360 Chem.setSMILESOutputCanonicalFormParameter(writer, True)
361
362 if args.disc_file:
363 # create writer for sorted out molecules (format specified by file extension)
364 disc_writer = Chem.MolecularGraphWriter(args.disc_file)
365 else:
366 disc_writer = None
367
368 # create instances of the default implementation of the Chem.Molecule interface for the input and output molecules
369 in_mol = Chem.BasicMolecule()
370
371 i = 0
372 num_changed = 0
373 num_disc = 0
374
375 try:
376 # read and process molecules one after the other until the end of input has been reached (or a severe error occurs)
377 while reader.read(in_mol):
378 # compose a molecule identifier
379 mol_id = Chem.getName(in_mol).strip()
380
381 if mol_id == '':
382 mol_id = '#' + str(i + 1) # fallback if name is empty or not available
383 else:
384 mol_id = '\'%s\' (#%s)' % (mol_id, str(i + 1))
385
386 try:
387 # process input molecule
388 out_mol, log_msg = processMolecule(in_mol, args)
389
390 if out_mol is None: # check whether the molecule has been sorted out
391 if args.verb_level > 1 and log_msg:
392 print('- Molecule %s: discarded, %s' % (mol_id, log_msg))
393
394 num_disc += 1
395
396 if not disc_writer: # check whether discarded molecules should be saved to a separate file
397 i += 1
398 continue
399
400 # write discarded molecule to the specified target file
401 out_mol = in_mol
402 out_mol_writer = disc_writer
403 else:
404 if log_msg:
405 if args.verb_level > 1 :
406 print('- Molecule %s: %s' % (mol_id, log_msg))
407
408 num_changed += 1
409 else:
410 if args.verb_level > 2:
411 print('- Molecule %s: left unchanged' % mol_id)
412
413 # molecule passed all checks and needs to be written to the regular output file
414 out_mol_writer = writer
415
416 try:
417 # calculate (if not already present) some basic properties of the output molecule
418 # that might be required for writing (output format dependent)
419 Chem.calcImplicitHydrogenCounts(out_mol, False)
420 Chem.perceiveHybridizationStates(out_mol, False)
421 Chem.perceiveSSSR(out_mol, False)
422 Chem.setRingFlags(out_mol, False)
423 Chem.setAromaticityFlags(out_mol, False)
424 Chem.perceiveComponents(out_mol, False)
425
426 # output molecule
427 if not out_mol_writer.write(out_mol):
428 sys.exit('Error: writing molecule %s failed' % mol_id)
429
430 except Exception as e: # handle exception raised in case of severe write errors
431 sys.exit('Error: writing molecule %s failed:\n%s' % (mol_id, str(e)))
432
433 i += 1
434
435 except Exception as e: # handle exception raised in case of severe processing errors
436 sys.exit('Error: processing of molecule %s failed:\n%s' % (mol_id, str(e)))
437
438 except Exception as e: # handle exception raised in case of severe read errors
439 sys.exit('Error: reading of molecule %s failed:\n%s' % (str(i), str(e)))
440
441 if args.verb_level > 0:
442 print('Summary:')
443 print(' - %s molecules processed' % str(i))
444 print(' - %s modified ' % str(num_changed))
445 print(' - %s discarded' % str(num_disc))
446
447 writer.close()
448 sys.exit(0)
449
450if __name__ == '__main__':
451 main()