Re: fixed tm values

2015-06-08 Thread Edward d'Auvergne
Hi Christina,

Welcome to the relax mailing lists!  As Troels mentioned, I have been
on holidays so can only now reply.  For details, please see below:

On 3 June 2015 at 14:32, Christina Möller  wrote:
> Dear Edward and relax users,
>
> I successfully performed the dauvergne_protocol.py analysis scripts to
> determine the ps - ns dynamics of a 14 kDa protein. The global
> correlation time tm is 9.6 ns.

This value does not seem unreasonable.

> Since other methods suggested a smaller
> global correlation time tm,

Which other methods have you used?  You should note the following issues:

- It is well known that MD underestimates the global correlation time
by about half and David Case is actively researching this area.
- Stoke's law and the related hydrodynamic beads model from Garcia de
la Tor are reasonable for estimating tm for isolated molecules.
However with the super high concentration in NMR samples, these tend
to underestimate the tm value by about half as they do not take
micro-viscosity into consideration.
- ORD measurements as well, as the tm changes between the lower
concentration optical spectroscopic sample and the NMR sample by about
a factor of 2 (micro-viscosity effects being a major factor again).


> I would like to know whether it is possible
> to fix the tm value in all rounds of individual model-free
> optimisations?

As Troels pointed out, this is of course possible.  You can take, for
example, the sample_scripts/model_free/diff_min.py script.  This
already implements a lot of what Troels described for you.  Note
however that this needs to be executed iteratively until the global
chi-squared value between iterations is identical.  See figure 2 of:

d'Auvergne, E. J. and Gooley, P. R. (2008). Optimisation of NMR
dynamic models II. A new methodology for the dual optimisation of the
model-free parameters and the Brownian rotational diffusion tensor. J.
Biomol. NMR, 40(2), 121-133. (
http://dx.doi.org/10.1007/s10858-007-9213-3 ).

This global iterative optimisation of the diffusion tensor is
essential, and you'll find it documented in Mandel et al., 1995 as
well as in the papers by the Dasha authors.  It will take between 5-20
iterations to properly converge (in rare cases it can be 2 iterations,
in others >50).  For more details, you really should read:

How to get model free parameters from output files (
http://thread.gmane.org/gmane.science.nmr.relax.user/1375/focus=1378
).

This requires an initial diffusion tensor estimate.  And as I
demonstrated in the above paper - as well as first shown in the
Korzhnev et al., 1999 bacteriorhodopsin fragment paper (
http://dx.doi.org/10.1023/a:1008356809071 ) - if this is too far off,
the global minimum will never be reached.  You can however directly
compare the two results using AIC values (not chi-squared values as
the individual model-free models for each residue will be different
and the effects of parsimony will not be taken into account).  An
additional thread which might be of interest is:

AIC to select diffusion model (
http://thread.gmane.org/gmane.science.nmr.relax.user/885/focus=891 ).

There are plenty of other relax-users mailing list threads on the
subject which you can search for at:

http://dir.gmane.org/gmane.science.nmr.relax.user


> The corresponding chi2 values might then be useful to
> evaluate the global correlation times that I obtained by different methods.

Note that you should compare the chi-squared values from the same
program, just to be sure.  Also note that the spherical angle and
Euler angle notations in Modelfree4, Dasha, Tensor2 and relax are not
compatible.  The problem is that the definitions of these angles are
not documented (except in relax) so if you take a diffusion tensor
from one and input it into another, you will see the angles swing
around wildly with the global iterative diffusion tensor estimate
until it converges to the same tensor but with the different angles
(except when a local minimum is hit).  There are 2406 Euler angle
conventions and symmetries for diffusion tensors!

One last thing to note is that relax is extremely flexible in what it
can do.  Using specially designed scripts, relax can replicate the
results of Modelfree4, Dasha, Tensor2, or DYNAMICS.  One exception is
that relax uses a real optimisation constraint algorithm (the
augmented Lagrangian or method of multipliers,
https://en.wikipedia.org/wiki/Augmented_Lagrangian_method ,
https://gna.org/projects/minfx/ ) which the other model-free softwares
do not, hence there can be cases where relax does not find exactly the
same result as the other softwares.

I hope all this information helps.  You should also consider Troels'
suggestion of the dx.map user function to see the diffusion tensor
parameter space (or 3D subsets of it).  I used this in figure 6 of:

d'Auvergne, E. J. and Gooley, P. R. (2008). Optimisation of NMR
dynamic models I. Minimisation algorithms and their performance within
the model-free and B

Re: fixed tm values

2015-06-03 Thread Troels Emtekær Linnet
Dear Christina Möller.

I am sure Edward can give you an answer to this.

I know that he is on holiday until the 7th, so he will probably first
respond after this.

Until then, I can try to give you some temporary idea's you can work with.
Beforehand, I would like to say, that I no experience with the whole setup
of modelfree, but
have only worked with the dispersion module of relax.

What you want to play with, is the function "value.set" and "minimise" in
relax.

Try this.
# Start relax in terminal mode
> relax
# Then see help functions
help(minimise)

# Function to look into
help(minimise.grid_search)
help(minimise.execute)
help(minimise.calculate)

"minimise.grid_search" is the function to find initial values, before
minimise.execute
is performed.
Grid searching is "expensive", and particular if clustering of residues is
part of this.

I have tried to redo some published results. Often these results for
dispersion only gives the global fitted paramter
of exchange rate or population, but neglect to report the individual
chemical shift change (dw).
See http://wiki.nmr-relax.com/CR72.

What I can do is then:

value.set(val=2000.0, param='kex', error=False, force=True)
value.set(val=0.98, param='pA', error=False, force=True)
minimise.grid_search()

minimise.grid_search will detect that kex and pA has already been set, and
therefore then lock these, and
will then only search for the optimal "dw" and "r2" parameter.

After this, I will do:
minimise.execute()

But here, all parameters are "set free", when minimising, and I don't think
you can "lock" a parameter in this.
And the results are often "close" to the published values, and within
experimental error.

What you can also do, is just manually setting a parameter value and do:
minimise.calculate
minimise.calculate simply just takes the stored parameter values and
calculate the chi2 value.
So here, you could test a range of your own tm, and do minimise.calculate
to see the chi2 values.

This goes fine for maybe testing one or two ranges of values of a parameter.
You would probably like to make tm/chi2 graph, to map the change in chi2
value.

What is more powerful is to use the dx.map function.
http://wiki.nmr-relax.com/Install_dx
http://wiki.nmr-relax.com/Dx_map

This maps 3 parameters, to visualize the chi2 hypersurface.
See Edwards paper on this:
http://www.nmr-relax.com/refs.html

Model-free model elimination: A new step in the model-free dynamic analysis
of NMR relaxation data
http://link.springer.com/article/10.1007%2Fs10858-006-9007-z

This will give you a VERY clear idea, how the chi2 surface is mapped for
the tm values.

relax has a lot of user functions to help with analysis, and a lot more
"plumbing functions", to do exactly as you wish.

Here are some ideas for a start, but adjust to your example.
(This is for dispersion)

--
state.load("final_state_from_your_analysis.bz2")

# A plumping function
from pipe_control.mol_res_spin import spin_loop

# Just test some print things
for spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True,
return_id=True, skip_desel=True):
 print dir(spin), spin.params

# Make an array
kex_arr = [1000., 2000., 3000.,]

for kex in kex_arr:
value.set(val=kex, param='kex', error=False, force=True)
minimise.calculate()
print kex, cdp.mol[0].res[0].spin[0].chi

-

If you need specific help with a saved state for practicing, you can ask
for support.
http://gna.org/support/?group=relax

Here you can describe the problem, and add a saved state from your final
minimization.
Then I can write a little script to help you out.

Note, you should only add a state, where you have deleted all residues,
except 2-3.
This is make sure, not to compromise any confidentiality of your data.

Best
Troels













2015-06-03 14:32 GMT+02:00 Christina Möller :

> Dear Edward and relax users,
>
> I successfully performed the dauvergne_protocol.py analysis scripts to
> determine the ps - ns dynamics of a 14 kDa protein. The global
> correlation time tm is 9.6 ns. Since other methods suggested a smaller
> global correlation time tm, I would like to know whether it is possible
> to fix the tm value in all rounds of individual model-free
> optimisations? The corresponding chi2 values might then be useful to
> evaluate the global correlation times that I obtained by different methods.
>
> Best regards,
> Christina
>
> --
> 
> Christina Möller
> ICS-6 / Structural Biochemistry
> Forschungszentrum Jülich
> D-52425 Jülich
> GERMANY
>
> E-mail: c.moel...@fz-juelich.de
> Tel.: +49-2461-619387
> Fax +49-2461-619387
>
>
>
>
> 
>
> 
> Forschungszentrum Juelich GmbH
> 52425 Juelich
> Sitz der Gesellschaft: Juelich
> Eingetragen im Handelsregister des Amtsgerichts Dueren Nr