idTarget: a web server for identifying protein targets of small chemical molecules with robust scoring functions and a divide-and-conquer docking approach.

Identification of possible protein targets of small chemical molecules is an important step for unravelling their underlying causes of actions at the molecular level. To this end, we construct a web server, idTarget, which can predict possible binding targets of a small chemical molecule via a divide-and-conquer docking approach, in combination with our recently developed scoring functions based on robust regression analysis and quantum chemical charge models. Affinity profiles of the protein targets are used to provide the confidence levels of prediction. The divide-and-conquer docking approach uses adaptively constructed small overlapping grids to constrain the searching space, thereby achieving better docking efficiency. Unlike previous approaches that screen against a specific class of targets or a limited number of targets, idTarget screen against nearly all protein structures deposited in the Protein Data Bank (PDB). We show that idTarget is able to reproduce known off-targets of drugs or drug-like compounds, and the suggested new targets could be prioritized for further investigation. idTarget is freely available as a web-based server at http://idtarget.rcas.sinica.edu.tw.

INTRODUCTION

Identification of targets of small chemical molecules is essential for unravelling the underlying molecular causes of actions. Often, natural products, i.e. compounds discovered from plants, animals, marine lives or other living organism, exhibit useful pharmaceutical effects, e.g. anti-inflammatory, anti-cancer and anti-viral effects, yet their molecular mechanisms remain elusive. On the other hand, many drugs are known to be accompanied with unpleasant adverse effects, but the molecular targets of such effects are largely unknown. On the contrary, there are also some old drugs whose additional beneficiary effects are discovered only recently. For example, the epigenetic mechanism of the anticancer effect of cholesterol-lowering drugs, statins, was uncovered rather recently (1).

Conventional virtual screening of chemical libraries has been used widely to search for new leads in drug development for a protein target (2). As the deposited structures of biomolecules in the Protein Data Bank (PDB) increase substantially in the past decades, searching for the targets of a given drug or small compounds (also known as inverse screening, target fishing, off-target prediction, etc.) has become a useful approach (3–7).

One of the major hurdles for target identification is the effectiveness of scoring functions (7,8). To evaluate the binding affinity of the small ligand and a protein target, an accurate yet generally applicable scoring function is essential. We recently developed three robust scoring functions, AutoDock4RRP, AutoDock4RAP and AutoDock4RGG (9) based on the energetic terms and the formulation of AutoDock4 (10). These scoring functions report the binding free energy in the experimental scale, which allows direct comparison among different protein–ligand systems. Two of these three robust scoring functions were constructed using atomic charges from quantum chemical calculations, namely, RESP (11) and AM1-BCC (12), and the robust regression analysis (13) was employed to mitigate the influence of outliers for the calibration of the scoring functions. These robust AutoDock4 scoring functions have been benchmarked for their capability in binding affinity prediction and binding pose prediction (9). For the assessment of binding affinity prediction with a large external set of 1427 complexes from PDBbind v2009, AutoDock4RAP obtained root-mean-square errors of 2.176 kcal/mol, while the size of the training set is only 147. Benchmarked by using two decoy sets (14,15), the robust AutoDock4 scoring functions outperformed most of other scoring functions for the binding pose prediction (9).

Here, we utilize an efficient docking approach to screen the protein targets. Evaluation of potential targets is carried out by using the AutoDock4 robust scoring functions and the affinity profile analysis to enhance the confidence level of prediction.

MATERIALS AND METHODS

Docking and scoring

The search engine of idTarget web server is MEDock (16), which generates initial docking poses of the small ligand. The global search algorithm used in MEDock has also been tested recently by random mathematical functions simulating rugged free energy landscapes with different dimensionalities (17). It was shown that this global search algorithm maintained very high searching efficiency even at the dimensionality of 30 (17), which should be sufficient for applying to most protein–ligand systems. It was also shown that the traditional genetic algorithm failed to deliver good searching efficiency as long as the dimensionality is beyond 20 (17).

Since the binding sites are generally not known, the search space (i.e. the docking grid), should cover the entire protein. However, the searching efficiency also depends on the size of the search space. To achieve better searching efficiency, as well as to overcome the memory limitation when docking to very large proteins, a divide-and-conquer docking approach, similar to the partial box approach in BDT (18), was adopted for efficient blind docking. It should be emphasized that the entire receptor surface is searched, instead of only a part of the surface or just the active site. This also allows the search for possible allosteric binding sites. The divide-and-conquer docking approach can easily take the advantage of the parallel computing facilities and reduces the searching time by limiting the size of grids. The large box that covers the entire protein surface is subsequently divided into smaller boxes, where the size of the smaller boxes is dynamically determined according to the size of query ligand. If a grid box is far away from the receptor, namely, no atom is within 1.42 times the length of the grid box to the center, it will be eliminated to reduce the computational cost. The pre-calculated energetic grid maps for the protein targets are stored in the idTarget web server, and therefore the time for grid construction can be saved.

To analyze the results of target screening, the affinity profile of each binding site is constructed. Our affinity profile analysis was inspired by the structure-based maximal affinity model recently proposed (19). Here, instead of maximum affinity, the affinity profile of a given pocket is modelled by a Gaussian function, whose width is determined by the range of the predicted binding affinities of the complexes of different ligands in the same binding pocket of the same protein. A Z-score of ligand j to a given protein pocket i is calculated based on above profile using the following equation: where E is the dock energy of ligand j docking to the protein pocket i. E represents the energy directly evaluated with the crystal pose. E and sd are the center and the width of the affinity profile of protein pocket i, respectively. A large negative Z value signifies an important target of the query compound.

The contraction-and-expansion strategy

Since direct screening for a large set of targets is computational intensive, a contraction-and- expansion strategy is adopted to perform target screening efficiently. In the contraction stage, the set of representative pockets for the protein–ligand complexes in the PDB, 3046 mean points of scPDB (20,21), were clustered by CD-Hit (22) under the cutoff of 40% sequence identity. After clustering of protein sequences in the PDB, a contracted list of targets (2091 targets) was obtained. The targets screening starts from this contracted list, subsequently, half of targets with lower docking energy are kept. In the expansion stage, targets that are homologous or contain similar binding sites to the top representatives are also selected for screening. The idTarget server is periodically updated to include new protein structures from PDB. The pockets of a new PDB entry will be compared with the set of representative pockets by SiteEngine (23) and subsequently assigned to the most similar pocket group.

In addition to the contraction-and-expansion strategy, idTarget provides two modes for searching binding poses. In the ‘scanning mode’, conventional docking jobs are performed for each protein structure. In the ‘fast mode’, the ligand will first be mapped to the binding sites of the homologous proteins by superposition of homologous protein structures with the pre-calculated structural alignment information by using CE (24), and then its binding poses are optimized by adaptive local searches to remove too close contacts with protein atoms. To extend the coverage of PDB structures, idTarget also performs a ligand similarity search to find structurally similar co-crystallized ligands. The Tanimoto coefficients calculated by OpenBabel 2.3.0 (25) are used to rank the results of the similarity search by matching molecular fingerprints.

INPUT

The input window of the idTarget web server requires an input ligand file for the target screening. Users can upload a ligand file with the popular molecular file formats (pdb, mol2, pdbqt or cif). Internally, idTarget uses the ‘pdbqt’ file format that contains the information of atomic charges and atom types. Users can select the Gasteiger, AM1-BCC or RESP model to calculate the atomic charges for the uploaded ligand. Otherwise, the default atomic charges will be calculated by the AM1-BCC charge model. idTarget accepts an ligand file with predetermined protonated state by checking the option, ‘Yes, no further addition of polar hydrogens is needed’, or the protonation state of the ligand will be assigned by the idTarget server. The idTarget web server also provides a graphical user interface using the Marvin Sketch applet for directly drawing a molecule on the web interface. As the default setting, the larger scale screening is performed with the idTarget data set. Users can also query a small set of proteins with known structures by editing a PDB ID list.

OUTPUT

Tracking the status of a job and retrieving the results can be achieved by input the Task ID. If the e-mail address is provided, the notification will be sent to users as soon as the calculation is finished. Although it usually takes hours to finish a job, depending on the loading of the server, users can always check the results while some targets have been screened. During the calculation of a job, the screening results are sorted immediately for showing the potential targets in the higher ranks. As depicted in Figure 1, two hierarchical tables and a molecular graphic visualizer are shown in the result page. The two-layer table (grouped by protein name) lists all the results sorted by dock energies. The PDB entries with the same target name will be grouped, and one can click on the button to expand/collapse the next layer or all results at once. The homologous proteins are grouped in the three-layer table (grouped by the binding pocket), where the representative one is shown in the first layer and the proteins within the same binding pocket are shown in the second layer. All the PDB entries with the same target name are collected in the third layer. Users can also use selected ranges of predicted binding free energy and Z-score as the filtering criteria to obtain a more condensed target list.

Figure 1.

An example of target screening for the HDAC2 inhibitor (LLX). The left panel gives two kinds of representations: grouped by protein name and grouped by homology. The right-hand side contains a viewer for visualizing the docking poses online and a panel showing 2D-similar ligands in PDB.

WEB SERVER

The web interface of idTarget was written in PHP (5.1.6). Marvin Sketch applet (http://www.chemaxon.com/products/marvin/marvinsketch/) was embedded in the idTarget job submission page so that users can draw their interested ligand online. In the result page, users can inspect the binding poses and manipulate the resulting 3D structures interactively via the embedded OpenAstexViewer applet (http://openastexviewer.net/web/).

EXAMPLES

To demonstrate the usage of the idTarget web server, three examples are presented here. The first example is the HIV-1 protease inhibitor, darunavir (DRV), which has been approved by FDA on 2006 and has been validated the inhibition against multidrug-resistant viruses (26). The second example is one of indirubin derivatives, 6-bromo-indirubin-3′ oxime (6BIO), which has been experimentally validated as an inhibitor of several kinds of protein kinases (27). The third example is an inhibitor of histone deacetylase 2 (HDAC2), namely, N-(4-aminobiphenyl-3-yl)-benzamide (dubbed ‘LLX’ from its PDB component ID). The target screening results of these three examples are archived in the web server with the Task IDs DRV_screening, 6BIO_screening and LLX_screening, respectively.

DRV

Mutations of HIV-1 protease (PR) residues usually reduce susceptibility to the PR inhibitors but not significantly to darunavir. There have been many structural and kinetic analyses of darunavir against several drug-resistant mutants (26,28–30). Table 1 gives the comparison of predicted and experimental binding affinities for the wild-type and mutant HIV-1 PR. These results indicate the comparable predicted binding affinity to the experimental data. Although the scoring function was not trained with these PR-inhibitor complexes, it still can discriminate the fine difference of affinity against mutants, except for the case of 3CYW. It should be noted that in the docking, we have excluded a well-known flap water molecule, which forms a hydrogen bond network between the flaps of PR and ligand. If that flap water was included in docking of DRV, the predicted free energy will be lowered by ∼1 kcal/mol. This flap water was removed because in general it may cause steric clashes for most ligands. Figure 2 depicts the binding pose of DRV against the wild-type PR with the rmsd of 0.85 Å compared to the crystal pose.

Table 1.

Darunavir screening results for wild-type and mutant HIV-1 PR

PDBID	ΔGpred.	Kipred.(nM)	Z-score	mutant	Kiexp. (nM)
3CYW	−10.99	8.8	0.26	G48V	17a
2F81	−10.82	11.7	0.32	L90M	0.03b
2IEN	−10.79	12.3	0.33	wild	0.22–1.0a,b,c
2IDW	−10.65	15.6	0.38	V82A	0.8–1.3b,c
2IEO	−10.55	18.5	0.41	I84V	3.2b,c
1T3R	−10.53	19.1	0.42	wild	0.06d
2HS1	−10.47	21.1	0.44	V32I	3.3e
2F80	−10.4	23.8	0.46	D30N	6.6b
2F8G	−10.38	24.6	0.47	I50V	2.0–18a,b
3D1Z	−10.38	24.6	0.47	I54M	1.6a
3D20	−10.38	24.6	0.47	I54V	5a

aExperimental data from J. Mol. Biol. (2008) 381, 102–115.

bExperimental data from J. Med. Chem. (2006) 49, 1379–1387.

cExperimental data from J. Mol. Biol. (2004) 338, 341–352.

dExperimental data from J. Med. Chem. (2005) 48, 1813–1822.

eExperimental data from J. Mol. Biol. (2006) 363, 161–173.

Figure 2.

The docking pose of darunavir (in magenta) compared to the x-ray pose (in cyan) in the wild-type protease (PDB ID: 2IEN) gives a RMSD of 0.846 Å.

6BIO

This example shows that not only the well-known kinase targets (CDK2, CDK5 and GSK-3β) could be identified by idTarget, but also the kinase PDK1 (27) and other proteins. More detailed processes are described in ‘Example 1’ in our idTarget web server. The predicted best targets are as follows: CDK2 (PDB ID: 1FVV, predicted Ki: 51.7 nM); CDK5 (PDB ID: 1UNH, predicted Ki: 71.2 nM); GSK-3β (PDB ID: 1Q41, predicted Ki: 63.3 nM); PDK1 (PDB ID: 2PE1, predicted Ki: 96.5 nM). Figure 3 depicts the best docking pose of 6-bromo-indirubin-3′ oxime (6BIO) in PDK1 that confirms the proposed pose in Zahler et al. (27).

Figure 3.

The best docking pose of 6-bromo-indirubin-3’ oxime (6BIO) in PDK1.

To compare directly with the results of Zahler et al. (27), we performed the screening with the same set of 5821 PDB entries by using the ‘User-edited list’ feature. The predicted binding free energy was used to rank the screening results, and the Z-score (<0) was the filter. For the known protein kinase targets (CDK2, CDK5 and GSK-3β), the top 1% and top 5% enrichment factors in our screening results were 6.45 and 11.94, respectively, while the top 1% and top 5% enrichment factors calculated from the results of Zahler et al. were 19.35 and 7.42, respectively. It should be noted that protein kinases are the only targets of interest in the work of Zahler et al. However, we found that two proteins, matrix metalloproteinase and glycogen phosphorylase, constitute the major part of the top 1% of our results. Interestingly, we found that a recent experimental work (31) showed that an indirubin derivative inhibits rabbit muscle glycogen phosphorylase b. The enrichment factor for glycogen phosphorylase (PYGM) in the top 1% was 18.18 in our results, but none of PYGM was identified in the top 5% of the results of Zahler et al. We further performed the ROC curve analysis (32) to assess of the performance of the prediction. Figure 4 shows the ROC curves and numerical values of the areas under the curve (AUC), where CDK2, CDK5 and GSK-3β are considered as known targets of 6BIO (in Figure 4a), and these three kinases plus PYGM are considered as known targets (in Figure 4b), respectively. Supplementary Figures S1–S4 show the ROC curves with different selections of known targets and Z-score filters. The more the stringent of the Z-score filter, the higher the AUC was obtained, which could be due to the nature of conformational selectivity, the general applicability and the universal physical scale of our robust scoring function (9). These analyses also indicate that our predictions are comparable with the prediction of Zahler et al.

Figure 4.

The red curve is obtained from the prediction of Zahler et al., while the black curve is obtained from the prediction of idTarget with Z-score < 0 set as the selection filter (1161 protein targets selected). (a) CDK2, CDK5 and GSK-3β are considered as known targets and others are considered as decoys. (b) CDK2, CDK5, GSK-3β and PYGM are considered as known targets and others are considered as decoys.

LLX

N-(4-aminobiphenyl-3-yl)-benzamide is a compound recently designed and synthesized to optimize the fitting of the HDAC2 binding pocket (33). It has been long recognized that there is a 14 -Å long internal cavity, also called ‘foot pocket’, adjacent to the catalytic zinc ion of HDAC. This foot pocket connects to a lipophilic tube with an angle of about 110°, and this compound sits at the kink between the foot pocket and the tube. As shown in Figure 1, the screenshot of the result page given by idTarget, the compound LLX recognized its co-crystalized protein structure (PDB ID: 3MAX) with the highest affinity.

DISCUSSION

Although darunavir is a relatively new drug, its adverse effects and interactions with other drugs have been reported (34,35). It may be possible to further investigate possible molecular causes of actions of these adverse effects and drug–drug interactions by scrutinizing the list of proteins identified by idTarget. Our prediction for possible targets of 6BIO is generally consistent with the previous work. (27). It is interesting to note that the top prediction of LLX is 3MAX, which is the HDAC2 structure crystallized by the same group that designs and synthesizes this N-(4-aminobiphenyl-3-yl)-benzamide. A possible explanation could be that this compound sits at the kink between the foot pocket and the lipophilic tube, which altogether forms a very curved channel structure and may be rather unique among protein structures.

Drugs can bind to a myriad of particles in the blood, such as red blood cell, leukocytes, platelets, as well as to some proteins, e.g. serum albumin, lipoproteins, α1-acid glycoproteins, α, β, γ-globulins, etc. (36). Usually, the metabolites of a drug (or a small molecule) can also bind to other human proteins to exert unwanted effects, or even bind to the original target protein and deliver similar actions. In view of this, target identification should be employed in tandem with metabolite prediction algorithms, e.g. MetaSite (37). Many drugs, e.g. vinblastine, verapamil, taxol, etc., are substrates of P-glycoprotein, whose low resolution structure has been available recently (38), and higher resolution structures should be on the way (G. Chang, personal communication). Due to the tremendous cost of large-scale protein binding assays, the comprehensive list of the proteins that a given compound binds is generally not available. One should also bear in mind that proteins are localized in compartmented organelles, and a given compound may not be able to be distributed to a specific compartment where a certain protein exists, and therefore the prediction of target identification should also be guided by the knowledge of protein subcellular locations to remove the false positives if only the physiological scenario is of interest.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Figures 1–4.

FUNDING

The National Science Council of Taiwan [NSC 100-2627-B-001-001, 99-2627-B-001-002, 98-2627-B-001-002, 98-2323- B-002-001, 98-2323-B-077-001, 97-2323-B-002-015, and 97-2923-M-001- 001-MY3]. Funding for open access charge: Research Center for Applied Sciences, Academia Sinica.

Conflict of interest statement. None declared.