Biskit.ReduceCoordinates

1 ## 2 ## Biskit, a toolkit for the manipulation of macromolecular structures 3 ## Copyright (C) 2004-2005 Raik Gruenberg & Johan Leckner 4 ## 5 ## This program is free software; you can redistribute it and/or 6 ## modify it under the terms of the GNU General Public License as 7 ## published by the Free Software Foundation; either version 2 of the 8 ## License, or any later version. 9 ## 10 ## This program is distributed in the hope that it will be useful, 11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of 12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13 ## General Public License for more details. 14 ## 15 ## You find a copy of the GNU General Public License in the file 16 ## license.txt along with this program; if not, write to the Free 17 ## Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18 ## 19 ## 20 ## $Revision: 2.4 $ 21 ## last $Date: 2006/09/12 06:34:36 $ 22 ## last $Author: leckner $ 23 24 """ 25 Create structures with reduced number of atoms. 26 """ 27 28 from PDBModel import PDBModel 29 import Numeric as N 30 import tools as T 31 import molUtils as MU 32 33

34 -class ReduceCoordinates:

35 """ 36 ReduceCoordinates 37 ================= 38 39 Translate a PDBModel or frames from a trajectory to structure(s) with 40 only one backbone and up to 2 side chain atoms per residue. 41 The new atoms are the centers of mass of several atoms and carry the 42 weight of the pooled atoms in an atom profile called 'mass'. 43 44 Examples 45 -------- 46 >>> ## create with reference PDBModel 47 >>> reducer = ReduceCoordinates( m_ref ) 48 >>> ## creates reduced PDBModel from m_ref 49 >>> m_red = reducer.reduceToModel() 50 51 OR: 52 >>> m_red_1 = reducer.reduceToModel( m1.getXyz() ) ## reduce many models 53 >>> m_red_2 = reducer.reduceToModel( m2.getXyz() ) ## with identical atoms 54 55 OR: 56 >>> ## reduce a complete Trajectory 57 >>> reducer = ReduceCoordinates( traj.ref ) 58 >>> red_ref= reducer.reduceToModel() 59 >>> frames = reducer.reduceXyz( traj.frames ) 60 >>> traj_red = Trajectory( ref=red_ref ) 61 >>> traj_red.frames = frames 62 """ 63 ## modify order of TYR/PHE ring atoms to move centers away from ring axis 64 aaAtoms = MU.aaAtoms 65 aaAtoms['TYR'] = ['N','CA','C','O','CB','CG','CD1','CE1','CD2', 66 'CE2','CZ','OH', 'OXT'] 67 aaAtoms['PHE'] = ['N','CA','C','O','CB','CG','CD1','CE1','CD2', 68 'CE2','CZ', 'OXT'] 69

70 - def __init__( self, model, maxPerCenter=4 ):

71 """ 72 Prepare reduction of coordinates from a given model. 73 74 @param model: reference model defining atom content and order 75 @type model: PDBModel 76 @param maxPerCenter: max number of atoms per side chain center atom 77 (default: 4) 78 @type maxPerCenter: int 79 """ 80 self.m = model 81 self.__addMassProfile( self.m ) 82 83 ## sort atoms within residues into standard order 84 def cmpAtoms( a1, a2 ): 85 """ 86 Comparison function for bringing atoms into standard order 87 within residues as defined by L{ aaAtoms }. 88 89 @param a1: model 90 @type a1: PDBModel 91 @param a2: model 92 @type a2: PDBModel 93 94 @return: int or list of matching positions 95 @rtype: [-1|0|1] 96 """ 97 res = a1['residue_name'] 98 target = self.aaAtoms[ res ] 99 try: 100 return cmp(target.index( a1['name'] ), 101 target.index( a2['name'] )) 102 except ValueError, why: 103 return cmp( a1['name'], a2['name'] )

104 ## s = "Unknown atom for %s %i: %s or %s" % \ 105 ## (res, a1['residue_number'], a1['name'], a2['name'] ) 106 ## raise PDBError( s ) 107 108 self.a_indices = self.m.argsort( cmpAtoms ) 109 self.m_sorted = self.m.sort( self.a_indices ) 110 111 ## remove H from internal model and from list of atom positions 112 maskH = self.m_sorted.remove( self.m_sorted.maskH() ) 113 self.a_indices = N.compress( maskH, self.a_indices ) 114 115 self.makeMap( maxPerCenter )

116 117

118 - def __addMassProfile( self, model ):

119 """ 120 Add a mass profile to the model. 121 122 @param model: model 123 @type model: PDBModel 124 """ 125 h = MU.atomMasses['H'] 126 masses = [] 127 for a in model.atoms: 128 if a['element'] == 'H': 129 masses += [ 0 ] 130 else: 131 masses += [ MU.atomMasses[ a['element'] ] + 132 h * MU.aaAtomsH[a['residue_name']][a['name']] ] 133 134 model.setAtomProfile( 'mass', masses )

135 136

137 - def group( self, a_indices, maxPerCenter ):

138 """ 139 Group a bunch of integers (atom indices in PDBModel) so that each 140 group has at most maxPerCenter items. 141 142 @param a_indices: atom indices 143 @type a_indices: [int] 144 @param maxPerCenter: max entries per group 145 @type maxPerCenter: int 146 147 @return: list of lists of int 148 @rtype: [[int],[int]..] 149 """ 150 ## how many groups are necessary? 151 n_centers = len( a_indices ) / maxPerCenter 152 if len( a_indices ) % maxPerCenter: 153 n_centers += 1 154 155 ## how many items/atoms go into each group? 156 nAtoms = N.ones(n_centers, 'i') * int(len( a_indices ) / n_centers) 157 i=0 158 while N.sum(nAtoms) != len( a_indices ): 159 nAtoms[i] += 1 160 i += 1 161 162 ## distribute atom indices into groups 163 result = [] 164 pos = 0 165 for n in nAtoms: 166 result += [ N.take( a_indices, N.arange(n) + pos) ] 167 pos += n 168 169 return result

170 171

172 - def nextAtom( self, resName, resNumber, name, chainId, segid):

173 """ 174 Create an atom dictionary. 175 176 @param resName: residue name 177 @type resName: str 178 @param resNumber: residue number 179 @type resNumber: int 180 @param name: atom name 181 @type name: str 182 @param chainId: chain identifier 183 @type chainId: str 184 @param segid: segnemt identifier 185 @type segid: str 186 187 @return: atom dictionary 188 @rtype: dict 189 """ 190 self.currentAtom += 1 191 192 return {'residue_name':resName, 'name':name, 'type':'ATOM', 193 'residue_number':resNumber, 194 'serial_number': self.currentAtom, 195 'segment_id': segid, 'chain_id':chainId, 196 'name_original':name, 'element':'X'}

197 198

199 - def makeMap( self, maxPerCenter=4 ):

200 """ 201 Calculate mapping between complete and reduced atom list. 202 Creates a (list of lists of int, list of atom dictionaries) 203 containing groups of atom indices into original model, new center atoms 204 205 @param maxPerCenter: max number of atoms per side chain center atom 206 (default: 4) 207 @type maxPerCenter: int 208 """ 209 210 resIndex = self.m_sorted.resIndex() 211 resModels= self.m_sorted.resModels() 212 m = self.m_sorted 213 214 self.currentAtom = 0 215 216 groups = [] 217 atoms = [] 218 219 for i in range( len( resIndex ) ): 220 221 first_atom = resIndex[ i ] 222 223 if i < len( resIndex )-1: 224 last_atom = resIndex[ i+1 ] - 1 225 else: 226 last_atom = len( self.a_indices ) - 1 227 228 res_name = m.atoms[ first_atom ]['residue_name'] 229 segid = m.atoms[ first_atom ]['segment_id'] 230 chainId = m.atoms[ first_atom ]['chain_id'] 231 res_number = m.atoms[first_atom]['serial_number'] 232 233 ## position of this residue's atoms in original PDBModel (unsorted) 234 a_indices = self.a_indices[ first_atom : last_atom+1 ] 235 236 ## for each center create list of atom indices and a center atom 237 if res_name != 'GLY' and res_name != 'ALA': 238 239 bb_a_indices = N.compress( resModels[i].maskBB(), a_indices) 240 sc_a_indices = N.compress( 241 N.logical_not( resModels[i].maskBB()), a_indices ) 242 243 sc_groups = self.group( sc_a_indices, maxPerCenter ) 244 245 else: 246 bb_a_indices = a_indices 247 sc_groups = [] 248 249 groups += [ bb_a_indices ] 250 atoms += [ self.nextAtom( res_name, res_number, 'BB', 251 chainId, segid ) ] 252 253 i = 0 254 for g in sc_groups: 255 groups += [ g ] 256 atoms += [ self.nextAtom( res_name, res_number, 'SC%i'%i, 257 chainId, segid) ] 258 i += 1 259 260 self.groups = groups 261 self.atoms = atoms

262 263

264 - def reduceXyz( self, xyz, axis=0 ):

265 """ 266 Reduce the number of atoms in the given coordinate set. The set must 267 have the same length and order as the reference model. It may have 268 an additional (time) dimension as first axis. 269 270 @param xyz: coordinates (N_atoms x 3) or (N_frames x N_atoms x 3) 271 @type xyz: array 272 @param axis: axis with atoms (default: 0) 273 @type axis: int 274 275 @return: coordinate array (N_less_atoms x 3) or 276 (N_frames x N_less_atoms x 3) 277 @rtype: array 278 """ 279 masses = self.m.atomProfile('mass') 280 r_xyz = None 281 282 for atom_indices in self.groups: 283 284 x = N.take( xyz, atom_indices, axis ) 285 m = N.take( masses, atom_indices ) 286 287 center = N.sum( x * N.transpose([m,]), axis=axis) / N.sum( m ) 288 289 if axis == 0: 290 center = center[N.NewAxis, :] 291 292 if axis == 1: 293 center = center[:, N.NewAxis, :] 294 295 if r_xyz == None: 296 r_xyz = center 297 298 else: 299 r_xyz = N.concatenate( (r_xyz, center), axis ) 300 301 return r_xyz

302 303

304 - def reduceToModel( self, xyz=None, reduce_profiles=1 ):

305 """ 306 Create a reduced PDBModel from coordinates. Atom profiles the source 307 PDBModel are reduced by averaging over the grouped atoms. 308 309 @param xyz: coordinte array (N_atoms x 3) or 310 None (->use reference coordinates) 311 @type xyz: array OR None 312 313 @return: PDBModel with reduced atom set and profile 'mass' 314 @rtype: PDBModel 315 """ 316 317 mass = self.m.atomProfile('mass') 318 xyz = xyz or self.m.getXyz() 319 320 mProf = [ N.sum( N.take( mass, group ) ) for group in self.groups ] 321 xyz = self.reduceXyz( xyz ) 322 323 result = PDBModel() 324 result.setAtoms( self.atoms ) 325 result.setXyz( xyz ) 326 result.setAtomProfile( 'mass', mProf ) 327 328 if reduce_profiles: 329 self.reduceAtomProfiles( self.m, result ) 330 331 result.rProfiles = self.m.rProfiles 332 333 return result

334 335

336 - def reduceAtomProfiles( self, from_model, to_model ):

337 """ 338 reduce all atom profiles according to the calculated map by calculating 339 the average over the grouped atoms. 340 341 @param from_model: model 342 @type from_model: PDBModel 343 @param to_model: model 344 @type to_model: PDBModel 345 """ 346 for profname in from_model.aProfiles: 347 348 p0 = from_model.atomProfile(profname) 349 info = from_model.profileInfo( profname ) 350 351 pr = [ N.average( N.take( p0, group ) ) for group in self.groups ] 352 353 to_model.setAtomProfile( profname, pr ) 354 to_model.setProfileInfo( profname, **info )

355 356 357 358 ############# 359 ## TESTING 360 ############# 361

362 -class Test:

363 """ 364 Test class 365 """ 366

367 - def run( self, local=0 ):

368 """ 369 run function test 370 371 @param local: transfer local variables to global and perform 372 other tasks only when run locally 373 @type local: 1|0 374 375 @return: atoms after reduction 376 @rtype: int 377 """ 378 m = PDBModel( T.testRoot()+'/com/1BGS.pdb' ) 379 m = m.compress( N.logical_not( m.maskH2O() ) ) 380 381 m.setAtomProfile('test', range(len(m))) 382 383 red = ReduceCoordinates( m, 4 ) 384 385 mred = red.reduceToModel() 386 387 if local: 388 print 'Atoms before reduction %i'%m.lenAtoms() 389 print 'Atoms After reduction %i'%mred.lenAtoms() 390 globals().update( locals() ) 391 392 return mred.lenAtoms()

393 394

395 - def expected_result( self ):

396 """ 397 Precalculated result to check for consistent performance. 398 399 @return: atoms after reduction 400 @rtype: int 401 """ 402 return 445

403 404 405 406 407 if __name__ == '__main__': 408 409 test = Test() 410 411 assert test.run( local=1 ) == test.expected_result() 412

Source Code for Module Biskit.ReduceCoordinates