Difference between revisions of "Tutorial for Relaxation dispersion analysis r1rho fixed time recorded on varian as sequential spectra"

From relax wiki
Jump to navigation Jump to search
(Using the {{caution}} template.)
 
(21 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 +
{{caution|This tutorial is incomplete.}}
 +
 
= Intro =
 
= Intro =
This tutorial presently cover the [http://svn.gna.org/svn/relax/branches/ relax_disp branch].<br>
+
 
This branch is under development, for testing it out, you need to use the source code. See [[Installation_linux#Checking_out_a_relax_branch]].
 
  
 
This tutorial is based on the analysis of R1rho data, analysed in a master thesis.
 
This tutorial is based on the analysis of R1rho data, analysed in a master thesis.
Line 219: Line 220:
  
 
== Preparation ==
 
== Preparation ==
=== making a spin file from SPARKY list ===
 
relax does not yet has the possibility to read spins from a sparky file. [https://gna.org/support/?3044 See support request].
 
 
So we create one.
 
 
<source lang="bash">
 
set ATOMS=`tail -n+4 peaks_list.tab | awk '{print $7}'`
 
set SCRIPT=relax_2_spins.py
 
 
foreach I (`seq 1 ${#ATOMS}`)
 
set ATOM=${ATOMS[$I]}; set SPIN=`echo $ATOM | sed -e "s/N-HN//g"`; set RESN=`echo $SPIN | sed -e "s/[0-9]*//g"`; set RESI=`echo $SPIN | sed -e "s/[A-Za-z]//g"`
 
echo $ATOM $SPIN $RESN $RESI
 
echo "spin.create(spin_name='N', spin_num=$I, res_name='$RESN', res_num=$RESI, mol_name=None)" >> $SCRIPT
 
end
 
 
cat $SCRIPT
 
</source>
 
 
 
=== Prepare directory for relax run ===
 
=== Prepare directory for relax run ===
 
Then we make a directory ready for relax
 
Then we make a directory ready for relax
Line 245: Line 228:
 
cp -r peak_lists* ../relax
 
cp -r peak_lists* ../relax
 
cp peaks_list.tab ../relax
 
cp peaks_list.tab ../relax
cp relax_2_spins.py ../relax
 
 
cd ../relax
 
cd ../relax
 
</source>
 
</source>
Line 259: Line 241:
 
</source>
 
</source>
  
== Analyse rx ==
+
== Scripts ==
 +
=== 1_setup_r1rho.py ===
 +
This a script file to be able to call the setup.
  
=== Setup rx data ===
+
file: '''1_setup_r1rho.py'''.
'''rx_1_ini.py'''
+
<source lang="Python">
<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.
 
# Create the data pipe.
pipe_name = 'base pipe'
+
pipe.create(pipe_name=ds.pipe_name, bundle=ds.pipe_bundle, pipe_type=ds.pipe_type)
pipe_bundle = 'relax_fit'
+
 
pipe.create(pipe_name=pipe_name, bundle=pipe_bundle, pipe_type='relax_fit')
+
# 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)
  
# Create the spins
+
        # Calculating the spin-lock offset in ppm, from offsets values provided in Hz.
script(file='relax_2_spins.py', dir=None)
+
        frq_N15_Hz = yOBS * 1E6
 +
        offset_ppm_N15 = float(deltadof2) / frq_N15_Hz * 1E6
 +
        omega_rf_ppm = centerPPM_N15 + offset_ppm_N15
  
# Set the spectra experimental properties/settings.
+
        # Set The spin-lock offset, omega_rf, in ppm.
script(file='rx_3_spectra_settings.py', dir=None)
+
        relax_disp.spin_lock_offset(spectrum_id=sp_id, offset=omega_rf_ppm)
  
# Save the program state before run.
+
        # Set the relaxation times (in s).
# This state file will also be used for loading, before a later cluster/global fit analysis.
+
        relax_disp.relax_time(spectrum_id=sp_id, time=time_sl)
state.save('ini_setup_rx', force=True)
+
 
 +
        # 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>
 
</source>
  
'''rx_3_spectra_settings.py'''
+
=== 2_pre_run_r2eff.py ===
<source lang="python">
+
This a script file to run the R2eff values only, with a high number of Monte Carlo simulations.
# Set settings for experiment
+
 
import os
+
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()
  
spin_locks = [0.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0]
+
# The models to analyse.
lock_powers = [35.0, 39.0, 41.0, 43.0, 46.0, 48.0]
+
if not hasattr(ds, 'models'):
ncycs = [0, 4, 10, 14, 20, 40]
+
    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 = 2000
 +
 
 +
# 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.
 
# Load the experiments settings file.
expfile = open('exp_parameters_sort.txt','r')
+
residues = open(ds.data_path+sep+'global_fit_residues.txt', 'r')
expfileslines = expfile.readlines()
+
residueslines = residues.readlines()
fpipe = open('pipe_names_rx.txt','w')
+
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, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05
 +
 
 +
#########################################
 +
#### 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, MODEL_MP05, MODEL_DPL94, MODEL_TP02, MODEL_TAP03]
 +
    ds.models = [MODEL_NOREX_R1RHO, MODEL_DPL94]
 +
 
 +
# 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)
 +
 
 +
# If R1 is not measured, then do R1 fitting.
 +
r1_fit=True
 +
 
 +
# 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, r1_fit=r1_fit)
 +
</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, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05
 +
 
 +
#########################################
 +
#### 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 MODEL_MP05, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05]
 +
MODELS = [MODEL_NOREX_R1RHO, MODEL_DPL94]
 +
 
 +
# Print results for each model.
 +
print("\n################")
 +
print("Printing results")
 +
print("################\n")
 +
 
 +
# Store all the pipe names.
 +
pipes = []
 +
 
 +
for model in MODELS:
 +
    # Skip R2eff model.
 +
    if model == MODEL_R2EFF:
 +
        continue
 +
 
 +
    # Switch to pipe.
 +
    pipe_name = '%s - relax_disp' % (model)
 +
    pipes.append(pipe_name)
 +
    pipe.switch(pipe_name=pipe_name)
 +
    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 spin_lock in spin_locks:
+
             # For all other parameters.
    for lock_power in lock_powers:
 
        for i in range(len(expfileslines)):
 
             line=expfileslines[i]
 
            if line[0] == "#":
 
                continue
 
 
             else:
 
             else:
                 # DIRN I deltadof2 dpwr2slock ncyc trim ss sfrq
+
                 # Get the value.
                DIRN = line.split()[0]
+
                 value = getattr(cur_spin, param)
                I = int(line.split()[1])
 
                deltadof2 = float(line.split()[2])
 
                dpwr2slock = float(line.split()[3])
 
                 ncyc = int(line.split()[4])
 
                trim = float(line.split()[5])
 
                ss = int(line.split()[6])
 
                set_sfrq = float(line.split()[7])
 
  
                 # Calculate spin_lock time
+
                 # Print value.
                 time_sl = 2*ncyc*trim
+
                 print("%-10s %-6s %-6s %3.8f" % ("Parameter:", param, "Value:", value))
  
                # Define file name for peak list.
 
                FNAME = "%s_%s_%s_%s_max_standard.ser"%(I, int(deltadof2), int(dpwr2slock), ncyc)
 
                sp_id = "%s_%s_%s_%s"%(I, int(deltadof2), int(dpwr2slock), ncyc)
 
  
                if deltadof2 == spin_lock and dpwr2slock == lock_power and ncyc == ncycs[0]:
+
# Print the final pipe.
                    pipename = "%s_%s"%(int(deltadof2), int(dpwr2slock))
+
pipe.switch(pipe_name='%s - relax_disp' % ('final'))
                    pipebundle = 'relax_fit'
+
print("\nFinal pipe")
                    pipe.copy(pipe_from='base pipe', pipe_to=pipename, bundle_to=pipebundle)
 
                    pipe.switch(pipe_name=pipename)
 
                    fpipe.write("%s %s"%(pipename, pipebundle)+"\n")
 
  
                #print spin_lock, deltadof2, lock_power, dpwr2slock
+
# Loop over the spins.
                if deltadof2 == spin_lock and dpwr2slock == lock_power:
+
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
                    # Load the peak intensities (first the backbone NH, then the tryptophan indole NH).
+
    # Generate spin string.
                    spectrum.read_intensities(file=FNAME, dir=os.path.join(os.getcwd(),"peak_lists"), spectrum_id=sp_id, int_method='height')
+
    spin_string = generate_spin_string(spin=cur_spin, mol_name=mol_name, res_num=resi, res_name=resn)
  
                    # Set the relaxation times.
+
    # Loop over the parameters.
                    relax_fit.relax_time(time=time_sl, spectrum_id=sp_id)
+
    print("\nOptimised model for spin: %s" % (spin_string))
 +
    param = 'model'
  
                    # Set the peak intensity errors, as defined as the baseplane RMSD.
+
    # Get the value.
                    spectrum.baseplane_rmsd(error=1.33e+03, spectrum_id=sp_id)
+
    value = getattr(cur_spin, param)
                   
+
    print("%-10s %-6s %-6s %6s" % ("Parameter:", param, "Value:", value))
  
expfile.close()
+
 
fpipe.close()
+
# Print the model selection
 +
print("Printing the model selection")
 +
model_selection(method='AIC', modsel_pipe='test', pipes=pipes)
 +
pipe.display()
 
</source>
 
</source>
  
Run by
+
=== 5_clustered_analyses.py ===
<source lang="bash">
+
This a script file to do a clustered analysis.
relax rx_1_ini.py -t log_rx_1_ini.log
+
 
</source>
+
file: '''5_clustered_analyses.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 pipe_control.mol_res_spin import spin_loop
 +
from specific_analyses.relax_disp.variables import MODEL_R2EFF, MODEL_NOREX_R1RHO, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05
 +
 
 +
#########################################
 +
#### Setup
 +
# The pipe names.
 +
if not (hasattr(ds, 'pipe_name') and hasattr(ds, 'pipe_bundle') and hasattr(ds, 'pipe_type') and hasattr(ds, 'pipe_bundle_cluster')):
 +
    # Set pipe name, bundle and type.
 +
    ds.pipe_name = 'base pipe'
 +
    ds.pipe_bundle = 'relax_disp'
 +
    ds.pipe_type = 'relax_disp'
 +
    ds.pipe_bundle_cluster = 'cluster'
 +
 
 +
# 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, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05]
 +
    ds.models = [MODEL_DPL94]
 +
 
 +
# The number of increments per parameter, to split up the search interval in grid search.
 +
# This is not used, when pointing to a previous result directory.
 +
# Then an average of the previous values will be used.
 +
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_models' + sep + ds.models[0]
 +
 
 +
# The result directory.
 +
if not hasattr(ds, 'results_dir'):
 +
    ds.results_dir = getcwd() + sep + 'results_clustering'
 +
 
 +
## 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.
 +
ini_pipe_name = '%s - %s' % (ds.models[0], ds.pipe_bundle)
 +
pipe.create(pipe_name=ini_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)
 +
 
 +
# Create a new pipe, where the clustering analysis will happen.
 +
# We will copy the pipe to get all information.
 +
pipe.copy(pipe_from=ini_pipe_name, pipe_to=ds.pipe_name, bundle_to=ds.pipe_bundle_cluster)
 +
pipe.switch(ds.pipe_name)
 +
 
 +
pipe.display()
  
=== analyse rx data ===
+
# Now cluster spins.
'''rx_4_run.py'''
+
#relax_disp.cluster('model_cluster', ":1-100")
<source lang="python">
+
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
import os
+
    # Here one could write some advanced selecting rules.
from auto_analyses.relax_fit import Relax_fit
+
    relax_disp.cluster('model_cluster', spin_id)
  
# Set settings for run.
+
# See the clustering in the current data pipe "cdp".
GRID_INC = 11; MC_NUM = 3
+
for key, value in cdp.clustering.iteritems():
 +
    print key, value
  
# Load the initial state setup
+
# Print parameter kex before copying.
state.load(state='ini_setup_rx.bz2')
+
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
 +
    print(cur_spin.kex)
  
# Load the pipe names.
+
## Make advanced parameter copy.
fpipe = open('pipe_names_rx.txt','r')
+
# It is more advanced than the value.copy user function, in that clustering is taken into account. 
pipenames = []
+
# When the destination data pipe has spin clusters defined, then the new parameter values, when required, will be taken as the median value.
for line in fpipe:
+
relax_disp.parameter_copy(pipe_from=ini_pipe_name, pipe_to=ds.pipe_name)
    pipenames.append([line.split()[0], line.split()[1]])
 
fpipe.close()
 
  
for pipename, pipebundle in pipenames:
+
# Print parameter kex after copying.
    pipe.switch(pipe_name=pipename)
+
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
     results_directory = os.path.join("rx", pipename)
+
     print(cur_spin.kex)
  
    Relax_fit(pipe_name=pipename, pipe_bundle=pipebundle, file_root='rx', results_dir=results_directory, grid_inc=GRID_INC, mc_sim_num=MC_NUM, view_plots=False)
+
pipe.display()
</source>
 
  
Run by
+
# Run the analysis.
<source lang="bash">
+
Relax_disp(pipe_name=ds.pipe_name, pipe_bundle=ds.pipe_bundle_cluster, results_dir=ds.results_dir, models=ds.models, grid_inc=ds.grid_inc, mc_sim_num=ds.mc_sim_num, modsel=ds.modsel)
relax rx_4_run.py -t log_rx_4_run.log
 
 
</source>
 
</source>
  
 
= See also =
 
= See also =
 
[[Category:Tutorials]]
 
[[Category:Tutorials]]

Latest revision as of 17:13, 6 November 2015

Caution  This tutorial is incomplete.

Intro

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

peak_lists
spectrometer_data
scripts

You can create the folders by

mkdir peak_lists spectrometer_data scripts

In the folder peak_lists should contain SPARKY list in SPARKY list format.
In the folder scripts we put scripts which help us processing the files.
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.

Establish file-overview

Make file with paths to fid files

We make a file list of filepaths to fid files.

ls -v -d -1 */fid > fid_files.ls
cat fid_files.ls

Do something in each directory

With the fid_files.ls, we can do something in each directory.
For example do a list files in each direcory.

set FIDS=`cat fid_files.ls`

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

Extract the spectra settings from Varian procpar file

Now we want to make a settings file we can read in relax.

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

You can check, if you have repetitions of experiments, by sorting the parameters, and see if they are dublicated.
We do this, by numerical sort columns 3,4 and 5 with the values for "deltadof2, dpwr2slock, ncyc".

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

Spectral process files

Copy data

We first copy the data

# Copy data
cp -r spectrometer_data spectrometer_data_processed
cd spectrometer_data_processed

Change format to NMRPipe

set CWD=$PWD
set DIRS=`cat fid_files.ls | sed 's/\/fid//g'`
cd ${DIRS[1]}
varian

Now we make a file to convert from binary format of Varian to NMRPipe.

  1. Now click, 'read parameters', check 'Rance-Kay'
  2. Remember to set Y-'Observe Freq MHz' to N15
  3. Remove sleep 5 from the script.
  4. 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.
Start nmrDraw by command

nmrDraw

Convert and spectral processing all

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

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

Convert NMRPipe to Sparky

Next we also want to convert them to SPARKY format.

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

Working with peaks

Check the peak list matches

Check that your peak list matches your spectrum.
Read the section in Check the peak list matches.

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

The final peak list is expected to be in:

/peak_lists/peaks_corr_final.list

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.
Read the section in Check for peak movement.

sparky ${DIRS[1]}/test.ucsf ${DIRS[10]}/test.ucsf ${DIRS[25]}/test.ucsf ${DIRS[50]}/test.ucsf

Measuring peak heights

We will use the program NMRPipe seriesTab to measure the intensities.

seriesTab needs a input file, where the ppm values from a SPARKY list has been converted to spectral points.
The spectral points value depends on the spectral processing parameters.

Generate spectral point file

Create a file with spectral point information with script stPeakList.pl .

stPeakList.pl ${DIRS[1]}/test.ft2 ../peak_lists/peaks_corr_final.list > peaks_list.tab
cat peaks_list.tab

Make a file name of .ft2 fil

echo "test.ft2" > ft2_file.ls

Measure the height or sum in a spectral point box

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

Analyse in relax

Preparation

Prepare directory for relax run

Then we make a directory ready for relax

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

See unique parameters

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

Scripts

1_setup_r1rho.py

This a script file to be able to call the setup.

file: 1_setup_r1rho.py.

# 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)

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.

# 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 = 2000

# 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)

3_analyse_models.py

This a script file to analyse the models.

file: 3_analyse_models.py.

# 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, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05

#########################################
#### 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, MODEL_MP05, MODEL_DPL94, MODEL_TP02, MODEL_TAP03]
    ds.models = [MODEL_NOREX_R1RHO, MODEL_DPL94]

# 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)

# If R1 is not measured, then do R1 fitting.
r1_fit=True

# 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, r1_fit=r1_fit)

4_inspect_results.py

This a script file to inspect results in relax.

file: 4_inspect_results.py.

# 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, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05

#########################################
#### 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 MODEL_MP05, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05]
MODELS = [MODEL_NOREX_R1RHO, MODEL_DPL94]

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

# Store all the pipe names.
pipes = []

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

    # Switch to pipe.
    pipe_name = '%s - relax_disp' % (model)
    pipes.append(pipe_name)
    pipe.switch(pipe_name=pipe_name)
    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))


# Print the model selection
print("Printing the model selection")
model_selection(method='AIC', modsel_pipe='test', pipes=pipes)
pipe.display()

5_clustered_analyses.py

This a script file to do a clustered analysis.

file: 5_clustered_analyses.py.

# 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 pipe_control.mol_res_spin import spin_loop
from specific_analyses.relax_disp.variables import MODEL_R2EFF, MODEL_NOREX_R1RHO, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05

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

# 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, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05]
    ds.models = [MODEL_DPL94]

# The number of increments per parameter, to split up the search interval in grid search.
# This is not used, when pointing to a previous result directory.
# Then an average of the previous values will be used.
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_models' + sep + ds.models[0]

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

## 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.
ini_pipe_name = '%s - %s' % (ds.models[0], ds.pipe_bundle)
pipe.create(pipe_name=ini_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)

# Create a new pipe, where the clustering analysis will happen.
# We will copy the pipe to get all information.
pipe.copy(pipe_from=ini_pipe_name, pipe_to=ds.pipe_name, bundle_to=ds.pipe_bundle_cluster)
pipe.switch(ds.pipe_name)

pipe.display()

# Now cluster spins.
#relax_disp.cluster('model_cluster', ":1-100")
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
    # Here one could write some advanced selecting rules.
    relax_disp.cluster('model_cluster', spin_id)

# See the clustering in the current data pipe "cdp".
for key, value in cdp.clustering.iteritems():
    print key, value

# Print parameter kex before copying.
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
    print(cur_spin.kex)

## Make advanced parameter copy.
# It is more advanced than the value.copy user function, in that clustering is taken into account.  
# When the destination data pipe has spin clusters defined, then the new parameter values, when required, will be taken as the median value.
relax_disp.parameter_copy(pipe_from=ini_pipe_name, pipe_to=ds.pipe_name)

# Print parameter kex after copying.
for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True):
    print(cur_spin.kex)

pipe.display()

# Run the analysis.
Relax_disp(pipe_name=ds.pipe_name, pipe_bundle=ds.pipe_bundle_cluster, results_dir=ds.results_dir, models=ds.models, grid_inc=ds.grid_inc, mc_sim_num=ds.mc_sim_num, modsel=ds.modsel)

See also