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, typical use: 
 21   
 22      >>> surface = get_surface("1FAT.pdb") 
 23   
 24  The surface is a Numeric array with all the surface vertices. 
 25   
 26  Distance to surface: 
 27   
 28      >>> dist = min_dist(coord, surface) 
 29   
 30  where coord is the coord of an atom within the volume bound by 
 31  the surface (ie. atom depth). 
 32   
 33  To calculate the residue depth (average atom depth of the atoms 
 34  in a residue): 
 35   
 36      >>> rd = residue_depth(residue, surface) 
 37   
 38  """ 
 39   
 40  from __future__ import print_function 
 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 with open(filename, "r") as fp: 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 = [float(x) for x in sl[0:3]] 64 vertex_list.append(vl) 65 return numpy.array(vertex_list)
66 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
111 -def residue_depth(residue, surface):
112 """ 113 Return average distance to surface for all 114 atoms in a residue, ie. the residue depth. 115 """ 116 atom_list=residue.get_unpacked_list() 117 length=len(atom_list) 118 d=0 119 for atom in atom_list: 120 coord=atom.get_coord() 121 d=d+min_dist(coord, surface) 122 return d/length
123 124
125 -def ca_depth(residue, surface):
126 if not residue.has_id("CA"): 127 return None 128 ca=residue["CA"] 129 coord=ca.get_coord() 130 return min_dist(coord, surface)
131 132
133 -class ResidueDepth(AbstractPropertyMap):
134 """ 135 Calculate residue and CA depth for all residues. 136 """
137 - def __init__(self, model, pdb_file):
138 depth_dict={} 139 depth_list=[] 140 depth_keys=[] 141 # get_residue 142 residue_list=Selection.unfold_entities(model, 'R') 143 # make surface from PDB file 144 surface=get_surface(pdb_file) 145 # calculate rdepth for each residue 146 for residue in residue_list: 147 if not is_aa(residue): 148 continue 149 rd=residue_depth(residue, surface) 150 ca_rd=ca_depth(residue, surface) 151 # Get the key 152 res_id=residue.get_id() 153 chain_id=residue.get_parent().get_id() 154 depth_dict[(chain_id, res_id)]=(rd, ca_rd) 155 depth_list.append((residue, (rd, ca_rd))) 156 depth_keys.append((chain_id, res_id)) 157 # Update xtra information 158 residue.xtra['EXP_RD']=rd 159 residue.xtra['EXP_RD_CA']=ca_rd 160 AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list)
161 162 163 if __name__=="__main__": 164 165 import sys 166 from Bio.PDB import PDBParser 167 168 p=PDBParser() 169 s=p.get_structure("X", sys.argv[1]) 170 model=s[0] 171 172 rd=ResidueDepth(model, sys.argv[1]) 173 174 for item in rd: 175 print(item) 176