Peaks and assignments

Peak{N}

A peak in an N-dimensional frequency-domain spectrum: where it is, how big it is, how broad it is, and what it is called.

A peak records what was picked, not what a fit produced. Positions are in ppm and linewidths are full widths at half maximum in Hz, both ordered like the dimensions of the spectrum the peak came from, so for a 2D spectrum the first coordinate is the direct dimension (x) and the second the indirect dimension (y). Linewidths are kept in Hz because, unlike ppm, they mean the same thing on every nucleus and at every field; the corresponding relaxation rate is R2 = π * linewidth.

Unknown quantities are NaN: a peak read from a Sparky list has a position and label but no intensity or linewidth.

A peak list is simply a Vector{Peak{N}}, so sort, filter, push!, deleteat! and indexing all work as usual. Peak lists are also Tables.jl tables, so DataFrame(pks) works directly.

Examples

julia> p = Peak(8.21, 120.4; label="G10")      # a 2D peak with only a position and label
Peak{2} G10
  position (ppm): 8.21, 120.4
  linewidth (Hz): NA, NA
  intensity:      NA

julia> ppm(p)
(8.21, 120.4)

julia> label(Peak(p; label="G10N"))            # a copy with a new label
"G10N"

pks = detectpeaks(spec)                   # Vector{Peak{2}}
sort(pks; by=intensity, rev=true)

See also detectpeaks, findpeak, findnearest, ppm, intensity, linewidth, loadpeaks and savepeaks.

intensity(peak)

The height of a peak at its maximum, in the units of the spectrum it was picked from, or NaN if unknown.

linewidth(peak)
linewidth(peak, n)

The full width at half maximum of a peak in Hz: a tuple ordered like the spectrum's dimensions, or the nth dimension's on its own. NaN where unknown. The corresponding transverse relaxation rate is π * linewidth(peak, n).

detectpeaks(spec; snr=10, threshold=nothing, sign=:both) -> Vector{Peak{N}}
detectpeaks(; kwargs...)

Find the peaks in a frequency-domain spectrum of any dimensionality.

The forms without the data argument return a function for use with |>.

A peak is a local maximum (or, for negative peaks, minimum) over its immediate neighbours in every dimension whose height exceeds threshold, which by default is snr times the noise level (estimated with estimatenoise if it has not been already). Positions and heights are refined to better than the digital resolution by parabolic interpolation, and linewidths are estimated from the half-height crossings in each dimension (NaN where an overlapping peak or the edge of the spectrum gets in the way). Heights are measured from the baseline, taken as the median of the spectrum, so a small offset does not produce spurious peaks; the reported intensities are heights above that baseline.

sign may be :both, :positive or :negative. Peaks are returned strongest first and are labelled X1, X2, … to match NMRAnalysis.jl.

Every dimension must be a FrequencyDimension. To pick peaks in one plane of a pseudo-2D experiment, index it first, e.g. detectpeaks(spec[:, 1]).

Examples

pks = detectpeaks(spec)                     # positive and negative peaks above 10σ
pks = detectpeaks(spec; snr=20, sign=:positive)
pks = detectpeaks(spec; threshold=1e6)      # an absolute threshold

See also findpeak, Peak.

findpeak(spec, position; label="") -> Peak
findpeak(position; label="")

The peak nearest to position (a tuple in ppm, ordered like the dimensions of spec), found by climbing from the closest data point to the top of the hill it sits on, or to the bottom of the valley if the spectrum is negative there. The position, height and linewidths are refined just as in detectpeaks.

The forms without the data argument return a function for use with |>.

This is what snapping a mouse click onto a peak means.

Examples

p = findpeak(spec, (8.2, 120.0); label="G10")
p = findpeak(spec1d, 4.7)

findnearest(peaks, position; scale=1) -> Union{Int,Nothing}

The index of the peak in peaks nearest to position (a tuple in ppm, or a number for 1D peaks), or nothing if the list is empty. This is what selecting a peak with a mouse click means; the index can be passed straight to deleteat! or used to replace the peak.

scale weights the dimensions before distances are compared, as a tuple or a single number, since a ppm of ¹⁵N is not a ppm of ¹H. A typical choice for an ¹H,¹⁵N spectrum is scale=(1, 0.2); a GUI can pass the inverse of the visible range in each dimension so that distances are measured as they appear on screen.

Examples

i = findnearest(pks, (8.2, 120.0); scale=(1, 0.2))
deleteat!(pks, i)

register!(peaks, spec; maxshift=nothing) -> peaks
register(peaks, spec; maxshift=nothing) -> Vector{Peak{N}}
register!(spec; maxshift=nothing)
register(spec; maxshift=nothing)

Translate a peak list onto a spectrum: find the single offset, applied to every peak, that best matches the list to the signals observed in spec, and move the peaks by it. register returns a moved copy and leaves peaks unchanged. The offset found is reported with @info.

The forms without the data argument return a function for use with |>.

This corrects for referencing differences, e.g. between a list predicted from BMRB shifts (see loadbmrb) and a spectrum recorded under different conditions. It does not move peaks individually.

The match is scored as $\sum_i \log(1 + |S(\delta_i + \Delta)| / \sigma)$: the absolute intensity of the spectrum $S$ at each shifted peak position, so that negative peaks count as much as positive ones, in units of the noise $\sigma$. The logarithm gives diminishing returns for very intense signals, so that many peaks landing on signal outweigh one peak landing on a huge one. To tolerate small disagreements between individual peaks, each point of the spectrum takes the largest absolute intensity among its immediate neighbours.

Every whole-point offset within maxshift is tried, so the match is global rather than the nearest local optimum, and the best is then refined to a fraction of a point. maxshift is the largest offset considered in each dimension in ppm, as a tuple or a single number, and defaults to 10% of the spectral width. Set it to 0 in a dimension to hold that dimension fixed.

Examples

pks = loadbmrb(15000, H1, N15)
register!(pks, spec)                        # shift the BMRB list onto the HSQC
register!(pks, spec; maxshift=(0.3, 2.0))   # search at most ±0.3 ppm ¹H, ±2 ppm ¹⁵N
register!(pks, spec; maxshift=(0.3, 0))     # shift ¹H only

register(pks, spec; maxshift=nothing) -> Vector{Peak}
register(spec; maxshift=nothing)

A copy of the peak list pks moved by the offset that best matches it to spec; pks is not changed. See register!, which moves the list in place. The second form returns a function for use with |>.

Examples

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

loadpeaks(filename) -> Vector{Peak{N}}
loadpeaks(filename, spec) -> Vector{Peak{N}}

Read a peak list: either a peaklist.csv as written by savepeaks and by NMRAnalysis.jl, a hand-made list of label x y lines, or a Sparky .list file.

CSV files are matched by column name, so column order and extra columns do not matter. Linewidths and intensities are read where present and are otherwise NaN. A list in which peaks were tracked plane by plane (with a plane column filled in) is read as one position per peak, taking each peak's first row.

Sparky lists name their dimensions w1, w2, … with the indirect dimension first, which is the reverse of NMRTools' order. Passing the spectrum spec lets the axis order be checked against where the shifts actually fall, so a list written the other way round is still read correctly.

Examples

pks = loadpeaks("out/peaklist.csv")
pks = loadpeaks("hsqc.list", spec)

savepeaks(filename, peaks)
savepeaks(filename)

Write a peak list. A filename ending in .list is written in Sparky format; anything else is written as a peaklist.csv that NMRAnalysis.jl reads directly, with positions in ppm, linewidths in Hz and intensities at full precision.

The forms without the data argument return a function for use with |>.

Peak labels identify peaks in NMRAnalysis, so they must be unique and may not contain commas. Unlabelled peaks are written as X1, X2, … by their position in the list.

Examples

pks = detectpeaks(spec)
savepeaks("peaklist.csv", pks)
savepeaks("hsqc.list", pks)

ChemicalShift

One assigned chemical shift: an atom of a residue, its nucleus, and its shift.

A shift list is a Vector{ChemicalShift}, read with loadbmrb and turned into the peaks expected in a given experiment with peaklist. Use filter to select residues before building peaks, e.g. filter(s -> s.residuename in ("ILE", "LEU", "VAL"), shifts).

loadbmrb(id::Integer; list=nothing) -> Vector{ChemicalShift}
loadbmrb(filename; list=nothing) -> Vector{ChemicalShift}
loadbmrb(id_or_filename, dims...; list=nothing) -> Vector{Peak{N}}

Download the assigned chemical shifts of BMRB entry id, or read them from a local NMR-STAR (v3) file.

Given dimensions as well, the shifts are assembled straight into the peak list of an experiment, exactly as peaklist does: H1, N15 for a ¹⁵N HSQC, H1, C13 for a ¹³C HSQC, or atom names such as "H", "N", "CO-1" for an HNCO.

An entry may hold several shift lists, e.g. under different sample conditions. The first is used unless list selects another by its ID. Shifts of nuclei NMRTools does not know about (such as ¹⁷O) are skipped.

Examples

loadbmrb(15000, H1, N15)                 # ¹⁵N HSQC peak list, side chains included
loadbmrb(15000, H1, C13)                 # ¹³C HSQC peak list
loadbmrb(15000, "H", "N")                # backbone amides only
loadbmrb(15000, "CO-1", "N")             # CON

shifts = loadbmrb(15000)                 # the shifts themselves, to filter or reuse
shifts = loadbmrb("bmr15000_3.str")      # a local NMR-STAR file

peaklist(shifts, atoms...) -> Vector{Peak{N}}
peaklist(shifts, H1, X) -> Vector{Peak{2}}
peaklist(shifts, X) -> Vector{Peak{1}}

The peaks expected in an N-dimensional experiment, assembled from a list of assigned chemical shifts (see loadbmrb). Dimensions are given in the order of the spectrum, direct dimension first.

By atom name. Each argument names an atom, optionally with a residue offset relative to the residue that the peak belongs to: "CA-1" is the Cα of the preceding residue. A peak is produced for each residue where every atom has an assigned shift. "CO" and "HN" are accepted for the BMRB names "C" and "H".

By nucleus. Passing H1 and another nucleus gives one peak for every directly bonded pair of atoms with assigned shifts: every N–H for an ¹⁵N HSQC, including side-chain NH and NH₂ groups, or every C–H for a ¹³C HSQC. Methyl and methylene protons with identical shifts give a single peak.

A 1D list. Passing a single nucleus gives one peak for every shift of that nucleus, e.g. the expected lines of a ¹³C or ¹⁹F 1D spectrum.

Peaks are labelled by residue and atom, e.g. L23H-N, K22C-L23N or L23HD1-CD1. Positions come from the shift list, while intensities and linewidths are NaN.

Examples

shifts = loadbmrb(15000)              # download BMRB entry 15000

peaklist(shifts, "H", "N")              # backbone amide ¹H,¹⁵N HSQC
peaklist(shifts, H1, N15)               # all NH groups, including side chains
peaklist(shifts, H1, C13)               # ¹H,¹³C HSQC
peaklist(shifts, C13)                   # every ¹³C shift, as a 1D list
peaklist(shifts, "C-1", "N")            # CON: C′ of residue i-1 with N of residue i
peaklist(shifts, "H", "N", "C-1")       # HNCO
peaklist(shifts, "H", "N", "CA")        # the i peaks of an HNCA
peaklist(shifts, "H", "N", "CA-1")      # and the i-1 peaks

ilv = filter(s -> s.residuename in ("ILE", "LEU", "VAL"), shifts)
peaklist(ilv, H1, C13)                  # ILV methyl (and other CH) peaks

pcls(A, y)

Compute the phase-constrained least squares solution:

\[y = A x e^{i\phi}\]

following the algorithm of Bydder (2010) Lin Alg & Apps.

Returns the tuple (x, ϕ), containing the component amplitudes and global phase.

If passed a matrix Y, the function will return a matrix of component amplitudes and a list of global phases corresponding to each row.

Arguments

• A: (m,n) complex matrix with component spectra
• y: (m,) complex vector containing the observed spectrum