Re: [Chimera-users] chimera and hydrophobic potentials
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)
participants (1)
-
Thomas Goddard