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