Re: fixed tm values
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
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