1import sys
2import argparse
3
4import CDPL.Chem as Chem
5import CDPL.Biomol as Biomol
6import CDPL.Pharm as Pharm
7import CDPL.MolProp as MolProp
8
9
10# reads and preprocesses the specified receptor structure
11def readAndPrepareReceptorStructure(args: argparse.Namespace) -> Chem.Molecule:
12 # create reader for receptor structure (format specified by file extension)
13 reader = Chem.MoleculeReader(args.receptor_file)
14
15 sup_fmts = [ Chem.DataFormat.MOL2,
16 Biomol.DataFormat.PDB,
17 Biomol.DataFormat.MMTF,
18 Biomol.DataFormat.MMCIF ]
19
20 if reader.getDataFormat() not in sup_fmts: # check if the format is supported by this script
21 sys.exit('Error: receptor input file format not supported')
22
23 rec_mol = Chem.BasicMolecule() # create an instance of the default implementation of the
24 # Chem.Molecule interface that will store the receptor struct.
25 try:
26 if not reader.read(rec_mol): # read receptor structure
27 sys.exit('Error: reading receptor structure failed')
28
29 except Exception as e:
30 sys.exit('Error: reading receptor structure failed:\n' + str(e))
31
32 # preprocess the receptor structure (removal of residues and
33 # calculation of properties required by the pharm. generation procedure)
34 try:
35 rem_atoms = False
36
37 # delete atoms belonging to residues that should be stripped
38 if args.strip_res_list:
39 atoms_to_rem = Chem.Fragment() # will store the atoms to delete
40 res_to_strip = { tlc.upper() for tlc in args.strip_res_list }
41
42 for atom in rec_mol.atoms: # identify and note atoms belonging to the stripped residues
43 if Biomol.getResidueCode(atom).upper() in res_to_strip:
44 atoms_to_rem.addAtom(atom)
45
46 if atoms_to_rem.numAtoms > 0:
47 rec_mol -= atoms_to_rem # delete atoms from the receptor structure
48 rem_atoms = True
49
50 if not args.quiet:
51 print('! Removed %s atoms from the receptor structure' % str(atoms_to_rem.numAtoms))
52
53 # prepares the receptor structure for pharmacophore generation
54 Chem.perceiveSSSR(rec_mol, rem_atoms)
55 Chem.setRingFlags(rec_mol, rem_atoms)
56 Chem.calcImplicitHydrogenCounts(rec_mol, rem_atoms);
57 Chem.perceiveHybridizationStates(rec_mol, rem_atoms);
58 Chem.setAromaticityFlags(rec_mol, rem_atoms);
59
60 if Chem.makeHydrogenComplete(rec_mol): # make implicit hydrogens (if any) explicit
61 Chem.calcHydrogen3DCoordinates(rec_mol) # calculate 3D coordinates for the added expl. hydrogens
62 Biomol.setHydrogenResidueSequenceInfo(rec_mol, False) # set residue information for the added expl. hydrogens
63
64 MolProp.calcAtomHydrophobicities(rec_mol, False) # calculate atom hydrophobicity values (needed for hydrophobic
65 # pharm. feature generation)
66 except Exception as e:
67 sys.exit('Error: processing of receptor structure failed: ' + str(e))
68
69 return rec_mol
70
71def parseArgs() -> argparse.Namespace:
72 parser = argparse.ArgumentParser(description='Generates pharmacophores describing the interactions between a given receptor structure and a set of ligand molecules.')
73
74 parser.add_argument('-r',
75 dest='receptor_file',
76 required=True,
77 metavar='<file>',
78 help='Receptor structure input file (*.mol2, *.pdb, *.mmtf, *.cif, *.mmcif)')
79 parser.add_argument('-l',
80 dest='ligands_file',
81 required=True,
82 metavar='<file>',
83 help='Ligand structure input file (*.sdf, *.mol2, *.cdf)')
84 parser.add_argument('-o',
85 dest='out_file',
86 required=True,
87 metavar='<file>',
88 help='Pharmacophore output file (*.pml, *.cdf)')
89 parser.add_argument('-s',
90 dest='strip_res_list',
91 required=False,
92 metavar='<three-letter code(s)>',
93 nargs='+',
94 help='Whitespace separated list of PDB three-letter codes specifying residues to remove from the receptor structure (e.g. an existing ligand)')
95 parser.add_argument('-q',
96 dest='quiet',
97 required=False,
98 action='store_true',
99 default=False,
100 help='Disable progress output (default: false)')
101 parser.add_argument('-x',
102 dest='gen_x_vols',
103 required=False,
104 action='store_true',
105 default=False,
106 help='Generate exclusion volumes (default: false)')
107
108 return parser.parse_args()
109
110def main() -> None:
111 args = parseArgs()
112
113 rec_mol = readAndPrepareReceptorStructure(args) # read and preprocess the receptor structure
114 lig_reader = Chem.MoleculeReader(args.ligands_file) # create reader for the ligand input file (format specified by file extension)
115 ph4_writer = Pharm.FeatureContainerWriter(args.out_file) # create writer for the generated pharmacophores (format specified by file extension)
116
117 lig_mol = Chem.BasicMolecule() # create an instance of the default implementation of the
118 # Chem.Molecule interface that will store the ligand structures
119 ia_ph4 = Pharm.BasicPharmacophore() # create an instance of the default implementation of the Pharm.Pharmacophore
120 # interface that will store the generated pharmacophores
121
122 ph4_gen = Pharm.InteractionPharmacophoreGenerator() # create an instance of the pharmacophore generator
123
124 ph4_gen.addExclusionVolumes(args.gen_x_vols) # specify whether to generate exclusion volume spheres
125 # on pharm. feature atoms of interacting residues
126 try:
127 i = 1
128
129 # read and process ligand molecules one after the other until the end of input has been reached (or a severe error occurs)
130 while lig_reader.read(lig_mol):
131 mol_id = Chem.getName(lig_mol).strip() # compose a simple ligand identifier for messages
132
133 if mol_id == '':
134 mol_id = '#' + str(i) # fallback if name is empty or not available
135 else:
136 mol_id = '\'%s\' (#%s)' % (mol_id, str(i))
137
138 if not args.quiet:
139 print('- Generating interaction pharmacophore of molecule %s...' % mol_id)
140
141 try:
142 Pharm.prepareForPharmacophoreGeneration(lig_mol) # make ligand ready for pharm. generation
143
144 ph4_gen.generate(lig_mol, rec_mol, ia_ph4, True) # generate the pharmacophore (True = extract ligand environment residues on-the-fly)
145
146 if not args.quiet:
147 print(' -> Generated %s features: %s' % (str(ia_ph4.numFeatures), Pharm.generateFeatureTypeHistogramString(ia_ph4)))
148
149 try:
150 if not ph4_writer.write(ia_ph4): # output pharmacophore
151 sys.exit('Error: writing interaction pharmacophore of molecule %s failed' % mol_id)
152
153 except Exception as e: # handle exception raised in case of severe write errors
154 sys.exit('Error: writing interaction pharmacophore of molecule %s failed: %s' % (mol_id, str(e)))
155
156 except Exception as e: # handle exception raised in case of severe processing errors
157 sys.exit('Error: interaction pharmacophore generation for molecule %s failed: %s' % (mol_id, str(e)))
158
159 i += 1
160
161 except Exception as e: # handle exception raised in case of severe read errors
162 sys.exit('Error: reading molecule %s failed: %s' % (str(i), str(e)))
163
164 if not args.quiet:
165 print('Done!')
166
167 ph4_writer.close()
168 sys.exit(0)
169
170if __name__ == '__main__':
171 main()