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