Trusted, Bureaucrats
4,223
edits
Changes
Jump to navigation
Jump to search
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]].
=== 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>
cp relax_2_spins.py ../relax
=== Setup rx data ===file: '''rx_1_ini1_setup_r1rho.py'''.<source lang="pythonPython"># 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 pathif not hasattr(ds, 'data_path'): ds.data_path = getcwd() ############################################ Start setup
spin_locks # The models to analyse.if not hasattr(ds, 'models'): ds.models = [0MODEL_R2EFF] # The number of increments per parameter, to split up the search interval in grid search.0if not hasattr(ds, 'grid_inc'): ds.grid_inc = 21 # The number of Monte-Carlo simulations, 500for the error analysis in the 'R2eff' model when exponential curves are fitted.0# For estimating the error of the fitted R2eff values, 1000# a high number should be provided.0Later 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.0if not hasattr(ds, 5000'results_dir'): ds.results_dir = getcwd() + sep + 'results_R2eff' ## The optimisation function tolerance.0## This is set to the standard value, 10000and should not be changed.0]lock_powers #if not hasattr(ds, 'opt_func_tol'):# ds.opt_func_tol = [351e-25#Relax_disp.0, 39opt_func_tol = ds.0opt_func_tol #if not hasattr(ds, 41'opt_max_iterations'):# ds.0, 43opt_max_iterations = int(1e7)#Relax_disp.0, 46opt_max_iterations = ds.0, 48opt_max_iteration ############################################ Run script with setup.0]ncycs script(file= [0'1_setup_r1rho.py', 4dir=ds.data_path) # To speed up the analysis, 10, 14, 20, 40]only select a few spins.deselect.all()
expfile residues = open(ds.data_path+sep+'exp_parameters_sortglobal_fit_residues.txt','r')expfileslines residueslines = expfileresidues.readlines()fpipe 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, sepimport re # relax module imports.from auto_analyses.relax_disp import Relax_dispfrom 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 pathif 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 = opengetcwd() + 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, sepimport re # relax module imports.from pipe_control.mol_res_spin import generate_spin_string, return_spin, spin_loopfrom specific_analyses.relax_disp.data import generate_r20_key, loop_exp_frqfrom specific_analyses.relax_disp.variables import MODEL_R2EFF, MODEL_NOREX_R1RHO, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05 ############################################# Setupresults_dir = getcwd() + sep + 'results_models' # Load the previous statestate.load(state='final_state.bz2', dir=results_dir) # Display all pipespipe.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'pipe_names_rx]: # Get the value.txt if param in ['r1_fit','wr2']: 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 lock_power in lock_powers: for i in range(len(expfileslines)): line=expfileslines[i] if line[0] == "#": continueFor all other parameters.
# 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]: pipename = "%s_%s"%(int(deltadof2), int(dpwr2slock)) pipebundle = 'relax_fit' # Print the final pipe.copy(pipe_from='base pipe', pipe_to=pipename, bundle_to=pipebundle) pipe.switch(pipe_name=pipename) fpipe.write("'%s %s"- relax_disp' %(pipename, pipebundle'final'))+print("\nnFinal pipe")
#print spin_lockLoop over the spins.for cur_spin, deltadof2mol_name, lock_powerresi, dpwr2slock if deltadof2 resn, spin_id in spin_loop(full_info=True, return_id= spin_lock and dpwr2slock True, skip_desel== lock_powerTrue): # Load the peak intensities (first the backbone NH, then the tryptophan indole NH)Generate spin string. spectrum.read_intensities spin_string = generate_spin_string(filespin=FNAMEcur_spin, dirmol_name=os.path.join(os.getcwd(),"peak_lists")mol_name, spectrum_idres_num=sp_idresi, int_methodres_name='height'resn)
# Set Loop over the relaxation timesparameters. relax_fit.relax_time print("\nOptimised model for spin: %s" % (timespin_string)) param =time_sl, spectrum_id=sp_id)'model'
# Set Get the peak intensity errors, as defined as the baseplane RMSDvalue. spectrum.baseplane_rmsd value = getattr(error=1.33e+03cur_spin, spectrum_id=sp_idparam) print("%-10s %-6s %-6s %6s" % ("Parameter:", param, "Value:", value))
expfile.close# Print the model selectionprint("Printing the model selection")model_selection(method='AIC', modsel_pipe='test', pipes=pipes)fpipepipe.closedisplay()
Run by=== 5_clustered_analyses.py ===This a script file to do a clustered analysis. file: '''5_clustered_analyses.py'''.<source lang="bashPython"># Python module imports.from os import getcwd, sepimport re # relax rx_1_inimodule imports.from auto_analyses.relax_disp import Relax_dispfrom data_store import Relax_data_store; ds = Relax_data_store()from pipe_control.mol_res_spin import spin_loopfrom 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 pathif 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.py ini_pipe_name = '%s -t log_rx_1_ini%s' % (ds.models[0], ds.pipe_bundle)pipe.create(pipe_name=ini_pipe_name, bundle=ds.pipe_bundle, pipe_type=ds.logpipe_type)</source># 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#relax_disp.py'cluster('model_cluster', ":1-100")<source langfor cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel="python">True):import os # Here one could write some advanced selecting rules.from auto_analyses relax_disp.relax_fit import Relax_fitcluster('model_cluster', spin_id)
Relax_fitpipe.display(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)</source>
Using the {{caution}} template.
{{caution|This tutorial is incomplete.}}
= Intro =
This tutorial is based on the analysis of R1rho data, analysed in a master thesis.
set OUT=$PWD/exp_parameters.txt
echo "# DIRN I deltadof2 dpwr2slock ncyc trim ss sfrqapod_rmsd" > $OUT
foreach I (`seq 1 ${#FIDS}`)
set FID=${FIDS[$I]}; set DIRN=`dirname $FID`
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
== Preparation ==
=== Prepare directory for relax run ===
Then we make a directory ready for relax
cp -r peak_lists* ../relax
cp peaks_list.tab ../relax
cd ../relax
</source>
</source>
== Analyse rx Scripts ===== 1_setup_r1rho.py ===This a script file to be able to call the 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 MHzyOBS = 81.050# In ppmyCAR = 118.078centerPPM_N15 = yCAR ## Read the chemical shift data.chemical_shift.read(file= 'base pipe1_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}pipe_bundle ## 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= 'relax_fitheight') # Set the peak intensity errors, as defined as the baseplane RMSD.pipe spectrum.createbaseplane_rmsd(pipe_nameerror=pipe_nameapod_rmsd, bundlespectrum_id=sp_id) # Set the relaxation dispersion experiment type. relax_disp.exp_type(spectrum_id=pipe_bundlesp_id, pipe_typeexp_type='relax_fitR1rho') # Set The spin-lock field strength, nu1, in Hz relax_disp.spin_lock_field(spectrum_id=sp_id, field=spin_lock_field_strength)
# Create Calculating the spinsspin-lock offset in ppm, from offsets values provided in Hz.script(file frq_N15_Hz ='relax_2_spins.py', diryOBS * 1E6 offset_ppm_N15 =Nonefloat(deltadof2)/ frq_N15_Hz * 1E6 omega_rf_ppm = centerPPM_N15 + offset_ppm_N15
# Set the spectra experimental properties/settingsThe spin-lock offset, omega_rf, in ppm.script relax_disp.spin_lock_offset(filespectrum_id='rx_3_spectra_settings.py'sp_id, diroffset=Noneomega_rf_ppm)
# Save Set the program state before runrelaxation times (in s). relax_disp.relax_time(spectrum_id=sp_id, time=time_sl) # This state file will also be used for loadingSet the spectrometer frequency. spectrometer.frequency(id=sp_id, frq=set_sfrq, units='MHz') # Read the R1 data# We do not read the R1 data, before a later cluster/global fit analysisbut rather with R1.state# relax_data.saveread(ri_id='R1', ri_type='ini_setup_rxR1', forcefrq=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=True7)
</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: '''rx_3_spectra_settings2_pre_run_r2eff.py'''.<source lang="pythonPython"># Set settings for experimentPython module imports.from os import getcwd, sepimport osre # relax module imports.from auto_analyses.relax_disp import Relax_dispfrom data_store import Relax_data_store; ds = Relax_data_store()from specific_analyses.relax_disp.variables import MODEL_R2EFF ############################################# Setup# The data pathif not hasattr(ds, 'data_path'): ds.data_path = getcwd()
# Load the experiments settings file.
else:
# DIRN I deltadof2 dpwr2slock ncyc trim ss sfrq DIRN = line.split()[0] I = int(line.split()[1]) deltadof2 = float(line.split()[2]) dpwr2slock = float(lineGet the value.split()[3]) ncyc value = int(line.split()[4]) trim = float(line.split()[5]) ss = int(line.split()[6]) set_sfrq = float(line.splitgetattr()[7]cur_spin, param)
# Calculate spin_lock timePrint value. time_sl = 2*ncyc*trimprint("%-10s %-6s %-6s %3.8f" % ("Parameter:", param, "Value:", value))
</source>
# Set settings See the clustering in the current data pipe "cdp".for runkey, value in cdp.clustering.iteritems():GRID_INC = 11; MC_NUM = 3 print key, value
# Load the initial state setupPrint parameter kex before copying.state.loadfor cur_spin, mol_name, resi, resn, spin_id in spin_loop(statefull_info=True, return_id=True, skip_desel='ini_setup_rxTrue): print(cur_spin.bz2'kex)
# Load the pipe names# Make advanced parameter copy.fpipe = open('pipe_names_rx# It is more advanced than the value.txt'copy user function,'r')pipenames = []for line in fpipe: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.append([linerelax_disp.splitparameter_copy()[0]pipe_from=ini_pipe_name, line.split()[1]])fpipepipe_to=ds.close(pipe_name)
# Print parameter kex after copying.for pipenamecur_spin, mol_name, resi, resn, pipebundle spin_id in pipenames: pipe.switchspin_loop(pipe_namefull_info=pipenameTrue, return_id=True, skip_desel=True): results_directory = os.pathprint(cur_spin.join("rx", pipenamekex)
# Run bythe analysis.<source langRelax_disp(pipe_name=ds.pipe_name, pipe_bundle="bash">relax rx_4_runds.pipe_bundle_cluster, results_dir=ds.results_dir, models=ds.models, grid_inc=ds.grid_inc, mc_sim_num=ds.py -t log_rx_4_runmc_sim_num, modsel=ds.logmodsel)
</source>
= See also =
[[Category:Tutorials]]
