Tutorial for sorting data stored as numpy to on-resonance R1rho analysis

From relax wiki
Jump to: navigation, search

Data background

This is data recorded at 600 and 950 MHz.
This should follow relax-4.0.1.Darwin.dmg installation on a Mac, only with GUI.

For each spectrometer frequency, the data is saved in np.arrays

  1. one for the residue number,
  2. one for the rates,
  3. one for the errorbars,
  4. one for the RF field strength.

They can be retrieved also with scipy's loadmat command.

The experiments are on-resonance R, and the rates are already corrected for the (small) offset effect, using the experimentally determined R1.

Specifically, the numpy shapes of the data is:

  1. For 600 MHz
    1. residues (1, 60)
    2. rates (60, 10)
    3. errorbars_rate (60, 10)
    4. RFfields (1, 10)
  1. For 950 Mhz
    1. residues (1, 61)
    2. rates (61, 19)
    3. errorbars_rate (61, 19)
    4. RFfields (1, 19)

For the 950 MHz, the 4000 Hz dispersion point is replicated.


An example of the data at the 2 fields is:

Create data files for relax

First prepare data, by running in Python:

python 1_prepare_data.py

The script is:

Run analysis in relax GUI

  • Start relax
  • Then click: new analysis or Cmd+n
  • Then click icon for Relaxation dispersionNext
  • Just accept name for the pipe
  • Open the interpreter: Cmd+p
  • Paste in
import os; os.chdir(os.getenv('HOME') + os.sep + 'Desktop' + os.sep + 'temp'); pwd()
  • Then do
script(file='2_load_data.py')
  • ( You can scroll through earlier commands with: Cmd+ )
  • Close the Grace Results viewer window
  • Close the interpreter window

Now save the current state! This saved state can now be loaded just before an analysis, and contain all setup and data.

  • File  → Save as (Shift+Ctrl+S) → Save as before_analysis.bz2

Now Quit relax, and start relax again.
Go to: File  → Open relax state (Ctrl+O) → Locate before_analysis.bz2

In the GUI, select:

  • For models, only select: No Rex and M61
  • R1 parameter optimisation to False, as this on-resonance data does not use R1.
  • The number of Monte-Carlo simulations to 3. (The minimum number before relax refuse to run).
  • Let other things be standard, and click execute


Note  The number of Monte-Carlo simulations is set to 3. This means that relax will do the analysis, and the fitted parameters will be correct, but the standard error of the parameters will be wrong.

The reason for this is, that the data "naturally"? for some spins contains measurements which are "doubtful".
When relax is performing Monte-Carlo simulations, all the data is first copied x times to the number of Monte-Carlo simulations, and each R2,eff point on the graph is modified randomly by a gaussian noise with a width described by the associated errors.

For spins which contains measurements which are "doubtful", relax will be spending "very long time" trying to fit a meaningful model.
This "time" is poorly spend on "bad data".

Rather, one should first try to

  • analyse the data quickly
  • make the graphs
  • examine which spins should be deselected
  • and then rerun the analysis with a higher number of monte-carlo simulations

After this comes the part with global fit/clustered analysis.

  • analyse the data again
  • Select which spins should be included in a global fit/clustered analysis
  • and then rerun with analysis

The script

Run the same analysis in relax through terminal

The analysis can be performed by

python 1_prepare_data.py
 
relax 2_load_data.py -t log.txt
# Or
mpirun -np 8 relax --multi='mpi4py' 2_load_data.py -t log.txt

Inspect and make graphs

  • When this is done, quit relax
  • Then to convert all xmgrace files to png

First make graphs

cd grace
./grace2images.py -t EPS,PNG
 
cd ..
cd No_Rex/
./grace2images.py -t EPS,PNG
cat r1rho_prime.out
 
cd ..
cd M61
./grace2images.py -t PNG
cat r1rho_prime.out
cat kex.out
cat phi_ex.out 
 
cat kex.out | grep "None          12"
cat phi_ex.out | grep "None          12"

This should give the images of the data

For model M61, spin 12 is:

Further analysis steps

Load the final state

  • Inspect which residues should be included in analysis
# Now simulate that all spins are first deselected, and then selected.
deselect.all()
sel_ids = [
":12",
":51",
]
for sel_spin in sel_ids:
    print("Selecting spin %s"%sel_spin)
    select.spin(spin_id=sel_spin, change_all=False)
 
# Inspect which residues should be analysed together for a clustered/global fit.
cluster_ids = sel_ids
 
# Cluster spins
for curspin in cluster_ids:
    print("Adding spin %s to cluster"%curspin)
    relax_disp.cluster('model_cluster', curspin)
 
# Show the pipe
print("\nPrinting all the available pipes.")
pipe.display()
 
# Get the selected models
print("\nChecking which model is stored per spin.")
for curspin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=False):
    print("For spin_id '%s the model is '%s''"%(spin_id, curspin.model))
 
# Copy pipe and switch to it.
pipe.copy(pipe_from="final - relax_disp", pipe_to="relax_disp_cluster", bundle_to="relax_disp_cluster")
pipe.switch(pipe_name="relax_disp_cluster")
pipe.display()
 
# Go again with clustered spins.
relax_disp.Relax_disp(pipe_name="relax_disp_cluster", pipe_bundle="relax_disp_cluster", results_dir=outdir+sep+"cluster", models=MODELS, grid_inc=GRID_INC, mc_sim_num=MC_NUM, exp_mc_sim_num=None, modsel=MODSEL,  pre_run_dir=None, optimise_r2eff=False, insignificance=0.0, numeric_only=False, mc_sim_all_models=False, eliminate=True, set_grid_r20=set_grid_r20, r1_fit=r1_fit)
 
# Get the clustered fitted values
print("\nChecking which value is stored per spin.")
kex = None
for curspin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=False):
    if kex == None:
        kex = curspin.kex
    self.assertEqual(curspin.kex, kex)
    print("For spin_id %s the kex is %.3f"%(spin_id, kex))
  • Run analysis again

See also