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