Nmrglue: An Open Source Python Package for the Analysis of Multidimensional NMR Data

Reading, Writing and Converting between Common File Formats

First and foremost, nmrglue can be used to read and write data from and to a number of common NMR file formats. At the present time, nmrglue can access and save data in formats recorded on Bruker and Agilent/Varian spectrometers, as well as files used by NMRPipe (Delaglio et al. 1995), Sparky (Goddard and Kneller 2008), SIMPSON (Bak et al. 2000), and the Rowland NMR Toolkit (Hoch and Stern 1996). Importantly, support for additional data formats can be readily implemented owing to the flexibility of the platform and is planned for future releases of nmrglue. Using the read function appropriate for the particular file format, nmrglue returns to the user the NMR data and any spectral parameters or other meta-data contained in the file. One- or multidimensional NMR data are read into memory and contained in an ndarray, a robust multidimensional array object provided by NumPy with which the scientific routines in the NumPy and SciPy libraries can subsequently interact. The dimensionality, size, direct dimension quadrature, and data type (integer or floating point values) of this object are set and appropriately converted by nmrglue, based on parameters contained in the file being read. Since NumPy provides no mechanism for the representation of hypercomplex data (States et al. 1982; Delsuc 1988), quadrature in the indirect dimensions must be manipulated as required when the data are transposed. Any meta-data and spectral parameters present in the file, such as the carrier frequencies and sweep widths, are stored in a dictionary using key-value pairs and returned to the user. Since each NMR file format contains different numbers and formatting styles of the spectral parameters, the structure and size of these dictionaries are unique to each file format. Nmrglue has a limited ability to convert between the parameter dictionaries corresponding to different file formats as discussed in more detail below.

Oftentimes it is neither necessary nor desirable to read an entire NMR data set into memory, especially when working with large multi-gigabyte 3D or 4D data sets. In many cases, only a small subset of the data is actually required at a given time, and loading the entire data set into memory can be a time and resource consuming process (or, indeed, impossible in cases where memory is limited). For instance, when visualizing a 3D spectrum as a series of 2D planes, only single 2D slices along a given frequency axis need to be accessed at any given time. Nmrglue provides methods for loading limited regions of NMR data into memory only when the data are needed to perform a calculation. This is accomplished via the read_lowmem functions, available for most file formats. The resulting array-like object which is returned behaves similarly to a NumPy ndarray object. It can be transposed, but no data are actually loaded into memory until a specific region of the data set is requested through a slicing operation, at which point the requested data are read and returned as a NumPy array.

In addition to being able to read NMR data and parameters in a number of file formats, nmrglue can write to most of the supported file formats using write functions present in the modules. In the case of large data sets, for which the low memory reading functions are employed, users can request a region of the data and write that region to disk using the write function or write the entire data set to a file trace by trace using the write_lowmem function. By using this functionality, only a single 1D data trace must be stored in memory at any given moment. Additional methods, which allow for writing to specific regions of a file, are under development and will be available in the future versions of the package.

Nmrglue can also be used to convert data between different formats. The convert module provides tools for this functionality. An example of conversion from Sparky to NMRPipe file format is given in Listing S1 available in the Supporting Information. As demonstrated by this example, the mechanism for the conversion is to read in the NMR data and parameters (line 4), create and load a conversion object (lines 7 and 8), and request the NMR data and parameters in the output format (line 9) which are then written to disk (line 12). Any manipulation of the data necessary for the conversion (sign changes, data type modifications, etc.) is performed internally so that the resulting array is correct for the requested output format. By using the read_lowmem and write_lowmem functions, large NMR data sets can be converted in this manner from one format to another in such a way that only a single trace is stored in memory at any given time.

Given that the supported file formats store differing numbers of spectral parameters, many of which are not specifically linked to a given dimension, it is not always possible to determine all of the parameters required for a particular format. During the conversion, nmrglue will attempt to find the correct parameters in the initial parameter dictionary or will fill in default values if this information is not available. Internally nmrglue uses a “universal dictionary” for these conversions, which contains the most fundamental parameters necessary to describe NMR data. Users can fill in corrected values of the various parameters by updating the resulting dictionary or by providing a custom universal dictionary when loading in the data. An example of this procedure is given in Listing S8, which will be discussed in detail in the Example Applications section. A more advanced file conversion system that is focused on ease of use is currently under development and will be included in future versions of the software.

Often locations in NMR spectra are referenced not by points, but in more convenient units such as Hz or ppm. To allow the use of these units, the file input/output modules contain functions for creating unit conversion objects, which convert to and from points in a given dimension to more common unit types including Hz, ppm, percent, microseconds, milliseconds, and seconds. The use of these objects to extract and display a spectrum in units of ppm will be shown in the Example Applications section.

Data Processing

In addition to providing the capabilities to read, write, and convert between various NMR file formats, nmrglue also contains a number of modules for processing NMR data. Since most of the data processing methods used in NMR spectroscopy are identical or highly analogous to signal processing techniques used in other scientific fields, the numerical routines in the NumPy and SciPy libraries provide efficient implementations of these algorithms. The proc_base module contains a number of common NMR processing functions including apodizations, spectral shifts, transforms, and filters. The proc_bl module provides functions for filtering, smoothing and flattening spectral baselines, and the proc_lp module has functions for linear prediction (Ni and Scheraga 1986), modeling and extrapolation using different algorithms including singular value decomposition (SVD) and total least squares (TLS). Given the widespread use and popularity of the NMRPipe software in the NMR community, the pipe_proc module has also been created to provide processing functions with names and parameters similar to their NMRPipe counterparts. In addition, the functions in this module also update the dictionary of spectral parameters as the data are processed, so that the file written after processing is nearly identical to files processed using NMRPipe. Overall, the simplicity and code readability of the nmrglue processing modules—which, together with the open distribution of the source code provide developers with example implementations of many common NMR processing methods—comes at some expense of speed. However, the reduced speed of these functions relative to highly optimized and specialized code is typically not a major hindrance, as Fourier transforms and other operations employed in NMR data processing are no longer time intensive tasks for modern computers. For example, the processing of a test 1500* × 166* 2D NMR data matrix with zero-filling to 4096* and 2048* points, respectively, which required 2.04 seconds using NMRPipe on a Dell Precision Workstation 470 with two Intel Xeon processors and 6 GB of RAM, could be achieved in 5.85 seconds using nmrglue. This difference in processing times can be partially attributed to nmrglue’s limitation of running on a single processor that is not present in NMRPipe.

Spectral Analysis

Extracting the wealth of information present in multidimensional NMR spectra can be a time and labor intensive task. In biomolecular NMR this analysis is oftentimes highly repetitive, with each residue in the protein giving rise to one or more cross-peaks which are analyzed in a similar manner. The process thus lends itself to streamlining and automation. Nmrglue contains basic analysis tools that can be used to create scripts and programs to tease out structural and dynamic information from NMR data. The peakpick module contains a function which can perform peak picking of NMR spectra of arbitrary dimensionality using several different algorithms. Estimation of peak positions and linewidths and the clustering of nearby peaks, which are all required for further analysis, can be performed within this module. One commonly used procedure involves fitting the experimental peaks to a model lineshape, and the linesh module enables such fitting of a peak or cluster of peaks to be performed in an arbitrary number of dimensions using the Levenberg-Marquardt algorithm (Marquardt 1963). Lorentzian, Gaussian, and Voigt lineshapes most commonly encountered in the analysis of NMR spectra are built into nmrglue, and users may also supply their own lineshape functions to be used in peak fitting.

In addition to the analysis capabilities described above, nmrglue can be used to rapidly develop new data analysis approaches. Such developments are aided by the underlying spectral segmentation and peak picking methods, a Levenberg-Marquardt least squares optimization algorithm which allows fitting parameters to be constrained and a spectral simulation function applicable to an arbitrary number of dimensions, as well as the extensive collections of fast and efficient scientific and computational routines available in NumPy and SciPy. Furthermore, C, C++ or Fortran code can be easily interfaced with Python using the C-API module in NumPy, F2PY, SWIG, or Cython, which allows the analysis routines written in these languages to be included in nmrglue scripts.

Interactivity

In addition to enabling the development of Python scripts that work with NMR data, nmrglue can also be used to examine NMR data in an interactive mode. This is most readily achieved using an enhanced Python shell, such as IPython, that allows data to be viewed in graphical form and includes a number of additional advanced features. Figure 1 shows a sample IPython notebook, in which nmrglue was used to examine and visualize 1D and 2D NMR data from within a web browser.

Installation and Documentation

Nmrglue was designed to provide a powerful environment in which NMR data can be converted, processed and analyzed, with particular attention paid to ease-of-use and substantial flexibility that allows users to develop and test new analysis methods. With these goals in mind, creation of quality documentation was given a high priority. Consequently, all user-facing functions and classes in the package have been thoroughly documented and are available using Python’s built-in help function or online at the nmrglue website at http://nmrglue.com. In addition, a tutorial introducing many of the features of the package, multiple examples of nmrglue scripts, as well as instructions for software installation, are also available at the website.

Visualization of 1D Time and Frequency Domain NMR Data

Visualization of NMR spectra is an essential step in the assignment of the resonances and is helpful in evaluating various processing methods. Nmrglue can be used in conjunction with matplotlib or another Python plotting package to visualize NMR data stored in any of the supported formats. As an example, Listings S2 and S3 provide Python scripts which plot the free induction decay (Listing S2) and corresponding ¹³C MAS solid-state NMR spectrum (Listing S3) recorded for a sample of amyloid fibrils formed from ¹³C,¹⁵N-enriched Y145Stop variant of the human prion protein (Helmus et al. 2008b; Helmus et al. 2010; Helmus et al. 2011). In both scripts the NMR data are read in from an NMRPipe file (line 5) and a unit conversion object is created (line 8), which is used in line 13 to provide a scale for the plot in milliseconds or ppm. Lines 11–19 plot the time and frequency domain spectra, set the scales of the axes and add appropriate labels using the matplotlib library. The last line in each script creates the output files, “fid.eps” and “spectrum.eps”, which are shown in Figures 2 and 3, respectively. Features can be added to or removed from the figures by modifying the Python scripts or using image editing software, and other graphics output formats, including png and PDF, can be created by appropriately modifying the extension of the filename on the last line.

Visualization of 2D NMR Spectra

2D NMR spectra can be plotted in a similar manner as shown in Listing S4 and accompanying Figure 4. In this script, a 2D ¹⁵N-¹³C spectrum is read from an NMRPipe file (line 5), the limits of the ppm scales for both axes are determined from unit conversion objects (lines 8–11), and a contour plot of the spectrum is created (lines 14–18). Three representative 1D ¹³C slices were also added to the plot in lines 21–27 to demonstrate the spectral resolution and signal-to-noise ratio, followed by setting the limits and labeling the frequency axes (lines 30–34).

Matplotlib and other Python plotting libraries can also be used to create plots that can be examined interactively. This is achieved by changing the last line of the Python scripts in Listings S2, S3 or S4 to fig.show(). Running the scripts with this modification will open an interactive window containing the plots in Figures 2, 3 or 4. In this environment users can zoom, pan and save the plot, and even more advanced interactive environments can be created using these and other Python tools.

Separation of Interleaved Pseudo 3D NMR Data Sets

The next two examples, given in Listings S5 and S6, show how nmrglue can be used to prepare NMR data for use by other software. In both of these examples pseudo 3D data sets, recorded on an Agilent/Varian VNMRS-500 spectrometer and consisting of series of 2D chemical shift correlation spectra acquired in an interleaved fashion as a function of a relaxation or dipolar evolution delay, are separated into their constituent 2D’s for processing with NMRPipe or another software package. These types of interleaved experiments, where the relaxation or dipolar evolution period is incremented as the innermost loop, are very common in solution and solid-state NMR since the effects of long-term instrument instabilities on the extracted relaxation rates or dipolar coupling constants are minimized for data collected in this manner (as opposed to recording complete 2D’s back-to-back with the relaxation or dipolar evolution period incremented as the outermost loop of the pseudo 3D data set). While other tools exist for separating these types of arrayed experiments, they often require the data to be collected in particular order and fail otherwise. By providing direct access to the NMR data, nmrglue can be readily used to create the required scripts that separate these types of interleaved experiments regardless of the order in which the parameters were arrayed.

In the first example (Listing S5), a relaxation experiment consisting of six 2D ¹⁵N-¹³C solid-state NMR spectra was recorded. The experiment was set up so that relaxation delay, labeled by the parameter “techo”, was the innermost arrayed parameter, followed by the phase and then the delay for the indirect ¹⁵N chemical shift dimension. The Python script separates these data set into six directories with names “techo_X.fid” where X is the value of the “techo” parameter. Each directory contains the corresponding 2D time-domain data in Agilent/Varian format for processing with NMRPipe or another software package. In the script, the interleaved data set is read using the read function of nmrglue’s varian module (line 4). For data collected with the arrayed parameter as the innermost loop, the array containing the NMR data will automatically be shaped by nmrglue so that the next-to-last dimension varies with the arrayed parameter. The size of the individual 2D spectra will be smaller than the full interleaved pseudo 3D data set, so the parameter dictionary is updated with this new size on line 7. Line 10 begins a loop over the relaxation times “techo”. Within each iteration of the loop, the appropriate region of the data set is extracted and written to the corresponding directory on line 13, which completes the script. Note, finally, that this script will work to separate other interleaved pseudo 3D NMR data sets, collected with the arrayed parameter as the innermost loop by modifying lines 10 and 11 to account for the name of the parameter being arrayed.

For interleaved experiments with a different order of the arrayed parameters the script in Listing S5 will fail, as will other programs that separate experiments collected in this manner. However, nmrglue provides the facilities for the rapid development of new scripts to deal with these cases. One such Python script, which separates a pseudo 3D data set consisting of a series of 2D z-filtered TEDOR spectra (Jaroniec et al. 2002) collected with the quadrature phase as the innermost loop, is given in Listing S6. In this case the number of applied TEDOR dipolar mixing cycles, set by the parameter “nredor”, corresponds to the next-to-last loop. This script has a similar layout to Listing S5, with a few minor modifications as required by the difference in the data ordering. Specifically, the Agilent/Varian data are read on line 5 with the keyword as_2d assigned a value of True which results in the NMR data being returned as a two-dimensional array. Lines 8–10 calculate and set the new size of the outputted data. Line 13 begins the loop over the “nredor” values. For each iteration of the loop a directory name is created (line 14) and the traces which make up the individual 2D spectra are extracted and added in the correct order to the sdata array (lines 16–18), which are then written to disk (line 19).

Processing of 2D NMR Spectra Collected with a S³E Filter

Recently it has been demonstrated that the spin-state selective excitation (S³E) approach (Meissner et al. 1997) can be used in the context of protein solid-state NMR to obtain high-resolution ¹⁵N-¹³C chemical shift correlation spectra for which the effects of ¹³C-¹³C J-coupling evolution during signal acquisition have been suppressed (Laage et al. 2009). We have recently employed this methodology for proteins containing covalently attached paramagnetic tags to rapidly collect 2D and 3D NMR spectra with high resolution and sensitivity (Nadaud et al. 2010; Nadaud et al. 2011; Sengupta et al. 2012). Processing data collected with a S³E filter involves separating the ‘A’ and ‘B’ blocks of the experiment, recombining them as the sum and difference data sets, processing these two data sets independently, and finally coadding the two resulting spectra. While this type of processing can be carried out entirely with existing software, such as NMRPipe, creating the necessary processing scripts requires the use of NMRPipe functions that are not common in routine processing scripts. Nmrglue can be used to accomplish such data processing in a very straight forward manner as illustrated by the Python script in Listing S7. This script takes an Agilent/Varian 2D data set recorded with a S³E filter, extracts the ‘A’ and ‘B’ blocks which were collected as the even and odd traces of the experiment (lines 4 and 5), and saves the sum and difference data sets as binary Agilent/Varian files (lines 6 and 7). These files can be converted and processed as conventional 2D ¹⁵N-¹³C NMR spectra using NMRPipe or other software, or processed further using nmrglue as demonstrated below.

Prior to processing the two data sets they are converted to NMRPipe files using the script in Listing S8. Here the Agilent/Varian “fid_sum” file is read (line 4), spectral parameters for the 2D spectra are recorded in the universal parameter dictionary (lines 8–14), a converter object is created and loaded with the Agilent/Varian data and the universal dictionary (lines 17 and 18) and NMRPipe formatted data are requested (lines 19). Finally, these data are written to the “test_sum.fid” file. The entire process is then repeated for the difference data set (lines 25–29), but without setting the spectral parameters which are identical to the sum data set. Listing S9 provides a script in which the direct dimensions of the sum (lines 4–12) and difference (lines 15–24) data sets are processed using nmrglue’s pipe_proc module. Note that the requisite 27.5 Hz shift of the spectra, corresponding to one half of the one-bond ¹³CO-¹³Cα J-coupling constant, is converted to points using a unit conversion object (lines 8–9). With the direct dimension processing completed, the two data sets can be coadded (lines 27) and the indirect dimension can be processed (lines 31–37) before saving the resulting spectra (line 40).

Covariance Processing of 2D NMR Spectra

In addition to conventional NMR data processing, nmrglue provides a convenient platform for the facile development of new processing methods that are not included in the common software packages. One example is covariance processing, which is an alternative to the Fourier transform for the processing of 2D NMR spectra (Brüschweiler and Zhang 2004). Recently, the Covariance NMR Toolbox, which enables covariance processing of NMR data within MATLAB or OCTAVE (Short et al. 2011), was released; to the authors’ best knowledge, this is currently the only available software package for such data processing. At the root of covariance NMR is the computation of the covariance of a matrix. The NumPy library contains an efficient algorithm to perform such a calculation and can be used for covariance processing, provided that the data are appropriately prepared before and after the calculation. The script in Listing S10 uses nmrglue to read in data from a 2D solid-state NMR ¹³C-¹³C chemical shift correlation experiment—recorded with the dipolar assisted rotational resonance (DARR) mixing scheme (Takegoshi et al. 2001) for a sample of nanotubes formed by the bolaamphiphilic self-assembly of 1,4,5,8-naphthalenetetracarboxylic acid diimide with L-lysine headgroups (Shao et al. 2010)—which have been Fourier transformed along the direct dimension (line 5). Subsequently, the script computes the covariance of the data using NumPy’s efficient algorithm (line 8), updates the necessary spectral parameters (line 11–14), and writes out the resulting covariance processed spectrum (line 17), which is shown in Figure 5.

Preparation of Strip Plots from Multiple 3D NMR Spectra

Triple resonance 3D chemical shift correlation experiments have become the standard methodology for determining the sequential backbone resonance assignments for ¹³C,¹⁵N enriched proteins. In solid-state NMR these assignments can frequently be established within Sparky or other software by using a combination of several complementary 3D ¹⁵N-¹³C-¹³C spectra (e.g., NCACX and NCOCX) (Baldus 2002). For publications and presentations, resonance assignment data are most conveniently summarized by showing aligned strips from multiple 3D data sets. This process, which can be rather tedious and time consuming, can be readily automated in nmrglue. This is illustrated in Figure 6, which shows representative strip plots from 3D CONCA, NCACX, and NCOCX spectra of the B3 immunoglobulin binding domain of protein G (GB3) in the microcrystalline solid phase (Nadaud et al. 2007). The Python script used to generate Figure 6 is presented in Listing S11. This script can be readily edited to alter the various aspects of the figure such contour colors, fonts, tick placement and labeling.

Analysis of NMR Relaxation Data

A multitude of NMR methods have been developed to extract structural and dynamic data in biological macromolecules in site-specific fashion. These methods typically rely on the quantification of an observable, such as a dipolar coupling or relaxation rate constant, in one of the indirect dimensions of a multidimensional NMR experiment. While existing software packages, such as NMRPipe, provide facilities for the rapid analysis of such data sets based on peak heights or volumes extracted using lineshape fitting routines, they often yield suboptimal results in the context of solid-state NMR spectra which exhibit irregular lineshapes and/or limited sensitivity and resolution. In such cases a simple summation of spectral intensities within a small rectangular region around each peak provides a robust approach for estimating the experimental peak volumes, and also, importantly, their uncertainties based on the intrinsic spectral signal-to-noise ratio. Our group has made extensive use of this approach to analyze multiple dipolar (Helmus et al. 2008a; Helmus et al. 2010) and relaxation (Nadaud et al. 2009; Helmus et al. 2010; Nadaud et al. 2010; Nadaud et al. 2011; Sengupta et al. 2012) trajectories in solid proteins in an automated manner.

This type of analysis in nmrglue is demonstrated here for a set of residue-specific longitudinal relaxation trajectories extracted from a series of 2D chemical shift correlation spectra recorded for the model B1 immunoglobulin binding domain of protein G (GB1). The Python script in Listing S12 reads in the integration limits and a list of spectra to analyze from text files (examples of which are available at the nmrglue website) on lines 5 and 6. An array that will hold the trajectories is created on line 9. Line 12 begins a loop in which a spectrum stored in an NMRPipe formatted file is read (line 16). An inner loop (lines 19–27) determines the correct order of the integration limits for each peak (lines 21–24) and computes the summation of all points within the spectrum that fall within these integration limits (line 27). Lines 30–36 create a series of text files ending with a “.dat” extension that contain the normalized relaxation trajectories. Trajectories which have not been normalized can be obtained by omitting line 31. Accurate locations of the integration regions around the cross-peaks are crucial for the reliable extraction of individual trajectories, especially for clusters of peaks exhibiting partial overlap. To aid in determining these locations the Python script in Listing S13 creates a series of 2D contour plots for each peak of interest (see Figure 7 for a representative example). These plots display the region defined by the current set of integration limits as well as somewhat smaller and larger integration rectangles.

The relaxation trajectories extracted using the script in Listing S12 can be further analyzed with a variety of commercial or home-built software. Here this analysis, consisting of a fit of the experimental trajectory to a decaying single exponential, is illustrated using the constrained least-squares optimization method in nmrglue. The Python script in Listing S14 defines the fitting function (lines 7–9) as well as a residuals function that calculates the difference between the experimental and simulated trajectories (lines 12–14). The latter function is required for the least squares optimization procedure on line 33. Lines 17–19 read in the relaxation times and set the initial values and constraints on the fit parameters. The bounds defined on line 19 restrict the amplitude scaling factor A to values between 0.98 and 1.02, a reasonable range for the normalized experimental relaxation trajectories. The script then creates an output file “fits.txt” for storing the fitting results (line 22–23). Line 26 begins a loop over all the peak trajectory files which have the “.dat” extension. The peak assignment is determined (line 28) and the corresponding trajectory loaded (line 32) based on the filename. Nmrglue’s least-squares optimization algorithm is then utilized to find the best-fit values for the amplitude scaling and relaxation rate parameters (lines 33 and 34), and the fitting results are stored in the output file. The quality of the fitting can be examined using the Python script in Listing S15, which creates plots of the experimental and simulated relaxation trajectories for each cross-peak. An example of such a plot, for GB1 residue D40, is shown Figure 8.

Collectively, the set of nmrglue scripts shown in Listings S12–15 enables the extraction and fitting of relaxation trajectories from a series of 2D NMR spectra in an automated and interactive manner. With relatively minor modifications these scripts can be used to analyze relaxation data from a series of 3D NMR chemical shift correlation spectra, or fit data from other types of arrayed NMR experiments by changing the fitting function in Listings S14 and S15.