Source Code for Module Biskit.AmberCrdParser

  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  ## $Revision: 2.9 $ 
 23  ## last $Date: 2007/04/06 10:16:06 $ 
 24  ## last $Author: graik $ 
 25  """ 
 26  Convert single amber crd into Trajectory object 
 27  """ 
 28   
 29  import re 
 30  import numpy.oldnumeric as N 
 31  import sys 
 32   
 33  import tools as T 
 34  from Trajectory import Trajectory 
 35  from PDBModel import PDBModel 
 36  from LogFile import StdLog 
 37   
 38 -def _use(): 
 39   
 40      print """ 
 41  Convert single amber crd into Trajectory object 
 42   
 43  amber2traj.py -i sim.crd -o traj_0.dat -r ref.pdb [-b -wat -hyd -rnres 
 44                -code PDBC ] 
 45   
 46      -b     traj has box info (3 additional coordinates per frame) 
 47      -wat   delete WAT, Cl-, Na+ residues (after parsing) 
 48      -hyd   delete all hydrogens (after parsing) 
 49      -rnres rename amber residues HIE/HID/HIP, CYX to HIS and CYS 
 50      -code  PDB code of molecule (otherwise first 4 letters of ref file name) 
 51   
 52      ref.pdb must have identical atom content as sim.crd 
 53      """ 
 54      sys.exit( 0 ) 
 55   
 56 -class ParseError( Exception ): 
 57      pass 
 58   
 59   
 60 -class AmberCrdParser: 
 61      """ 
 62      Convert an Amber-generated crd file into a Trajectory object. 
 63      """ 
 64   
 65 -    def __init__( self, fcrd, fref, box=0, rnAmber=0, pdbCode=None, 
 66                    log=StdLog(), verbose=0 ): 
 67          """ 
 68          @param fcrd: path to input coordinate file 
 69          @type  fcrd: str 
 70          @param fref: PDB or pickled PDBModel with same atom content and order 
 71          @type  fref: str 
 72          @param box: expect line with box info at the end of each frame 
 73                      (default: 0) 
 74          @type  box: 1|0 
 75          @param rnAmber: rename amber style residues into standard (default: 0) 
 76          @type  rnAmber: 1|0 
 77          @param pdbCode: pdb code to be put into the model (default: None) 
 78          @type  pdbCode: str 
 79          @param log: LogFile instance [Biskit.StdLog] 
 80          @type  log: Biskit.LogFile 
 81          @param verbose: print progress to log [0] 
 82          @type  verbose: int 
 83          """ 
 84          self.fcrd = T.absfile( fcrd ) 
 85          self.crd  = T.gzopen( self.fcrd ) 
 86   
 87          self.ref  = PDBModel( T.absfile(fref), pdbCode=pdbCode ) 
 88          self.box  = box 
 89   
 90          self.n = self.ref.lenAtoms() 
 91   
 92          self.log = log 
 93          self.verbose = verbose 
 94   
 95          if rnAmber: 
 96              self.ref.renameAmberRes() 
 97   
 98          ## pre-compile pattern for line2numbers 
 99          xnumber = "-*\d+\.\d+"              # optionally negtive number 
100          xspace  = ' *'                      # one or more space char 
101          self.xnumbers = re.compile('('+xspace+xnumber+')') 
102   
103          ## pre-compute lines expected per frame 
104          self.lines_per_frame = self.n * 3 / 10 
105   
106          if self.n % 10 != 0:  self.lines_per_frame += 1 
107          if self.box:          self.lines_per_frame += 1 
108   
109          ## mark chains (the TER position might get lost when deleting atoms) 
110          ## should not be necessary any longer 
111  ##        if not self.ref.getAtoms()[0].get('chain_id',''): 
112  ##            self.ref.addChainId() 
113   
114 -    def line2numbers( self, l ): 
115          """ 
116          convert a line from crd/vel file to list of float numbers 
117           
118          @param l: line 
119          @type  l: str 
120           
121          @return: list of floats 
122          @rtype: [float] 
123          """ 
124          match = self.xnumbers.findall( l ) 
125   
126          return [ round( float(strCrd),3) for strCrd in match ]  
127   
128   
129 -    def nextLine( self ): 
130          """ 
131          extract next 10 coordinates from crd file 
132   
133          @return: coordinates 
134          @rtype: [float]     
135          """ 
136          l = self.crd.readline() 
137          if l == '': 
138              raise EOFError('EOF') 
139   
140          return self.line2numbers( l ) 
141   
142   
143 -    def nextFrame( self ): 
144          """ 
145          Collect next complete coordinate frame 
146   
147          @return: coordinate frame 
148          @rtype: array 
149          """ 
150   
151          i = 0 
152          xyz = [] 
153          while i != self.lines_per_frame: 
154   
155              if self.box and i+1 == self.lines_per_frame: 
156   
157                  ## skip box info 
158                  if len( self.nextLine() ) != 3: 
159                      raise ParseError( "BoxInfo must consist of 3 numbers." ) 
160   
161              else: 
162                  xyz += self.nextLine() 
163   
164              i += 1 
165   
166          return N.reshape( xyz, ( len(xyz) / 3, 3 ) ).astype(N.Float32) 
167   
168   
169 -    def crd2traj( self ): 
170          """ 
171          Convert coordinates into a Trajectory object. 
172   
173          @return: trajectory object 
174          @rtype: Trajectory 
175          """ 
176          ## skip first empty line 
177          self.crd.readline() 
178   
179          xyz = [] 
180          i = 0 
181   
182          if self.verbose: self.log.write( "Reading frames .." ) 
183   
184          try: 
185              while 1==1: 
186   
187                  xyz += [ self.nextFrame() ] 
188                  i += 1 
189   
190                  if i % 100 == 0 and self.verbose: 
191                      self.log.write( '#' ) 
192   
193          except EOFError: 
194              if self.verbose: self.log.add("Read %i frames." % i) 
195   
196          t = Trajectory( refpdb=self.ref ) 
197   
198          t.frames = N.array( xyz ).astype(N.Float32) 
199   
200          t.setRef( self.ref ) 
201          t.ref.disconnect() 
202   
203          return t 
204   
205  import Biskit.test as BT 
206  import tempfile 
207   
208 -class Test( BT.BiskitTest ): 
209      """Test AmberCrdParser""" 
210   
211 -    def prepare(self): 
212          root = T.testRoot() + '/amber/' 
213          self.finp = root + 'sim.crd.gz' 
214          self.fref = root + '1HPT_0.pdb' 
215          self.fout = tempfile.mktemp('.dat', 'traj_') 
216   
217 -    def cleanUp(self): 
218          T.tryRemove( self.fout ) 
219           
220   
221 -    def test_AmberCrdParser(self): 
222          """AmberCrdParser test""" 
223           
224          self.p = AmberCrdParser( self.finp, self.fref, box=True, rnAmber=True, 
225                                   log=self.log, verbose=self.local ) 
226          self.t = self.p.crd2traj() 
227   
228          self.t.removeAtoms(lambda a: a['residue_name'] in ['WAT','Na+','Cl-'] ) 
229          self.t.removeAtoms(lambda a: a['element'] == 'H' ) 
230   
231          if self.local: 
232              print "Dumping result to ", self.fout 
233   
234          T.Dump( self.t, T.absfile(self.fout) ) 
235   
236          if self.local: 
237              print "Dumped Trajectory with %i frames and %i atoms." % \ 
238                    (len(self.t), self.t.lenAtoms() ) 
239   
240          self.assertEqual( len(self.t), 10 ) 
241          self.assertEqual( self.t.lenAtoms(), 440 ) 
242   
243  if __name__ == '__main__': 
244   
245      BT.localTest() 
246