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