Tuesday, 19 July 2011

Displaying pharmacophores from Pharmer in Avogadro

Avogadro has a plugin architecture that allows you to drop a Python script into a particular folder to add additional functionality, whether it be a new menu item, a new display method or whatever. I've previously used this to write an extension to read comp chem log files using cclib. This time I need to add a new display method (or Engine in the lingo of Avogadro) that will depict a pharmacophore.

So what's Pharmer? Well, until the release of Pharao by Silicos last year there were no open source pharmacophore programs. And the recent release of Pharmer by David Koes bumps the number up to two. Pharmer focuses on pharmacophore search - ultra-fast pharmacophore search, that is. It does this by precomputing an index that allows very efficient searching. The common denominator of Pharmer and Pharao is not only "Phar", but also that both use the Open Babel library.

The file format used by Pharmer is simply JSON, which is easily edited by a text editor. Here's an example that was created using the ZINC Pharmer interface: HSP90.txt. If you install Avogadro 1.0.1 with Python support on Windows, also have Python 2.6 installed in C:\Python26, and copy the Python script shown below to C:\Program Files (x86)\Avogadro\bin\engineScripts, then you can add the Pharmacophore engine (click "Add" in the "Display Types" panel; note: you can add more than one Pharmacophore engine if you want to display more than one pharmacophore at the same time). Once added, click on Settings and Browse to open the Pharmer JSON file. If nothing appears on the screen, try rotating the display slightly.

Note: The image above was created by joining the output of the following Pharmer command using "obabel --join" and depicting with Avogadro:
pharmer dbsearch -dbdir db_712_actives -in HSP90.txt -out results.sdf -reduceconfs 1 -sort-rmsd

pharmacophore.py:
import os
from PyQt4.Qt import *
import Avogadro as avo
from numpy import array

# Get JSON from my own Python installation
import sys
sys.path.append("C:\\Python26\\Lib\\site-packages")
sys.path.append("C:\\Python26\\lib")
import json

colorlookup = {
  'Hydrophobic': [0, 139, 69],
  'HydrogenAcceptor': [255, 165, 0],
  'HydrogenDonor': [211, 211, 211],
  'Aromatic': [205, 0, 205],
  'PositiveIon': [0, 0, 238],
  'NegativeIon': [255, 0, 0],
  }

class Engine(QObject):
  # declare the changed() signal
  __pyqtSignals__ = ("changed()",)

  # constructor
  def __init__(self):
    QObject.__init__(self)
    self.widget = None
    self.textbox = None
    self.filename = ""
    self.pharma = []

  def name(self):
    return "Pharmacophore"

  def flags(self):
    return avo.EngineFlags.NoFlags

  def debug(self, text):
    return
    filename = os.path.join("C:\\", "Users", "Noel", "Desktop", "tmp.txt")
    output = open(filename, "a")
    print >> output, text
    output.close()

  @pyqtSignature("")
  def browse(self):
    self.filename = str(QFileDialog.getOpenFileName())
    if not self.filename:
      return
    
    self.textbox.setText(self.filename)
    with open(self.filename) as f:
      query = json.load(f)
      self.pharma = [x for x in query['points'] if x['enabled']]
    if len(self.pharma) > 0:
      for p in self.pharma:
        x += p['x']
        y += p['y']
        z += p['z']
      x /= len(self.pharma)
      y /= len(self.pharma)
      z /= len(self.pharma)
      for p in self.pharma:
        p['x'] -= x
        p['y'] -= y
        p['z'] -= z
    self.debug("Pharma: %s" % self.pharma)
    self.emit(SIGNAL("changed()"))

  @pyqtSignature("")
  def clear(self):
    self.filename = ""
    self.textbox.setText(self.filename)
    self.pharma = []
    self.emit(SIGNAL("changed()"))
    
  def settingsWidget(self):
    self.widget = QWidget()

    layout = QVBoxLayout(self.widget)

    self.textbox = QLineEdit(self.filename)
    self.textbox.setReadOnly(True)
    layout.addWidget(self.textbox)
    
    buttons = QHBoxLayout(self.widget)
    layout.addLayout(buttons)

    self.browsebtn = QPushButton("Browse")
    self.clearbtn = QPushButton("Clear")
    buttons.addWidget(self.browsebtn)
    buttons.addWidget(self.clearbtn)
    
    QObject.connect(self.browsebtn, SIGNAL("clicked()"), self, SLOT("browse()"))
    QObject.connect(self.clearbtn, SIGNAL("clicked()"), self, SLOT("clear()"))

    return self.widget

  def renderOpaque(self, pd):
    self.debug("Len of pharma: %d" % len(self.pharma))
    if len(self.pharma) == 0:
      return
    
    # Painter 
    painter = pd.painter
    # Molecule
    molecule = pd.molecule

    # Color
    color = pd.colorMap

    for point in self.pharma:
      if point['name'] in colorlookup:
        colvals = [float(x) / 255. for x in colorlookup[point['name']]] + [1.0]
        color.setFromRgba(*colvals)
      else:
        color.setFromRgba(0.3, 0.6, 1.0, 1.0)
      painter.setColor(color)        
      
      r = point['radius']
      self.debug("Radius: %.1f" % r)
      x, y, z = point['x'], point['y'], point['z']
      painter.drawSphere(array([x, y, z]), r)

8 comments:

  1. Nice engine, Noel! One comment:

    "If nothing appears on the screen, try rotating the display slightly."

    Calling molecule.update() at the end of the renderOpaque method should fix that, although it seems like something the rendering engine should take care of automatically.

    ReplyDelete
  2. The changed() SIGNAL is not being detected by Avogadro. I've filed a bug.

    ReplyDelete
  3. Another problem (which I described on the Avo mailing list) is that it's not possible to draw transparent spheres from Python. If only the Python bindings were made 1st class citizens of Avo (and if the docs were tidied up), then I think you would see a lot more contributions here.

    ReplyDelete
  4. Please note also the third open source, pharmacophore software: Open3DQSAR: "An open-source software aimed at high-throughput
    chemometric analysis of molecular interaction fields" http://open3dqsar.sourceforge.net/

    ReplyDelete
  5. @filips: That's an interesting suggestion. Traditionally the 3D QSAR (MIF) field is seen as separate from pharmacophores, but there certainly is some overlap.

    ReplyDelete
  6. Hi, trying to figure out engines with Avogadro - for an ubuntu build where would I put this script in order for it to appear as an engine from within Avogadro?

    Cheers

    Clyde

    ReplyDelete
  7. I don't know, and the info doesn't seem to be on the Avo wiki. You'll haae to email the list to find out...

    ReplyDelete
  8. Figured it out by looking at what the package manager was doing when it installed avogadro in case it's of interest it's:

    /usr/share/libavogadro/engineScripts

    Clyde

    ReplyDelete