BACKGROUND
the study of changes in the levels of metabolites and lipids has become essential for the comprehensive understanding of human health
<cit> .
chromatography-coupled mass spectrometry  techniques have become the standard method for characterizing the human metabolome
<cit>  and lipidome
<cit> .
the technique generates a spectrum of peaks describing the sample in the plane defined by the retention time from the chromatograph and the mass-to-charge ratio from the mass spectrometer.
each peak in this plane is either generated by an ion arising from one of the compounds present in the sample, or is an artifact of the measurement without association to any of the compounds.
the association between the peaks and compounds is unknown a priori.
the produced peak data are noisy: first, sample preparation introduces sources of uncertainty that propagate to the analysis
<cit> .
second, the accuracy of the device is limited
<cit>  and it produces biases.
third, peak identification, annotation and pre-processing steps produce additional layers of uncertainty
<cit> .
the effect of errors at all these levels is exacerbated by the “small n, large p” problem: experiments cover a very limited number of samples, n, while the set of compounds measured, p, is potentially large.
however, there also is strong informative structure in the data: first, each compound may generate multiple adduct peaks
<cit>   and isotope peaks
<cit>  , whose positions and shapes provide information about the identity of the compound.
second, the concentrations of different compounds generated by or participating in similar biological processes may be highly correlated
<cit> .
an increasing number of machine learning algorithms are being developed for inferring such structure either from raw spectral data
<cit>  or from processed intensity data
<cit> .
the inference of covariate effects—the differences between sample groups determined by the controlled covariates of the experiment, such as an intervention—is in the core of the comparative analysis of spectral profiles
<cit> .
in addition to the controlled covariates, confounding factors may affect the observations and are subject to the experiment design.
in this work, we focus on inferring effects of the controlled covariates from the data.
<dig>  even though the distinct isotope peaks are not visible to the eye here, they are clearly separable by the mass spectrometer.
in the figure, adduct types and compounds are marked by colors and characters, respectively.
the existence of additional peaks in the spectrum is usually seen as a problem and only the main peak of each identified compound is used for further analysis.
all peaks are a result of the ionisation process where a charged particle is attached to or detached from a compound.
each such compound-ion pair produces a distinct adduct peak.
random variation in the ionisation process leads to inconsistencies between batches of samples, perceived as variation in the ratio of intensities of the peaks associated with one compound.
this is a major source of error for all existing analysis approaches regardless of the choice of the peak used for the analysis.
on the other hand, the distribution of the intensities of isotope peaks is by nature well preserved across both samples and compounds.
moreover, the natural isotopic distribution of a compound is known and can be used to make peak annotation more precise.
in this way, isotope peaks provide reliable additional information about the differences in compound concentrations between sample groups.we propose a probabilistic approach for extending statistical analysis to all available peaks and demonstrate that the additional peaks can provide a real benefit to the inference of covariate effects .
the approach is used to cluster the peaks that are likely to arise from a single compound together and to infer the changes in concentrations of the compounds more accurately based on all these peaks.
by this approach, we are addressing the problem of inadequate sample-size by introducing additional data describing the compounds behind the noisy measurements.
to solve the problem we introduce the following assumptions about the generative process of the data within a bayesian model: first, samples carry between-group differences in their compound concentrations and the differences arise from responses to controlled covariates.
second, multiple observed spectral peaks follow an identical generative process and their heights are a noisy reflection of the true concentration level of the compound.
third, shapes of the peaks from one compound are generated through an identical process following the properties of the measurement device, and thus these shapes are highly similar.
the approach presented in this paper consists of two stages of computational inference:  peaks that share a compound as their generative source are clustered together, and  the responses to controlled covariates of the experiment are inferred on these clusters of peaks.
the clustering part of the approach is based on a nonparametric bayesian dirichlet process model
<cit> .
to improve the performance of this model, we have redefined the prior distributions from a normal distribution to a beta distribution to improve the match to the peak shape similarity observations.
the model for inferring the responses to covariates operates on clusters inferred in the first part.
a bayesian multi-way model
<cit>  is suitable for this task.
this model itself could be used for clustering summarized mass spectral intensity data, but in this work, we demonstrate that the clustering can be done more accurately based upon the similarity of chromatographic peak shapes.
methods
this section describing the models consists of two parts: clustering of spectral peaks and inference of covariate effects.
to maintain the mathematical rigor in the section, we use the terms “samples,” “variables” and “clusters” to refer to the experimental runs of the mass spectrometer, the peaks in the mass spectrometry data, and the yet unknown compounds in the experimental runs, respectively.
in the equations, we denote them by the indices
i= <dig> …,n,j= <dig> …,p,k= <dig> …,k,
respectively, where n, p and k are their respective total numbers.
we use bold capital, bold non-capital and non-bold non-capital symbols to refer to matrices, vectors and scalars, respectively .
clusters of peaks based on the similarity
following earlier work
<cit> , we measure the similarity between the shapes of two peaks by their pearson correlation computed over a window of retention time after a standard peak alignment
<cit>  across the samples.
truncating negative values to zero, this leads to a distinct similarity matrix q
i··∈ <cit>
p×p
for each sample i. in the notation, the operator “ ·” indicates that the entire dimension of the array is included, not only the single item j. since a peak is not necessarily observed in every sample, there may be missing values in the matrices.
therefore, we construct an additional mask r∈{ <dig> }
n×p×p
with binary values
rijj′ indicating whether the peak pair  in sample i appears together within the window where the similarity is measured and whether both of the peaks are observed.
an unidentified peak may still be present in the sample below the limit of detection of the mass spectrometer.
however, then it is not useful for the inference of covariate effects and, thus, is treated as missing.
model
we assume that the peaks are generated through a dirichlet process
<cit> : there is an unknown number of clusters and an unknown and variable number of peaks that arise from each of the clusters.
peaks are assumed to have a one-out-of-many association: each peak is associated with only one of the unknown clusters.
with these basic assumptions, we can infer the p-by-k clustering matrix v from the data q. value v
j
k
= <dig> in the clustering matrix v assigns the peak j to the cluster k. to make the following equations more compact, we use an additional variable,
εjj′=vj·vj′·t∈{ <dig> }, which is an inner product of the cluster indicator vectors of the peaks j and j′, to denote whether the two peaks come from the same or different clusters .
we set a spike-and-slab prior
<cit>  for the peak shape similarity to model the inherent sparse structure of the data.
the similarity between any pair of observed peaks  is assumed to follow a beta distribution, but the shape of the distribution is assumed to depend on whether the pair comes from the same cluster or from different clusters  or , when
εjj′= <dig> or  <dig>  respectively).
also the probability of a missing similarity value is assumed to depend on the cluster assignment of the pair .
the distributional assumptions are
qijj′|εjj′∼rijj′1-p0inbetaqijj′|ain,bin+p0inδrijj′,εjj′= <dig> rijj′1-p0outbetaqijj′|aout,bout+p0outδrijj′,εjj′= <dig>
with the first and the second row of the equation stating the distributions of a peak pair from the same cluster and different clusters, respectively.
the likelihood of the entire peak shape data,
ℒq,r|v=∏i=1n∏j=1p-1∏j′=j+1ppqijj′,rijj′|εjj′,
becomes a product over all peak pairs and samples following the distributional assumption of equation  <dig>
we further assume that the observed peaks are generated from an unknown finite subset of an infinite set of clusters with an equal prior probability,
pεjj′=1=1p-1+αdp,
for any pair of peaks to be generated from the same cluster.
these assumptions define the dirichlet process, controlled by the concentration parameter αdp, which determines the prior probability mass outside the existing clusters.
following from this prior assumption, the probability of assigning peak j to an existing cluster k,
pvjk=1|q,r,v-j,·∝skℒq,r|v-j,·,vjk= <dig>
becomes weighted by the current size of the cluster,
sk=v-j,ktv-j,k.
in the notation, matrices v·,-k
and v-j,· correspond to the matrix v with the column k and the row j omitted, respectively.
alternatively, with probability
pvj,k+1=1|q,r,v∝αdpℒq,r|v-j,·,vj,k+1= <dig>
the process may create a new cluster with the index k+ <dig> and only the peak j inside.
then, the likelihood term is weighted by the dirichlet process concentration parameter αdp, which can be seen as a pseudo-count for the number of peaks outside the current k clusters.
inference
we infer the posterior distribution of the clustering via gibbs sampling, which results in a set of s samples of the clustering v, s= <dig> …,s, from the true posterior distribution p. the computational complexity of a gibbs iteration is
okp <dig>  further analysis can operate on the entire posterior distribution of the clustering through integration, or on a point estimate of the distribution.
we follow earlier work
<cit>  and acquire a point estimate of the posterior distribution of the clustering through finding the least-squares clustering .
covariate effects based on peak heights
having inferred the grouping of similar peaks into clusters that each correspond to a compound, we infer the differences in concentrations between sample groups for each cluster given the peak height data
x∈rp×n and the clustering v. again, some values in the data may be missing.
model
after a peak-specific centering based on the control group, the observed peak heights for each sample i are assumed to be normally distributed with a cluster-specific mean
x·ilat:
x·i|v,xlat,σ2∼nvx·ilat,Λ,
where the diagonal matrix
Λ
contains the peak-specific variance parameters
σ2∈r+p.
the cluster-specific means are assumed to be normally distributed with a sample group-specific prior
α
,
x·ilat|α,ai∼nα·ai,i,
where a
i
∈{ <dig> …,l
a
} is an indicator of group membership  for sample i and i is a k-by-k identity matrix.
the corresponding covariate effects are arranged into an k-by- l
a
matrix
α
and the effects are assumed to be independent and normally distributed,
α·l∼δα·l,l=1n <dig> i,l= <dig> …,la,
except for the first level, l= <dig>  which is defined as the baseline level and thus is fixed to zero.
the model is not limited to one covariate: the cluster-specific mean
xi·lat can be expressed as a sum of effects of multiple covariates and their interaction effects .
further, the model is readily extensible for dependent covariate effects
<cit> .
the peak-specific variance parameter,
σj2∼scale-inv-χ2n <dig> σ <dig>
follows a scaled inverse- χ <dig> distribution with n <dig> prior samples and a scale
σ <dig>
inference and analysis
we infer the covariate effects via gibbs sampling.
now the clustering matrix v has been learned earlier, and is thus taken as known in the model.
computational complexity of a gibbs iteration is
onpk <dig>  the clustering and the covariate effects can be inferred overnight on a standard desktop computer for a typical-sized data set.
the posterior distributions of the covariate effects
α
are descriptive of the differences between the sample groups and, thus, interesting from the analysis point of view.
to assess the significance of the difference between a sample group, c=l> <dig>  and the control group, c= <dig>  for a cluster k, we can study the posterior probability of the effect α
k
l
being greater or less than zero.
comparison methods
we call the method described above model  <dig>  we compared the performance of the following approaches and refer to them as models  <dig>   <dig> and 3:
<dig>  the multi-peak approach using both peak shape and height information, as proposed in this work ,
<dig>  the multi-peak approach using peak height information only
<cit>  ,
<dig>  the typical single-peak approach .
for the studied real data sets, we discovered that peak height information alone is not enough for clustering the peaks into individual compounds due to the high level of noise and strong correlations between compounds.
thus, for real data we compared model  <dig> to model  <dig> and highlight the benefit gained by using peak shape information.
model  <dig> assumes the generative gaussian latent variable model of the equations 7– <dig> for the intensity observations x and a uniform multinomial prior for the clustering of the peaks.
the clustering is inferred by gibbs sampling together with the covariate effects.
model  <dig> quantifies the difference between the covariate level, c=l, and the control level, c= <dig>  as the difference of their means based on the main peak j,
αj,l=1∑i=1nδai,l∑i=1nδai,lxj,i-1∑i=1nδai,1∑i=1nδai,1xj,i.
the kronecker delta function
δai,l selects the samples that have the covariate level l by receiving the value  <dig>  when a
i
=l, and  <dig>  otherwise.
when the data are log-transformed, the mean difference corresponds to the fold change computed in many analysis platforms such as mzmine
<cit>  and xcms
<cit> .
experiments
we demonstrate the performance of the proposed method through three experiments: a simulated data experiment, a spike-in benchmark experiment with known changes in concentrations, and a gene silencing experiment with measurements of the lipidome of cancer cells.
evaluation measures
evaluation of the performance on real data sets is not a trivial task, as there is no ground truth available: neither the identity of the peaks nor the true effect sizes are known.
thus, we also used spike-in data, where the true covariate effects are known, although only a small number of the peaks are annotated.
for the simulated and benchmark experiments, we computed the mean squared error  between inferred and true covariate effects as an evaluation metric.
as a result of the log-transformation of the intensity data, we were quantifying relative changes between sample groups, independent of the average height of each peak.
in the model, we thus assumed that the change is preserved across the peaks of one compound, in relative terms.
the significance of the difference in the mse of the proposed approach and the comparison method was tested by the paired one-sided t-test.
the false discovery rate was controlled by the benjamini-hochberg step-up procedure
<cit> .
additionally for the simulated experiment, we studied the inference of the statistical significance of effects, since the true distribution of the data was known.
to assess the sensitivity of the approaches to noise in natural lipidomic data lacking a ground truth, we used two types of indirect evaluation: first, we studied the consistency of the inferred covariate effects given a prior assumption about their similarity.
second, we examined the robustness of the inferred covariate effects to noise.
finally, we demonstrated differences between the multi-peak and single-peak approaches through examples of qualitative analysis of annotated peak clusters.
simulated data
we started by investigating the performance of the proposed approach on synthetic data, where the true covariate effects are known.
we focused on a usual task in exploratory analysis of biological data: the detection of significant non-zero covariate effects.
we measured the performance by the accuracy at this task—the ratio of true positive and true negative significant differences among all effects.
we used the 95% posterior quantiles to determine the significance.
additionally, we compared the approaches by the mse to the true effects and studied the performance of the two clustering models by computing the normalized information distance
<cit>  to the true clustering.
the approaches were tested across a set of potential experimental settings to study how the observation of additional peaks and samples affects the performance.
simulated data were generated by assuming the latent structure of model  <dig>  the following data parameters were varied on a grid: sample-size n=2×{ <dig> ,15} and peak-specific noise σ2={ <dig> }.
additionally, the number of peaks per cluster was varied between  <dig>   <dig> and  <dig>  covariate effects
α
·2= were generated for each data set.
the experiment was repeated  <dig> times with independent data sets.
the results are reported in the results and discussion section.
benchmark data with known changes in concentrations
the benchmark data set of apple samples
<cit>  includes a set of annotated spike-in compounds with increases of  <dig>   <dig> or 100% in concentrations.
we started with the raw spectral data
<cit>  in order to acquire the shapes of the peaks in addition to their heights.
the mass spectra were pre-processed using mzmine 2
<cit>  .
we used standard pre-processing methodology also used in the original publications of the data sets, thus maintaining the focus of the work on the potential benefit gained from the multiple peaks.
the compared approaches were on the same line in terms of the data.
we evaluated the approaches by the mse between inferred and true covariate effects.
if the cluster contained multiple annotated peaks, the effect of each annotated peak was evaluated separately for the single-peak approach.
clusters with no annotated peaks were considered to have a 0% true effect and the effect of the single-peak approach was inferred based on the strongest peak of the cluster.
lipidomic data from a gene silencing study
the data come from a recent experiment
<cit>  to study the effects of gene silencing treatments on lipidomic profiles and growth of breast cancer tissue.
driven by the lack of ground truth about the covariate effects, we evaluated the inferred effects indirectly in two ways:  by quantifying the consistency of the effects within a lipid family and  by quantifying the robustness of the magnitudes of the inferred effects across the lipidome in presence of additional noise.
additionally, we investigated the stability of the inferred clustering on the data and qualitatively analyzed differences between the covariate effects of single peaks and the effects inferred on clusters of peaks by model  <dig>
the data included  <dig> lipidomic profiles of breast cancer cells from the zr-75- <dig> cell line.
we inferred the effects of seven distinct silencing interventions on metabolism- regulating genes  at two time points.
the raw spectra were pre-processed with mzmine  <dig> as described in the original publication
<cit> , in addition to which the shape similarities of the peaks were computed.
consistency of effect signs
in the first task, we quantified the consistency as the accuracy at predicting the covariate effect of a test lipid given the model on the covariate effects of other lipids of the same family.
for instance, we predicted the effect of a gene silencing treatment on the sphingomyelin sm based on the sphingomyelin compounds in the training set.
we examined the sign of the effect instead of the absolute effect, since even within a family of lipids the changes have a high variance and thus cannot be reliably predicted without imposing additional information about the biological system.
we predicted the signs of the covariate effects for test lipids in a three-fold cross-validation setting with  <dig> randomizations.
the examined lipids included the annotated members from the three most abundant families of lipids that had two or more peaks identified with the clustering model .
further, we studied the influence of noise to the consistency by adding independent normally distributed noise  on the peak intensity observations.
added noise variance σ= <dig> was equal to the existing original variance in the data, and the upper bound for the signal-to-noise ratio then was  <dig>  .
robustness of effect magnitudes
to evaluate the inferred effects at the scale of the entire observed lipidome, we examined the consistency of inferred covariate effects between the original and noise-added data sets.
this experiment simulated the situation where the true effects are known , but the data based on which the effects are inferred are noisy .
to measure the consistency, we computed the spearman correlation between the covariate effects inferred from the original and the added-noise data sets.
we studied all clusters with two or more peaks, constituting 20% of the clusters.
RESULTS
simulated data
on a normal level of noise , the multi-peak approaches  always performed better at detecting significant covariate effects than the single-peak approach  and only with enough samples the performance of model  <dig> became comparable to model  <dig>  the inferred clustering of model  <dig> was perfect while the clustering performance of model  <dig> heavily depended on the number of samples available .
on a high level of noise , only model  <dig> worked .
the reason for the failure of model  <dig> was the imperfectly inferred clustering .
the good performance of model  <dig> resulted from the clustering step, which is robust to noise in the peak heights.
the peak shape similarity gave strong evidence for inferring the clusters already from a single sample.
the mse between the inferred and true covariate effects for model  <dig> was smaller compared to model  <dig> in all the  <dig> setups of the experimental grid .
the difference was statistically significant in  <dig> setups and in all setups at the high level of noise.the performance of model  <dig> clearly improved, when more peaks from a cluster were present in the data .
this was pronounced at a high level of noise, when the observation of a single peak is unreliable for inferring the covariate effects.
in a similar way as in averaging over samples, the model is able to overcome peak-specific noise also by averaging over multiple peaks.
benchmark data with known changes in concentrations
in the first demonstration on real uplc-ms data
<cit> , we show that model  <dig> can infer the artificial perturbations in a spike-in experiment more accurately than the single-peak approach.
in the positive ion mode, the model inferred  <dig> clusters, among which  <dig> clusters included more than one peak.
seven clusters included annotated peaks from the spike-in compounds, four of which included more than one annotated peak .
peaks from two compounds were distributed to two and four clusters, respectively.
in the negative ion mode, the model inferred  <dig> clusters, among which  <dig> clusters were non-singletons.
three clusters included annotated peaks from the spike-in compounds, all of these clusters included more than one annotated peak and all peaks from one compound were clustered together.
in both the ion modes, all clusters with annotated peaks were specific to one compound.
model  <dig> had a lower error than model  <dig> at all magnitudes of the true effect with the strongest relative improvement occurring at the small magnitudes .
the difference was statistically significant for covariate effects from  <dig> to 40% .
lipidomic data from a gene silencing study
in the second demonstration on real uplc-ms data
<cit> , we show that model  <dig> can infer more consistent covariate effects in two ways even though the true effects are unknown.
consistency and robustness of effects
when examining the consistency of effects within a lipid family, model  <dig> was more consistent than model  <dig> at all levels of noise .
when no noise was added and also at moderate levels of noise, both approaches performed clearly better than expected by random chance.
when noise was added, model  <dig> suffered more and its performance reduced to the random level more rapidly.
given the assumption about the general similarity of lipids within a family is true, model  <dig> inferred the covariate effects more consistently.when examining the robustness of effect magnitudes, model  <dig> was more consistent than model  <dig> when noise was added to the data .
the confidence intervals from the  <dig> randomizations did not overlap at all at moderate levels of noise.
stability
since the proposed approach is sensitive to the inferred clustering of the data, we analyzed the stability of the inferred clustering on biological data, using the lipidomic gene silencing data as a case study.
we tested the influence of the concentration parameter αdp in the dirichlet process clustering model.
the clustering result for the lipidomic gene silencing data was robust to changes in the magnitude of the concentration parameter .
as expected, the number of clusters increased, when the preset value of the concentration parameter increased, but the relative change was small.
qualitative analysis
finally, we give concrete examples of potential findings that the approaches can uncover and demonstrate how analysis based on a single peak may lead to a different outcome depending on the choice of the peak.the intervention-driven changes of individual peaks from two lipids along with the covariate effects inferred by models  <dig> and  <dig> are shown in figure
<dig>  in the case of the sphingomyelin sm, there were strong covariate effects inferred by model  <dig> but many of these effects became weaker when inferred based on multiple peaks by model  <dig>  on the contrary, model  <dig> inferred weak covariate effects for the ceramide cer but based on multiple peaks and model  <dig>  one of the effects was actually among the top-5% strongest effects across the observed lipidome.
CONCLUSIONS
we have empirically demonstrated that a model-based integration of multiple peaks can lead to an improved accuracy in the inference of covariate effects, and we introduced a novel method for this task.
while the sample-size is always restricted by external constraints such as the experiment budget or the availability of suitable patients, the inference based on multiple peaks gives a shortcut to extracting more information from the limited set of samples, thereby directly addressing the “small n, large p” problem.
however, some types of systematic measurement error, such as some batch effects, escape this treatment and can only be reduced by introducing independent replicates.
based on the results presented in this work, we argue that additional peaks are especially useful when the signal-to-noise ratio is low and the differences between sample groups are small.
we suggest that all the detected peaks that can be associated with a compound should be taken into account in the comparative analysis.
this is possible through the two-step generative modeling approach presented in this work:  by identifying the peaks that can be associated with one compound through clustering the peaks based on their shape similarity and  by the inference of covariate effects on the clusters, each representing one compound.
abbreviations
anova: analysis of variance; cer: ceramide; sm: sphingomyelin; uplc-ms: ultra performance liquid chromatography-mass spectrometry.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
the method was developed jointly by ts, sk and sr. ts had a lead role at implementing the model, designing and implementing the experiments, and at preparing the manuscript.
all authors read and approved the manuscript.
supplementary material
additional file 1
supplementary material.
more details of the experiments.
click here for file
acknowledgements
the authors would like to thank sandra castillo, peddinti v. gopalacharyulu, mika hilvo and matej orešič for providing data and for useful discussions.
the authors would also like to thank rónán daly and joe wandy for useful discussions.
this work was supported by the academy of finland , the finnish foundation for technology promotion  and the finnish doctoral programme in computational sciences fics .
the calculations presented in the work were performed using computer resources within the aalto university school of science "science-it" project.

 the study of changes in the levels of metabolites and lipids has become essential for the comprehensive understanding of human health. chromatography-coupled mass spectrometry techniques have become the standard method for characterizing the human metabolome and lipidome. the technique generates a spectrum of peaks describing the sample in the plane defined by the retention time from the chromatograph and the mass-to-charge ratio from the mass spectrometer. each peak in this plane is either generated by an ion arising from one of the compounds present in the sample, or is an artifact of the measurement without association to any of the compounds. the association between the peaks and compounds is unknown a priori. the produced peak data are noisy: first, sample preparation introduces sources of uncertainty that propagate to the analysis. second, the accuracy of the device is limited and it produces biases. third, peak identification, annotation and pre-processing steps produce additional layers of uncertainty. the effect of errors at all these levels is exacerbated by the “small n, large p” problem: experiments cover a very limited number of samples, n, while the set of compounds measured, p, is potentially large. however, there also is strong informative structure in the data: first, each compound may generate multiple adduct peaks and isotope peaks, whose positions and shapes provide information about the identity of the compound. second, the concentrations of different compounds generated by or participating in similar biological processes may be highly correlated. an increasing number of machine learning algorithms are being developed for inferring such structure either from raw spectral data or from processed intensity data. the inference of covariate effects—the differences between sample groups determined by the controlled covariates of the experiment, such as an intervention—is in the core of the comparative analysis of spectral profiles. in addition to the controlled covariates, confounding factors may affect the observations and are subject to the experiment design. in this work, we focus on inferring effects of the controlled covariates from the data. even though the distinct isotope peaks are not visible to the eye here, they are clearly separable by the mass spectrometer. in the figure, adduct types and compounds are marked by colors and characters, respectively. the existence of additional peaks in the spectrum is usually seen as a problem and only the main peak of each identified compound is used for further analysis. all peaks are a result of the ionisation process where a charged particle is attached to or detached from a compound. each such compound-ion pair produces a distinct adduct peak. random variation in the ionisation process leads to inconsistencies between batches of samples, perceived as variation in the ratio of intensities of the peaks associated with one compound. this is a major source of error for all existing analysis approaches regardless of the choice of the peak used for the analysis. on the other hand, the distribution of the intensities of isotope peaks is by nature well preserved across both samples and compounds. moreover, the natural isotopic distribution of a compound is known and can be used to make peak annotation more precise. in this way, isotope peaks provide reliable additional information about the differences in compound concentrations between sample groups.we propose a probabilistic approach for extending statistical analysis to all available peaks and demonstrate that the additional peaks can provide a real benefit to the inference of covariate effects. the approach is used to cluster the peaks that are likely to arise from a single compound together and to infer the changes in concentrations of the compounds more accurately based on all these peaks. by this approach, we are addressing the problem of inadequate sample-size by introducing additional data describing the compounds behind the noisy measurements. to solve the problem we introduce the following assumptions about the generative process of the data within a bayesian model: first, samples carry between-group differences in their compound concentrations and the differences arise from responses to controlled covariates. second, multiple observed spectral peaks follow an identical generative process and their heights are a noisy reflection of the true concentration level of the compound. third, shapes of the peaks from one compound are generated through an identical process following the properties of the measurement device, and thus these shapes are highly similar. the approach presented in this paper consists of two stages of computational inference: peaks that share a compound as their generative source are clustered together, and the responses to controlled covariates of the experiment are inferred on these clusters of peaks. the clustering part of the approach is based on a nonparametric bayesian dirichlet process model. to improve the performance of this model, we have redefined the prior distributions from a normal distribution to a beta distribution to improve the match to the peak shape similarity observations. the model for inferring the responses to covariates operates on clusters inferred in the first part. a bayesian multi-way model is suitable for this task. this model itself could be used for clustering summarized mass spectral intensity data, but in this work, we demonstrate that the clustering can be done more accurately based upon the similarity of chromatographic peak shapes. methods this section describing the models consists of two parts: clustering of spectral peaks and inference of covariate effects. to maintain the mathematical rigor in the section, we use the terms “samples,” “variables” and “clusters” to refer to the experimental runs of the mass spectrometer, the peaks in the mass spectrometry data, and the yet unknown compounds in the experimental runs, respectively. in the equations, we denote them by the indices i= …,n,j= …,p,k= …,k, respectively, where n, p and k are their respective total numbers. we use bold capital, bold non-capital and non-bold non-capital symbols to refer to matrices, vectors and scalars, respectively. clusters of peaks based on the similarity following earlier work, we measure the similarity between the shapes of two peaks by their pearson correlation computed over a window of retention time after a standard peak alignment across the samples. truncating negative values to zero, this leads to a distinct similarity matrix q i··∈ p×p for each sample i. in the notation, the operator “ ·” indicates that the entire dimension of the array is included, not only the single item j. since a peak is not necessarily observed in every sample, there may be missing values in the matrices. therefore, we construct an additional mask r∈{ } n×p×p with binary values rijj′ indicating whether the peak pair in sample i appears together within the window where the similarity is measured and whether both of the peaks are observed. an unidentified peak may still be present in the sample below the limit of detection of the mass spectrometer. however, then it is not useful for the inference of covariate effects and, thus, is treated as missing. model we assume that the peaks are generated through a dirichlet process: there is an unknown number of clusters and an unknown and variable number of peaks that arise from each of the clusters. peaks are assumed to have a one-out-of-many association: each peak is associated with only one of the unknown clusters. with these basic assumptions, we can infer the p-by-k clustering matrix v from the data q. value v j k = in the clustering matrix v assigns the peak j to the cluster k. to make the following equations more compact, we use an additional variable, εjj′=vj·vj′·t∈{ }, which is an inner product of the cluster indicator vectors of the peaks j and j′, to denote whether the two peaks come from the same or different clusters. we set a spike-and-slab prior for the peak shape similarity to model the inherent sparse structure of the data. the similarity between any pair of observed peaks is assumed to follow a beta distribution, but the shape of the distribution is assumed to depend on whether the pair comes from the same cluster or from different clusters or, when εjj′= or respectively). also the probability of a missing similarity value is assumed to depend on the cluster assignment of the pair. the distributional assumptions are qijj′|εjj′∼rijj′1-p0inbetaqijj′|ain,bin+p0inδrijj′,εjj′= rijj′1-p0outbetaqijj′|aout,bout+p0outδrijj′,εjj′= with the first and the second row of the equation stating the distributions of a peak pair from the same cluster and different clusters, respectively. the likelihood of the entire peak shape data, ℒq,r|v=∏i=1n∏j=1p-1∏j′=j+1ppqijj′,rijj′|εjj′, becomes a product over all peak pairs and samples following the distributional assumption of equation we further assume that the observed peaks are generated from an unknown finite subset of an infinite set of clusters with an equal prior probability, pεjj′=1=1p-1+αdp, for any pair of peaks to be generated from the same cluster. these assumptions define the dirichlet process, controlled by the concentration parameter αdp, which determines the prior probability mass outside the existing clusters. following from this prior assumption, the probability of assigning peak j to an existing cluster k, pvjk=1|q,r,v-j,·∝skℒq,r|v-j,·,vjk= becomes weighted by the current size of the cluster, sk=v-j,ktv-j,k. in the notation, matrices v·,-k and v-j,· correspond to the matrix v with the column k and the row j omitted, respectively. alternatively, with probability pvj,k+1=1|q,r,v∝αdpℒq,r|v-j,·,vj,k+1= the process may create a new cluster with the index k+ and only the peak j inside. then, the likelihood term is weighted by the dirichlet process concentration parameter αdp, which can be seen as a pseudo-count for the number of peaks outside the current k clusters. inference we infer the posterior distribution of the clustering via gibbs sampling, which results in a set of s samples of the clustering v, s= …,s, from the true posterior distribution p. the computational complexity of a gibbs iteration is okp further analysis can operate on the entire posterior distribution of the clustering through integration, or on a point estimate of the distribution. we follow earlier work and acquire a point estimate of the posterior distribution of the clustering through finding the least-squares clustering. covariate effects based on peak heights having inferred the grouping of similar peaks into clusters that each correspond to a compound, we infer the differences in concentrations between sample groups for each cluster given the peak height data x∈rp×n and the clustering v. again, some values in the data may be missing. model after a peak-specific centering based on the control group, the observed peak heights for each sample i are assumed to be normally distributed with a cluster-specific mean x·ilat: x·i|v,xlat,σ2∼nvx·ilat,Λ, where the diagonal matrix Λ contains the peak-specific variance parameters σ2∈r+p. the cluster-specific means are assumed to be normally distributed with a sample group-specific prior α, x·ilat|α,ai∼nα·ai,i, where a i ∈{ …,l a } is an indicator of group membership for sample i and i is a k-by-k identity matrix. the corresponding covariate effects are arranged into an k-by- l a matrix α and the effects are assumed to be independent and normally distributed, α·l∼δα·l,l=1n i,l= …,la, except for the first level, l= which is defined as the baseline level and thus is fixed to zero. the model is not limited to one covariate: the cluster-specific mean xi·lat can be expressed as a sum of effects of multiple covariates and their interaction effects. further, the model is readily extensible for dependent covariate effects. the peak-specific variance parameter, σj2∼scale-inv-χ2n σ follows a scaled inverse- χ distribution with n prior samples and a scale σ inference and analysis we infer the covariate effects via gibbs sampling. now the clustering matrix v has been learned earlier, and is thus taken as known in the model. computational complexity of a gibbs iteration is onpk the clustering and the covariate effects can be inferred overnight on a standard desktop computer for a typical-sized data set. the posterior distributions of the covariate effects α are descriptive of the differences between the sample groups and, thus, interesting from the analysis point of view. to assess the significance of the difference between a sample group, c=l> and the control group, c= for a cluster k, we can study the posterior probability of the effect α k l being greater or less than zero. comparison methods we call the method described above model we compared the performance of the following approaches and refer to them as models and 3: the multi-peak approach using both peak shape and height information, as proposed in this work, the multi-peak approach using peak height information only, the typical single-peak approach. for the studied real data sets, we discovered that peak height information alone is not enough for clustering the peaks into individual compounds due to the high level of noise and strong correlations between compounds. thus, for real data we compared model to model and highlight the benefit gained by using peak shape information. model assumes the generative gaussian latent variable model of the equations 7– for the intensity observations x and a uniform multinomial prior for the clustering of the peaks. the clustering is inferred by gibbs sampling together with the covariate effects. model quantifies the difference between the covariate level, c=l, and the control level, c= as the difference of their means based on the main peak j, αj,l=1∑i=1nδai,l∑i=1nδai,lxj,i-1∑i=1nδai,1∑i=1nδai,1xj,i. the kronecker delta function δai,l selects the samples that have the covariate level l by receiving the value when a i =l, and otherwise. when the data are log-transformed, the mean difference corresponds to the fold change computed in many analysis platforms such as mzmine and xcms. experiments we demonstrate the performance of the proposed method through three experiments: a simulated data experiment, a spike-in benchmark experiment with known changes in concentrations, and a gene silencing experiment with measurements of the lipidome of cancer cells. evaluation measures evaluation of the performance on real data sets is not a trivial task, as there is no ground truth available: neither the identity of the peaks nor the true effect sizes are known. thus, we also used spike-in data, where the true covariate effects are known, although only a small number of the peaks are annotated. for the simulated and benchmark experiments, we computed the mean squared error between inferred and true covariate effects as an evaluation metric. as a result of the log-transformation of the intensity data, we were quantifying relative changes between sample groups, independent of the average height of each peak. in the model, we thus assumed that the change is preserved across the peaks of one compound, in relative terms. the significance of the difference in the mse of the proposed approach and the comparison method was tested by the paired one-sided t-test. the false discovery rate was controlled by the benjamini-hochberg step-up procedure. additionally for the simulated experiment, we studied the inference of the statistical significance of effects, since the true distribution of the data was known. to assess the sensitivity of the approaches to noise in natural lipidomic data lacking a ground truth, we used two types of indirect evaluation: first, we studied the consistency of the inferred covariate effects given a prior assumption about their similarity. second, we examined the robustness of the inferred covariate effects to noise. finally, we demonstrated differences between the multi-peak and single-peak approaches through examples of qualitative analysis of annotated peak clusters. simulated data we started by investigating the performance of the proposed approach on synthetic data, where the true covariate effects are known. we focused on a usual task in exploratory analysis of biological data: the detection of significant non-zero covariate effects. we measured the performance by the accuracy at this task—the ratio of true positive and true negative significant differences among all effects. we used the 95% posterior quantiles to determine the significance. additionally, we compared the approaches by the mse to the true effects and studied the performance of the two clustering models by computing the normalized information distance to the true clustering. the approaches were tested across a set of potential experimental settings to study how the observation of additional peaks and samples affects the performance. simulated data were generated by assuming the latent structure of model the following data parameters were varied on a grid: sample-size n=2×{,15} and peak-specific noise σ2={ }. additionally, the number of peaks per cluster was varied between and covariate effects α ·2= were generated for each data set. the experiment was repeated times with independent data sets. the results are reported in the results and discussion section. benchmark data with known changes in concentrations the benchmark data set of apple samples includes a set of annotated spike-in compounds with increases of or 100% in concentrations. we started with the raw spectral data in order to acquire the shapes of the peaks in addition to their heights. the mass spectra were pre-processed using mzmine 2. we used standard pre-processing methodology also used in the original publications of the data sets, thus maintaining the focus of the work on the potential benefit gained from the multiple peaks. the compared approaches were on the same line in terms of the data. we evaluated the approaches by the mse between inferred and true covariate effects. if the cluster contained multiple annotated peaks, the effect of each annotated peak was evaluated separately for the single-peak approach. clusters with no annotated peaks were considered to have a 0% true effect and the effect of the single-peak approach was inferred based on the strongest peak of the cluster. lipidomic data from a gene silencing study the data come from a recent experiment to study the effects of gene silencing treatments on lipidomic profiles and growth of breast cancer tissue. driven by the lack of ground truth about the covariate effects, we evaluated the inferred effects indirectly in two ways: by quantifying the consistency of the effects within a lipid family and by quantifying the robustness of the magnitudes of the inferred effects across the lipidome in presence of additional noise. additionally, we investigated the stability of the inferred clustering on the data and qualitatively analyzed differences between the covariate effects of single peaks and the effects inferred on clusters of peaks by model the data included lipidomic profiles of breast cancer cells from the zr-75- cell line. we inferred the effects of seven distinct silencing interventions on metabolism- regulating genes at two time points. the raw spectra were pre-processed with mzmine as described in the original publication, in addition to which the shape similarities of the peaks were computed. consistency of effect signs in the first task, we quantified the consistency as the accuracy at predicting the covariate effect of a test lipid given the model on the covariate effects of other lipids of the same family. for instance, we predicted the effect of a gene silencing treatment on the sphingomyelin sm based on the sphingomyelin compounds in the training set. we examined the sign of the effect instead of the absolute effect, since even within a family of lipids the changes have a high variance and thus cannot be reliably predicted without imposing additional information about the biological system. we predicted the signs of the covariate effects for test lipids in a three-fold cross-validation setting with randomizations. the examined lipids included the annotated members from the three most abundant families of lipids that had two or more peaks identified with the clustering model. further, we studied the influence of noise to the consistency by adding independent normally distributed noise on the peak intensity observations. added noise variance σ= was equal to the existing original variance in the data, and the upper bound for the signal-to-noise ratio then was. robustness of effect magnitudes to evaluate the inferred effects at the scale of the entire observed lipidome, we examined the consistency of inferred covariate effects between the original and noise-added data sets. this experiment simulated the situation where the true effects are known, but the data based on which the effects are inferred are noisy. to measure the consistency, we computed the spearman correlation between the covariate effects inferred from the original and the added-noise data sets. we studied all clusters with two or more peaks, constituting 20% of the clusters.  simulated data on a normal level of noise, the multi-peak approaches always performed better at detecting significant covariate effects than the single-peak approach and only with enough samples the performance of model became comparable to model the inferred clustering of model was perfect while the clustering performance of model heavily depended on the number of samples available. on a high level of noise, only model worked. the reason for the failure of model was the imperfectly inferred clustering. the good performance of model resulted from the clustering step, which is robust to noise in the peak heights. the peak shape similarity gave strong evidence for inferring the clusters already from a single sample. the mse between the inferred and true covariate effects for model was smaller compared to model in all the setups of the experimental grid. the difference was statistically significant in setups and in all setups at the high level of noise.the performance of model clearly improved, when more peaks from a cluster were present in the data. this was pronounced at a high level of noise, when the observation of a single peak is unreliable for inferring the covariate effects. in a similar way as in averaging over samples, the model is able to overcome peak-specific noise also by averaging over multiple peaks. benchmark data with known changes in concentrations in the first demonstration on real uplc-ms data, we show that model can infer the artificial perturbations in a spike-in experiment more accurately than the single-peak approach. in the positive ion mode, the model inferred clusters, among which clusters included more than one peak. seven clusters included annotated peaks from the spike-in compounds, four of which included more than one annotated peak. peaks from two compounds were distributed to two and four clusters, respectively. in the negative ion mode, the model inferred clusters, among which clusters were non-singletons. three clusters included annotated peaks from the spike-in compounds, all of these clusters included more than one annotated peak and all peaks from one compound were clustered together. in both the ion modes, all clusters with annotated peaks were specific to one compound. model had a lower error than model at all magnitudes of the true effect with the strongest relative improvement occurring at the small magnitudes. the difference was statistically significant for covariate effects from to 40%. lipidomic data from a gene silencing study in the second demonstration on real uplc-ms data, we show that model can infer more consistent covariate effects in two ways even though the true effects are unknown. consistency and robustness of effects when examining the consistency of effects within a lipid family, model was more consistent than model at all levels of noise. when no noise was added and also at moderate levels of noise, both approaches performed clearly better than expected by random chance. when noise was added, model suffered more and its performance reduced to the random level more rapidly. given the assumption about the general similarity of lipids within a family is true, model inferred the covariate effects more consistently.when examining the robustness of effect magnitudes, model was more consistent than model when noise was added to the data. the confidence intervals from the randomizations did not overlap at all at moderate levels of noise. stability since the proposed approach is sensitive to the inferred clustering of the data, we analyzed the stability of the inferred clustering on biological data, using the lipidomic gene silencing data as a case study. we tested the influence of the concentration parameter αdp in the dirichlet process clustering model. the clustering result for the lipidomic gene silencing data was robust to changes in the magnitude of the concentration parameter. as expected, the number of clusters increased, when the preset value of the concentration parameter increased, but the relative change was small. qualitative analysis finally, we give concrete examples of potential findings that the approaches can uncover and demonstrate how analysis based on a single peak may lead to a different outcome depending on the choice of the peak.the intervention-driven changes of individual peaks from two lipids along with the covariate effects inferred by models and are shown in figure in the case of the sphingomyelin sm, there were strong covariate effects inferred by model but many of these effects became weaker when inferred based on multiple peaks by model on the contrary, model inferred weak covariate effects for the ceramide cer but based on multiple peaks and model one of the effects was actually among the top-5% strongest effects across the observed lipidome.  we have empirically demonstrated that a model-based integration of multiple peaks can lead to an improved accuracy in the inference of covariate effects, and we introduced a novel method for this task. while the sample-size is always restricted by external constraints such as the experiment budget or the availability of suitable patients, the inference based on multiple peaks gives a shortcut to extracting more information from the limited set of samples, thereby directly addressing the “small n, large p” problem. however, some types of systematic measurement error, such as some batch effects, escape this treatment and can only be reduced by introducing independent replicates. based on the results presented in this work, we argue that additional peaks are especially useful when the signal-to-noise ratio is low and the differences between sample groups are small. we suggest that all the detected peaks that can be associated with a compound should be taken into account in the comparative analysis. this is possible through the two-step generative modeling approach presented in this work: by identifying the peaks that can be associated with one compound through clustering the peaks based on their shape similarity and by the inference of covariate effects on the clusters, each representing one compound. abbreviations anova: analysis of variance; cer: ceramide; sm: sphingomyelin; uplc-ms: ultra performance liquid chromatography-mass spectrometry. competing interests the authors declare that they have no competing interests. authors’ contributions the method was developed jointly by ts, sk and sr. ts had a lead role at implementing the model, designing and implementing the experiments, and at preparing the manuscript. all authors read and approved the manuscript. supplementary material additional file 1 supplementary material. more details of the experiments. click here for file acknowledgements the authors would like to thank sandra castillo, peddinti v. gopalacharyulu, mika hilvo and matej orešič for providing data and for useful discussions. the authors would also like to thank rónán daly and joe wandy for useful discussions. this work was supported by the academy of finland, the finnish foundation for technology promotion and the finnish doctoral programme in computational sciences fics. the calculations presented in the work were performed using computer resources within the aalto university school of science "science-it" project.
mass spectrometry-based metabolomic analysis depends upon the identification of spectral peaks by their mass and retention time. changes of concentration between sample groups can be inferred more accurately when multiple peaks are available. statistical analysis that follows the identification currently relies on one main peak of each compound. however, a compound present in the sample typically produces several spectral peaks due to its isotopic properties and the ionization process of the mass spectrometer device. in this work, we investigate the extent to which these additional peaks can be used to increase the statistical strength of differential analysis.
 mass spectrometry-based metabolomic analysis depends upon the identification of spectral peaks by their mass and retention time. statistical analysis that follows the identification currently relies on one main peak of each compound. however, a compound present in the sample typically produces several spectral peaks due to its isotopic properties and the ionization process of the mass spectrometer device. in this work, we investigate the extent to which these additional peaks can be used to increase the statistical strength of differential analysis.  we present a bayesian approach for integrating data of multiple detected peaks that come from one compound. we demonstrate the approach through a simulated experiment and validate it on ultra performance liquid chromatography-mass spectrometry experiments for metabolomics and lipidomics. peaks that are likely to be associated with one compound can be clustered by the similarity of their chromatographic shape. changes of concentration between sample groups can be inferred more accurately when multiple peaks are available.  when the sample-size is limited, the proposed multi-peak approach improves the accuracy at inferring covariate effects. an r implementation and data are available at http://research.ics.aalto.fi/mi/software/peakanova/. anova-type modelingbayesian modelingclusteringmass spectrometrymetabolomicslipidomicsnonparametric bayes 