1
2
3
4
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
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 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
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
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
92 surface=_read_vertex_array(surface_file)
93
94
95
96
97
98 return surface
99
100
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
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
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
131 """
132 Calculate residue and CA depth for all residues.
133 """
135 depth_dict={}
136 depth_list=[]
137 depth_keys=[]
138
139 residue_list=Selection.unfold_entities(model, 'R')
140
141 surface=get_surface(pdb_file)
142
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
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
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