Source Code for Module Bio.PDB.ResidueDepth'

  1  # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk) 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  5   
  6  """Calculation of residue depth using command line tool MSMS. 
  7   
  8  This module uses Michel Sanner's MSMS program for the surface calculation 
  9  (specifically commands msms and pdb_to_xyzr). See: 
 10  http://mgltools.scripps.edu/packages/MSMS 
 11   
 12  Residue depth is the average distance of the atoms of a residue from 
 13  the solvent accessible surface. 
 14   
 15  Residue Depth: 
 16   
 17      >>> rd = ResidueDepth(model, pdb_file) 
 18      >>> print rd[(chain_id, res_id)] 
 19   
 20  Direct MSMS interface: 
 21   
 22      Typical use: 
 23   
 24          >>> surface = get_surface("1FAT.pdb") 
 25   
 26      Surface is a Numeric array with all the surface 
 27      vertices. 
 28   
 29      Distance to surface: 
 30   
 31          >>> dist = min_dist(coord, surface) 
 32   
 33      where coord is the coord of an atom within the volume 
 34      bound by the surface (ie. atom depth). 
 35   
 36      To calculate the residue depth (average atom depth 
 37      of the atoms in a residue): 
 38   
 39          >>> rd = residue_depth(residue, surface) 
 40  """ 
 41   
 42  import os 
 43  import tempfile 
 44   
 45  import numpy 
 46   
 47  from Bio.PDB import Selection 
 48  from Bio.PDB.AbstractPropertyMap import AbstractPropertyMap 
 49  from Bio.PDB.Polypeptide import is_aa 
 50   
 51   
 52 -def _read_vertex_array(filename): 
 53      """ 
 54      Read the vertex list into a Numeric array. 
 55      """ 
 56      fp=open(filename, "r") 
 57      vertex_list=[] 
 58      for l in fp.readlines(): 
 59          sl=l.split() 
 60          if not len(sl)==9: 
 61              # skip header 
 62              continue 
 63          vl=map(float, sl[0:3]) 
 64          vertex_list.append(vl) 
 65      fp.close() 
 66      return numpy.array(vertex_list) 
 67   
 68   
 69 -def get_surface(pdb_file, PDB_TO_XYZR="pdb_to_xyzr", MSMS="msms"): 
 70      """ 
 71      Return a Numeric array that represents 
 72      the vertex list of the molecular surface. 
 73   
 74      PDB_TO_XYZR --- pdb_to_xyzr executable (arg. to os.system) 
 75      MSMS --- msms executable (arg. to os.system) 
 76      """ 
 77      # extract xyz and set radii 
 78      xyz_tmp=tempfile.mktemp() 
 79      PDB_TO_XYZR=PDB_TO_XYZR+" %s > %s" 
 80      make_xyz=PDB_TO_XYZR % (pdb_file, xyz_tmp) 
 81      os.system(make_xyz) 
 82      assert os.path.isfile(xyz_tmp), \ 
 83          "Failed to generate XYZR file using command:\n%s" % make_xyz 
 84      # make surface 
 85      surface_tmp=tempfile.mktemp() 
 86      MSMS=MSMS+" -probe_radius 1.5 -if %s -of %s > "+tempfile.mktemp() 
 87      make_surface=MSMS % (xyz_tmp, surface_tmp) 
 88      os.system(make_surface) 
 89      surface_file=surface_tmp+".vert" 
 90      assert os.path.isfile(surface_file), \ 
 91          "Failed to generate surface file using command:\n%s" % make_surface 
 92      # read surface vertices from vertex file 
 93      surface=_read_vertex_array(surface_file) 
 94      # clean up tmp files 
 95      # ...this is dangerous 
 96      #os.system("rm "+xyz_tmp) 
 97      #os.system("rm "+surface_tmp+".vert") 
 98      #os.system("rm "+surface_tmp+".face") 
 99      return surface 
100   
101   
102 -def min_dist(coord, surface): 
103      """ 
104      Return minimum distance between coord 
105      and surface. 
106      """ 
107      d=surface-coord 
108      d2=numpy.sum(d*d, 1) 
109      return numpy.sqrt(min(d2)) 
110   
111   
112 -def residue_depth(residue, surface): 
113      """ 
114      Return average distance to surface for all 
115      atoms in a residue, ie. the residue depth. 
116      """ 
117      atom_list=residue.get_unpacked_list() 
118      length=len(atom_list) 
119      d=0 
120      for atom in atom_list: 
121          coord=atom.get_coord() 
122          d=d+min_dist(coord, surface) 
123      return d/length 
124   
125   
126 -def ca_depth(residue, surface): 
127      if not residue.has_id("CA"): 
128          return None 
129      ca=residue["CA"] 
130      coord=ca.get_coord() 
131      return min_dist(coord, surface) 
132   
133   
134 -class ResidueDepth(AbstractPropertyMap): 
135      """ 
136      Calculate residue and CA depth for all residues. 
137      """ 
138 -    def __init__(self, model, pdb_file): 
139          depth_dict={} 
140          depth_list=[] 
141          depth_keys=[] 
142          # get_residue 
143          residue_list=Selection.unfold_entities(model, 'R') 
144          # make surface from PDB file 
145          surface=get_surface(pdb_file) 
146          # calculate rdepth for each residue 
147          for residue in residue_list: 
148              if not is_aa(residue): 
149                  continue 
150              rd=residue_depth(residue, surface) 
151              ca_rd=ca_depth(residue, surface) 
152              # Get the key 
153              res_id=residue.get_id() 
154              chain_id=residue.get_parent().get_id() 
155              depth_dict[(chain_id, res_id)]=(rd, ca_rd) 
156              depth_list.append((residue, (rd, ca_rd))) 
157              depth_keys.append((chain_id, res_id)) 
158              # Update xtra information 
159              residue.xtra['EXP_RD']=rd 
160              residue.xtra['EXP_RD_CA']=ca_rd 
161          AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list) 
162   
163   
164  if __name__=="__main__": 
165   
166      import sys 
167      from Bio.PDB import PDBParser 
168   
169      p=PDBParser() 
170      s=p.get_structure("X", sys.argv[1]) 
171      model=s[0] 
172   
173      rd=ResidueDepth(model, sys.argv[1]) 
174   
175      for item in rd: 
176          print item 
177