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  import warnings 
  9  from math import pi 
 10   
 11  from Bio.PDB.AbstractPropertyMap import AbstractPropertyMap 
 12  from Bio.PDB.PDBParser import PDBParser 
 13  from Bio.PDB.Polypeptide import CaPPBuilder, is_aa 
 14  from Bio.PDB.Vector import rotaxis 
 15   
 16   
17 -class _AbstractHSExposure(AbstractPropertyMap):
18 """ 19 Abstract class to calculate Half-Sphere Exposure (HSE). 20 21 The HSE can be calculated based on the CA-CB vector, or the pseudo CB-CA 22 vector based on three consecutive CA atoms. This is done by two separate 23 subclasses. 24 """
25 - def __init__(self, model, radius, offset, hse_up_key, hse_down_key, 26 angle_key=None):
27 """ 28 @param model: model 29 @type model: L{Model} 30 31 @param radius: HSE radius 32 @type radius: float 33 34 @param offset: number of flanking residues that are ignored in the calculation 35 of the number of neighbors 36 @type offset: int 37 38 @param hse_up_key: key used to store HSEup in the entity.xtra attribute 39 @type hse_up_key: string 40 41 @param hse_down_key: key used to store HSEdown in the entity.xtra attribute 42 @type hse_down_key: string 43 44 @param angle_key: key used to store the angle between CA-CB and CA-pCB in 45 the entity.xtra attribute 46 @type angle_key: string 47 """ 48 assert(offset>=0) 49 # For PyMOL visualization 50 self.ca_cb_list=[] 51 ppb=CaPPBuilder() 52 ppl=ppb.build_peptides(model) 53 hse_map={} 54 hse_list=[] 55 hse_keys=[] 56 for pp1 in ppl: 57 for i in range(0, len(pp1)): 58 if i==0: 59 r1=None 60 else: 61 r1=pp1[i-1] 62 r2=pp1[i] 63 if i==len(pp1)-1: 64 r3=None 65 else: 66 r3=pp1[i+1] 67 # This method is provided by the subclasses to calculate HSE 68 result=self._get_cb(r1, r2, r3) 69 if result is None: 70 # Missing atoms, or i==0, or i==len(pp1)-1 71 continue 72 pcb, angle=result 73 hse_u=0 74 hse_d=0 75 ca2=r2['CA'].get_vector() 76 for pp2 in ppl: 77 for j in range(0, len(pp2)): 78 if pp1 is pp2 and abs(i-j)<=offset: 79 # neighboring residues in the chain are ignored 80 continue 81 ro=pp2[j] 82 if not is_aa(ro) or not ro.has_id('CA'): 83 continue 84 cao=ro['CA'].get_vector() 85 d=(cao-ca2) 86 if d.norm()<radius: 87 if d.angle(pcb)<(pi/2): 88 hse_u+=1 89 else: 90 hse_d+=1 91 res_id=r2.get_id() 92 chain_id=r2.get_parent().get_id() 93 # Fill the 3 data structures 94 hse_map[(chain_id, res_id)]=(hse_u, hse_d, angle) 95 hse_list.append((r2, (hse_u, hse_d, angle))) 96 hse_keys.append((chain_id, res_id)) 97 # Add to xtra 98 r2.xtra[hse_up_key]=hse_u 99 r2.xtra[hse_down_key]=hse_d 100 if angle_key: 101 r2.xtra[angle_key]=angle 102 AbstractPropertyMap.__init__(self, hse_map, hse_keys, hse_list)
103
104 - def _get_cb(self, r1, r2, r3):
105 """This method is provided by the subclasses to calculate HSE.""" 106 return NotImplemented
107
108 - def _get_gly_cb_vector(self, residue):
109 """ 110 Return a pseudo CB vector for a Gly residue. 111 The pseudoCB vector is centered at the origin. 112 113 CB coord=N coord rotated over -120 degrees 114 along the CA-C axis. 115 """ 116 try: 117 n_v=residue["N"].get_vector() 118 c_v=residue["C"].get_vector() 119 ca_v=residue["CA"].get_vector() 120 except: 121 return None 122 # center at origin 123 n_v=n_v-ca_v 124 c_v=c_v-ca_v 125 # rotation around c-ca over -120 deg 126 rot=rotaxis(-pi*120.0/180.0, c_v) 127 cb_at_origin_v=n_v.left_multiply(rot) 128 # move back to ca position 129 cb_v=cb_at_origin_v+ca_v 130 # This is for PyMol visualization 131 self.ca_cb_list.append((ca_v, cb_v)) 132 return cb_at_origin_v
133 134 135
136 -class HSExposureCA(_AbstractHSExposure):
137 """ 138 Class to calculate HSE based on the approximate CA-CB vectors, 139 using three consecutive CA positions. 140 """
141 - def __init__(self, model, radius=12, offset=0):
142 """ 143 @param model: the model that contains the residues 144 @type model: L{Model} 145 146 @param radius: radius of the sphere (centred at the CA atom) 147 @type radius: float 148 149 @param offset: number of flanking residues that are ignored in the calculation of the number of neighbors 150 @type offset: int 151 """ 152 _AbstractHSExposure.__init__(self, model, radius, offset, 153 'EXP_HSE_A_U', 'EXP_HSE_A_D', 'EXP_CB_PCB_ANGLE')
154
155 - def _get_cb(self, r1, r2, r3):
156 """ 157 Calculate the approximate CA-CB direction for a central 158 CA atom based on the two flanking CA positions, and the angle 159 with the real CA-CB vector. 160 161 The CA-CB vector is centered at the origin. 162 163 @param r1, r2, r3: three consecutive residues 164 @type r1, r2, r3: L{Residue} 165 """ 166 if r1 is None or r3 is None: 167 return None 168 try: 169 ca1=r1['CA'].get_vector() 170 ca2=r2['CA'].get_vector() 171 ca3=r3['CA'].get_vector() 172 except: 173 return None 174 # center 175 d1=ca2-ca1 176 d3=ca2-ca3 177 d1.normalize() 178 d3.normalize() 179 # bisection 180 b=(d1+d3) 181 b.normalize() 182 # Add to ca_cb_list for drawing 183 self.ca_cb_list.append((ca2, b+ca2)) 184 if r2.has_id('CB'): 185 cb=r2['CB'].get_vector() 186 cb_ca=cb-ca2 187 cb_ca.normalize() 188 angle=cb_ca.angle(b) 189 elif r2.get_resname()=='GLY': 190 cb_ca=self._get_gly_cb_vector(r2) 191 if cb_ca is None: 192 angle=None 193 else: 194 angle=cb_ca.angle(b) 195 else: 196 angle=None 197 # vector b is centered at the origin! 198 return b, angle
199
200 - def pcb_vectors_pymol(self, filename="hs_exp.py"):
201 """ 202 Write a PyMol script that visualizes the pseudo CB-CA directions 203 at the CA coordinates. 204 205 @param filename: the name of the pymol script file 206 @type filename: string 207 """ 208 if len(self.ca_cb_list)==0: 209 warnings.warn("Nothing to draw.", RuntimeWarning) 210 return 211 fp=open(filename, "w") 212 fp.write("from pymol.cgo import *\n") 213 fp.write("from pymol import cmd\n") 214 fp.write("obj=[\n") 215 fp.write("BEGIN, LINES,\n") 216 fp.write("COLOR, %.2f, %.2f, %.2f,\n" % (1.0, 1.0, 1.0)) 217 for (ca, cb) in self.ca_cb_list: 218 x,y,z=ca.get_array() 219 fp.write("VERTEX, %.2f, %.2f, %.2f,\n" % (x,y,z)) 220 x,y,z=cb.get_array() 221 fp.write("VERTEX, %.2f, %.2f, %.2f,\n" % (x,y,z)) 222 fp.write("END]\n") 223 fp.write("cmd.load_cgo(obj, 'HS')\n") 224 fp.close()
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