2.3.3. Pharmacophore to Pharmacophore Alignment

The script align_ph4s_to_ph4.py overlays a set of input pharmacophores with a given reference pharmacophore and then writes the calculated alignment pose(s) to a specified output file.

Synopsis

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

Mandatory options

-r <file>

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

-i <file>

Pharmacophore input file (*.pml, *.cdf)

-o <file>

Aligned pharmacophore output file (*.pml, *.cdf)

Other options

-h, --help

Show help message and exit

-s <file>

Pharmacophore alignment score output file

-n <integer>

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

-x

Perform an exhaustive alignment search (default: false)

-d <float>

Minimum required score difference between two consecutively output pharmacophore alignment poses (default: 0.0)

-Q

If specified, only alignments where the positions of the features of the input pharmacophores lie strictly within the tolerance spheres of the reference pharmacophore features will be considered as being valid. Otherwise, alignments where the position of at least one feature of the aligned pairs lies within the tolerance sphere of the other feature are also valid (default: false)

-q

Disable progress output (default: false)

-p

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

Code

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

Download source file