jjhelmus / nmrglue

A module for working with NMR data in Python
BSD 3-Clause "New" or "Revised" License
211 stars 86 forks source link

bruker pdata: serWarning: Failed to determine udic parameters for dim: 0 #38

Closed dmaziuk closed 8 years ago

dmaziuk commented 8 years ago

(Submitting this as an issue since I don't see a user mailing list or anything like that)

I'm playing with the basic tutorial code trying to plot a topspin 1H spectrum and I'm having hard time figuring nmrglue out. What I'd like to get is this: figure_topspin

The best I managed so far is the attached figure_nmrglue, with a

nmrglue-0.6_dev-py2.7.egg/nmrglue/fileio/bruker.py:84: UserWarning: Failed to determine udic parameters for dim: 0

thrown in the process (I'm calling ng.bruker.guess_udic( dic, data ) in order to construct a ng.fileiobase.uc_from_udic( udic ) so I can uc.ppm_scale() the plot). Am I doing everything wrong or is ppm scaling not working?

Thanks in advance

atomman commented 8 years ago

Hi While it wont solve the gues_udic problem I usually construct the ppm scale from a bruker spectrum like this: dic, data = ng.bruker.read_pdata(pathFinal) SF=dic['procs']['SF'] SW_p=dic['procs']['SW_p'] SI=dic['procs']['SI'] OFFSET=dic['procs']['OFFSET'] ppm_scale=np.linspace(OFFSET,OFFSET-SW_p/SF,SI)

Den lørdag den 3. oktober 2015 skrev dmaziuk notifications@github.com:

(Submitting this as an issue since I don't see a user mailing list or anything like that)

I'm playing with the basic tutorial code trying to plot a topspin 1H spectrum and I'm having hard time figuring nmrglue out. What I'd like to get is this: http://www.bmrb.wisc.edu/ftp/pub/bmrb/metabolomics/entry_directories/bmse001105/nmr/spectra_png/1D_1H/1D_1H_0.png

The best I managed so far is the attached [image: figure_nmrglue] https://cloud.githubusercontent.com/assets/3136401/10258955/9c1fb2c8-6928-11e5-8134-d390acaf157c.png, with a

nmrglue-0.6_dev-py2.7.egg/nmrglue/fileio/bruker.py:84: UserWarning: Failed to determine udic parameters for dim: 0

thrown in the process (I'm calling ng.bruker.guess_udic( dic, data ) in order to construct a ng.fileiobase.uc_from_udic( udic ) so I can uc.ppm_scale() the plot). Am I doing everything wrong or is ppm scaling not working?

Thanks in advance

— Reply to this email directly or view it on GitHub https://github.com/jjhelmus/nmrglue/issues/38.

Sent from Gmail Mobile

dmaziuk commented 8 years ago

Thank you,

(dic, data) = ng.bruker.read_pdata("/share/subedit/metabolomics/bmse001105/nmr/set01/1/pdata/1")
data = ng.proc_base.di( data )
fig = plt.figure()
ax = fig.add_subplot(111)
SF = dic["procs"]["SF"]
SW_p = dic["procs"]["SW_p"]
SI = dic["procs"]["SI"]
OFFSET = dic["procs"]["OFFSET"]
ppm_scale = numpy.linspace( OFFSET, OFFSET - SW_p / SF, SI )
ax.plot( ppm_scale, data, "k-", linewidth = 0.3 )
fig.savefig('figure_nmrglue.png')

results in figure_nmrglue which is much better. Any tips as to plotting the Y axis in relative intensity?

atomman commented 8 years ago

Hi You should not use ng.proc_base.rev on processed data without also reversing the ppm scale before you plot it. Alternative use plt.gca().invert_xaxis() on the non-reversed data when you plot it

Also no need to delete imaginaries from pdata since there are none.

Den mandag den 5. oktober 2015 skrev dmaziuk notifications@github.com:

Thank you, but

(dic, data) = ng.bruker.read_pdata("/share/subedit/metabolomics/bmse001105/nmr/set01/1/pdata/1") data = ng.proc_base.di( data ) data = ng.proc_base.rev( data ) fig = plt.figure() ax = fig.add_subplot(111) SF = dic["procs"]["SF"] SW_p = dic["procs"]["SW_p"] SI = dic["procs"]["SI"] OFFSET = dic["procs"]["OFFSET"] ppm_scale = numpy.linspace( OFFSET, OFFSET - SW_p / SF, SI ) ax.plot( ppm_scale, data, "k-", linewidth = 0.3 ) fig.savefig('figure_nmrglue.png')

results in [image: figure_nmrglue] https://cloud.githubusercontent.com/assets/3136401/10270821/c6d9da6a-6ac3-11e5-83f7-60ad9a565c1b.png which unfortunately isn't any closer to the image generated by the spectroscopist using topspin.

— Reply to this email directly or view it on GitHub https://github.com/jjhelmus/nmrglue/issues/38#issuecomment-145400523.

Sent from Gmail Mobile

dmaziuk commented 8 years ago

I can figure out about reversing the x-axis, thank you, my problem is the examples/docs are sorely lacking.

Is there a link between data and dic?

What data point corresponds to my scaled x=0?

Just copy-pasting stuff from the examples doesn't work: with processed bruker data guess_udic() crashes. Trying to use process_and_plot_nmrglue.py example on unprocessed data instead, guess_udic() works and I can make a uc object:

udic = ng.bruker.guess_udic( dic, data )
uc = ng.fileiobase.uc_from_udic( udic )
...
ax.plot( uc.ppm_scale(), data, "k-", linewidth = 0.3 )

but then matplotlib croaks with

Traceback (most recent call last):
  File "./process_and_plot_nmrglue.py", line 31, in <module>
    ax.plot(uc.ppm_scale(), data, 'k-')
  File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 3996, in plot
    for line in self._get_lines(*args, **kwargs):
  File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 330, in _grab_next_args
    for seg in self._plot_args(remaining, kwargs):
  File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 308, in _plot_args
    x, y = self._xy_from_xy(x, y)
  File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 248, in _xy_from_xy
    raise ValueError("x and y must have same first dimension")
ValueError: x and y must have same first dimension

I'd like to make use of nmrglue but at this point I really don't see how. :(

atomman commented 8 years ago

If you have calculated the ppm scale and call I AX the the following definition lets you find a datapoint correponding to a given ppm value

definition used to convert units

def ppm2dat(AX,ppm): return (np.abs(AX-ppm)).argmin()

I use nmrglue extensively for processed bruker data and find it very useful

Den mandag den 5. oktober 2015 skrev dmaziuk notifications@github.com:

I can figure out about reversing the x-axis, thank you, my problem is the examples/docs are sorely lacking.

Is there a link between data and dic?

What data point corresponds to my scaled x=0?

Just copy-pasting stuff from the examples doesn't work: with processed bruker data guess_udic() crashes. Trying to use process_and_plot_nmrglue.py example on unprocessed data instead, guess_udic() works and I can make a uc object:

udic = ng.bruker.guess_udic( dic, data ) uc = ng.fileiobase.uc_from_udic( udic ) ... ax.plot( uc.ppm_scale(), data, "k-", linewidth = 0.3 )

but then matplotlib croaks with

Traceback (most recent call last): File "./process_and_plot_nmrglue.py", line 31, in ax.plot(uc.ppm_scale(), data, 'k-') File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 3996, in plot for line in self._get_lines(_args, *_kwargs): File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 330, in _grab_next_args for seg in self._plot_args(remaining, kwargs): File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 308, in _plot_args x, y = self._xy_from_xy(x, y) File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 248, in _xy_from_xy raise ValueError("x and y must have same first dimension") ValueError: x and y must have same first dimension

I'd like to make use of nmrglue but at this point I really don't see how. :(

— Reply to this email directly or view it on GitHub https://github.com/jjhelmus/nmrglue/issues/38#issuecomment-145627973.

Sent from Gmail Mobile

dmaziuk commented 8 years ago

Sorry, I just program computers here, I don't understand what you mean. What is AX?

atomman commented 8 years ago

AX is a calculated ppm scale and ppm2dat is a custom function that converts a ppm value to a corresponding data point.

I guess I misspelled in the earlier message...

Den mandag den 5. oktober 2015 skrev dmaziuk notifications@github.com:

Sorry, I just program computers here, I don't understand what you mean. What is AX?

— Reply to this email directly or view it on GitHub https://github.com/jjhelmus/nmrglue/issues/38#issuecomment-145652666.

Sent from Gmail Mobile

dmaziuk commented 8 years ago

Above, ppm_scale was calculated as numpy.linspace( OFFSET, OFFSET - SW_p / SF, SI ) which as far as I can tell returns an array. So you're saying that is the AX?

Let me put it this way: I'd like to a) limit the X axis to 0:8 and b) take the tallest peak, make its Y value "15" and limit the Y axis to 0:20. From looking at the plotting examples and bruker examples, I've no clue where to even begin. If I ever get this to work, this should go into the wiki... and/or the examples directory...

atomman commented 8 years ago

This will do what you want. Why would you scale the data to 20? if you want the exact scaling as in topspin you should read in the data with the scale_data option set to True.

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

def ppm2dat(AX,ppm):

function used to convert from ppm to datapoint

idx = (np.abs(AX-ppm)).argmin(axis=0)
return idx

dic, data = ng.bruker.read_pdata('./10/pdata/1') SF=dic['procs']['SF'] SW_p=dic['procs']['SW_p'] SI=dic['procs']['SI'] OFFSET=dic['procs']['OFFSET'] ppm_scale=np.linspace(OFFSET,OFFSET-SW_p/SF,SI) data_scaled=data/np.float(data[data.argmax(axis=0)])*20 min = ppm2dat(ppm_scale,0) max = ppm2dat(ppm_scale,8) if min>max: min,max = max,min

Plot and reverse the x axis

plt.plot(ppm_scale[min:max+1],data_scaled[min:max+1]) plt.gca().invert_xaxis()

On Tue, Oct 6, 2015 at 6:54 AM, dmaziuk notifications@github.com wrote:

Above, ppm_scale was calculated as numpy.linspace( OFFSET, OFFSET - SW_p / SF, SI ) which as far as I can tell returns an array. So you're saying that is the AX?

Let me put it this way: I'd like to a) limit the X axis to 0:8 and b) take the tallest peak, make its Y value "15" and limit the Y axis to 0:20. From looking at the plotting examples and bruker examples, I've no clue where to even begin. If I ever get this to work, this should go into the wiki...

— Reply to this email directly or view it on GitHub https://github.com/jjhelmus/nmrglue/issues/38#issuecomment-145741898.

jjhelmus commented 8 years ago

Just chiming in with a few additional comment.

Thanks @dmaziuk for the question and interest in nmrglue and thank you @atomman for taking the time to provide help on this topic.

There is a nmrglue mailing list which is a great place to ask question like this. The list is low traffic but does have a large reach. Neither the mailing list or the package's landing page are mentioned in the README. I'll update this shortly to include links to both.

Many of the questions you are asking are covered in the nmrglue tutorial. This document explains the data and dic variables, the use and creation of universal dictionaries, unit conversion objects, and using matplotlib to create plots.

Reading processed Bruker data is a relatively new feature in nmrglue and there may still be some issues which need to be addressed. The guess_udic function is particularly brittle. As the name states the function tries to guess a universal dictionary, which implies that it might fail. The warning you are seeing is an indication that the automatic determination has failed and the parameters will need to be provided manually.

Unfortunately, I have not had much experience with Bruker data, when I was doing NMR research I worked almost exclusively on Agilent/Varian systems, and so much of the format and variables contained in the metadata files is unfamiliar to me. I have relied on other in the community to provide expert advice on the best way to handle this data. I would welcome a contribution that improved the guess_udic function in the bruker module as improving this function would be a great help to many users.

@atomman has provided a script which can create the plot which should meet your specification. Below is another which used the unit conversion objects to determine the limits and attempt to produce a figure like the first one you provided:

import nmrglue as ng
import matplotlib.pyplot as plt

dic, data = ng.bruker.read_pdata('./set01/1/pdata/1/', scale_data=True)

# scale the data to be similar to the data in the figure provided
data *= 15. / data.max()  # scale so that the largest value is 15

# create a universal dictionary
# I am not familar enough with Bruker data to know if this is a typically how
# spectral parameters are stored in processed data
sw = dic['procs']['SW_p']  # spectal width in Hz
obs = dic['procs']['SF']   # observation frequences in MHz
carrier = dic['procs']['OFFSET'] * obs - (sw / 2.) # carrier freq in Hz

udic = ng.bruker.guess_udic(dic, data)
udic[0]['size']     = data.shape[0]
udic[0]['complex']  = True
udic[0]['encoding'] = 'direct'
udic[0]['sw']       = sw
udic[0]['obs']      = obs
udic[0]['car']      = carrier
udic[0]['label']    = '1H'
udic[0]['freq']     = True
udic[0]['time']     = False

# create a unit conversion object for the direct dimension
uc = ng.fileiobase.uc_from_udic(udic)
ppm_scale = uc.ppm_scale()

# plot the data
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ppm_scale, data, 'k-')

ax.set_xlim(6.5, 0.4)
ax.set_ylim(-0.5, 18.5)
ax.grid()

plt.savefig('figure.png')

Results:

figure

dmaziuk commented 8 years ago

OK, thank you guys.

@jjhelmus : so, why do you pick xlim = 6.5? What I mean is, there's a chloroform peak at ~7.24 ppm. So far I'm unable to figure out how to find its index as most of the data points are "nonzero" as in 0.000835907636396. Is there a way to find it from the data and dic?

jjhelmus commented 8 years ago

I choose 6.5 because that is the limit used in the figure you linked.

The raw spectral data does not contain information about peak locations. Rather the data must be analyzed to determine the locations of the peaks. nmrglue has some rudimentary peak picking methods which can be used through nmrglue.peakpick.pick.

Here is an example where the peaks in the spectrum are found, their locations in ppm determined and then plotted. The value for MINIMUM_PEAK_HEIGHT was set by visual inspection of the spectra.

from __future__ import print_function
import nmrglue as ng
import matplotlib.pyplot as plt

dic, data = ng.bruker.read_pdata('./set01/1/pdata/1/', scale_data=True)

# scale the data to be similar to the data in the figure provided
data *= 15. / data.max()  # scale so that the largest value is 15

# create a universal dictionary
sw = dic['procs']['SW_p']  # spectal width in Hz
obs = dic['procs']['SF']   # observation frequences in MHz
carrier = dic['procs']['OFFSET'] * obs - (sw / 2.) # carrier freq in Hz

udic = ng.bruker.guess_udic(dic, data)
udic[0]['size']     = data.shape[0]
udic[0]['complex']  = True
udic[0]['encoding'] = 'direct'
udic[0]['sw']       = sw
udic[0]['obs']      = obs
udic[0]['car']      = carrier
udic[0]['label']    = '1H'
udic[0]['freq']     = True
udic[0]['time']     = False

# create a unit conversion object for the direct dimension
uc = ng.fileiobase.uc_from_udic(udic)
ppm_scale = uc.ppm_scale()

# find all peaks in the spectrum
MINIMUM_PEAK_HEIGHT = 0.25
peak_table = ng.peakpick.pick(data, MINIMUM_PEAK_HEIGHT)
peak_locations = peak_table['X_AXIS'].astype('int')
peak_heights = data[peak_locations]
peak_locations_ppm = [uc.ppm(p) for p in peak_locations]

# print out the peak table
print("---------Peaks detected--------")
print("peak_location(points)\tpeak_location(ppm)\tpeak_height")
npeaks = len(peak_locations)
for peak_num in range(npeaks):
    loc_pts = peak_locations[peak_num]
    loc_ppm = peak_locations_ppm[peak_num]
    height = peak_heights[peak_num]
    print('%i\t%f\t%f' % (loc_pts, loc_ppm, height))

# plot the data
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ppm_scale, data, 'k-')  # plot the spectrum as a black line
# mark detected peaks with red circles
ax.plot(peak_locations_ppm, peak_heights, 'r.')
ax.set_xlim(8.0, -0.4)
ax.set_ylim(-0.5, 18.5)
ax.grid()

Output:

---------Peaks detected--------
peak_location(points)   peak_location(ppm)  peak_height
11200   7.262703    2.253766
15092   5.358657    0.309782
15137   5.336642    2.859392
15189   5.311203    0.290627
15367   5.224122    0.318989
18027   3.922796    0.551951
21329   2.307390    1.399205
21934   2.011412    5.079666
22814   1.580898    1.549221
23455   1.267308    15.000000
24247   0.879846    8.937476
26045   0.000228    5.542150

peak_figure

dmaziuk commented 8 years ago

Yeah, that's what I though you'd say. Processed data actually has picked peaks in peaklist.xml and also some stuff in peakrng and peaks -- if pp was run on the spectrum -- but nmrglue isn't reading those I take it. Oh well...

jjhelmus commented 8 years ago

Nope, nmrglue does not have the ability to read in these files. Only the 1r, 1i and proc/procs files.

axz91 commented 1 year ago

Here my data does not have SW_p,

what i can do, thanks.

(32, 64, 2048) {'procs': {'_coreheader': ['##TITLE=Parameter List, ParaVision 6.0.1', '##JCAMPDX=4.24', '##DATATYPE=Parameter Values', '##ORIGIN=Bruker BioSpin MRI GmbH', '##OWNER=QHeLab'], '_comments': ['$$ 2023-08-05 20:59:01.500 -0400 QHeLab@wobble', '$$ /opt/cs', '$$ process /opt/PV6.0.1/prog/mod/parx'], 'PPARMOD': 2, 'SI': 2048, 'DC': 1, 'OFFSET': 9.69458184838613, 'SF': 300.325441834007, 'YMAX_p': 32766, 'YMIN_p': 0},