Find and fit peaks in a spectrum
This example shows how to find and fit peaks in a spectrum. This example fits the two groups of peaks in the spectrum of 1,3 diaminopropane which is available from the Madison Metabolomics Consortium Database as expnmr_00001_1.tar.
The following script will process this data creating a spectrum, perform peak picking using the pick function in nmrglue, fit the spectrum using fit_spectrum, simulate the spectrum using sim_NDregion and finally plot the results.
#! /usr/bin/env python
import nmrglue as ng
import matplotlib.pyplot as plt
# read in the bruker formatted data
dic, data = ng.bruker.read('expnmr_00001_1')
# remove the digital filter
dic['acqus']['GRPDLY'] = 72 # set the group delay to 72
data = ng.bruker.remove_digital_filter(dic, data, )
# process the spectrum
data = ng.proc_base.zf_size(data, 32768) # zero fill to 32768 points
data = ng.proc_base.fft(data) # Fourier transform
data = ng.proc_base.ps(data, p0=-50.0) # phase correction
data = ng.proc_base.di(data) # discard the imaginaries
data = ng.proc_base.rev(data) # reverse the data
data = data[20000:25000] # extract out the peaks of interest
# perform peak picking
peaks = ng.peakpick.pick(data, 25000)
npeaks = len(peaks) # 9 in this case
# fit the spectrum
lineshapes = ['g']
params = [[(x, lw)] for x, lw in zip(peaks['X_AXIS'], peaks['X_LW'])]
amps = peaks['VOL']
# bounds must be a list of length npeaks (9) with each elements of the list having a length
# of the dimensionalty of the data (1), each element of this sub-list must be a list
# of length equal to the number of lineshape parameters (2). The elements of this
# sub-sub-list must be 2-tuples containing the minimum and maximum values for the parameter.
# In this case the peak center is un-constrained and the linewidth is
# constrained to be greater than 0.
bounds = [[[(None, None), (0, None)]] for i in range(npeaks)]
ampbounds = [None for i in range(len(peaks))]
centers = [(i, ) for i in peaks['X_AXIS']]
rIDs = list(peaks['cID'])
box_width = (25,)
error_flag = False
params_best, amp_best, iers = ng.linesh.fit_spectrum(
data, lineshapes, params, amps, bounds, ampbounds, centers, rIDs,
box_width, error_flag)
# simulate the spectrum
sdata = ng.linesh.sim_NDregion(data.shape, lineshapes, params_best, amp_best)
# plot
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(data, 'k-')
ax.plot(peaks['X_AXIS'], data[peaks['X_AXIS'].astype('int')], 'ro')
ax.plot(sdata, 'b-')
ax.axhline(y=25000, color='k', ls='--')
plt.show()Note in this example all the peaks at fit individually in separate clusters which limited the size of box_width (try setting it to 50). To cluster the peaks into two sets change the peak picking line to:
peaks = ng.peakpick.pick(data, 25000, c_ndil=100)Which results in fitting two cluster of peaks, the first containing three peaks, the second five.
