The automated comparison of protein-ligand binding sites provides useful insights into yet unexplored site similarities. Various stages of computational and chemical biology research can benefit from this knowledge. The search for putative off-targets and the establishment of polypharmacological effects by comparing binding sites led to promising results for numerous projects. Although many cavity comparison methods are available, a comprehensive analysis to guide the choice of a tool for a specific application is wanting. Moreover, the broad variety of binding site modeling approaches, comparison algorithms, and scoring metrics impedes this choice. Herein, we aim to elucidate strengths and weaknesses of binding site comparison methodologies. A detailed benchmark study is the only possibility to rationalize the selection of appropriate tools for different scenarios. Specific evaluation data sets were developed to shed light on multiple aspects of binding site comparison. An assembly of all applied benchmark sets (ProSPECCTs–Protein Site Pairs for the Evaluation of Cavity Comparison Tools) is made available for the evaluation and optimization of further and still emerging methods. The results indicate the importance of such analyses to facilitate the choice of a methodology that complies with the requirements of a specific scientific challenge.
Due to the huge number of available binding site comparison algorithms[5], a large-scale analysis of all methods is beyond the scope of this article. We have therefore restricted the evaluation to a small, but still diverse subset of promising algorithms. These were derived from an analysis of successful applications within medicinal chemistry projects[5]. A comparison of web server-based tools became infeasible due to the number and size of analyzed data sets, so we restricted our evaluation to standalone tools. The methods analyzed herein and their fields of utilization are summarized in Table 1. Intriguingly, the success of nearly all of those studies resulted from the use of binding site comparison as part of a workflow combining different computational methods, e. g. MD simulations or molecular docking studies. A recent, impressive example shows that the combination of various tools of structure-based modeling and binding site comparison delivers insight into putative mechanisms of drug action[17]. Herein, we analyzed binding site comparison methods which are based on fingerprints (e. g. , SiteAlign[18] and TIFP (Interaction Fingerprint Triplets) [19]), graphs (e. g. , Cavbase[20,21] and IsoMIF[22]), grids (VolSite/Shaper[23]), and those that make use of alternative approaches (PocketMatch[24] and SiteHopper[25]). KRIPO (Key Representation of Interaction in Pockets), that was originally designed for sub-pocket matching to facilitate bioisosteric replacement, can also be applied for ligand binding site comparison[26] and was included in our analysis. Additionally, TM-align[27], which was developed to compare protein structures based on their overall structure, was evaluated as it was successfully applied in various medicinal chemistry scenarios. The freely available tools CMASA[28], COFACTOR[29], and PocketFEATURE[30] could not be analyzed. Both CMASA and COFACTOR can only be applied with pre-prepared data sets; they enable the user to compare the binding site of interest against sets of precomputed binding sites. PocketFEATURE is not publicly available. Grim (Graph Interaction Matching) [19] and RAPMAD (Rapid Pocket Matching using Distances) [31] were included as they make use of binding site representations highly similar to TIFP and Cavbase. The impact of the underlying data structures was evaluated (fingerprints vs. graph models for TIFP and Grim, histograms vs. graph models for RAPMAD and Cavbase, respectively). While fingerprint- and histogram-based methods are usually characterized by low computational demands, graph models enable a more accurate binding site characterization accompanied by higher run times. The representation of the binding site for the subsequent comparison algorithm has a large impact on the outcome of the investigations. We therefore classified all approaches as depicted in Fig 1. While most tools encode the binding site features based on the underlying ligand-interacting residues, other approaches use surface-based binding site representations, e. g. through a projection of physicochemical site features onto the respective surface patches. Probably the highest level of abstraction is achieved by programs that encode protein-ligand interactions, i. e. they do not depend on a distinct set of residue classes or functional groups, but on different types of interactions. While the residue type is crucial for most residue- and surface-based methods, interaction-based methods often rely on the nature of the bound ligand. However, this classification only holds true for binding site modeling. The final scoring scheme applied to a binding site pair match might also include other properties, for example, surface similarity. Moreover, the data structure for comparison differs for methods within a category. While residue-based comparisons are achieved using versatile approaches, trends can be derived for surface- and interaction-based methods. Interaction-based methods make use of graph and fingerprint representations. Graph and 3D point approaches are often used for surface-based comparisons. Brief descriptions of the methods and explanations of the scoring schemes are given in the SI (S1 Text) to outline the general concept underlying each method. This information on different scoring measures and comparison approaches is included to support the understanding of the outcomes presented herein. The explanations are by no means complete, but are focused on the most essential points. Although these descriptions are sufficient for the scope of this study, we encourage readers to refer to the methods’ publications to gain further insight into usage and parameter details. In most cases, default settings were applied in our work, hence parameter optimization might lead to improved performance.
As for many other computational methods, the success of binding site comparison methods has to be evaluated using benchmark data sets. Some evaluation sets previously created to test the tools analyzed herein are summarized in the SI (S1 Table). Although we do not claim that our collection is comprehensive, there are no commonly used, state-of-the-art evaluation data sets available. Standardized benchmark sets are accessible for broadly applied modeling approaches such as pharmacophore searches[61] and molecular docking[16,62]. In contrast, the high diversity of the applied benchmark sets for binding site comparison makes it difficult to draw definitive conclusions in comparing the different tools. There is a need for a common evaluation scheme to assess the applicability of highly diverse methodologies. Even knowledge of the underlying concept is not sufficient to make a confident choice of an appropriate tool. For example, Cavbase and SiteAlign make use of similar binding site representations, but the underlying comparison algorithms and scoring schemes differ. It is not obvious whether a certain method should be used preferentially. The evaluation of SiteAlign (S1 Table) would suggest that the method is able to detect evolutionary relationships by identifying similarities between binding sites with similar biological functions. However, it has been successfully applied for the investigation of protein-ligand interactions (Table 1). TM-align was never evaluated as a cavity comparison tool, but it has given hints on interesting binding site relationships (Table 1). The high number of application domains encouraged us to develop novel, objective data sets with respect to specific aspects of binding site comparison. Pairs of similar and dissimilar binding sites enable an objective and detailed analysis of available tools. Various benchmark studies were performed for the methods mentioned above. Based on the results, we attempted to discern whether a suitable method, both in general and for specific application domains, can be selected. The first data set (structures with identical sequences) was designed to evaluate the tools’ sensitivity to the binding site definition. This definition often depends on the size and location of bound ligands. Different ligands can address various regions of the binding site of interest (sub-pockets). Because of these different site definitions, similar sub-pockets are more difficult to match. Although binding sites can interact with a broad variety of ligands, they share common properties and distinct similarities. Thus, the scoring scheme has to be optimized for this scenario to enrich similar binding site pairs with different cavity definitions in a ranked list. Data set 1 contains structures with identical sequences, which bind to chemically different ligands located at identical sites, leading to diverse binding site definitions used for the comparisons. The second data set assesses the tools’ performance with respect to the binding site flexibility. This is an important factor when comparing two similar pockets. The data set is based on protein models extracted from ligand-bound solution NMR structures with more than one model in the structure ensemble. The next two data sets were used to elucidate the scoring metrics’ discrimination performance on nearly identical binding sites differing by a single substitution or by multiple mutations. Artificial protein structures were created by randomly picking binding site residues and substituting them with physicochemically different residues (decoy set 1) and residues that lead to a change in the binding sites’ size and physicochemical properties (decoy set 2). A comparison of binding site pairs of proteins on a data set containing both the original sequences and their artificially generated counterparts should lead to an enrichment of binding site pairs of proteins with the original sequences. Pairs of original and modified binding sites should obtain lower similarity scores. Two pre-existing data sets were used to evaluate whether tools are able to differentiate between binding sites that bind different ligands and to identify similarities between binding sites occupied by identical or highly similar ligands. One such data set was described by Kahraman et al. [63]. This was originally designed to evaluate whether binding site shape and ligand shape are related. The authors state that the variability of binding site shape cannot solely be explained by the conformational variability of the ligand. Although the data set structures are derived from unrelated proteins, it was not investigated whether there are local similarities between the binding sites with similar ligands. In contrast, the data set of Barelier et al. [64] comprises 62 pairs of unrelated proteins binding to similar ligands. It includes 19 pairs of binding sites that show local similarities (as “observed” by the ligand atoms) whereas the remaining pairs do not display any obvious resemblance. Finally, we established a data set of binding site pairs whose similarity was correctly identified by at least one binding site comparison tool as described in a recent perspective[5]. We combined those similar binding site pairs with a diverse data set of sc-PDB[65] derived binding sites. A comparison of the query structures (binding sites with known similarities) to the complete data set (data set of successful applications) can be performed. This analysis allows the assessment of whether the tools are able to enrich similar site pairs within the high-scoring hits. This benchmark set contains the most interesting pairs as their similarity proved useful in various medicinal chemistry projects. Table 2 gives an overview of the benchmark data sets, their main goals, composition, and experimental quality (if applicable). For some data sets, we generated two versions (e. g. data set 1 and data set 1. 2) to focus on specific aspects of binding site comparison. Details regarding the individual data sets can be found in the Methods section and S1 Fig and S2, S3 and S5–S11 Tables. As all data sets include active (i. e. similar) and inactive (i. e. dissimilar) pairs, the ROC (receiver operating characteristics) curves as well as EFs (enrichment factors) at different percentages of the screened data set for all analyses can be calculated. A detailed summary of all AUC (area under the ROC curve) values and EFs obtained for the applied tools is given in S14, S16, S18, S20, S22, S24, S26, S28, S30 and S32 Tables. The significance of the AUC values and the differences between the methods’ AUC values for the different data sets are provided in S15, S17, S19, S21, S23, S25, S27, S29, S31 and S33 Tables. A final aspect for choosing an appropriate tool is the computational cost required by the different methods. The interplay of cavity preparation, number of modeled binding site properties, the implementation of the comparison algorithm, and filtering steps determines the final CPU time. Therefore, the run time of binding site preparation as well as comparison per method was analyzed. Table 6 summarizes the different algorithms with respect to the run time for binding site pre-processing of all 100 structures of data set 5 and for 10,000 comparisons (all-against-all). As PocketMatch, SiteAlign, and TM-align rely on user-defined cavities, which were prepared in an automated manner using a Python script, i. e. exemplary run times are given that might change for a different site extraction method. Interestingly, the tools were not as robust as expected and failed at different steps of preparation and comparison. The failure rates of all methods for the analyzed data sets are given in S10 Fig. The most interesting finding was the importance of PDB to MOL2 conversion to decrease the failure rates of Grim and TIFP. Both methods should be preferentially applied with MOL2 files of the protein-ligand complexes. The calculation of the average run times per comparison takes the number of omitted comparisons into consideration. For Cavbase and RAPMAD, all structures could be processed, but the ligand-defined binding site was not identified for nine structures. Therefore, the comparison run time is given for 8,281 comparisons. For FuzCav and Shaper, only 96 out of 100 binding sites could be prepared. A pre-processing of the prepared pockets (with the PDB file parser pdbconv of IChem) with VolSite for a subsequent Shaper comparison (VolSite/Shaper) led to only 57 pockets derived from PDB files. In contrast, 76 cavities were extracted with MOL2 files as input for pdbconv. TIFP calculations led to 47 fingerprints for PDB files and 77 for MOL2 files. In Table 6, the numbers in brackets summarize the numbers of prepared structures and comparisons for each tool. A clear correlation between the comparison method used and the run time could not be observed. Nevertheless, it is possible to differentiate between very fast methods (several microseconds per comparison), moderately fast methods (several milliseconds per comparison), and comparably slow methods (several seconds per comparison). Depending on the desired outcome and the size of the data set, the computational cost might be a limiting factor and the use of some methods becomes infeasible. For a comparison of a single binding site of interest against all known pockets as stored in the sc-PDB[65] (9,283 entries) on a single CPU, three days will be necessary using SiteAlign while PocketMatch will perform these comparisons within 0. 28 seconds. For an all-against-all comparison, RAPMAD is approximately 9,200 times faster than the graph-based method Cavbase. In contrast, TIFP is only twice as fast as the graph-based method Grim. Such minor differences suggest the use of the better-performing method. The run time also depends on the type of comparisons. Pairwise comparisons were performed for Grim, IsoMIF, VolSite/Shaper, SiteAlign, SMAP, TIFP, and TM-align, i. e. the tool was invoked for each comparison separately. This might become the time-limiting step for some methods. Cavbase, FuzCav, PocketMatch, ProBiS, RAPMAD, SiteEngine, and SiteHopper allow the comparison of one query against a list of targets in one run. This causes a speedup as the tool has to be invoked only once. For ProBiS, the conversion of PDB files to a surface file format is necessary. Thus the cavity preparation time reflects this preliminary step. Additionally, KRIPO and RAPMAD enable even faster comparisons by providing SQLite databases of the modeled sites and the automated calculation of similarity matrices. These factors illustrate the difficulty in analyzing and comparing run times. In this analysis, we chose the fastest option available for each tool. The computational cost will therefore differ depending on the tool’s implementation and the type of approach employed (comparing distinct binding site pairs, query-based analyses, generation of distance matrices, etc.). Finally, the nature of the data set has to be taken into consideration. Changing the binding site size might have a huge impact on the run time per comparison. While the run times of fingerprint-based methods are often not affected by this factor, graph-based methods might become significantly slower with an increasing graph size. Cavbase, for example, fails to compare very large binding sites. ProBiS uses pre-filtering steps to speed up comparisons.
Plenty of factors influence the decision for a suitable binding site comparison tool. S13 Table summarizes the tools’ characteristics with respect to the most important criteria and provides a detailed assessment of the tools’ performance. First of all, the necessary pre-processing of structures might have an impact on the choice. The preparation of large data sets, e. g. for the prediction of potential off-targets using all known binding sites, requires tools that enable an automated and flawless binding site processing and annotation. We divided this criterion into two parts. Firstly, it is important that binding sites can be prepared in an automated fashion. For example, the KRIPO developers generated scripts to automatically process the binding sites of interest. The same holds true for FuzCav, Grim, TIFP, Shaper, and VolSite/Shaper (although for some PDB structures ligand MOL2 templates have to be provided). Cavbase and RAPMAD require XML-formatted cavity descriptions as input which can be generated with the help of the CSD Python API[90] in an automated fashion. Secondly, the inclusion of all available binding sites might be crucial (S10 Fig). The prediction of potential off-targets and identification of novel binding sites require comprehensive cavity data sets. It is also vital that the chosen tools are reasonably fast with respect to such applications (run time). Sometimes, it is of interest to investigate the similarity between predicted binding sites and already known ones, e. g. for the identification of druggable binding sites of novel targets. Many interaction-based methods will not be able to compare predicted pockets (Grim, KRIPO, TIFP). SiteEngine will also fail as the tool relies on protein- and ligand-surface construction to compare the resulting sites’ surfaces. Furthermore, SMAP does not provide the possibility to process detected binding sites. The only possibility is the generation of putative ligand binding modes via docking or pharmacophore searches to obtain structures with bound ligands for these methods. The tools Cavbase, FuzCav, RAPMAD, Shaper, and IsoMIF come along with implemented binding site identification approaches. Additionally, PocketMatch, ProBiS, SiteAlign, and TM-align process user-defined binding sites based on residues. Therefore, it is possible to apply them to externally predicted binding sites. SiteHopper was shown to compare predicted pockets as long as fpocket[91]-derived cavities are utilized as a pseudo-ligand[53]. Based on the analyzed data sets we can also draw conclusions with respect to a suitable input for binding site comparison. Tools that showed a poor performance for data set 1 should not be used if only a small number of protein-ligand complex structures are known. These tools are very sensitive toward the nature of the ligand. The problem may possibly be circumvented by using docking-derived poses of ligands in the binding site of interest. Some methods’ scoring schemes suffer from a small similarity score difference window for similar and dissimilar binding sites. In many cases, it is advisable to use all available crystal structures or multiple NMR models to retrieve promising results. Although the accuracy of NMR structures is poorly validated[92], they seem to be an acceptable choice for different structural biology and computational modeling challenges[93,94]. The G-factors calculated with PROCHECK-NMR (Table 2) indicate a lower quality as compared to the X-ray data sets. Nevertheless, the average of -0. 28 point toward usual dihedral angles and main-chain covalent forces. Several tools were developed to facilitate the selection of a diverse set of representative protein structures[95,96]. In the absence of experimental structural ensembles, the use of MD-derived binding site conformations as input for a comparison might be beneficial. With respect to the decoy data sets, there are only a few tools that are able to distinguish between similar pockets with minor dissimilarities and highly dissimilar sites. SMAP, SiteHopper, and KRIPO are best suited to convincingly penalize differences between protein structures and their artificially created decoy variants with different numbers of substitutions. Other tools failed to score dissimilarities appropriately. This can be partially attributed to the fact that the binding site’s shape was retained in data set 3; yet, the introduction of major shape differences (data set 4) did not cause substantial improvements. Tools which do not reliably score minor differences in structures might provide better results if they are applied to elucidate similarities between unrelated cavities. In such cases, residue changes and geometric deviations have to be tolerated to identify such partial similarities within the best-ranked binding site pairs. The evaluation of the tools with regard to the data set of Barelier et al. [64] shows that a final visual inspection is unavoidable if the user is interested in a drug repurposing strategy or the establishment of polypharmacology. Therefore, we have to differentiate between the methods depending on the possibility to visually analyze the resulting binding site alignment. The data set of similar binding site pairs extracted from lifelike analyses gives some insight into the basic applicability of the tools for different scenarios. The observation that all tools showed a very high early enrichment underlines the fact that all tools perform well for relevant examples. The individual strengths and weaknesses of all tools finally level out to considerably high early enrichment of similar binding site pairs for all methods. Whenever possible, it is advisable to use more than one tool as they might complement one another. The comparison methods Grim and TIFP are characterized by a comparably low applicability toward the identification of binding site similarities. The original application domain both tools were designed for was the rescoring of docking poses and pose selection, e. g. in a structure-based virtual screening campaign. For these approaches they perform well as shown in the original publication[19] and quite recently for the D3R grand challenge 2015[97]. We could not use the pre-processed structures as stored in the sc-PDB[65] as our data sets make use of all available PDB structures. Restricting the binding site comparison with TIFP and Grim to sc-PDB derived structures might provide a different perspective for both tools as their performance depends highly on the pre-processing of the binding sites. Both tools assign interaction patterns depending on the nature of the ligand which impairs their success for some of our chosen data sets. Nevertheless, their performance on the data set of lifelike examples was good. Finally, the application of score thresholds to discriminate between similar and dissimilar binding sites can be discussed. Similarity scoring for binding sites is challenging. Often, score cut-off values have to be optimized for specific studies[11]. Our investigations also show the strong context dependency of such thresholds (S12 Table). If a clear classification of similar and dissimilar pairs is necessary, the appropriate cut-off has to be evaluated with a data set that reflects the main purpose of the analysis. Score thresholds for the decoy data sets (data set 3 and 4) should be taken into account if binding sites with minor dissimilarities are compared. In contrast to the elucidation of local binding site similarities or similarities of pockets of unrelated proteins, the less restrictive cut-offs determined for the data set of Kahraman and colleagues[63] or data set 7 are more appropriate. The methods SMAP and ProBiS offer additional statistical measures to estimate the significance of a binding site match. Nevertheless, we want to emphasize that the application of thresholds often means disregarding potentially interesting and yet unexplored binding site similarities. The Venn diagram in Fig 15A summarizes six potentially meaningful categories which influence the decision on a useful tool. Depending on the problem to be addressed, other criteria out of a multitude as given in S13 Table might also facilitate the choice. We strongly recommend this comprehensive table for further details. Nevertheless, we wanted to give an example for guiding the way to the most appropriate tool according to the criteria selected. FuzCav seems to fulfil all criteria in the Venn diagram and can be safely applied for evolutionary analyses. However, the binding site matches cannot be visualized. Hence, the method is not suited for query-based drug repurposing projects which require a detailed examination of the identified similarities. The use of KRIPO enables the generation of site alignments, but the method is not applicable for predicted sites and consequently not suited for function prediction. ProBiS failed to rank the similarity according to the number of residue substitutions and will probably fail to accurately score minor cavity dissimilarities. As this is decisive in the analysis of evolutionary binding site relationships, such scenarios should be analyzed with SiteAlign, SiteEngine, or SMAP. A potential user can thereby iteratively exclude tools to arrive at a final choice. Fig 15B puts all results into context and can be used to discuss possible application areas for the binding site analysis tools. The applicability suggestions for the different types of tools are based on previous studies illustrating different uses in drug discovery[5]. However, the diagram also hints at possible fields of applications yet unexploited, for example, the use of interaction-based methods for selectivity prediction or surface-based methods for target elucidation. Of course, the general aim of a study, the basis of data, and the number of comparisons necessary will influence the choice of a method. The embedding of a tool in an elaborate in silico workflow might serve to diminish certain weaknesses. Methods never applied to date for a certain field of research might provide useful results within specialized projects. Some obvious new applications are the use of interaction-based methods for the analysis of protein-protein and protein-ligand interactions or drug repurposing approaches exploiting the outcome of surface-based site comparisons. The prediction of off-targets with comparably fast residue-based approaches is another unexplored application. Furthermore, some tools have never been applied in independent studies. Several projects which made use of Grim[97], RAPMAD[98], and SiteHopper[99] demonstrate their unique capabilities. Their strengths as highlighted herein might encourage researchers to actively use these tools. Although it is not possible to identify any single tool which fits all needs and performs well for all data sets, we can provide some guidance regarding different aspects of binding site comparison. The impact of the ease and completeness of binding site pre-processing steps increases with an increasing number of proteins used for comparison. While both factors are crucial for elaborate projects, e. g. off-target prediction or the identification of novel binding sites, their importance decreases for relatively small data sets used in selectivity profiling or polypharmacology prediction. An application to predicted binding sites is relevant for off-target prediction, drug repurposing, or the identification of potential targets. For polypharmacology elucidation, one generally refers to known druggable binding sites. The analysis of similar binding sites to elucidate minor structural differences is highly influenced by the binding site’s definition and flexibility, as is the ranking with respect to binding site properties. In contrast, the elucidation of similarities between the cavities of unrelated proteins does not necessarily depend on accurate scoring. MD simulations or the use of NMR ensembles might help to circumvent a potential failure of binding site comparison due to insufficient consideration of protein flexibility. It is necessary to unravel the similarity between functionally unrelated binding sites which bind similar ligands when understanding polypharmacology. Visualization of identified similarities, which is essential when dealing with non-obvious similarities, is not necessarily crucial for selectivity profiling or analyzing evolutionary relationships. These criteria should become the focus of future benchmark analyses of other promising site comparison approaches as such analyses can guide the rational choice of a method. Ultimately, the choice of a comparison method depends on the focus of the study. These investigations can help to ease the choice of a suitable tool, though restricted to a limited subset of available approaches. The publication of the generated data sets and the benchmark results can assist in the assessment of tools and the establishment of reliable workflows that consider individual strengths and weaknesses. We hope that the assembly of benchmark sets (ProSPECCTs) and the conclusions drawn from the evaluation encourage researchers to objectively assess the advantages and drawbacks of individual approaches. Finally, this guide could facilitate the final choice of a suitable method and enable researchers to derive an advantage from these–as far as our experience goes—widely underemployed binding site comparison approaches.
The structures of all data sets were prepared in the same way to ensure an equivalent basis for all binding site comparison tools. First, modified residues were identified and the respective HETATM record names in the PDB files were changed to ATOM. This modification did not affect the tools’ performance, but was essential for the site processing with SiteEngine which detects ligands based on the HETATM record names. The binding site’s defining ligand was identified and renamed to LIG for further processing steps. Other HETATM entries were deleted to ensure the exclusion of buffer ions, cofactors, and prosthetic groups. The final steps were realized with the help of the pdbcur tool of the CCP4[101] software package (version 6. 5). Alternative locations with the highest occupancy were retained, while for locations with identical occupancy values the first one was retained. Finally, ANISOU entries were removed. The resolution and R-factors (R work) for all structures were downloaded as a report from the PDB[1]. The mean values, standard deviations, minima, and maxima of these parameters were calculated for all data sets including X-ray structures. G-factors were calculated using PROCHECK[66] and PROCHECK-NMR[67]. These values measure the degree of unusual dihedral angles and covalent forces of the main chain. G-factors below -0. 5 hint at unusual structure properties. For the groups of structures with identical sequences and the NMR ensembles, all residues were renumbered according to the sequence alignment calculated with default settings in MOE2013[102]. Binding site residues were assigned based on a 5 Å radius of all ligand atoms. Subsequently, the Cα atoms and all atoms of the binding site-defining residues were aligned using the “match” command of UCSF Chimera[60]. The mean RMSD values, standard deviations, minima and maxima of all pairwise comparison were calculated to characterize the binding site flexibility. For the data set of identical structures, Tanimoto coefficients based on the ECFP4 fingerprints were calculated in a pairwise manner for all groups using KNIME[103]. Binding site descriptors were calculated using DoGSite. The binding site ligand was chosen as the reference ligand and the pocket was defined by the ligand. Apart from these changes, default settings were used. The resulting pocket descriptors were analyzed for each NMR ensemble of data set 2, the sequence-culled set of the sc-PDB and the structures of the data set of Barelier and co-workers[64]. Unless stated differently, default parameters of the analyzed tools were used. The exception that holds true for all methods is that the scoring measure used for ranking was selected based on the early enrichment for the data sets of structures with identical sequences, NMR, and decoy structures with descending priority (Table 9). Regardless, all scores were calculated for further evaluation. Additional information for the different scoring schemes used is given in the SI (S1 Text). Most binding site comparison tools applied in this study have some major limitations that have to be taken into account when applying them to different types of data sets. They will be discussed in the following. PocketMatch (version 2. 0), SiteAlign (version 4. 0), and TM-align (version 20170708) rely on pre-processed ligand binding sites. The respective binding sites were extracted by means of a Python script which creates an output PDB file that contains all protein atoms within a given radius of the ligand atoms. This excised binding site was used for the comparison. A 5 Å radius was applied for the creation of input cavities for PocketMatch and SiteAlign. Residues within 10 Å of the ligand were used as the binding site definition for TM-align to guarantee a sufficient number of residues to yield a reliable alignment. The programs Cavbase and RAPMAD require XML-formatted pockets as input for binding site comparison. The extraction of the pockets was achieved using the CSD Python API 1. 3[90] from the CCDC. Residues with missing backbone atoms which were not part of the binding site were excluded to obtain pre-processed cavities as these residues were not properly processed. The cavity extraction is based on pockets detected by LIGSITE[76]. Thus, the binding site of interest could not be found for some PDB entries. We used only cavities including ligand atoms whereas other detected pockets were excluded from the analysis. If more than one pocket for the ligand of interest was identified, the pocket with the largest cavity volume was used. Some very large cavities could not be processed using Cavbase. The use of RAPMAD is restricted to similarity scoring of binding sites as it does not generate a binding site alignment. FuzCav comparisons can be performed for binding sites extracted with the pdbconv tool of IChem, i. e. the complete site or residues within 4,6, 8,10, or 12 Å of the ligand can be used for the comparison. A 6Å radius was used in our analyses. IsoMIF is a molecular interaction-based method that generally relies on bound ligands. A cut-off radius can be used to restrict the binding site dimensions depending on the ligand. The use of the additional tool GetCleft[108] allows for the inclusion of predicted cavities. KRIPO fingerprint databases can be individually created by users for their own data sets. A ligand database is provided for all PDB entries. For proprietary or modified structures, a ligand database has to be generated. Fingerprints can be prepared for ligand fragments as well as for complete ligands with the command line tool “kripo”. Subsequently, the KRIPO DB package (“kripodb”) was used to perform all-against-all comparisons for the derived fingerprint databases for complete ligands. The final similarity matrix did not include similarities between identical structures. The identity scores were automatically set to 1. PocketMatch does not take modified residues into account. The tool is suitable for similarity ranking but the alignments are not available as output. The PDB convention of 80 characters per line has to be fulfilled for all input pockets. Binding sites are defined by the ligand structures for ProBiS comparisons. A distance threshold can be applied to modify the cavity definition. By default, ProBiS does not provide scores for all binding site pairs of interest. Cliques with poor scores and/or low z-scores are deleted. The “noprune” and “z-score” options offer the output of insignificant matches together with the significant ones. The use of Shaper relies on cavity definitions by VolSite. The program predicts the druggability of binding sites and excludes those that are denoted non-druggable. Not all cavities can therefore be processed with Shaper. We additionally applied Shaper to binding sites extracted by the pdbconv tool alone to allow for a more complete processing of binding sites although this is not recommended by the developers. The complete detected cavity as well as cavities defined within a 4,6, 8,10, or 12 Å radius of the ligand can be used. According to the recommendations, we used a 6 Å radius. Modified residues cannot be processed by SiteAlign. A binding site definition based on residue names and numbers is necessary. Insertion numbers are not supported requiring a preliminary renumbering of residues with insertion codes in the PDB files. The comparison of one query against a list of targets frequently fails with “Segmentation Error”. It is unavoidable to do comparisons of one query against one target each to avoid high failure rates. Alternate atom locations are not supported by SiteEngine. A single conformation has to be retained in the PDB file. Moreover, the HETATM entries for modified residues have to be changed to ATOM entries for correct surface construction. The tool is highly sensitive toward PDB files with more than 80 characters per line. Some large protein structures were not properly processed and had to be excluded. SiteEngine comparisons depend on ligand-defined cavities and its use is restricted to protein binding sites with bound ligands. The distance cut-off can be manually adjusted to define the binding site. SiteHopper initially creates binding site patches. This step failed for some protein-ligand complexes. Furthermore, the tool was not able to process residues with missing backbone atoms. These residues were therefore excluded for the comparisons after ensuring that they are not part of the binding site of interest. SiteHopper relies on ligand defined binding sites, but the cut-off radius can be adjusted. For small ligands, an additional flag had to be set during protein-ligand splitting (-min_atoms). This flag was set to 0 to process phosphate binding sites and to 1 for the remaining data sets. SMAP (version 2. 0) comparisons worked for nearly all structures after providing the PDB files in an appropriate manner. In some cases, the tool was highly sensitive with respect to the provided structures (including REMARK lines). The binding site definition can be modified using a distance threshold. The tools TIFP and Grim (as tools of IChem version 5. 2. 6) are based on interaction fingerprints. Only ligand-occupied cavities can therefore be compared. The tool pdbconv of the IChem toolbox must be applied to extract proper binding sites and interaction fingerprints. Problems also arose when the ligand of interest giben in the PDB format was not in the predefined template files. Input structure preparation therefore included a PDB to MOL2 conversion with the CSD Python API 1. 3[90] to obtain more reliable fingerprints. Consequently, we applied two analyses (with PDB and MOL2 files). The same holds true for FuzCav, Shaper, and VolSite/Shaper. The algorithm underlying TM-align relies on a given residue sequence alignment. Binding site atoms within 10 Å of the bound ligand were used to ensure a meaningful comparison. As the tool relies on an initial sequence alignment of cavity residues, a small number of residues prevents the alignment of the residues of interest. The ranking lists of the binding site comparison tools were used to investigate their performance with respect to AUC and EF. The ROC curves were plotted with the help of the KNIME[103] ROC Curve node to analyze sensitivity (true positive rate, Eq 1) and specificity (true negative rate, Eq 2) of the tools. The AUC values for the resulting ROC curves were also calculated using KNIME. Statistical analyses of the AUC differences were performed according to DeLong et al. [109] as implemented in the R[110] package pROC[111]. The EF describes the enrichment of similar (active) binding site pairs opposed to the number of similar pairs identified in a random screening (Eq 3). The EF for x% of the screened data set is calculated based on the number of true actives at this percentage (N (actives) x%) and the number of all pairs at this percentage (Nx%) in the list of pairs with ranked similarity/distance score, the number of true active pairs (N (actives) 100%), and the number of all pairs in the complete data set (N100%). For tools with more than one scoring scheme, we analyzed the early EFs. The score that led to the highest early enrichment for the identical structures was taken into account. If no distinction was possible, the results for the data set of NMR structures were used in a similar manner. Finally, the data set of decoy structures was taken into account. The applied scores can be found in Table 9. Notched box plots of scores for active and inactive pairs in the data set of structures with identical sequences were generated using the software package R[110]. The Welch’s two-sample t-test[112] for the similar and dissimilar pair score distributions was performed using the software package R[110] The Spearman’s rank correlation coefficient (Spearman’s Rho, rS) was calculated for the decoy binding site data sets according to Eq 4. The raw scores (Xi) and the number of mutations (Yi) are converted to ranks (rank X, rank Y). The covariance of the rank variables divided by the product of the standard deviations of both rank variables gives the final correlation coefficient. In this study, the correlation between binding site similarity score and number of binding site mutations was calculated. The general expectation is that the higher the number of binding site mutations the lower the score. Thus, if the Spearman’s Rho equals -1, it denotes a perfect correlation. Optimum similarity score cut-off values for each method were determined using the R[110] package pROC[111]. Youden’s J statistic[85] was applied without weights to derive a score threshold that optimizes both sensitivity and specificity for the corresponding data set. The thresholds based on data set 1 and data set 7 were used to determine the methods’ sensitivity and specificity for data set 7. All-against-all comparisons with the structures of data set 5 were performed on an Intel Xeon workstation (E5-2690 with 2. 90GHz and 32 GB RAM) in a serial manner (single core). The time for the pocket preparation was disregarded for SiteAlign, PocketMatch, and TM-align as this was realized separately. It has to be considered that for the remaining binding site comparison software, this is achieved within the given run time for preparation. For ProBiS, preparation and comparison are performed on the fly. Additionally, the time for summarizing the results is not included. The run time was assessed with the Linux “time” command (user time).
Binding site similarities are useful in the context of promiscuity prediction, drug repurposing, the analysis of protein-ligand and protein-protein complexes, function prediction, and further fields of general interest in chemical biology and biochemistry. Many years of research have led to the development of a multitude of methods for binding site analysis and comparison. On the one hand, their availability supports research. On the other hand, the huge number of methods hampers the efficient selection of a specific tool. Our research is dedicated to the analysis of different cavity comparison tools. We use several binding site data sets to establish guidelines which can be applied to ensure a successful application of comparison methods by circumventing potential pitfalls.
infographics
chemical characterization
protein interactions
social sciences
archaeology
protein structure
research and analysis methods
sequence analysis
computer and information sciences
sequence alignment
bioinformatics
proteins
biological databases
proteomics
binding analysis
molecular biology
protein structure comparison
biochemistry
proteomic databases
data visualization
database and informatics methods
graphs
biology and life sciences
macromolecular structure analysis
A benchmark driven guide to binding site comparison: An exhaustive evaluation using tailor-made data sets (ProSPECCTs)