Hi Edward,
ok I uploaded a tar file with synthetic peaks for one spin system:
file #22989
The saturation recovery formula is:
I0*(1 - exp(−R1*t))
The parameters I used were
I0 = 1000000000000000.00
and
R1 = 0.5
Obviously, the same numbers also work for the inversion recovery experiment:
I(t) = I∞ − I0*exp(−R1*t)
with I∞ set to 1000000000000000.00 as well. Not sure how much that helps
though.
If "inv" is not implemented, then that would explain the weird results :)
Unfortunately, I have no primary reference for the saturation recovery
experiment either.
Best regards,
Andras
________________________________________
From: edward.dauvergne@xxxxxxxxx [edward.dauvergne@xxxxxxxxx] On Behalf Of
Edward d'Auvergne [edward@xxxxxxxxxxxxx]
Sent: Thursday, November 27, 2014 11:32 AM
To: Boeszoermenyi, Andras
Cc: relax-devel@xxxxxxx
Subject: Re: [sr #3195] Fitting of saturation recovery experiment
Hi Andras,
Thank you for your help.
No problems. Hopefully we can get this running without too much effort.
Are Monte Carlo simulations used for the error analysis in CCPNMR
(http://www.ccpn.ac.uk/software/analysis)? What do you mean that it
is not robust enough? From the equation you sent, it does not look
like your experiment is either inversion recovery or saturation
recovery.
I don't know how CCPNMR does it's fitting (it's a bit of black box), but I
don't believe that it does Monte Carlo simulations. What I mean by not
robust enough is, that the fit fails instead of producing a result. The
fitting routine in ccpnmr takes seconds or less, and therefore I guess that
it is not very sophisticated.
I thought it did carry out real Monte Carlo simulations, unlike a lot
of the other software out there (excluding curvefit and relax), but I
no longer remember. I think I asked this once on the relax mailing
lists and received a response, but I cannot find the old posts.
Anyway, the fitting is also very quick in relax as this is an
incredibly simple problem. The optimisation part for a single spin
can be measured in milliseconds.
Do you have a good original reference for your experiment? I.e. one
with the pulse sequence and also the exponential equation used? You
said that the formula for your data is:
I(t) = I0*(1 - exp(−R1*t)) (1)
In the relax implementation, the formula is:
I(t) = I0*exp(-Rx*t) (2)
These are clearly not the same as (2) heads to zero as t->inf, whereas
(1) heads to I0 as t->inf.
I am not comparing the saturation recovery experiment to a standard
relaxation decay, but to an inversion recovery experiment. The inversion
recovery experiment is implemented in relax "according to the manual" as
I'm not sure what a standard R1 decay is. I would have guessed a 180
degree pulse on the spin in the ground state and then watch the peak
intensity return to the ground state - i.e. standard inversion
recovery.
I(t) = I∞ − I0*exp(−R1*t), which is a lot more similar to the formula for
saturation recovery.
This is currently not implemented in relax. This is where the
inversion-recovery branch stalled. So selecting the model 'inv' has
no effect. This will require someone with motivation to help complete
this ;) The hardest part is the synthetic data and creation of the
associated system test.
The reference that I am using for saturation recovery is:
Respondek et al., 2007, Mapping the Orientation of Helices in Micelle-Bound
Peptides by Paramagnetic Relaxation Waves. JACS.
For permanent reference, the DOI link is
http://dx.doi.org/10.1021/ja069004f . Equation (1) in that reference
is:
I(t) = I0*(1 - exp(−R1*t)) (1)
So this implies that you start at time 0 with zero intensity and then
watch the magnitization return to I0. Is this correct? Is this what
'saturation recovery' means? From Figure 3, and your grace plots
attached to https://gna.org/task/?7415, I would guess so. It's a pity
this paper doesn't have the reference for this experiment type. Would
you happen to know the primary reference for 'saturation recovery' or
'signal pre-saturation recovery'?
Anyway, with test data this formula should be very easy to add support for.
If you look at the plots in the grace folder of the project that I
uploaded, you can see that everything looks exactly as it should, except
that the Rx values are negative.
This is very strange as the only formula currently implemented in relax is:
I(t) = I0*exp(-Rx*t) (2)
In the saved state file, I can see that the chi-squared values are
huge though. This implies that the fit is terrible! So I would not
trust these numbers.
I will need to look into the generation of synthetic peak lists some more,
as I have no experience with that. But once we have those, I can load
synthetic peaks (in e.g. sparky format) into ccpnmr and fit them there to
double check results.
It can be done with a calculator, spreadsheet, or script. Just
calculate new numbers and replace the old using a text editor. For
example take your 800mhz-303k-*.list files, delete all lines except
for the header and first data line (G17), and then replace its numbers
with your hand calculated ones. That's it. You could also run it
through a relax script to see what happens. If you upload the new
files to https://gna.org/task/?7415, this relax script, and say what
the I0 and Rx numbers are in the comment section, all in one post,
then I could create a system test with this in ~5 min.
Cheers,
Edward
P. S. Looking at the grace/intensities_norm.agr file attached to
https://gna.org/task/?7415, I think we will need to rethink how
normalisation occurs for this experiment type! Dividing all points by
I0 might be better. Then all curves terminate at the same intensity
of 1, and you can better differentiate the rates, as setting the first
point to 1 in these graphs to normalise is a bit silly.
Best wishes,
Andras
________________________________________
From: edward.dauvergne@xxxxxxxxx [edward.dauvergne@xxxxxxxxx] On Behalf Of
Edward d'Auvergne [edward@xxxxxxxxxxxxx]
Sent: Thursday, November 27, 2014 4:07 AM
To: Boeszoermenyi, Andras
Cc: relax-devel@xxxxxxx
Subject: Re: [sr #3195] Fitting of saturation recovery experiment
Hi Andras,
Please see below:
I would be very glad to help out with that.
Great!
My input data format is sparky peak lists (I export them from ccpnmr, but
that doesn't make a difference))
relax should handle this no problem. There are system tests set up to
check this compatibility:
https://gna.org/support/?3044
http://www.nmr-relax.com/api/3.3/test_suite.system_tests.peak_lists.Peak_lists-class.html#test_ccpn_analysis
http://www.nmr-relax.com/api/3.3/test_suite.system_tests.peak_lists-pysrc.html#Peak_lists.test_ccpn_analysis
The other part of the input data is a sequence file. I have problems
loading a cyana generated structure. Not sure how it would behave with a
structure outputted from ccpnmr or other programs. The structure is not in
the pdb yet.
This should be no problem in relax. Any column formatted file will be
sufficient. You can even use the spectrum.read_spins user function to
load the sequence out of a peak list. Anyway, it is better to discuss
such topics on the relax-users mailing list. I would recommend
signing up for the relax-users and relax-devel mailing list using the
links:
https://mail.gna.org/listinfo/relax-users/
https://mail.gna.org/listinfo/relax-devel/
Using these same links, you can very quickly terminate your
subscription later on. The advantage of subscribing is that your
messages will be instantly accepted. Otherwise your message will have
to wait for the mailing list maintainer, Chris MacRaild
(http://www.monash.edu.au/pharm/research/researchers/profile.html?sid=47804&pid=5391).
As he is based in Melbourne, Australia, there will be a time zone
delay between messages sent, and messages accepted.
I have two spin systems in there, but for testing, we can throw out the
NHe from my tryptophan indole rings.
For testing it is best to use synthetic data, see below.
For my test set, I can even calculate rates for inversion recovery and
saturation recovery in ccpnmr. The reason that I changed to relax is, that
the ccpnmr fitting routine is not robust enough for other sets.
Are Monte Carlo simulations used for the error analysis in CCPNMR
(http://www.ccpn.ac.uk/software/analysis)? What do you mean that it
is not robust enough? From the equation you sent, it does not look
like your experiment is either inversion recovery or saturation
recovery.
I'm very happy to send you real peak lists, or do they have to be
synthetic for some reason? I could make real peak lists without tryptophan
NHes if that helps. For synthetic peak lists, I have no experience in
making those. But I think I could figure that out.
For testing an implementation, real data is to be avoided. The reason
is because you do not know beforehand what the real answer is. In
software development, this can often lead to an incorrect
implementation.
A synthetic set is easy. Take a set of Sparky formatted peak lists
for a single spin. Then pick Rx and I0 values that are reasonable.
For R1, a value of 2.0 might be a good choice. You could have a
second spin with a different R1 value if you wish. For the I0, pick a
huge value in the order of 1e15. The reason is because the peak
heights in a Sparky peak list are stored as integer values and not
floating point values. By starting with a large intensity, you
thereby minimise truncation artefacts. For the R1 and I0 value,
together with the time points used in the experiment, the formula of
interest can be used to then calculate the intensity for each time
point. These can then replace the values in the Sparky peak lists to
create the synthetic test set. Do not truncate the numbers by using
the format 9.62E+07, instead include the full integer number into the
text file.
Using this, a system test can then be set up in relax. It would
consist of a script similar to your wr1043_0N_relax_fit.py script.
Optimisation would be performed, and then the values of R1 and I0
checked to see if the exact values can be found (to a very high degree
of accuracy). For example as in the
Relax_fit.test_zooming_grid_search system test which can be run as:
$ relax -sd Relax_fit.test_zooming_grid_search system
This assumes you are not using the Mac OS X dmg distribution file.
The test itself is:
http://www.nmr-relax.com/api/3.3/test_suite.system_tests.relax_fit.Relax_fit-class.html#test_zooming_grid_search
http://www.nmr-relax.com/api/3.3/test_suite.system_tests.relax_fit-pysrc.html#Relax_fit.test_zooming_grid_search
This uses the script located at
test_suite/system_tests/scripts/relax_fit_zooming_grid.py.
I am attaching one set of relax that I did on this system with inversion
recovery to the task and if this is what you are looking for, then I can
also send you R1 values for inversion and saturation recovery fitting of
this set in ccpnmr.
For reference, the file is attached to task #7415
(https://gna.org/task/?7415) as the file Relax.tar.gz
(https://gna.org/task/download.php?file_id=22980). The synthetic data
would be far better. Could such synthetic data be put through CCPNMR?
Once the synthetic cases are implemented and returning the expected
result, then it could be applied to real data and comparisons between
software made. However comparing the results of the synthetic data
for different software is a far better method, as it is a real
positive control.
I did have a problem with the experiment that I ran on relax though. It
does converge and produces beautiful values, but the Rx values are
negative. Which doesn't make any sense to me. The pattern is exactly as it
should be, but the values are simply wrong. Do you have an idea what could
be going on there?
Do you have a good original reference for your experiment? I.e. one
with the pulse sequence and also the exponential equation used? You
said that the formula for your data is:
I(t) = I0*(1 - exp(−R1*t)) (1)
In the relax implementation, the formula is:
I(t) = I0*exp(-Rx*t) (2)
These are clearly not the same as (2) heads to zero as t->inf, whereas
(1) heads to I0 as t->inf. To determine exactly what is happening,
the synthetic data would be best. Are you sure you reported the
correct equation? Equation (1) is zero when t is zero, which does not
seem to make sense. Should the saturation recover not start at I0 and
return to zero? Equation (2) is supposed to be the one used for
saturation recovery. Maybe your data is actually in the form of:
I(t) = -I0*exp(-Rx*t) (3)
Cheers,
Edward
P. S. Note that there is a 6 hour time zone difference, so I may only
be able to answer you in your mornings.
Best regards,
Andras
________________________________________
From: edward.dauvergne@xxxxxxxxx [edward.dauvergne@xxxxxxxxx] On Behalf Of
Edward d'Auvergne [edward@xxxxxxxxxxxxx]
Sent: Wednesday, November 26, 2014 1:04 PM
To: Andras Boeszoermenyi
Cc: Boeszoermenyi, Andras; relax-devel@xxxxxxx
Subject: Re: [sr #3195] Fitting of saturation recovery experiment
Hi Andras,
I was wondering if you were interested in helping out with the
inversion-recovery task at https://gna.org/task/index.php?7415 ? What
is your input data format, i.e. which peak list format are you using?
Would it be possible to create synthetic peak lists containing with a
single spin system whereby you have calculated what the peak
intensities are based on your equation and some I0 and Rx parameter
values? If this could be done, it could be added to the relax test
suite directories and form the basis of an implementation. It could
easily then be used to create a system test which checks that the
fitting occurs as expected. Then implementing this should be easy.
Cheers,
Edward
On 26 November 2014 at 16:59, Andras Boeszoermenyi
<NO-REPLY.INVALID-ADDRESS@xxxxxxx> wrote:
URL:
<http://gna.org/support/?3195>
Summary: Fitting of saturation recovery experiment
Project: relax
Submitted by: andras
Submitted on: Wed 26 Nov 2014 03:59:50 PM UTC
Category: Feature request
Priority: 5 - Normal
Severity: 3 - Normal
Status: None
Assigned to: None
Originator Email:
Open/Closed: Open
Discussion Lock: Any
Operating System: Mac OS
_______________________________________________________
Details:
Hi,
would it be possible to add an option to fit saturation recovery
experiments
to the existing 'exp' and 'inv' options in relax_fit?
The formula should look like this:
I(t) = I0*(1 - exp(−R1*t))
Best regards,
Andras Boeszoermenyi
_______________________________________________________
Reply to this item at:
<http://gna.org/support/?3195>
_______________________________________________
Message sent via/by Gna!
http://gna.org/
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-devel mailing list
relax-devel@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-devel
RE: [sr #3195] Fitting of saturation recovery experiment