API: NMRAnalysis

NMRAnalysis

Interactive analysis of NMR relaxation, diffusion, exchange and 2D peak series. Set your working directory to the folder holding your experiments, call one of the routines below, and use ?name for help on any of them. Bruker experiment numbers work wherever a path does, and any parameter that is not given is read from the pulse-sequence annotations or acquisition parameters where possible, and otherwise asked for.

Automatic analysis

• analyse(filename), analyse([filenames...]): choose a routine from the annotations

1D experiments

Each opens a window to pick integration and noise regions, fits live, and returns the results when the window closes.

• relaxation1d(filename): R₁, R₂ and inversion recovery
• diffusion1d(filename): diffusion coefficient and hydrodynamic radius
• tract(trosy, antitrosy): rotational correlation time from TRACT
• calibration1d(filename): pulse-length calibration from a nutation experiment
• kinetics1d(filename, times): intensities of named regions against time
• peaktrack1d(filenames): positions of peaks and multiplets across a series of spectra
• r1rho([directory]): R₁ρ relaxation dispersion
• exchange1d([filenames]): CEST and R₁ρ fitted by Bloch-McConnell simulation

2D experiments

Each opens a window to place and fit peaks, and returns the analysis when the window closes: results(expt) gives one row per peak and planeresults(expt) one per peak per plane. All take skipplanes to leave planes out of the fitting, and peaklist to load a saved peak list as the window opens.

• fit2d(files): positions, linewidths and amplitudes, with no model
• relaxation2d(files): R₁ or R₂ from an exponential decay
• recovery2d(files): inversion or saturation recovery
• modelfit2d(files, x, equation, parameters): a model of your own
• hetnoe2d(reference, saturated): heteronuclear NOE
• ccr2d(decay, buildup, T): cross-correlated relaxation rates
• methylccr2d(buildup, decay, T): methyl S²τc
• cest2d(file): CEST profiles
• cpmg2d(file): CPMG relaxation dispersion
• pre2d(diamagnetic, paramagnetic; type): paramagnetic relaxation enhancement Γ₂
• lpre2d(diamagnetic, paramagnetic; type): Γ₁ and Γ₂ from a series of recycle delays
• peaktrack2d(files): peaks that move from plane to plane
• titration2d(files): binding isotherms and a global Kd
• rdc2d(; isotropic, aligned): one-bond couplings and RDCs
• summaryplot(expt): a fitted parameter against residue number

Other tools

• viscosity(solvent, T): solvent viscosity
• B1Calibration: B₁ calibrations, from calibration1d or a reference pulse

DEFAULTINHOMOGENEITY

Fractional B₁ inhomogeneity assumed when nothing has been measured: 5%, typical of a standard probe. Only a nutation calibration can do better than this guess.

AnalysisCall

The call that produced an analysis, recorded by the entry point so summary.txt can print a line that repeats it. reproducible is false where one of the positional arguments has no Julia literal - a spectrum passed as an NMRData rather than as a path - in which case no such line is printed at all rather than a misleading one.

B1Calibration(power, ν1; nuc=nothing, inhomogeneity=DEFAULTINHOMOGENEITY, source="")
B1Calibration(spec, nuc; inhomogeneity=DEFAULTINHOMOGENEITY)

The B₁ field a power level produces, and how widely that field is distributed across the sample.

The field follows the ideal square-root power law,

ν₁(p) = ν₁ref · 10^(−L·(dB(p) − dB(pref))/20),

in which the linearity L is 1 for a perfectly linear amplifier. One measurement fixes ν₁ref and assumes L = 1; two or more fit both, so L is a direct diagnostic of amplifier compression, and departures from 1 are what a multi-power calibration is for.

The second form is the nominal calibration read from a spectrum's own reference pulse (p1/pl1, via referencepulse): one point, L = 1, and an assumed inhomogeneity. It reproduces hz(p, refpower, refpulse, 90) exactly, so an analysis given no calibration behaves as it always did.

Use hz(power, cal) to get a field from a power and Power(ν1, cal) to invert it.

Arguments

• power: the power levels calibrated at, as NMRTools Powers.
• ν1: the field measured at each, in Hz, with uncertainties.
• nuc: the nucleus (channel) this calibration applies to, for checking it is applied to the right one.
• inhomogeneity: fractional width of the B₁ distribution (0.05 for 5%).
• source: where the calibration came from, for provenance in saved results.

Example

cal = B1Calibration([Power(18.0, :dB), Power(12.0, :dB)], [1250.0 ± 5.0, 2480.0 ± 9.0])
hz(Power(15.0, :dB), cal)    # field at an intermediate power, Hz
Power(2000.0, cal)           # power needed for a 2 kHz field

B1Distribution(σ)
B1Distribution(calibration)
B1Distribution(scaling, weight)

The distribution of B₁ across the sample, as the quadrature a simulation actually needs: a set of field scalings relative to nominal and the weight each carries.

σ is the fractional width of a Gaussian distribution, sampled at the three Gauss-Hermite nodes 1 ± √3σ with weights 1/6, 2/3, 1/6. σ = 0 gives the single nominal field, so a simulation with no inhomogeneity costs exactly what it did before.

Holding the quadrature rather than the distribution is what keeps this general: another family, or a B₁ profile measured on the probe, is a different set of scalings and weights and nothing downstream changes. The third form takes them directly.

Example

d = B1Distribution(0.05)          # 91.3%, 100% and 108.7% of nominal
b1average(ν -> ν^2, d, 1000.0)    # ⟨ν²⟩ over the distribution

Power(ν1, calibration) -> Power

The power level giving a field strength of ν1 Hz, inverting hz(power, cal). This is what turns a list of wanted spin-lock strengths into a list of powers to set.

acqusvalue(spec, keys...) -> value or nothing

Acquisition-parameter lookup returning nothing rather than throwing when the parameter is absent, or present but empty. The counterpart of annotation.

analysiscall(func, args...; kwargs...) -> AnalysisCall

Record a call. Keyword arguments whose values have no literal (a SeriesModel object, say) are left out, since the point of the record is that it can be pasted back into Julia.

Annotation lookup returning nothing rather than throwing when absent.

ask(label, default; unit="", note="", type=Float64, prompt=true) -> value
ask(label; unit="", note="", type=Float64, prompt=true) -> value

Ask the user for a parameter value, on one line.

Given a default, the value is offered for confirmation and pressing enter accepts it; without one (or with default === nothing) a value must be typed. unit and note are shown alongside to say what is wanted and where the default came from. Re-asks until the answer parses as type.

With prompt=false nothing is printed: the default is returned, or an ArgumentError thrown if there is none.

askchoice(question, options; default=1, prompt=true) -> Int

Ask the user to pick one of options, returning its index. With prompt=false, or if the menu is cancelled, default is returned.

askpath(label) -> String

Ask for the path of an experiment folder, re-asking until it exists. Used by the entry points that can be called with no arguments at all.

askvector(label, n; unit="", prompt=true) -> Vector{Float64}

Ask for the n values of an arrayed acquisition variable - a list of relaxation delays, say. They can be typed directly, separated by spaces or commas, or the name given of a file holding one per line, which is how Bruker writes a vdlist.

b1average(f, distribution, ν1)

Average f(ν) over the B₁ distribution of a nominal field ν1.

This is the right average for anything linear in the observable, which is to say for a directly measured intensity: a CEST profile is the weighted sum of the profiles the separate parts of the sample give. For a quantity that reaches the data through a fit, see b1rate.

b1rate(f, distribution, ν1, T)

Apparent relaxation rate over a B₁ distribution, where f(ν) is the rate at field ν and T the evolution time the rate was measured over.

A spread of B₁ gives a spread of rates, so the decay is a sum of exponentials, and what a monoexponential fit to the data reports is

R = −ln(Σᵢ wᵢ exp(−Rᵢ T)) / T

rather than the weighted mean of the rates ⟨R⟩. The two agree as T → 0 and diverge by about ½T·var(R) beyond that, which at the low spin-lock strengths of an R₁ρ dispersion is comparable with the uncertainty on the measurement. T is the time the experiment actually sampled, so the simulated quantity is the one the data reduction produced.

backupfile(filepath)

Rename an existing file to <name>.bak, so a save never silently destroys the last one.

backupfolder(folder) -> String

Move an existing output folder aside to <folder>_previous and return the (now empty) path, so that a save never silently destroys the last one and never leaves stale files from it behind. An earlier _previous is replaced.

Folder-level rather than file-level, because the stale files are the problem: a region or a peak deleted between one save and the next would otherwise leave its plot and its data in place, looking like part of the current result.

calibrationfit(dB, ν1, dBref) -> (ν1ref, linearity)

Weighted straight-line fit of log₁₀ν₁ against power, returning the field at dBref and the linearity (1 for an ideal amplifier).

The fit is linear because the ideal power law is a straight line in these coordinates, so there is nothing to iterate and no initial guess to get wrong. Written as the five weighted sums of the standard treatment (Numerical Recipes §15.2), with the abscissa centred on dBref so that the intercept is log₁₀ν₁ref and its variance comes straight from the covariance matrix. Uncertainties on ν₁ transform as σ(log₁₀ν) = σ(ν)/(ν ln10).

callstring(call, extra=[]) -> String or nothing

call as Julia code to paste back, with the keyword arguments in extra (names and already-rendered values) appended, one keyword to a line. nothing where the call has no literal form (see AnalysisCall).

callvalue(x) -> String or nothing

x as a Julia literal, or nothing where it has none. Vectors of numbers are written out in full: a gradient ramp or a delay list is exactly what has to be repeated, so abbreviating it would defeat the purpose.

csvcolumn(name, unit) -> String

A CSV column header: the ASCII parameter or coordinate name, with its unit in parentheses where it has one ("R (s-1)", "A"). Units are ASCII so that a file survives being opened on any machine, and a reader locates a column by stripping the parenthesised part.

csvcolumns(name, unit) -> (String, String)

Value and uncertainty column headers for one parameter. The uncertainty repeats the unit, so the pair is symmetrical for anything reading them.

csvvalue(x) -> String

One cell, at full precision - these are machine files, and rounding is irreversible. "NA" for anything absent or not applicable, which is distinct from a blank key (a blank key means the row applies to every value of it).

The field strength measured at each power level, with its uncertainty.

inhomogeneity(calibration) -> Float64

Fractional width of the B₁ distribution (0.05 for 5%).

isexpt(directory) -> Bool

Check if the given directory contains an NMR experiment by verifying the existence of a 'pdata/1/title' file, which is typical for Bruker NMR data structure.

Return true if the directory has the expected title file structure, false otherwise.

linearity(calibration)

Fitted power-law exponent relative to the ideal, so 1 means the amplifier follows ν₁ ∝ √W. Below 1 the field falls short of the ideal law as the power is raised, the usual signature of amplifier compression. Exactly 1, without uncertainty, where a single measurement left nothing to fit.

A name as it reads mid-sentence: a lower-case first letter, unless it begins an acronym.

parsevector(response) -> Vector{Float64} or nothing

Read a whitespace- or comma-separated list of numbers, or the contents of a file named by response. Returns nothing if the response is neither.

pick_expt_folders(directory_path; experiment_type_filter="") -> Vector{String}

The power levels a calibration was measured at.

The calibration's reference power.

safename(label) -> String

A label made safe to use as a filename, so that an entity the user named anything at all still writes to <name>.csv beside its plot without escaping its folder.

sanitizelabel(label) -> String

Strip commas from a user-entered region or peak label. Labels are written unescaped into CSV cells (see writetable), so a comma in one would corrupt the row.

select_expts(directory_path=""; experiment_type_filter="")

Interactively select NMR experiment folders from a given directory path.

Arguments

• directory_path::String: Path to start folder selection from. If empty, uses last selected folder or current directory.
• experiment_type_filter: Optional filter to apply for experiment types.

Returns

• Vector{String}: Array of selected experiment folder paths. Returns empty array if no selection is made.

shortpath(path) -> String

The last two components of path, with a Bruker pdata/<n> tail dropped first: /Users/chris/NMR/crick-701/sophia_990_260823/10/pdata/1 → sophia_990_260823/10.

What identifies a spectrum to the person who recorded it is the dataset folder and the experiment number, not the fifty characters above them, and a full path overflows a GUI panel and pads every line of a summary. The call that produced an analysis is still recorded in full (see AnalysisCall), so nothing needed to repeat it is lost.

Anything that is not a path is returned as it stands.

Where the calibration came from, for provenance in saved results.

stderrors(fit) -> Vector{Float64}

Standard errors of a least-squares fit, or NaNs where the covariance matrix cannot be inverted.

LsqFit.stderror inverts JᵀJ, which is singular whenever the data leave a parameter unconstrained: an integration region dragged onto the baseline, an amplitude that has gone to zero, two coincident points on the fit axis. Throwing is the right answer in a script and the wrong one in an interactive fit, where the same exception arrives inside a Makie callback and takes the window's refit down with it. A NaN uncertainty instead propagates into the results and the plots, where it is visible and harmless.

The warning is capped: a fit re-run on every drag of a region would otherwise fill the console with the same line.

viscosity(solvent, T)

Calculate the viscosity of a solvent (in mPa s) at a given temperature.

Arguments

• solvent: Should be either :h2o or :d2o.
• T: The temperature (in K) at which the viscosity is calculated.

Notes

The calculation is based on the formula proposed by Cho et al in the paper "J Phys Chem B (1999) 103 1991-1994".

windowname(analysis) -> String

The title of an analysis window, naming the package and the analysis: windowname("1D peak tracking") == "NMRAnalysis.jl - 1D peak tracking".

writetable(filepath, comments, header, rows) -> String

Write one CSV: each line of comments as a # comment (the experiment description), then the header row, then rows. Any existing file is backed up first. Returns the path.

Field strength (Hz) at the calibration's reference power, with its uncertainty.

hz(power, calibration) -> Float64

The B₁ field strength (Hz) a power level produces, per calibration.

The calibration's own uncertainty is not propagated: it is far smaller than the B₁ inhomogeneity that B1Distribution accounts for, and the analyses downstream treat the nominal field as exact.

The nucleus (channel) a calibration was measured on, or nothing if not recorded.