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  """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 -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 assert os.path.isfile(xyz_tmp), \ 82 "Failed to generate XYZR file using command:\n%s" % make_xyz 83 # make surface 84 surface_tmp=tempfile.mktemp() 85 MSMS=MSMS+" -probe_radius 1.5 -if %s -of %s > "+tempfile.mktemp() 86 make_surface=MSMS % (xyz_tmp, surface_tmp) 87 os.system(make_surface) 88 surface_file=surface_tmp+".vert" 89 assert os.path.isfile(surface_file), \ 90 "Failed to generate surface file using command:\n%s" % make_surface 91 # read surface vertices from vertex file 92 surface=_read_vertex_array(surface_file) 93 # clean up tmp files 94 # ...this is dangerous 95 #os.system("rm "+xyz_tmp) 96 #os.system("rm "+surface_tmp+".vert") 97 #os.system("rm "+surface_tmp+".face") 98 return surface
99 100
101 -def min_dist(coord, surface):
102 """ 103 Return minimum distance between coord 104 and surface. 105 """ 106 d=surface-coord 107 d2=numpy.sum(d*d, 1) 108 return numpy.sqrt(min(d2))
109
110 -def residue_depth(residue, surface):
111 """ 112 Return average distance to surface for all 113 atoms in a residue, ie. the residue depth. 114 """ 115 atom_list=residue.get_unpacked_list() 116 length=len(atom_list) 117 d=0 118 for atom in atom_list: 119 coord=atom.get_coord() 120 d=d+min_dist(coord, surface) 121 return d/length
122
123 -def ca_depth(residue, surface):
124 if not residue.has_id("CA"): 125 return None 126 ca=residue["CA"] 127 coord=ca.get_coord() 128 return min_dist(coord, surface)
129
130 -class ResidueDepth(AbstractPropertyMap):
131 """ 132 Calculate residue and CA depth for all residues. 133 """
134 - def __init__(self, model, pdb_file):
135 depth_dict={} 136 depth_list=[] 137 depth_keys=[] 138 # get_residue 139 residue_list=Selection.unfold_entities(model, 'R') 140 # make surface from PDB file 141 surface=get_surface(pdb_file) 142 # calculate rdepth for each residue 143 for residue in residue_list: 144 if not is_aa(residue): 145 continue 146 rd=residue_depth(residue, surface) 147 ca_rd=ca_depth(residue, surface) 148 # Get the key 149 res_id=residue.get_id() 150 chain_id=residue.get_parent().get_id() 151 depth_dict[(chain_id, res_id)]=(rd, ca_rd) 152 depth_list.append((residue, (rd, ca_rd))) 153 depth_keys.append((chain_id, res_id)) 154 # Update xtra information 155 residue.xtra['EXP_RD']=rd 156 residue.xtra['EXP_RD_CA']=ca_rd 157 AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list)
158 159 160 if __name__=="__main__": 161 162 import sys 163 from Bio.PDB import PDBParser 164 165 p=PDBParser() 166 s=p.get_structure("X", sys.argv[1]) 167 model=s[0] 168 169 rd=ResidueDepth(model, sys.argv[1]) 170 171 172 for item in rd: 173 print item 174