mailRe: Help with settings script for R1rho analysis.


Others Months | Index by Date | Thread Index
>>   [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Header


Content

Posted by Troels Emtekær Linnet on October 22, 2013 - 19:17:
Hi Edward.

I prefer to do the calculations straight away in my scripts, from the information which is ex-tractable from the procpar files.
But I guess that is a matter of taste.

But now I think I got it. :-) Thanks!
I guess that relax calculates "omega_rf_ppm - chemShift{peakName} ", but I could not locate this functions call.
( I do load the chemical shifts from a SPARKY list, (the seriesTab format is not yet supported). )


I now do the settings script, by settting variables:
# In MHz
yOBS = 81.050
# In ppm
yCAR = 118.078
centerPPM_N15 = yCAR

And then I calculate the offset in ppm for each spectrum.

# Calculating the spin-lock offset in ppm, from offsets values provided in Hz.
#frq_N15_Hz = set_sfrq * 1E6 * gyro15N / gyro1H
frq_N15_Hz = yOBS * 1E6
offset_ppm_N15 = float(deltadof2) / frq_N15_Hz * 1E6
omega_rf_ppm = centerPPM_N15 + offset_ppm_N15

And the range is now 118.078 ppm to 241.45 ppm.

I tried to locate the corresponding calculation of OMEGA in relax, but i was not successful.
I looked in:

lib/dispersion/dpl94.py
target_functions/relax_disp.py

I wonder where how to locate the calculation of theta?

Best
Troels






2013/10/22 Edward d'Auvergne <edward@xxxxxxxxxxxxx>
Hi,

You could look at the relax code for how omega_eff is calculated.
However you will never use this - it is never input into relax.  If
you have a look at the sample_scripts/relax_disp/R1rho_analysis.py
script, you will see that all is needed is an equivalent of the table:

# The spectral data - spectrum ID, peak list file name, spin-lock
field strength (Hz), the spin-lock offset (ppm), the relaxation time
(s), spectrometer frequency (Hz), and experimental error (RMSD of the
base plane noise for each spectrum).
data = "">     ['ref_500MHz',       'ref_500MHz.list',     ,   None, 110.0, 0.1,
500e6, 200000.0]
    ['nu_1000.0_500MHz', 'nu_1000.0_500MHz.list', 1000.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_1500.0_500MHz', 'nu_1500.0_500MHz.list', 1500.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_2000.0_500MHz', 'nu_2000.0_500MHz.list', 2000.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_2500.0_500MHz', 'nu_2500.0_500MHz.list', 2500.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_3000.0_500MHz', 'nu_3000.0_500MHz.list', 3000.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_3500.0_500MHz', 'nu_3500.0_500MHz.list', 3500.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_4000.0_500MHz', 'nu_4000.0_500MHz.list', 4000.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_4500.0_500MHz', 'nu_4500.0_500MHz.list', 4500.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_5000.0_500MHz', 'nu_5000.0_500MHz.list', 5000.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_5500.0_500MHz', 'nu_5500.0_500MHz.list', 5500.0, 110.0, 0.1,
500e6, 200000.0]
    ['nu_6000.0_500MHz', 'nu_6000.0_500MHz.list', 6000.0, 110.0, 0.1,
500e6, 200000.0]
    ['ref_800MHz',       'ref_800MHz.list',     ,   None, 110.0, 0.1,
800e6, 200000.0]
    ['nu_1000.0_800MHz', 'nu_1000.0_800MHz.list', 1000.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_1500.0_800MHz', 'nu_1500.0_800MHz.list', 1500.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_2000.0_800MHz', 'nu_2000.0_800MHz.list', 2000.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_2500.0_800MHz', 'nu_2500.0_800MHz.list', 2500.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_3000.0_800MHz', 'nu_3000.0_800MHz.list', 3000.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_3500.0_800MHz', 'nu_3500.0_800MHz.list', 3500.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_4000.0_800MHz', 'nu_4000.0_800MHz.list', 4000.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_4500.0_800MHz', 'nu_4500.0_800MHz.list', 4500.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_5000.0_800MHz', 'nu_5000.0_800MHz.list', 5000.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_5500.0_800MHz', 'nu_5500.0_800MHz.list', 5500.0, 110.0, 0.1,
800e6, 200000.0]
    ['nu_6000.0_800MHz', 'nu_6000.0_800MHz.list', 6000.0, 110.0, 0.1,
800e6, 200000.0]
]

These values can be hardcoded into a script.  There is no need to do
this programatically - they will not change.  Though you could
programmatically generate such a table in a separate script, if you
wish.  You almost have this table anyway with the
exp_parameters_sort.txt file.  It is good, for sanity's sake, to have
such a complete summary table in the standard units.  As for
omega_eff, it uses the data in this table together with the chemical
shifts loaded via the chemical_shift.read user function to
automatically determine the values.  It would be best read these
shifts from one of your *.ser files.

I would recommend to use the
sample_scripts/relax_disp/R1rho_analysis.py script as a template for
your whole analysis.  For the test suite, I would recommend copying
test_suite/system_tests/scripts/relax_disp/r1rho_off_res_tp02.py and
making a few minor modifications.  For example using a table as in the
sample script.  This r1rho_off_res_tp02.py script does not have a
table as it deals with perfect synthetic data.

Regards,

Edward

On 22 October 2013 14:54, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
> Hi Edward.
>
> I am digging through old perl scripts, which is preparing data for our
> analysis in IgorPro.
>
> So I guess my challenge is now to match this to the procedures of relax, and
> establish the difference.
>
> -----------
> $expList="expList.txt";
> $peakFile='peaks.dat';
> $centerPPM=118.085;
> $frq=81.050;
> -----------
>
> expList.txt contains lines with: "fit_R1_filename"  "spin_lock_offset_HZ"
> "spin_lock_field_Hz"
>
> ----------------------
> open IN, "$expList" or die "Cannot open $expList for read";
> while (<IN>){
> @process = split (/\s+/, $_);
> if ($process[0] ne "#"){
> $fileName{$i}=$process[0];
> $offset{$i}=$process[1];
> $omega1{$i}=$process[2];
> $i++;
> };
> };
> close (IN);
> -------------------
>
> And then it read the peak file to extract ppm for N15.
>
> -------------
> open IN , "$peakFile" or die "Cannot open $peakFile for read";
> while (<IN>) {
> @process = split (/\s+/, $_);
> if ($process[0] eq ""){splice (@process, 0, 1)};
> if ($process[0] == 1){ $read = 1};
> if ($read == 1){
> $chemShift{$process[6]}=$process[4];
> };
> };
> close (IN);
> -----------
> The chemical shifts are stored in a dictionary under the residue name
> ($process[6]) with the ppm values of N15 ($process[4]).
>
> Then R1r and R1r_err are read in from external exponential fit files.
> ------------
> for ($i = 1; $i <= $numFiles; $i++){
> open IN, "$fileName{$i}" or die "Cannot open $fileName{$i} for read";
> while (<IN>){
> @process = split (/\s+/, $_);
> if ($process[0] ne "#"){
> $R1r{$process[0]}{$i}=$process[1];
> $R1r_err{$process[0]}{$i}=$process[2];
> };
> };
> close(IN);
> };
> -------------------
>
>
> Then omega is calculated.
> ------------
> foreach $peakName (keys %chemShift) {
> open (OUT, ">../residueFiles/$peakName.dat");
> for ($i=1;$i<=$numFiles;$i++){
> $OMEGA=($centerPPM-$chemShift{$peakName})*$frq+$offset{$i};
> $omegaEFF=sqrt($OMEGA**2+$omega1{$i}**2);
> if (($omega1{$i}/$OMEGA) > 0){
> $theta=180/$PI*abs(atan($omega1{$i}/$OMEGA));
> }else{
> $theta=180-180/$PI*abs(atan($omega1{$i}/$OMEGA));
> };
> printf OUT "%s %s %s %s %s
> %s\n",$OMEGA,$omega1{$i},$omegaEFF,$theta,$R1r{$peakName}{$i},$R1r_err{$peakName}{$i};
> };
> close (OUT);
> };
> --------------
>
> I guess I would need to find out how relax calculates omegaEFF, the
> effective field in the rotating frame ?
>
> Best
> Troels
>
>
>
>
> 2013/10/22 Edward d'Auvergne <edward@xxxxxxxxxxxxx>
>>
>> Hi,
>>
>> Ok, I made the assumption that deltadof2 was not the spin-lock offset
>> (in ppm) but rather the spin-lock field strength (in Hz).  The ppm
>> values as printed out from the r1rho_1_ini.py script are far too low
>> for 15N.  They should be in the range of 100-120 ppm!  I.e. around the
>> 15N chemical shifts.  They can be much higher and much lower for
>> off-resonance data, but having many at an offset of 0.0 ppm seems a
>> little strange.
>>
>> Regards,
>>
>> Edward
>>
>> On 22 October 2013 13:18, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
>> wrote:
>> > Hi Edward.
>> >
>> > Thanks for looking at this.
>> >
>> > The offset for N15 spinlock is given in Hz as parameter deltadof2
>> > for our pulse sequence.
>> >
>> > So I guess that this conversions should be appropriate:
>> > omega_rf_ppm = float(deltadof2) / (set_sfrq * 1E6) * 1E6
>> >
>> > where set_sfrq=799.7773991 is the spectrometer frequency.
>> >
>> > Best
>> > Troels
>> >
>> >
>> > 2013/10/22 Edward d'Auvergne <edward@xxxxxxxxxxxxx>
>> >>
>> >> Hi,
>> >>
>> >> In the relax prompt, have a look at:
>> >>
>> >> relax> help(relax_disp.spin_lock_offset)
>> >>
>> >> The unit you need is ppm.  You can find this number directly from the
>> >> Bruker acqus or Varian procpar files.  It is the ppm offset of the
>> >> spin-lock pulse which is manually set by the person recording for each
>> >> spectrum.  You may need to go to the original pulse sequence to know
>> >> which pulse this is.  It will be a fixed value for each spectrum
>> >> collected and you can use that value directly.  The calculations in
>> >> your scripts seem far to complicated and it looks like the internal
>> >> conversions performed within relax.  You should never need
>> >> gyromagnetic ratios, factors of pi, etc.  If you do, you should
>> >> probably go back to the original experiments and pull out the correct
>> >> offset number in ppm for each spectrum.  I hope this helps.
>> >>
>> >> Regards,
>> >>
>> >> Edward
>> >>
>> >>
>> >>
>> >>
>> >> On 22 October 2013 10:46, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
>> >> wrote:
>> >> > Hi Edward.
>> >> >
>> >> > I wonder if would have time to look at the settings script in:
>> >> > test_suite/shared_data/dispersion/Kjaergaard_et_al_2013
>> >> >
>> >> > It is the scripts:
>> >> > r1rho_1_ini.py
>> >> > r1rho_3_spectra_settings.py
>> >> >
>> >> > In r1rho, I am unsure if I do the correct conversion from Hz to ppm
>> >> > for
>> >> > omega_rf.
>> >> >
>> >> > Again, if there is hidden radian units?
>> >> >
>> >> > It is when I set:
>> >> > relax_disp.spin_lock_offset
>> >> >
>> >> > And I wonder, how to start modifying, so R1 rates can be read per
>> >> > spectra?
>> >> >
>> >> > Best
>> >> > Troels
>> >> >
>> >> >
>> >> >
>> >> >
>> >> >
>> >> > _______________________________________________
>> >> > 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
>> >> >
>> >
>> >
>
>


Related Messages


Powered by MHonArc, Updated Tue Oct 22 19:40:07 2013