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]].
ls -v -d -1 */*echo "test.ft2 " > ft2_files.lscat ft2_filesft2_file.ls
OR make the sum in a box:=== See unique parameters ===
seriesTab tail -in peaks_listn +2 exp_parameters.txt | awk '{print $3}' | sort -k1,1n | uniqtail -n +2 exp_parameters.txt | awk '{print $4}' | sort -k1,1n | uniqtail -n +2 exp_parameters.txt | awk '{print $5}' | sort -k1,1n | uniqtail -n +2 exp_parameters.txt | awk '{print $6}' | sort -k1,1n | uniqtail -n +2 exp_parameters.txt | awk '{print $7}' | sort -k1,1n | uniqtail -n +2 exp_parameters.tab txt | awk '{print $8}' | sort -out peaks_list_sum_dx1_dy1k1,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 pathif 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 MHzyOBS = 81.050# In ppmyCAR = 118.078centerPPM_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 ft2_files.ls 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-sum lock field strength, nu1, in Hz relax_disp.spin_lock_field(spectrum_id=sp_id, field=spin_lock_field_strength) # Calculating the spin-dx 1 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-dy 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)
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 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
<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>
<source lang="bash">
set CWD=$PWD
set FIDSDIRS=`cat fid_files.ls| sed 's/\/fid//g'`set DIRN=`dirname cd ${FIDSDIRS[1]}`cd $DIRN
varian
</source>
<source lang="bash">
set DIRS=`cat fid_files.ls | sed 's/\/fid//g'`sparky ${FIDSDIRS[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 ==
<source lang="bash">
stPeakList.pl 0.fid${DIRS[1]}/test.ft2 ../peak_lists/peaks_corr_final.list > peaks_list.tab
cat peaks_list.tab
</source>
=== Make a file with paths to name of .ft2 files fil ===Then we make a file list of filepaths to .ft2 files.
<source lang="bash">
</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 peaks_list_max_standard${FNAME}_max_standard.ser -list ft2_files../ft2_file.ls -max seriesTab -in ../peaks_list.tab -out peaks_list_max_dx1_dy1${FNAME}_max_dx1_dy1.ser -list ft2_files../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 ../relaxcp exp_parameters.txt ../relaxcp exp_parameters_sort.txt ../relaxcp -r peak_lists* ../relaxcp peaks_list.tab ../relaxcd ../relax
</source>
<source lang="bash">
</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 = 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)
</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 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()
</source>
=== 5_clustered_analyses.py ===
This a script file to do a clustered analysis.
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()
# 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)
</source>
= See also =
[[Category:Tutorials]]
