peak picking in a 2D plane of a 4D spectrum

[Thomas Evangelidis]

Greetings,

nmrglue is great project! Unfortunately seems to be dead. I hope that I am proven wrong.

I am trying to extract a 2D plane of a 4D HCNH NOESY spectrum, do automatic peak picking, and plot the picked peaks as in your example. Below is my script.

import matplotlib.pyplot as plt
import nmrglue as ng

# Path to the Sparky file
file_path = "4D_HCNH_NOESY_NUS_reconstructed.ucsf"

# Read data from Sparky UCSF file
dic, data = ng.sparky.read(file_path)

# The 'dic' variable contains the metadata including the universal dictionary (udic)
# for each dimension. We need to extract this udic to create unit conversion objects.
udic = ng.sparky.guess_udic(dic, data)

# Create unit conversion objects for each dimension
# Assuming the dimensions are in the same order: 0: HC, 1: C, 2: N, 3: HN
uc = {
"hc": ng.fileiobase.uc_from_udic(udic, dim=0), # HC dimension
"c": ng.fileiobase.uc_from_udic(udic, dim=1), # C dimension
"n": ng.fileiobase.uc_from_udic(udic, dim=2), # N dimension
"hn": ng.fileiobase.uc_from_udic(udic, dim=3), # HN dimension
}

# Specific N and HN values for filtering
n_ppm = 115.438
hn_ppm = 8.671

# Convert N and HN ppm values to points; this must be wrong!
n_point = uc["n"].ppm(n_ppm)
hn_point = uc["hn"].ppm(hn_ppm)

# Extract the slice at N=120.86 ppm and HN=8.934 ppm
data_slice = data[:, :, int(n_point), int(hn_point)].real

# Define a threshold for peak detection
threshold = 8.5e4

# Detect all peaks within the extracted slice with a threshold
peaks = ng.peakpick.pick(data_slice, pthres=threshold, algorithm="thres", msep=[1, 1])

# Contour levels
clevs = [threshold * 1.4 ** i for i in range(10)]

# Plotting
fig, ax = plt.subplots(figsize=(8, 6), constrained_layout=True)

# Set the extent for the plot. Adjust the order of limits based on your preference
ext = [*uc["hc"].ppm_limits(), *uc["c"].ppm_limits()]
ax.contour(data_slice, levels=clevs, extent=ext, cmap="bone", linewidths=0.5)

# Adjust xlim and ylim based on your data range in HC and C dimensions
ax.set_xlim(uc["hc"].ppm_limits()[1], uc["hc"].ppm_limits()[0]) # Adjust according to your preference
ax.set_ylim(uc["c"].ppm_limits()[1], uc["c"].ppm_limits()[0]) # Adjust according to your preference

# Add markers for peak positions
x = uc["hc"].ppm(peaks['X_AXIS']) # Adjust axis labels if necessary
y = uc["c"].ppm(peaks['Y_AXIS']) # Adjust axis labels if necessary
ax.scatter(x, y, marker=".", s=10, color="red")

ax.set_xlabel("$^{1}$H (ppm)")
ax.set_ylabel("$^{13}$C (ppm)")

plt.show()

However, it raises the following exception. Obviously, the ng.peakpick.pick doesn't like the way the 4D spectrum was sliced. Could someone please troubleshoot me?

/home/thomas/.conda/envs/nmrglue/lib/python3.12/site-packages/nmrglue/fileio/sparky.py:517: UserWarning: Bad file size in header 14417920692 vs 0
  warn(f'Bad file size in header {seek_pos} vs {expected_seek_pos}')
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[44], line 38
     35 threshold = 8.5e4
     37 # Detect all peaks within the extracted slice with a threshold
---> 38 peaks = ng.peakpick.pick(data_slice, pthres=threshold, algorithm="thres", msep=[1, 1])
     40 # Contour levels
     41 clevs = [threshold * 1.4 ** i for i in range(10)]

File ~/.conda/envs/nmrglue/lib/python3.12/site-packages/nmrglue/analysis/peakpick.py:214, in pick(data, pthres, nthres, msep, algorithm, est_params, lineshapes, edge, diag, c_struc, c_ndil, cluster, table, axis_names)
    210 # scales = np.zeros(np.array(locations).shape,dtype=float)
    211 # amps = np.zeros(len(locations),dtype=float)
    213 for i, (l, seg_slice) in enumerate(zip(locations, seg_slices)):
--> 214     null, scales[i], amps[i] = guess_params_slice(data, l, seg_slice,
    215                                                   ls_classes)
    217 ########################################################
    218 # return locations, scales and amplitudes as requested #
    219 ########################################################
    220 if cluster:

File ~/.conda/envs/nmrglue/lib/python3.12/site-packages/nmrglue/analysis/peakpick.py:370, in guess_params_slice(data, location, seg_slice, ls_classes)
    345 """
    346 Guess the parameter of a peak in a segment.
    347
   (...)
    367
    368 """
    369 # find the rectangular region around the segment
--> 370 region = data[seg_slice]
    371 edge = [s.start for s in seg_slice]
    372 rlocation = [l - s.start for l, s in zip(location, seg_slice)]

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

[Jonathan Helmus]

What version of NumPy and nmrglue are you using to run this script?  There was a change [1] merged into nmrglue last month to fix a bug in the peak picking module that occurs when using recent versions of NumPy. You may need to update nmrglue to include this fix or downgrade NumPy to an earlier version.  

No release of nmrglue has been made with this fix but you can install the latest development version using:

python -m pip install -e nmrglue @ git+https://github.com/jjhelmus/nmrglue.git

Cheers,

 - Jonathan Hemus

[1] https://github.com/jjhelmus/nmrglue/pull/216

--
You received this message because you are subscribed to the Google Groups "nmrglue-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nmrglue-discu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nmrglue-discuss/eb8736d1-351f-4875-ac25-5a735c890667n%40googlegroups.com.

[Thomas Evangelidis]

Hi Jonathan,

Thank you for your reply. I am glad that this mailing list is still active. :)

I use numpy version 1.26.4 and nmrglue 0.10. After installing your development version, my script finished and created the attached image, which is nonsense. I am sure I am doing something wrong, e.g. in the following two lines:

n_point = uc["n"].ppm(n_ppm)
hn_point = uc["hn"].ppm(hn_ppm)

Would you be so kind as to review it, please?

Kind regards,

Thomas

You received this message because you are subscribed to a topic in the Google Groups "nmrglue-discuss" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/nmrglue-discuss/AP_s3JZn72c/unsubscribe.
To unsubscribe from this group and all its topics, send an email to nmrglue-discu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nmrglue-discuss/CAEBOqQzBXqP3wP0itfM5uK_gOBMe4mpJEAESReUzXFOx__Guuw%40mail.gmail.com.

[Kaustubh Mote]

It seems to me that the threshold is too low. Can you plot the spectrum without picking peaks with an increased threshold and ensure that the contours are present only in regions where there are peaks. Once you are sure that the threshold is optimal, you can pick the peaks again.

[Thomas Evangelidis]

Thank you Kaustubh. I will try what you suggested, but I am afraid that the issue arises at an earlier stage when selecting the correct 2D slice. Is there a built-in function in NMRglue that, given a value in ppm, can tell you the bin index it corresponds to? Or should I program this function myself, one that reads uc["n"] and uc["hn"], and based on their dimensions and minimum and maximum values, will predict the bin index of the given peak? Because I am afraid I am plotting a 2D slice without any peaks.

Thanks,

Thomas

[Jonathan Helmus]

Thomas,

 The unit conversion object can convert between PPM and the index in the array.  For your example you would use:

n_point = uc["n"](n_ppm, "ppm")

hn_point = uc["hn"](h_ppm,  "ppm")

or 

n_point = uc["n"]("115.438 ppm")
hn_point = uc["hn"]("8.671 ppm")

The ppm method on the class does the inverse conversion, from points to ppm. The tutorial has a section with examples on using these objects, https://nmrglue.readthedocs.io/en/latest/tutorial.html#the-unit-conversion-object. 

Cheers,

  - Jonathan Helmus

--
You received this message because you are subscribed to the Google Groups "nmrglue-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nmrglue-discu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nmrglue-discuss/CAACvdx24r0G6NG%3DEpKrT71D7Hd4%2Bof%2BVki%3DZiPZE6GQFeJHxUA%40mail.gmail.com.