1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
# make yield compatible with Python 2.2
from __future__ import generators
from Numeric import array, sum, sqrt
import tempfile
import os
import sys
from Bio.PDB import *
from AbstractPropertyMap import AbstractPropertyMap
__doc__="""
Calculation of residue depth (using Michel Sanner's MSMS program for the
surface calculation).
Residue depth is the average distance of the atoms of a residue from
the solvent accessible surface.
Residue Depth:
rd=ResidueDepth(model, pdb_file)
print rd[(chain_id, res_id)]
Direct MSMS interface:
Typical use:
surface=get_surface("1FAT.pdb")
Surface is a Numpy array with all the surface
vertices.
Distance to surface:
dist=min_dist(coord, surface)
where coord is the coord of an atom within the volume
bound by the surface (ie. atom depth).
To calculate the residue depth (average atom depth
of the atoms in a residue):
rd=residue_depth(residue, surface)
"""
def _read_vertex_array(filename):
"""
Read the vertex list into a Numpy array.
"""
fp=open(filename, "r")
vertex_list=[]
for l in fp.readlines():
sl=l.split()
if not len(sl)==9:
# skip header
continue
vl=map(float, sl[0:3])
vertex_list.append(vl)
fp.close()
return array(vertex_list)
def get_surface(pdb_file, PDB_TO_XYZR="pdb_to_xyzr", MSMS="msms"):
"""
Return a Numpy array that represents
the vertex list of the molecular surface.
PDB_TO_XYZR --- pdb_to_xyzr executable (arg. to os.system)
MSMS --- msms executable (arg. to os.system)
"""
# extract xyz and set radii
xyz_tmp=tempfile.mktemp()
PDB_TO_XYZR=PDB_TO_XYZR+" %s > %s"
make_xyz=PDB_TO_XYZR % (pdb_file, xyz_tmp)
os.system(make_xyz)
# make surface
surface_tmp=tempfile.mktemp()
MSMS=MSMS+" -probe_radius 1.5 -if %s -of %s > "+tempfile.mktemp()
make_surface=MSMS % (xyz_tmp, surface_tmp)
os.system(make_surface)
surface_file=surface_tmp+".vert"
# read surface vertices from vertex file
surface=_read_vertex_array(surface_file)
# clean up tmp files
# ...this is dangerous
#os.system("rm "+xyz_tmp)
#os.system("rm "+surface_tmp+".vert")
#os.system("rm "+surface_tmp+".face")
return surface
def min_dist(coord, surface):
"""
Return minimum distance between coord
and surface.
"""
d=surface-coord
d2=sum(d*d, 1)
return sqrt(min(d2))
def residue_depth(residue, surface):
"""
Return average distance to surface for all
atoms in a residue, ie. the residue depth.
"""
atom_list=residue.get_unpacked_list()
length=len(atom_list)
d=0
for atom in atom_list:
coord=atom.get_coord()
d=d+min_dist(coord, surface)
return d/length
def ca_depth(residue, surface):
if not residue.has_id("CA"):
return None
ca=residue["CA"]
coord=ca.get_coord()
return min_dist(coord, surface)
class ResidueDepth(AbstractPropertyMap):
"""
Calculate residue and CA depth for all residues.
"""
def __init__(self, model, pdb_file):
depth_dict={}
depth_list=[]
depth_keys=[]
# get_residue
residue_list=Selection.unfold_entities(model, 'R')
# make surface from PDB file
surface=get_surface(pdb_file)
# calculate rdepth for each residue
for residue in residue_list:
if not is_aa(residue):
continue
rd=residue_depth(residue, surface)
ca_rd=ca_depth(residue, surface)
# Get the key
res_id=residue.get_id()
chain_id=residue.get_parent().get_id()
depth_dict[(chain_id, res_id)]=(rd, ca_rd)
depth_list.append((residue, (rd, ca_rd)))
depth_keys.append((chain_id, res_id))
# Update xtra information
residue.xtra['EXP_RD']=rd
residue.xtra['EXP_RD_CA']=ca_rd
AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list)
if __name__=="__main__":
import sys
p=PDBParser()
s=p.get_structure("X", sys.argv[1])
model=s[0]
rd=ResidueDepth(model, sys.argv[1])
for item in rd:
print item
|
