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