Source Code for Module Biskit.Dock.Analyzer

  1  ## Automatically adapted for numpy.oldnumeric Mar 26, 2007 by alter_code1.py 
  2   
  3  #!/usr/bin/env python 
  4  ## 
  5  ## Biskit, a toolkit for the manipulation of macromolecular structures 
  6  ## Copyright (C) 2004-2006 Raik Gruenberg & Johan Leckner 
  7  ## 
  8  ## This program is free software; you can redistribute it and/or 
  9  ## modify it under the terms of the GNU General Public License as 
 10  ## published by the Free Software Foundation; either version 2 of the 
 11  ## License, or any later version. 
 12  ## 
 13  ## This program is distributed in the hope that it will be useful, 
 14  ## but WITHOUT ANY WARRANTY; without even the implied warranty of 
 15  ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
 16  ## General Public License for more details. 
 17  ## 
 18  ## You find a copy of the GNU General Public License in the file 
 19  ## license.txt along with this program; if not, write to the Free 
 20  ## Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 
 21  ## 
 22  ## 
 23  ## $Revision: 2.13 $ 
 24  ## last $Date: 2007/04/06 10:16:08 $ 
 25  ## last $Author: graik $ 
 26   
 27  """ 
 28  Analyze HEX docking result. 
 29  """ 
 30   
 31  import Biskit.tools as t 
 32  from Biskit import Trajectory, mathUtils,  molUtils 
 33   
 34   
 35  import numpy.oldnumeric as N 
 36  import numpy.oldnumeric.random_array as RandomArray 
 37  import copy 
 38   
 39  import biggles 
 40   
 41 -class AnalyzeError( Exception ): 
 42      pass 
 43   
 44 -class Analyzer: 
 45   
 46 -    def __init__( self, verbose=1, **options ): 
 47          """ 
 48          @param verbose: verbosity level (default: 1) 
 49          @type  verbose: 1|0         
 50          @param options: needs:: 
 51                           rec,lig - file name, receptor, ligand trajectories 
 52                           ref     - file name, pickled reference complex 
 53          @type  options: any 
 54   
 55          @raise AnalyzeError: if atoms are not aligned 
 56          """ 
 57          self.options = options 
 58   
 59          if verbose: t.flushPrint("\nLoading...") 
 60          self.t_lig = t.Load( options['lig'] ) 
 61          self.t_rec = t.Load( options['rec'] ) 
 62          self.com= t.Load( options['ref'] ) 
 63   
 64          ## delete H from all players 
 65          self.t_lig.removeAtoms( self.t_lig.getRef().maskH() ) 
 66          self.t_rec.removeAtoms( self.t_rec.getRef().maskH() ) 
 67   
 68          self.com.rec_model.remove( self.com.rec_model.maskH() ) 
 69          self.com.lig_model.remove( self.com.lig_model.maskH() ) 
 70          self.com.lig_transformed = None 
 71   
 72          ## equalize atom content of free (trajectory) and bound models 
 73          if verbose: t.flushPrint('\nCasting...') 
 74   
 75          bnd_rec = self.com.rec() 
 76          bnd_lig = self.com.lig_model 
 77   
 78          self.t_rec.sortAtoms() 
 79          self.t_lig.sortAtoms() 
 80   
 81          bnd_rec = bnd_rec.sort() 
 82          bnd_lig = bnd_lig.sort() 
 83   
 84          #m_bnd_rec, m_t_rec = bnd_rec.equalAtoms( self.t_rec.getRef() ) 
 85          #m_bnd_lig, m_t_lig = bnd_lig.equalAtoms( self.t_lig.getRef() ) 
 86   
 87          i_bnd_rec, i_t_rec = bnd_rec.compareAtoms( self.t_rec.getRef() ) 
 88          i_bnd_lig, i_t_lig = bnd_lig.compareAtoms( self.t_lig.getRef() ) 
 89   
 90          #self.t_rec.removeAtoms( N.logical_not( m_t_rec ) ) 
 91          #self.t_lig.removeAtoms( N.logical_not( m_t_lig ) ) 
 92   
 93          self.t_rec = self.t_rec.takeAtoms( i_t_rec ) 
 94          self.t_lig = self.t_lig.takeAtoms( i_t_lig ) 
 95   
 96          #self.mask_free_lig = m_t_lig 
 97          #self.mask_free_rec = m_t_rec 
 98          self.i_free_lig = i_t_lig 
 99          self.i_free_rec = i_t_rec 
100   
101          #bnd_rec.remove( N.logical_not( m_bnd_rec ) ) 
102          #bnd_lig.remove( N.logical_not( m_bnd_lig ) ) 
103          bnd_rec = bnd_rec.take( i_bnd_rec ) 
104          bnd_lig = bnd_lig.take( i_bnd_lig ) 
105   
106          #self.mask_bnd_rec = m_bnd_rec 
107          #self.mask_bnd_lig = m_bnd_lig 
108          self.i_bnd_rec = i_bnd_rec 
109          self.i_bnd_lig = i_bnd_lig 
110   
111          ## put 'equalized' models back into ref complex 
112          self.com.rec_model = bnd_rec 
113          self.com.lig_model = bnd_lig 
114          self.com.lig_transformed = None 
115   
116          ## check 
117          if not self.t_rec.getRef().equals( self.com.rec() )[1] or \ 
118             not self.t_lig.getRef().equals( self.com.lig() )[1]: 
119   
120              raise AnalyzeError('Atoms are not aligned.') 
121   
122          ## native contact matrix 
123          self.contacts = self.com.atomContacts() 
124   
125          self.hexContacts = None 
126   
127   
128 -    def __categorizeHexSurf(self, cutoff=0.1): 
129          """ 
130          Compare complexes of list to native complex to see if 
131          their contact surfaces overlapp with the native complex. 
132           
133          @param cutoff: fraction cutoff for defining a overlap (default: 0.1) 
134          @type  cutoff: float 
135           
136          @return: list of len(self.hexContacts) overlapping with 
137                   native contact surface of lig and rec (0 - no overlap, 
138                   1 - rec OR lig overlapps, 2- rec AND lig overlapps) 
139          @rtype: [0|1|2] 
140          """ 
141          result = [ self.com.fractionNativeSurface( c, self.contacts ) 
142                     for c in self.hexContacts ] 
143   
144          result = [ N.sum( N.greater( o, cutoff ) ) for o in result ] 
145          return result 
146   
147   
148 -    def setHexComplexes(self, com_lst, n=20 ): 
149          """ 
150          add contact matrices of hex-generated (wrong) complexes for comparison 
151           
152          @param com_lst: ComplexList with contacts calculated 
153          @type  com_lst: ComplexList 
154          """ 
155          t.flushPrint('adding hex-generated complexes') 
156   
157          self.hexContacts = [] 
158   
159  #        f = lambda a: molUtils.elementType(a['element']) == 'p' 
160   
161          i = 0 
162          while i < n: 
163   
164              com = com_lst[i] 
165              t.flushPrint('#') 
166   
167              try: 
168                  if com['fractNatCont'] == 0.0:#com['fnac_10'] == 0.0: 
169                      com.rec_model.remove( com.rec().maskH() ) 
170                      com.lig_model.remove( com.lig_model.maskH() ) 
171   
172                      com.rec_model = com.rec_model.sort() 
173                      com.lig_model = com.lig_model.sort() 
174   
175                      com.rec_model.keep( self.i_free_rec ) 
176                      com.lig_model.keep( self.i_free_lig ) 
177                      com.lig_transformed = None 
178   
179                      self.hexContacts += [ com.atomContacts() ] 
180   
181                  else: 
182                      n += 1 
183                  i+= 1 
184              except: 
185                  print t.lastError() 
186   
187          self.hexSurfaces = self.__categorizeHexSurf( 0.2 ) 
188   
189   
190 -    def random_contacts( self, contMat, n, maskRec=None, maskLig=None ): 
191          """ 
192          Create randomized surface contact matrix with same number of 
193          contacts and same shape as given contact matrix. 
194           
195          @param contMat: template contact matrix 
196          @type  contMat: matrix 
197          @param n: number of matrices to generate 
198          @type  n: int 
199          @param maskRec: surface masks (or something similar) 
200          @type  maskRec: [1|0] 
201          @param maskLig: surface masks (or something similar) 
202          @type  maskLig: [1|0] 
203           
204          @return: list of [n] random contact matricies 
205          @rtype: [matrix] 
206          """ 
207          a,b = N.shape( contMat ) 
208          nContacts = N.sum( N.sum( contMat )) 
209   
210          if not maskLig: 
211              r_size, l_size = N.shape( contMat ) 
212              maskLig = N.ones( l_size ) 
213              maskRec = N.ones( r_size ) 
214   
215          c_mask = N.ravel( N.outerproduct( maskRec, maskLig ) ) 
216          c_pos = N.nonzero( c_mask ) 
217   
218          # get array with surface positions from complex 
219          cont = N.take( N.ravel(contMat), c_pos ) 
220          length = len( cont ) 
221   
222          result = [] 
223   
224          for i in range( n ): 
225              # create random array 
226              ranCont = mathUtils.randomMask( nContacts,length ) 
227   
228              # blow up to size of original matrix 
229              r = N.zeros(a*b) 
230              N.put( r, c_pos, ranCont) 
231   
232              result += [ N.reshape( r, (a,b) ) ] 
233   
234          return result 
235   
236   
237 -    def __shuffleList(self, lst ): 
238          """ 
239          shuffle order of lst 
240   
241          @param lst: list to shuffle 
242          @type  lst: [any] 
243           
244          @return: shuffeled list 
245          @rtype: [any] 
246          """ 
247          pos = RandomArray.permutation( len( lst )) 
248          return N.take( lst, pos ) 
249   
250   
251 -    def shuffledLists( self, n, lst, mask=None ): 
252          """ 
253          shuffle order of a list n times, leaving masked(0) elements untouched 
254   
255          @param n: number of times to shuffle the list 
256          @type  n: int 
257          @param lst: list to shuffle 
258          @type  lst: [any] 
259          @param mask: mask to be applied to lst 
260          @type  mask: [1|0] 
261   
262          @return: list of shuffeled lists 
263          @rtype: [[any]]         
264          """ 
265          if not mask: 
266              mask = N.ones( len(lst)  ) 
267   
268          if type( lst ) == list: 
269              lst = N.array( lst ) 
270           
271          pos = N.nonzero( mask ) 
272   
273          rand_pos = N.array( [ self.__shuffleList( pos ) for i in range(n) ] ) 
274   
275          result = [] 
276          for p in rand_pos: 
277   
278              r = copy.copy( lst ) 
279              N.put( r, p, N.take( lst, pos ) ) 
280              result += [r] 
281   
282          return result 
283   
284   
285 -    def __plotSequence(self, list, **arg): 
286          """ 
287          add curve to current plot, with x1..xn = 0..n and y1..yn = |list| 
288   
289          @param list: list to plot 
290          @type  list: [any]         
291          """ 
292   
293          self.plot.add( biggles.Curve( range( len(list) ), list, **arg ) ) 
294   
295          if arg.has_key('label'): 
296              self.plot.lastY -= .05 
297   
298              self.plot.add( biggles.PlotLabel( self.plot.lastX, self.plot.lastY, 
299                                                arg['label'], **arg )) 
300   
301 -    def report(self): 
302          """ 
303          override for actual plotting 
304          """ 
305          pass 
306   
307   
308 -    def initPlot(self): 
309          """ 
310          override for plot creation 
311          """ 
312          self.page = biggles.FramedPlot() 
313          self.plot = self.page 
314   
315   
316 -    def plotAll( self ): 
317          """ 
318          Show plot 
319          """ 
320          self.initPlot() 
321          self.report() 
322          self.page.show() 
323   
324   
325   
326  ############# 
327  ##  TESTING         
328  ############# 
329  import Biskit.test as BT 
330           
331 -class Test(BT.BiskitTest): 
332      """Hmmer test""" 
333   
334       
335 -    def prepare(self): 
336          import tempfile 
337          self.f_out = tempfile.mktemp( '_test_rec.traj' ) 
338   
339 -    def cleanUp(self): 
340          t.tryRemove( self.f_out ) 
341   
342 -    def test_Analyzer( self): 
343          """Dock.Analyzer test """ 
344          from Biskit import Trajectory 
345          from Biskit.Dock import ComplexList 
346   
347          ## create a minimal 1-frame receptor trajectory from a pdb file 
348          self.t_rec = Trajectory( [t.testRoot()+'/rec/1A2P.pdb'], 
349                                   verbose=self.local) 
350          t.Dump( self.t_rec, self.f_out ) 
351   
352          ## load a complex list 
353          cl = t.Load( t.testRoot() + '/dock/hex/complexes.cl') 
354   
355          self.a= Analyzer( rec = self.f_out, 
356                            lig = t.testRoot()+'/lig_pcr_00/traj.dat', 
357                            ref = t.testRoot()+'/com/ref.complex', 
358                            verbose = self.local) 
359   
360          ## shuffle this list five times 
361          shuff_lst = self.a.shuffledLists( 5, range(8) ) 
362   
363          ## create two random contact matrices 
364          rand_mat = self.a.random_contacts( cl[0].atomContacts(), 2 ) 
365           
366          self.assertEqual( N.shape(rand_mat[1]), (1075, 876) )  
367       
368           
369  if __name__ == '__main__': 
370      
371      BT.localTest() 
372