View source for Tutorial for Relaxation dispersion analysis r1rho fixed time recorded on varian as sequential spectra

= Intro =
This tutorial is not complete.<br>

This tutorial is based on the analysis of R1rho data, analysed in a master thesis.

The spectra is not recorded interleaved, but as a series of spectra with experimental changes.

= Preparation =
You want to make a working dir, with different folders

<source lang="text">
peak_lists
spectrometer_data
scripts
</source>

You can create the folders by
<source lang="text">
mkdir peak_lists spectrometer_data scripts
</source>

In the folder '''peak_lists''' should contain SPARKY list in [[SPARKY_list | SPARKY list format]]. <br>
In the folder '''scripts''' we put scripts which help us processing the files. <br>
In the folder '''spectrometer_data''' should be a directory containing directories with all experiments where each directory contain files: '''fid''' and '''procpar''' as the output from recording on Varian.<br>

= Establish file-overview =

== Make file with paths to fid files ==
We make a file list of filepaths to '''fid''' files.
<source lang="bash">
ls -v -d -1 */fid > fid_files.ls
cat fid_files.ls
</source>

== Do something in each directory ==
With the '''fid_files.ls''', we can do ''something'' in each directory. <br>
For example do a list files in each direcory.

<source lang="bash">
set FIDS=`cat fid_files.ls`

foreach I (`seq 1 ${#FIDS}`)
set FID=${FIDS[$I]}; set DIRN=`dirname $FID`
cd $DIRN
ls 
cd ..
end
</source>

== Extract the spectra settings from Varian procpar file ==
Now we want to make a settings file we can read in relax.

<source lang="bash">
set FIDS=`cat fid_files.ls`
set OUT=$PWD/exp_parameters.txt

echo "# DIRN I deltadof2 dpwr2slock ncyc trim ss sfrq apod_rmsd" > $OUT
foreach I (`seq 1 ${#FIDS}`)
set FID=${FIDS[$I]}; set DIRN=`dirname $FID`
cd $DIRN
set deltadof2=`awk '/^deltadof2 /{f=1;next}f{print $2;exit}' procpar`
set dpwr2slock=`awk '/^dpwr2slock /{f=1;next}f{print $2;exit}' procpar`
set ncyc=`awk '/^ncyc /{f=1;next}f{print $2;exit}' procpar`
set trim=`awk '/^trim /{f=1;next}f{print $2;exit}' procpar`
set ss=`awk '/^ss /{f=1;next}f{print $2;exit}' procpar`
set sfrq=`awk '/^sfrq /{f=1;next}f{print $2;exit}' procpar`
set apodrmsd=`showApod test.ft2 | grep "REMARK Automated Noise Std Dev in Processed Data:" | awk '{print $9}'`
echo "$DIRN $I $deltadof2 $dpwr2slock $ncyc $trim $ss $sfrq $apodrmsd" >> $OUT
cd ..
end

cat $OUT
</source>

You can check, if you have repetitions of experiments, by sorting the parameters, and see if they are dublicated.<br>
We do this, by numerical sort columns 3,4 and 5 with the values for "deltadof2, dpwr2slock, ncyc".
<source lang="bash">
sort -b -k 3,3n -k 4,4n -k 5,5n exp_parameters.txt | awk '{print $3, $4, $5}'
sort -b -k 3,3n -k 4,4n -k 5,5n exp_parameters.txt > exp_parameters_sort.txt
</source>

= Spectral process files =
== Copy data ==
We first copy the data
<source lang="bash">
# Copy data
cp -r spectrometer_data spectrometer_data_processed
cd spectrometer_data_processed
</source>

== Change format to NMRPipe ==
<source lang="bash">
set CWD=$PWD
set DIRS=`cat fid_files.ls | sed 's/\/fid//g'`
cd ${DIRS[1]}
varian
</source>

Now we make a file to convert from binary format of Varian to NMRPipe.
# Now click, 'read parameters', check 'Rance-Kay'
# Remember to set Y-'Observe Freq MHz' to N15
# Remove '''sleep 5''' from the script.
# Click 'Save script' to make '''fid.com''' file, and 'Quit', and run the script.

== Spectral processing ==
This step can be done by following wiki page [[Spectral_processing]].<br>
Start '''nmrDraw''' by command
 nmrDraw

== Convert and spectral processing all ==
Now we want to convert all spectra.<br>
You should have a '''fid.com''' and '''nmrproc.com''' in the first FID folder.<br>
We now copy these script into all of the experimental folders, and execute them.

<source lang="bash">
cd $CWD

set FIDS=`cat fid_files.ls`
set DIRN1=`dirname $PWD/${FIDS[1]}`
 
foreach I (`seq 2 ${#FIDS}`)
set FID=${FIDS[$I]}; set DIRN=`dirname $FID`
cd $DIRN
echo $DIRN
cp -f $DIRN1/fid.com .
cp -f $DIRN1/nmrproc.com .
./fid.com
./nmrproc.com
cd ..
end
</source>

== Convert NMRPipe to Sparky ==
Next we also want to convert them to SPARKY format.
<source lang="bash">
set FTS=`ls -v -d -1 */*.ft2`

foreach FT ($FTS)
    set DNAME=`dirname $FT`
    set BNAME=`basename $FT`
    set FNAME=`echo $BNAME | cut -d'.' -f1`
    echo $FT $DNAME $BNAME $FNAME
    pipe2ucsf $FT ${DNAME}/${FNAME}.ucsf
end
</source>

= Working with peaks =

== Check the peak list matches ==
Check that your peak list matches your spectrum.<br>
Read the section in [[Tutorial_for_Relaxation_dispersion_analysis_cpmg_fixed_time_recorded_on_varian_as_fid_interleaved#Check_the_peak_list_matches | Check the peak list matches]].

<source lang="bash">
set DIRS=`cat fid_files.ls | sed 's/\/fid//g'`
sparky ${DIRS[1]}/test.ucsf
</source>

The final peak list is expected to be in:
<source lang="bash">
/peak_lists/peaks_corr_final.list
</source>
And have been saved by SPARKY, so it is in this format [[SPARKY_list]].

== Check for peak movement ==
Your should check, that the peaks do not move at the different experiments. Try opening some random spectra, and overlay them in SPARKY.<br>
Read the section in [[Tutorial_for_Relaxation_dispersion_analysis_cpmg_fixed_time_recorded_on_varian_as_fid_interleaved#Check_for_peak_movement | Check for peak movement]].

<source lang="bash">
sparky ${DIRS[1]}/test.ucsf ${DIRS[10]}/test.ucsf ${DIRS[25]}/test.ucsf ${DIRS[50]}/test.ucsf
</source>

== Measuring peak heights ==
We will use the program [[NMRPipe_seriesTab | NMRPipe seriesTab]] to measure the intensities.

'''seriesTab''' needs a input file, where the ppm values from a [[SPARKY_list | SPARKY list]] has been converted to spectral points.<br>
The spectral points value depends on the spectral processing parameters.

=== Generate spectral point file ===
Create a file with spectral point information with script 
[[Tutorial_for_Relaxation_dispersion_analysis_cpmg_fixed_time_recorded_on_varian_as_fid_interleaved_scripts#stPeakList.pl | stPeakList.pl]] .

<source lang="bash">
stPeakList.pl ${DIRS[1]}/test.ft2 ../peak_lists/peaks_corr_final.list > peaks_list.tab
cat peaks_list.tab
</source>

=== Make a file name of .ft2 fil ===
<source lang="bash">
echo "test.ft2" > ft2_file.ls
</source>

=== Measure the height or sum in a spectral point box ===
<source lang="bash">
mkdir peak_lists
 
foreach line ("`tail -n+2 exp_parameters.txt`")
  set argv=( $line )
  set DIRN=$1
  set I=$2
  set deltadof2=$3
  set dpwr2slock=$4
  set ncyc=$5
  set trim=$6
  set ss=$7
  set sfrq=$8
  echo $I
  set FNAME=${I}_${deltadof2}_${dpwr2slock}_${ncyc}
  cd $DIRN
  seriesTab -in ../peaks_list.tab -out ${FNAME}_max_standard.ser -list ../ft2_file.ls -max
  seriesTab -in ../peaks_list.tab -out ${FNAME}_max_dx1_dy1.ser -list ../ft2_file.ls -max -dx 1 -dy 1
  seriesTab -in ../peaks_list.tab -out ${FNAME}_sum_dx1_dy1.ser -list ../ft2_file.ls -sum -dx 1 -dy 1
  cp ${FNAME}_max_standard.ser ../peak_lists
  cd ..
end
</source>

= Analyse in relax =

== Preparation ==
=== Prepare directory for relax run ===
Then we make a directory ready for relax
<source lang="bash">
mkdir ../relax
cp exp_parameters.txt ../relax
cp exp_parameters_sort.txt ../relax
cp -r peak_lists* ../relax
cp peaks_list.tab ../relax
cd ../relax
</source>

=== See unique parameters ===
<source lang="bash">
tail -n +2 exp_parameters.txt | awk '{print $3}' | sort -k1,1n | uniq
tail -n +2 exp_parameters.txt | awk '{print $4}' | sort -k1,1n | uniq
tail -n +2 exp_parameters.txt | awk '{print $5}' | sort -k1,1n | uniq
tail -n +2 exp_parameters.txt | awk '{print $6}' | sort -k1,1n | uniq
tail -n +2 exp_parameters.txt | awk '{print $7}' | sort -k1,1n | uniq
tail -n +2 exp_parameters.txt | awk '{print $8}' | sort -k1,1n | uniq
</source>

== Scripts ==
=== 1_setup_r1rho.py ===
This a script file to be able to call the setup.

file: '''1_setup_r1rho.py'''.
<source lang="Python">
# Python module imports.
from os import getcwd, sep

# relax module imports.
from data_store import Relax_data_store; ds = Relax_data_store()

#########################################
#### Setup
# The pipe names.
if not (hasattr(ds, 'pipe_name') and hasattr(ds, 'pipe_bundle') and hasattr(ds, 'pipe_type')):
    # Set pipe name, bundle and type.
    ds.pipe_name = 'base pipe'
    ds.pipe_bundle = 'relax_disp'
    ds.pipe_type = 'relax_disp'

# The data path
if not hasattr(ds, 'data_path'):
    ds.data_path = getcwd()

#########################################
### Start setup
# Create the data pipe.
pipe.create(pipe_name=ds.pipe_name, bundle=ds.pipe_bundle, pipe_type=ds.pipe_type)

# Read the spins.
spectrum.read_spins(file='1_0_46_0_max_standard.ser', dir=ds.data_path+sep+'peak_lists')

# Name the isotope for field strength scaling.
spin.isotope(isotope='15N')

# Load the experiments settings file.
expfile = open(ds.data_path+sep+'exp_parameters_sort.txt', 'r')
expfileslines = expfile.readlines()
expfile.close()

# In MHz
yOBS = 81.050
# In ppm
yCAR = 118.078
centerPPM_N15 = yCAR

## Read the chemical shift data.
chemical_shift.read(file='1_0_46_0_max_standard.ser', dir=ds.data_path+sep+'peak_lists')

## The lock power to field, has been found in an calibration experiment.
spin_lock_field_strengths_Hz = {'35': 431.0, '39': 651.2, '41': 800.5, '43': 984.0, '46': 1341.11, '48': 1648.5}

## Apply spectra settings.
for i in range(len(expfileslines)):
    line = expfileslines[i]
    if line[0] == "#":
        continue
    else:
        # DIRN I deltadof2 dpwr2slock ncyc trim ss sfrq
        DIRN = line.split()[0]
        I = int(line.split()[1])
        deltadof2 = line.split()[2]
        dpwr2slock = line.split()[3]
        ncyc = int(line.split()[4])
        trim = float(line.split()[5])
        ss = int(line.split()[6])
        set_sfrq = float(line.split()[7])
        apod_rmsd = float(line.split()[8])
        spin_lock_field_strength = spin_lock_field_strengths_Hz[dpwr2slock]

        # Calculate spin_lock time
        time_sl = 2*ncyc*trim

        # Define file name for peak list.
        FNAME = "%s_%s_%s_%s_max_standard.ser"%(I, deltadof2, dpwr2slock, ncyc)
        sp_id = "%s_%s_%s_%s"%(I, deltadof2, dpwr2slock, ncyc)

        # Load the peak intensities.
        spectrum.read_intensities(file=FNAME, dir=ds.data_path+sep+'peak_lists', spectrum_id=sp_id, int_method='height')

        # Set the peak intensity errors, as defined as the baseplane RMSD.
        spectrum.baseplane_rmsd(error=apod_rmsd, spectrum_id=sp_id)

        # Set the relaxation dispersion experiment type.
        relax_disp.exp_type(spectrum_id=sp_id, exp_type='R1rho')

        # Set The spin-lock field strength, nu1, in Hz
        relax_disp.spin_lock_field(spectrum_id=sp_id, field=spin_lock_field_strength)

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

        # Set The spin-lock offset, omega_rf, in ppm.
        relax_disp.spin_lock_offset(spectrum_id=sp_id, offset=omega_rf_ppm)

        # Set the relaxation times (in s).
        relax_disp.relax_time(spectrum_id=sp_id, time=time_sl)

        # Set the spectrometer frequency.
        spectrometer.frequency(id=sp_id, frq=set_sfrq, units='MHz')


# Read the R1 data
# We do not read the R1 data, but rather with R1.
# relax_data.read(ri_id='R1', ri_type='R1', frq=cdp.spectrometer_frq_list[0], file='R1_fitted_values.txt', dir=data_path, mol_name_col=1, res_num_col=2, res_name_col=3, spin_num_col=4, spin_name_col=5, data_col=6, error_col=7)
</source>

=== 2_pre_run_r2eff.py ===
This a script file to run the R2eff values only, with a high number of Monte Carlo simulations.

file: '''2_pre_run_r2eff.py'''.
<source lang="Python">
# Python module imports.
from os import getcwd, sep
import re

# relax module imports.
from auto_analyses.relax_disp import Relax_disp
from data_store import Relax_data_store; ds = Relax_data_store()
from specific_analyses.relax_disp.variables import MODEL_R2EFF


#########################################
#### Setup
# The data path
if not hasattr(ds, 'data_path'):
    ds.data_path = getcwd()

# The models to analyse.
if not hasattr(ds, 'models'):
    ds.models = [MODEL_R2EFF]

# The number of increments per parameter, to split up the search interval in grid search.
if not hasattr(ds, 'grid_inc'):
    ds.grid_inc = 21

# The number of Monte-Carlo simulations, for the error analysis in the 'R2eff' model when exponential curves are fitted.
# For estimating the error of the fitted R2eff values,
# a high number should be provided. Later the high quality R2eff values will be read for subsequent model analyses.
if not hasattr(ds, 'exp_mc_sim_num'):
    ds.exp_mc_sim_num = 10

# The result directory.
if not hasattr(ds, 'results_dir'):
    ds.results_dir = getcwd() + sep + 'results_R2eff'

## The optimisation function tolerance.
## This is set to the standard value, and should not be changed.
#if not hasattr(ds, 'opt_func_tol'):
#    ds.opt_func_tol = 1e-25
#Relax_disp.opt_func_tol = ds.opt_func_tol

#if not hasattr(ds, 'opt_max_iterations'):
#    ds.opt_max_iterations = int(1e7)
#Relax_disp.opt_max_iterations = ds.opt_max_iteration


#########################################
### Run script with setup.
script(file='1_setup_r1rho.py', dir=ds.data_path)

# To speed up the analysis, only select a few spins.
deselect.all()

# Load the experiments settings file.
residues = open(ds.data_path+sep+'global_fit_residues.txt', 'r')
residueslines = residues.readlines()
residues.close()

# Split the line string into number and text.
r = re.compile("([a-zA-Z]+)([0-9]+)([a-zA-Z]+)(-)([a-zA-Z]+)")

for i, line in enumerate(residueslines):
    if line[0] == "#":
        continue
    else:
        re_split = r.match(line)
        #print re_split.groups()
        resn = re_split.group(1)
        resi = int(re_split.group(2))
        isotope = re_split.group(3)

        select.spin(spin_id=':%i@%s'%(resi, isotope), change_all=False)

# Run the analysis.
Relax_disp(pipe_name=ds.pipe_name, pipe_bundle=ds.pipe_bundle, results_dir=ds.results_dir, models=ds.models, grid_inc=ds.grid_inc, exp_mc_sim_num=ds.exp_mc_sim_num)
</source>

=== 3_analyse_models.py ===
This a script file to analyse the models.

file: '''3_analyse_models.py'''.
<source lang="Python">
# Python module imports.
from os import getcwd, sep
import re

# relax module imports.
from auto_analyses.relax_disp import Relax_disp
from data_store import Relax_data_store; ds = Relax_data_store()
from specific_analyses.relax_disp.variables import MODEL_R2EFF, MODEL_NOREX_R1RHO_FIT_R1, MODEL_DPL94_FIT_R1, MODEL_TP02_FIT_R1, MODEL_TAP03_FIT_R1, MODEL_MP05_FIT_R1

#########################################
#### Setup
# The pipe names.
if not (hasattr(ds, 'pipe_name') and hasattr(ds, 'pipe_bundle') and hasattr(ds, 'pipe_type')):
    # Set pipe name, bundle and type.
    ds.pipe_name = 'base pipe'
    ds.pipe_bundle = 'relax_disp'
    ds.pipe_type = 'relax_disp'

# The data path
if not hasattr(ds, 'data_path'):
    ds.data_path = getcwd()

# The models to analyse.
if not hasattr(ds, 'models'):
    #ds.models = [MODEL_NOREX_R1RHO_FIT_R1, MODEL_DPL94_FIT_R1, MODEL_TP02_FIT_R1, MODEL_TAP03_FIT_R1, MODEL_MP05_FIT_R1]
    ds.models = [MODEL_NOREX_R1RHO_FIT_R1, MODEL_DPL94_FIT_R1]

# The number of increments per parameter, to split up the search interval in grid search.
if not hasattr(ds, 'grid_inc'):
    ds.grid_inc = 10

# The number of Monte-Carlo simulations for estimating the error of the parameters of the fitted models.
if not hasattr(ds, 'mc_sim_num'):
    ds.mc_sim_num = 10

# The model selection technique. Either: 'AIC', 'AICc', 'BIC'
if not hasattr(ds, 'modsel'):
    ds.modsel = 'AIC'

# The previous result directory with R2eff values.
if not hasattr(ds, 'pre_run_dir'):
    ds.pre_run_dir = getcwd() + sep + 'results_R2eff' + sep + 'R2eff'

# The result directory.
if not hasattr(ds, 'results_dir'):
    ds.results_dir = getcwd() + sep + 'results_models'

## The optimisation function tolerance.
## This is set to the standard value, and should not be changed.
#if not hasattr(ds, 'opt_func_tol'):
#    ds.opt_func_tol = 1e-25
#Relax_disp.opt_func_tol = ds.opt_func_tol

#if not hasattr(ds, 'opt_max_iterations'):
#    ds.opt_max_iterations = int(1e7)
#Relax_disp.opt_max_iterations = ds.opt_max_iteration

#########################################
# Create the data pipe.
pipe.create(pipe_name=ds.pipe_name, bundle=ds.pipe_bundle, pipe_type=ds.pipe_type)

# Load the previous results into the base pipe.
results.read(file='results', dir=ds.pre_run_dir)

# Run the analysis.
Relax_disp(pipe_name=ds.pipe_name, pipe_bundle=ds.pipe_bundle, results_dir=ds.results_dir, models=ds.models, grid_inc=ds.grid_inc, mc_sim_num=ds.mc_sim_num, modsel=ds.modsel)
</source>

=== 4_inspect_results.py ===
This a script file to inspect results in relax.

file: '''4_inspect_results.py'''.
<source lang="Python">
# Python module imports.
from os import getcwd, sep
import re

# relax module imports.
from pipe_control.mol_res_spin import generate_spin_string, return_spin, spin_loop
from specific_analyses.relax_disp.data import generate_r20_key, loop_exp_frq
from specific_analyses.relax_disp.variables import MODEL_R2EFF, MODEL_NOREX_R1RHO_FIT_R1, MODEL_DPL94_FIT_R1, MODEL_TP02_FIT_R1, MODEL_TAP03_FIT_R1, MODEL_MP05_FIT_R1

#########################################
#### Setup
results_dir = getcwd() + sep + 'results_models'

# Load the previous state
state.load(state='final_state.bz2', dir=results_dir)

# Display all pipes
pipe.display()

# Define models which have been analysed.
#MODELS = [MODEL_NOREX_R1RHO_FIT_R1, MODEL_DPL94_FIT_R1, MODEL_TP02_FIT_R1, MODEL_TAP03_FIT_R1, MODEL_MP05_FIT_R1]
MODELS = [MODEL_NOREX_R1RHO_FIT_R1, MODEL_DPL94_FIT_R1]

# Print results for each model.
print("\n################")
print("Printing results")
print("################\n")

for model in MODELS:
    # Skip R2eff model.
    if model == MODEL_R2EFF:
        continue

    # Switch to pipe.
    pipe.switch(pipe_name='%s - relax_disp' % (model))
    print("\nModel: %s" % (model))

    # Loop over the spins.
    for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
        # Generate spin string.
        spin_string = generate_spin_string(spin=cur_spin, mol_name=mol_name, res_num=resi, res_name=resn)

        # Loop over the parameters.
        print("\nOptimised parameters for spin: %s" % (spin_string))
        for param in cur_spin.params + ['chi2']:
            # Get the value.
            if param in ['r1_fit', 'r2']:
                for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True):
                    # Generate the R20 key.
                    r20_key = generate_r20_key(exp_type=exp_type, frq=frq)

                    # Get the value.
                    value = getattr(cur_spin, param)[r20_key]

                    # Print value.
                    print("%-10s %-6s %-6s %3.8f" % ("Parameter:", param, "Value:", value))

            # For all other parameters.
            else:
                # Get the value.
                value = getattr(cur_spin, param)

                # Print value.
                print("%-10s %-6s %-6s %3.8f" % ("Parameter:", param, "Value:", value))


# Print the final pipe.
pipe.switch(pipe_name='%s - relax_disp' % ('final'))
print("\nFinal pipe")

# Loop over the spins.
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
    # Generate spin string.
    spin_string = generate_spin_string(spin=cur_spin, mol_name=mol_name, res_num=resi, res_name=resn)

    # Loop over the parameters.
    print("\nOptimised model for spin: %s" % (spin_string))
    param = 'model'

    # Get the value.
    value = getattr(cur_spin, param)
    print("%-10s %-6s %-6s %6s" % ("Parameter:", param, "Value:", value))
</source>

= See also =
[[Category:Tutorials]]