mailRe: R1rho RD analysis


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

Header


Content

Posted by Troels Emtekær Linnet on June 03, 2014 - 18:35:
Hi Peixiang.

I am glad to see your interest in relax.

I joined the project around a year ago, since I know that fitting
CPMG and R1rho data in our lab, can be quite a struggle for our students.

So relax has come as a relief, where the GUI can get most things done.
But it seems very much, that you would like to get to the power of the
scripting.

The combined experimental setup for an R1rho experiment can pile up.
The combined combination of:
1) relax_disp.spin_lock_field
2) relax_disp.spin_lock_offset
3) relax_disp.relax_time

can quickly extend to 70-100+ experiments.

What I intended with the tutorial, was to imagine that you had such a
directory with 50+ experiments.
http://wiki.nmr-relax.com/Tutorial_for_Relaxation_dispersion_analysis_r1rho_fixed_time_recorded_on_varian_as_sequential_spectra#Extract_the_spectra_settings_from_Varian_procpar_file

That little bash script, will extract the different settings for each
experiment and put the settings into one column formatted settings file.
-----
# DIRN I deltadof2 dpwr2slock ncyc trim ss sfrq apod_rmsd
ldn_20120917_15N_NCBD_20mMPi_pH6c5_31C_t1rho46.fid 46 0 35 0 0.005 16
799.7773991 745.77
ldn_20120917_15N_NCBD_20mMPi_pH6c5_31C_t1rho48.fid 48 0 35 4 0.005 128
799.7773991 1274.17
ldn_20120917_15N_NCBD_20mMPi_pH6c5_31C_t1rho47.fid 47 0 35 10 0.005 16
799.7773991 732.333
....
-------

Then the idea is to measure the height of each experiment, with NMRPipe,
using seriesTab.

And as an error estimator, using just:
`showApod test.ft2 | grep "REMARK Automated Noise Std Dev in Processed
Data:" | awk '{print $9}'`.
This is the last column in the above example.

That should cover it.

This following would be a script I would suggest.

There are many ways how to organise this.
You could also make a list of lists in a script file, and just loop through
this.

If you look at the following example script.
Then the models are:
['R2eff', 'DPL94']

So, relax would exponential fit R2eff data from the intensities.

You already have all the R2eff data.
So here you could skip the intensities reading and the R2eff model.

You still need to set all the meta-data settings, to allow relax to match
the combinations of
data to a spectra id.

Then you need to read in the R2eff data, instead of the intensities.

I am actually a little unsure how that should be performed.

I think I will try to make to scripts, and see how far I get.
But haven't tried this before, and don't have test data for it.


Suggestion for Normal
------------------
from auto_analyses.relax_disp import Relax_disp

# Set pipe name, bundle and type.
pipe_name = 'base pipe'
pipe_bundle = 'relax_disp'
pipe_type= 'relax_disp'

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

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

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

# Load the experiments settings file.
expfile = open(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=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.
# Count settings
j = 0
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=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')

        # Add to counter
        j += 1

# Read the R1 data
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)

# Cluster residues
cluster_ids = [
":13@N",
":15@N",
":16@N",
":25@N",
":26@N",
":28@N",
":39@N",
":40@N",
":41@N",
":43@N",
":44@N",
":45@N",
":49@N",
":52@N",
":53@N"]

# Cluster spins
for curspin in cluster_ids:
    print("Adding spin %s to cluster"%curspin)
    relax_disp.cluster('model_cluster', curspin)

# De-select for analysis those spins who have not been clustered
for free_spin in cdp.clustering['free spins']:
    print("Deselecting free spin %s"%free_spin)
    deselect.spin(spin_id=free_spin, change_all=False)

# The dispersion models.
MODELS = ['R2eff', 'DPL94']

# The grid search size (the number of increments per dimension).
GRID_INC = 21

# The number of Monte Carlo simulations to be used for error analysis at
the end of the analysis.
MC_NUM = 3

# Model selection technique.
MODSEL = 'AIC'

# Execute the auto-analysis (fast).
# Standard parameters are: func_tol=1e-25, grad_tol=None, max_iter=10000000,
OPT_FUNC_TOL = 1e-1
relax_disp.Relax_disp.opt_func_tol = OPT_FUNC_TOL
OPT_MAX_ITERATIONS = 1000
relax_disp.Relax_disp.opt_max_iterations = OPT_MAX_ITERATIONS


# Run the analysis.
Relax_disp(pipe_name=pipe_name, pipe_bundle=pipe_bundle, results_dir=None,
models=MODELS, grid_inc=GRID_INC, mc_sim_num=MC_NUM, modsel=MODSEL)
--------------------

Suggestion for r2eff_read


#.... SAME as before
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

        # 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

        # Construct ID
        part_id = = "%s_%s_%s"%(set_sfrq, time_sl, omega_rf_ppm)
        sp_id = "%s_%s"%(part_id, spin_lock_field_strength)

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

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

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

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

        FNAME = "%s.txt"%(sp_id)
        relax_disp.r2eff_read(id=part_id, file=sp_id+'.txt', dir=None,
disp_frq=spin_lock_field_strength, spin_id_col=1, mol_name_col=2,
res_num_col=3, res_name_col=4, spin_num_col=5, spin_name_col=6, data_col=7,
error_col=8)

        # Add to counter
        j += 1

#...

# The dispersion models.
MODELS = ['DPL94']

#....

# Run the analysis.
relax_disp.Relax_disp(pipe_name=pipe_name, pipe_bundle=pipe_bundle,
results_dir=None, models=MODELS, grid_inc=GRID_INC, mc_sim_num=MC_NUM,
modsel=MODSEL)







2014-06-03 15:27 GMT+02:00 pma <Peixiang.Ma@xxxxxx>:

Hi, Edward,

Thanks, as your instruction, I tried Gui, but I could not quit selection
relaxation dispersion models. (I just report as a bug)

1) In r2eff_read, I am not quite clear about the pop window. You said I
should load partial experiment ID string as ('%s_%s'%(id,dis_point)).
Then you have the second window to load the file. what does the file mean?
The file is the fitted R2eff or a list for all the R2eff file names?

Is that possible to load a series of files? I guess so, but do you have an
example somewhere?

2) actually, I found script is more powerful, so I just use Troels
mentioned trick. type help(relax_disp.r2eff), then I found the description
about this function.
I made one script, it is included in the log file.

I attached the log file after I run the script.You can see that it reads
all the R2eff value perfectly, but did not do the calculation there.

It is a bit strange to me. . Do you have an idea why?

Do not care too much about the value in the script, such as Grid_inc,
MC_value, they are only for test now.



Peixiang



On 05/30/2014 04:45 PM, Edward d'Auvergne wrote:

Hi Peixiang,

If you do only have R2eff data, i.e. there is no original peak lists,
I'll assume that you would fit into the category of 'power user'.  In
this case in the GUI, instead of using the peak intensity wizard via
the spectra list 'Add' button, try the "User
functions->relax_disp->r2eff_read" or "User
functions->relax_disp->r2eff_read_spin" menu entries.  You will have
to deselect the 'R2eff' model.  Clicking on the 'Execute' button will
then run the analysis.

Regards,

Edward






--

Peixiang Ma, PhD.
Institut de Biologie Structurale (IBS)
UMR 5075 - CEA - CNRS - UJF
41 rue Jules Horowitz
38027 Grenoble Cedex 1, France
Email: Peixiang.Ma@xxxxxx


_______________________________________________
relax (http://www.nmr-relax.com)

This is the relax-users mailing list
relax-users@xxxxxxx

To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-users




Related Messages


Powered by MHonArc, Updated Fri Jun 06 18:00:08 2014