1
2
3
4
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
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
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
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
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
85 surface=_read_vertex_array(surface_file)
86
87
88
89
90
91 return surface
92
93
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
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
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
124 """
125 Calculate residue and CA depth for all residues.
126 """
128 depth_dict={}
129 depth_list=[]
130 depth_keys=[]
131
132 residue_list=Selection.unfold_entities(model, 'R')
133
134 surface=get_surface(pdb_file)
135
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
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
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