1
2
3
4
5
6 """Half-sphere exposure and coordination number calculation."""
7
8 import warnings
9 from math import pi
10
11 from Bio.PDB.AbstractPropertyMap import AbstractPropertyMap
12 from Bio.PDB.PDBParser import PDBParser
13 from Bio.PDB.Polypeptide import CaPPBuilder, is_aa
14 from Bio.PDB.Vector import rotaxis
15
16
18 """
19 Abstract class to calculate Half-Sphere Exposure (HSE).
20
21 The HSE can be calculated based on the CA-CB vector, or the pseudo CB-CA
22 vector based on three consecutive CA atoms. This is done by two separate
23 subclasses.
24 """
25 - def __init__(self, model, radius, offset, hse_up_key, hse_down_key,
26 angle_key=None):
27 """
28 @param model: model
29 @type model: L{Model}
30
31 @param radius: HSE radius
32 @type radius: float
33
34 @param offset: number of flanking residues that are ignored in the calculation
35 of the number of neighbors
36 @type offset: int
37
38 @param hse_up_key: key used to store HSEup in the entity.xtra attribute
39 @type hse_up_key: string
40
41 @param hse_down_key: key used to store HSEdown in the entity.xtra attribute
42 @type hse_down_key: string
43
44 @param angle_key: key used to store the angle between CA-CB and CA-pCB in
45 the entity.xtra attribute
46 @type angle_key: string
47 """
48 assert(offset>=0)
49
50 self.ca_cb_list=[]
51 ppb=CaPPBuilder()
52 ppl=ppb.build_peptides(model)
53 hse_map={}
54 hse_list=[]
55 hse_keys=[]
56 for pp1 in ppl:
57 for i in range(0, len(pp1)):
58 if i==0:
59 r1=None
60 else:
61 r1=pp1[i-1]
62 r2=pp1[i]
63 if i==len(pp1)-1:
64 r3=None
65 else:
66 r3=pp1[i+1]
67
68 result=self._get_cb(r1, r2, r3)
69 if result is None:
70
71 continue
72 pcb, angle=result
73 hse_u=0
74 hse_d=0
75 ca2=r2['CA'].get_vector()
76 for pp2 in ppl:
77 for j in range(0, len(pp2)):
78 if pp1 is pp2 and abs(i-j)<=offset:
79
80 continue
81 ro=pp2[j]
82 if not is_aa(ro) or not ro.has_id('CA'):
83 continue
84 cao=ro['CA'].get_vector()
85 d=(cao-ca2)
86 if d.norm()<radius:
87 if d.angle(pcb)<(pi/2):
88 hse_u+=1
89 else:
90 hse_d+=1
91 res_id=r2.get_id()
92 chain_id=r2.get_parent().get_id()
93
94 hse_map[(chain_id, res_id)]=(hse_u, hse_d, angle)
95 hse_list.append((r2, (hse_u, hse_d, angle)))
96 hse_keys.append((chain_id, res_id))
97
98 r2.xtra[hse_up_key]=hse_u
99 r2.xtra[hse_down_key]=hse_d
100 if angle_key:
101 r2.xtra[angle_key]=angle
102 AbstractPropertyMap.__init__(self, hse_map, hse_keys, hse_list)
103
105 """This method is provided by the subclasses to calculate HSE."""
106 return NotImplemented
107
109 """
110 Return a pseudo CB vector for a Gly residue.
111 The pseudoCB vector is centered at the origin.
112
113 CB coord=N coord rotated over -120 degrees
114 along the CA-C axis.
115 """
116 try:
117 n_v=residue["N"].get_vector()
118 c_v=residue["C"].get_vector()
119 ca_v=residue["CA"].get_vector()
120 except:
121 return None
122
123 n_v=n_v-ca_v
124 c_v=c_v-ca_v
125
126 rot=rotaxis(-pi*120.0/180.0, c_v)
127 cb_at_origin_v=n_v.left_multiply(rot)
128
129 cb_v=cb_at_origin_v+ca_v
130
131 self.ca_cb_list.append((ca_v, cb_v))
132 return cb_at_origin_v
133
134
135
137 """
138 Class to calculate HSE based on the approximate CA-CB vectors,
139 using three consecutive CA positions.
140 """
141 - def __init__(self, model, radius=12, offset=0):
142 """
143 @param model: the model that contains the residues
144 @type model: L{Model}
145
146 @param radius: radius of the sphere (centred at the CA atom)
147 @type radius: float
148
149 @param offset: number of flanking residues that are ignored in the calculation of the number of neighbors
150 @type offset: int
151 """
152 _AbstractHSExposure.__init__(self, model, radius, offset,
153 'EXP_HSE_A_U', 'EXP_HSE_A_D', 'EXP_CB_PCB_ANGLE')
154
156 """
157 Calculate the approximate CA-CB direction for a central
158 CA atom based on the two flanking CA positions, and the angle
159 with the real CA-CB vector.
160
161 The CA-CB vector is centered at the origin.
162
163 @param r1, r2, r3: three consecutive residues
164 @type r1, r2, r3: L{Residue}
165 """
166 if r1 is None or r3 is None:
167 return None
168 try:
169 ca1=r1['CA'].get_vector()
170 ca2=r2['CA'].get_vector()
171 ca3=r3['CA'].get_vector()
172 except:
173 return None
174
175 d1=ca2-ca1
176 d3=ca2-ca3
177 d1.normalize()
178 d3.normalize()
179
180 b=(d1+d3)
181 b.normalize()
182
183 self.ca_cb_list.append((ca2, b+ca2))
184 if r2.has_id('CB'):
185 cb=r2['CB'].get_vector()
186 cb_ca=cb-ca2
187 cb_ca.normalize()
188 angle=cb_ca.angle(b)
189 elif r2.get_resname()=='GLY':
190 cb_ca=self._get_gly_cb_vector(r2)
191 if cb_ca is None:
192 angle=None
193 else:
194 angle=cb_ca.angle(b)
195 else:
196 angle=None
197
198 return b, angle
199
201 """
202 Write a PyMol script that visualizes the pseudo CB-CA directions
203 at the CA coordinates.
204
205 @param filename: the name of the pymol script file
206 @type filename: string
207 """
208 if len(self.ca_cb_list)==0:
209 warnings.warn("Nothing to draw.", RuntimeWarning)
210 return
211 fp=open(filename, "w")
212 fp.write("from pymol.cgo import *\n")
213 fp.write("from pymol import cmd\n")
214 fp.write("obj=[\n")
215 fp.write("BEGIN, LINES,\n")
216 fp.write("COLOR, %.2f, %.2f, %.2f,\n" % (1.0, 1.0, 1.0))
217 for (ca, cb) in self.ca_cb_list:
218 x,y,z=ca.get_array()
219 fp.write("VERTEX, %.2f, %.2f, %.2f,\n" % (x,y,z))
220 x,y,z=cb.get_array()
221 fp.write("VERTEX, %.2f, %.2f, %.2f,\n" % (x,y,z))
222 fp.write("END]\n")
223 fp.write("cmd.load_cgo(obj, 'HS')\n")
224 fp.close()
225
226
228 """
229 Class to calculate HSE based on the real CA-CB vectors.
230 """
231 - def __init__(self, model, radius=12, offset=0):
232 """
233 @param model: the model that contains the residues
234 @type model: L{Model}
235
236 @param radius: radius of the sphere (centred at the CA atom)
237 @type radius: float
238
239 @param offset: number of flanking residues that are ignored in the calculation of the number of neighbors
240 @type offset: int
241 """
242 _AbstractHSExposure.__init__(self, model, radius, offset,
243 'EXP_HSE_B_U', 'EXP_HSE_B_D')
244
246 """
247 Method to calculate CB-CA vector.
248
249 @param r1, r2, r3: three consecutive residues (only r2 is used)
250 @type r1, r2, r3: L{Residue}
251 """
252 if r2.get_resname()=='GLY':
253 return self._get_gly_cb_vector(r2), 0.0
254 else:
255 if r2.has_id('CB') and r2.has_id('CA'):
256 vcb=r2['CB'].get_vector()
257 vca=r2['CA'].get_vector()
258 return (vcb-vca), 0.0
259 return None
260
261
263 - def __init__(self, model, radius=12.0, offset=0):
264 """
265 A residue's exposure is defined as the number of CA atoms around
266 that residues CA atom. A dictionary is returned that uses a L{Residue}
267 object as key, and the residue exposure as corresponding value.
268
269 @param model: the model that contains the residues
270 @type model: L{Model}
271
272 @param radius: radius of the sphere (centred at the CA atom)
273 @type radius: float
274
275 @param offset: number of flanking residues that are ignored in the calculation of the number of neighbors
276 @type offset: int
277
278 """
279 assert(offset>=0)
280 ppb=CaPPBuilder()
281 ppl=ppb.build_peptides(model)
282 fs_map={}
283 fs_list=[]
284 fs_keys=[]
285 for pp1 in ppl:
286 for i in range(0, len(pp1)):
287 fs=0
288 r1=pp1[i]
289 if not is_aa(r1) or not r1.has_id('CA'):
290 continue
291 ca1=r1['CA']
292 for pp2 in ppl:
293 for j in range(0, len(pp2)):
294 if pp1 is pp2 and abs(i-j)<=offset:
295 continue
296 r2=pp2[j]
297 if not is_aa(r2) or not r2.has_id('CA'):
298 continue
299 ca2=r2['CA']
300 d=(ca2-ca1)
301 if d<radius:
302 fs+=1
303 res_id=r1.get_id()
304 chain_id=r1.get_parent().get_id()
305
306 fs_map[(chain_id, res_id)]=fs
307 fs_list.append((r1, fs))
308 fs_keys.append((chain_id, res_id))
309
310 r1.xtra['EXP_CN']=fs
311 AbstractPropertyMap.__init__(self, fs_map, fs_keys, fs_list)
312
313
314 if __name__=="__main__":
315
316 import sys
317
318 p=PDBParser()
319 s=p.get_structure('X', sys.argv[1])
320 model=s[0]
321
322
323 RADIUS=13.0
324 OFFSET=0
325
326 hse=HSExposureCA(model, radius=RADIUS, offset=OFFSET)
327 for l in hse:
328 print l
329 print
330
331 hse=HSExposureCB(model, radius=RADIUS, offset=OFFSET)
332 for l in hse:
333 print l
334 print
335
336 hse=ExposureCN(model, radius=RADIUS, offset=OFFSET)
337 for l in hse:
338 print l
339 print
340
341 for c in model:
342 for r in c:
343 try:
344 print r.xtra['PCB_CB_ANGLE']
345 except:
346 pass
347