2.3.2. Molecule to Reference Pharmacophore Alignment

The script align_mols_to_ph4.py overlays a set of input molecules with a given reference pharmacophore and outputs the molecules translated/rotated to the calculated alignment pose(s).

Synopsis

python align_mols_to_ph4.py [-h] -r <file> -i <file> -o <file> [-n <integer>] [-x] [-d <float>] [-q] [-p]

Mandatory options

-r <file>

Reference pharmacophore input file (*.pml, *.cdf)

-i <file>

Molecule input file

-o <file>

Aligned molecule output file

Other options

-h, --help

Show help message and exit

-n <integer>

Number of top-ranked alignment solutions to output per molecule (default: best alignment solution only)

-x

Perform an exhaustive alignment search (default: false)

-d <float>

Minimum required RMSD between two consecutively output molecule alignment poses (default: 0.0)

-q

Disable progress output (default: false)

-p

Ignore feature orientations, feature position matching only (default: false)

Code

  1import sys
  2import argparse
  3
  4import CDPL.Chem as Chem
  5import CDPL.Pharm as Pharm
  6import CDPL.Math as Math
  7
  8
  9# reads and returns the specified alignment reference pharmacophore
 10def readRefPharmacophore(filename: str) -> Pharm.Pharmacophore:
 11    # create pharmacophore reader instance
 12    reader = Pharm.PharmacophoreReader(filename)
 13
 14    # create an instance of the default implementation of the Pharm.Pharmacophore interface
 15    ph4 = Pharm.BasicPharmacophore()
 16
 17    try:
 18        if not reader.read(ph4): # read reference pharmacophore
 19            sys.exit('Error: reading reference pharmacophore failed')
 20                
 21    except Exception as e: # handle exception raised in case of severe read errors
 22        sys.exit('Error: reading reference pharmacophore failed:\n' + str(e))
 23
 24    # remove exclusion volumes
 25    Pharm.removeFeaturesWithType(ph4, Pharm.FeatureType.EXCLUSION_VOLUME)
 26
 27    return ph4
 28
 29# generates and returns the pharmacophore of the specified molecule
 30def genPharmacophore(mol: Chem.Molecule) -> Pharm.Pharmacophore:
 31    Pharm.prepareForPharmacophoreGeneration(mol)       # call utility function preparing the molecule for pharmacophore generation
 32        
 33    ph4_gen = Pharm.DefaultPharmacophoreGenerator()    # create an instance of the pharmacophore generator default implementation
 34    ph4 = Pharm.BasicPharmacophore()                   # create an instance of the default implementation of the Pharm.Pharmacophore interface
 35
 36    ph4_gen.generate(mol, ph4)                         # generate the pharmacophore
 37
 38    return ph4
 39
 40def parseArgs() -> argparse.Namespace:
 41    parser = argparse.ArgumentParser(description='Aligns a set of input molecules onto a given reference pharmacophore.')
 42
 43    parser.add_argument('-r',
 44                        dest='ref_ph4_file',
 45                        required=True,
 46                        metavar='<file>',
 47                        help='Reference pharmacophore input file (*.pml, *.cdf)')
 48    parser.add_argument('-i',
 49                        dest='in_file',
 50                        required=True,
 51                        metavar='<file>',
 52                        help='Molecule input file')
 53    parser.add_argument('-o',
 54                        dest='out_file',
 55                        required=True,
 56                        metavar='<file>',
 57                        help='Aligned molecule output file')
 58    parser.add_argument('-n',
 59                        dest='num_out_almnts',
 60                        required=False,
 61                        metavar='<integer>',
 62                        default=1,
 63                        help='Number of top-ranked alignment solutions to output per molecule (default: best alignment solution only)',
 64                        type=int)
 65    parser.add_argument('-x',
 66                        dest='exhaustive',
 67                        required=False,
 68                        action='store_true',
 69                        default=False,
 70                        help='Perform an exhaustive alignment search (default: false)')
 71    parser.add_argument('-d',
 72                        dest='min_pose_rmsd',
 73                        required=False,
 74                        metavar='<float>',
 75                        default=0.0,
 76                        help='Minimum required RMSD between two consecutively output molecule alignment poses (default: 0.0)',
 77                        type=float)
 78    parser.add_argument('-q',
 79                        dest='quiet',
 80                        required=False,
 81                        action='store_true',
 82                        default=False,
 83                        help='Disable progress output (default: false)')
 84    parser.add_argument('-p',
 85                        dest='pos_only',
 86                        required=False,
 87                        action='store_true',
 88                        default=False,
 89                        help='Ignore feature orientations, feature position matching only (default: false)')
 90
 91    return parser.parse_args()
 92
 93def main() -> None:
 94    args = parseArgs()
 95
 96    # read the reference pharmacophore
 97    ref_ph4 = readRefPharmacophore(args.ref_ph4_file) 
 98
 99    # create reader for input molecules (format specified by file extension)
100    mol_reader = Chem.MoleculeReader(args.in_file) 
101
102    Chem.setMultiConfImportParameter(mol_reader, False) # treat conformers as individual molecules
103    
104    # create writer for aligned molecules (format specified by file extension)
105    mol_writer = Chem.MolecularGraphWriter(args.out_file) 
106
107    # create an instance of the default implementation of the Chem.Molecule interface
108    mol = Chem.BasicMolecule()
109
110    # create an instance of the class implementing the pharmacophore alignment algorithm
111    almnt = Pharm.PharmacophoreAlignment(True) # True = aligned features have to be within the tolerance spheres of the ref. features
112
113    if args.pos_only:                          # clear feature orientation information
114        Pharm.clearOrientations(ref_ph4)
115        Pharm.removePositionalDuplicates(ref_ph4)
116    
117    almnt.addFeatures(ref_ph4, True)               # set reference features (True = first set = reference)
118    almnt.performExhaustiveSearch(args.exhaustive) # set minimum number of top. mapped feature pairs
119    
120    # create pharmacophore fit score calculator instance
121    almnt_score = Pharm.PharmacophoreFitScore(True) # True = aligned features have to be within the tolerance spheres of the ref. features
122    
123    # read and process molecules one after the other until the end of input has been reached
124    try:
125        i = 1
126
127        while mol_reader.read(mol):
128            # compose a simple molecule identifier
129            mol_id = Chem.getName(mol).strip() 
130
131            if mol_id == '':
132                mol_id = '#' + str(i)  # fallback if name is empty
133            else:
134                mol_id = '\'%s\' (#%s)' % (mol_id, str(i))
135
136            if not args.quiet:
137                print('- Aligning molecule %s...' % mol_id)
138
139            try:
140                mol_ph4 = genPharmacophore(mol)    # generate input molecule pharmacophore
141
142                if args.pos_only:                  # clear feature orientation information
143                    Pharm.clearOrientations(mol_ph4)
144                    Pharm.removePositionalDuplicates(mol_ph4)
145                    
146                almnt.clearEntities(False)         # clear features of previously aligned pharmacophore
147                almnt.addFeatures(mol_ph4, False)  # specify features of the pharmacophore to align
148
149                almnt_solutions = []               # stores the found alignment solutions
150                
151                while almnt.nextAlignment():                                     # iterate over all alignment solutions that can be found
152                    score = almnt_score(ref_ph4, mol_ph4, almnt.getTransform())  # calculate alignment score
153                    xform = Math.Matrix4D(almnt.getTransform())                  # make a copy of the alignment transformation (mol. ph4 -> ref. ph4) 
154
155                    almnt_solutions.append((score, xform))
156
157                if not args.quiet:
158                    print(' -> Found %s alignment solution(s)' % str(len(almnt_solutions)))
159                
160                saved_coords = Math.Vector3DArray()      # create data structure for storing 3D coordinates
161
162                Chem.get3DCoordinates(mol, saved_coords) # save the original atom coordinates
163
164                struct_data = None
165
166                if Chem.hasStructureData(mol):           # get existing structure data if available
167                    struct_data = Chem.getStructureData(mol)
168                else:                                    # otherwise create and set new structure data
169                    struct_data = Chem.StringDataBlock()
170
171                    Chem.setStructureData(mol, struct_data)
172
173                # add alignment score entry to struct. data
174                struct_data.addEntry('<PharmFitScore>', '') 
175                
176                output_cnt = 0
177                last_pose = None
178                
179                # order solutions by desc. alignment score
180                almnt_solutions = sorted(almnt_solutions, key=lambda entry: entry[0], reverse=True)
181
182                # output molecule alignment poses until the max. number of best output solutions has been reached
183                for solution in almnt_solutions:
184                    if output_cnt == args.num_out_almnts:
185                        break
186
187                    curr_pose = Math.Vector3DArray(saved_coords)
188
189                    Math.transform(curr_pose, solution[1])  # transform atom coordinates
190
191                    # check whether the current pose is 'different enough' from
192                    # the last pose to qualify for output
193                    if args.min_pose_rmsd > 0.0 and last_pose and Math.calcRMSD(last_pose, curr_pose) < args.min_pose_rmsd:
194                        continue
195
196                    # apply the transformed atom coordinates
197                    Chem.set3DCoordinates(mol, curr_pose)  
198
199                    # store alignment score in the struct. data entry
200                    struct_data[len(struct_data) - 1].setData(format(solution[0], '.4f'))     
201                    
202                    try:
203                        if not mol_writer.write(mol): # output the alignment pose of the molecule
204                            sys.exit('Error: writing alignment pose of molecule %s failed' % mol_id)
205
206                    except Exception as e: # handle exception raised in case of severe write errors
207                        sys.exit('Error: writing alignment pose of molecule %s failed:\n%s' % (mol_id, str(e)))
208
209                    last_pose = curr_pose
210                    output_cnt += 1
211
212                if not args.quiet:
213                    print(' -> %s alignment pose(s) saved' % str(output_cnt))
214
215            except Exception as e:
216                sys.exit('Error: pharmacophore alignment of molecule %s failed:\n%s' % (mol_id, str(e)))
217
218            i += 1
219                
220    except Exception as e: # handle exception raised in case of severe read errors
221        sys.exit('Error: reading input molecule failed:\n' + str(e))
222
223    mol_writer.close()
224    sys.exit(0)
225        
226if __name__ == '__main__':
227    main()

Download source file