## "Fossies" - the Fresh Open Source Software Archive

### Member "relax-5.0.0/docs/latex/curvefit.tex" (18 Apr 2019, 44183 Bytes) of package /linux/privat/relax-5.0.0.src.tar.bz2:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) TeX and LaTeX source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %                                                                             %
3 % Copyright (C) 2006,2008,2012-2014 Edward d'Auvergne                         %
4 %                                                                             %
5 % This file is part of the program relax (http://www.nmr-relax.com).          %
6 %                                                                             %
7 % This program is free software: you can redistribute it and/or modify        %
9 % the Free Software Foundation, either version 3 of the License, or           %
10 % (at your option) 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               %
15 % GNU General Public License for more details.                                %
16 %                                                                             %
17 % You should have received a copy of the GNU General Public License           %
18 % along with this program.  If not, see <http://www.gnu.org/licenses/>.       %
19 %                                                                             %
20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21
22
23 % Relaxation curve-fitting.
24 %%%%%%%%%%%%%%%%%%%%%%%%%%%
25
26 \chapter[Relaxation curve-fitting]{The $\Rone$ and $\Rtwo$ relaxation rates -- relaxation curve-fitting} \label{ch: relax-fit}
27 \index{relaxation curve-fitting|textbf}
28
29
30 \begin{figure*}[h]
31   \includegraphics[width=5cm, bb=0 0 1701 1701]{graphics/analyses/r1_600x600} \hfill \includegraphics[width=5cm, bb=0 0 1701 1701]{graphics/analyses/r2_600x600}
32 \end{figure*}
33
34
35 % Introduction.
36 %%%%%%%%%%%%%%%
37
38 \section{Introduction to relaxation curve-fitting}
39
40 The fitting of exponentials to relaxation curves (relaxation curve-fitting or as used throughout this chapter abbreviated simply as relax-fit) involves a number of steps including the loading of data, the calculation of both the average peak intensity\index{peak!intensity} across replicated spectra and the standard deviations\index{standard deviation} of those peak intensities, selection of the experiment type, optimisation of the parameters of the exponential curves during the fit for each observed spin, Monte Carlo simulations\index{Monte Carlo simulation} to find the parameter errors, and saving and viewing the results.
41 To simplify the process a sample script will be followed step by step as was done with the NOE calculation.
42
43
44
45 % The models.
46 %%%%%%%%%%%%%
47
48 \section{The exponential curve models}
49 \label{sect: exponential curve models}
50
51 A number of different models are supported in this analysis.
52 These include the two parameter exponential decay to zero, the inversion recovery experiment, and the saturation recovery experiment.
53 These can be selected using the \uf{relax\ufus{}fit\ufsep{}select\ufus{}model} user function.
54
55 The default is the two parameter exponential decay whereby the magnetisation starts at $I_0$ and decays to zero.
56 It has the parameters \{$\mathrm{R}_x$, $I_0$\}.
57 The formula of this function is
58 \begin{equation}
59   I(t) = I_0 e^{-\mathrm{R}_x \cdot t} ,
60 \end{equation}
61
62 \noindent where $I(t)$ is the peak intensity at any time point $t$, $I_0$ is the initial intensity, and $\mathrm{R}_x$ is the relaxation rate (either the $\Rone$ or $\Rtwo$).
63
64 In the inversion recovery experiment, the magnetisation starts at a negative value at $-I_0$ and relaxes to a positive $I_{\infty}$ value.
65 This curve consists of three parameters \{$\mathrm{R}_x$, $I_0$, $I_{\infty}$\}.
66 The formula is
67 \begin{equation}
68   I(t) = I_{\infty} - I_0 e^{-\mathrm{R}_x \cdot t} .
69 \end{equation}
70
71 In the saturation recovery experiment, the magnetisation starts at zero and relaxes to a positive $I_{\infty}$ value.
72 The model consists of the two parameters \{$\mathrm{R}_x$, $I_{\infty}$\} and has the formula
73 \begin{equation}
74   I(t) = I_{\infty} \left( 1 - e^{-\mathrm{R}_x \cdot t} \right) .
75 \end{equation}
76
77
78
79 % From spectra to peak intensities.
80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81
82 \section{From spectra to peak intensities for the relaxation rates} \label{sect: spectra to intensities}
83
84 The following subsections simply contain advice on how to go from the recorded FIDs to the peak lists ready to be input into relax.
85 This need not be followed -- it is simply a set of recommendations for obtaining the highest quality relaxation rates.
86
87
88 % Temperature control and calibration.
89 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
90
91 \subsection{Temperature control and calibration} \label{sect: temperature control and calibration}
92
93 \includegraphics[bb=0 0 18 18]{graphics/oxygen_icons/128x128/status/weather-clear}
94
95 Before starting with the spectral processing, it should be noted that proper temperature control and calibration are essential for relaxation data.
96 Small temperature changes can have an effect on the viscosity and hence global tumbling of the molecule being studied and, as the molecular diffusion tensor is the major contributor to relaxation, any non-consistent data will likely lead to artificial motions appearing in subsequent model-free analyses.
97
98 Per-experiment temperature calibration is essential and the technique used will need to be specified for BMRB data deposition.
99 Note that the standard MeOH/ethylene glycol calibration of a spectrometer is of no use when you are running experiments which pump in large amounts of power into the probe head.
100 Although the R1 experiment should be about the same temperature as a HSQC and hence be close to the standard MeOH/ethylene glycol spectrometer calibration, the R2 CPMG or spin lock and, to a lesser extent, the NOE pre-saturation pump a lot more power into the probe head.
101 The power differences can either cause the temperature in the sample to be too high or too low.
102 This is unpredictable as the thermometer used by the VT unit is next to the coils in the probe head and not inside the NMR sample.
103 So the VT unit tries to control the temperature inside the probe head rather than in the NMR sample.
104 However between the thermometer and the sample is the water of the sample, the glass of the NMR tube, the air gap where the VT unit controls air flow and the outside components of the probe head protecting the electronics.
105 If the sample, the probe head or the VT unit is changed, this will have a different affect on the per-experiment temperature.
106 The VT unit responds differently under different conditions and may sometimes over or under compensate by a couple of degrees.
107 Therefore each relaxation data set from each spectrometer requires a per-experiment calibration.
108
109 Explicit temperature control techniques are also essential for relaxation data collection.
110 Again the technique used will be asked for by relax for BMRB data deposition.
111 A number of factors can cause significant temperature fluctuations between individual relaxation experiments.
112 This includes the daily temperature cycle of the room housing the spectrometer, different amounts of power for the individual experiments, etc.
113 The best methods for eliminating such problems are single scan interleaving and temperature compensation block.
114 Single scan interleaving is the most powerful technique for averaging the temperature fluctuations not only across different experiments, but also across the entire measurement time.
115 The application of off-resonance temperature compensation blocks at the start of the experiment is useful for the R2 and will normalise the temperature between the individual experiments, but single scan or single fid interleaving is nevertheless required for normalising the temperature across the entire measurement.
116
117
118 % Spectral processing.
119 %~~~~~~~~~~~~~~~~~~~~~
120
121 \subsection{Spectral processing}
122
123 For the best measurement of peak heights across the myriad of NMR spectral analysis software, it is recommend to zero fill a lot -- 8k to 16k would give the best results.
124 This does not increase the information content of the spectrum or decrease the errors, it simply interpolates.
125 Even if the NMR spectral software performs 3-point quadratic interpolation between the highest points to determine the peak height, the additional free interpolation will make the estimation more accurate.
126
127 Additionally, care must be taken to properly scale the first point as this can cause a baseline roll which will affect peak heights.
128 A very useful description comes directly from the \href{http://spin.niddk.nih.gov/NMRPipe/doc1/}{NMRPipe manual}:
129
130 \begin{quotation}
131 Depending on the delay, the first point of the FID should be adjusted before Fourier Transform.  The first point scaling factor is selected by the window function argument \prompt{-c}.
132
133 If the required first order phase P1 for the given dimension is 0.0, the first point scaling factor should be 0.5.  This is because the discrete Fourier transform does the equivalent of counting the point at t=0 twice.  If the first point is not scaled properly in this case, ridge-line baseline offsets in the spectrum will result.
134
135 In all other cases (P1 is not zero), this scale factor should be 1.0. This is because the first point of the FID no longer corresponds to t=0, and so it shouldn't be scaled. If the scale factor is not set correctly, it will introduce a baseline distortion which is either zero-order or sinusoidal, depending on what first-order phase is required. When possible, it is best to set up experiments with either exactly 0, 1/2, or 1-point delay.  There are several reasons:
136 \begin{itemize}
137   \item Phase correction values can be determined easily.
138   \item If the delay is not a multiple of 1/2 point, the phase of folded peaks will be distorted.
139   \item The Hilbert transform (HT) is used, sometimes automatically, to reconstruct previously deleted imaginary data for interactive rephasing or inverse processing. But, the HT can only reconstruct imaginary data perfectly if the phase is a multiple of 1/2 point.
140   \item Data with P1 = 360 have the first point t=0 missing (i.e.\ 1 point delay). Since the first point of the FID corresponds to the sum of points in the corresponding spectrum, this missing first point can be restored'' by adding a constant to the phased spectrum.  This can be done conveniently by automated zero-order baseline correction, as shown in table~\ref{table: NMRPipe -c}.
141
142     \begin{table}
143     \begin{center}
144     \caption{Summary, First Point Scaling and Phase Correction}
145     \begin{tabular}{llll}
146     \toprule
147     Delay & P1 & FID & Spectrum\\
148     \midrule
149     0   point &   0 & Scale -c 0.5 \\
150     1/2 point & 180 & Scale -c 1.0 & Folded peaks have opposite sign \\
151     1   point & 360 & Scale -c 1.0 & Use POLY -auto -ord 0'' \\
152     \bottomrule
153     \label{table: NMRPipe -c}
154     \end{tabular}
155     \end{center}
156     \end{table}
157
158 \end{itemize}
159 \end{quotation}
160
161
162 Here is an example NMRPipe script designed for optimal relaxation rate extraction:
163 \begin{lstlisting}[language=csh]
164 #!/bin/csh
165
166 setenv FILEROOT $1 167 set PHASE=81.4 168 169 echo "\n# Fourier Transform (nmrPipe fid/*.fid to ft/*.dat)" 170 echo "# t2 phase is set to$PHASE"
171 echo "# t1 phase is set to 0.0\n"
172
173 nmrPipe -in fid/$FILEROOT.fid \ 174 | nmrPipe -fn SOL \ 175 | nmrPipe -fn GM -g1 15 -g2 20 -c 0.5 \ 176 | nmrPipe -fn ZF -size 8192 \ 177 | nmrPipe -fn FT -auto \ 178 | nmrPipe -fn PS -p0$PHASE -p1 0.0 -di -verb \
179 | nmrPipe -fn TP \
180 | nmrPipe -fn SP -off 0.5 -end 0.98 -pow 2 -c 0.5 \
181 | nmrPipe -fn ZF -size 8192 \
182 | nmrPipe -fn FT -auto \
183 | nmrPipe -fn PS -p0 0.0 -p1 0.0 -di -verb \
184 | nmrPipe -fn TP \
185 | nmrPipe -fn POLY -auto \
186 | nmrPipe -fn EXT -left -sw \
187 | nmrPipe -out ft/$FILEROOT.dat -ov 188 \end{lstlisting} 189 190 The script is run by suppling the FILEROOT value as a command line option so if the script is called \file{nmrpipe.sh} and the \file{var2pipe} or \file{bruk2pipe} processed file \file{R1\osus{}ncyc4.fid} is in the \directory{fid} directory, you would run: 191 \begin{lstlisting}[language=bash,numbers=none] 192$ ./nmrpipe.sh R1_ncyc4
193 \end{lstlisting}
194
195 The \directory{ft} directory must exist for this script to execute.
196 Different experiment specific options may be needed such as:
197 \begin{lstlisting}[language=csh,numbers=none]
198 | nmrPipe -fn REV \
199 | nmrPipe -fn FT -neg \
200 | nmrPipe -fn PS -rs 2.5 \
201 \end{lstlisting}
202
203 The script should be changed for different phasing, first point scaling, a polynomial baseline correction added in the direct dimension or removed from the indirect dimension, solvent suppression removed or changed, and the window functions modified for optimal spectral quality.
204 Each system and spectrum is different, so it is recommended that to find the optimal processing that each part of the script be removed and re-added one-by-one between processing and checking of the resultant spectrum.
205 Note that the extraction at the end after the polynomial baseline correction in the indirect dimension is important as the baseline correction often displays a much better performance when the empty part of the spectrum is used in the calculation.
206
207
208
209 % Measuring peak intensities.
210 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211
212 \subsection{Measuring peak intensities}
213
214 For the measurement of peak intensities, again care must be taken.
215 A read of the paper:
216 \begin{itemize}
217   \item \bibentry{Viles01}
218 \end{itemize}
219
220 is highly recommended.
221 Despite the recommendations in the discussion of this paper, a different methodology using peak heights can be used to solve the same problems.
222 This will be discussed in a paper which is currently in preparation from the Gooley group.
223 The steps involved are:
224 \begin{itemize}
225   \item For the first spectrum in the time series, shift the peak list to the tops of the peaks (for example using \gui{pc} in Sparky on subsets of peaks).
226   \item Copy this \nth{1} spectrum list onto all spectra, shifting the peaks to the top as in the previous step.
227   \item When the peak disappears into the noise, leave it at its current position and do not type \gui{pc} or equivalent.
228     This will add weight to the first point in the subsequent step.
229   \item Once all spectra are shifted, calculate an average peak list.
230   \item Copy this average peak list onto fresh copies of all spectra.
231   \item Measure peak heights using this averaged peak list.
232 \end{itemize}
233
234 This will produce the most accurate peak intensity measurements until better, more robust peak shape integration comes along.
235 This is a special technique which is designed to minimise the white-noise bias talked about in the \citet{Viles01} paper.
236 As the noise often decreases with the decrease in total spectral power, using the tops of the peaks means that you are actually measuring the real peak height plus positive noise in all cases.
237 This non-constant additional positive noise contribution can result in a double exponential in the measured data.
238 The technique above eliminates this as you then measure close to real peak height with the addition of white noise centred at zero -- it is both negative and positive to equal amounts -- rather than the peak high with noise contribution strongly biased towards the positive.
239 Where the peaks disappear, you then are measuring the pure baseplane noise.
240 This is fine as these white-noise data points centred at zero will help in the subsequent exponential fit in relax.
241
242 If using Sparky then, to be sure that the peak heights are properly updated, for each spectrum type \gui{pa} to select all peaks, \gui{ph} to update all selected peak heights, \gui{lt} to show the spectrum peaks window, make sure \gui{data height} is selected in the options, and then save the peak list.
243
244
245
246 % Script UI.
247 %%%%%%%%%%%%
248 \section{Relaxation curve-fitting in the prompt/script UI mode}
249
250
251 % The sample script.
252 %~~~~~~~~~~~~~~~~~~~
253
254 \subsection{Relax-fit script mode -- the sample script}
255
256 The following is a verbatim copy of the contents of the \file{sample\osus{}scripts\ossep{}relax\osus{}fit.py} file.
257 If your copy of the sample script is different than that below, please send an email to the relax-devel mailing list\index{mailing list!relax-devel} to tell the relax developers that the manual is out of date (see section~\ref{sect: relax-devel mailing list} on page~\pageref{sect: relax-devel mailing list}).
258 You will need to first copy this script to a dedicated analysis directory containing peak lists, a PDB file and a file listing unresolved spin systems, and then modify its contents to suit your specific analysis.
259 The script contents are:
260
261 \begin{lstlisting}
262 # Script for relaxation curve-fitting.
263
264 # Create the 'rx' data pipe.
265 pipe.create('rx', 'relax_fit')
266
267 # Load the backbone amide 15N spins from a PDB file.
271
272 # Spectrum names.
273 names = [
274     'T2_ncyc1_ave',
275     'T2_ncyc1b_ave',
276     'T2_ncyc2_ave',
277     'T2_ncyc4_ave',
278     'T2_ncyc4b_ave',
279     'T2_ncyc6_ave',
280     'T2_ncyc9_ave',
281     'T2_ncyc9b_ave',
282     'T2_ncyc11_ave',
283     'T2_ncyc11b_ave'
284 ]
285
286 # Relaxation times (in seconds).
287 times = [
288     0.0176,
289     0.0176,
290     0.0352,
291     0.0704,
292     0.0704,
293     0.1056,
294     0.1584,
295     0.1584,
296     0.1936,
297     0.1936
298 ]
299
300 # Loop over the spectra.
301 for i in range(len(names)):
302     # Load the peak intensities.
304
305     # Set the relaxation times.
306     relax_fit.relax_time(time=times[i], spectrum_id=names[i])
307
308 # Specify the duplicated spectra.
309 spectrum.replicated(spectrum_ids=['T2_ncyc1_ave', 'T2_ncyc1b_ave'])
310 spectrum.replicated(spectrum_ids=['T2_ncyc4_ave', 'T2_ncyc4b_ave'])
311 spectrum.replicated(spectrum_ids=['T2_ncyc9_ave', 'T2_ncyc9b_ave'])
312 spectrum.replicated(spectrum_ids=['T2_ncyc11_ave', 'T2_ncyc11b_ave'])
313
314 # Peak intensity error analysis.
315 spectrum.error_analysis()
316
317 # Deselect unresolved spins.
318 deselect.read(file='unresolved', mol_name_col=1, res_num_col=2, res_name_col=3, spin_num_col=4, spin_name_col=5)
319
320 # Set the relaxation curve type.
321 relax_fit.select_model('exp')
322
323 # Grid search.
324 minimise.grid_search(inc=11)
325
326 # Minimise.
327 minimise.execute('newton', constraints=False)
328
329 # Monte Carlo simulations.
330 monte_carlo.setup(number=500)
331 monte_carlo.create_data()
332 monte_carlo.initial_values()
333 minimise.execute('newton', constraints=False)
334 monte_carlo.error_analysis()
335
336 # Save the relaxation rates.
337 value.write(param='rx', file='rx.out', force=True)
338
339 # Save the results.
340 results.write(file='results', force=True)
341
342 # Create Grace plots of the data.
343 grace.write(y_data_type='chi2', file='chi2.agr', force=True)    # Minimised chi-squared value.
344 grace.write(y_data_type='i0', file='i0.agr', force=True)    # Initial peak intensity.
345 grace.write(y_data_type='rx', file='rx.agr', force=True)    # Relaxation rate.
346 grace.write(x_data_type='relax_times', y_data_type='peak_intensity', file='intensities.agr', force=True)    # Average peak intensities.
347 grace.write(x_data_type='relax_times', y_data_type='peak_intensity', norm=True, file='intensities_norm.agr', force=True)    # Average peak intensities (normalised).
348
349 # Display the Grace plots.
350 grace.view(file='chi2.agr')
351 grace.view(file='i0.agr')
352 grace.view(file='rx.agr')
353 grace.view(file='intensities.agr')
354 grace.view(file='intensities_norm.agr')
355
356 # Save the program state.
357 state.save('rx.save', force=True)
358 \end{lstlisting}
359
360 The next sections will break this script down into its logical components and explain how these parts will be interpreted by relax.
361 To execute this script, please see section~\ref{sect: scripting} on page~\pageref{sect: scripting} for details.
362
363
364 % Initialisation of the data pipe.
365 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
366
367 \subsection{Relax-fit script mode -- initialisation of the data pipe} \label{Rx initialisation}
368
369 The data pipe is simply created by the command
370
371 \begin{lstlisting}[firstnumber=3]
372 # Create the 'rx' data pipe.
373 pipe.create('rx', 'relax_fit')
374 \end{lstlisting}
375
376 This user function will then create a relaxation exponential curve-fitting specific data pipe labelled \promptstring{rx}.
377 The second argument sets the pipe type to that of the relaxation curve-fitting.
378 Setting the pipe type is important so that the program knows which user functions are compatible with the data pipe, for example in the steady-state NOE analysis the function \uf{minimise.execute} (see page~\pageref{uf: minimise.execute}) is meaningless as the NOE values are calculated directly rather than optimised.
379
380
381 % Spin systems.
382 %~~~~~~~~~~~~~~
383
384 \subsection{Relax-fit script mode -- setting up the spin systems}
385
386 The first thing which needs to be completed prior to any spin specific command is to generate the molecule, residue and spin data structures for storing the spin specific data.
387 In the sample script above this is generated from a PDB file, however a plain text file with the sequence information can be used instead (see the \uf{sequence\ufsep{}read} user function on page~\pageref{uf: sequence.read} for more details).
388 In the case of the sample script, the command
389
390 \begin{lstlisting}[firstnumber=6]
391 # Load the backbone amide 15N spins from a PDB file.
393 \end{lstlisting}
394 \index{PDB}
395
396 will load the PDB file \file{Ap4Aase\osus{}new\osus{}3.pdb} into relax.
397 Then
398
399 \begin{lstlisting}[firstnumber=8]
402 \end{lstlisting}
403
404 will generate the molecule, residue, and spin sequence for the current data pipe.
405 In this situation there will be a single spin system per residue generated corresponding to the backbone amide nitrogens as well as $^{15}$N spins set up for the tryptophan indole nitrogens.
406 Although the 3D coordinates have been loaded into the program from the PDB file, this structural information serves no purpose when calculating $\Rone$ and $\Rtwo$ values.
407
408
410 %~~~~~~~~~~~~~~~~~~
411
413
414 To load the peak intensities\index{peak!intensity} into relax the \uf{spectrum\ufsep{}read\ufus{}intensities} and \uf{relax\ufus{}fit\ufsep{}relax\ufus{}times} user functions are executed.
415 Important keyword arguments for these user functions are the file name and directory, the spectrum identification string and the relaxation time period of the experiment in seconds.
416 By default the file format will be automatically detected.
417 Currently Sparky\index{software!Sparky}, XEasy\index{software!XEasy}, NMRView\index{software!NMRView}, and generic columnar formatted peak lists are supported.
418 To be able to import any other type of format please send an email to the relax development mailing list\index{mailing list!relax-devel} with the details of the format.
419 Adding support for new formats is trivial.
420 The following series of commands -- an expansion of the \prompt{for} loop in the sample script -- will load peak intensities from six different relaxation periods, four of which have been duplicated, from Sparky peak lists with the peak heights in the 10$^{\textrm{th}}$ column.
421
422 \begin{lstlisting}[numbers=none]
424 relax_fit.relax_time(spectrum_id='1',   time=0.0176)
426 relax_fit.relax_time(spectrum_id='1b',  time=0.0176)
428 relax_fit.relax_time(spectrum_id='2',   time=0.0352)
430 relax_fit.relax_time(spectrum_id='4',   time=0.0704)
432 relax_fit.relax_time(spectrum_id='4b',  time=0.0704)
434 relax_fit.relax_time(spectrum_id='6',   time=0.1056)
436 relax_fit.relax_time(spectrum_id='9',   time=0.1584)
438 relax_fit.relax_time(spectrum_id='9b',  time=0.1584)
440 relax_fit.relax_time(spectrum_id='11',  time=0.1936)
442 relax_fit.relax_time(spectrum_id='11b', time=0.1936)
443 \end{lstlisting}
444
445 The replicated spectra a set up with the commands
446
447 \begin{lstlisting}[firstnumber=47]
448 # Specify the duplicated spectra.
449 spectrum.replicated(spectrum_ids=['T2_ncyc1_ave', 'T2_ncyc1b_ave'])
450 spectrum.replicated(spectrum_ids=['T2_ncyc4_ave', 'T2_ncyc4b_ave'])
451 spectrum.replicated(spectrum_ids=['T2_ncyc9_ave', 'T2_ncyc9b_ave'])
452 spectrum.replicated(spectrum_ids=['T2_ncyc11_ave', 'T2_ncyc11b_ave'])
453 \end{lstlisting}
454
455 Note that the relaxation time period should be calculated directly from the pulse sequence (as the sum of delays and pulses for the period), as the estimated time may not match the real time.
456 For the Sparky peak lists, by default relax assumes that the intensity value is in the 4$^{\textrm{th}}$ column.
457 A typical file looks like:
458
459 {\scriptsize \begin{verbatim}
460      Assignment         w1         w2   Data Height
461
462         LEU3N-HN    122.454      8.397       129722
463         GLY4N-HN    111.999      8.719       422375
464         SER5N-HN    115.085      8.176       384180
465         MET6N-HN    120.934      8.812       272100
466         ASP7N-HN    122.394      8.750       174970
467         SER8N-HN    113.916      7.836       218762
468        GLU11N-HN    122.194      8.604        30412
469        GLY12N-HN    110.525      9.028        90144
470 \end{verbatim}}
471
472 By supplying the \keyword{int\_col} argument to the \uf{spectrum\ufsep{}read\ufus{}intensities} user function, this can be changed.
473 A typical XEasy file will look like:
474
475 {\scriptsize \begin{verbatim}
476  No.  Color    w1      w2     ass. in w1     ass. in w2    Volume     Vol. Err.  Method  Comment
477
478    2    2    10.014 134.221   HN  21 LEU      N  21 LEU    7.919e+03  0.00e+00     m
479    3    2    10.481 132.592  HE1  79 TRP    NE1  79 TRP    1.532e+04  0.00e+00     m
480   17    2     9.882 129.041   HN 110 PHE      N 110 PHE    9.962e+03  0.00e+00     m
481   18    2     8.757 128.278   HN  52 ASP      N  52 ASP    2.041e+04  0.00e+00     m
482   19    2    10.086 128.297   HN  69 SER      N  69 SER    9.305e+03  0.00e+00     m
483   20    3     9.111 127.707   HN  15 ARG      N  15 ARG    9.714e+03  0.00e+00     m
484 \end{verbatim}}
485
486 where the peak height is in the \prompt{Volume} column.
487 And for an NMRView file:
488
489 {\tiny \begin{verbatim}
490 label dataset sw sf
491 H1 N15
492 cNTnC_noe0.nv
493 2505.63354492 1369.33557129
494 499.875 50.658000946
495 H1.L H1.P H1.W H1.B H1.E H1.J H1.U N15.L N15.P N15.W N15.B N15.E N15.J N15.U vol int stat comment flag0
496 0 {70.HN} 10.75274 0.02954 0.05379 ++ 0.0 {} {70.N} 116.37241 0.23155 0.35387 ++ 0.0 {} -6.88333129883 -0.1694 0 {} 0
497 1 {72.HN} 9.67752 0.03308 0.05448 ++ 0.0 {} {72.N} 126.41302 0.27417 0.37217 ++ 0.0 {} -5.49038267136 -0.1142 0 {} 0
498 2 {} 8.4532 0.02331 0.05439 ++ 0.0 {} {} 122.20137 0.38205 0.33221 ++ 0.0 {} -2.58034267191 -0.1320 0 {} 0
499 \end{verbatim}}
500
501
502 % The rest of the setup.
503 %~~~~~~~~~~~~~~~~~~~~~~~
504
505 \subsection{Relax-fit script mode -- the rest of the setup} \label{sect: Rx setup fin}
506
507 Once all the peak intensity data has been loaded a few calculations are required prior to optimisation.
508 Firstly the peak intensities for individual spins needs to be averaged across replicated spectra.
509 The peak intensity errors also have to be calculated using the standard deviation formula.
510 These two operations are executed by the user function
511
512 \begin{lstlisting}[firstnumber=53]
513 # Peak intensity error analysis.
514 spectrum.error_analysis()
515 \end{lstlisting}
516
517 Any spins which cannot be resolved due to peak overlap were included in a file called \file{unresolved}.
518 This file can consist of optional columns of the molecule name, the residue name and number, and the spin name and number.
519 The matching spins are excluded from the analysis by the user function
520
521 \begin{lstlisting}[firstnumber=56]
522 # Deselect unresolved spins.
523 deselect.read(file='unresolved', mol_name_col=1, res_num_col=2, res_name_col=3, spin_num_col=4, spin_name_col=5)
524 \end{lstlisting}
525
526 Finally the experiment type is specified by the command
527
528 \begin{lstlisting}[firstnumber=59]
529 # Set the relaxation curve type.
530 relax_fit.select_model('exp')
531 \end{lstlisting}
532
533 The argument \promptstring{exp} sets the relaxation curve to a two parameter \{$\mathrm{R}_x$, $I_0$\} exponential which decays to zero.
534 Changing the user function argument to \promptstring{inv} will select the inversion recovery experiment, and changing it to \promptstring{sat} will select the saturation recovery experiment (see section~\ref{sect: exponential curve models} on page~\pageref{sect: exponential curve models}).
535
536
537
538 % Optimisation.
539 %~~~~~~~~~~~~~~
540
541 \subsection{Relax-fit script mode -- optimisation of exponential curves}
542
543 Now that everything has been setup minimisation can be used to optimise the parameter values.
544 Firstly a grid search is applied to find a rough starting position for the subsequent optimisation algorithm.
545 Eleven increments per dimension of the model (in this case the two dimensions \{$\mathrm{R}_x$, $I_0$\}) is sufficient.
546 The user function for executing the grid search is
547
548 \begin{lstlisting}[firstnumber=62]
549 # Grid search.
550 minimise.grid_search(inc=11)
551 \end{lstlisting}
552
553 The next step is to select one of the minimisation algorithms to optimise the model parameters
554
555 \begin{lstlisting}[firstnumber=65]
556 # Minimise.
557 minimise.execute('newton', constraints=False)
558 \end{lstlisting}
559
560
561
562 % Error analysis.
563 %~~~~~~~~~~~~~~~~
564
565 \subsection{Relax-fit script mode -- error analysis}
566
567 Only one technique adequately estimates parameter errors when the parameter values where found by optimisation -- Monte Carlo simulations\index{Monte Carlo simulation}.
568 In relax this can be implemented by using a series of functions from the \uf{monte\ufus{}carlo} user function class.
569 Firstly the number of simulations needs to be set
570
571 \begin{lstlisting}[firstnumber=68]
572 # Monte Carlo simulations.
573 monte_carlo.setup(number=500)
574 \end{lstlisting}
575
576 For each simulation, randomised relaxation curves will be fit using exactly the same methodology as the original exponential curves.
577 These randomised curves are created by back calculation from the fitted model parameter values and then each point on the curve randomised using the error values set earlier in the script
578
579 \begin{lstlisting}[firstnumber=70]
580 monte_carlo.create_data()
581 \end{lstlisting}
582
583 As a grid search for each simulation would be too computationally expensive, the starting point for optimisation for each simulation can be set to the position of the optimised parameter values of the model
584
585 \begin{lstlisting}[firstnumber=71]
586 monte_carlo.initial_values()
587 \end{lstlisting}
588
589 Then exactly the same optimisation as was used for the model can be performed
590
591 \begin{lstlisting}[firstnumber=72]
592 minimise.execute('newton', constraints=False)
593 \end{lstlisting}
594
595 The parameter errors are then determined as the standard deviation of the optimised parameter values of the simulations
596
597 \begin{lstlisting}[firstnumber=73]
598 monte_carlo.error_analysis()
599 \end{lstlisting}
600
601
602 % Finishing off.
603 %~~~~~~~~~~~~~~~
604
605 \subsection{Relax-fit script mode -- finishing off}
606
607 To finish off, the script first saves the relaxation rates together with their errors in a simple text file
608
609 \begin{lstlisting}[firstnumber=75]
610 # Save the relaxation rates.
611 value.write(param='rx', file='rx.out', force=True)
612 \end{lstlisting}
613
614 Grace plots are created and viewed
615
616 \begin{lstlisting}[firstnumber=81]
617 # Create Grace plots of the data.
618 grace.write(y_data_type='chi2', file='chi2.agr', force=True)    # Minimised chi-squared value.
619 grace.write(y_data_type='i0', file='i0.agr', force=True)    # Initial peak intensity.
620 grace.write(y_data_type='rx', file='rx.agr', force=True)    # Relaxation rate.
621 grace.write(x_data_type='relax_times', y_data_type='peak_intensity', file='intensities.agr', force=True)    # Average peak intensities.
622
623 grace.write(x_data_type='relax_times', y_data_type='peak_intensity', norm=True, file='intensities_norm.agr', force=True)    # Average peak intensities (normalised).
624 \end{lstlisting}
625
626 and viewed
627
628 \begin{lstlisting}[firstnumber=88]
629 # Display the Grace plots.
630 grace.view(file='chi2.agr')
631 grace.view(file='i0.agr')
632 grace.view(file='rx.agr')
633 grace.view(file='intensities.agr')
634 grace.view(file='intensities_norm.agr')
635 \end{lstlisting}
636
637 and finally the program state is saved for future reference
638
639 \begin{lstlisting}[firstnumber=95]
640 # Save the program state.
641 state.save(file='rx.save', force=True)
642 \end{lstlisting}
643
644
645
646 % GUI.
647 %%%%%%
648
649 \newpage
650 \section{The relaxation curve-fitting auto-analysis in the GUI}
651
652 The $\Rone$ and $\Rtwo$ relaxation rates can be calculated using the relax GUI (see Figures~\ref{fig: screenshot: R1 analysis} and~\ref{fig: screenshot: R2 analysis}).
653 These auto-analyses can be selected using the analysis selection wizard (Figure~\ref{fig: screenshot: analysis wizard} on page~\pageref{fig: screenshot: analysis wizard}).
654 Just as with the steady-state NOE in the next chapter, these auto-analyses are very similar in spirit to the sample script described in this chapter, though the Grace 2D visualisation is more advanced.
655 If you have read this chapter, the usage of these analyses should be self explanatory.
656
657 As in the script/prompt UI section above, the example of protein $^{15}$N $\Rone$ relaxation analysis will be performed in the following sections.
658 To keep track of all the messages relax produces for future reference, you can run the relax GUI with the following command line arguments:
659
660 \example{\$relax --log log --gui} 661 662 The messages will then appear both in the relax controller window (see Figure~\ref{fig: screenshot: relax controller} on page~\pageref{fig: screenshot: relax controller}) and in the \file{log} file. 663 664 665 % Initialisation of the data pipe. 666 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 667 668 \subsection{Relax-fit GUI mode -- initialisation of the data pipe} 669 670 To begin the analysis, launch the analysis selection wizard (see Figure~\ref{fig: screenshot: analysis wizard} on page \pageref{fig: screenshot: analysis wizard}). 671 Select either the$\Rone$or$\Rtwo$analyses, and change the name of the analysis if you plan on running multiple analyses from different field strengths in one relax instance. 672 673 \begin{minipage}[h]{\linewidth} 674 \centerline{ 675 \includegraphics[ 676 width=0.8\textwidth, 677 bb=14 14 1415 1019 678 ] 679 {graphics/screenshots/r1_analysis/analysis_wizard1} 680 } 681 \end{minipage} 682 683 Then click on the \guibutton{Next} button. 684 On the second page click on \guibutton{Start} to commence the analysis -- this second part of the wizard does not need to be changed. 685 For the$\Rone$and$\Rtwo$analyses in the GUI, a data pipe bundle containing only a single data pipe for that analysis will be created. 686 This data pipe bundle can be safely ignored. 687 688 \begin{minipage}[h]{\linewidth} 689 \centerline{ 690 \includegraphics[ 691 width=0.8\textwidth, 692 bb=14 14 1415 1019 693 ] 694 {graphics/screenshots/r1_analysis/analysis_wizard2} 695 } 696 \end{minipage} 697 698 699 % General setup. 700 %~~~~~~~~~~~~~~~ 701 702 \subsection{Relax-fit GUI mode -- general setup} 703 704 You will now be presented with a blank analysis tab: 705 706 \begin{minipage}[h]{\linewidth} 707 \centerline{ 708 \includegraphics[ 709 width=0.8\textwidth, 710 bb=14 14 1415 1019 711 ] 712 {graphics/screenshots/r1_analysis/blank} 713 } 714 \end{minipage} 715 716 Here there are two things unique to the GUI which need to be preformed: 717 \begin{description} 718 \item[NMR frequency label:] First set the NMR frequency label. 719 This is only used for the name of the output file. 720 For example if you set the label to \guistring{1200}, the file \file{r1.1200.out} will be created at the end of the analysis. 721 \item[Results directory:] All of the automatically created results and Grace files will be placed into this directory. 722 The \gui{Results directory} can now be changed. 723 \end{description} 724 725 726 % Spin systems. 727 %~~~~~~~~~~~~~~ 728 729 \subsection{Relax-fit GUI mode -- setting up the spin systems} 730 731 As the relaxation data is at the level of the spins, the molecule, residue and spin data structures need to be set up. 732 In the$\Rone$and$\Rtwo$GUI analysis tabs, there is a special \gui{Spin systems} GUI element designed for this. 733 This will initially say \gui{0 spins loaded and selected}. 734 Click on the \guibutton{Spin editor} button to launch the spin viewer window. 735 The steps for setting up the spin containers using PDB files are described in section~\ref{sect: GUI - structural data} on page~\pageref{sect: GUI - structural data} or for sequence files in section~\ref{sect: GUI - sequence file} on page~\pageref{sect: GUI - sequence file}. 736 737 738 % Unresolved spins. 739 %~~~~~~~~~~~~~~~~~~ 740 741 \subsection{Relax-fit GUI mode -- unresolved spins} 742 743 As in the prompt/script UI section~\ref{sect: Rx setup fin}, the spins can be deselected at this point using the same \file{unresolved} file. 744 This is described in detail in section~\ref{sect: GUI - deselect spins} on page~\pageref{sect: GUI - deselect spins}. 745 746 747 748 % Loading the data. 749 %~~~~~~~~~~~~~~~~~~ 750 751 \subsection{Relax-fit GUI mode -- loading the data} 752 753 At this stage, the peak intensity data needs to be loaded. 754 In both the$\Rone$and$\Rtwo$analysis tabs is a \gui{Spectra list} GUI element. 755 Click on the \guibutton{Add} button to launch the peak intensity loading wizard: 756 757 \begin{minipage}[h]{\linewidth} 758 \centerline{ 759 \includegraphics[ 760 width=0.8\textwidth, 761 bb=14 14 1415 1019 762 ] 763 {graphics/screenshots/r1_analysis/peak_intensity_bb_peaks} 764 } 765 \end{minipage} 766 767 In this example, a Sparky peak list containing the peak heights determined from the averaged chemical shift positions for all spectra will be loaded. 768 Set the spectrum ID string to a unique value. 769 Click on \guibutton{Next}. 770 This will most likely cause a \prompt{RelaxWarning} message to appear for all peak list elements which do not correspond to any spins loaded into the relax data store: 771 772 \begin{minipage}[h]{\linewidth} 773 \centerline{ 774 \includegraphics[ 775 width=0.8\textwidth, 776 bb=14 14 1415 1019 777 ] 778 {graphics/screenshots/r1_analysis/peak_intensity_warnings} 779 } 780 \end{minipage} 781 782 These messages must be carefully checked to be sure that the correct data has been loaded. 783 A \prompt{RelaxError} might be thrown if the peak list is corrupted or if the dimension has been incorrectly given. 784 In this case check the message, go \guibutton{Back}, fix the problem, and click on \guibutton{Next} again. 785 Then click on \guibutton{Next}. 786 You should now see the error type page: 787 788 \begin{minipage}[h]{\linewidth} 789 \centerline{ 790 \includegraphics[ 791 width=0.8\textwidth, 792 bb=14 14 1415 1019 793 ] 794 {graphics/screenshots/r1_analysis/peak_intensity_err_type} 795 } 796 \end{minipage} 797 798 The description for this wizard page should be very carefully read -- it will tell you about all of the error analysis options available and how these are implemented in relax. 799 For the protein relaxation example, replicated spectra have been collected. 800 Therefore the option \gui{Replicated spectra} will be chosen. 801 The \gui{Baseplane RMSD} option is documented in the NOE chapter. 802 After clicking on \guibutton{Next} you will see: 803 804 \begin{minipage}[h]{\linewidth} 805 \centerline{ 806 \includegraphics[ 807 width=0.8\textwidth, 808 bb=14 14 1415 1019 809 ] 810 {graphics/screenshots/r1_analysis/peak_intensity_replicates1} 811 } 812 \end{minipage} 813 814 For the first of the duplicate spectra, or any spectrum without a duplicate, you can click on the \guibutton{Skip} button. 815 If this is the second spectrum you have loaded from a duplicated set, select the two replicated spectra and then click on \guibutton{Next}: 816 817 \begin{minipage}[h]{\linewidth} 818 \centerline{ 819 \includegraphics[ 820 width=0.8\textwidth, 821 bb=14 14 1415 1019 822 ] 823 {graphics/screenshots/r1_analysis/peak_intensity_replicates2} 824 } 825 \end{minipage} 826 827 Finally set the relaxation time period for this experiment in seconds: 828 829 \begin{minipage}[h]{\linewidth} 830 \centerline{ 831 \includegraphics[ 832 width=0.8\textwidth, 833 bb=14 14 1415 1019 834 ] 835 {graphics/screenshots/r1_analysis/peak_intensity_times} 836 } 837 \end{minipage} 838 839 All delays and pulse lengths in the pulse sequence should be carefully checked to be sure that the time is exactly what you would expect -- the estimated time may not match the real time. 840 To set the time and close the wizard, click on the \guibutton{Finish} button. 841 842 This procedure should be repeated for every experiment you have collected (you could, as an alternative, load all at the same time using the \guibutton{Apply} button at each stage). 843 In the end you should see something such as: 844 845 \begin{minipage}[h]{\linewidth} 846 \centerline{ 847 \includegraphics[ 848 width=0.8\textwidth, 849 bb=14 14 1415 1019 850 ] 851 {graphics/screenshots/r1_analysis/analysis_tab_full} 852 } 853 \end{minipage} 854 855 856 % Optimisation and error analysis. 857 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 858 859 \subsection{Relax-fit GUI mode -- optimisation and error analysis} 860 861 Back in the main$\Rone\$ analysis tab, the grid search increments and number of Monte Carlo simulations can be changed.
862 The default values of 21 grid search increments and 500 MC simulations are optimal -- lower values are not recommended.
863 To perform the optimisation and error analysis, click on the \guibutton{Execute relax} button.
864 The relax controller will open to show you the progress of the optimisation and simulations:
865
866 \begin{minipage}[h]{\linewidth}
867   \centerline{
868     \includegraphics[
869       width=0.8\textwidth,
870       bb=14 14 1415 1019
871     ]
872     {graphics/screenshots/r1_analysis/mc_sim}
873   }
874 \end{minipage}
875
876 Once finished, the \gui{Results viewer} window will also appear:
877
878 \begin{minipage}[h]{\linewidth}
879   \centerline{
880     \includegraphics[
881       width=0.8\textwidth,
882       bb=14 14 1415 1019
883     ]
884     {graphics/screenshots/r1_analysis/fin}
885   }
886 \end{minipage}
887
888 This window can be used to open the text files in the default text editor for your operating system or the 2D Grace plots in \prompt{xmgrace} if available on your system.
889
890
891
892 % Final checks.
893 %%%%%%%%%%%%%%%
894
895 \section{Final checks of the curve-fitting}
896
897 To be sure that the data has been properly collected and that no instrumentation or pulse sequence timing errors have occurred, it is essential to carefully check the \file{intensities.agr} and \file{intensities\osus{}norm.agr} 2D Grace\index{software!Grace|textbf} files.
898 These are plots of the decay curves for each spin system analysed, and any non-exponential behaviour should be clearly visible (see Figure~\ref{fig: screenshot: xmgrace peak intensities}).
899 If Xmgrace or a compatible program is not available for your operating system, the Grace files contain text representations of the curves at the end which can opened, edited and visualised in any another 2D graphing software package.
900
901 % Xmgrace screenshot
902 \begin{figure}
903   \centerline{
904     \includegraphics[
905       width=\textwidth,
906       bb=14 14 923 724
907     ]
908     {graphics/screenshots/xmgrace_peak_intensities}
909   }
910   \caption[Peak intensity 2D plot xmgrace screenshot]{
911     Screenshot of the 2D peak intensity plots for the exponential relaxation curves in Xmgrace.
912   }
913   \label{fig: screenshot: xmgrace peak intensities}
914 \end{figure}
915
916 Note that errors resulting in systematic bias in the data -- for example if temperature control (single-scan interleaving or temperature compensation blocks) or per-experiment/per-spectrometer temperature calibration on MeOH or ethylene glycol have not been performed -- will not be detected by looking at the decay curves.
917 See section~\ref{sect: temperature control and calibration} or the \uf{relax\ufus{}data\ufsep{}temp\ufus{}calibration} user function documentation on page~\pageref{uf: relax_data.temp_calibration} and the \uf{relax\ufus{}data\ufsep{}temp\ufus{}control} user function documentation on page~\pageref{uf: relax_data.temp_control} for more details.