Evaluation of DOCK 6 as a Pose Generation and Database Enrichment Tool

Scott R. Brozell,§, 1 Sudipto Mukherjee,§, 2 Trent E. Balius,2 Daniel R. Roe,1 David A. Case,1 and Robert C. Rizzo*, 2, 3

1 BioMaPS Institute and Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, New Jersey 08854, USA

2 Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, New York 11794, USA

3 Institute of Chemical Biology Drug Discovery, Stony Brook University, Stony Brook, New York 11794, USA

§ These authors contributed equally to this work

* Corresponding author e-mail: , phone: 631-944-2891, fax: 631-632-8490

SUPPLEMENTAL INFORMATION

Ligand Atom Mapping Algorithm. To address issues with organizer supplied ligand data not being in the same atom order and not having the same atom name (which would make mapping trivial), an atom mapping routine was implemented into the program cpptraj (part of the AmberTools suite of programs). The basic algorithm uses an iterative procedure that consists of four steps, illustrated in Figure S1. The first step is to identify unique atoms in both the reference molecule and the target molecule (Figure S1, blue atoms). An atom is considered unique if its local environment (the atom, plus the atoms it is bonded to, plus the atoms they are bonded to; i.e., up to an adjacency of two) is not repeated anywhere else in the structure. The second step is to identify atoms that are unique by virtue of being the only ones of their kind bonded to unique atoms, e.g. the only hydrogen atom bonded to a unique carbon atom (Figure S1, red atoms). Any unique target atoms that match unique reference atoms are then mapped. The third step is to identify atoms in the target and reference that are bonded to partially mapped chiral centers (i.e., the chiral atom plus two atoms it is bonded to have been mapped) using improper dihedral calculations (Figure S1, orange atoms). This algorithm therefore preserves chirality when mapping two molecules.

At this point, the algorithm returns to the first step to see if any more atoms have become unique given the current atom mapping. If no new atoms have been mapped after the third step, the final step is to match remaining unmapped atoms based on their element and bonding to already-mapped atoms (Figure S1, green atoms); if any atoms are mapped this way the algorithm returns to the first step (Figure S1, cyan atoms), otherwise mapping is finished. If mapping is incomplete (which can result if the target molecule does not have the same number of atoms as the reference molecule) the algorithm simply maps what atoms it can and reports the discrepancy to the user. In this study, several targets were missing a hydrogen atom compared to the reference. However, this does not impact success rate calculations since no hydrogens are included in rmsd calculations. One issue that can arise with this method is that it may not be possible to identify any unique atoms in the first step; this can happen if, for example, the molecule to be mapped is highly symmetric (such was the case with the ligand from 1tz8). When this happens, a list of atoms is made which are the “most unique” (i.e., have the smallest number of matches) and mapping is tried using these atoms as a first guess. Although the algorithm is quite simple compared to more complex graph-matching algorithms (e.g., see references [1-3]), all ligands from the organizer supplied set were able to be mapped.

Figure S1. Figure illustrating the procedure followed by the atom-mapping algorithm in cpptraj for the ligand from 1g9v.

ROC Curve Examples with Missing Data. As discussed in Computational Details it is often necessary to account for systems in which some active or decoy molecules do not produce a final answer [4,5]. DOCK completion statistics for this study are shown in Tables S6-S7. From a practical standpoint there are several possible approaches to deal with such systems. Figure S2a demonstrates a situation whereby a ROC curve is generated ignoring molecules for which a final answer was not obtained (black line). In this scenario the number of positives/actives (P) and negatives/decoys (N) employed to compute the TPR and FPR rates become Pdocked and Ndocked and the ROC curves reaches TPR=1, FPR=1. In contrast, the gray line in Figure S2a uses the underlying data but number of actives and decoys used in the calculation are based on the initial number (Pinitial and Ninitial) regardless of whether they actually produced a final result. Note that in this case the ROC curve does not reach TPR=1, FPR=1 which could be problematic for making AUC comparisons among different systems. Figure S2b demonstrates three different ways to force a ROC curve to achieve TPR=1, FPR=1 in which missing data can be assumed as yielding perfect enrichment (blue upper line), no enrichment (lower red line), or random enrichment (purple middle line). For systems in which not all molecules were successfully docked, ROC curves in this report were generated using the reasonable assumption that there would be random enrichment.

Figure S2. (a) Hypothetical ROC curves computed using two different values for the total number of molecules classified as active (P, positive) or decoy (N, negative). The gray curve was computed using Pinitial and Ninitial and the black curve was computed using Pdocked and Ndocked. (b) Missing data can be assumed as yielding perfect enrichment (blue upper line), no enrichment (lower red line), or random enrichment (purple middle line) to ensure the ROC curve will reach TPR=1, FPR=1. The dashed line is the random ROC curve.


DETAILED SYSTEM PREPARATION PROTOCOLS

(1) ASTEX PDB. For consistency with the recently described SB2010 [6] database, systems in Astex dataset were built from scratch using raw coordinates downloaded from the Protein Data Bank. Improvements over the earlier SB2010 protocol include MM energy minimization of receptor hydrogens in the presence of the ligand, and support for small molecule co-factors. However, the 1T9B structure was excluded because insufficient electron density was observed for the cofactor in the binding site. Only the first occurrence of each binding site was included in the setups. Protonation states for ligands were assigned by hand by consulting the corresponding crystal structure reference. The total number of DOCK-ready structures using ASTEX PDB protocol was 84.

(2) ASTEX SUP. The organizers requested that participants use the supplied re-refined structures "as is". Although there were numerous problems with some of the supplied structures, the spirit, if not the letter of that request was followed as much as possible. In contrast to the ASTEX PDB protocol which only included the first occurrence of a binding site, the ASTEX SUP protocol incorporated all viable binding sites (N=147). The DockPrep tool in Chimera 1.2470 [7-9] was used to add AM1BCC [10,11] partial charges to the ligand and assign Amber FF99 charges for the receptor. Four of the 85 Astex receptors failed during automated preparation due to missing TER cards: for 1of1 SO4 was removed; for 1w1p SO4 and GOL were removed; for 1w2g and 1yv3 TER cards were added. Four Astex ligands also failed during automated preparation due to incorrect reading of the sdf charge cards. For 1mmv, 1uou, and 1yv3, Chimera did not estimate the charges correctly resulting in odd numbers of electrons, so these were rerun manually using antechamber with the correct charges. For 1q41, Chimera added a hydrogen, so this was removed and antechamber was rerun manually with the correct charge. In general, protonation states for the supplied ligands remained unchanged.

(3) ASTEX AMB. Use of Amber score requires many additional input files beyond those of the other DOCK scoring functions. The general Amber score procedure has been previously described [12] so only some context and new developments for DOCK are presented. Amber score preparation in DOCK uses three amberization scripts, one each for the ligand, receptor, and complex. These scripts invoke the LEaP and antechamber programs to generate the input files containing the topology, force field parameters, and initial coordinates. DOCK Amber score now provides four libraries that are automatically sourced in the amberization steps: cofactors, ions, modified residues, and a user definable library. The cofactor library was substantially expanded for this work to forty molecules, and it now uses a scripted protocol, for consistency and robustness, that downloads structures from the PDBeChem Ligand Dictionary [13] and executes antechamber to assign GAFF parameters and AM1BCC charges. The ions library is a corrected and expanded version of Amber’s ions94 library. The modified residues library was created for this study from the four nonstandard residues in the Astex data set; it is under development. The user definable library exists to override the other libraries, but it was extensively used in this investigation to correct nonstandard atom names. The 1994 parameter set was used for the protein atoms [14], and the Generalized Amber Force Field (GAFF) parameter set for ligand and cofactor atoms was used [15]. Upon processing, 65 of the 85 Astex receptor files supplied by the organizers initially failed mostly due to atom naming issues, missing TER cards, bad CONECT records, and replacing HETATM with ATOM. For expediency in these 65 all hydrogens were removed and replaced using LEaP with its defaults (i.e., histidine with hydrogen on the epsilon nitrogen), TERs were added manually as needed, and CONECTs were removed manually as needed.

(4) DUD PDB. This protocol is the same as the ASTEX PDB protocol in which SB2010 [6] procedures were employed to set up the receptors. Thirty eight out of forty systems were successfully prepared using this protocol. Two systems (1ah3 and 1l21) could not be added owing to issues with the coordinates deposited in the PDB. For this preparation, active and decoy ligands were used directly as downloaded from the DUD website and for the 10 systems re-run using WOMBAT actives, the WOMBAT ligands were assigned AM1BCC partial charges with the MOE program and combined with the appropriate DUD decoys.

(5) DUD SUP. The DUD dataset provided by the organizers, contained protonated versions of the original DUD receptors, but in many cases the partial charge assignments appeared problematic. For example, many atoms had missing (0.000) partial charges in the supplied mol2 files, and for many proteins the sums of the charges for a given residue were non-integral. In addition, the atom records had nonstandard names meant for small molecules (e.g., c1, c2), not the standard atom names which allows recognition as protein atoms. Naming is especially important for proteins because force field assignments depend wholly on the atom names. Furthermore, even if Sybyl atom types are fully specified in the mol2 files, the combination of incomplete partial charges and nonstandard atom names is problematic for completion of the partial charge assignments because the Sybyl atom types cannot be mapped back to the protein atom names. The MOE program was used in an attempt to assign partial charges, but because of naming issues default Gasteiger [16] charges were assigned. Identical procedures used to setup active and decoy ligand sets for the DUD PDB preparation were used here for the DUD SUP protocols.

SUPPLEMENTAL RESULTS

Multi site Breakdown. As an extension of the multi site docking analysis (Table 3) a breakdown based on site index and the number of binding sites in each protein was also performed. Table S1 contains statistics for the RU/ASTEX SUP multi site runs (one run each DOCK random seed 0) categorized by the number of sites per receptor. The organizers' pdbid_lig.sdf files define both the total number of sites and the order of sites: 1/1 denotes the set of the singly sited receptors, 1/ ≥ 2 denotes the set of the first sites of the receptors with 2 or more binding sites, 2/ ≥ 2 denotes the set of the second sites of the receptors with 2 or more binding sites, etc. Although there is not too large a variation across these sets, considering the relatively small number of sites in each set, Table S1 reveals that docking to a given site of a multiple site receptor does not guarantee success in docking to another of its sites.

Table S1. Astex Pose Reproduction Statistics of Rmsd values of the Top Grid Scored Poses by Sets of Site Index.
Site index /
No. sitesa / Subset
Size / Symmetry rmsd (Å) / Success%
Min / Max / Mean / Std Dev / Median
1/1 / 48 / 0.24 / 8.60 / 2.11 / 2.28 / 0.99 / 64.6
1/ ≥ 2 / 37 / 0.28 / 15.25 / 2.52 / 3.19 / 1.03 / 67.6
2/ ≥ 2 / 37 / 0.24 / 10.37 / 2.20 / 2.58 / 0.89 / 70.3
1/ ≥ 3 / 9 / 0.51 / 4.63 / 1.61 / 1.35 / 0.77 / 66.7
2/ ≥ 3 / 9 / 0.24 / 5.42 / 1.50 / 1.56 / 0.89 / 77.8
3/ ≥ 3 / 9 / 0.37 / 4.72 / 1.89 / 1.52 / 0.86 / 55.6
1/ ≥ 4 / 8 / 0.51 / 4.63 / 1.72 / 1.40 / 0.93 / 62.5
2/ ≥ 4 / 8 / 0.24 / 5.42 / 1.58 / 1.64 / 0.94 / 75.0
3/ ≥ 4 / 8 / 0.37 / 4.72 / 1.71 / 1.51 / 0.79 / 62.5
4/ ≥ 4 / 8 / 0.32 / 10.84 / 2.27 / 3.29 / 1.02 / 75.0
aData from RU/ASTEX SUP docking runs. Multi site set.

Potentially Problematic Sites. The organizers noted that 22/85 systems had significant structural problems, including crystal contacts with the ligand, alternate side chain atom contacts with the ligand, and little or no density for residues or cofactors in the active site. Also three ligands had alternate conformations. The total number of binding sites affected is 41/147. The organizers requested an analysis of the impact of these problems on the results. For this set of sites, Table S2 shows the statistics of the top scored rmsds for the grid score RU/ASTEX SUP protocol. For each set, statistics are based on only a single DOCK run (random seed 0). Surprisingly, the success rate obtained of 73.2% for this set of 41 is higher than the 68.0% obtained for the total ASTEX SUP/RU set of 147. Thus, if structural problems are present, the overall docking results, in this case, do not appear to be adversely affected. In addition, the organizers reported that there was a mistake in the stereochemistry for three of the ligands in the data set they provided (1gpk, 1hvy, and 1s3v). The effect of this on the reported results is unclear.