Source Code for Module Biskit.Dock.hexTools

  1  ## Automatically adapted for numpy.oldnumeric Mar 26, 2007 by alter_code1.py 
  2   
  3  ## 
  4  ## Biskit, a toolkit for the manipulation of macromolecular structures 
  5  ## Copyright (C) 2004-2006 Raik Gruenberg & Johan Leckner 
  6  ## 
  7  ## This program is free software; you can redistribute it and/or 
  8  ## modify it under the terms of the GNU General Public License as 
  9  ## published by the Free Software Foundation; either version 2 of the 
 10  ## License, or any later version. 
 11  ## 
 12  ## This program is distributed in the hope that it will be useful, 
 13  ## but WITHOUT ANY WARRANTY; without even the implied warranty of 
 14  ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
 15  ## General Public License for more details. 
 16  ## 
 17  ## You find a copy of the GNU General Public License in the file 
 18  ## license.txt along with this program; if not, write to the Free 
 19  ## Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 
 20  ## 
 21  ## 
 22  ## last $Author: graik $ 
 23  ## last $Date: 2007/03/26 18:40:42 $ 
 24  ## $Revision: 2.9 $ 
 25   
 26  """ 
 27  Various common functions used by the docking modeles. 
 28  """ 
 29   
 30  import numpy.oldnumeric as N 
 31  import os 
 32   
 33  import Biskit.tools as t 
 34  import tempfile 
 35  from Biskit import PDBDope, molUtils 
 36   
 37   
 38 -def createHexPdb_single( model, fout=None ): 
 39      """ 
 40      Write PDB of one structure for hex. 
 41       
 42      @param model: model 
 43      @type  model: PDBModel 
 44      @param fout: out file name (default: pdbCode + _hex.pdb) 
 45      @type  fout: str 
 46       
 47      @return: file name of result PDB 
 48      @rtype: str 
 49      """ 
 50      fout = fout or model.pdbCode + '_hex.pdb' 
 51      fout = t.absfile( fout ) 
 52      model.writePdb( fout ) 
 53      return fout 
 54   
 55  ##     dic = { 1:model } 
 56  ##     return createHexPdb( dic, fout ) 
 57   
 58   
 59 -def createHexPdb( modelDic, fout=None ): 
 60      """ 
 61      Write pdb for hex with models separated by MODEL%i/ENDMODEL 
 62       
 63      @param modelDic: dictionary mapping an integer to each model 
 64      @type  modelDic: dict {int:PCRModel} 
 65      @param fout: output name, default is pdbCode + _hex.pdb 
 66      @type  fout: str 
 67       
 68      @return: file name of result PDB 
 69      @rtype: str 
 70      """ 
 71      fout = fout or modelDic[1].pdbCode + '_hex.pdb' 
 72      fout = t.absfile( fout ) 
 73   
 74      ## open new file 
 75      out = open(fout, 'w') 
 76      ## fetch name for temporary file 
 77      tmpFile = tempfile.mktemp('_pcr2hex.pdb')   
 78   
 79      numbers = modelDic.keys() 
 80      numbers.sort() 
 81      for n in numbers: 
 82   
 83          ## write temporary pdb and open it 
 84          modelDic[ n ].writePdb( tmpFile ) 
 85   
 86          pdb_temp = open( tmpFile ) 
 87   
 88          out.write("MODEL%6i\n" % n) 
 89          lines = pdb_temp.readlines()    # get all lines 
 90          for line in lines: 
 91              if line[:3] <> 'END':       # filter out END 
 92                  out.write(line) 
 93          out.write("ENDMDL\n") 
 94   
 95          ## remove temporary file 
 96          pdb_temp.close() 
 97          os.remove( tmpFile ) 
 98   
 99      out.write("END") 
100      out.close() 
101      return fout 
102   
103   
104 -def centerSurfDist( model, surf_mask, mask=None ): 
105      """ 
106      Calculate the longest and shortest distance from 
107      the center of the molecule to the surface. 
108       
109      @param mask: atoms not to be considerd (default: None) 
110      @type  mask: [1|0] 
111      @param surf_mask: atom surface mask, needed for minimum surface distance 
112      @type  surf_mask: [1|0] 
113   
114      @return: max distance, min distance 
115      @rtype: float, float 
116      """ 
117      if mask is None: 
118              mask = model.maskHeavy() 
119   
120      ## calculate center of mass 
121      center = model.centerOfMass() 
122   
123      ## surface atom coordinates 
124      surf_xyz = N.compress( mask*surf_mask, model.getXyz(), 0 ) 
125   
126      ## find the atom closest and furthest away from center 
127      dist = N.sqrt( N.sum( (surf_xyz-center)**2 , 1 ) ) 
128      minDist = min(dist) 
129      maxDist = max(dist) 
130   
131      return maxDist, minDist 
132   
133   
134 -def createHexInp( recPdb, recModel, ligPdb, ligModel, comPdb=None, 
135                    outFile=None, macDock=None, silent=0, sol=512 ): 
136      """ 
137      Prepare a Hex macro file for the docking of the receptor(s) 
138      against ligand(s). 
139   
140      @param recPdb: hex-formatted PDB 
141      @type  recPdb: str 
142      @param recModel: hex-formatted PDB 
143      @type  recModel: str 
144      @param ligPdb: PDBModel, get distances from this one 
145      @type  ligPdb: PDBModel 
146      @param ligModel: PDBModel, getdistances from this one 
147      @type  ligModel: PDBModel 
148      @param comPdb: reference PDB 
149      @type  comPdb: str 
150      @param outFile: base of file name for mac and out 
151      @type  outFile: str 
152       
153      @param macDock: None -> hex decides (from the size of the molecule), 
154                      1 -> force macroDock, 0-> force off (default: None) 
155      @type  macDock: None|1|0 
156      @param silent: don't print distances and macro warnings (default: 0) 
157      @type  silent: 0|1 
158      @param sol: number of solutions that HEx should save (default: 512) 
159      @type  sol: int 
160       
161      @return: HEX macro file name, HEX out generated bu the macro, 
162               macro docking status 
163      @rtype: str, str, boolean 
164      """ 
165      ## files and names 
166      recCode = t.stripFilename( recPdb )[0:4]           
167      ligCode = t.stripFilename( ligPdb )[0:4] 
168   
169      outFile = outFile or recCode + '-' + ligCode 
170   
171      ## hex macro name 
172      macName = t.absfile( outFile + '_hex.mac' ) 
173   
174      ## hex rotation matrix output name 
175      outName_all = t.absfile( outFile + '_hex.out'  ) 
176      outName_clust = t.absfile( outFile + '_hex_cluster.out') 
177   
178      ## add surface profiles if not there 
179      if not recModel.aProfiles.has_key('relAS'): 
180          #t.flushPrint('\nCalculating receptor surface profile') 
181          rec_asa = PDBDope( recModel ) 
182          rec_asa.addSurfaceRacer() 
183      if not ligModel.aProfiles.has_key('relAS'): 
184          #t.flushPrint('\nCalculating ligand surface profile') 
185          lig_asa = PDBDope( ligModel ) 
186          lig_asa.addSurfaceRacer() 
187   
188      ## surface masks, > 95% exposed 
189      rec_surf_mask = N.greater( recModel.profile('relAS'), 95 ) 
190      lig_surf_mask = N.greater( ligModel.profile('relAS'), 95 ) 
191   
192      ## maximun and medisn distance from centre of mass to any surface atom 
193      recMax, recMin = centerSurfDist( recModel, rec_surf_mask ) 
194      ligMax, ligMin = centerSurfDist( ligModel, lig_surf_mask ) 
195   
196      ## approxinate max and min center to centre distance 
197      maxDist = recMax + ligMax  
198      minDist = recMin + ligMin 
199   
200      ## molecular separation and search range to be used in the docking 
201      molSep = ( maxDist + minDist ) / 2 
202      molRange = 2 * ( maxDist - molSep ) 
203   
204      if not silent: 
205          print 'Docking setup: %s\nRecMax: %.1f RecMin: %.1f\nLigMax: %.1f LigMin: %.1f\nMaxDist: %.1f MinDist: %.1f\nmolecular_separation: %.1f r12_range: %.1f\n'%(outFile, recMax, recMin, ligMax, ligMin, maxDist, minDist, molSep, molRange) 
206   
207      if recMax > 30 and ligMax > 30 and not silent: 
208          print '\nWARNING! Both the receptor and ligand radius is ', 
209          print 'greater than 30A.\n'      
210   
211      ## determine docking mode to use 
212      macroDocking = 0 
213   
214      if macDock==None: 
215          if recMax > 35 and not silent: 
216              print '\nReceptor has a radius that exceeds 35A ', 
217              print '-> Macro docking will be used' 
218              macroDocking = 1 
219      else: 
220          macroDocking = macDock 
221   
222      ##################### 
223      ## write macro file 
224   
225      macOpen= open( macName, 'w') 
226   
227      macOpen.write('# -- ' + macName + ' --\n') 
228      macOpen.write(' \n') 
229      macOpen.write('open_receptor '+ t.absfile(recPdb) +'\n') 
230      macOpen.write('open_ligand '+ t.absfile(ligPdb) +'\n') 
231   
232      if comPdb and comPdb[-4:] == '.pdb': 
233          macOpen.write('open_complex '+comPdb+'\n') 
234   
235      macOpen.write('\n') 
236   
237      head = """ 
238  # -------------- general settings ---------------- 
239  disc_cache 1                   # disc cache on (0 off) 
240  docking_sort_mode 1            # Sort solutions by cluster (0 by energy) 
241  docking_cluster_mode 1         # Display all clusters (0 display best) 
242  docking_cluster_threshold 2.00 
243  # docking_cluster_bumps  number 
244   
245  # ------------ molecule orientation -------------- 
246  molecule_separation %(separation)i 
247  commit_view """%({'separation': round(molSep)} ) 
248   
249   
250      macro =""" 
251  # -------------- macro docking ------------------- 
252  macro_min_coverage 25 
253  macro_sphere_radius 15 
254  macro_docking_separation 25 
255  activate_macro_model""" 
256   
257   
258      tail = """ 
259  # -------------- docking setup ------------------- 
260  docking_search_mode 0          # full rotational search 
261   
262  receptor_range_angle  180      # 0, 15, 30, 45, 60, 75, 90, 180 
263  docking_receptor_samples 720   # 362, 492, 642, 720, 980, 1280 
264   
265  ligand_range_angle  180 
266  docking_ligand_samples 720 
267   
268  twist_range_angle 360          # 0, 15, 30, 60, 90, 180, 360 
269  docking_alpha_samples 128      # 64, 128, 256 
270   
271  r12_step 0.500000              # 0.1, 0.2, 0.25, 0.5, 0.75, 1, 1.5, 2 
272  r12_range %(range)i 
273   
274  docking_radial_filter 0        # Radial Envelope Filter - None 
275   
276  grid_size 0.600                # 0.4, 0.5, 0.6, 0.75, 1.0 
277  # docking_electrostatics 0       # use only surface complimentarity 
278  docking_electrostatics 1      # use electrostatic term for scoring clusters 
279   
280  docking_main_scan 16     #  
281  docking_main_search 26 
282   
283  max_docking_solutions %(nr_sol)i # number of solutions to save 
284   
285  # -------------- post-processing ---------------- 
286  docking_refine 0    # None 
287  #  docking_refine 1    # Backbone Bumps 
288  #  docking_refine 2    # MM energies 
289  #  docking_refine 3    # MM minimization 
290   
291  # ---------------- run docking ------------------ 
292  activate_docking 
293  #  save_docking %(output_clust)s 
294  #  save_range 1 512 ./ dock .pdb 
295   
296  # ------------ also save all solutions ---------- 
297  docking_sort_mode 0            # Sort solutions by energy (1 by cluster) 
298  save_docking %(output_all)s""" \ 
299      %({'range':round(molRange), 'output_all':outName_all, 
300         'nr_sol':int(sol), 'output_clust':outName_clust} ) 
301   
302      macOpen.writelines( head ) 
303   
304      ## macro docking will not work with multiple models, if both are added to 
305      ## the hex macro file - macrodocking will be skipped during the docking run 
306      if macroDocking: 
307          macOpen.writelines( macro ) 
308   
309      macOpen.writelines( tail ) 
310   
311      macOpen.close() 
312   
313      return macName, outName_all, macroDocking 
314   
315   
316   
317  ############# 
318  ##  TESTING         
319  ############# 
320  import Biskit.test as BT 
321           
322 -class Test(BT.BiskitTest): 
323      """Test case""" 
324   
325 -    def test_hexTools(self): 
326          """Dock.hexTools test""" 
327          from Biskit import PDBModel 
328          self.m = PDBModel( t.testRoot() + '/com/1BGS.pdb' ) 
329          dist = centerSurfDist( self.m , self.m.maskCA() ) 
330   
331          self.assertAlmostEqual( dist[0], 26.880979538, 7) 
332   
333   
334  if __name__ == '__main__': 
335   
336      BT.localTest() 
337