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  from __future__ import print_function 
 43   
 44  import os 
 45  import tempfile 
 46   
 47  import numpy 
 48   
 49  from Bio.PDB import Selection 
 50  from Bio.PDB.AbstractPropertyMap import AbstractPropertyMap 
 51  from Bio.PDB.Polypeptide import is_aa 
 52   
 53   
54 -def _read_vertex_array(filename):
55 """ 56 Read the vertex list into a Numeric array. 57 """ 58 with open(filename, "r") as fp: 59 vertex_list=[] 60 for l in fp.readlines(): 61 sl=l.split() 62 if not len(sl)==9: 63 # skip header 64 continue 65 vl = [float(x) for x in sl[0:3]] 66 vertex_list.append(vl) 67 return numpy.array(vertex_list)
68 69
70 -def get_surface(pdb_file, PDB_TO_XYZR="pdb_to_xyzr", MSMS="msms"):
71 """ 72 Return a Numeric array that represents 73 the vertex list of the molecular surface. 74 75 PDB_TO_XYZR --- pdb_to_xyzr executable (arg. to os.system) 76 MSMS --- msms executable (arg. to os.system) 77 """ 78 # extract xyz and set radii 79 xyz_tmp=tempfile.mktemp() 80 PDB_TO_XYZR=PDB_TO_XYZR+" %s > %s" 81 make_xyz=PDB_TO_XYZR % (pdb_file, xyz_tmp) 82 os.system(make_xyz) 83 assert os.path.isfile(xyz_tmp), \ 84 "Failed to generate XYZR file using command:\n%s" % make_xyz 85 # make surface 86 surface_tmp=tempfile.mktemp() 87 MSMS=MSMS+" -probe_radius 1.5 -if %s -of %s > "+tempfile.mktemp() 88 make_surface=MSMS % (xyz_tmp, surface_tmp) 89 os.system(make_surface) 90 surface_file=surface_tmp+".vert" 91 assert os.path.isfile(surface_file), \ 92 "Failed to generate surface file using command:\n%s" % make_surface 93 # read surface vertices from vertex file 94 surface=_read_vertex_array(surface_file) 95 # clean up tmp files 96 # ...this is dangerous 97 #os.system("rm "+xyz_tmp) 98 #os.system("rm "+surface_tmp+".vert") 99 #os.system("rm "+surface_tmp+".face") 100 return surface
101 102
103 -def min_dist(coord, surface):
104 """ 105 Return minimum distance between coord 106 and surface. 107 """ 108 d=surface-coord 109 d2=numpy.sum(d*d, 1) 110 return numpy.sqrt(min(d2))
111 112
113 -def residue_depth(residue, surface):
114 """ 115 Return average distance to surface for all 116 atoms in a residue, ie. the residue depth. 117 """ 118 atom_list=residue.get_unpacked_list() 119 length=len(atom_list) 120 d=0 121 for atom in atom_list: 122 coord=atom.get_coord() 123 d=d+min_dist(coord, surface) 124 return d/length
125 126
127 -def ca_depth(residue, surface):
128 if not residue.has_id("CA"): 129 return None 130 ca=residue["CA"] 131 coord=ca.get_coord() 132 return min_dist(coord, surface)
133 134
135 -class ResidueDepth(AbstractPropertyMap):
136 """ 137 Calculate residue and CA depth for all residues. 138 """
139 - def __init__(self, model, pdb_file):
140 depth_dict={} 141 depth_list=[] 142 depth_keys=[] 143 # get_residue 144 residue_list=Selection.unfold_entities(model, 'R') 145 # make surface from PDB file 146 surface=get_surface(pdb_file) 147 # calculate rdepth for each residue 148 for residue in residue_list: 149 if not is_aa(residue): 150 continue 151 rd=residue_depth(residue, surface) 152 ca_rd=ca_depth(residue, surface) 153 # Get the key 154 res_id=residue.get_id() 155 chain_id=residue.get_parent().get_id() 156 depth_dict[(chain_id, res_id)]=(rd, ca_rd) 157 depth_list.append((residue, (rd, ca_rd))) 158 depth_keys.append((chain_id, res_id)) 159 # Update xtra information 160 residue.xtra['EXP_RD']=rd 161 residue.xtra['EXP_RD_CA']=ca_rd 162 AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list)
163 164 165 if __name__=="__main__": 166 167 import sys 168 from Bio.PDB import PDBParser 169 170 p=PDBParser() 171 s=p.get_structure("X", sys.argv[1]) 172 model=s[0] 173 174 rd=ResidueDepth(model, sys.argv[1]) 175 176 for item in rd: 177 print(item) 178