Relaxation trajectory analysis example

Introduction

This example is taken from Listing S12 - S15 in the 2013 JBNMR nmrglue paper. In this example a series of 3D NMRPipe files containing relaxation trajectories for a solid state NMR experment and analyzed.

Instructions

Execute python extract_traj.py to extract the relaxation trajectories from the data set. ‘XXX.dat’ files are created for the peaks defined in the boxes.in file. The spectra.in file defines which spectra the trajectories will be extracted from.

Execute python plot_boxes.py to create plots showing the peak and integration limits for all peaks defined in boxes.in. peak_XXX_spectrum_X.png files are created for all peaks and spectra. A example plot is provided as Figure 7 of the paper, which corresponds to peak_D40_spectrum_0.png

Execute python fit_exp.py to fit all relaxation trajectories. The fitting results are provided in the fits.txt file. The relaxation_times.in file defines the relaxation times for each spectra.

Execute python plot_trajectories to create plots of all experimental and fit relaxation tranjectories. This script creates a series of XXX_plot.png files. An example plot is provided as Figure 8 of the paper, which corresponds to D40_plot.png.

The data used in this example is available for download part1, [177.28 MB] part2, [177.27 MB] part3, [177.27 MB] part4. [59.09 MB]

Listing S12

[extract_traj.py]

import nmrglue as ng
import numpy as np

# read the integration limits and list of spectra
peak_list = np.recfromtxt("boxes.in", names=True)
spectra_list = np.recfromtxt("spectra.in")

# create an array to hold the trajectories
trajectories = np.empty((peak_list.size, spectra_list.size), dtype='float')

# loop over the spectra
for sn, spectra in enumerate(spectra_list):

    # read in the spectra data
    print "Extracting peak intensities from:", spectra
    dic, data = ng.pipe.read(spectra)

    # loop over the integration limits
    for i, (name, x0, y0, x1, y1) in enumerate(peak_list):

        if x0 > x1:
            x0, x1 = x1, x0
        if y0 > y1:
            y0, y1 = y1, y0

        # integrate the region and save in trajectories array
        trajectories[i][sn] = data[y0:y1 + 1, x0:x1 + 1].sum()

# write out the trajectories for each peak
for itraj, peak_traj in enumerate(trajectories):
    peak_traj /= peak_traj.max()    # normalize the peak's trajectory
    fname = peak_list.peak_label[itraj] + '.dat'
    f = open(fname, 'w')
    for v in peak_traj:
        f.write(str(v) + '\n')
    f.close()

Listing S13

[plot_boxes.py]

import numpy as np
import nmrglue as ng
import matplotlib.pyplot as plt
import matplotlib.cm

# plot parameters
xpad = 5                        # padding around peak box on x-axis
ypad = 5                        # padding around peak box on y-axis
cmap = matplotlib.cm.Blues_r    # contour map (colors to use for contours)

# contour levels
cl = 30000 * 1.20 ** np.arange(20)

# read in the box limits and list of spectra
peak_list = np.recfromtxt("boxes.in", names=True)
spectra_list = np.recfromtxt("spectra.in")

# loop over the spectra 
for spec_number, spectra in enumerate(spectra_list):

    # read in the spectral data
    dic, data = ng.pipe.read(spectra)

    # loop over the peaks
    for peak, x0, y0, x1, y1 in peak_list:

        if x0 > x1:
            x0, x1 = x1, x0
        if y0 > y1:
            y0, y1 = y1, y0

        # slice the data around the peak
        slice = data[y0 - ypad:y1 + 1 + ypad, x0 - xpad:x1 + 1 + xpad]

        # create the figure
        fig = plt.figure()
        ax = fig.add_subplot(111)

        # plot the contours
        print "Plotting:", peak, spec_number
        extent = (x0 - xpad + 1, x1 + xpad - 1, y0 - ypad + 1, y1 + ypad - 1)
        ax.contour(slice, cl, cmap=cmap, extent=extent)

        # draw a box around the peak
        ax.plot([x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0], 'k--')

        # draw lighter boxes at +/- 1 point
        ax.plot([x0 - 1, x1 + 1, x1 + 1, x0 - 1, x0 - 1],
                [y0 - 1, y0 - 1, y1 + 1, y1 + 1, y0 - 1], 'k--', alpha=0.35)
        ax.plot([x0 + 1, x1 - 1, x1 - 1, x0 + 1, x0 + 1],
                [y0 + 1, y0 + 1, y1 - 1, y1 - 1, y0 + 1], 'k--', alpha=0.35)

        # set the title, save the figure
        ax.set_title('Peak: %s Spectrum: %i'%(peak, spec_number))
        fig.savefig('peak_%s_spectrum_%i'%(peak, spec_number))
        del(fig)

Listing S14

[fit_exp.py]

import glob

import numpy as np
from nmrglue.analysis.leastsqbound import leastsqbound

# exponential function to fit data to.
def fit_func(p, x):
    A, R2 = p
    return A * np.exp(-1.0 * np.array(x) * R2 / 1.0e6)

# residuals between fit and experimental data.
def residuals(p, y, x):
    err = y - fit_func(p, x)
    return err

# prepare fitting parameters
relaxation_times = np.loadtxt("relaxation_times.in")
x0 = [1.0, 0.10]  # initial fitting parameter
bounds = [(0.98, 1.02), (None, None)] # fitting constraints

# create an output file to record the fitting results
output = open('fits.txt', 'w')
output.write("#Peak\tA\t\tR2\t\tier\n")

# loop over the trajecory files
for filename in glob.glob('*.dat'):
    
    peak = filename[:3]
    print "Fitting Peak:", peak

    # fit the trajectory using contrainted least squares optimization
    trajectory = np.loadtxt(filename)
    x, ier = leastsqbound(residuals, x0, bounds=bounds, 
                            args=(trajectory, relaxation_times))
    
    # write fitting results to output file
    output.write('%s\t%.6f\t%.6f\t%i\n' % (peak, x[0], x[1], ier))

output.close()  # close the output file

Listing S15

[plot_trajectories.py]

import numpy as np
import matplotlib.pyplot as plt

# exponential function used to fit the data
def fit_func(p,x):
    A, R2 = p
    return A * np.exp(-1.0 * np.array(x) * R2 / 1.0e6)

fitting_results = np.recfromtxt('fits.txt')
experimental_relaxation_times = np.loadtxt("relaxation_times.in")
simulated_relaxation_times = np.linspace(0,4000000,2000)

# loop over the fitting results
for peak, A, R2, ier in fitting_results:

    print "Plotting:", peak
    
    # load the experimental and simulated relaxation trajectories
    experimental_trajectory = np.loadtxt(peak + '.dat')
    simulated_trajectory = fit_func((A, R2), simulated_relaxation_times)
    
    # create the figure
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(experimental_relaxation_times, experimental_trajectory, 'or')
    ax.plot(simulated_relaxation_times, simulated_trajectory, '-k')
    ax.set_title(peak)
    fig.savefig(peak+"_plot.png")

Example input files

[spectra.in]

data/Ytau_100.fid/test.ft2
data/Ytau_100000.fid/test.ft2
data/Ytau_250000.fid/test.ft2
data/Ytau_500000.fid/test.ft2
data/Ytau_750000.fid/test.ft2
data/Ytau_1000000.fid/test.ft2
data/Ytau_1500000.fid/test.ft2
data/Ytau_2000000.fid/test.ft2
data/Ytau_3000000.fid/test.ft2
data/Ytau_4000000.fid/test.ft2

[boxes.in]

#peak_label    X0     Y0      X0      Y1     
A20     4068    938     4079    913
A24     3992    1013    4000    997
A26     4065    962     4075    940
A34     4009    985     4018    958
A48     4028    1034    4036    1010
C28     4035    1115    4044    1092
D36     3994    987     4003    973
D40     4076    802     4085    774
D46     4155    899     4163    883
D47     4053    967     4062    941
E15     4162    1022    4170    996
E19     4176    902     4185    875
E27     4036    1084    4044    1054
E42     4136    1055    4142    1026
E56     4107    821     4115    794
F30     4013    1060    4023    1031
F52     4097    828     4105    799
G09     4054    1249    4063    1220
G14     4068    1331    4077    1304
G38     4098    1254    4106    1227
G41     4091    1283    4099    1259
I06     4087    903     4096    884

[relaxation_times.in]

# time in us
100
100000
250000
500000
750000
1000000
1500000
2000000
3000000
4000000

Example output

[fits.txt]

#Peak	A		R2		ier
F30	0.981006	0.259152	1
D40	0.980906	0.328635	1
G14	0.994899	0.108697	1
A48	0.984907	0.125793	1
E15	0.996023	0.088139	1
D47	0.993891	0.101194	1
G38	1.005751	0.127826	1
F52	1.014378	0.063087	1
G09	0.993838	0.097431	1
A34	1.002384	0.152056	1
I06	0.997428	0.060497	1
E56	0.987760	0.157588	1
D36	0.985922	0.153331	1
A20	0.999985	0.102821	1
E27	0.993052	0.298782	1
D46	0.980000	0.074505	1
C28	0.980000	0.680741	1
E42	0.996588	0.210449	1
A26	0.989291	0.258348	1
E19	0.992583	0.165457	1
G41	0.980000	0.528399	1
A24	0.992186	0.417901	1
../_images/peak_D40_spectrum_0.png ../_images/D40_plot.png