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