Source Code for Module Biskit.WhatIf

  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.6 $ 
 21  ## last $Date: 2006/09/13 10:28:14 $ 
 22  ## last $Author: leckner $ 
 23  """ 
 24  Run WhatIf, calculate and collect accessibility data 
 25  """ 
 26   
 27  import Numeric as N 
 28  import string 
 29  import tempfile 
 30   
 31  from Biskit import Executor, TemplateError 
 32  ## import Biskit.settings as S 
 33  import Biskit.tools as T 
 34   
 35   
 36 -class WhatIf_Error( Exception ): 
 37      pass 
 38   
 39 -class WhatIf( Executor ): 
 40      """ 
 41      WhatIf 
 42      ====== 
 43      Run WhatIf accessibility calculation 
 44       
 45      Result 
 46      ------ 
 47      A WhatIf log file containing the accessibilities (in square Angstrom) 
 48      per residue. The ASA is specified for the entire residue as well for 
 49      the backbone and side chain. 
 50   
 51      Reference 
 52      --------- 
 53      U{ http://swift.cmbi.kun.nl/whatif/}     
 54   
 55      @note: WhatIf is not free software, therefore the calculations performed 
 56             by this class can now also be performed by the freeware program 
 57             SurfaceRacer, see L{ Biskit.SurfaceRacer }.  
 58      """ 
 59   
 60      whatif_script ="""SOUP  \nINISOU 
 61   
 62  GETMOL  \n%(f_pdb)s  \nTEST 
 63   
 64  ACCESS  \nINIACC 
 65  SETACC  \nALL  \n0   \n      
 66   
 67  DOLOG  \n%(f_relativeASA)s  \n0 
 68  VACACC  \nN  \nALL 
 69  NOLOG 
 70   
 71  DOLOG  \n%(f_residueASA)s  \n0 
 72  SHOACC  \nALL  \n0 
 73  NOLOG"""  
 74   
 75   
 76 -    def __init__( self, model, **kw ): 
 77          """ 
 78          @param model: PDBModel 
 79          @type  model:  
 80   
 81          @param kw: additional key=value parameters for Executor: 
 82          @type  kw: key=value pairs 
 83          :: 
 84            debug    - 0|1, keep all temporary files (default: 0) 
 85            verbose  - 0|1, print progress messages to log (log != STDOUT) 
 86            node     - str, host for calculation (None->local) NOT TESTED 
 87                            (default: None) 
 88            nice     - int, nice level (default: 0) 
 89            log      - Biskit.LogFile, program log (None->STOUT) (default: None) 
 90          """ 
 91          Executor.__init__( self, 'whatif', template=self.whatif_script, 
 92                             f_out='/dev/null', **kw ) 
 93   
 94          self.f_pdb = tempfile.mktemp('_whatif.pdb') 
 95          self.f_relativeASA = tempfile.mktemp('_whatif_relative.log') 
 96          self.f_residueASA = tempfile.mktemp('_whatif_residue.log') 
 97   
 98          self.model = model.clone( deepcopy=1 ) 
 99   
100   
101 -    def prepare( self ): 
102          """ 
103          Overrides Executor method. 
104          """ 
105          self.__prepareModel( self.model, self.f_pdb )  
106   
107   
108 -    def __prepareModel( self, m, f_pdb_out ): 
109          """ 
110          Prepare a model that WhatIf likes. 
111           - consecutive numbering of residues 
112           - remove hydrogens (will be deleted anyway) 
113           
114          @param m: model  
115          @type  m: PDBModel 
116          @param f_pdb_out: name of pdb file to write 
117          @type  f_pdb_out: str 
118          """ 
119          ## Whatif deletes hydrogens anyway 
120          m = m.compress( m.maskHeavy() ) 
121          m.renumberResidues() 
122          m.writePdb( f_pdb_out, ter=0 ) 
123   
124   
125 -    def cleanup( self ): 
126          """ 
127          Tidy up the mess you created. 
128          """ 
129          Executor.cleanup( self ) 
130   
131          if not self.debug: 
132              T.tryRemove( self.f_pdb ) 
133              T.tryRemove( self.f_relativeASA ) 
134              T.tryRemove( self.f_residueASA ) 
135   
136              T.tryRemove('FOR???.DAT', wildcard=1) 
137              T.tryRemove('pdbout.tex') 
138              T.tryRemove('pdbout.txt') 
139              T.tryRemove('TEXSTORE.DAT') 
140              T.tryRemove('TEXTABLE.DAT') 
141   
142   
143 -    def isFailed( self ): 
144          """ 
145          Overrides Executor method 
146          """ 
147          return self.error is None  
148   
149   
150 -    def finish( self ): 
151          """ 
152          Overrides Executor method 
153          """ 
154          Executor.finish( self ) 
155          self.result = self.parse_whatif( ) 
156   
157   
158 -    def __read_residueASA( self ): 
159          """ 
160          Read solvent accessibility calculated with WHATIF and return array 
161          of ASA-values. First column is the total ASA, second column the ASA 
162          of the backbone, the third column is the ASA of the side-chain. 
163   
164          @return: array (3 x len(residues)) with ASA values 
165          @rtype: array 
166          """ 
167          ## [14:-27] skip header and tail lines line 
168   
169          lines = open(self.f_residueASA).readlines()[14:-27] 
170          ASA_values = map(lambda line: string.split(line)[-3:], lines) 
171          ASA_values = N.array(map(lambda row: map(string.atof, row), 
172                                   ASA_values)) 
173   
174          return ASA_values 
175   
176   
177 -    def __read_relativeASA( self ): 
178          """ 
179          Read the relative ASA data. Return dic indexed by resNumber 
180          with atom and relative accessibility (%) 
181   
182          @return: relative atom ASA dictionary. e.g:: 
183                   { 1: { 'N':34.6, 'CA':81.5 .. }, 2 : {.. }.... } 
184          @rtype: dict of dict 
185          """ 
186          AtomASA = {}  
187   
188          logFile = open( self.f_relativeASA  ) 
189   
190          while 1: 
191              line = logFile.readline() 
192   
193              # first line of new residue 
194              if line[:8] == 'Residue:': 
195                  # get residue number 
196                  residue = int( string.split(line)[1] ) 
197                  AtomASA[residue] = {} 
198   
199                  # skip two lines 
200                  line = logFile.readline()  
201                  line = logFile.readline()  
202                  # read first line of atom list 
203                  line = logFile.readline() 
204                  while line[:20] != '                    ': 
205   
206                      atom = string.split(line)[0] 
207                      value = float( string.split(line)[-1] ) 
208                      AtomASA[ residue ][ atom ] = value 
209   
210                      line = logFile.readline() 
211   
212              if line == '': #eof 
213                  break 
214   
215          return AtomASA 
216   
217   
218 -    def __exposedResidues( self, ASA_values, sidechainCut=0.0, 
219                           backboneCut=0.0, totalCut=0.0  ): 
220          """ 
221          Decide what is a surface exposed residue and what is not. 
222          sidechainCut, backboneCut, totalCut - float, cutoff value 
223          for what will be considered as a exposed residue. All 
224          three values have to pass the test. 
225   
226          @param ASA_values: array with ASA values for side chains, backbone 
227                             and total calculated in L{__read_residueASA}. 
228          @type  ASA_values: array 
229          @param sidechainCut: cutoff ASA value for considering the side chain 
230                               to consider thew residue being exposed 
231                               (default: 0.0)  
232          @type  sidechainCut: float 
233          @param backboneCut: cutoffvalue for back bone ASA 
234          @type  backboneCut: float  
235          @param totalCut: cutoff for total ASA 
236          @type  totalCut: float    
237   
238          @return: residue mask, where 0 = burried 
239          @rtype: [1|0] 
240          """ 
241          col_0 = N.greater( N.transpose(ASA_values)[0], totalCut ) 
242          col_1 = N.greater( N.transpose(ASA_values)[1], backboneCut ) 
243          col_2 = N.greater( N.transpose(ASA_values)[2], sidechainCut ) 
244   
245          col_012 = N.concatenate( ([col_0],[col_1],[col_2]) )  
246   
247          exposedList = N.greater(N.sum(col_012), 0) 
248   
249          return exposedList 
250   
251   
252 -    def parse_whatif( self ): 
253          """ 
254          Takes a PDBModel and returns three lists:: 
255   
256            relativeASA - a list of atom relative accessibilities 
257                          (in %  - compared to the resideu in a 
258                          Gly-XXX-Gly tripeptide, an approximation 
259                          of the unfolded state) 
260   
261                        -> N.array( 1 x N_atoms ) 
262   
263            residueASA  - array 3 x N_residues of float 
264                          acc.surf.area per res: sidechain, backbone, total 
265   
266   
267            burriedResidues - a burried atom mask 
268                              note: cutoff values set above 
269   
270                            ->  [0, 1, 0, 0, 1, 1, 1, 0 ...] 
271   
272          @return: relative accessibility, residue ASA and burried atom mask 
273          @rtype: array, array, [1|0] 
274          """ 
275          # read ASA per residue log 
276          resASA = self.__read_residueASA( ) 
277   
278          # and decide what is a burried residue 
279          burriedResidues = self.__exposedResidues( resASA ) 
280   
281          # read the relative assesibility per atom log 
282          atomRelASA = self.__read_relativeASA( ) 
283   
284          return self.__asaAtomList( atomRelASA ), resASA, burriedResidues 
285   
286   
287 -    def __asaAtomList( self, relativeASA ): 
288          """ 
289          convert ASA dict to list of values with len == atom number 
290           
291          @param relativeASA: see __read_relativeASA() 
292          @type  relativeASA: dict 
293   
294          @return: relative accessibilities e.g. [ 98.0, 0.0, 12.0, .. ], 
295                   Numpy array 1 x N(atoms) 
296          @rtype:  [float] 
297          """ 
298          result = [] 
299   
300          resIndex = 1 
301          for res in self.model.resList(): 
302   
303              for atom in res: 
304   
305                  if relativeASA.has_key( resIndex ) \ 
306                         and relativeASA[ resIndex ].has_key( atom['name'] ): 
307   
308                      result.append( relativeASA[ resIndex ][ atom['name' ]] ) 
309                  else: 
310                      result.append( 0.0 ) 
311   
312              resIndex += 1 
313   
314          return result 
315   
316   
317  ############# 
318  ##  TESTING         
319  ############# 
320           
321 -class Test: 
322      """ 
323      Test class 
324      """ 
325       
326 -    def run( self, local=0 ): 
327          """ 
328          run function test 
329           
330          @param local: transfer local variables to global and perform 
331                        other tasks only when run locally 
332          @type  local: 1|0 
333           
334          @return: sum of the total ASA for each residue 
335          @rtype:  float 
336          """ 
337          from Biskit import PDBModel 
338   
339          ## Loading PDB... 
340   
341          f = T.testRoot()+"/com/1BGS.pdb" 
342          m = PDBModel(f) 
343   
344          m = m.compress( m.maskProtein() ) 
345          m = m.compress( m.maskHeavy() ) 
346   
347          ## Starting WhatIf 
348          x = WhatIf( m, debug=0, verbose=0 ) 
349   
350          ## Running 
351          atomAcc, resAcc, resMask = x.run() 
352   
353          if local: 
354              ## check that model hasn't changed 
355              m_ref = PDBModel(f) 
356              m_ref = m.compress( m.maskProtein() ) 
357              for k in m_ref.atoms[0].keys(): 
358                  ref = [ m_ref.atoms[0][k] for i in range( m_ref.lenAtoms() ) ] 
359                  mod = [ m.atoms[0][k] for i in range( m.lenAtoms() ) ] 
360                  if not ref == mod: 
361                      print 'Not equal ', k 
362                  else: 
363                      print 'Equal ', k 
364   
365              ## display exposed residues in PyMol 
366              from Pymoler import Pymoler 
367              pm = Pymoler() 
368              model = pm.addPdb( m, '1' ) 
369              pm.colorRes( '1', resAcc[:,0] ) 
370              pm.show() 
371   
372              print "\nResult for first 10 atoms/residues: " 
373              print '\nAccessability (A^2):\n', atomAcc[:10] 
374              print '\nResidue accessability (A^2)' 
375              print '[total, backbone, sidechain]:\n', resAcc[:10] 
376              print '\nExposed residue mask:\n',resMask[:10] 
377              print '\nTotal atom    accessability (A^2): %.2f'%sum(atomAcc)  
378              print '      residue accessability (A^2): %.2f'%sum(resAcc)[0] 
379   
380              globals().update( locals() ) 
381                                 
382          return N.sum(resAcc[:,0]) 
383   
384   
385 -    def expected_result( self ): 
386          """ 
387          Precalculated result to check for consistent performance. 
388   
389          @return: sum of the total ASA for each residue 
390          @rtype:  float 
391          """ 
392          return 2814.6903 
393       
394   
395  if __name__ == '__main__': 
396   
397      test = Test() 
398   
399      assert abs( test.run( local=1 ) - test.expected_result() ) < 1e-8 
400