2056
J. Chem. Inf. Model. 2009, 49, 2056–2066
LigMatch: A Multiple Structure-Based Ligand Matching Method for 3D Virtual Screening Sarah L. Kinnings and Richard M. Jackson* Institute of Molecular and Cellular Biology and Astbury Centre for Structural Molecular Biology, Faculty of Biological Sciences, University of Leeds, Leeds LS2 9JT, United Kingdom Received June 10, 2009
We have developed a new virtual screening (VS) method called LigMatch and evaluated its performance on 13 protein targets using a filtered and clustered version of the directory of useful decoys (DUD). The method uses 3D structural comparison to a crystallographically determined ligand in a bioactive ‘template’ conformation, using a geometric hashing method, in order to prioritize each database compound. We show that LigMatch outperforms several other widely used VS methods on the 13 DUD targets. We go on to demonstrate that improved VS performance can be gained from using multiple, structurally diverse templates rather than a single template ligand for a particular protein target. In this case, a 2D fingerprint-based method is used to select a ligand template from a set of known bioactive conformations. Furthermore, we show that LigMatch performs well even in the absence of 2D similarity to the template ligands, thereby demonstrating its robustness with respect to purely 2D methods and its potential for scaffold hopping. INTRODUCTION
Virtual screening (VS) is used to rank databases of chemical compounds according to decreasing probability of activity against a particular target. This prioritization means that only those molecules with significant probabilities of activity are screened experimentally, therefore decreasing the costs associated with the drug discovery process.1 With a number of success stories already reported,2 VS methods are increasingly playing an important role in drug discovery.3 When a crystal structure of the target protein is available, a structure-based approach will commonly be used.4 Structurebased methods often involve docking small molecules into the structure of the target protein and measuring the resulting interactions.3 However, it is well-known that current docking scoring methods are generally poor predictors of binding affinity,5 and even their ability to correctly rank candidates has been questioned.6 When the structure of the target protein is unavailable, a ligand-based approach must be followed. Ligand-based approaches require one or more known active molecules against which the database compounds can be compared.4 There are numerous examples in the literature of where ligand-based methods have been shown to outperform docking algorithms.7-10 Alternatively, when the structure of the target protein in complex with a ligand is available, the 3D bioactive conformation of the ligand will commonly be used as the ‘template’, i.e., reference molecule against which the database compounds are compared. Various different structure-based, ligand-centric methods often use the bioactive conformations of template ligands, including shape-matching algorithms,11 molecular field descriptors,12-14 and pharmacophore fingerprints.15-17 Alternatively, low-energy conformations,8 or even 2D representations,18-20 of the known * Corresponding author. Telephone: +44 (0)113 343 2592. Fax: +44 (0)113 343 3167. E-mail:
[email protected].
actives may be used as template ligands when no structural information is available. One of the main challenges when using any ligand-centric approach is to find the appropriate balance between search specificity, which is required for target selectivity, and search flexibility, which is necessary for removing any dependency on the template ligand.3 Indeed, ligand-based methods, particularly those based on simple molecular topology, are often criticized for finding close structural analogues instead of discovering novel structures, a concept known as ‘scaffold hopping’. Scaffold hopping is a highly desirable outcome of VS because it not only provides ways of improving the pharmacological properties of known lead compounds, but it also allows the exploration of unpatented regions of chemical space. VS methods are only truly useful for scaffold hopping if they are able to retrieve actives that are both leadlike and structurally novel, and a major consideration when demonstrating this in a retrospective validation is the data set of actives and decoys used. It is important that the decoys have similar physicochemical properties to the actives, otherwise the actives may simply be distinguished from the decoys based on trivial differences, such as molecular weight or logP.4 Huang et al.,21 recently constructed the directory of useful decoys (DUD) in which the decoys were specifically selected to be physicochemically similar yet topologically dissimilar to the actives, thereby providing a more realistic scenario for the validation of ligand-based VS. In this study, we use a version of the DUD that has been modified by Good and Oprea22 to ensure its suitability for benchmarking ligand-centric VS methods. The vast majority of similarity methods reported to date have made use of only a single bioactive template ligand. However, by focusing on only a single bioactive molecule, it is likely that actives with dissimilar structures will be overlooked, leading to a high false-negative rate.23 With the increasing availability of crystal structures of multiple,
10.1021/ci900204y CCC: $40.75 2009 American Chemical Society Published on Web 08/17/2009
STRUCTURE-BASED LIGAND MATCHING METHOD
structurally diverse, bioactive ligands for a given protein target, there has been recent interest in the use of multiple template ligands in VS.1 Indeed, fingerprint searching is often more effective when multiple template ligands are used due to an increase in the amount of information available for the calculations.24,25 Accordingly, several different approaches have been applied to multiple template-based fingerprint searching, including consensus fingerprints,26 centroid fingerprints,24,25,27 bit-scaling techniques28 and nearest-neighbor methods.24,25,28-31 A number of studies comparing these different search strategies have shown that the nearest-neighbor methods and the centroid methods often give the best performance.27,29 Consensus fingerprints, centroid fingerprints, and bit-scaling techniques emphasize bit positions that are shared by all template ligands and are, therefore, likely to account for their observed bioactivity.25 Although they are powerful methods, they may be of less use when applied to the structurally diverse molecules that are typical of high-throughput screening experiments.29 Indeed, they are more likely to rediscover the common features of known bioactive molecules rather than discover new chemotypes, which is a desirable feature of scaffoldhopping approaches. Conversely, the nearest-neighbor method individually calculates the similarity of each compound in the database to all of the template ligands, and then either the largest observed similarity score is selected to represent each database compound (MAX fusion rule) or the mean of the similarity scores is taken for the k nearest neighbors (SUM fusion rule).25 Both Hert et al.29 and Whittle et al.31 showed that the best recovery of active compounds was achieved when ranking the database molecules according to the MAX fusion rule. In this study we have extended this 2D approach by using the Tanimoto coefficient and the MAX fusion rule to select a template ligand against which each database compound can be compared in terms of its 3D structure, using a geometric hashing method.32 Here we report a new VS method called LigMatch, which uses a geometric hashing method32 in order to compare database compounds with the 3D bioactive conformation of a template ligand. We evaluate the performance of LigMatch on protein targets specifically selected from the DUD to allow a direct comparison with a recent study carried out by Cheeseright et al.4 For the 13 DUD single template targets, we show that the performance of LigMatch is better for six of them than for any of the four other popular VS methods, i.e., FieldScreen,4 DOCK,4,21 ROCS,33 and a 2D fingerprintbased method.34,35 Subsequently, we show that performance is highly dependent on the choice of ligand used as the template, and we demonstrate that the use of a 2D method to select the template ligand from a set of multiple template ligands can improve the performance of LigMatch. Finally, we systematically remove similarity between the template ligands and the database molecules and show that LigMatch is able to perform well even in the absence of 2D similarity to the template ligands, thereby demonstrating its potential for scaffold hopping. COMPUTATIONAL METHODS
Active and Decoy Molecules. The DUD (http://dud. docking.org/r2/, accessed 11/21/2008) provides a welldefined, unbiased set of active and decoy compounds that is
J. Chem. Inf. Model., Vol. 49, No. 9, 2009 2057
publicly available for the validation of docking algorithms. It contains sets of known actives for 40 target proteins, and sets of decoys that have been specifically selected to be physicochemically similar yet topologically dissimilar to these actives. Five molecular descriptors (molecular weight, logP, number of hydrogen-bond acceptors, number of hydrogen-bond donors, and number of rotatable bonds) were used to select decoys with similar physicochemical properties to the actives, while topological dissimilarities were calculated based on fingerprints. The strict criteria used in the decoy selection process mean that the DUD provides a stringent test with which to evaluate VS performance. As one of the best and most comprehensive collections of actives and decoys, the DUD is currently considered to be the industry standard for the validation of docking algorithms.23 However, since the DUD was designed as a benchmarking set for docking, it requires further processing to make it suitable for use in the validation of ligand-based methods.4 Indeed, in its raw state, there are many targets for which nearly all actives are trivial analogues of a central structure, causing it to be an inappropriate benchmark for 3D molecular similarity methods.36 In a recent study, Good and Oprea22 filtered the DUD actives for lead-likeness (molecular weight (Mw) < 450 and AlogP < 4.5) before clustering them according to chemotype using a reduced graph representation. In this way, they not only removed large molecules with unsuitable physicochemical properties, but they also addressed the problem of analogue bias. This filtered and clustered set of actives was used in the validation of FieldScreen, a field-based method reported recently by Cheeseright et al.4 In order to restrict their analyses to targets with a sufficient number of chemotypes, the authors selected a subset of 13 targets from the DUD, which each possessed at least 15 different clusters of actives. They then filtered the sets of decoys that corresponded to each of these targets using the same Mw and AlogP thresholds that were used by Good and Oprea when filtering the actives. The sets of active and decoy molecules corresponding to each of the 13 selected target proteins (ace, ache, cdk2, cox2, egfr, fxa, hivrt, inha, p38, pde5, pdgfrb, src, and vegfr2) were downloaded directly from the DUD (http://dud.docking.org/ r2/; DUD release 2, 10/22/2006) in mol2 file format. In order to allow a direct comparison with the results presented in the paper by Cheeseright et al., ranked lists of actives and decoys were obtained from the authors, and actives and decoys that were absent from their results were removed from our data set. It is worth noting that where multiple tautomeric forms of a molecule were present in the DUD, only a single tautomer was retained for our data set. Conformer Generation. Since our geometric hashing alignment and scoring method uses rigid conformations, multiple conformers of each active and decoy molecule were generated using Openeye’s OMEGA.37 OMEGA, which uses a deterministic rule-based method for conformer searching, was chosen because it has been shown to be very reliable and effective at reproducing bioactive conformations.38 An energy window of 25 kcal mol-1 was used since it was shown to be essential for best conformer generation performance in a study by Kirchmair et al.39 Indeed, a lower energy threshold could result in valuable conformations being discarded, whereas a higher energy threshold may produce high-energy conformations that are not likely to represent
2058
J. Chem. Inf. Model., Vol. 49, No. 9, 2009
KINNINGS
AND JACKSON
Table 1. Summary of the Data Sets Used in This Study
target angiotensin-converting enzyme (ace) acetylcholinesterase (ache) cyclin-dependent kinase 2 (cdk2) cyclooxygenase-2 (cox2) epidermal growth factor receptor (egfr) factor Xa (fxa) HIV reverse transcriptase (hivrt) enoyl ACP reductase (inha) P38 mitogen activated protein (p38) phosphodiesterase (pde5) platelet derived growth factor receptor kinase (pdgfrb) tyrosine kinase SRC (src) vascular endothelial growth factor receptor (vegfr2)
PDB code no. of actives (and ligand code) no. of possible (and active of single template useda conformers)c templatesb
no. of decoys (and decoy conformers)d
no. of clusters of activese (and no. of clusters represented by the possible templates)f
1o86 (LPR)
6
46 (2 310)
1,705 (119 737)
19 (3)
1eve (E20) 1ckp (PVB) 1cx2 (558) 1m17 (AQ4)
26 131 4 10
100 (5 162) 47 (1 525) 207 (1 675) 365 (10 980)
3,652 (166 556) 1,545 (54 966) 9,335 (336 371) 12,626 (526 281)
19 (4) 32 (14) 43 (4) 40 (1)
1f0r (815) 1rt1 (MKC) 1p44 (GEQ) 1kv2 (B96)
71 63 11 53
63 (5 298) 32 (605) 57 (1 696) 137 (2 545)
1,809 (143 852) 1,353 (37 853) 2,401 (84 357) 5,220 (141 564)
19 (6) 15 (9) 23 (1) 20 (4)
1xp0 (VDN) 1t46 (STI)g
8 n/a
26 (616) 120 (2 126)
1,504 (97 937) 4,529 (251 313)
22 (1) 22 (n/a)
2src (ANP) 1fgi (SU1)h
13 16
98 (1 395) 48 (1 631)
4,546 (312 785) 2,376 (132 311)
21 (1) 31 (0)
a The crystallographic ligand used as the single template in our study. It is identical to those used in the studies by Huang et al. and Cheeseright et al. b The total number of different ligands retrieved from all crystal structures of the target protein present in the RCSB PDB. These ligands were used for the multiple template method, and include the ligand used in the single template method (except in the case of vegfr2). c The number of active molecules taken from the DUD, and the number of conformers of these actives generated using OMEGA. d The number of decoy molecules taken from the DUD, and the number of conformers of these decoys generated using OMEGA. e Good and Oprea assigned each active to a cluster based on chemotype, therefore allowing an arithmetic weighting procedure to be used to remove analogue bias from the results. f The most similar active to each template ligand (greater than or equal to a Tanimoto coefficient threshold of 0.85) was selected, and the total number of clusters to which the selected actives belonged was then calculated. g Since no crystal structure of pdgfrb was available, the ligand coordinates for the template were taken from a homology model based on a c-Kit kinase-gleevec complex. h Ligand coordinates of SU5402 (ligand code; SU1) from the closely related fgfr1 protein kinase (PDB code; 1fgi) were used as the template.
bioactive conformations. A diversity threshold (root-meansquare deviation, rmsd) of 1 Å was also used to discard similar conformations, and a maximum of 100 conformers were generated per molecule. These parameters were used in the work by Sperandio et al.3 because they represented an optimal balance between speed and performance and agreed well with the study carried out by Kirchmair et al. For each target, a set of actives comprising of the original coordinates of the active molecules from the DUD, plus all generated conformers of these actives, was compiled. A set of decoys was compiled for each target in the same way. The numbers of actives and decoys for each target are given in Table 1. Selection of Single-Template Ligands. For each target protein, we used the coordinates of the ligand used in the studies by Huang et al.21 and subsequently by Cheeseright et al., which are derived from those in the RCSB Protein Data Bank (PDB)40 but with added hydrogen atoms. Since no crystal structure of pdgfrb was available, the ligand coordinates for the template were taken from a homology model based on a c-Kit kinase-gleevec complex (PDB code, 1t46; ligand code, STI). For vegfr2, the ligand coordinates of SU5402 (ligand code, SU1) from the closely related fgfr1 protein kinase (PDB code, 1fgi) were used as the template. Details of each of the 13 templates used (referred to as the ‘DUD ligands’) are given in Table 1. Scoring Using the Geometric Hashing Method. A geometric hashing method described previously32 was used to compare all actives and decoys for a particular protein target with the template ligand. The geometric hashing method proceeds as follows: A series of possible triplets consisting of three atoms forming a triangle is generated for both the database molecule and the template ligand. All possible
triplets between the two molecules are then compared, and those with the same atom types at triangle vertices (all atoms of the same element are treated as equivalent) and similar interatomic distances are defined as a match. A least-squares fitting routine is then used to perform a transformation for each match, resulting in a translational and rotational matrix that maps one 3D triplet onto another. Each match is then given a score corresponding to the number of coincident atoms in the mapping, and only the highest score, known as the ‘atom-atom score’, is retained. The atom-atom score was determined for every conformer of each active or decoy molecule, with both the best and the mean atom-atom score of all conformers being assigned to that particular molecule. The best and the mean atom-atom scores were then normalized to take into account the maximum possible number of matching atoms, by dividing by the number of atoms present in the active or decoy molecule. The resulting list of actives and decoys was subsequently ranked by both the descending normalized best score and descending normalized mean score. Since the normalized mean score gave slightly better performance overall (see the Supporting Information, Table SI 1), it was used in the studies reported herein. In terms of computational time, the geometric hashing method is able to process around 130 conformers per minute on a single processor, which equates to almost four million conformers in one day on a Linux cluster of 20 CPUs. Scoring Using Other Methods. In order to provide a comparison with LigMatch, four other methods were compared, using exactly the same sets of actives and decoys for each target. Both the FieldScreen and the DOCK rankings were obtained from Cheeseright et al.4,14 In addition, we ran ROCS and a 2D fingerprint-based method. ROCS is a shape-
STRUCTURE-BASED LIGAND MATCHING METHOD
based method that assesses the volume overlap of two molecules, using Gaussians parametrized according to the hard-sphere volume of heavy atoms.11 We used the ROCS combo score (default) to rank the lists of actives and decoys. The combo score is the sum of the shape Tanimoto and the scaled color score, where the former is a measure of the shape similarity and the latter is a measure of the chemical match between two molecules. Finally, in order to compare the performance of LigMatch with a 2D fingerprint-based method, Tanimoto coefficients were used to prioritize the actives and the decoys for each target. Open Babel34,35 was used to calculate the Tanimoto coefficients, using Daylight fingerprints,41 between the sets of actives and decoys and their corresponding templates. The actives and decoys for each target were then ranked by descending Tanimoto coefficient. Comparison of Methods on Single Templates. Since receiver operator characteristic (ROC) curve analysis is considered to be the best approach for characterizing the performance of VS methods,23 ROC curves were used to compare the overall performance of the five different methods. An arithmetic weighting procedure4,42 was applied in order to remove any analogue bias from the results. As mentioned previously, all active molecules for each target had been assigned by Good and Oprea to 15 or more different clusters, based on their chemotype (http://dud.docking.org/clusters/, accessed 01/14/2009). Each active molecule was given a weighting that was equal to the inverse of the size of the cluster to which it belonged. For instance, if five active molecules with similar chemotypes had been assigned to the same cluster, each active within that cluster would be given a weighting of one-fifth. Where there were actives that were absent from the Cheeseright et al. results (due to a failure to be processed by both FieldScreen and DOCK), the weightings of the remaining actives were modified accordingly. The number of clusters of actives for each target is given in Table 1. LigMatch Multiple-Template Method. A protein sequence corresponding to the relevant chain of each of 12 different target proteins was downloaded from the RCSB PDB. Since the X-ray crystal structure of pdgfrb had not been solved at the time of writing, no further analysis was conducted on this target. A protein-protein BLAST43 search with an e-value cutoff of e-40 and a sequence identity cutoff of 70% was carried out against the sequences of all proteins present in the RCSB PDB (http://www.rcsb.org/, accessed 09/16/2008 (kinases) and 02/10/2009 (nonkinases)). Protein structures with bound ligands were downloaded, and where there were multiple structures with the same ligand, only one representative structure was retained. Ligands were inspected to ensure that they were biologically relevant (resulting in the removal of small molecules, such as lipids and ions) and bound in the correct site, i.e., the same site as the DUD ligand. Ligands meeting these criteria were selected, so that for each target there was a set of different ligands with known crystallographic coordinates that would serve as templates against which actives and decoys could be scored. Note that the crystal structure of the DUD ligand was also included in each set of templates. The total number of templates for each target is given in Table 1 (for further details refer to the Supporting Information). Simply using
J. Chem. Inf. Model., Vol. 49, No. 9, 2009 2059
LigMatch to compare the active and decoy molecules to every individual template molecule for each target would be a timeconsuming process, especially for those targets with a large number of templates, e.g., cdk2, which has a total of 131 templates. Therefore, for each active or decoy molecule, a single template was selected from all of the templates available for a particular target by using Open Babel to calculate Tanimoto coefficients between that active or decoy molecule and each of the available template ligands. In accordance with the MAX fusion rule, the template with the highest Tanimoto coefficient was chosen as the template against which that active or decoy molecule and all of its conformers would be scored using the geometric hashing method. In order to quantify any overlap between the template ligands and the actives in the data set, the active with the highest Tanimoto coefficient, greater than or equal to a threshold of 0.85, to each template was selected. This threshold was chosen because it is generally accepted that molecules with Tanimoto coefficients of greater than 0.85 will exhibit similar biological activity.44-46 For each target, the total number of different clusters to which the selected actives belonged was then calculated. From Table 1 it can clearly be seen that the overlap between the templates and the actives is relatively low, i.e., the templates represent only a small number of clusters of actives in each case. Comparison of Single and Multiple Template Methods Using Ligmatch. A common measure of the performance of VS methods is the enrichment factor (EF), which measures the number of actives found within a compound database at a particular stage of VS. ROC enrichment (ROCE) differs slightly from the EF because it expresses the percentage of actives observed as a proportion of the percentage of the decoys observed rather than as a proportion of the percentage of the entire data set observed. Thus, ROCE does not show a dependence on the ratio of actives to decoys in the data set.47 Since early enrichment is an important measure of performance in VS, arithmetic weighted ROCE (awROCE) was calculated at the 1% level for all 12 targets in order to compare the LigMatch single template method with the LigMatch multiple template method. awROCE was calculated at the 1% level because it is likely that, in the drug discovery process, a similarly small proportion of the ranked database would be subject to further experimental screening. As with the EF, ROCE values of greater than 1.0 signify enrichment with respect to random. To provide a measure of VS performance across the entire data set, the arithmetic weighted area under the ROC curve (awAUC) was also calculated for all targets, using both the LigMatch single and multiple template methods. In contrast to the EF, the AUC shows no dependence on the ratio of actives to decoys in a data set. Ideal distributions of actives and decoys result in an AUC value approaching 1.0, whereas random distributions result in a value of 0.5. Performance of LigMatch in the Absence of Template Similarity. The method used to select the template ligand in the LigMatch multiple template method was modified in order to investigate the effects of reduced 2D similarity between the templates and the active and decoy molecules. Accordingly, the ligand with the highest Tanimoto coefficient to an active or decoy molecule, below or equal to a certain threshold, was chosen as the template for that particular molecule. Thresholds of 1.0 (i.e., no threshold),
2060
J. Chem. Inf. Model., Vol. 49, No. 9, 2009
0.8, 0.6, 0.4, and 0.2 were used, and where there were no templates with Tanimoto coefficients below or equal to a certain threshold, the relevant active or decoy molecules were given a normalized mean score and a weighting of zero. The mean enrichment across all 12 target proteins was calculated for each threshold, before being plotted as a ROC curve to give a generalized view of the overall performance. The above procedure was also performed using the 2D similarity method where the Tanimoto coefficient was used as the basis for prioritizing the actives and the decoys. In this way, it was possible to compare the performance of LigMatch with the 2D similarity method in an unbiased way at low template similarity. RESULTS AND DISCUSSION
Comparison of Methods on Single Templates. The performance of LigMatch on the benchmark of 13 selected target protein templates used in the study by Cheesewright et al. is shown by the awROC curves in Figure 1. The performance of FieldScreen, DOCK, ROCS, and the 2D fingerprint-based method on the same data set is also given for comparison. LigMatch gives the best performance out of all of the methods for six of the targets, in particular for ache, inha, p38, src, and vegfr2. For the targets where other methods outperform LigMatch, such as ace and egfr, it is never the worst performing method. For these two targets, ROCS and the 2D fingerprint-based method give superior performance. However, they appear unable to maintain this level of performance across all targets, particularly for fxa, p38, and vegfr2, where their performance is close to random. The performance of the methods in terms of enrichment (awROCE at 1, 2, and 5%) is given in Table 2. Although ROCS outperforms LigMatch at the 1% level for six of the 13 targets, it gives very poor performance for many of the other targets, in particular for src and vegfr2, where it fails to retrieve any actives at the 1% level. On average, across all 13 targets, LigMatch outperforms each of the other methods at all three enrichment levels. In addition, the standard deviations across the targets are lower than those of the 2D fingerprint method and ROCS (which are the two most competitive methods in terms of their performance). One target family that would appear to be particularly troublesome for many of the methods is the protein kinase family, which makes up six of the 13 targets, cdk2, egfr, p38, pdgfrb, src, and vegfr2. Interestingly, Huang et al. commented on the high failure rates observed for the protein kinase targets when using DOCK. Furthermore, in a recent study where they studied the performance of ROCS against all 40 DUD targets, Kirchmair et al.23 also showed that, when compared to other target families, the lowest performance, on average, was achieved for the protein kinases. A possible explanation of this failure could lie within the characteristics of the highly conserved protein kinase ATP binding site. Indeed, the ATP binding site is highly flexible and can alter its shape and size in response to ligand binding. Therefore, the use of a single ligand as a template ligand, or indeed a single protein structure as a template for docking, could be insufficient to represent the wide diversity of binding modes exhibited by the different ligands that bind a single protein kinase target (for a recent paper examining the structural diversity of protein kinase binding sites see Kinnings and
KINNINGS
AND JACKSON
Jackson (2009)48). With the exception of egfr, LigMatch gives the best performance overall for the protein kinase targets. For p38, src, and vegfr2 in particular, LigMatch performed well, while several of the other methods failed to perform better than random. Comparison of Single and Multiple Template Methods Using LigMatch. We have compared the LigMatch single template method (described above) with using multiple, structurally diverse templates for a target protein cocrystallized with different bioactive ligands. In this case, the 2D fingerprint was used to select a ligand template based on the Tanimoto coefficient. The effect of using multiple templates, rather than single templates, on the performance of LigMatch is shown in Figures 2 and 3. All awROCE values at 1% (Figure 2) are well above 1.0 (ranging from ∼10-70), and all awAUC values (Figure 3) are above 0.5, indicating that both single and multiple template methods perform much better than random across all targets. The multiple template method generally outperforms the single template method in terms of both early enrichment (awROCE at 1%) and enrichment across the entire data set (awAUC). As shown in Figure 2, the multiple template method improves the early enrichment for seven out of 12 targets. For one target, ace, there is no change, and for the remaining four targets, early enrichment is better using the single template method. In general, the greatest improvements in early enrichment from using the multiple template method can be observed where early enrichment is relatively low using the single template method. The awAUC values in Figure 3 show a very similar trend, suggesting that the results achieved for early enrichment are consistent with enrichment across the entire data set. For instance, pde5, p38, and cdk2, which have low early enrichment with the single template method, show the greatest percentage increases with the multiple template method (300, 185, and 79% improvement, respectively). Achieving good enrichment does not simply depend upon analogs of the template ligands being present in the data set of actives. Indeed, for pde5, the target showing the greatest improvement in enrichment, only one out of a total of 22 clusters of actives were represented by the eight different ligand templates (see Table 1). The targets that show a decrease in performance with the multiple template method already have the greatest early enrichment with the single template method, i.e. src, followed by inha and then vegfr2. In all cases except src, the early enrichment is never significantly reduced by moving to a multiple template model, and in many cases, it is significantly improved, especially for pde5, p38, fxa, hivrt, and cdk2. It would appear that the multiple template method removes some of the randomness from choosing an appropriate template, leading to a more robust method overall. For those targets that performed poorly using the single template method, but then showed a significant improvement upon using the multiple template method, the unsuitability of the DUD ligand may have been a contributing factor. For instance, p38 is a mitogen-activated protein kinase (PDB code, 1kv2), with the template ligand 1-(5-tert-butyl-2-ptolyl-2h-pyrazol-3-yl)- 3-[4-(2-morpholin-4-yl-ethoxy)-naphthalen-1-yl]-urea (ligand code, B96). The 185% increase in early enrichment when using the multiple template method on this target is possibly due to the unsuitability of B96 as
STRUCTURE-BASED LIGAND MATCHING METHOD
J. Chem. Inf. Model., Vol. 49, No. 9, 2009 2061
Figure 1. Arithmetic weighted ROC (awROC) curves to compare the performance of five different VS methods across the 13 selected DUD targets. The lines are colored as follows: LigMatch (green), FieldScreen (blue), DOCK (red), ROCS (purple), 2D fingerprint-based method (yellow), and random performance (black).
a template ligand. Indeed, B96 is a large ligand (MW > 500) that is bound in the rarely seen ‘DFG-out’ conformation of the protein.49 There are many known actives which receive low scores when using B96 in the single template method. For instance, the lowest ranked active (ZINC03815765), which was only retrieved after observing >99% of the data
set, is shown in Figure 4. The predicted binding pose is wrong with the molecule positioned in the back of the binding cleft and a major steric clash with a leucine residue (Leu167). However, when using the multiple template method, 4-(4-fluorophenyl)-1-(4-piperidinyl)-5-(2- amino-4pyrimidinyl)-imidazole (SB4) is selected as the template, and
2062
J. Chem. Inf. Model., Vol. 49, No. 9, 2009
KINNINGS
AND JACKSON
Table 2. Arithmetic Weighted ROC Enrichment (awROCE)a Values at 1, 2, and 5% for Five Different VS Methods Across the 13 Selected DUD Targets LigMatch
FieldScreen
DOCK
ROCS
2D fingerprint-based method
target
1%
2%
5%
1%
2%
5%
1%
2%
5%
1%
2%
5%
1%
2%
5%
ace ache cdk2 cox2 egfr fxa hivrt inha p38 pde5 pdgfrb src vegfr2
20.0 20.8 15.9 41.5 23.1 26.0 31.1 51.1 24.2 11.4 18.2 60.0 39.1
10.0 19.1 8.3 28.1 15.1 16.4 15.6 28.3 15.1 6.8 9.5 32.6 19.6
5.4 8.1 5.5 13.8 7.5 11.0 8.7 12.2 7.6 3.6 8.4 15.2 10.0
12.6 20.4 3.8 29.5 29.5 2.8 20.0 31.2 1.8 4.5 13.6 7.0 8.1
8.9 13.5 1.9 17.8 18.1 5.4 11.7 15.6 0.9 9.7 9.1 3.7 7.3
4.7 7.3 0.8 10.4 9.5 5.4 5.1 6.5 0.5 4.8 3.8 2.5 3.5
14.2 0.0 9.9 5.3 6.8 13.7 2.2 0.0 0.0 6.8 0.0 4.8 3.2
8.4 0.0 7.0 7.1 7.0 7.0 2.2 0.0 0.0 8.0 0.0 2.6 3.2
4.6 0.8 2.8 5.5 4.5 6.2 2.2 0.0 0.4 4.3 0.0 1.0 1.3
44.5 23.7 20.3 48.1 77.6 10.5 20.0 25.6 4.9 4.5 22.7 0.0 0.0
22.9 15.4 10.2 29.3 39.6 5.3 10.0 12.8 5.4 2.8 11.4 0.0 0.8
10.2 7.2 4.4 12.8 17.1 2.1 4.0 5.4 3.2 1.4 4.5 0.0 0.3
48.2 11.9 15.9 11.0 71.7 5.3 26.7 39.1 2.1 13.6 29.5 0.4 9.7
29.6 6.0 8.8 6.1 41.0 5.3 13.3 19.6 1.1 6.8 16.2 4.9 6.5
14.0 2.9 5.2 4.3 17.5 2.1 10.2 9.6 0.4 3.6 7.1 8.0 4.0
mean sd
29.4 14.5
17.3 8.2
9.0 3.4
14.2 10.8
9.5 5.7
5.0 3.0
5.1 5.0
4.0 3.4
2.6 2.2
23.3 22.3
12.8 11.7
5.6 5.1
21.9 20.7
12.7 11.5
6.8 4.9
a
Equivalent awROC values are provided in the Supporting Information.
Figure 2. Arithmetic weighted ROC enrichment (awROCE) at 1% for both the single (light gray) and the multiple template methods (dark gray) using LigMatch for each of the 12 DUD targets.
Figure 4. The best scoring conformation of a p38 active molecule (ZINC03815765, colored by element) superimposed onto the single template ligand B96 (shown in yellow). The protein structure (PDB code, 1kv2) is shown as a ribbon and colored by secondary structure (Leu167 is shown in green stick representation). N.B. the protein structure was not considered in this study.
Figure 3. Arithmetic weighted area under the curve (awAUC) for both the single (light gray) and the multiple template methods (dark gray) using LigMatch for each of the twelve DUD targets.
the resulting binding pose is much more realistic (Figure 5), causing the active to be given a much higher score and to be retrieved after observing
