A QuantumBio QuickStudy is a concise example, tutorial, or “mini-study” which answers a specific scientific or technical question using QuantumBio methods and tools. The following QuickStudy was suggested by Dr. Sirish Kaushik Lakkaraju. If you have a question you would like to see answered in a future QuickStudy please email info@quantumbioinc.com and suggest it.
- Question: How well does MovableType (MTScoreE and MTScoreES) perform with the CatS set from D3R?
- Download the input/output for each step and follow along.
Introduction
On February 6, 2020, the Cathepsin S (CatS) Dataset from Grand Challenge 3 was provided by the Drug Design Data Resource (D3R). The CatS set consists both of a set of X-ray structures corresponding to 19 reported binding affinities shown in the following table along with a set of ~117 additional ligands provided in SMILES format.
Table 1: Structure ID’s, Ligand ID’s, and pIC50 data as provided by D3R.
| StrucID | Ligand | IC50[um] | IC50[M] | pIC50 |
|---|---|---|---|---|
| gmad | B8J | 0.055 | 5.50E-08 | -7.26 |
| gosb | B8S | 0.056 | 5.60E-08 | -7.25 |
| haan | BC7 | 0.16 | 1.60E-07 | -6.80 |
| hhcf | B8Y | 0.11 | 1.10E-07 | -6.96 |
| jjod | BJS | 0.015 | 1.50E-08 | -7.82 |
| mekm | BCJ | 0.058 | 5.80E-08 | -7.24 |
| pbwz | BJJ | 0.12 | 1.20E-07 | -6.92 |
| pmgn | FRG | 0.2 | 2.00E-07 | -6.70 |
| qjan | BFV | 0.29 | 2.90E-07 | -6.54 |
| tjyg | BJY | 0.2 | 2.00E-07 | -6.70 |
| ttlc | B8V | 0.2 | 2.00E-07 | -6.70 |
| tuhd | BHV | 0.17 | 1.70E-07 | -6.77 |
| uavh | BHJ | 0.03 | 3.00E-08 | -7.52 |
| wcgq | B9Y | 0.23 | 2.30E-07 | -6.64 |
| yqqd | BJV | 0.14 | 1.40E-07 | -6.85 |
| yquj | BGJ | 0.064 | 6.40E-08 | -7.19 |
| yrpk | BG7 | 0.26 | 2.60E-07 | -6.59 |
| znpb | BAJ | 0.03 | 3.00E-08 | -7.52 |
| zrqq | B9S | 0.074 | 7.40E-08 | -7.13 |
Step 1: Provided structures and MTScore
Beginning with the 24 PDB files found in GC3_CatS/CatS/CatS_poses/chainA of the CatS zip file, we removed the 5 cases which did not have provided IC50 data listed in the GC3_CatS/CatS/final_CatS_score_compounds_D3R_GC3.csv.csv file. This IC50 data was converted to pIC50 for use in subsequent analysis as reported in Table 1 above.
Each of these 19 structures was treated with the following two routines:
$ moebatch -run /path/to/DivConSuite/svl/run/qbDockPair.svl
-rec ${targetbasename}.pdb -delwat -protonate
$ qmechanic pro_<ligand>_predock.pdb –ligand lig_<ligand>_predock.mol2
-h garf –mtdock <ligand>_dock.sdf –mtscore ensemble –np 2 -v 2
When moebatch is executed, qbDockPair.svl reads the PDB files provided in GC3_CatS and docks any included ligand within the active site (note: extraneous species like SO4, GOL, etcetera were removed prior to docking). The -protonate command line option performs a Protonate3D step and the -delwat option removes all water molecules (after protonation is performed). Once these steps are complete, MOE’s built in docker generates 25 ligand binding modes or poses using a rigid-receptor/flexible-ligand docking protocol. This qbDockPair.svl workflow encapsulates a number of modifications to MOE-defaults. Specifically, 100 ligand poses are selected and optimized within the pocket using the Amber10:EHT potential (instead of DEFAULT: 30), and the top 25 poses are selected using MOE GBVI/WSA dG score (instead of DEFAULT: 5). These settings increase time requirements for docking, but they also increase the amount of sampling performed during docking. These settings may be increased further by using the -maxpose option for increasing the number of poses to optimize and the -remaxpose option for increasing the number of poses to return.
After moebatch completed docking each ligand, qmechanic iswasexecuted on the output files generated from the qbDockPair.svl script. The pro_<ligand>_predock.pdb and lig_<ligand>_predock.mol2 files were those provided by qbDockPair.svl and generated based on the receptor and ligand respectively. Finally, the <ligand>_dock.sdf file includes 25 docked poses for this run. The MT-GARF pair potential was selected for this study, and moebatch-generated poses were supplied to MTScore for an ensemble MT score.
The MTScoreES analysis yields a Pearson-R = 0.58 (while the Ensemble score with 25 poses has a Pearson-R = 0.51). While this is a “good start” it is by no means a limit to how well the method is expected to perform. We could certainly explore the provided X-ray models and perform an outlier analysis, but often we have seen that initial X-ray structures have “hot spots” and other problems which can impact results. With that in mind, we reran the docking with some additional settings to see if we could clean up the prediction.
Step 2: Impact of improved structure minimization on MTScore
Again, beginning with the 19 structure files noted in Table 1, we added the -preopt flag to the call to qbDockPair.svl. As the argument keyword implies, the -preopt flag optimizes the protein ligand complex (along with any water molecules, ligands, and so on) using the Amber10:EHT potential as implemented in MOE.
$ moebatch -run /path/to/DivConSuite/svl/run/qbDockPair.svl
-rec ${targetbasename}.pdb -delwat -protonate -preopt
$ qmechanic pro_<ligand>_predock.pdb –ligand lig_<ligand>_predock.mol2
-h garf –mtdock <ligand>_dock.sdf –mtscore ensemble –np 2 -v 2
In Step 1, the pro_<ligand>_predock.pdb and lig_<ligand>_predock.mol2 files were effectively recapitulations of the original input PDB files plus protons. In Step 2 however, the pro_<ligand>_predock.pdb and lig_<ligand>_predock.mol2 files are not only protonated versions but they are optimized as well and the <ligand>_dock.sdf file includes the 25 poses for these optimized versions. No other changes were made in this run.
For Step 2, we observe some improvement to MTScoreES (from Pearson-R=0.58 to Pearson-R = 0.61). This suggests we are likely on the right track. When considering the ensemble score (MTScoreE) however, this additional minimization step improves our prediction substantially and the Pearson-R moved from 0.51 to 0.64 with fewer outliers than in Step 1.
Step 3: Outlier Cleanup and Conclusions
Finally, we can further increase the ensemble Pearson-R = 0.73 (from R=0.64) with a slope of 0.87 by removing a single outlier ligand (BAJ). The end-state Pearson-R remains unchanged with removal of the outlier and has a slope of 0.64. Given that 0.73 is significantly better than 0.61 and the slope is likewise approaching 1, this observation would suggest that for this set, the crystal structures are not the best/only representations and don’t provide the most sampling opportunities for the method.
Next: Follow up analyses could be performed to explore why BAJ exhibited problems with MT and hurt the overall prediction. It could be that this is a limitation of the GARF potential, or the perhaps the docker had more trouble with this particular ligand leading to questionable poses. The observation that the ensemble score is significantly more predictive than the end state score suggests that additional sampling for this set would be worth exploring. Specifically, would molecular dynamics (MD) snapshots lead to greater predictive performance? Finally, there are 117 additional ligands listed in the final_CatS_score_compounds_D3R_GC3.csv.csv which could be docked/scored.