Source Code for Module gromacs.analysis.plugins.cysaccessibility

  1  # GromacsWrapper plugin: cysaccessibility.py 
  2  # Copyright (c) 2009 Oliver Beckstein <orbeckst@gmail.com> 
  3  # Released under the GNU Public License 3 (or higher, your choice) 
  4  # See the file COPYING for details. 
  5   
  6  """ 
  7  CysAccessibility plugin 
  8  ======================= 
  9   
 10  Cysteine accessibility is analyzed by histogramming the distance of 
 11  shortest approach of water molecules to the sulfhydryl group of Cys. 
 12   
 13  See class docs for more details. 
 14   
 15  .. note:: 
 16   
 17     This plugin is the canonical example for how to structure plugins that 
 18     conform to the plugin API (see docs :mod:`gromacs.analysis.core` for 
 19     details). 
 20   
 21  Plugin class 
 22  ------------ 
 23   
 24  .. autoclass:: CysAccessibility 
 25     :members: worker_class 
 26     :undoc-members: 
 27   
 28  Worker class 
 29  ------------ 
 30   
 31  The worker class performs the analysis. 
 32   
 33  .. autoclass:: _CysAccessibility 
 34     :members: 
 35   
 36   
 37  """ 
 38  from __future__ import with_statement 
 39   
 40  __docformat__ = "restructuredtext en" 
 41   
 42  import sys 
 43  import os.path 
 44  import warnings 
 45  import subprocess 
 46   
 47  import gromacs 
 48  from gromacs.utilities import AttributeDict 
 49  from gromacs.analysis.core import Worker, Plugin 
 50  import dist 
 51   
 52   
 53  # Worker classes that are registered via Plugins (see below) 
 54  # ---------------------------------------------------------- 
 55  # These must be defined before the plugins. 
 56   
 57 -class _CysAccessibility(Worker): 
 58      """Analysis of Cysteine accessibility.""" 
 59   
 60 -    def __init__(self,**kwargs): 
 61          """Set up  customized Cysteine accessibility analysis. 
 62   
 63          :Arguments: 
 64             cysteines : list 
 65                list of *all* resids (eg from the sequence) that are used as 
 66                labels or in the form 'Cys<resid>'. (**required**) 
 67             cys_cutoff : number 
 68                cutoff in nm for the minimum S-OW distance [1.0] 
 69   
 70          Note that *all* Cys residues in the protein are analyzed. Therefore, 
 71          the list of cysteine labels *must* contain as many entries as there are 
 72          cysteines in the protein. There are no sanity checks. 
 73          """ 
 74          # specific arguments: take them before calling the super class that 
 75          # does not know what to do with them 
 76          cysteines = kwargs.pop('cysteines',None)     # sequence resids as labels (NOT necessarily Gromacs itp) 
 77          cys_cutoff = kwargs.pop('cys_cutoff', 1.0)   # nm 
 78   
 79          # super class init: do this before doing anything else 
 80          # (also sets up self.parameters and self.results) 
 81          super(_CysAccessibility,self).__init__(**kwargs) 
 82   
 83          # process specific parameters now 
 84          try: 
 85              self.parameters.cysteines = map(int, cysteines)  # sequence resids 
 86          except (TypeError,ValueError): 
 87              raise ValueError("Keyword argument cysteines MUST be set to sequence of resids.") 
 88          self.parameters.cysteines.sort()                 # sorted because make_ndx returns sorted order 
 89          self.parameters.cutoff = cys_cutoff 
 90   
 91          if not self.simulation is None: 
 92              self._register_hook() 
 93   
 94 -    def _register_hook(self, **kwargs): 
 95          """Run when registering; requires simulation.""" 
 96   
 97          super(_CysAccessibility, self)._register_hook(**kwargs) 
 98          assert not self.simulation is None 
 99   
100          # filename of the index file that we generate for the cysteines 
101          self.parameters.ndx = self.plugindir('cys.ndx') 
102          # output filenames for g_dist, indexed by Cys resid 
103          self.parameters.filenames = dict(\ 
104              [(resid, self.plugindir('Cys%d_OW_dist.txt.bz2' % resid)) 
105               for resid in self.parameters.cysteines]) 
106          # default filename for the combined plot 
107          self.parameters.figname = self.figdir('mindist_S_OW') 
108   
109   
110      # override 'API' methods of base class 
111           
112 -    def run(self, cutoff=None, force=False, **gmxargs): 
113          """Run ``g_dist -dist cutoff`` for each cysteine and save output for further analysis. 
114   
115          By default existing data files are *not* overwritten. If this 
116          behaviour is desired then set the *force* = ``True`` keyword 
117          argument. 
118   
119          All other parameters are passed on to :func:`gromacs.g_dist`. 
120          """ 
121   
122          if cutoff is None: 
123              cutoff = self.parameters.cutoff 
124          else: 
125              self.parameters.cutoff = cutoff    # record cutoff used 
126   
127          ndx = self.parameters.ndx     # could move this into _register_hook? 
128          if not os.path.isfile(ndx): 
129              warnings.warn("Cysteine index file %r missing: running 'make_index_cys'." % ndx) 
130              self.make_index_cys() 
131   
132          for resid in self.parameters.cysteines: 
133              groupname = 'Cys%(resid)d' % vars() 
134              commands = [groupname, 'OW'] 
135              filename = self.parameters.filenames[resid] 
136              if not force and self.check_file_exists(filename, resolve='warning'): 
137                  continue 
138              print "run_g_dist: %(groupname)s --> %(filename)r" % vars() 
139              sys.stdout.flush() 
140              with open(filename, 'w') as datafile: 
141                  p = gromacs.g_dist.Popen( 
142                      s=self.simulation.tpr, f=self.simulation.xtc, n=ndx, dist=cutoff, input=commands,  
143                      stderr=None, stdout=subprocess.PIPE, **gmxargs) 
144                  compressor = subprocess.Popen(['bzip2', '-c'], stdin=p.stdout, stdout=datafile) 
145                  p.communicate() 
146   
147 -    def analyze(self,**kwargs): 
148          """Mindist analysis for all cysteines. Returns results for interactive analysis."""         
149   
150          results = AttributeDict() 
151          for resid in self.parameters.cysteines: 
152              groupname = 'Cys%(resid)d' % vars()    # identifier should be a valid python variable name 
153              results[groupname] = self._mindist(resid) 
154          self.results = results 
155          return results 
156   
157 -    def plot(self, **kwargs): 
158          """Plot all results in one graph, labelled by the result keys. 
159   
160          :Keywords: 
161             figure 
162                 - ``True``: save figures in the given formats 
163                 - "name.ext": save figure under this filename (``ext`` -> format) 
164                 - ``False``: only show on screen 
165             formats : sequence 
166                 sequence of all formats that should be saved [('png', 'pdf')] 
167             plotargs     
168                 keyword arguments for pylab.plot() 
169          """ 
170   
171          import pylab 
172          figure = kwargs.pop('figure', False) 
173          extensions = kwargs.pop('formats', ('pdf','png')) 
174          for name,result in self.results.items(): 
175              kwargs['label'] = name 
176              result.plot(**kwargs) 
177          pylab.legend(loc='best') 
178          if figure is True: 
179              for ext in extensions: 
180                  self.savefig(ext=ext) 
181          elif figure: 
182              self.savefig(filename=figure) 
183   
184       
185      # specific methods 
186   
187 -    def make_index_cys(self): 
188          """Make index file for all cysteines and water oxygens.  
189   
190          **NO SANITY CHECKS**: The SH atoms are simply labelled consecutively 
191          with the resids from the cysteines parameter. 
192          """ 
193          commands_1 = ['keep 0', 'del 0', 'r CYSH & t S', 'splitres 0', 'del 0']  # CYS-S sorted by resid 
194          commands_2 = ['t OW', 'q']                                               # water oxygens 
195          commands = commands_1[:] 
196          for groupid, resid in enumerate(self.parameters.cysteines): 
197              commands.append('name %(groupid)d Cys%(resid)d'  % vars())           # name CYS-S groups canonically 
198          commands.extend(commands_2) 
199          return gromacs.make_ndx(f=self.simulation.tpr, o=self.parameters.ndx,  
200                                  input=commands, stdout=None) 
201   
202 -    def _mindist(self,resid): 
203          """Analyze minimum distance for resid.""" 
204          filename = self.parameters.filenames[resid] 
205          return dist.Mindist(filename,cutoff=self.parameters.cutoff) 
206   
207   
208  # Public classes that register the worker classes 
209  #------------------------------------------------ 
210   
211 -class CysAccessibility(Plugin): 
212      """*CysAccessibility* plugin. 
213       
214      For each frame of a trajectory, the shortest distance of all water oxygens 
215      to all cysteine sulphur atoms is computed. For computational efficiency, 
216      only distances smaller than a cutoff are taken into account. A histogram of 
217      the distances shows how close water molecules can get to cysteines. The 
218      closest approach distance should be indicative of the reactivity of the SH 
219      group with crosslinking agents. 
220   
221      .. class:: CysAccessibility(cysteines, [cys_cutoff[, name[, simulation]]]) 
222       
223      :Arguments: 
224          name : string 
225              plugin name (used to access it) 
226          simulation : instance 
227              The :class:`gromacs.analysis.Simulation` instance that owns the plugin. 
228          cysteines : list 
229              list of *all* resids (eg from the sequence) that are used as 
230              labels or in the form 'Cys<resid>'. (**required**) 
231          cys_cutoff : number 
232              cutoff in nm for the minimum S-OW distance [1.0] 
233   
234      Note that *all* Cys residues in the protein are analyzed. Therefore, the 
235      list of cysteine labels *must* contain as many entries as there are 
236      cysteines in the protein. There are no sanity checks. 
237   
238      """ 
239      worker_class = _CysAccessibility 
240