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