Back to index

python-biopython  1.60
HSExposure.py
Go to the documentation of this file.
00001 # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
00002 # This code is part of the Biopython distribution and governed by its
00003 # license.  Please see the LICENSE file that should have been included
00004 # as part of this package.
00005 
00006 """Half-sphere exposure and coordination number calculation."""
00007 
00008 import warnings
00009 from math import pi
00010 
00011 from Bio.PDB.AbstractPropertyMap import AbstractPropertyMap
00012 from Bio.PDB.PDBParser import PDBParser
00013 from Bio.PDB.Polypeptide import CaPPBuilder, is_aa
00014 from Bio.PDB.Vector import rotaxis
00015 
00016 
00017 class _AbstractHSExposure(AbstractPropertyMap):
00018     """
00019     Abstract class to calculate Half-Sphere Exposure (HSE).
00020 
00021     The HSE can be calculated based on the CA-CB vector, or the pseudo CB-CA
00022     vector based on three consecutive CA atoms. This is done by two separate 
00023     subclasses. 
00024     """
00025     def __init__(self, model, radius, offset, hse_up_key, hse_down_key, 
00026             angle_key=None):
00027         """
00028         @param model: model
00029         @type model: L{Model}
00030 
00031         @param radius: HSE radius
00032         @type radius: float
00033 
00034         @param offset: number of flanking residues that are ignored in the calculation
00035         of the number of neighbors
00036         @type offset: int
00037 
00038         @param hse_up_key: key used to store HSEup in the entity.xtra attribute
00039         @type hse_up_key: string
00040 
00041         @param hse_down_key: key used to store HSEdown in the entity.xtra attribute
00042         @type hse_down_key: string
00043 
00044         @param angle_key: key used to store the angle between CA-CB and CA-pCB in 
00045             the entity.xtra attribute
00046         @type angle_key: string
00047         """
00048         assert(offset>=0)
00049         # For PyMOL visualization
00050         self.ca_cb_list=[]
00051         ppb=CaPPBuilder()
00052         ppl=ppb.build_peptides(model)
00053         hse_map={}
00054         hse_list=[]
00055         hse_keys=[]
00056         for pp1 in ppl:
00057             for i in range(0, len(pp1)):
00058                 if i==0:
00059                     r1=None
00060                 else:
00061                     r1=pp1[i-1]
00062                 r2=pp1[i]
00063                 if i==len(pp1)-1:
00064                     r3=None
00065                 else:
00066                     r3=pp1[i+1]
00067                 # This method is provided by the subclasses to calculate HSE
00068                 result=self._get_cb(r1, r2, r3)
00069                 if result is None:
00070                     # Missing atoms, or i==0, or i==len(pp1)-1
00071                     continue
00072                 pcb, angle=result
00073                 hse_u=0
00074                 hse_d=0
00075                 ca2=r2['CA'].get_vector()
00076                 for pp2 in ppl:
00077                     for j in range(0, len(pp2)):
00078                         if pp1 is pp2 and abs(i-j)<=offset:
00079                             # neighboring residues in the chain are ignored 
00080                             continue
00081                         ro=pp2[j]
00082                         if not is_aa(ro) or not ro.has_id('CA'):
00083                             continue
00084                         cao=ro['CA'].get_vector()
00085                         d=(cao-ca2)
00086                         if d.norm()<radius:
00087                             if d.angle(pcb)<(pi/2):
00088                                 hse_u+=1
00089                             else:
00090                                 hse_d+=1
00091                 res_id=r2.get_id()
00092                 chain_id=r2.get_parent().get_id()
00093                 # Fill the 3 data structures
00094                 hse_map[(chain_id, res_id)]=(hse_u, hse_d, angle)
00095                 hse_list.append((r2, (hse_u, hse_d, angle)))
00096                 hse_keys.append((chain_id, res_id))
00097                 # Add to xtra
00098                 r2.xtra[hse_up_key]=hse_u
00099                 r2.xtra[hse_down_key]=hse_d
00100                 if angle_key:
00101                     r2.xtra[angle_key]=angle
00102         AbstractPropertyMap.__init__(self, hse_map, hse_keys, hse_list)
00103 
00104     def _get_cb(self, r1, r2, r3):
00105         """This method is provided by the subclasses to calculate HSE."""
00106         return NotImplemented
00107 
00108     def _get_gly_cb_vector(self, residue):
00109         """
00110         Return a pseudo CB vector for a Gly residue.
00111         The pseudoCB vector is centered at the origin.
00112 
00113         CB coord=N coord rotated over -120 degrees 
00114         along the CA-C axis.
00115         """
00116         try:
00117             n_v=residue["N"].get_vector()
00118             c_v=residue["C"].get_vector()
00119             ca_v=residue["CA"].get_vector()
00120         except:
00121             return None
00122         # center at origin
00123         n_v=n_v-ca_v
00124         c_v=c_v-ca_v
00125         # rotation around c-ca over -120 deg
00126         rot=rotaxis(-pi*120.0/180.0, c_v)
00127         cb_at_origin_v=n_v.left_multiply(rot)
00128         # move back to ca position
00129         cb_v=cb_at_origin_v+ca_v
00130         # This is for PyMol visualization
00131         self.ca_cb_list.append((ca_v, cb_v))
00132         return cb_at_origin_v
00133 
00134 
00135 
00136 class HSExposureCA(_AbstractHSExposure):
00137     """
00138     Class to calculate HSE based on the approximate CA-CB vectors,
00139     using three consecutive CA positions.
00140     """
00141     def __init__(self, model, radius=12, offset=0):
00142         """
00143         @param model: the model that contains the residues
00144         @type model: L{Model}
00145 
00146         @param radius: radius of the sphere (centred at the CA atom)
00147         @type radius: float
00148 
00149         @param offset: number of flanking residues that are ignored in the calculation            of the number of neighbors
00150         @type offset: int
00151         """
00152         _AbstractHSExposure.__init__(self, model, radius, offset, 
00153                 'EXP_HSE_A_U', 'EXP_HSE_A_D', 'EXP_CB_PCB_ANGLE')
00154 
00155     def _get_cb(self, r1, r2, r3):
00156         """
00157         Calculate the approximate CA-CB direction for a central
00158         CA atom based on the two flanking CA positions, and the angle
00159         with the real CA-CB vector. 
00160         
00161         The CA-CB vector is centered at the origin.
00162 
00163         @param r1, r2, r3: three consecutive residues
00164         @type r1, r2, r3: L{Residue}
00165         """
00166         if r1 is None or r3 is None:
00167             return None
00168         try:
00169             ca1=r1['CA'].get_vector()
00170             ca2=r2['CA'].get_vector()
00171             ca3=r3['CA'].get_vector()
00172         except:
00173             return None
00174         # center
00175         d1=ca2-ca1
00176         d3=ca2-ca3
00177         d1.normalize()
00178         d3.normalize()
00179         # bisection
00180         b=(d1+d3)
00181         b.normalize()
00182         # Add to ca_cb_list for drawing
00183         self.ca_cb_list.append((ca2, b+ca2))
00184         if r2.has_id('CB'):
00185             cb=r2['CB'].get_vector()
00186             cb_ca=cb-ca2
00187             cb_ca.normalize()
00188             angle=cb_ca.angle(b)
00189         elif r2.get_resname()=='GLY':
00190             cb_ca=self._get_gly_cb_vector(r2)
00191             if cb_ca is None:
00192                 angle=None
00193             else:
00194                 angle=cb_ca.angle(b)
00195         else:
00196             angle=None
00197         # vector b is centered at the origin!
00198         return b, angle
00199 
00200     def pcb_vectors_pymol(self, filename="hs_exp.py"):
00201         """
00202         Write a PyMol script that visualizes the pseudo CB-CA directions 
00203         at the CA coordinates.
00204 
00205         @param filename: the name of the pymol script file
00206         @type filename: string
00207         """
00208         if len(self.ca_cb_list)==0:
00209             warnings.warn("Nothing to draw.", RuntimeWarning)
00210             return
00211         fp=open(filename, "w")
00212         fp.write("from pymol.cgo import *\n")
00213         fp.write("from pymol import cmd\n")
00214         fp.write("obj=[\n")
00215         fp.write("BEGIN, LINES,\n")
00216         fp.write("COLOR, %.2f, %.2f, %.2f,\n" % (1.0, 1.0, 1.0))
00217         for (ca, cb) in self.ca_cb_list:
00218             x,y,z=ca.get_array()
00219             fp.write("VERTEX, %.2f, %.2f, %.2f,\n" % (x,y,z))
00220             x,y,z=cb.get_array()
00221             fp.write("VERTEX, %.2f, %.2f, %.2f,\n" % (x,y,z))
00222         fp.write("END]\n")
00223         fp.write("cmd.load_cgo(obj, 'HS')\n")
00224         fp.close()
00225 
00226 
00227 class HSExposureCB(_AbstractHSExposure):
00228     """
00229     Class to calculate HSE based on the real CA-CB vectors.
00230     """
00231     def __init__(self, model, radius=12, offset=0):
00232         """
00233         @param model: the model that contains the residues
00234         @type model: L{Model}
00235 
00236         @param radius: radius of the sphere (centred at the CA atom)
00237         @type radius: float
00238 
00239         @param offset: number of flanking residues that are ignored in the calculation            of the number of neighbors
00240         @type offset: int
00241         """
00242         _AbstractHSExposure.__init__(self, model, radius, offset,
00243                 'EXP_HSE_B_U', 'EXP_HSE_B_D')
00244 
00245     def _get_cb(self, r1, r2, r3):
00246         """
00247         Method to calculate CB-CA vector.
00248 
00249         @param r1, r2, r3: three consecutive residues (only r2 is used)
00250         @type r1, r2, r3: L{Residue}
00251         """
00252         if r2.get_resname()=='GLY':
00253             return self._get_gly_cb_vector(r2), 0.0
00254         else:
00255             if r2.has_id('CB') and r2.has_id('CA'):
00256                 vcb=r2['CB'].get_vector()
00257                 vca=r2['CA'].get_vector()
00258                 return (vcb-vca), 0.0
00259         return None
00260 
00261 
00262 class ExposureCN(AbstractPropertyMap):
00263     def __init__(self, model, radius=12.0, offset=0):
00264         """
00265         A residue's exposure is defined as the number of CA atoms around 
00266         that residues CA atom. A dictionary is returned that uses a L{Residue}
00267         object as key, and the residue exposure as corresponding value.
00268 
00269         @param model: the model that contains the residues
00270         @type model: L{Model}
00271 
00272         @param radius: radius of the sphere (centred at the CA atom)
00273         @type radius: float
00274 
00275         @param offset: number of flanking residues that are ignored in the calculation            of the number of neighbors
00276         @type offset: int
00277 
00278         """
00279         assert(offset>=0)
00280         ppb=CaPPBuilder()
00281         ppl=ppb.build_peptides(model)
00282         fs_map={}
00283         fs_list=[]
00284         fs_keys=[]
00285         for pp1 in ppl:
00286             for i in range(0, len(pp1)):
00287                 fs=0
00288                 r1=pp1[i]
00289                 if not is_aa(r1) or not r1.has_id('CA'):
00290                     continue
00291                 ca1=r1['CA']
00292                 for pp2 in ppl:
00293                     for j in range(0, len(pp2)):
00294                         if pp1 is pp2 and abs(i-j)<=offset:
00295                             continue
00296                         r2=pp2[j]
00297                         if not is_aa(r2) or not r2.has_id('CA'):
00298                             continue
00299                         ca2=r2['CA']
00300                         d=(ca2-ca1)
00301                         if d<radius:
00302                             fs+=1
00303                 res_id=r1.get_id()
00304                 chain_id=r1.get_parent().get_id()
00305                 # Fill the 3 data structures
00306                 fs_map[(chain_id, res_id)]=fs
00307                 fs_list.append((r1, fs))
00308                 fs_keys.append((chain_id, res_id))
00309                 # Add to xtra
00310                 r1.xtra['EXP_CN']=fs
00311         AbstractPropertyMap.__init__(self, fs_map, fs_keys, fs_list)
00312 
00313 
00314 if __name__=="__main__":
00315 
00316     import sys
00317 
00318     p=PDBParser()
00319     s=p.get_structure('X', sys.argv[1])
00320     model=s[0]
00321 
00322     # Neighbor sphere radius
00323     RADIUS=13.0
00324     OFFSET=0
00325 
00326     hse=HSExposureCA(model, radius=RADIUS, offset=OFFSET)
00327     for l in hse:
00328         print l
00329     print
00330 
00331     hse=HSExposureCB(model, radius=RADIUS, offset=OFFSET)
00332     for l in hse:
00333         print l
00334     print
00335 
00336     hse=ExposureCN(model, radius=RADIUS, offset=OFFSET)
00337     for l in hse:
00338         print l
00339     print
00340 
00341     for c in model:
00342         for r in c:
00343             try:
00344                 print r.xtra['PCB_CB_ANGLE']
00345             except:
00346                 pass
00347 
00348