Package Bio :: Package PDB :: Module ResidueDepth'
[hide private]
[frames] | no frames]

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  # make yield compatible with Python 2.2 
  7  from __future__ import generators 
  8   
  9  from Numeric import array, sum, sqrt 
 10  import tempfile 
 11  import os 
 12  import sys 
 13   
 14  from Bio.PDB import * 
 15  from AbstractPropertyMap import AbstractPropertyMap 
 16   
 17  __doc__=""" 
 18  Calculation of residue depth (using Michel Sanner's MSMS program for the 
 19  surface calculation). 
 20   
 21  Residue depth is the average distance of the atoms of a residue from  
 22  the solvent accessible surface. 
 23   
 24  Residue Depth: 
 25   
 26      rd=ResidueDepth(model, pdb_file) 
 27   
 28      print rd[(chain_id, res_id)] 
 29   
 30  Direct MSMS interface: 
 31   
 32      Typical use: 
 33   
 34          surface=get_surface("1FAT.pdb") 
 35   
 36      Surface is a Numeric array with all the surface  
 37      vertices.   
 38   
 39      Distance to surface: 
 40   
 41          dist=min_dist(coord, surface) 
 42   
 43      where coord is the coord of an atom within the volume 
 44      bound by the surface (ie. atom depth). 
 45   
 46      To calculate the residue depth (average atom depth 
 47      of the atoms in a residue): 
 48   
 49      rd=residue_depth(residue, surface) 
 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 array(vertex_list)
67
68 -def get_surface(pdb_file, PDB_TO_XYZR="pdb_to_xyzr", MSMS="msms"):
69 """ 70 Return a Numeric array that represents 71 the vertex list of the molecular surface. 72 73 PDB_TO_XYZR --- pdb_to_xyzr executable (arg. to os.system) 74 MSMS --- msms executable (arg. to os.system) 75 """ 76 # extract xyz and set radii 77 xyz_tmp=tempfile.mktemp() 78 PDB_TO_XYZR=PDB_TO_XYZR+" %s > %s" 79 make_xyz=PDB_TO_XYZR % (pdb_file, xyz_tmp) 80 os.system(make_xyz) 81 # make surface 82 surface_tmp=tempfile.mktemp() 83 MSMS=MSMS+" -probe_radius 1.5 -if %s -of %s > "+tempfile.mktemp() 84 make_surface=MSMS % (xyz_tmp, surface_tmp) 85 os.system(make_surface) 86 surface_file=surface_tmp+".vert" 87 # read surface vertices from vertex file 88 surface=_read_vertex_array(surface_file) 89 # clean up tmp files 90 # ...this is dangerous 91 #os.system("rm "+xyz_tmp) 92 #os.system("rm "+surface_tmp+".vert") 93 #os.system("rm "+surface_tmp+".face") 94 return surface
95 96
97 -def min_dist(coord, surface):
98 """ 99 Return minimum distance between coord 100 and surface. 101 """ 102 d=surface-coord 103 d2=sum(d*d, 1) 104 return sqrt(min(d2))
105
106 -def residue_depth(residue, surface):
107 """ 108 Return average distance to surface for all 109 atoms in a residue, ie. the residue depth. 110 """ 111 atom_list=residue.get_unpacked_list() 112 length=len(atom_list) 113 d=0 114 for atom in atom_list: 115 coord=atom.get_coord() 116 d=d+min_dist(coord, surface) 117 return d/length
118
119 -def ca_depth(residue, surface):
120 if not residue.has_id("CA"): 121 return None 122 ca=residue["CA"] 123 coord=ca.get_coord() 124 return min_dist(coord, surface)
125
126 -class ResidueDepth(AbstractPropertyMap):
127 """ 128 Calculate residue and CA depth for all residues. 129 """
130 - def __init__(self, model, pdb_file):
131 depth_dict={} 132 depth_list=[] 133 depth_keys=[] 134 # get_residue 135 residue_list=Selection.unfold_entities(model, 'R') 136 # make surface from PDB file 137 surface=get_surface(pdb_file) 138 # calculate rdepth for each residue 139 for residue in residue_list: 140 if not is_aa(residue): 141 continue 142 rd=residue_depth(residue, surface) 143 ca_rd=ca_depth(residue, surface) 144 # Get the key 145 res_id=residue.get_id() 146 chain_id=residue.get_parent().get_id() 147 depth_dict[(chain_id, res_id)]=(rd, ca_rd) 148 depth_list.append((residue, (rd, ca_rd))) 149 depth_keys.append((chain_id, res_id)) 150 # Update xtra information 151 residue.xtra['EXP_RD']=rd 152 residue.xtra['EXP_RD_CA']=ca_rd 153 AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list)
154 155 156 if __name__=="__main__": 157 158 import sys 159 160 p=PDBParser() 161 s=p.get_structure("X", sys.argv[1]) 162 model=s[0] 163 164 rd=ResidueDepth(model, sys.argv[1]) 165 166 167 for item in rd: 168 print item 169