2.3.4. Pharmacophore Feature Mapping Analysis

The script report_ph4_ftr_mpg.py reports information about detected matches between the features of a reference pharmacophore and a set of input pharmacophores.

Synopsis

python report_ph4_ftr_mpg.py [-h] -r <file> -i <file> -o <file> [-x] [-Q] [-q] [-p]

Mandatory options

-r <file>

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

-i <file>

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

-o <file>

Feature mapping info output file

Other options

-h, --help

Show help message and exit

-x

Do not remove exclusion volumes before processing (default: false)

-d <float>

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

-Q

If specified, only matches 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, matches where the position of at least one feature of a pair lies within the tolerance sphere of the other feature are also considered (default: false)

-q

Disable progress output (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 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    return ph4
 24
 25def parseArgs() -> argparse.Namespace:
 26    parser = argparse.ArgumentParser(description='Reports information about the detected matches between the features of a reference pharmacophore and a set of input pharmacophores.')
 27
 28    parser.add_argument('-r',
 29                        dest='ref_ph4_file',
 30                        required=True,
 31                        metavar='<file>',
 32                        help='Reference pharmacophore input file (*.pml, *.cdf)')
 33    parser.add_argument('-i',
 34                        dest='in_file',
 35                        required=True,
 36                        metavar='<file>',
 37                        help='Pharmacophore input file (*.pml, *.cdf)')
 38    parser.add_argument('-o',
 39                        dest='out_file',
 40                        required=True,
 41                        metavar='<file>',
 42                        help='Feature mapping info output file')
 43    parser.add_argument('-Q',
 44                        dest='query_mode',
 45                        required=False,
 46                        action='store_true',
 47                        default=False,
 48                        help='If specified, only matches where the positions of the features of the input pharmacophores lie strictly '\
 49                        'within the tolerance spheres of the reference pharmacophore features will be considered as being valid. Otherwise, '\
 50                        'matches where the position of at least one feature of a pair lies within the tolerance sphere of the '\
 51                        'other feature are also considered (default: false)')
 52    parser.add_argument('-x',
 53                        dest='keep_x_vols',
 54                        required=False,
 55                        action='store_true',
 56                        default=False,
 57                        help='Do not remove exclusion volumes before processing (default: false)')
 58    parser.add_argument('-q',
 59                        dest='quiet',
 60                        required=False,
 61                        action='store_true',
 62                        default=False,
 63                        help='Disable progress output (default: false)')
 64 
 65    return parser.parse_args()
 66
 67def main() -> None:
 68    args = parseArgs()
 69
 70    # read the reference pharmacophore
 71    ref_ph4 = readRefPharmacophore(args.ref_ph4_file) 
 72
 73    # will store the processed reference features
 74    working_ref_ph4 = Pharm.FeatureSet()
 75
 76    if args.keep_x_vols:
 77        working_ref_ph4.assign(ref_ph4)
 78    else:
 79        Pharm.removeFeaturesWithType(ref_ph4, working_ref_ph4, Pharm.FeatureType.EXCLUSION_VOLUME)
 80
 81    # will store the processed input features
 82    working_ipt_ph4 = Pharm.FeatureSet()
 83
 84    # create reader for input pharmacophores (format specified by file extension)
 85    ph4_reader = Pharm.PharmacophoreReader(args.in_file) 
 86
 87    # open output file for writing
 88    out_file = open(args.out_file, 'w')
 89
 90    # write output file CSV column headers
 91    out_file.write('Input Pharm. Index, Ref. Feature Index, Ref. Feature Type, Inp. Feature Index, Inp. Feature Type, Distance Score, Angle Score\n')
 92    
 93    # create an instances of the default implementation of the Pharm.Pharmacophore interface
 94    ph4 = Pharm.BasicPharmacophore()
 95
 96    # create an instance of the class implementing the perception and storage of pharmacophore feature mappings 
 97    ftr_mpg = Pharm.SpatialFeatureMapping(args.query_mode) # arg = True -> input features must be within the tolerance spheres of the ref. features
 98
 99    # create identity transformation for input pharmacophore features (already aligned)
100    id_xform = Math.Matrix4D(Math.DIdentityMatrix(4, 4))
101   
102    # read and process pharmacphores one after the other until the end of input has been reached
103    try:
104        i = 1
105        
106        while ph4_reader.read(ph4):
107            # compose a simple parmacophore identifier
108            ph4_id = Pharm.getName(ph4).strip() 
109
110            if ph4_id == '':
111                ph4_id = f'#{i}'       # fallback if name is empty
112            else:
113                ph4_id = f'\'{ph4_id}\' (#{i})'
114
115            if not args.quiet:
116                print(f'- Processing input pharmacophore {ph4_id}...')
117
118            try:
119                if args.keep_x_vols:
120                    working_ipt_ph4.assign(ph4)
121                else:
122                    Pharm.removeFeaturesWithType(ph4, working_ipt_ph4, Pharm.FeatureType.EXCLUSION_VOLUME)
123                
124                unmpd_ref_ftrs = set(range(0, working_ref_ph4.numFeatures))
125                unmpd_ipt_ftrs = set(range(0, working_ipt_ph4.numFeatures))
126            
127                ftr_mpg.perceive(working_ref_ph4, working_ipt_ph4, id_xform)
128
129                # output found feature mappings where the features are of the same type
130                for ftr_pair in ftr_mpg.items():
131                    out_file.write(f'{i - 1}, {ref_ph4.getFeatureIndex(ftr_pair[0])}, {Pharm.getType(ftr_pair[0])}, {ph4.getFeatureIndex(ftr_pair[1])}, '\
132                                   f'{Pharm.getType(ftr_pair[1])}, {ftr_mpg.getPositionMatchScore(ftr_pair[0], ftr_pair[1]):.4g}, '\
133                                   f'{ftr_mpg.getGeometryMatchScore(ftr_pair[0], ftr_pair[1]):.4g}\n')
134                    
135                    unmpd_ref_ftrs.discard(working_ref_ph4.getFeatureIndex(ftr_pair[0]))
136                    unmpd_ipt_ftrs.discard(working_ipt_ph4.getFeatureIndex(ftr_pair[1]))
137
138                if not args.quiet:
139                    print(f' -> Num. reference features matching input features of the same type: {working_ref_ph4.numFeatures - len(unmpd_ref_ftrs)}')
140                    print(f' -> Num. input features matching reference features of the same type: {working_ipt_ph4.numFeatures - len(unmpd_ipt_ftrs)}')
141                    
142                het_mpd_ref_ftrs = set()
143                het_mpd_ipt_ftrs = set()
144                
145                # output valid (in terms of distance) matches of not yet mapped reference and
146                # input features (ignoring the type)
147                for ref_ftr_idx in list(unmpd_ref_ftrs):
148                    ref_ftr = working_ref_ph4.getFeature(ref_ftr_idx)
149                    
150                    for ipt_ftr_idx in unmpd_ipt_ftrs:
151                        ipt_ftr = working_ipt_ph4.getFeature(ipt_ftr_idx)
152                        dist_score = ftr_mpg.getPositionMatchFunction()(ref_ftr, ipt_ftr, id_xform)
153                        
154                        if dist_score > 0:
155                            unmpd_ref_ftrs.discard(ref_ftr_idx)
156                            het_mpd_ref_ftrs.add(ref_ftr_idx)
157                            het_mpd_ipt_ftrs.add(ipt_ftr_idx)
158                            
159                            out_file.write(f'{i - 1}, {ref_ph4.getFeatureIndex(ref_ftr)}, {Pharm.getType(ref_ftr)}, {ph4.getFeatureIndex(ipt_ftr)}, '\
160                                           f'{Pharm.getType(ipt_ftr)}, {dist_score:.4g}, -1\n')
161
162                if not args.quiet:
163                    print(f' -> Num. remaining reference features matching input features of any type: {len(het_mpd_ref_ftrs)}')
164                    print(f' -> Num. remaining input features matching reference features of any type: {len(het_mpd_ipt_ftrs)}')
165                 
166                unmpd_ipt_ftrs -= het_mpd_ipt_ftrs
167
168                # output unmapped reference features
169                for ref_ftr_idx in unmpd_ref_ftrs:
170                    ref_ftr = working_ref_ph4.getFeature(ref_ftr_idx)
171
172                    out_file.write(f'{i - 1}, {ref_ph4.getFeatureIndex(ref_ftr)}, {Pharm.getType(ref_ftr)}, -1, -1, -1, -1\n')
173
174                # output unmapped input features
175                for ipt_ftr_idx in unmpd_ipt_ftrs:
176                    ipt_ftr = working_ipt_ph4.getFeature(ipt_ftr_idx)
177
178                    out_file.write(f'{i - 1}, -1, -1, {ph4.getFeatureIndex(ipt_ftr)}, {Pharm.getType(ipt_ftr)}, -1, -1\n')
179
180                if not args.quiet:
181                    print(f' -> Num. unmatched reference features: {len(unmpd_ref_ftrs)}')
182                    print(f' -> Num. unmatched input features: {len(unmpd_ipt_ftrs)}')
183             
184            except Exception as e:
185                sys.exit(f'Error: processing of input parmacophore {ph4_id} failed:\n{str(e)}')
186
187            i += 1
188                
189    except Exception as e: # handle exception raised in case of severe read errors
190        sys.exit(f'Error: reading input parmacophore failed:\n{str(e)}')
191
192    out_file.close()
193    sys.exit(0)
194        
195if __name__ == '__main__':
196    main()

Download source file