1import sys
2import os
3import numpy
4
5import CDPL.Chem as Chem
6
7
8# calculates the FAME descriptor for a given atom of the provided molecule
9def genFAMEDescriptor(ctr_atom: Chem.Atom, molgraph: Chem.MolecularGraph, radius: int) -> numpy.array:
10 env = Chem.Fragment() # for storing of extracted environment atoms
11 descr = numpy.zeros((Chem.SybylAtomType.MAX_TYPE + 1) * (radius + 1))
12
13 Chem.getEnvironment(ctr_atom, molgraph, radius, env) # extract environment of center atom reaching
14 # out up to 'radius' bonds
15 for atom in env.atoms: # iterate over extracted environment atoms
16 sybyl_type = Chem.getSybylType(atom) # retrieve Sybyl type of environment atom
17 top_dist = Chem.getTopologicalDistance(ctr_atom, atom, molgraph) # get top. distance between center atom and environment atom
18 descr[top_dist * (Chem.SybylAtomType.MAX_TYPE + 1) + sybyl_type] += 1 # instead of 1 (= Sybyl type presence) also any other numeric atom
19 # property could be summed up here
20 return descr
21
22# function called for each read molecule
23def procMolecule(molgraph: Chem.MolecularGraph) -> None:
24 Chem.calcImplicitHydrogenCounts(molgraph, False) # calculate implicit hydrogen counts and set corresponding property for all atoms
25 Chem.perceiveHybridizationStates(molgraph, False) # perceive atom hybridization states and set corresponding property for all atoms
26 Chem.perceiveSSSR(molgraph, False) # perceive smallest set of smallest rings and store as Chem.MolecularGraph property
27 Chem.setRingFlags(molgraph, False) # perceive cycles and set corresponding atom and bond properties
28 Chem.setAromaticityFlags(molgraph, False) # perceive aromaticity and set corresponding atom and bond properties
29 Chem.perceiveSybylAtomTypes(molgraph, False) # perceive Sybyl atom types and set corresponding property for all atoms
30 Chem.calcTopologicalDistanceMatrix(molgraph, False) # calculate topological distance matrix and store as Chem.MolecularGraph property
31
32 for atom in molgraph.atoms:
33 descr = genFAMEDescriptor(atom, molgraph, 5) # generate atom environment descriptor using a radius of five bonds
34
35 print(descr)
36
37def main() -> None:
38 if len(sys.argv) < 2:
39 sys.exit('Usage: %s <input mol. file>' % sys.argv[0])
40
41 # create reader for input molecules (format specified by file extension)
42 reader = Chem.MoleculeReader(sys.argv[1])
43
44 # create an instance of the default implementation of the Chem.Molecule interface
45 mol = Chem.BasicMolecule()
46
47 # read and process molecules one after the other until the end of input has been reached
48 try:
49 while reader.read(mol):
50 try:
51 procMolecule(mol)
52 except Exception as e:
53 sys.exit('Error: processing of molecule failed: ' + str(e))
54
55 except Exception as e: # handle exception raised in case of severe read errors
56 sys.exit('Error: reading molecule failed: ' + str(e))
57
58 sys.exit(0)
59
60if __name__ == '__main__':
61 main()