Biskit.Ramachandran

45

46 - def __init__( self, models, name=None, profileName='relAS', verbose=1 ):

47 """ 48 @param models: List of models display a Ramachandran plot for 49 @type models: [ PDBModel ] OR PDBModel 50 @param name: model name, will show up in plot 51 @type name: str 52 @param profileName: name of profile to use for coloring 53 (default: 'relAS') 54 @type profileName: str 55 @param verbose: verbosity level (default: 1) 56 @type verbose: 1|0 57 """ 58 if not biggles: 59 raise ImportError, 'biggles module could not be imported.' 60 61 if type(models) != type([]): 62 models = [ models ] 63 64 self.psi = [] 65 self.phi = [] 66 67 self.gly = [] 68 self.pro = [] 69 70 self.prof=[] 71 self.profileName = profileName 72 73 self.name=name 74 75 self.verbose = verbose 76 77 # calculate angles, profiles ... 78 self.calc( models ) 79 80 self.prof = N.ravel(self.prof) 81 self.gly = N.ravel(self.gly) 82 self.pro = N.ravel(self.pro)

83 84

85 - def calc( self, models ):

86 """ 87 Calculate angles, profiles and other things needed. 88 89 @param models: List of models 90 @type models: [ PDBModel ] 91 """ 92 res_count = 0 93 for m in models: 94 95 ## add profile if not there 96 if self.profileName: 97 self.prof += [ self.calcProfiles( m ) ] 98 99 ## calclate phi and psi angles for model 100 self.phi_and_psi( m ) 101 102 ## get list with GLY and PRO residue indices 103 gly_atomInd = m.indices(lambda a: a['residue_name']=='GLY') 104 gly_resInd = N.array( m.atom2resIndices( gly_atomInd ) ) 105 pro_atomInd = m.indices(lambda a: a['residue_name']=='PRO') 106 pro_resInd = N.array( m.atom2resIndices( pro_atomInd ) ) 107 self.gly.append( gly_resInd + res_count ) 108 self.pro.append( pro_resInd + res_count ) 109 res_count += m.lenResidues()

110 111

112 - def calcProfiles( self, m ):

113 """ 114 Calculate needed profiles. 115 116 @param m: PDBModel to calculate data for 117 @type m: PDBModel 118 """ 119 if self.verbose: print "Initiating PDBDope..." 120 d = PDBDope( m ) 121 122 if not self.profileName in m.aProfiles.keys(): 123 124 if self.profileName in ['MS', 'AS', 'curvature', 'relAS', 'relMS']: 125 if self.verbose: print "Adding SurfaceRacer profile...", 126 d.addSurfaceRacer() 127 128 if self.profileName in ['relASA']: 129 if self.verbose: print "Adding WhatIf ASA...", 130 d.addASA() 131 132 if self.profileName in ['density']: 133 if self.verbose: print "Adding surface density...", 134 d.addDensity() 135 136 if not self.profileName in m.rProfiles.keys(): 137 138 if self.profileName in ['cons_abs', 'cons_max', 'cons_ent']: 139 if self.verbose: print "Adding conservation data...", 140 d.addConservation() 141 142 if self.profileName in ['ASA_total', 'ASA_sc', 'ASA_bb']: 143 if self.verbose: print "Adding WhatIf ASA...", 144 d.addASA() 145 146 if self.verbose: print 'Done.' 147 148 ## convert atom profiles to average residue profile 149 if self.profileName in m.aProfiles.keys(): 150 prof = [] 151 aProfile = m.profile( self.profileName ) 152 resIdx = m.resIndex() 153 resIdx += [ m.lenAtoms()] 154 for i in range(len(resIdx)-1): 155 prof += [ N.average( N.take(aProfile, range(resIdx[i], 156 resIdx[i+1]) ) )] 157 else: 158 prof = m.profile( self.profileName ) 159 160 return prof

161 162

163 - def phi_and_psi( self, model ):

164 """ 165 Calculate phi and psi torsion angles for all 166 residues in model:: 167 168 phi - rotation about the N-CA bond 169 - last position in a chain = None 170 psi - rotation about CA-C 171 - first position in a chain = None 172 173 @param model: PDBModel 174 @type model: PDBModel 175 """ 176 for c in range( model.lenChains(breaks=1) ): 177 cModel = model.takeChains( [c], breaks=1 ) 178 179 xyz = cModel.xyz 180 181 xyz_CA = N.compress( cModel.maskCA(), xyz, 0 ) 182 xyz_N = N.compress( cModel.mask( ['N'] ), xyz, 0 ) 183 xyz_C = N.compress( cModel.mask( ['C'] ), xyz, 0 ) 184 185 ## phi: c1 - N 186 ## c2 - CA 187 ## c3 - C 188 ## c4 - N of next residue 189 for i in range( len(xyz_N)-1 ): 190 self.phi += [self.dihedral( xyz_N[i], xyz_CA[i], 191 xyz_C[i], xyz_N[i+1] )] 192 self.phi += [None] 193 194 ## psi: c1 - C of previous residue 195 ## c2 - N 196 ## c3 - CA 197 ## c4 - C 198 self.psi += [None] 199 for i in range( 1, len(xyz_N) ): 200 self.psi += [self.dihedral( xyz_C[i-1], xyz_N[i], 201 xyz_CA[i], xyz_C[i] )]

202 203

204 - def dihedral( self, coor1, coor2, coor3, coor4 ):

205 """ 206 Calculates the torsion angle of a set of four atom coordinates. 207 The dihedral angle returned is the angle between the projection 208 of i1-i2 and the projection of i4-i3 onto a plane normal to i2-i3. 209 210 @param coor1: coordinates 211 @type coor1: [float] 212 @param coor2: coordinates 213 @type coor2: [float] 214 @param coor3: coordinates 215 @type coor3: [float] 216 @param coor4: coordinates 217 @type coor4: [float] 218 """ 219 vec21 = coor2 - coor1 220 vec32 = coor3 - coor2 221 L = MU.cross( vec21, vec32 ) 222 L_norm = N.sqrt(sum(L**2)) 223 224 vec43 = coor4 - coor3 225 vec23 = coor2 - coor3 226 R = MU.cross( vec43, vec23 ) 227 R_norm = N.sqrt(sum(R**2)) 228 229 S = MU.cross( L, R ) 230 angle = sum( L*R ) / ( L_norm * R_norm ) 231 232 ## sometimes the value turns out to be ever so little greater than 233 ## one, to prevent N.arccos errors for this, set angle = 1.0 234 if angle > 1.0: angle = 1.0 235 236 if angle < -1.0: angle = -1.0 237 238 angle = N.arccos(angle) *180/N.pi 239 if sum(S*vec32) < 0.0: 240 angle = -angle 241 242 return angle

243 244

245 - def ramachandran( self ):

246 """ 247 Create all the ramachandran plot points. 248 249 @return: list of biggles.Point objects (all the points of the 250 plot)and a biggles.Inset object (property scale). 251 @rtype: [ biggles.Point ], biggles.Inset 252 """ 253 p = [] 254 255 ## calculate colors and create a legend if a property is given 256 if self.profileName: 257 palette = CS('plasma', 0, 100) 258 col = palette.color_array( self.prof ) 259 260 legend = Legend( palette.legend() ) 261 inset = biggles.Inset((1.1, 0.60), (1.2, .97), legend) 262 263 else: 264 col = ['black']*len(self.phi) 265 inset = None 266 267 ## add data points to plot 268 for i in range(len(self.phi)): 269 ## don't add termini - has missing angles 270 if self.phi[i] and self.psi[i]: 271 if i in self.gly: 272 p += [biggles.Point( self.psi[i], self.phi[i], 273 type="star", size=1, color=col[i] )] 274 elif i in self.pro: 275 p += [biggles.Point( self.psi[i], self.phi[i], 276 type="filled square", size=1, 277 color=col[i] )] 278 else: 279 p += [biggles.Point( self.psi[i], self.phi[i], 280 type="filled circle", size=1, 281 color=col[i] )] 282 return p, inset

283 284

285 - def ramachandran_background( self ):

286 """ 287 Creates a background (favoured regions) for a ramachandran plot. 288 289 @return: list of biggles.Point objects 290 @rtype: [ biggles.Point ] 291 """ 292 bg = [] 293 mat = biggles.read_matrix( T.projectRoot() + 294 '/external/biggles/ramachandran_bg.dat') 295 x, y = N.shape(mat) 296 for i in range(x): 297 for j in range(y): 298 if mat[i,j] < 200: 299 a = (360./y)*j - 180 300 b = (360./x)*(x-i)- 180 301 bg += [ biggles.Point( a, b, type="dot" )] 302 return bg

303 304

305 - def show( self, fileName=None ):

306 """ 307 Show ramachandran plot. 308 """ 309 plot = biggles.FramedPlot() 310 plot.xrange = (-180., 180.) 311 plot.yrange = (-180., 180.) 312 plot.xlabel = "$\Phi$" 313 plot.ylabel = "$\Psi$" 314 315 if self.name: 316 plot.title = self.name 317 318 ## add allowed regions 319 bg_plot = self.ramachandran_background( ) 320 for p in bg_plot: 321 plot.add( p ) 322 323 ## add ramachandran phi, psi valies 324 points, inset = self.ramachandran( ) 325 for p in points: 326 plot.add(p) 327 if inset: 328 plot.add( inset ) 329 330 plot.add( biggles.PlotLabel( 1.14, 0.55, self.profileName, size=2) ) 331 plot.add( biggles.PlotLabel( 1.1, 0.45, "GLY star", size=2) ) 332 plot.add( biggles.PlotLabel( 1.12, 0.40, "PRO square", size=2) ) 333 334 plot.show() 335 336 if fileName: 337 plot.write_eps( fileName )

Source Code for Module Biskit.Ramachandran