[Chimera-users] chimera and hydrophobic potentials

[Thomas Goddard]

Hi Miguel,

  I wrote a conversion program cns2delphi.py that takes the k2f output
and produces a delphi file and tried it on the data you sent.  The
delphi file is read properly by Chimera DelphiViewer and can be
used to color an MSMS surface.

  A problem with the data you sent is that the hydrophobicity map is
zero over the entire MSMS surface.  There are blobs of hydrophobicity
a few angstroms away from the surface but they don't touch the
surface.  I observed this by showing the map with the Chimera volume
viewer (just use File/Open...).  You can change the probe radius for
generating the MSMS surface by going to Favorites/Model Panel,
selecting the MSMS surface, clicking the Attributes button, changing
the probe radius setting, and pressing Enter.  (You need to set the
surface color to a solid color before changing the radius or an error
message is generated.)  Unfortunately that does not push the molecular
surface out uniformly -- only the surface in the cavities is pushed
further from the molecule.  Apparently the probe radius is the radius
of a ball rolled over the molecule with standard atom radii.  The MSMS
surface is the surface of the region that cannot be touched by the
ball.  So a bigger ball does not change the surface location at
protrusions.  Maybe a better approach is to specify larger atom radii.

  Below is the cns2delphi.py Python program.  It requires the Numeric
Python module which is not a standard part of Python.  The Python that
comes with Chimera includes Numeric.  You can convert k2f output to
delphi format as follows.

  % /usr/local/chimera/bin/python2.3 cns2delphi.py 1ltlA.k2f 1ltlA.phi

This creates the 1ltlA.phi file.

  To use the Chimera DelphiViewer extension, start Chimera, open
1ltlA.pdb, make a molecular surface with Actions/Surface/Show, start
Extensions/Delphi/DelphiViewer, enter Potential File 1ltlA.phi, click
the Browse button next to MSMS Surface in the Molecular surface pane
so the surface is listed, then press Apply in the Molecular surface
pane.  This will color the surface.  You can also just use
File/Open... on 1ltlA.phi to see an isosurface of the hydrophobicity
map.  You can drag the vertical bar on the histogram shown in the
volume viewer dialog to change the isosurface threshold level.

       Tom

-----
cns2delphi.py follows.  Preserve indentation since it is important in Python.

#!/usr/bin/env python
# ----------------------------------------------------------------------------
# Convert CNS volume to Delphi format.
#

# ----------------------------------------------------------------------------
# Returns matrix, xyz_origin, and xyz_step
#
def read_cns_file(path):

    f = open(path, 'r')

    f.readline()                          # First line is blank
    ntitle_line = f.readline()            # integer number of comment lines
    try:
      ntitle = int(ntitle_line.split()[0])
    except ValueError:
      raise SyntaxError, ('Invalid CNS comment line count on line 2: %s'
                          % ntitle_line[:80])

    for t in range(ntitle):
      f.readline()

    extent_line = f.readline()
    extent = split_fields(extent_line, 8, 9)
    try:
      na, amin, amax, nb, bmin, bmax, nc, cmin, cmax = map(int, extent)
    except ValueError:
      raise SyntaxError, 'Invalid CNS grid size line: %s' % extent_line[:80]
        
    cell_line = f.readline()
    cell = split_fields(cell_line, 12, 6)
    try:
      cell_params = map(float, cell)
    except ValueError:
      raise SyntaxError, ('Invalid CNS cell parameters line: %s'
                          % cell_line[:80])
    cell_size = cell_params[0:3]
    cell_angles = cell_params[3:6]

    axis_order_line = f.readline()        # Must be ZYX

    asize, bsize, csize = (amax - amin + 1, bmax - bmin + 1, cmax - cmin + 1)
    if asize <= 0 or bsize <= 0 or csize <=0:
      raise SyntaxError, ('Bad CNS grid size (%d,%d,%d)'
                          % (asize, bsize, csize))

    import Numeric
    data = Numeric.zeros((csize, bsize, asize), Numeric.Float32)
    for c in range(csize):
      read_floats(f, 1)                        # c section number
      values = read_floats(f, asize * bsize)
      array = Numeric.array(values, Numeric.Float32)
      plane = Numeric.reshape(array, (bsize, asize))
      data[c,:,:] = plane

    f.readline()                          # footer - int value must be -9999
    f.readline()                          # density average and standard dev

    f.close()

    xyz_step = (cell_size[0]/na, cell_size[1]/nb, cell_size[2]/nc)
    xyz_origin = map(lambda gmin, step: gmin * step,
                     (amin, bmin, cmin), xyz_step)

    return data, xyz_origin, xyz_step
  
# -----------------------------------------------------------------------------
# Read ascii float values on as many lines as needed to get count values.
#
def read_floats(f, count):

  values = [0.0] * count
  c = 0
  while c < count:
    line = f.readline()
    if line == '':
      raise FileFormatError('File is shorter than expected')
    fields = split_fields(line, 12, 6)
    if c + len(fields) > count:
      fields = fields[:count-c]
    for v in map(float, fields):
      values[c] = v
      c = c + 1

  return values
  
# -----------------------------------------------------------------------------
#
def split_fields(line, field_size, max_fields):

  fields = []
  for k in range(0, len(line), field_size):
    f = line[k:k+field_size].strip()
    if f:
      fields.append(f)
    else:
      break
  return fields[:max_fields]


# ----------------------------------------------------------------------------
#
def write_delphi_file(path, matrix, xyz_center, scale):

    zsize, ysize, xsize = matrix.shape
    if xsize != ysize or xsize != zsize:
      raise ValueError, "Matrix size (%d,%d,%d) has to be the same in all 3 dimensions" % (zsize, ysize, xsize)
    size = xsize

    
    uplbl = delphi_record('Converted from ' + path)
    morelabels = delphi_record('%70s' % '')
    import Numeric
    data_flat = Numeric.ravel(matrix)
    data = delphi_record(floats_as_string(data_flat))
    botlbl = delphi_record('')
    params = delphi_record(floats_as_string((scale,) + tuple(xyz_center)) +
                           integer_as_string(size))

    f = open(path, 'wb')
    f.write(uplbl + morelabels + data + botlbl + params)
    f.close()

# ----------------------------------------------------------------------------
#
def delphi_record(data):

    size_str = integer_as_string(len(data))
    r = size_str + data + size_str
    return r

# ----------------------------------------------------------------------------
#
def integer_as_string(i):

    import Numeric
    return Numeric.array((i,), Numeric.Int32).tostring()

# ----------------------------------------------------------------------------
#
def floats_as_string(f):

    import Numeric
    return Numeric.array(f, Numeric.Float32).tostring()

# ----------------------------------------------------------------------------
#
def cns_to_delphi_file(cns_path, delphi_path):

    data, xyz_origin, xyz_step = read_cns_file(cns_path)
    if xyz_step[0] != xyz_step[1] or xyz_step[0] != xyz_step[2]:
      raise Value_Error, 'Step sizes (%.2g, %.2g, %.2g) are not the same along all axes' % xyz_step
    scale = 1.0 / xyz_step[0]

    size = list(data.shape)
    size.reverse()
    xyz_center = map(lambda o, st, sz: o + st * ((sz-1)/2),
                     xyz_origin, xyz_step, size)

    write_delphi_file(delphi_path, data, xyz_center, scale)

# ----------------------------------------------------------------------------
#
if __name__ == '__main__':

  import sys
  if len(sys.argv) != 3:
    sys.stderr.write('Syntax: cns2delphi <cnsfile> <delphifile>\n')
    sys.exit(1)

  cns_path = sys.argv[1]
  delphi_path = sys.argv[2]
  cns_to_delphi_file(cns_path, delphi_path)