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.Polypeptide import CaPPBuilder, is_aa 
 15  from Bio.PDB.Vector import rotaxis 
 16   
 17   
18 -class _AbstractHSExposure(AbstractPropertyMap):
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
26 - def __init__(self, model, radius, offset, hse_up_key, hse_down_key, 27 angle_key=None):
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 """Method is provided by the subclasses to calculate HSE (PRIVATE).""" 106 return NotImplemented
107
108 - def _get_gly_cb_vector(self, residue):
109 """Return a pseudo CB vector for a Gly residue (PRIVATE). 110 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 Exception: 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 -class HSExposureCA(_AbstractHSExposure):
136 """Class to calculate HSE based on the approximate CA-CB vectors. 137 138 Uses three consecutive CA positions. 139 """ 140
141 - def __init__(self, model, radius=12, offset=0):
142 """Initialse class. 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 """Calculate approx CA-CB direction (PRIVATE). 158 159 Calculate the approximate CA-CB direction for a central 160 CA atom based on the two flanking CA positions, and the angle 161 with the real CA-CB vector. 162 163 The CA-CB vector is centered at the origin. 164 165 @param r1, r2, r3: three consecutive residues 166 @type r1, r2, r3: L{Residue} 167 """ 168 if r1 is None or r3 is None: 169 return None 170 try: 171 ca1 = r1['CA'].get_vector() 172 ca2 = r2['CA'].get_vector() 173 ca3 = r3['CA'].get_vector() 174 except Exception: 175 return None 176 # center 177 d1 = ca2 - ca1 178 d3 = ca2 - ca3 179 d1.normalize() 180 d3.normalize() 181 # bisection 182 b = (d1 + d3) 183 b.normalize() 184 # Add to ca_cb_list for drawing 185 self.ca_cb_list.append((ca2, b + ca2)) 186 if r2.has_id('CB'): 187 cb = r2['CB'].get_vector() 188 cb_ca = cb - ca2 189 cb_ca.normalize() 190 angle = cb_ca.angle(b) 191 elif r2.get_resname() == 'GLY': 192 cb_ca = self._get_gly_cb_vector(r2) 193 if cb_ca is None: 194 angle = None 195 else: 196 angle = cb_ca.angle(b) 197 else: 198 angle = None 199 # vector b is centered at the origin! 200 return b, angle
201
202 - def pcb_vectors_pymol(self, filename="hs_exp.py"):
203 """Write PyMol script for visualization. 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 """Class to calculate HSE based on the real CA-CB vectors.""" 231
232 - def __init__(self, model, radius=12, offset=0):
233 """Initialize class. 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 """Method to calculate CB-CA vector. 249 250 @param r1, r2, r3: three consecutive residues (only r2 is used) 251 @type r1, r2, r3: L{Residue} 252 """ 253 if r2.get_resname() == 'GLY': 254 return self._get_gly_cb_vector(r2), 0.0 255 else: 256 if r2.has_id('CB') and r2.has_id('CA'): 257 vcb = r2['CB'].get_vector() 258 vca = r2['CA'].get_vector() 259 return (vcb - vca), 0.0 260 return None
261 262
263 -class ExposureCN(AbstractPropertyMap):
264 """Residue exposure as number of CA atoms around its CA atom.""" 265
266 - def __init__(self, model, radius=12.0, offset=0):
267 """Initialize. 268 269 A residue's exposure is defined as the number of CA atoms around 270 that residues CA atom. A dictionary is returned that uses a L{Residue} 271 object as key, and the residue exposure as corresponding value. 272 273 @param model: the model that contains the residues 274 @type model: L{Model} 275 276 @param radius: radius of the sphere (centred at the CA atom) 277 @type radius: float 278 279 @param offset: number of flanking residues that are ignored in the calculation of the number of neighbors 280 @type offset: int 281 282 """ 283 assert(offset >= 0) 284 ppb = CaPPBuilder() 285 ppl = ppb.build_peptides(model) 286 fs_map = {} 287 fs_list = [] 288 fs_keys = [] 289 for pp1 in ppl: 290 for i in range(0, len(pp1)): 291 fs = 0 292 r1 = pp1[i] 293 if not is_aa(r1) or not r1.has_id('CA'): 294 continue 295 ca1 = r1['CA'] 296 for pp2 in ppl: 297 for j in range(0, len(pp2)): 298 if pp1 is pp2 and abs(i - j) <= offset: 299 continue 300 r2 = pp2[j] 301 if not is_aa(r2) or not r2.has_id('CA'): 302 continue 303 ca2 = r2['CA'] 304 d = (ca2 - ca1) 305 if d < radius: 306 fs += 1 307 res_id = r1.get_id() 308 chain_id = r1.get_parent().get_id() 309 # Fill the 3 data structures 310 fs_map[(chain_id, res_id)] = fs 311 fs_list.append((r1, fs)) 312 fs_keys.append((chain_id, res_id)) 313 # Add to xtra 314 r1.xtra['EXP_CN'] = fs 315 AbstractPropertyMap.__init__(self, fs_map, fs_keys, fs_list)
316