Dear Atul.
Please use the support tracker when you upload files.
The email list is fragile, and cannot handle email attachment.
I made small modifications to your script.
I normally split them up.
1st script for setup, and printing stuff to check.
2nd script, which just call the previous script, and run analysis.
So, I would:
relax 1_setup_r1rho.py
to check everything.
When I am happy, I just:
relax 2_pre_run_r2eff.py
which call script 1, setup, and then run.
Did you know you can also:
relax 1_setup_r1rho.py -p
This runs the script, and enters terminal. Like "python -i script.py"
Any variable saved in the "cdp" can be accessed.
Try:
relax 1_setup_r1rho.py -p
print cdp.temp
You can also run relax with the debug option "-d"
relax 2_pre_run_r2eff.py -d
This gives more information about where the error is.
Anyway.
The problem seems related to, that you only have one reference
intensity per offset.
----1_setup_r1rho.py
# Python module imports.
from os import sep, getcwd
# relax module imports.
from pipe_control.mol_res_spin import generate_spin_string,
return_spin, spin_loop
from specific_analyses.relax_disp.data import average_intensity,
generate_r20_key, get_curve_type, has_exponential_exp_type,
has_r1rho_exp_type, loop_exp_frq, loop_exp_frq_offset_point,
loop_exp_frq_offset_point_time, loop_time, return_grace_file_name_ini,
return_param_key_from_data
# Set up the data pipe.
#######################
# The path to the data files.
DATA_PATH = getcwd()
# Create the data pipe.
pipe_name = 'base pipe'
pipe_bundle = 'relax_disp'
pipe.create(pipe_name=pipe_name, bundle=pipe_bundle, pipe_type='relax_disp')
# Read the spins from a spectum file instead. This already contain the
spin name.
spectrum.read_spins(file='2100_off_reso_R1rho.list', dir=DATA_PATH)
# Set the isotope information.
spin.isotope(isotope='15N')
# Loop over spins, to see current setup.
for cur_spin, mol_name, resi, resn, spin_id in
spin_loop(full_info=True, return_id=True, skip_desel=True):
print mol_name, resi, resn, spin_id
print " "
print cur_spin
# The spectral data - spectrum ID, peak list file name, spin-lock
field strength (Hz), the spin-lock offset (ppm), the relaxation time
(s), spectrometer frequency (Hz), and experimental error (RMSD of the
base plane noise for each spectrum).
data = [
['ref_off_reso_R1rho', 'ref_off_reso_R1rho.list', None, 150.87,
0.320, 900.21422558574e6, 90000.0],
['2100_off_reso_R1rho', '2100_off_reso_R1rho.list', 1500.0,
141.01, 0.320, 900.21422558574e6, 90000.0],
['2728_off_reso_R1rho', '2728_off_reso_R1rho.list', 1500.0,
147.89, 0.320, 900.21422558574e6, 90000.0],
['3357_off_reso_R1rho', '3357_off_reso_R1rho.list', 1500.0,
154.78, 0.320, 900.21422558574e6, 90000.0],
['3985_off_reso_R1rho', '3985_off_reso_R1rho.list', 1500.0,
161.66, 0.320, 900.21422558574e6, 90000.0],
['4614_off_reso_R1rho', '4614_off_reso_R1rho.list', 1500.0,
168.55, 0.320, 900.21422558574e6, 90000.0],
['5242_off_reso_R1rho', '5242_off_reso_R1rho.list', 1500.0,
175.43, 0.320, 900.21422558574e6, 90000.0],
['5871_off_reso_R1rho', '5871_off_reso_R1rho.list', 1500.0,
182.32, 0.320, 900.21422558574e6, 90000.0],
['6500_off_reso_R1rho', '6500_off_reso_R1rho.list', 1500.0,
189.21, 0.320, 900.21422558574e6, 90000.0]
]
# Loop over the spectra.
for id, file, field, offset, relax_time, H_frq, rmsd in data:
print("\n###############")
print("id=%s, field=%s, offset=%3.2f, time=%1.3f, frq=%3.1f,
rmsd=%3.2f" % (id, field, offset, relax_time, H_frq/1e6, rmsd) )
# Load the peak intensities and set the errors.
spectrum.read_intensities(file=file, dir=DATA_PATH,
spectrum_id=id, int_method='height')
spectrum.baseplane_rmsd(spectrum_id=id, error=rmsd)
# Set the relaxation dispersion experiment type.
relax_disp.exp_type(spectrum_id=id, exp_type='R1rho')
# Set the relaxation dispersion spin-lock field strength (nu1).
relax_disp.spin_lock_field(spectrum_id=id, field=field)
# Set the spin-lock offset.
relax_disp.spin_lock_offset(spectrum_id=id, offset=offset)
# Set the relaxation times (in s).
relax_disp.relax_time(spectrum_id=id, time=relax_time)
# Set the NMR field strength of the spectrum.
spectrometer.frequency(id=id, frq=H_frq)
# Load the R1 data.
relax_data.read(ri_id='R1', ri_type='R1', frq=H_frq,
file='R1_relax.out', 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)
# Read the chemical shift data.
chemical_shift.read(file='2100_off_reso_R1rho.list', dir=DATA_PATH)
# Loop over spins, to see current setup.
for cur_spin, mol_name, resi, resn, spin_id in
spin_loop(full_info=True, return_id=True, skip_desel=True):
print mol_name, resi, resn, spin_id
print " "
print cur_spin
print " "
# Loop over setup.
for id in cdp.exp_type.keys():
print id, cdp.exp_type[id], cdp.spectrometer_frq[id],
cdp.spin_lock_offset[id], cdp.spin_lock_nu1[id]
print ""
cdp.temp = []
for exp_type, frq, offset, point, ei, mi, oi, di in
loop_exp_frq_offset_point(return_indices=True):
for time, ti in loop_time(exp_type=exp_type, frq=frq,
offset=offset, point=point, return_indices=True):
print exp_type, frq, offset, point, time
cdp.temp.append([exp_type, frq, offset, point, time])
------2_pre_run_r2eff.py
# Python module imports.
from os import sep, getcwd
# relax module imports.
from auto_analyses.relax_disp import Relax_disp
# Analysis variables.
#####################
# The dispersion models.
MODELS = ['R2eff']
# The grid search size (the number of increments per dimension).
GRID_INC = 11
# The number of Monte Carlo simulations to be used for error analysis
at the end of the analysis.
MC_NUM = 10
# A flag which if True will activate Monte Carlo simulations for all
models. Note this will hugely increase the computation time.
MC_SIM_ALL_MODELS = False
# The results directory.
RESULTS_DIR = getcwd() + sep + 'results_T'
# The directory of results of an earlier analysis without clustering.
PRE_RUN_DIR = None
# The model selection technique to use.
MODSEL = 'AIC'
# The flag for only using numeric models in the final model selection.
NUMERIC_ONLY = False
# The R1rho value in rad/s by which to judge insignificance. If the
maximum difference between two points on all dispersion curves for a
spin is less than this value, that spin will be deselected.
INSIGNIFICANCE = 1.0
# Run the setup script
script('1_setup_r1rho.py')
pipe_name = 'base pipe'
pipe_bundle = 'relax_disp'
# Auto-analysis execution.
##########################
# Do not change!
Relax_disp(pipe_name=pipe_name, pipe_bundle=pipe_bundle,
results_dir=RESULTS_DIR, models=MODELS, grid_inc=GRID_INC,
mc_sim_num=MC_NUM, modsel=MODSEL, insignificance=INSIGNIFICANCE,
numeric_only=NUMERIC_ONLY)
2014-08-22 23:01 GMT+02:00 Atul Srivastava <asrivast@xxxxxxx>:
Dear Edward,
Thanks for providing help and insights. I used an improved script. The
scripts is able to calculate the R2eff and the omegaeffective properly, but
it is not able to fit any model. It also sees dispersion point as a constant
1500 Hz (see attached file in results/R2eff/disp_55_N.out).
I have attached all sample data, script and results for your review. Please
see the last few lines in the log file for the error message.
Thanks,
Atul
On Fri, Aug 22, 2014 at 3:02 AM, Edward d'Auvergne <edward@xxxxxxxxxxxxx>
wrote:
Hi Atul,
Again continuing from Troels' post at
http://thread.gmane.org/gmane.science.nmr.relax.user/1718/focus=1732.
Well, actually, as Troels fully covered all of you questions in more
detail than I could have provided, I don't have anything to add. If
you have other questions, please don't hesitate to ask.
Regards,
Edward
On 21 August 2014 12:00, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:
[snip]
relax_disp.relax_time(spectrum_id=id, time=relax_time)
-> This is used for the initial R2eff calculation, for exponential
curve fitting.
# Set the NMR field strength of the spectrum.
spectrometer.frequency(id=id, frq=H_frq)
-> This is used for conversion between nucleus, etc.
# Load the R1 data.
relax_data.read(ri_id='500MHz', ri_type='R1', frq=500e6,
file='R1_500MHz.out', 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)
-> R1 needs to be loaded. This is because that R1 is part of the
equations. :-) The next release of relax, will implement feature where
R1 is fitted.
Note, that fitting introduces another variable in the equations.
# Read the chemical shift data.
chemical_shift.read(file='ref_500MHz.list', dir=DATA_PATH)
-> This is needed to get the ppm position of the nucleus.
-> See now:
http://www.nmr-relax.com/manual/Dispersion_model_summary.html
-> The "average resonance in the rotating frame" Omega(bar) = w_(bar) -
w_rf."
-> Here w_bar is the chemical shift, and w_rf is the offset.
-> Let is review figure Fig1_Palmer_Massi_2006.png
-> This is the S_z axis.
That should be it.
If you have any questions, please don't hesitate to write again.
If you need more help, consider writing a support request on the
homepage tracker.
- https://gna.org/support/?group=relax
Add following information:
# Please attach a system info file
relax -i -t relax_i.txt
# Please write up, which buttons you pushed, or attach your script.
# Consider adding your data in "sample" format. Meaning that you
delete all confidential information from the files, and only have 1-2
residues left for testing.
If you write such a support request, it it easier to share script files,
and
help other users.
The benefits from such a support request is:
- The information is available to all users, which can benefit
others in same situation.
- The information can be tracked back.
- The relax manual can be expanded, to help future users in same
situation.
Good luck!
Best
Troels Emtekær Linnet
PhD student
Copenhagen University
SBiNLab, 3-0-41
Ole Maaloes Vej 5
2200 Copenhagen N
Tlf: +45 353-22083
Lync Tlf: +45 353-30195
Re: field strength off-resonance R1rho constant relax time relaxation dispersion