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