1
2
3
4
5
6
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
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
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
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
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
88 surface=_read_vertex_array(surface_file)
89
90
91
92
93
94 return surface
95
96
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
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
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
127 """
128 Calculate residue and CA depth for all residues.
129 """
131 depth_dict={}
132 depth_list=[]
133 depth_keys=[]
134
135 residue_list=Selection.unfold_entities(model, 'R')
136
137 surface=get_surface(pdb_file)
138
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
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
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