Evaluation of Normalization Methods to Pave the Way Towards Large-Scale LC-MS-Based Metabolomics Profiling Experiments

Data

Background of the datasets

Two different datasets were included in this study. Set I was used to evaluate the different normalization methods, and Set II was used to validate the best performing normalization method (Fig. 1). Set I itself was composed of two subsets: (i) Subset A was obtained from LC-MS experiments on a freshly diluted mixture of 38 physiological amino acid standards [Product No. A9906, Sigma; preparation of this mixture is described elsewhere (Jankevics et al., 2012)], which were performed at three different time points, hereafter called time blocks: four LC-MS runs in time block 0 (T0), seven runs 2 months later (T2), and 11 runs 3 months later (T3), while (ii) Subset B was obtained from LC-MS experiments on a single, at −80°C stored, extract of Leishmania donovani (MHOM/NP/03/BPK282/0cl4), which was measured at two different time points: seven LC-MS runs at T0, and five runs at T2. It is important to note that biological variability is not present in Set I because it involves replicate measures of (i) a freshly diluted standard (Subset A) and (ii) one single parasite extract (Subset B). Set II consists of data acquired from LC-MS experiments on extracts of two clinical strains of Leishmania donovani with a different drug susceptibility towards antimonials: MHOM/NP/03/BPK282/0cl4 (sensitive) and MHOM/NP/03/BPK275/0cl18 (resistant) (Fig.1: Set II). Four extracts per strain per time block (T0 and T2) were prepared from independent cultures and measured soon after preparation. In this case, biological variability might be responsible for both intra- and inter-batch variability, as the data include biological replicates of the same strain (intra-batch) and freshly prepared extracts per time block (inter-batch). Methods for phenotyping and according metabolic differences related to antimonial-resistance for these strains have been reported elsewhere (t'Kindt et al., 2010b). Metabolite changes are regarded as biologically and significantly different if the average signal intensity differs at least a two-fold (e.g., fold change BPK275/BPK282 >2 or <0.5) and if the calculated p value of the t-test is smaller than 0.05. Set II will serve as a benchmark dataset with negative and positive controls used in the validation part of the manuscript.

OPEN IN VIEWER

FIG. 1.

Experimental design. Description of sample preparation, number of runs and source of possible batch variability for each experiment (T0, timeblock 0; T2, timeblock 2; T3, timeblock 3).

Employed normalization methods

Systematic changes in the average intensity levels across different experimental runs/timeblocks can obscure the biological information. As a consequence, the conclusions drawn from the non-normalized data may be incorrect. To remove the systematic variation across different experimental blocks observed in the original data, we have applied seven different data-driven normalization methods, which are commonly used for large-scale microarray datasets: (1) linear baseline normalization, (2) normalization by sum (analogous to total ion count normalization), (3) median normalization, (4) cyclic-Loess normalization, (5) quantile normalization, (6) probabilistic quotient normalization, and (7) cubic-spline normalization. These normalization methods were applied to both Subsets A and B of dataset I. More information about each of these normalization methods is available in the Supplementary Material. R scripts allowing applications of the described normalization methods to other data sources are also provided in the Supplementary Material (Supplementary Material is available online at www.liebertpub.com/omi).

Evaluation of normalization methods

In order to assess the performance of the different normalization methods, we assume that the majority of the metabolites are approximately similarly expressed across different experimental runs. As a result, the observed intensity measures for all metabolites are adjusted to satisfy this assumption. As described by Mar et al. (2009), Schmid et al. (2010), and Hill et al. (2001), this is a well-known assumption in gene-expression data normalization. Since the observed intensities are influenced by different factors, it is hard to predict the effect of this assumption on the assessment of normalization methods. By adopting the assumption is satisfied, the performance of the applied normalization methods was assessed using the variance and coefficient of variation between replicated measurements or time blocks. Both metrics are common statistical measures reflecting the spread (variation) of measurements across different experimental runs. Accordingly, the better the normalization method, the lower the coefficient of variation and the lower the variance between replicated measurements will be. More details on the calculation of these variability measures are given below.

Variance between replicated measurements

Let y_ij be the log intensity value for metabolite i (I=1,2,3..a) at run time j (j=1,2,3,…k), using the replicated observations due to different experimental runs. Our goal is to derive the variability measures of y_ij. It is expected that successful normalization should reduce the between-experimental run variability for each metabolite, as compared to the original data. A line plot for the variance was constructed to inspect visually which normalization method consistently reduces the variability across all the considered metabolites, as depicted in Figure 2 (top panel). The variance for the i^th metabolite \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sigma_i^2$$ \end{document} is computed as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\sigma_i^2 = \frac { \sum_{ j = 1 } ^k ( y_ { ij } - \bar { y } _i ) ^2 } { n - 1 } \tag { 1 } \end{align*} \end{document}

OPEN IN VIEWER

FIG. 2.

Line plots for the log-intensity variance for different metabolites of the amino acid standard across runs before and after normalization (a), and box-whisker plots for the coefficients of variation (b).

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sigma_i^2$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{y}_i$$ \end{document} are the variance and the mean of the i^th metabolite across different experimental runs, respectively; k represents the number of runs.

Coefficient of variation

The coefficient of variation (CV) is a measure that allows us to assess the degree of spread in a given data set relative to its mean value. It represents the ratio of the standard deviation to the mean. Since it is a unitless measure of variation, it is a useful statistic for comparing the degree of variation from one data series to another, even if the means are drastically different from each other. As a result, we used the CV (2) as an extra criterion to discriminate between different normalization methods. The better the normalization method, the lower the average metabolite specific coefficient of variation. The coefficient of variation for the i^th metabolite is defined as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}CV_i = \sigma_i / \bar{y}_i \tag{2}\end{align*} \end{document}

where CV_i, σ _i , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{y}_i$$ \end{document} are the coefficient of variation, standard deviation (square root of (1)), and mean for the i^th metabolite, respectively.

It is important to stress that, when the logarithmized data contains negative values, the CV is meaningless. Thus, for logarithmized data, the application of CV should be handled with caution. Since the log transformed data contains only positive values larger than 10 before (Fig. 3, lower channel) and after normalization (Fig. 4), we used equation (2) to compute CV.

OPEN IN VIEWER

FIG. 3.

(a) Principal component analyses (PCA) distinguishes replicate runs between timeblocks (T0, T2, and T3) for both amino acid standard (left, T0 with runs 1–4, T2 with runs 5–11, and T3 with runs 12–22) and Leishmania sample (right, T0 with runs 1–7 and T2 with runs 8–12). PCA is an unsupervised cluster method based on quantitative measurement of 28 amino acids and derivatives (left) and 189 metabolites (right). (b) Box-whisker plots for both amino acid standard (left) and Leishmania sample (right) before normalization.

OPEN IN VIEWER

FIG. 4.

Box-whisker plots for the log-intensity data (a) before and (b—f) after different normalization methods for the amino acid standard (Subset A). The top and bottom of the box represents the 75^th and 25^th percentiles of the distribution, respectively. The line inside the box represents the median value, and the dots indicate data points on the extreme of the distribution.

In addition to the numerical variability measures, box-plots and heat maps were prepared to visualize the effect of the employed normalization techniques. For the heat maps (Fig. 5), the intensities of each metabolite were rescaled between 0 and 100 to improve graphical interpretation. Moreover, an ANOVA model (3) was fitted to the data before and after normalization in order to test whether the batch effect was successfully removed from the data after normalization. ANOVA is a special case of regression model when all the explanatory variables are categorical/class variables (factors). The model was formulated as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}y_{ij} = \beta_0 + \beta_1 batch_2 + \beta_3 batch_3 + \varepsilon_{ij} \tag{3}\end{align*} \end{document}

OPEN IN VIEWER

FIG. 5.

Metabolic profiles of the 135 metabolites detected in dataset II in heat map format with hierarchical clustering before normalization (left) and after cyclic-Loess normalization (right). The samples are presented along the x-axis, the 135 detected metabolites are presented along the y-axis; the major classes of metabolites are color coded on the right. The intensity of each metabolite was rescaled between 0 (blue) and 100 (yellow). Before normalization, the clustering algorithm detects the measurement period as group. After normalization, the sensitive and resistant strains are properly separated by the clustering.

where y_ij is the log intensity value for metabolite i at run time j, and batch₂ and batch₃ are indicator variables for measurement taken from T2 and T3, respectively, for dataset I Subset A, and for Subset B only batch₂ kept in the model.

If the batch effects (time block) are successfully removed, their effects should become equal (β₁=β₂). This ANOVA model was also used to find differentially abundant metabolites for the validation data set.

Observed variation in the original dataset

Both principal component analysis and box-whisker plots (Fig. 3) distinguish a clear clustering per time block for the replicate runs of both the amino acid standard (Set I Subset A) and the Leishmania sample (Set I Subset B). While systematic variation between time blocks is present in all data (Sub-)sets, it is important to note that in each (Sub-)set, other plausible sources of variation inherent to the origin of the extracts or standards, can be present. As such, Subset A of dataset I likely also contains variation induced by random errors in sample dilution, while Subset B of Set I might contain variation induced by storage of the extract at −80°C between T0 and T2. Supplementary Figure 2S shows how the experiment-oriented biases influence the measured intensity for four randomly selected compounds from dataset I Subset A. Besides systematic LC-MS difference, Set II contains additional biological and technical components in its total variation induced by the independent cultures at the two different time blocks and the preparation of the extracts from these cultures at those time blocks, respectively.

For both subsets, the largest variability (presented by principal component 1 which shows the greatest variability in the dataset: PC1>80%) is observed between samples measured on T0 and T2 (Fig. 3a). The box-plots in Figure 3b further clarify that the signal intensity of the T0 dataset is much lower compared to the other time blocks (total ion count of T0 dataset is 7-fold lower). The second principal component presents the second largest variability in the dataset and shows the separation between T2 and T3 in the amino acid standard dataset (PC2=9.62%). An ‘intra-batch variability’ is observed for both datasets. For example, PC2=12.8% illustrates this variability for the Leishmania sample. It is unlikely that storage of the Leishmania sample at −80°C is responsible for the large ‘inter-batch variability’ because a similar variability (PC2=9.62%) is observed for the amino acid standard replicates that were freshly made for each time block, suggesting that the variability induced by storage is equivalent to the variability of pipetting out a fresh standard.

Variation after applying normalization techniques

To evaluate the performance of the different normalization methods, normalized data were presented in a box-plot (Fig. 4 and Supplementary Fig. 3S) and heat map format (Supplementary Figs. 4S and 5S). The variance and the coefficient of variation were calculated and presented in Table 1 and the line plot in Figure 2. For the purpose of direct visual comparison, only results from the normalization methods which have the same scale as the original data are presented. Figure 4 presents the box-whisker plots for the log-intensity of amino acid standard before (Fig. 4a) and after different normalization methods (Fig. 4b–4f). The same plots have been constructed for the Leishmania sample of Subset B of dataset I and can be found in the Supplementary Material (Supplementary Fig. 3S). Supplementary Figure 4S shows the heat map of the amino acid standard data before and after normalization. In the plotted heat map, the 22 LC-MS experiments are presented along the x-axis, and the 28 identified amino acids (and derivatives) are depicted along the y-axis. The original intensity of each metabolite in the sample set was rescaled to a range between 0 (blue) and 100 (yellow). A similar heat map was provided for the Leishmania sample (Subset B of dataset I, Supplementary Fig. 5S). In addition, Figure 2 (top panel) shows the variance of the replicated measurements for a given metabolite across the experimental runs, whereas Figure 2 (lower panel) shows the box-whisker plots for the coefficients of variation. Normalization by using a reference metabolite (i.e., an amino acid standard that remained constant across different runs) does not perform well as compared to the data-driven methods (Supplementary Fig. 6S). Thus, normalization on a limited set of internal standards is not easily extrapolated to the entire set of analytes in the data.

It becomes clear that all the employed normalization methods considerably reduce the variability observed in the original data (Fig. 3, Table 1, Supplementary Fig. 7S). First, Figure 4 clearly shows that normalization removes the block effects seen in the original data. After normalization, the mean intensities are similar across different experimental runs. For quantile normalization, the distribution of the normalized intensity is identical across experimental runs as expected for any type of datasets (Fig. 4f). According to the ANOVA model results (Supplementary Table 1S), and the box-plot of the CVs, showed a significant difference between normalized and non-normalized data [Fig. 2 (lower panel)]. The whiskers of the box-plots ARE lie on the same interval for the different normalization methods (Fig. 2b). Hence, there is no statistically significant difference between the different normalization methods (Fig. 2). These findings are also reflected in the heat maps of the data shown in Supplementary Figures 4S and 5S.

Nevertheless, the CV shows that the cyclic-Loess and cubic-spline normalization methods performed slightly better for the amino acid standards and the Leishmania sample, respectively, (Table 1) in terms of variance reduction. Hence, the cyclic-Loess normalization method was selected for the validation section, although, in general, all normalization methods performed well.

Validation of the selected normalization technique

In this Section, we evaluate the performance of the cyclic-Loess normalization on an additional dataset (Set II in Fig.1). In ideal circumstances, this dataset would be the result of a designed experiment, which contains positive and negative controls. The positive controls could be spiked-in metabolites with known differential concentrations, which could show that a given normalization removes unwanted variation while retaining the biologically interesting effects. In such a case, normalization should not attenuate nor amplify the measurements of the spiked-in concentration differences of the positive controls. The negative controls could be standards that are present at a constant concentration in all the samples of the experiment (cfr., Subset I of dataset I). In that case, normalization should not change the intensity measurements of the internal standards. For the purpose of the validation, we used the dataset described in the Material and Methods Section (Fig.1, Set II), which does not contain aforementioned controls. However, it should be pointed out that experiments similar to the one described above have been conducted in the past by t'Kindt et al. (2010b). In the experiment, metabolic differences correlated to antimonial-resistance have been reported and validated through biological assays for the Leishmania strains under scrutiny. The outcome of the previous study can serve as the set of positive and negative controls in order to evaluate to which extent the normalization influences the reported statistics or fold changes and whether the findings are reproducible. It should be mentioned that the dataset of t'Kindt et al. (2010b) was recorded in a different laboratory by using a different technology.

The distribution of the logarithmically transformed intensities of the 135 identified metabolites is visualized by the box-plots in Supplementary Figure 8S for the non-normalized and cyclic-Loess normalized data. The mean intensity level of the 135 metabolites measured on time block T2 are higher as compared to the measurements of T0 (Fig. 8bS). We can already state that the cyclic-Loess normalization was able to remove the systematic effect in a similar fashion as indicated by previous data sets (Fig. 8cS). Furthermore, it can be noted that, within a measurement period, the metabolite intensity distributions between the sensitive (BPK282) and resistant (BPK275) parasite isolates are comparable (Supplementary Fig. 8bS). The absence of such a pronounced systematic effect within a measurement period makes the data especially suited to study the impact of normalization on a downstream statistical analysis. In other words, by conducting the statistical analysis within the time block before and after normalization, we can evaluate if the number of differential findings and their consistency is influenced, since, in this case, there is no apparent need for data normalization.

A straightforward method to assess the quality of the normalization can be achieved by a hierarchical clustering approach (HCA) on the metabolite intensities (constructed with heat map.2 function from gplots package in R). The clustering algorithm looks for correlations between the metabolomics profiles and aims to group similar data together. The heat map in Figure 5 illustrates this principle. Clearly, the intensity differences between the measurement blocks are large, as indicated by the color code. These obvious differences are captured by the hierarchical clustering. As a result, metabolomics profiles within a time block are recognized as similar and consequently, are grouped together. The grouping is illustrated by the tree-like structure at the top of the figure. On the other hand, after normalization, the heat map seems more uniform and therefore the clustering algorithm groups the data according to other metabolite features. In this case, it can be seen that the resistant and sensitive strains are grouped together. Note that the grouping at the lower hierarchical levels in the treelike structure still reveals the time block structure, which can be due to technical differences in the preparation of the samples. This preliminary and simple analysis already indicates that normalization removes the unwanted systematic block effects while retaining the valuable metabolite information. In the next analysis, we evaluate the effect of normalization on the results of a statistical analysis.

A simple statistical analysis by means of volcano plots (Cui et al, 2003) (Fig. 6) is proposed. Metabolites are considered to be statistically and biologically significant when they have a p value below a 5% significance level based on a two-sample t-test and an average fold change larger than 2.

OPEN IN VIEWER

FIG. 6.

Volcano plot of the pooled data (a) before and (b) after normalization. The plots indicate that the systematic effect results in an increased variability, which yields lower −log10 (p values) (plot a).

In a first step, a Volcano plot analysis is conducted for each time period T0 and T2 separately using aforementioned thresholds. The number of differential findings, comparing the analysis before and after normalization, when comparing parasite isolate BPK275 with BPK282 in time period T0, is reported in Table 2a. Before normalization, 25 metabolites were upregulated, 48 metabolites were downregulated, and 62 metabolites were not differentially regulated. This information can be read from Table 2a in the row indicated with the label ‘Total Before’. The column with the label ‘Total After’ in Table 2a, illustrates the results after data normalization (i.e., 33, 32, and 70 metabolites found to be up-, down, and non-differentially-regulated, respectively). Because the aim of the study is to evaluate the effect of data normalization, the results are split into nine categories to illustrate the consistency of the differential findings before and after normalization. In ideal circumstances, when normalization would not influence the results of a statistical analysis at all, a complete agreement of the differential findings would be achieved. In such a situation, only the diagonal cells (highlighted in light gray) that indicate the number of metabolites which are found to be up-, not- or downregulated in the analysis before and after normalization, should contain values. However, in this study, 8 metabolites (highlighted in black in Table 2a) that were not differentially expressed using the data before normalization, were found to be upregulated after the normalization. Similarly, there are 16 metabolites (highlighted in black in Table 2a) that were downregulated before normalization, and were not differentially expressed after normalization. The same trend is also present in Table 2b (highlighted in black). A possible explanation is the presence of a slight systematic bias between the parasite isolates BPK275 and BPK282, as can be observed in Figure 8S. The figure illustrates that the log-intensity distribution of the BPK282 isolates are slightly higher than BPK275 when comparing within the time periods. Nevertheless, it is reassuring that the majority of the findings for the 135 metabolites are in the diagonal cells, highlighted by light gray in Table 2 (i.e., are consistently found before and after normalization). This result indicates indeed that within a time period, a good reproducibility of the measurements was achieved with negligible systematic effects. The effect of an unnecessary normalization does not alter the outcome of a statistical analysis severely. A similar observation can be made when looking at the findings in time period T2 displayed in Table 2b.

OPEN IN VIEWER

Table 2.

Result of Statistical Analysis Between BPK275 and BPK282 on 135 Metabolites with p<0.05 and Fold Change of >2, Before and After Normalization

The consistency of the findings before (columns) and after (rows) cyclic-Loess normalization are split into over nine categories. The values in a particular column illustrate how the findings before normalization are distributed over the rows, which represent the cases, up-, not-, or downregulated after normalization. Panel a) and Panel b) contain the result for time period T0 and T2, respectively.

Second, in a statistical analysis it is good practice to include all available data because a larger sample size enables a more powerful test to discern differential metabolites. Instead of performing an analysis for each time period separately, it is preferable to pool the collected data. Therefore, the same statistical test, as described earlier, was applied on the pooled data set, before and after normalization. Panel a of Table 3 contains the results of the analysis before data normalization and compares them with our benchmark dataset published by t'Kindt et al. (2010) and described in Materials and Methods. The table should be interpreted similarly as the consistency matrix presented in Table 2. It can be noted that pooling the data without accounting for systematic variation is undesirable, as only 1 out of the 31 upregulated metabolites from t'Kindt et al., is found to be upregulated in dataset II. Furthermore, 11 out of 19 downregulated metabolites from t'Kindt et al. are downregulated in dataset II. A possible explanation for this poor agreement is that, in this case, pooling the data sets increases the variability present in the data and makes it more difficult to detect statistically significant differences regardless the higher number of data points in the pooled analysis. In fact, there seems also to be a bias towards finding mostly downregulated metabolites. For these reasons, it is mandatory to remove the systematic bias from the data prior to a statistical analysis. Panel b of Table 3 contains the result of the statistical analysis after applying normalization. It can be observed that the agreement between the differential findings of t'Kindt et al. (2010b) and dataset II is improved (i.e., 16 out of 31 and 15 out of 19 are consistently found up- and downregulated, respectively). However, the agreement between both studies is still poor, regardless which normalization strategy is used. The cells highlighted in black from Panel b of Table 3 illustrate the disagreement. For example, in dataset II, there are 13 additional upregulated metabolites that could not be found by t'Kindt et al. However, it is reassuring that 9 of the 13 additional metabolites were below the limit of quantification in the study of t'Kindt et al. and consequentially were classified as not regulated. A similar conclusion can be drawn for the 12 downregulated metabolites. Closer inspection of the results confirm that, although the metabolites are not found differential according to the volcano plots, the same trends with respect to up- and downregulation were observed (Fig. 6). As depicted in Supplementary Figure 9S (Supplementary material are available online at www.liebertpub.com/omi), two-fold change increases or decreases are indicated by blue and red dots for each identified analyte. A possible reason for these inconsistent findings is due to inter-laboratory variation. It is worth keeping in mind that the results of t'Kindt et al. were obtained one year earlier (2010) on another LC-MS platform (2.1 mm HILIC column coupled to a Finnigan LTQ-Orbitrap XL mass spectrometer) and that the biological samples, although originating from the same Leishmania clone and prepared by the same protocol in the same laboratory settings, might be slightly different. A possible conclusion from these observations could be a more optimal validation data set is required.

OPEN IN VIEWER

Table 3.

Result of Statistical Analysis Between BPK275 and BPK282 on 135 Metabolites with p<0.05 and Fold Change of >2 for Dataset II and Findings of t'Kindt et al.

The table evaluates the consistency of the finding before (Panel a) and after (Panel b) cyclic-Loess normalization, results are split over nine categories. The rows contain the results of a statistical analysis on dataset II. The columns represent the results of t‘Kindt et al. (2010b).

However, when comparing the differential findings on dataset II before and after normalization across Panels a and b of Table 3, as highlighted by bold and italic numbers in the columns indicated by ‘Total’, one can conclude that pooling without taking into consideration the systematic effects is unfavorable for a downstream analysis: after normalization, 29 and 27 differential metabolites could be found opposed to the 17 differential findings before normalization. The numbers indicated by bold and italic in Panel b of Table 3 could be contrasted with the bold and italic numbers of the two panels in Table 2. Indeed the number of differential metabolites found when analyzing the separate time blocks before normalization corresponds to the findings presented in Panel b of Table 3.

Third, to put previous results in a better perspective and to illustrate that the inconsistent findings are not an artifact of the data normalization, we compared the statistical analysis on the time blocks T0 and T2 before normalization separately with the findings of t'Kindt et al. (2010b). Such a comparison is valid because only a slight systematic bias was present in the time blocks, as indicated by Table 2. The results of the comparison are displayed by the consistency matrix of Supplementary Table 2S (Supplementary Data are available online at www.liebertpub.com/omi) for time period T0 (Panel a) and T2 (Panel b). The cells highlighted in black indicate inconsistent findings across the comparison and confirms our hypothesis that these changes are primarily caused by inter-laboratory variation. More information on the consistency of findings is provided in Supplementary Figure 9S (Supplementary Material is available online at www.liebertpub.com/omi), which illustrates that in this example downregulated metabolites are consistently found disregarding the normalization.

Another way of illustrating the consequences of not properly accounting for the systematic effects is provided by the volcano plot presented in Figure 6. The plots indicate that the systematic variation results in an increased variability, which yield lower −log10 (p values). In addition, there seems to be a linear trend in the averaged ratios (resistant/sensitive strain) when compared to the normalized data, where the ratios are more scattered in the plot.