Package Bio :: Package PDB :: Module HSExposure
[hide private]
[frames] | no frames]

Source Code for Module Bio.PDB.HSExposure

  1  # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk) 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  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   
16 -class _AbstractHSExposure(AbstractPropertyMap):
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 # For PyMOL visualization 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 # This method is provided by the subclasses to calculate HSE 67 result=self._get_cb(r1, r2, r3) 68 if result is None: 69 # Missing atoms, or i==0, or i==len(pp1)-1 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 # neighboring residues in the chain are ignored 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 # Fill the 3 data structures 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 # Add to xtra 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
103 - def _get_gly_cb_vector(self, residue):
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 # center at origin 118 n_v=n_v-ca_v 119 c_v=c_v-ca_v 120 # rotation around c-ca over -120 deg 121 rot=rotaxis(-pi*120.0/180.0, c_v) 122 cb_at_origin_v=n_v.left_multiply(rot) 123 # move back to ca position 124 cb_v=cb_at_origin_v+ca_v 125 # This is for PyMol visualization 126 self.ca_cb_list.append((ca_v, cb_v)) 127 return cb_at_origin_v
128 129 130
131 -class HSExposureCA(_AbstractHSExposure):
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
150 - def _get_cb(self, r1, r2, r3):
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 # center 170 d1=ca2-ca1 171 d3=ca2-ca3 172 d1.normalize() 173 d3.normalize() 174 # bisection 175 b=(d1+d3) 176 b.normalize() 177 # Add to ca_cb_list for drawing 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 # vector b is centered at the origin! 193 return b, angle
194
195 - def pcb_vectors_pymol(self, filename="hs_exp.py"):
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
222 -class HSExposureCB(_AbstractHSExposure):
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
240 - def _get_cb(self, r1, r2, r3):
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
257 -class ExposureCN(AbstractPropertyMap):
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 # Fill the 3 data structures 301 fs_map[(chain_id, res_id)]=fs 302 fs_list.append((r1, fs)) 303 fs_keys.append((chain_id, res_id)) 304 # Add to xtra 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 # Neighbor sphere radius 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