PyRx – Python Prescription

Thanks to Mark Fortner for sharing this article. It showed up today on my virtual screening feed from PubMed, but it’s too good of an article to miss. Similar to AutoGrow, e-LEA3D is using fragments to built lead compounds. Building blocks are extracted from Comprehensive Medicinal Chemistry database by Symyx. It was interesting to find that besides DOCK Blaster there are also TarFisDock and SCFBIO servers that are using DOCK. Cloud Computing For Drug Discovery is an interesting subject and I knew that there are servers such as DockingServer, DrugFinder and GridSystems that are doing excellent job in this area. It’s interesting also to mention that current version of e-LEA3D is using PLANTS for docking and Frog to generate 3D conformations. The most exciting part is that I rediscovered ZODIAC which is also a GUI for PLANTS. ZODIAC is very similar in spirit to PyRx in that it’s an open source state-of-the-art tool for running molecular docking experiments. In summary, both e-LEA3D and ZODIAC are terrific tools for drug discovery.

There is a new article in J. Chem. Inf. Model. that talks about ParaDockS. This is the only other docking program, besides AutoDock 4 and Vina, that is open source. Congradulations to ParaDockS team for making such a great product available! This means that PyRx can some day use ParaDockS and vise versa, so here I’ll summirize main features of ParaDockS. First, it’s using particle-swarm optimizer, similar to PSO@Autodock. Unlike, AutoDock 4, which uses force field-derived (AMBER) objective function, ParaDockS is using knowledge-based potentials. I look forward for the new develoments in this area.

Here are the results of a new PubChem BioAssay AID: 2427 with 42 compounds (22 active, 2 inactive and 18 inconclusive).

AutoDock did a perfect job here giving favourable binding energy for active compounds compared to inactives. As a result, we have a prefect ROC curve here. I made the makers bigger so you can see that all actives are on a vertical x=0 line and all inactives lay on a horizontal y=1 line.

I’ve got another email from My NCBI (Subject: What’s new for ‘BioAssay – Limits: Protein 3D Structure). It was about AID: 2353 – qHTS Assay for Inhibitors of RecQ-Like Dna Helicase 1 with 23 active and 1272 tested compounds. There are 2 Links for Protein Structure corresponding to open and closed conformations. I run virtual screening for both of them, which produced similar results.

This looks more like an ROC curve you’ll see in publications. In contrast to my previous results, this has fewer active compounds. One would think that AutoDock did a better job here since ROC curves are above the diagonal. However, as you can see from the plots on the left, lowest binding energies are coming from inactive compounds. The image below shows molecular surface for 2WWY (chain A) in green and lowest binding energy conformation for CID: 5328760 (-6.03 kcal/mol) and CID: 11957460 (-10.37 kcal/mol). The first one being the most active in the assay shows a weaker (virtual) binding compared to the second one, which although has the best computed binding energy is not active in the assay. The favourable binding energy is partally because of the size; CID: 11957460 is one of the largest in this assay. So the question is: what kind of changes needs to be done to virtual screening to make it agree more with experiment?

When preparing ligands for docking PyRx users can select between PyBabel (MGLTools) and Open Bable partial charges using Edit -> Preferences menu. See PyRx Screencast – Open Bable for details. PyRx is using PyBabel charges, by default. I run the dataset from my previous blog post, using Open Bable, to see if I can get better results. The figure below show the new results.

This shows that PyBable charges perform better for this virtual screening. The reason might be that AutoDock is calibrated using PyBable charges. Another reason might be that I use Open Bable for ligands only and not for the target. Whatever the reason, the new database table plotting tool is a very useful tool when it comes to comparing different virtual screening. 

I’ve added instructions on how to update to the latest version from Source Code and how to submit a feature equest to Downloads page. The major new feature in 0.5 is database table plotting feature that can be used for BioAssays and ROC Curves. To compare ROC curves generated by two different docking algorithms, I run the dataset from my previous blog post through AutoDock Vina. There are third party packages available for AutoDock Vina such as PaDEL-ADV, VcPpt, MakeVinaCommand and PyMOL plugin. The new version of ADT also has features for working with output files generated by Vina. To explore an option for adding a new Wizard to PyRx for AutoDock Vina, I run AutoDock Vina on our cluster, using a script from Vina Manual. Here are the results:

This is just a proof of concept to show that you can use PyRx to analyse results from different docking programmes.

I’ve got an email from My NCBI about new BioAssay – Limit: Protein 3D Structure. AID: 2158 has Compounds Active: 287; Tested: 572; it’s exaclty that I was looking for. Based on this data, I’ve implemented new features that will available in upcoming 0.5 release. Interested users can check it out right away. The new feature allows you to plot Binding Energies to provide a bird’s-eye view of docking results. In addition, if you have corresponding BioAssay available, you can also plot ROC Curve next to Docking Results. ROC curves appear in many recent publications that prompted me to add this feature to PyRx. Google search on ROC Curve brought me to Receiver operating characteristic – Wikipedia. This page confused me more than it helped me, so I started searching for ROC curve further. I came across Let’s ROC that had the following article mentioned in the comments: Triballeau, N. et al. J Med Chem. 2005, 48, 2534-2547. This article is very well written and it helped me to better understand the use of ROC curves in Virtual Screenings.

This plot shows Docking Results and ROC Curve generated for AID: 2158 BioAssay. Note that results with lowest Binding Energy, as predicted by AutoDock, are active compounds. It might seem that overall AutoDock doesn’t do well because some of the points on ROC curve are below the diagonal. However, this all depends on the dataset chosen. If the dataset had only a few active compounds, then chances are, that it would have produced an ideal ROC curve. In that sence, ROC curve, by itself, is not informative, unless you have a plot of Docking Results like the way PyRx does.

I started searching PubChem BioAssays and limited results to the ones that link to Protein 3D Structure. There are currently only 61 BioAssays available. I wanted to select a BioAssay that would have around 100 compounds tested, but only handful that are active. In order to test PyRx and see what new features to add, I have been working with qHTS Assay for Promiscuous and Specific Inhibitors of Cruzain (with detergent) (AID: 1478). There were 617 compounds with the 3D records and I selected 1ME4 as a target for docking. After reading corresponding references from PubChem BioAssay and PDB, I prepared input pdbqt files for ligands and macromolecule. Since there were more that 6 ligand atom types, I needed AutoDock version 4.2 or higher. Luckily there is new Autodock4.2.1 service available.To make PyRx aware of this new service, I’ve added a new page for AutoDock Preferences, where users can change the name of the service, should a new version of AutoDock become available.

I’ve also made the following changes that will be available in the upcoming release of PyRx:

• Modified Remote Jobs Query in webSerives; it now updates the GUI after checking each job individually.
• Made changes to AutoDockPage to make parsing and displaying docking results faster.

You can access the latest version of the code from https://sourceforge.net/projects/pyrx/develop.

As a side note, I was trying to find PubChem Compound ID (CID) for T10. PubChem’s Chemical structure search did not recognize InChI provided by PDB. Fortunately, InChI Resolver from ChemSpider was able to read InChI from PDB (ChemSpider ID:  4451401) and I was able to find correpsonding PubChem ID (CID 5289424). It turns out that InChI versions in PDB and PubChem are not the same. It would be nice if there would be a link in Ligand Summary page that PDB displays that would point to a PubChem ID, whenever corresponding ligand is available in PubChem.

I’ve got a green light from our Office of Technology Development to distribute PyRx under Open Source BSD License. For users, this means that you are free to use PyRx for academic as well as commercial purposes. For potential developers, this means that you can now contribute code to PyRx. I’ve created a project for PyRx at http://sourceforge.net/projects/pyrx and you are welcome to join the project. Visit https://sourceforge.net/projects/pyrx/develop for more information.

I’m happy to announce that PyRx has made into Chemical & Engineering News. Please refer to the following publication if you’d like to cite PyRx in your article:

• L.K. Wolf, C&EN 87, 31 (2009)

Also, I’ll be presenting PyRx during Virtual Screening & Computer Aided Drug Design training session at the Summer Institute organized by National Biomedical Computation Resource (NBCR). Here is a link to the outline titled Getting Started with PyRx.