Register a peak list onto a spectrum

A peak list made for one spectrum rarely lines up exactly with another, because of differences in referencing, temperature or sample conditions. This applies to a list carried over from an earlier spectrum of the same sample, and even more to one predicted from assigned chemical shifts. register! finds the single offset that best matches the list to the observed peaks, and moves the list by it. register returns a moved copy instead.

using NMRTools, CairoMakie
spec = loadnmr(example("2d_hn_ubq"))

A list from the BMRB

predicted = peaklist(loadbmrb(17769), "H", "N")
moved = register(predicted, spec)

[ Info: Downloading BMRB entry 17769
[ Info: Peak list shifted by -0.231381, 0.278222 ppm

The offset applied is logged. See Assigning ubiquitin from the BMRB for a full example.

A list from another spectrum

A list picked in one spectrum can be registered onto another, for example the next point of a titration:

series = loadnmr(example("titration2d_fir_nbox"))
pks = detectpeaks(series[1]; snr=20, sign=:positive)
moved = register(pks, series[end])
fig, ax, _ = nmrplot(series[end][7.5 .. 8.7, 116 .. 126]; title="")
nmrplot!(ax, filter(p -> 7.5 < ppm(p, 1) < 8.7 && 116 < ppm(p, 2) < 126, moved))

[ Info: Peak list shifted by 0.00417994, 0.0440127 ppm

A saved list

Load the list with loadpeaks first:

pks = loadpeaks("peaklist.csv")
register!(pks, spec)

Limiting the search

Every offset within maxshift (by default 10% of the spectral width) is tried. Give maxshift in ppm, one value per dimension, to narrow the search, or 0 to hold a dimension fixed:

register!(pks, spec; maxshift=(0.3, 2.0))   # at most ±0.3 ppm ¹H, ±2 ppm ¹⁵N
register!(pks, spec; maxshift=(0.3, 0))     # shift ¹H only

How the offset is chosen

Each candidate offset is scored by how much signal lies under the shifted peaks. Absolute intensities are used, so negative peaks count as much as positive ones, and the score grows only logarithmically with intensity, so a single strong peak cannot outweigh many weaker matches. The best offset is refined to a fraction of a point.

Only a uniform translation is applied; individual peaks are not moved. Use findpeak on each registered peak to snap it to the nearest maximum:

snapped = [findpeak(series[end], ppm(p); label=label(p)) for p in moved[1:3]]

3-element Vector{Peak{2}}:
  label  x (ppm)  y (ppm)  lwx (Hz)  lwy (Hz)  intensity
     X1  7.70461  118.476   17.1524   8.37339   224358.0
     X2  6.85562   111.62   20.6546   17.5474   136587.0
     X3  6.79706  113.074   12.8805   11.2275   143071.0

Higher dimensions

Registration works in any number of dimensions, as long as the list and spectrum have the same number:

register!(peaklist(shifts, "H", "N", "CA"), loadnmr(example("3d_hnca_ubq")))

In a pipeline

loadbmrb(17769, H1, N15) |> register(spec) |> savepeaks("bmrb.csv")