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