Nowadays,NMR (Nuclear Magnetic Resonance) spectra as a powerful tool of metabolomics,which aims at studying the metabolites in organisms,tissues or biofluids,is widely used in many research fields such as pharmacology,toxicology,nutriology and disease diagnosis^{[1, 2, 3, 4, 5]}.

The advantages on little,or without,sample preparation,and high-throughput analysis as well,make proton NMR spectroscopy (¹H NMR) standing out in all analytical strategies^{[6, 7]}. Ideally,the NMR spectra should have consistent chemical shifts and stable baseline. However,it always drifts due to uncontrolled factors,such as variation of pH value and temperature,intermolecular interactions in the samples,instabilities of the NMR apparatus (temperature and vibration of the probe,etc.),and other insuperable physicochemical interactions^{[6, 8]}. Thus the NMR spectra post processing which mainly includes noise reduction,spectrum alignment (except for the microenvironment investigation) and peak extraction is an essential procedure for spectroscopy analysis^[9]. Due to the huge data sets,especially with a high-throughput strategy,of ¹H NMR spectra,the automatic batch spectrum alignment is a great challenge. Although a metabolic expert could carry out this tedious task manually,it has serious shortcomings because it is time consuming,highly subjective and often empirical on several fields,such as NMR spectroscopy,biochemistry and biology. The bias introduced by different users may affect the final experiment results. Hitherto,many methods have been employed in spectrum alignment,such as COW (correlation optimized warping),MSPA (multi-scale peak alignment) and GPA (Gaussian smoothing peak alignment). However,these algorithms have their limitations,and the utilizations are usually very complex,not suitable for common researchers^{[10, 11, 12]}.

Furthermore,the high-throughput data analysis has become an inevitable trend in the development of modern technology. Thus the variable extraction method of a parallel NMR spectra data is a crucial factor influencing the efficiency of the scientific research. Currently,there are two variable extraction methods commonly used in the high-throughput NMR data analysis process: subsection integration method and target analysis method^{[13, 14]}. In the subsection integration method,the whole NMR spectrum is divided into a set of variables which are also called ‘bins’ or ‘buckets’ corresponding to the spectra regions. The bins have fixed chemical shift range (such as  0.002 or 0.005),and the area of the bins are calculated for the variable selection. It is simple and fast for most NMR spectra. Furthermore,it only needs a little NMR experience for the investigators. However,the results usually lack biological significance and are easily influenced by pH,ionic concentration,peak overlap or other factors. In addition,it is based on the sacrifice of data resolution. To overcome these drawbacks,the target analysis method was proposed. This approach is based on matching the peaks with the spectrum of the reference compounds from databases or libraries which are defined as a set of ¹H NMR spectra,with one file per authentic compound. The variables have strong chemical meaning,since the resonance peaks are the finger prints of the chemical compounds. However,good prior knowledge with NMR spectra is required for the researcher,and databases or libraries for target chemical compounds are needed. Furthermore,this complex data analysis approach is very time consuming and highly demanding on researchers. To fulfill the aims of lower dependence on the prior knowledge concerning the composition of complex mixture and the reduction of time,the peak extraction approach should be a great choice after the spectrum is aligned accurately. Therefore,the present work proposes an efficient NMR spectra processing method for spectrum alignment and peak extraction.

Rat plasma and urine samples are widely used in animal metabolism research. Consequently,the efficiency of our approach was tested by applying to the NMR spectra of rat plasma and urine. The implementation Matlab code and the NMR spectra data of rat plasma and urine samples in this article are available upon request.

The following sections include two steps of the method,spectrum alignment and peak extraction. The steps of the spectrum alignment and peak extraction are summarized in the flowchart (Fig. 1). In order to show it clearly,the detailed theories are described as follows.

Fig. 1 The flowchart of NMR spectrum alignment and peak extraction

Figure options

Spectra preprocessing usually contains baseline correction and noise reduction,as well as artifacts removal in each spectrum. For baseline correction,the most commonly used method is MLR (multiple linear regression) provided by the commercial software,for example,Topspin,the software provided by Bruker. This step is always carried out by researchers manually,on each spectrum separately and independently from other spectra,which is efficient for a few spectra. In the current study,the spectra of both the urine and plasma samples were preprocessed in Topspin manually. Then the raw NMR data were analyzed in the following steps.

Ideally,a one-to-one correspondence for each of NMR spectra across all samples should be established. However,it is generally not the case because of uncontrolled changes in chemical shifts of NMR peaks due to slight differences in the detection environment (samples and the NMR apparatus). Thus,the alignment step is essential^[15]. Our processing method of spectrum alignment also offers a batch analysis for a set of NMR spectra. At first,a pre-alignment was completed according to the chemical shift of the internal standard substance or a characteristic peak of one target compound (such as NAA in the brain extracts). Then,a random spectrum of the spectra was chosen as the reference spectrum. It was divided into several consecutive searching intervals. In order to decrease the influence of the shift of the spectra during the alignment,the start and end points of the searching intervals were set to the flat region (very few NMR signals). The spectrum alignment was implemented in every interval.

In each interval,the highest peak value (Y_max) and the corresponding chemical shift of this peak (δ_ref) were found at first. This peak value was set to the shift target of the other spectrum. Due to the different chemical environment and the chemical components in various samples,this peak should not be the highest value in the other spectra. Thus a search of rage (Δδ) of the chemical shift was set at first. According to the experience,the default value of Δδ' was 0.2. All the peaks were found in the range of $[{\delta _{ref}} - \Delta \delta ',\;\,{\delta _{ref}} + \Delta \delta ']$ of the spectra.

For searching the peak values,the following algorithm was employed. Between the searching region $[{\delta _{ref}} - \Delta \delta ',\;\,{\delta _{ref}} + \Delta \delta ']$ in the other spectra,if Y_m > Y_m-1 and Y_m > Y_m+1 $({\delta _{m - 1}} \ge \;{\delta _{ref}} - \Delta \delta '$ and ${\delta _{m - 1}} \le \;{\delta _{ref}} + \Delta \delta ')$,Y_m was the peak value,and the corresponding chemical shift was δ_m. At the end,the number of peaks (pn) was found; all the peaks with the value ${Y_{p1}},\;{Y_{p2}},\; \cdots ,\;{Y_{pn}}$ and their reference chemical shift ${\delta _{p1}},\;{\delta _{p2}},\; \cdots ,\;{\delta _{pn}}$ were sought out.

Then,these peaks were moved to δ_ref in succession,and the correlation coefficient R_px (px =p1 - pn) between the spectrum and the reference spectrum in the search region $[{\delta _{ref}} - \Delta \delta ',\;\,{\delta _{ref}} + \Delta \delta ']$ was calculated. The correlation coefficient was described in the following section. Finally,the maximum of the correlation coefficient was chosen,R_max,Y_px was the corresponding peak value,and δ_px was the chemical shift for this peak. The step length of the shift was calculated from the difference δ_px and δ_ref. Tiny gap was generated after the shift of the spectrum,and it would be filled with the neighboring points. The same method was also applied to the other regions.

Here,the correlation coefficient R^[16] was introduced,which was widely used in the evaluation of the efficiency of spectrum alignment. The method of calculating R is shown in the following:

Peak extraction method provides a fast way of reducing the complexity of NMRS analysis task. It stands out compared to other two kinds of popular methods of analysis (segmented-direct-integration or target analysis) for enhancing sensitivity required in biologically meaningful areas and does not require the user to have strong NMR background^{[3, 17]}.

The peak extraction approach was applied to the whole spectrum. However,the former peak searching method would involve noise values. Thus a new peak-finding approach was introduced in this application. At first,the maximum value (Y_max1) was found in the whole spectrum,and δ_max1 was the corresponding chemical shift of this peak. Then the whole region of this peak was found according to the following method: If ${Y_{\max 1 - i}} > {Y_{\max 1 - (i - 1)}},\;{\delta _{\max 1 - (i - 1)}}$ was the chemical shift of the start point of this peak,and if ${Y_{\max 1 + j}} > {Y_{\max 1 + (j - 1)}},\;{\delta _{\max 1 + (j - 1)}}$ was the chemical shift of the end point of this peak $(i,j = 1,\cdots ,n)$,then the first peak was obtained between the start point and the end points which was deducted from the whole spectrum. The maximum value of the new spectrum and the new peak were found based on the former description.

In order to avoid the influence of the noise,a noise value (Y_noise) should be defined at first. If Y_max was smaller than Y_noise,the peak extraction method was completed. However,some meaningless peaks would be included with this whole spectrum searching method,which had only a few points. Thus these tiny peaks were merged with the neighboring ones. Using this method all the peaks should be interested,and stand for one or a few chemicals. Fewer variables were selected at the end,and it could be better adapted for many statistical analyses.

In order to validate the efficiency of the method,this new method was applied to align the NMR spectra of the wine samples,which were always used to evaluate the efficiency of the alignment^{[3, 4, 20]}. Furthermore,the NMR spectra of the most popular bio-fluids samples (urine and plasma) were employed. During the metabolomics analysis,the spectrum of urine is the most difficult for peak alignment,as the pH value varies in various samples,while the spectrum of plasma is much easier. In order to save the space,the urine spectra (of 60 samples) were employed for peak alignment,and the blood spectra (of 23 samples) were used in peak extraction.

During the preparation of the samples,100 μL of pure urine/plasma mixture was transferred to a 5 mm NMR tube containing 60 μL D2O and 440 μL phosphate buffer (pH=7.2,5 mmol/L,formate as the internal standard compound). The spectra were acquired on a Bruker Avance 600 spectrometer (Bruker Biospin) operating at 600.13 MHz.

Owing to pH variation and complex components,it is very difficult to perform spectrum alignment for urine samples. Our spectrum alignment theory and algorithm lead to a remarkable improvement as shown in the results in Fig. 2 and Fig. 3.

Fig. 2 Overview of the results of spectrum alignment in a typical local region for 60 rat urine samples: (a) and (b) show the result in the spectral region contains three reference peaks (from δ3.40 to 3.46). The x-axis is the chemical shift; the y-axis is the peak intensity

Figure options

Fig. 3 Overview of the whole 60 rat urine NMR spectra: (a) pre-alignment according to Formate δ-8.46; (b) the reconstruction after subsection alignment. The x-axis is the chemical shift; the y-axis is the peak intensity. The spectra cover 3.8 spectral region (from δ4.60 to 0.80)

Figure options

The correlation coefficient R only reached 0.634 4 after the pre-alignment based on the internal standard compound: Formate δ-8.46 [Fig. 3(a)]. According to the former description,one spectrum was selected as the reference spectrum. Then this spectrum was divided into several regions (Table 1). The peak alignment was completed in each region independently. If the default value Δδ' was not perfect for the alignment,it could be adjusted manually (Table 2). In order to demonstrate the efficiency,a typical region of the spectrum (from  3.384 2 to 3.538 3) was collected in Fig. 2 (zoom in from δ3.40 to 3.46). After the reconstruction of these 22 regions,the correlation coefficients R of these 60 spectra covering a 3.8 spectral region (from δ4.60 to 0.80) have been successfully achieved to 0.904 3 [Fig. 3(b)].

To generate relevant peaks with a meaningful characteristic corresponding to chemical compound or the overlapping of several compounds,the peak extraction processing method is extremely important^{[14, 18]}. The original data always include a great deal of meaningless information (noise). Thus the crux of this step is to remove this information or combine it to the meaningful peaks. In the rat plasma spectrum,a region of 4.0 spectral region (from δ4.50 to 0.50) was selected. Using the upper peak extracting method,there were 1 162 peaks or noise peaks (mostly noise) selected. This should cause a tremendous burden for further data analysis. However,a noise value was set to 15 000,the total number of variables was dramatically decreased to 395 (Fig. 4). Here the noise value was set manually according to the characteristic of the spectra. The whole peak dataset contains the signal intensities,the reference chemical shift and the peak region. All these variables were exported to Excel,and could be converted to other formats,for adapting all kinds of further statistical analysis^[19].

Fig. 4 Results of the peak extraction of the rat plasma NMR spectra. (a) a 4.0 spectral region (from δ4.50 to 0.50); (b) the enlargement of the peaks information mainly centered at 1.00 in (a). The x-axis is the chemical shift; the y-axis is the peak intensity. The circle means the acme of a peak. And the dotted lines represent the peak regions

Figure options

The spectrum alignment method was also applied to the spectra of the wine samples which were used to evaluate the efficiency of the spectrum alignment before^{[3, 4, 20]}. According to the literature retrieval,the best correlation coefficients for these spectra were 0.98 with the GPA (Gaussian smoothing algorithm) method^[20]. The correlation coefficient was improved to 0.99 with this new method (data are not shown). Thus this method can be better than other algorithms.

The current study provides a new approach for spectrum alignment and peak extraction for one dimensional NMR spectroscopy. It is fast and efficient. However,it also has a limitation. The searching region was defined manually. As the NMR spectra are diverse,and the influence of the environment is unpredictable,we hope this can be set automatically in the further study.

The algorithm presented in this work provides a computationally efficient and versatile tool for the one-dimensional metabolomics NMR spectroscopy. Although there is a limitation in the procedure of the spectrum alignment,it has better alignment efficiency in various samples. Furthermore,the operator needs only little experience in data analysis or computer skills. It also provides a new variable selection method: peak extraction approach for NMR spectroscopy analysis. This method dramatically decreases the number of the variables for metabolomics study,and each variable has their own chemical information. This should be useful for data interpretation during the NMR analysis. Owing to open source Matlab environment,other researchers can write their personalized code for their own requirement,and establish a set of suitable methods for their own research areas easily.