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