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