1
2
3
4
5
6 """
7 Helix bundle analysis
8 =====================
9
10 Analysis of helix bundles with :func:`gromacs.g_bundle`.
11
12 .. SeeAlso:: HELANAL does some things better than g_bundle, in
13 particular it can find kinks automatically.
14
15 Plugin class
16 ------------
17
18 .. autoclass:: HelixBundle
19 :members: worker_class
20 :undoc-members:
21
22 Example
23 ~~~~~~~
24
25 .. Warning:: Note that the helices should be sorted in ascending order
26 because all indices are sorted that way. If they are in
27 mixed order then the correspondence between the name
28 column and the helices is lost. This is a shortcoming of
29 :program:`make_ndx` that is not easy to rectify.
30
31 Define helices in a reST table (use column names as shown!), set up
32 the plugin, and add it to a simulation instance::
33
34 helixtable = '''
35 Table[helices]: helices top = N-term, bot = C-term; kink
36 ===== ================ ================ ===============
37 name top bottom kink
38 ===== ================ ================ ===============
39 TM1 G28 P29 F30 T53 S54 S55 I41 Q42 V43
40 TM2 Q57 V58 W59 I84 R85 W86 F78 T79 Q80
41 TM3 R101 G102 S103 L135 L136 T137 F116 W117 F118
42 TM4 L142 P143 L144 T157 F158 Y159 F149 G150 A151
43 TM6 F209 S210 T211 D229 I230 V231 G219 W220 I221
44 TM7 E243 G244 Q245 L277 V278 G279 V261 P262 A263
45 TM8 P297 M298 A299 C327 S328 T329 N314 P315 A316
46 TM9 F336 K337 T338 L348 L349 M350 V342 S343 A344
47 ===== ================ ================ ===============
48 '''
49 hb = HelixBundle(helixtable=helixtable, offset=-9)
50 S = Simulation(....)
51 S.add_plugin(hb)
52
53
54 Worker class
55 ------------
56
57 The worker class performs the analysis.
58
59 .. autoclass:: _HelixBundle
60 :members:
61
62
63 """
64 from __future__ import with_statement
65
66 __docformat__ = "restructuredtext en"
67
68 import os.path
69 import warnings
70
71 import recsql
72
73 import gromacs
74 from gromacs.utilities import AttributeDict
75 from gromacs.analysis.core import Worker, Plugin
76
77 import logging
78 logger = logging.getLogger('gromacs.analysis.plugins.helixbundle')
79
80
81
82
83
85 """HelixBundle worker class."""
86
88 """Set up HelixBundle analysis.
89
90 :Keywords:
91 *helixtable*
92 reST table with columns "name", "top", "bottom",
93 "kink"; see :mod:`gromacs.analysis.plugins.helixbundle`
94 for details
95 *offset*
96 add the *offset* to the residue numbers in *helixtable* [0]
97 *helixndx*
98 provide a index file with the appropriate groups
99 instead of the table; also requires *na*
100 *na*
101 number of helices
102 *with_kinks*
103 take kinks into account [True]
104 *name*
105 plugin name [HelixBundle]
106 *simulation*
107 The :class:`gromacs.analysis.Simulation` instance that
108 owns the plugin [None]
109 """
110 helixndx = kwargs.pop('helixndx', None)
111 helixtable = kwargs.pop('helixtable', None)
112 na = kwargs.pop('na', None)
113 offset = kwargs.pop('offset', 0)
114 with_kinks = kwargs.pop('with_kinks', True)
115
116 if helixndx is None:
117 if helixtable is None:
118 raise ValueError("HelixBundle: requires a table *helixtable* with residues "
119 "that indicate the helix boundaries and kink or a index "
120 "file *helixndx* appropriate for g_bundle. See the "
121 "docs for details.")
122 self.helices = recsql.rest_table.Table2array(helixtable, autoconvert=True).recarray()
123 na = len(self.helices)
124 elif na is None:
125 raise ValueError("When using *helixndx* one must also provide the total number of "
126 "helices *na*. See g_bundle docs for details.")
127
128 super(_HelixBundle, self).__init__(**kwargs)
129
130 self.parameters.na = na
131 self.parameters.offset = offset
132 self.parameters.with_kinks = with_kinks
133
134 if not self.simulation is None:
135 self._register_hook()
136
138 """Run when registering; requires simulation."""
139
140 super(_HelixBundle, self)._register_hook(**kwargs)
141 assert not self.simulation is None
142
143 self.helixndx = self.plugindir('helices.ndx')
144
145 self.parameters.filenames = {
146 'length': self.plugindir('len.xvg'),
147 'distance': self.plugindir('dist.xvg'),
148 'z': self.plugindir('z.xvg'),
149 'tilt': self.plugindir('tilt.xvg'),
150 'tilt_radial': self.plugindir('tilt_radial.xvg'),
151 'tilt_lateral': self.plugindir('tilt_lateral.xvg'),
152 }
153 if self.parameters.with_kinks:
154 self.parameters.filenames.update(
155 {'kink': self.plugindir('kink.xvg'),
156 'kink_radial': self.plugindir('kink_radial.xvg'),
157 'kink_lateral': self.plugindir('kink_lateral.xvg'),
158 })
159
160
161 self.parameters.fignames = {
162 'z': self.figdir('z'),
163 'tilt': self.figdir('tilt'),
164 'kink': self.figdir('kink'),
165 }
166
168 """Build g_bundle index file from a record array.
169
170 :Keywords:
171 *force*
172 - ``None`` does nothing if the index file already exists
173 - ``False`` raises exception if it exists
174 - ``True`` always rebuild index file
175
176 Format of the index table:
177 - must contain columns name, top, bottom, kink
178 - a row corresponds to one helis
179 - name contains the name of the helix e.g. 'TM5'
180 - a entry is a *string* of residues which is split on white space;
181 entries in the same column must have the same number of residues
182 (limitation of g_bundle)
183
184 """
185
186 if self.check_file_exists(self.helixndx, resolve="indicate", force=force):
187 logger.warning("make_index(): The index file %r already exists; "
188 "use force = True to rebuild", self.helixndx)
189 return self.helixndx
190
191 logger.info("Making index file for %d helices %r", len(self.helices), self.helixndx)
192 logger.debug("make_index: tpr = %r", self.simulation.tpr)
193 logger.debug("make_index: offset = %r", self.parameters.offset)
194
195
196 tops = []
197 bottoms = []
198 kinks = []
199
200 for h in self.helices:
201
202 tops.extend(h.top.split())
203 bottoms.extend(h.bottom.split())
204 kinks.extend(h.kink.split())
205
206 tpr = self.simulation.tpr
207 offset = self.parameters.offset
208
209
210
211 Itop = gromacs.cbook.IndexBuilder(tpr, tops, offset=offset, out_ndx='top.ndx')
212 name, topndx = Itop.combine(name_all='top', defaultgroups=False)
213
214 Ibot = gromacs.cbook.IndexBuilder(tpr, bottoms, offset=offset, out_ndx='bottom.ndx')
215 name, botndx = Ibot.combine(name_all='bottom', defaultgroups=False)
216
217 Ikink = gromacs.cbook.IndexBuilder(tpr, kinks, offset=offset, out_ndx='kinks.ndx')
218 name, kinkndx = Ikink.combine(name_all='kink', defaultgroups=False)
219
220 Iall = gromacs.cbook.IndexBuilder(tpr, ndx=[topndx, botndx, kinkndx])
221 Iall.cat(self.helixndx)
222
223 return self.helixndx
224
225
226
227 - def run(self, force=None, **gmxargs):
228 """Analyze trajectory and write HelixBundle files.
229
230 :Arguments:
231 - *force*: ``True`` does analysis and overwrites existing files
232 - *gmxargs*: additional keyword arguments for :func:`gromacs.g_bundle`
233
234 .. Note:: The plugin default is *z* = ``True``, i.e. the tilt is computed
235 relative to the box z-axis.
236 """
237 gmxargs.setdefault('z', True)
238 gmxargs['na'] = self.parameters.na
239
240 if self.check_file_exists(self.parameters.filenames['tilt'], resolve='warning', force=force):
241 return
242
243 self.make_index(force=force)
244
245 logger.info("Analyzing HelixBundle (with_kinks=%r)...", self.parameters.with_kinks)
246 f = self.parameters.filenames
247 if self.parameters.with_kinks:
248 gromacs.g_bundle(s=self.simulation.tpr, f=self.simulation.xtc, n=self.helixndx,
249 ol=f['length'], od=f['distance'], oz=f['z'],
250 ot=f['tilt'], otr=f['tilt_radial'], otl=f['tilt_lateral'],
251 ok=f['kink'], okr=f['kink_radial'], okl=f['kink_lateral'],
252 input=['top','bottom', 'kink'],
253 **gmxargs)
254 else:
255 gromacs.g_bundle(s=self.simulation.tpr, f=self.simulation.xtc, n=self.helixndx,
256 ol=f['length'], od=f['distance'], oz=f['z'],
257 ot=f['tilt'], otr=f['tilt_radial'], otl=f['tilt_lateral'],
258 input=['top','bottom'],
259 **gmxargs)
260
261
263 """Collect output xvg files as :class:`gromacs.formats.XVG` objects.
264
265 :Returns: a dictionary of the results and also sets ``self.results``.
266 """
267 from gromacs.formats import XVG
268
269 logger.info("Preparing HelixBundle graphs as XVG objects.")
270 results = AttributeDict( (k, XVG(fn)) for k,fn in self.parameters.filenames.items() )
271 self.results = results
272 return results
273
274 - def plot(self, **kwargs):
275 """Plot all results in one graph, labelled by the result keys.
276
277 :Keywords:
278 figure
279 - ``True``: save figures in the given formats
280 - "name.ext": save figure under this filename (``ext`` -> format)
281 - ``False``: only show on screen
282 formats : sequence
283 sequence of all formats that should be saved [('png', 'pdf')]
284 plotargs
285 keyword arguments for pylab.plot()
286 """
287
288 import pylab
289 figure = kwargs.pop('figure', False)
290 extensions = kwargs.pop('formats', ('pdf','png'))
291 for name,result in self.results.items():
292 kwargs['label'] = name
293 try:
294 result.plot(**kwargs)
295 except AttributeError:
296 warnings.warn("Sorry, plotting of result %(name)r is not implemented" % vars(),
297 category=UserWarning)
298 pylab.legend(loc='best')
299 if figure is True:
300 for ext in extensions:
301 self.savefig(ext=ext)
302 elif figure:
303 self.savefig(filename=figure)
304
305
306
307
308
309
310
312 """*HelixBundle* plugin.
313
314 :func:`gromacs.g_bundle` helix analysis
315
316 .. class:: HelixBundle([helixtable, offset, with_kinks, [name[, simulation]]])
317
318 :Arguments:
319 *helixtable*
320 reST table with columns "name", "top", "bottom",
321 "kink"; see :mod:`gromacs.analysis.plugins.helixbundle`
322 for details
323 *offset*
324 add the *offset* to the residue numbers in *helixtable* [0]
325 *helixndx*
326 provide a index file with the appropriate groups
327 instead of the table; also requires *na*
328 *na*
329 number of helices
330 *with_kinks*
331 take kinks into account [True]
332 *name*
333 plugin name [HelixBundle]
334 *simulation*
335 The :class:`gromacs.analysis.Simulation` instance that owns
336 the plugin. [None]
337 """
338 worker_class = _HelixBundle
339