BACKGROUND
nuclear magnetic resonance spectroscopy  is a powerful and widely applied analytical high-throughput technique to reveal and compare the quantitative metabolic profile of a given tissue in relation to various environmental and clinical parameters.
a typical nmr spectrum is composed out of an x-axis, which indicates the resonance frequencies of the observed molecule, and a y-axis, which denotes the corresponding intensities, i.e., abundance.
to analyse experimental nmr datasets, multivariate methods such as principle components analysis  or univariate techniques like student-t test are commonly applied.
however, chemical and physical sample variations due to, among others, differences in ph, temperature, ion content and the concentration of metabolites, make the analysis of the data challenging.
to address these challenges, several preprocessing steps are commonly applied, including noise reduction, normalization, baseline correction, peak picking and spectrum alignment, prior to statistical analysis.
a crucial and often depreciated aspect in this process is peak alignment, which aims to compensate for small variations in the position of corresponding peaks between spectra.
a number of spectral alignment approaches have previously been proposed.
however, most of them come with particular disadvantages.
for example, some methods use dynamic programming, like correlation optimized warping  and dynamic time warping   <cit> .
due to their computational complexity an alignment task based on these techniques may take hours.
several authors worked towards solutions to speed up this alignment process  <cit>  used a fast fourier transformation  cross-correlation engine to improve the alignment speed .
they also introduced an advanced extension, called recursive peak alignment by fft , which recursively divides the spectrum into meaningful segments and aligns them until a certain degree of goodness is obtained.
some advanced peak picking approaches are recursive segment-wise peak alignment   <cit>  and generalized fuzzy hought transform   <cit> .
other authors applied search algorithms to peak alignment, such as genetic algorithms in paga  <cit>  and beam searching in pabs  <cit> .
recently  <cit> , introduced the interval-correlation-shifting  algorithm, which aligns spectra by maximizing the cross-correlation between user-defined intervals.
another approach that is commonly employed for the peak alignment of mass spectral data is based on hierarchical clustering and could be applied as well on nmr spectral data  <cit> .
most of these methods apply hierarchical clustering to the entire collection of all peaks from the individual spectra and "cut off" the resulting dendrogram at a suitable height to produce a number of clusters used for alignment.
this approach works well on nmr data that is already calibrated to some extent.
however, in some datasets, the peak positions of chemical resonances are significantly shifted between the samples.
this strong shift could make the nmr spectra unclear to separate, which may lead to the wrong clustering, i.e.
alignment, of peaks.
the effect of strongly shifted spectra also challenges the methods based on spectral binning, like cow and icoshift, because peaks could mistakenly be assigned to the wrong bins.
to address the problems with misaligned spectra, we first focus on the development of a robust and highly confident alignment algorithm.
the method is based on a peak-picking approach for nmr spectra, called hierarchical cluster-based peak alignment .
the alignment is embedded in a workflow  that starts with deriving the peak list of each spectrum separately, followed by selecting the reference spectrum to which other spectra are aligned.
the high-quality of the alignment introduced here allows for the use of standard and stable components up- and downstream the processing pipeline.
for example, prior to the peak-picking, noise reduction is accomplished by wavelet filtering.
further more, differential regions in an nmr spectra between two conditions  are determined based on the ratio of the between-group to within-group sums of squares , as proposed by  <cit> .
the bw-ratio is calculated for each data point in the nmr spectra.
because the spectra are well aligned, no variability, due to shifted peaks, is added to the bw-statistic.
implementation
the complete workflow combines four major steps, illustrated in figure  <dig> and itemized below.
each of these major steps is discussed in more detail in separate subsections.
the core of the workflow, is a novel approach for spectral alignment.
<dig>  noise reduction and peak picking: remove instrument noise and collect peaks from the nmr spectra.
<dig>  reference determination: selecting a dominant and representative spectrum, which is used as template for the alignment.
<dig>  peak alignment: shifts are estimated by hierarchical clustering on the selected peaks and the template.
<dig>  quantification: determination of differential region in the nmr spectra based on the bw-ratio.
noise reduction and peak finding
although peak detection is a common step in the analysis of spectral data, it remains a challenging task, due to the noise inherent to many spectrometric techniques.
a typical peak detection approach consists of three steps: smoothing or denoising, baseline correction and the actual peak picking.
for each step, several methods have been developed, but there does not seem to be a "golden standard".
since the detected peaks are the input of the alignment procedure, selecting an appropriate technique is important.
in our implementation, the method of  <cit> , recommended in  <cit>  is chosen.
it employs continuous wavelet transform  to detect peaks in a nmr spectrum.
we used the package massspecwavelet  <cit>  which integrates wavelet-based filters in the smoothing step, continuous wavelet transform for baseline correction, and signal to noise ratio thresholding and ridge lines for the peak picking  <cit> .
in our approach, the method of du et al.
<cit>  could not be applied to the whole spectrum directly, due to computational complexity and possible insensitivity towards low intensity peaks in the nmr spectrum.
to address this, the function peakdetectioncwt of the library massspecwavelet is used only on partitions of the spectrum.
therefore the spectrum was divided into predefined equal window-sized segments.
next, the function peakdetectioncwt is applied to a four-segment sized window that slides along the spectrum, one segment at a time.
all detected peaks are subsequently merged into a final peaklist.
this strategy avoids missing peaks in the margin of a segment.
the segment size was chosen as a power of  <dig>  in order to boost the speed of the cwt  <cit> .
it should be noted that the segment size has to be wide enough to cover the widest peaks in spectra.
noisy peaks are removed by applying a baseline minimum intensity threshold .
the optimal threshold is dataset-dependent, and should be supplied a priori by the user.
for the selection of suitable parameters for the peak detection, we refer to  <cit> .
the user may also use alternative peak detection methods, but it is important to take into account that the quality of the alignment and the downstream analysis can suffer from counterfeit peak detection.
reference determination
the second step in the workflow consists of the selection of a reference spectrum.
the reference spectrum serves as a kind of template for other spectra to be aligned against.
the ideal reference contains the highest number of common chemical constituents  <cit> .
several methods exist to select an appropriate reference spectrum: either prior knowledge about the dataset can be used, or the reference can be selected according to a predefined criterion  <cit>  propose a method based on the product of the pearson correlation coefficients between spectra.
the approach of  <cit>  creates a virtual "average" spectrum as a reference.
alternatively, the spectrum that is most similar to the loading of the first principle component in a pca model can be chosen  <cit> .
it is difficult to define "common chemical constituents" when spectra are ill aligned.
to ensure robust reference selection, we introduce an alternative method that incorporates a heuristic strategy to find the optimal template.
it selects the spectrum with the largest goodness as the reference, which is defined by the sum of the inter-spectrum distances, according to formula  <dig> and 2:
goodness=-∑t∑pt∈tminps∈s
ref=argmaxs)
where ps and pt represent the resonance frequencies, i.e., location along the x-axis, of the corresponding peaks of spectrum s and t, respectively.
the distance is computed by the sum of the minimum distances between every peak pt of t and the peaks ps of s, because the closest peak is chosen as the best candidate to be matched with pt.
note, that the distance is not symmetric, that is the distance of s to t and the distance of t to s are not necessarily equal.
in principle, the method projects each high-dimensional nmr spectrum as a data-point in a one-dimensional space and selects the data-point in the middle of the data cloud as the reference.
a possible disadvantage of this method is that it favours the spectrum with a high number of peaks, which is not necessarily in the center of the data cloud.
if that is the case, a lower ranked spectrum could be selected instead.
ideally, the perfect reference candidate has enough peaks corresponding to chemical resonance areas, i.e.
frequencies and the peaks of the reference stay at the center of the distribution of peaks of other spectra in order to minimize the distance.
the computational complexity to find the reference spectrum includes o for finding the minimum value and the second sum in formula ; and o for the other sum of formula  and the maximum of formula ; where n is the maximum number of peaks of a spectrum and m is the number of spectra.
the total complexity is then o * o * o + o) or o or o. it should be noticed that good peak detection is required to restrain the computational complexity within acceptable limits.
in practice, finding the reference spectrum is fast, since the number of peaks n of a spectrum is small, given that peak selection up-stream of the workflow is adequate.
spectrum alignment
the alignment method is called hierarchical cluster-based peak alignment , and is based on shifting segments of the original spectra, first globally and then in smaller local segments.
in the first iteration, clupa does the alignment on the whole spectrum.
in subsequent iterations, clupa recursively splits the segment of previous steps into two smaller segments, and applies an individual alignment shift to each of these.
the boundaries of these segments are defined by applying a hierarchical clustering algorithm on the combined peak list of the reference and the target spectrum.
the clustering is based on the distance between peaks: peaks are hierarchically combined into a tree according to their distance.
after each round, the segment is split into two smaller segments by cutting at the second level of the hierarchical tree.
at each iteration, the present segments of the original spectra are aligned to the reference spectrum.
to this end, an optimal shift step is determined between the target and reference spectrum within the present segments.
this shift is subsequently applied to the corresponding segment in the target spectrum.
the process is repeated until a segment contains only peaks of a single spectrum, or until there are only two peaks left in the list.
a detailed outline of the algorithm is presented in the following pseudo-code.
input
• targetspec: original spectrum to be aligned
• refspec: original reference spectrum
• mergedpeaklist: merged list of all peaks detected in targetspec and refspec
output
• targetspec: spectrum after alignment
clupa
<dig>  startpoint, endpoint = findshiftregion
<dig>  shiftstep = findshiftstep
<dig>  targetspec, mergedpeaklist = shift
<dig>  hclustertree = buildhcluster
<dig>  leftpeaklist = cutlefttree
<dig>  if satifiedcondition then
<dig>  targetspec = clupa
<dig>  rightpeaklist = cutrighttree
<dig>  if satifiedcondition then
<dig>  targetspec = clupa
<dig>  return
step 1: findshiftregion: determines the region of the spectra for alignment, which must contain all peaks of the mergedpeaklist.
the left most of the region is the lowest intensity point in the spectral region on the left side of the list.
similarly, the end point of the region is the lowest intensity point in the spectral region on the right side of the list.
it takes linear time o to select the lowest intensity point, where d is the length of the spectral region.
this value rapidly decreases when the algorithm iterates with increasingly smaller regions.
step 2: findshiftstep: finds a shift step by using fft cross-correlation.
the step has a computational cost o. in implementation, the function is simulated from the function fftcorr in rafft  <cit> .
alternatives, for example pearson correlation and weighted cross correlation in  <cit>  could also be applied as an alternative, however, pearson correlation takes more time than using fft cross correlation.
step 3: shift: moves shiftstep - the number of points on the current region of targetspec spectrum and peak list.
this simple step takes linear time o to complete.
step 4: buildhcluster: creates a hierarchical cluster tree from the peaks of the mergedpeaklist.
in the first step, each cluster has a single peak.
the hierarchical cluster tree is built by in turn grouping the closest two clusters together.
there are three common ways to define the distance between two clusters: single linkage, complete linkage and average linkage  <cit> .
single linkage distance of two clusters is the minimum distance of two objects or peaks in each cluster, complete linkage distance is computed by maximum distance of the two objects, and average linkage is calculated by the average of all distances between any pairs of objects between the two clusters.
by default, the average linkage is selected as the distance for clustering, however, other types of distance could be applied.
this step takes o to o  <cit> , depending on which linkage distance is selected, where n is the number of elements in the peak list.
the calculation of single linkage distance is faster than its alternatives.
in the implementation, the standard function hclust() in the stats package of the r software ) is applied.
step 5: cutlefttree: gets the peaks corresponding to the left tree at the second level of hclustertree.
step 8: cutrighttree: gets the peaks corresponding to the right tree at the second level of hclustertree.
in the implementation, the function cutree() of the stats package is applied.
it takes o to be finished.
step  <dig> and 9: satifiedcondition verifies whether the condition that determines whether a tree should be further divided are met.
the tree has to contain peaks of both reference and target spectra, otherwise the spectral region does not need to be further aligned.
the maximum computational complexity in each stage is o + o + o + o +  <dig> * o) or o, and o for the recursive loop.
therefore, the total complexity of the algorithm is o * o.
quantitative analysis
the thorough alignment of the nmr spectra, outlined in previous section, allows for a very straightforward and simple quantification analysis.
the data returned by clupa are well aligned and presented in a rectangular format.
this means that the data are represented as a matrix in which the rows denote the spectra and the columns indicate the resonance frequency instances.
the quantification strategy is inspired by the method proposed by tan et al  <cit>  for annotating regions of significance in seldi-tof spectra.
in our approach, for each frequency instance in the matrix, we performed an analysis of the data based on the ratio of the between-group to within-group sums of squares  as proposed by dudoit et al  <cit> .
this ratio is related to the f-statistic or one-way anova, but we do not make assumptions about the distribution and the independency of the data.
for example, the bw-ratio of frequency instance j is calculated as
bw=ΣiΣki2ΣiΣki <dig>
where x.j- and xkj- denote the average intensity level of frequency instance j between all nmr spectra and within the nmr spectra belonging to group k only, respectively.
the variable i is an indicator for the group.
the variable xij denotes the i-th intensity values of frequency j belonging to group k. statistical inference on the bw-statistic is based on the sampled null distribution from the observed data.
for this purpose the average of the groups xkj- is subtracted from the data.
next, data points are sampled for each frequency instance j using the same grouping as in the original setting.
based on the sampled data, the bw-ratio is calculated, which is representative for data showing no difference in the group means.
this computation is iterated with replacement of the data until a well sampled null distribution is obtained.
for each frequency j, the value of the bw-statistic is determined for the desired alpha-level based on the sampled null distribution.
we are aware that nmr-spectra are correlated between consecutive frequency instances and cannot be treated as independent variables.
a model which takes into account such a correlation structure to discern noise peaks from bona fide molecule peaks is part of future research.
in addition, testing for every possible data entry can make the detection of differential frequency instances inefficient, because a valid nmr-peak can be composed out of multiple, redundant, data points.
therefore, correction of multiplicity is done by using the conservative bonferoni correction.
doing so, the alpha-level is adjusted for the average number of peaks returned by the peak picking algorithm used in the first step of the workflow.
we can already announce that this ad-hoc quantitative analysis is very fast and robust.
an interesting point is that the differential frequencies do not necessarily represent the apex of the nmr-peak, but are rather found in the shoulder of the peak.
we also want to stress that such an ad-hoc approach is only suitable for well aligned data.
misalignment would increase the variability, making it difficult to discriminate differential frequencies.
RESULTS
the evaluation of the workflow proposed in this paper is outlined in three sections.
first, we compare the alignment method against other available methods.
second, we evaluate the quantitative power of the approach.
third, we discuss the efficiency of the workflow.
spectral alignment evaluation
we compared the proposed method to two alternatives that use fft cross correlation to measure correlation: icoshift  <cit>  and rafft  <cit> , the comparison was made using two different datasets: 1) the wine dataset, which is a known benchmark nmr dataset  <cit> ; 2) a huntington dataset  <cit> .
details regarding the data and the methodology are outlined in the methods section at the end of this paper, and in the original papers that describe the data  <cit> .
as mentioned in  <cit> , most methods such as cow, rspa and icoshift perform well on certain regions of the spectra, but perform less in specific regions, such as the lactate and ethanol region.
figure  <dig> shows gray scale representations and spectrum plots of the unaligned spectra and the results of the different alignment methods within this region.
rafft did not properly align this region.
icoshift fails to properly align the lactic region in the automatic interval selection mode.
however, when we specified the intervals from  <cit>  using the user-defined interval mode, a satisfying alignment was obtained.
clupa aligns well on all main resonance peaks mentioned in  <cit> , without the requirement of user interaction.
in the case of the lactic acid region, clupa yields a good alignment, as long the peak detection algorithm successfully detects the peaks in this region.
the same methods were also applied to the huntington dataset.
we show two types of selected regions: a region containing the highest peaks and the largest region containing strongly overlapping peaks.
a) region with highest peaks: between  <dig> and  <dig> or  <dig>  to  <dig>  ppm
this region  is also the lactic acid region, similar to the wine dataset, in which there are many small peaks around a very intense peak.
considering the grey scale plots, rafft, icoshift and clupa align the highest peak region in an almost identical way.
the spectral plots, which boost the lower intensity scales, however show that clupa generally outperforms the other alignment options with regards to the lower intensity peaks surrounding the highest peaks.
nevertheless there are still some misalignments, potentially caused by the fact that the peak detection algorithm does not fully detect these peaks.
this problem may be overcome by adjusting the peak detection parameters.
the grey-scale plots also reveal spectra  <dig> and  <dig> as commonly misaligned, in accordance to their outlier status .
b) region with high overlap: from  <dig> to  <dig> index value or from  <dig>  until  <dig>  ppm
the large number of peaks in this region  makes the alignment of spectra by a binning method challenging, since the bin may contain the wrong peaks for the alignment.
icoshift does give a very good alignment for some spectra, but it is not always successful.
in some cases, the peaks are not fully aligned together and there are some artefact lines in the spectral plots of the icoshift results.
both rafft and clupa show better results in this region.
the grey scale pictures show that rafft gives a perfect alignment for the left side of the region and yields an adequate improvement in the right side region.
the best alignment is obtained by the clupa algorithm, which correctly matches the high peaks of the region.
the grey scale plots also reveals the noisy spectra  <dig>   <dig> and  <dig> as outliers.
a different way to assess the efficiency of the alignment algorithm consists of a correlation analysis between spectra.
to this aim, the pearson correlation coefficient was calculated between each pair of spectra within the  dataset.
a visual representation of all pairwise correlations for each method is shown in figure  <dig>  in general, all alignment methods yield an apparent improvement.
although manual-icoshift works well on all main resonance peaks , its correlation map shows lesser correlation for several spectra.
the reason may be found in the fact that manual-icoshift does alignment on some user-defined regions of the spectra, and does not fully consider all data points.
the correlation maps of rafft, auto-icoshift and clupa are similar.
we proposed a heuristic method to find the reference spectrum where other spectra are aligned against.
the method selects the spectrum with the maximum goodness, which is computed from the distance between the peaks of the reference candidate and the others.
to further evaluate the impact of reference selection, we used a simple test based on the average correlation values, obtained as explained above.
each sample was selected once as reference, and the others were aligned against it through the clupa alignment algorithm.
the distribution of the average correlation values for each run was drawn, and the correlation value corresponding to the use of the reference found by the heuristic method was marked by red lines, in figure  <dig> a and  <dig> b. in the case of the wine dataset, the correlations are high with a narrow distribution.
the reference is close to the center of this distribution.
in the case of the huntington dataset, the heuristic selects a reference in the center of the distribution of the non-outliers, clearly demonstrating that the heuristic method works.
this is further confirmed by the fact that this heuristic ranks the outlier spectra  lowest among all candidate reference spectra .
quantitative analysis
the quantitative analysis was performed as outlined in the workflow section.
the alpha-level of  <dig>  was corrected for the  <dig> peaks found above an intensity of  <dig> .
the analysis returned  <dig> regions of interest, which had a significant bw-ratio.
the top panel of figure  <dig> a displays the significant regions over the entire x-axis.
the blue line indicates the bw-statistic and the black line denotes the bw-statistic at a percentile of alpha/ <dig> using the null distribution.
the null distribution is composed out of  <dig> samples.
bw-values above the black line indicate a differential region.
the plot in the bottom panel displays the nmr-spectra for the red and white wine.
from the figures it becomes clear that differential regions do not necessarily correspond to regions with intense peak abundances because the bw-statistic takes into account the variability between the groups instead of looking at absolute differences.
in order to demonstrate the proper functioning of the approach, we present several close-ups of interesting regions.
for example, figure  <dig> b presents a close-up of the most intense peaks in the data .
from the top panel it is clear that this region is not selected as a differential region; the bw-statistics stays under the critical value corresponding to the significance level.
a different scenario is observed in figure  <dig> c and  <dig> d, where differential regions are detected.
a remarkable observation is that in figure  <dig> c and  <dig> d the maxmimum bw-statistic is not reached at the apex of the nmr-peaks but in the ascending and descending shoulder of the peak.
it seems that the due to lower variability in the shoulder, i.e., more robust measurements, it should be easier to detect differences.
an approach which focuses on finding differences in the peak heights might discard valuable information.
we compared these results with the information in the original research paper  <cit> , and found that most findings are in accordance to the published results.
in figure  <dig> b, the highest intensity region or the ethanol region ranging  <dig> - <dig>  ppm is not a good region for classification because the authors had to use a very highly complex ipls model with  <dig> components for it.
leave one out for rmsecv and r <dig>  are  <dig>  and 84%, however, the authors reported that they were overfitting due to the complex model  <cit> .
figure  <dig> c is the methanol region ranging  <dig> - <dig>  ppm.
the result from the paper showed that they used simpler ipls model with  <dig> components, achieve  <dig>  and 75% for rmsecv and r <dig>  respectively.
this means that this region is suitable for classification .
we do not have quantification information corresponding to figure  <dig> d, or the acid acic region.
algorithm speed and parameter settings
since the alignment algorithm  requires most of the computational resources of the workflow, in comparison with the statistical analysis, we briefly discuss its operational characteristics in this section.
on a macbook pro  <dig>  ghz  the whole process needs minutes to complete.
as other peak based alignment methods, the peak detection process is the most time consuming.
the alignment of the wine dataset takes totally  <dig>  minutes, consisting of  <dig>  minutes for peak detection,  <dig>  minutes for finding the reference and  <dig>  minutes for the alignment process.
the huntington dataset required nearly four times the time for the wine dataset:  <dig>  minutes.
most of the time cost is reserved for peak detection .
the remainder includes  <dig>  minutes for finding the reference spectrum and  <dig>  minutes for the alignment.
besides the parameters for the peak detection, the workflow does not require parameters for the initialization prior to running.
one extra optional parameter for clupa, but also for both rafft and icoshift, is the maximum shift point, which is inherited from the fft cross-correlation function.
this value is usually  <dig>  ppm, more or less  <dig> points, depending on the data.
if this parameter is not set, the fft cross-correlation function will check on whole spectral regions to get the optimal shift step.
furthermore, users need to set the parameters for the peak detection algorithm or use other methods instead.
even though default values are set in the software implementation, they are not guaranteed to work optimally for any dataset.
the selection of optimal parameters for peak detection algorithm is beyond the scope of this paper.
CONCLUSIONS
in this paper, we proposed a simple and efficient workflow which is centered around recursive hierarchical clustering.
this strategy, which is referred to as speaq  is feasible because prior to the alignment, a peak-picking and reference selection approach is used to reduce the complexity related to alignment.
the method aligns the target spectrum to the reference spectrum in a top-down fashion and makes use of fast fourier transformation  cross-correlation to reduce the computation time to find the shift step.
as conventional recursive strategy like rspa or rafft, the clupa algorithm aligns a larger segment and then recursively divides it into smaller segments to refine the alignment.
the differences between them lie in the method and criteria to divide the spectra.
rafft selects the dividing point based on the optimal shifts computed from fast fourier transform cross-correlation, rspa detects in the reference and target spectrum lists of segments, then tries to find the corresponding segments between them, and goes into segments for further alignment.
our proposed method builds a hierarchical cluster tree from peak lists of reference and target, and divides the spectra into smaller regions based on the most distant clusters of the tree.
however, instead of considering many clusters simultaneously, we cut off the dendrogram into just two sub trees  and then go further into each of them for finer alignment.
another feature of the workflow which is different from most other peak-picking based alignment methods is that instead of grouping peaks together to create segments for alignment, speaq does the alignment first, then it groups peaks and creates smaller segments.
in general, this strategy decreases the shift of peaks before the hierarchical cluster tree is built, thus, it significantly reduces overlapping peaks that yield wrong clustering.
despite the heuristic approach, selecting the perfect reference remains a non-trivial task, since one spectrum may be a good reference for a specific region but not for other regions.
within the software implementation, we therefore also allow users to intervene into the process by setting reference in the region they are confident about.
this also allows to skip certain regions during the alignment.
the maximum shift can also be set in a region-specific manner.
all this information can be readily provided to the software in an extra text file.
the speaq workflow was compared to two widely used methods icoshift and rafft.
the method shows several advantages over the others.
it is easy to use because users do not need to set many initial parameters.
experiments on both a public dataset and an huntington dataset show that it performs very well compared to the others.
all method use fft cross correlation for finding the shift because of its advantage in computational time.
however, speaq is slower than both icoshift and rafft, due to the time for peak detection and its recursive nature, but this is generally not a problem with the powerful computers we have nowadays.
in summary, the speaq workflow is unique in the way the subsequent steps are ordered.
because of rearranging the order of processing steps, simple and robust methods could be adopted for the workflow.
a case in point is the effective analysis based on the bw-ratio.
this could only be achieved since nmr spectra are well aligned.
if this were not the case, misalignment would introduce sources of variability into the statistic making it difficult to select differential regions.
i should be noticed that the proposed method is still open for improvements.
for example, from the top panel of figure  <dig> a it can be seen that the null hypothesis  is remarkably constant.
this constant line might suggest that the properties of the null distribution are applicable over the entire frequency range.
parametrization of the method can avoid the sampling of the null distribution and improve computational efficiency.
the method currently does not take into account covariation of the intensities of peaks of the same molecule across samples  <cit> .
the method was tested only on nmr datasets.
similar data types, such as chromatography data, may also benefit from the approach presented here.

 nuclear magnetic resonance spectroscopy is a powerful and widely applied analytical high-throughput technique to reveal and compare the quantitative metabolic profile of a given tissue in relation to various environmental and clinical parameters. a typical nmr spectrum is composed out of an x-axis, which indicates the resonance frequencies of the observed molecule, and a y-axis, which denotes the corresponding intensities, i.e., abundance. to analyse experimental nmr datasets, multivariate methods such as principle components analysis or univariate techniques like student-t test are commonly applied. however, chemical and physical sample variations due to, among others, differences in ph, temperature, ion content and the concentration of metabolites, make the analysis of the data challenging. to address these challenges, several preprocessing steps are commonly applied, including noise reduction, normalization, baseline correction, peak picking and spectrum alignment, prior to statistical analysis. a crucial and often depreciated aspect in this process is peak alignment, which aims to compensate for small variations in the position of corresponding peaks between spectra. a number of spectral alignment approaches have previously been proposed. however, most of them come with particular disadvantages. for example, some methods use dynamic programming, like correlation optimized warping and dynamic time warping. due to their computational complexity an alignment task based on these techniques may take hours. several authors worked towards solutions to speed up this alignment process used a fast fourier transformation cross-correlation engine to improve the alignment speed. they also introduced an advanced extension, called recursive peak alignment by fft, which recursively divides the spectrum into meaningful segments and aligns them until a certain degree of goodness is obtained. some advanced peak picking approaches are recursive segment-wise peak alignment and generalized fuzzy hought transform. other authors applied search algorithms to peak alignment, such as genetic algorithms in paga and beam searching in pabs. recently, introduced the interval-correlation-shifting algorithm, which aligns spectra by maximizing the cross-correlation between user-defined intervals. another approach that is commonly employed for the peak alignment of mass spectral data is based on hierarchical clustering and could be applied as well on nmr spectral data. most of these methods apply hierarchical clustering to the entire collection of all peaks from the individual spectra and "cut off" the resulting dendrogram at a suitable height to produce a number of clusters used for alignment. this approach works well on nmr data that is already calibrated to some extent. however, in some datasets, the peak positions of chemical resonances are significantly shifted between the samples. this strong shift could make the nmr spectra unclear to separate, which may lead to the wrong clustering, i.e. alignment, of peaks. the effect of strongly shifted spectra also challenges the methods based on spectral binning, like cow and icoshift, because peaks could mistakenly be assigned to the wrong bins. to address the problems with misaligned spectra, we first focus on the development of a robust and highly confident alignment algorithm. the method is based on a peak-picking approach for nmr spectra, called hierarchical cluster-based peak alignment. the alignment is embedded in a workflow that starts with deriving the peak list of each spectrum separately, followed by selecting the reference spectrum to which other spectra are aligned. the high-quality of the alignment introduced here allows for the use of standard and stable components up- and downstream the processing pipeline. for example, prior to the peak-picking, noise reduction is accomplished by wavelet filtering. further more, differential regions in an nmr spectra between two conditions are determined based on the ratio of the between-group to within-group sums of squares, as proposed by. the bw-ratio is calculated for each data point in the nmr spectra. because the spectra are well aligned, no variability, due to shifted peaks, is added to the bw-statistic. implementation the complete workflow combines four major steps, illustrated in figure and itemized below. each of these major steps is discussed in more detail in separate subsections. the core of the workflow, is a novel approach for spectral alignment. noise reduction and peak picking: remove instrument noise and collect peaks from the nmr spectra. reference determination: selecting a dominant and representative spectrum, which is used as template for the alignment. peak alignment: shifts are estimated by hierarchical clustering on the selected peaks and the template. quantification: determination of differential region in the nmr spectra based on the bw-ratio. noise reduction and peak finding although peak detection is a common step in the analysis of spectral data, it remains a challenging task, due to the noise inherent to many spectrometric techniques. a typical peak detection approach consists of three steps: smoothing or denoising, baseline correction and the actual peak picking. for each step, several methods have been developed, but there does not seem to be a "golden standard". since the detected peaks are the input of the alignment procedure, selecting an appropriate technique is important. in our implementation, the method of, recommended in is chosen. it employs continuous wavelet transform to detect peaks in a nmr spectrum. we used the package massspecwavelet which integrates wavelet-based filters in the smoothing step, continuous wavelet transform for baseline correction, and signal to noise ratio thresholding and ridge lines for the peak picking. in our approach, the method of du et al. could not be applied to the whole spectrum directly, due to computational complexity and possible insensitivity towards low intensity peaks in the nmr spectrum. to address this, the function peakdetectioncwt of the library massspecwavelet is used only on partitions of the spectrum. therefore the spectrum was divided into predefined equal window-sized segments. next, the function peakdetectioncwt is applied to a four-segment sized window that slides along the spectrum, one segment at a time. all detected peaks are subsequently merged into a final peaklist. this strategy avoids missing peaks in the margin of a segment. the segment size was chosen as a power of in order to boost the speed of the cwt. it should be noted that the segment size has to be wide enough to cover the widest peaks in spectra. noisy peaks are removed by applying a baseline minimum intensity threshold. the optimal threshold is dataset-dependent, and should be supplied a priori by the user. for the selection of suitable parameters for the peak detection, we refer to. the user may also use alternative peak detection methods, but it is important to take into account that the quality of the alignment and the downstream analysis can suffer from counterfeit peak detection. reference determination the second step in the workflow consists of the selection of a reference spectrum. the reference spectrum serves as a kind of template for other spectra to be aligned against. the ideal reference contains the highest number of common chemical constituents. several methods exist to select an appropriate reference spectrum: either prior knowledge about the dataset can be used, or the reference can be selected according to a predefined criterion propose a method based on the product of the pearson correlation coefficients between spectra. the approach of creates a virtual "average" spectrum as a reference. alternatively, the spectrum that is most similar to the loading of the first principle component in a pca model can be chosen. it is difficult to define "common chemical constituents" when spectra are ill aligned. to ensure robust reference selection, we introduce an alternative method that incorporates a heuristic strategy to find the optimal template. it selects the spectrum with the largest goodness as the reference, which is defined by the sum of the inter-spectrum distances, according to formula and 2: goodness=-∑t∑pt∈tminps∈s ref=argmaxs) where ps and pt represent the resonance frequencies, i.e., location along the x-axis, of the corresponding peaks of spectrum s and t, respectively. the distance is computed by the sum of the minimum distances between every peak pt of t and the peaks ps of s, because the closest peak is chosen as the best candidate to be matched with pt. note, that the distance is not symmetric, that is the distance of s to t and the distance of t to s are not necessarily equal. in principle, the method projects each high-dimensional nmr spectrum as a data-point in a one-dimensional space and selects the data-point in the middle of the data cloud as the reference. a possible disadvantage of this method is that it favours the spectrum with a high number of peaks, which is not necessarily in the center of the data cloud. if that is the case, a lower ranked spectrum could be selected instead. ideally, the perfect reference candidate has enough peaks corresponding to chemical resonance areas, i.e. frequencies and the peaks of the reference stay at the center of the distribution of peaks of other spectra in order to minimize the distance. the computational complexity to find the reference spectrum includes o for finding the minimum value and the second sum in formula; and o for the other sum of formula and the maximum of formula; where n is the maximum number of peaks of a spectrum and m is the number of spectra. the total complexity is then o * o * o + o) or o or o. it should be noticed that good peak detection is required to restrain the computational complexity within acceptable limits. in practice, finding the reference spectrum is fast, since the number of peaks n of a spectrum is small, given that peak selection up-stream of the workflow is adequate. spectrum alignment the alignment method is called hierarchical cluster-based peak alignment, and is based on shifting segments of the original spectra, first globally and then in smaller local segments. in the first iteration, clupa does the alignment on the whole spectrum. in subsequent iterations, clupa recursively splits the segment of previous steps into two smaller segments, and applies an individual alignment shift to each of these. the boundaries of these segments are defined by applying a hierarchical clustering algorithm on the combined peak list of the reference and the target spectrum. the clustering is based on the distance between peaks: peaks are hierarchically combined into a tree according to their distance. after each round, the segment is split into two smaller segments by cutting at the second level of the hierarchical tree. at each iteration, the present segments of the original spectra are aligned to the reference spectrum. to this end, an optimal shift step is determined between the target and reference spectrum within the present segments. this shift is subsequently applied to the corresponding segment in the target spectrum. the process is repeated until a segment contains only peaks of a single spectrum, or until there are only two peaks left in the list. a detailed outline of the algorithm is presented in the following pseudo-code. input • targetspec: original spectrum to be aligned • refspec: original reference spectrum • mergedpeaklist: merged list of all peaks detected in targetspec and refspec output • targetspec: spectrum after alignment clupa startpoint, endpoint = findshiftregion shiftstep = findshiftstep targetspec, mergedpeaklist = shift hclustertree = buildhcluster leftpeaklist = cutlefttree if satifiedcondition then targetspec = clupa rightpeaklist = cutrighttree if satifiedcondition then targetspec = clupa return step 1: findshiftregion: determines the region of the spectra for alignment, which must contain all peaks of the mergedpeaklist. the left most of the region is the lowest intensity point in the spectral region on the left side of the list. similarly, the end point of the region is the lowest intensity point in the spectral region on the right side of the list. it takes linear time o to select the lowest intensity point, where d is the length of the spectral region. this value rapidly decreases when the algorithm iterates with increasingly smaller regions. step 2: findshiftstep: finds a shift step by using fft cross-correlation. the step has a computational cost o. in implementation, the function is simulated from the function fftcorr in rafft. alternatives, for example pearson correlation and weighted cross correlation in could also be applied as an alternative, however, pearson correlation takes more time than using fft cross correlation. step 3: shift: moves shiftstep - the number of points on the current region of targetspec spectrum and peak list. this simple step takes linear time o to complete. step 4: buildhcluster: creates a hierarchical cluster tree from the peaks of the mergedpeaklist. in the first step, each cluster has a single peak. the hierarchical cluster tree is built by in turn grouping the closest two clusters together. there are three common ways to define the distance between two clusters: single linkage, complete linkage and average linkage. single linkage distance of two clusters is the minimum distance of two objects or peaks in each cluster, complete linkage distance is computed by maximum distance of the two objects, and average linkage is calculated by the average of all distances between any pairs of objects between the two clusters. by default, the average linkage is selected as the distance for clustering, however, other types of distance could be applied. this step takes o to o, depending on which linkage distance is selected, where n is the number of elements in the peak list. the calculation of single linkage distance is faster than its alternatives. in the implementation, the standard function hclust() in the stats package of the r software ) is applied. step 5: cutlefttree: gets the peaks corresponding to the left tree at the second level of hclustertree. step 8: cutrighttree: gets the peaks corresponding to the right tree at the second level of hclustertree. in the implementation, the function cutree() of the stats package is applied. it takes o to be finished. step and 9: satifiedcondition verifies whether the condition that determines whether a tree should be further divided are met. the tree has to contain peaks of both reference and target spectra, otherwise the spectral region does not need to be further aligned. the maximum computational complexity in each stage is o + o + o + o + * o) or o, and o for the recursive loop. therefore, the total complexity of the algorithm is o * o. quantitative analysis the thorough alignment of the nmr spectra, outlined in previous section, allows for a very straightforward and simple quantification analysis. the data returned by clupa are well aligned and presented in a rectangular format. this means that the data are represented as a matrix in which the rows denote the spectra and the columns indicate the resonance frequency instances. the quantification strategy is inspired by the method proposed by tan et al for annotating regions of significance in seldi-tof spectra. in our approach, for each frequency instance in the matrix, we performed an analysis of the data based on the ratio of the between-group to within-group sums of squares as proposed by dudoit et al. this ratio is related to the f-statistic or one-way anova, but we do not make assumptions about the distribution and the independency of the data. for example, the bw-ratio of frequency instance j is calculated as bw=ΣiΣki2ΣiΣki where x.j- and xkj- denote the average intensity level of frequency instance j between all nmr spectra and within the nmr spectra belonging to group k only, respectively. the variable i is an indicator for the group. the variable xij denotes the i-th intensity values of frequency j belonging to group k. statistical inference on the bw-statistic is based on the sampled null distribution from the observed data. for this purpose the average of the groups xkj- is subtracted from the data. next, data points are sampled for each frequency instance j using the same grouping as in the original setting. based on the sampled data, the bw-ratio is calculated, which is representative for data showing no difference in the group means. this computation is iterated with replacement of the data until a well sampled null distribution is obtained. for each frequency j, the value of the bw-statistic is determined for the desired alpha-level based on the sampled null distribution. we are aware that nmr-spectra are correlated between consecutive frequency instances and cannot be treated as independent variables. a model which takes into account such a correlation structure to discern noise peaks from bona fide molecule peaks is part of future research. in addition, testing for every possible data entry can make the detection of differential frequency instances inefficient, because a valid nmr-peak can be composed out of multiple, redundant, data points. therefore, correction of multiplicity is done by using the conservative bonferoni correction. doing so, the alpha-level is adjusted for the average number of peaks returned by the peak picking algorithm used in the first step of the workflow. we can already announce that this ad-hoc quantitative analysis is very fast and robust. an interesting point is that the differential frequencies do not necessarily represent the apex of the nmr-peak, but are rather found in the shoulder of the peak. we also want to stress that such an ad-hoc approach is only suitable for well aligned data. misalignment would increase the variability, making it difficult to discriminate differential frequencies.  the evaluation of the workflow proposed in this paper is outlined in three sections. first, we compare the alignment method against other available methods. second, we evaluate the quantitative power of the approach. third, we discuss the efficiency of the workflow. spectral alignment evaluation we compared the proposed method to two alternatives that use fft cross correlation to measure correlation: icoshift and rafft, the comparison was made using two different datasets: 1) the wine dataset, which is a known benchmark nmr dataset; 2) a huntington dataset. details regarding the data and the methodology are outlined in the methods section at the end of this paper, and in the original papers that describe the data. as mentioned in, most methods such as cow, rspa and icoshift perform well on certain regions of the spectra, but perform less in specific regions, such as the lactate and ethanol region. figure shows gray scale representations and spectrum plots of the unaligned spectra and the results of the different alignment methods within this region. rafft did not properly align this region. icoshift fails to properly align the lactic region in the automatic interval selection mode. however, when we specified the intervals from using the user-defined interval mode, a satisfying alignment was obtained. clupa aligns well on all main resonance peaks mentioned in, without the requirement of user interaction. in the case of the lactic acid region, clupa yields a good alignment, as long the peak detection algorithm successfully detects the peaks in this region. the same methods were also applied to the huntington dataset. we show two types of selected regions: a region containing the highest peaks and the largest region containing strongly overlapping peaks. a) region with highest peaks: between and or to ppm this region is also the lactic acid region, similar to the wine dataset, in which there are many small peaks around a very intense peak. considering the grey scale plots, rafft, icoshift and clupa align the highest peak region in an almost identical way. the spectral plots, which boost the lower intensity scales, however show that clupa generally outperforms the other alignment options with regards to the lower intensity peaks surrounding the highest peaks. nevertheless there are still some misalignments, potentially caused by the fact that the peak detection algorithm does not fully detect these peaks. this problem may be overcome by adjusting the peak detection parameters. the grey-scale plots also reveal spectra and as commonly misaligned, in accordance to their outlier status. b) region with high overlap: from to index value or from until ppm the large number of peaks in this region makes the alignment of spectra by a binning method challenging, since the bin may contain the wrong peaks for the alignment. icoshift does give a very good alignment for some spectra, but it is not always successful. in some cases, the peaks are not fully aligned together and there are some artefact lines in the spectral plots of the icoshift results. both rafft and clupa show better results in this region. the grey scale pictures show that rafft gives a perfect alignment for the left side of the region and yields an adequate improvement in the right side region. the best alignment is obtained by the clupa algorithm, which correctly matches the high peaks of the region. the grey scale plots also reveals the noisy spectra and as outliers. a different way to assess the efficiency of the alignment algorithm consists of a correlation analysis between spectra. to this aim, the pearson correlation coefficient was calculated between each pair of spectra within the dataset. a visual representation of all pairwise correlations for each method is shown in figure in general, all alignment methods yield an apparent improvement. although manual-icoshift works well on all main resonance peaks, its correlation map shows lesser correlation for several spectra. the reason may be found in the fact that manual-icoshift does alignment on some user-defined regions of the spectra, and does not fully consider all data points. the correlation maps of rafft, auto-icoshift and clupa are similar. we proposed a heuristic method to find the reference spectrum where other spectra are aligned against. the method selects the spectrum with the maximum goodness, which is computed from the distance between the peaks of the reference candidate and the others. to further evaluate the impact of reference selection, we used a simple test based on the average correlation values, obtained as explained above. each sample was selected once as reference, and the others were aligned against it through the clupa alignment algorithm. the distribution of the average correlation values for each run was drawn, and the correlation value corresponding to the use of the reference found by the heuristic method was marked by red lines, in figure a and b. in the case of the wine dataset, the correlations are high with a narrow distribution. the reference is close to the center of this distribution. in the case of the huntington dataset, the heuristic selects a reference in the center of the distribution of the non-outliers, clearly demonstrating that the heuristic method works. this is further confirmed by the fact that this heuristic ranks the outlier spectra lowest among all candidate reference spectra. quantitative analysis the quantitative analysis was performed as outlined in the workflow section. the alpha-level of was corrected for the peaks found above an intensity of. the analysis returned regions of interest, which had a significant bw-ratio. the top panel of figure a displays the significant regions over the entire x-axis. the blue line indicates the bw-statistic and the black line denotes the bw-statistic at a percentile of alpha/ using the null distribution. the null distribution is composed out of samples. bw-values above the black line indicate a differential region. the plot in the bottom panel displays the nmr-spectra for the red and white wine. from the figures it becomes clear that differential regions do not necessarily correspond to regions with intense peak abundances because the bw-statistic takes into account the variability between the groups instead of looking at absolute differences. in order to demonstrate the proper functioning of the approach, we present several close-ups of interesting regions. for example, figure b presents a close-up of the most intense peaks in the data. from the top panel it is clear that this region is not selected as a differential region; the bw-statistics stays under the critical value corresponding to the significance level. a different scenario is observed in figure c and d, where differential regions are detected. a remarkable observation is that in figure c and d the maxmimum bw-statistic is not reached at the apex of the nmr-peaks but in the ascending and descending shoulder of the peak. it seems that the due to lower variability in the shoulder, i.e., more robust measurements, it should be easier to detect differences. an approach which focuses on finding differences in the peak heights might discard valuable information. we compared these results with the information in the original research paper, and found that most findings are in accordance to the published results. in figure b, the highest intensity region or the ethanol region ranging - ppm is not a good region for classification because the authors had to use a very highly complex ipls model with components for it. leave one out for rmsecv and r are and 84%, however, the authors reported that they were overfitting due to the complex model. figure c is the methanol region ranging - ppm. the result from the paper showed that they used simpler ipls model with components, achieve and 75% for rmsecv and r respectively. this means that this region is suitable for classification. we do not have quantification information corresponding to figure d, or the acid acic region. algorithm speed and parameter settings since the alignment algorithm requires most of the computational resources of the workflow, in comparison with the statistical analysis, we briefly discuss its operational characteristics in this section. on a macbook pro ghz the whole process needs minutes to complete. as other peak based alignment methods, the peak detection process is the most time consuming. the alignment of the wine dataset takes totally minutes, consisting of minutes for peak detection, minutes for finding the reference and minutes for the alignment process. the huntington dataset required nearly four times the time for the wine dataset: minutes. most of the time cost is reserved for peak detection. the remainder includes minutes for finding the reference spectrum and minutes for the alignment. besides the parameters for the peak detection, the workflow does not require parameters for the initialization prior to running. one extra optional parameter for clupa, but also for both rafft and icoshift, is the maximum shift point, which is inherited from the fft cross-correlation function. this value is usually ppm, more or less points, depending on the data. if this parameter is not set, the fft cross-correlation function will check on whole spectral regions to get the optimal shift step. furthermore, users need to set the parameters for the peak detection algorithm or use other methods instead. even though default values are set in the software implementation, they are not guaranteed to work optimally for any dataset. the selection of optimal parameters for peak detection algorithm is beyond the scope of this paper.  in this paper, we proposed a simple and efficient workflow which is centered around recursive hierarchical clustering. this strategy, which is referred to as speaq is feasible because prior to the alignment, a peak-picking and reference selection approach is used to reduce the complexity related to alignment. the method aligns the target spectrum to the reference spectrum in a top-down fashion and makes use of fast fourier transformation cross-correlation to reduce the computation time to find the shift step. as conventional recursive strategy like rspa or rafft, the clupa algorithm aligns a larger segment and then recursively divides it into smaller segments to refine the alignment. the differences between them lie in the method and criteria to divide the spectra. rafft selects the dividing point based on the optimal shifts computed from fast fourier transform cross-correlation, rspa detects in the reference and target spectrum lists of segments, then tries to find the corresponding segments between them, and goes into segments for further alignment. our proposed method builds a hierarchical cluster tree from peak lists of reference and target, and divides the spectra into smaller regions based on the most distant clusters of the tree. however, instead of considering many clusters simultaneously, we cut off the dendrogram into just two sub trees and then go further into each of them for finer alignment. another feature of the workflow which is different from most other peak-picking based alignment methods is that instead of grouping peaks together to create segments for alignment, speaq does the alignment first, then it groups peaks and creates smaller segments. in general, this strategy decreases the shift of peaks before the hierarchical cluster tree is built, thus, it significantly reduces overlapping peaks that yield wrong clustering. despite the heuristic approach, selecting the perfect reference remains a non-trivial task, since one spectrum may be a good reference for a specific region but not for other regions. within the software implementation, we therefore also allow users to intervene into the process by setting reference in the region they are confident about. this also allows to skip certain regions during the alignment. the maximum shift can also be set in a region-specific manner. all this information can be readily provided to the software in an extra text file. the speaq workflow was compared to two widely used methods icoshift and rafft. the method shows several advantages over the others. it is easy to use because users do not need to set many initial parameters. experiments on both a public dataset and an huntington dataset show that it performs very well compared to the others. all method use fft cross correlation for finding the shift because of its advantage in computational time. however, speaq is slower than both icoshift and rafft, due to the time for peak detection and its recursive nature, but this is generally not a problem with the powerful computers we have nowadays. in summary, the speaq workflow is unique in the way the subsequent steps are ordered. because of rearranging the order of processing steps, simple and robust methods could be adopted for the workflow. a case in point is the effective analysis based on the bw-ratio. this could only be achieved since nmr spectra are well aligned. if this were not the case, misalignment would introduce sources of variability into the statistic making it difficult to select differential regions. i should be noticed that the proposed method is still open for improvements. for example, from the top panel of figure a it can be seen that the null hypothesis is remarkably constant. this constant line might suggest that the properties of the null distribution are applicable over the entire frequency range. parametrization of the method can avoid the sampling of the null distribution and improve computational efficiency. the method currently does not take into account covariation of the intensities of peaks of the same molecule across samples. the method was tested only on nmr datasets. similar data types, such as chromatography data, may also benefit from the approach presented here. 
the algorithm aligns a target spectrum to the reference spectrum in a top-down fashion by building a hierarchical cluster tree from peak lists of reference and target spectra and then dividing the spectra into smaller segments based on the most distant clusters of the tree. nuclear magnetic resonance spectroscopy is a powerful technique to reveal and compare quantitative metabolic profiles of biological tissues. to reduce the computational time to estimate the spectral misalignment, the method makes use of fast fourier transformation cross-correlation. for example, noise reduction, normalization, baseline correction, peak picking, spectrum alignment and statistical analysis are indispensable components in any nmr analysis pipeline. we introduce a novel suite of informatics tools for the quantitative analysis of nmr metabolomic profile data. correlation maps, spectral and grey scale plots show clear improvements in comparison to other methods, and the down-to-earth quantitative analysis works well for the clupa-aligned spectra. this differential analysis is related to the calculation of the f-statistic or a one-way anova, but without distributional assumptions. the core of the processing cascade is a novel peak alignment algorithm, called hierarchical cluster-based peak alignment. since the method returns a high-quality alignment, we can propose a simple methodology to study the variability of the nmr spectra. for each aligned nmr data point the ratio of the between-group and within-group sum of squares is calculated to quantify the difference in variability between and within predefined groups of nmr spectra. statistical inference based on the bw-ratio is achieved by bootstrapping the null distribution from the experimental data. however, chemical and physical sample variations make the analysis of the data challenging, and typically require the application of a number of preprocessing steps prior to data interpretation.
 nuclear magnetic resonance spectroscopy is a powerful technique to reveal and compare quantitative metabolic profiles of biological tissues. however, chemical and physical sample variations make the analysis of the data challenging, and typically require the application of a number of preprocessing steps prior to data interpretation. for example, noise reduction, normalization, baseline correction, peak picking, spectrum alignment and statistical analysis are indispensable components in any nmr analysis pipeline.  we introduce a novel suite of informatics tools for the quantitative analysis of nmr metabolomic profile data. the core of the processing cascade is a novel peak alignment algorithm, called hierarchical cluster-based peak alignment. the algorithm aligns a target spectrum to the reference spectrum in a top-down fashion by building a hierarchical cluster tree from peak lists of reference and target spectra and then dividing the spectra into smaller segments based on the most distant clusters of the tree. to reduce the computational time to estimate the spectral misalignment, the method makes use of fast fourier transformation cross-correlation. since the method returns a high-quality alignment, we can propose a simple methodology to study the variability of the nmr spectra. for each aligned nmr data point the ratio of the between-group and within-group sum of squares is calculated to quantify the difference in variability between and within predefined groups of nmr spectra. this differential analysis is related to the calculation of the f-statistic or a one-way anova, but without distributional assumptions. statistical inference based on the bw-ratio is achieved by bootstrapping the null distribution from the experimental data.  the workflow performance was evaluated using a previously published dataset. correlation maps, spectral and grey scale plots show clear improvements in comparison to other methods, and the down-to-earth quantitative analysis works well for the clupa-aligned spectra. the whole workflow is embedded into a modular and statistically sound framework that is implemented as an r package called "speaq", which is freely available from http://code.google.com/p/speaq/. 