Influence of Missing Values Substitutes on Multivariate Analysis of Metabolomics Data

Missing values are known to be problematic for the analysis of gas chromatography-mass spectrometry (GC-MS) metabolomics data. Typically these values cover about 10%–20% of all data and can originate from various backgrounds, including analytical, computational, as well as biological. Currently, the most well known substitute for missing values is a mean imputation. In fact, some researchers consider this aspect of data analysis in their metabolomics pipeline as so routine that they do not even mention using this replacement approach. However, this may have a significant influence on the data analysis output(s) and might be highly sensitive to the distribution of samples between different classes. Therefore, in this study we have analysed different substitutes of missing values namely: zero, mean, median, k-nearest neighbours (kNN) and random forest (RF) imputation, in terms of their influence on unsupervised and supervised learning and, thus, their impact on the final output(s) in terms of biological interpretation. These comparisons have been demonstrated both visually and computationally (classification rate) to support our findings. The results show that the selection of the replacement methods to impute missing values may have a considerable effect on the classification accuracy, if performed incorrectly this may negatively influence the biomarkers selected for an early disease diagnosis or identification of cancer related metabolites. In the case of GC-MS metabolomics data studied here our findings recommend that RF should be favored as an imputation of missing value over the other tested methods. This approach displayed excellent results in terms of classification rate for both supervised methods namely: principal components-linear discriminant analysis (PC-LDA) (98.02%) and partial least squares-discriminant analysis (PLS-DA) (97.96%) outperforming other imputation methods.


Introduction
Metabolomics is the study of all metabolites in a biological system under a given set of conditions [1]. The classical technologies for the analysis of metabolites (i.e., chemical entities) include: Gas chromatography-mass spectrometry (GC-MS), liquid chromatography-mass spectrometry (LC-MS), capillary electrophoresis-mass spectrometry (CE-MS), nuclear magnetic resonance (NMR) spectroscopy, Fourier transform-infrared (FT-IR) spectroscopy, and many more [2]. These methods have been widely applied in countless metabolomics research studies. These include for example: cell line studies probing intracellular metabolites [3]; animal models for cancer study, such as mouse tumours [4]; as well as general metabolic profiling for the analysis of human serum [5]. Unfortunately not all of these methods generate a complete data set, for instance GC-MS and LC-MS analyses employ chromatographic separation prior to MS and thus require a complex deconvolution step to transform these 3D data matrices into lists of annotated features (metabolites) with (relative) abundances. The result of this process is often a matrix that is not fully populated and thus presents a major problem in data processing due to these missing values [6].
We believe that one of the most important steps in GC-MS data pre-processing of metabolomics data, next to peak detection and alignment, is handling missing values [7]. In our experience, typically 10%-20% of GC-MS data from an experiment are missing, with the same level of missing values also reported in many other studies (e.g., Hrydziuszko and Viant (2011) [8] and Xia et al. (2009) [9]). Whilst it is obvious that a missing value is when a matrix contains an empty cells usually recorded as -Not a Number‖ (NaN) or more worryingly as a zero which cannot be readily distinguished from the real absence of a feature rather than a failure in the analysis. The obvious question that arises from these observations is: What are the roots of these missing values? Prior to any analysis it is good practice to identify the origins of missing values whether they are truly missing or not [10]. Missing values may arise due to various reasons, such as: (1) limits in computational detection; (2) imperfection of the algorithms whereby they fail in the identification of some of the signals from the background; (3) low intensity of the signals used; (4) measurement error; and finally (5) deconvolution that may result in false negative during separation of overlapping signals [6,8,[10][11][12][13][14][15][16][17].
Currently, the most well-known substitute for missing values is their substitution with a mean value [6]. In fact, some researchers do not specifically state how this aspect of data analysis in their metabolomics pipeline has been performed and use this replacement approach as a common practice. This is despite the fact that this problem has been well recognised in the literature as an important aspect and also appears within the minimum reporting standards for data analysis for metabolomics [7]. Methods which have been reported in the literature include: replacing missing values by half of the minimum value found in the data set [9]; missing value imputation using probabilistic principal component analysis (PPCA) [15], Bayesian PCA (BPCA) [18] or singular value decomposition imputation (SVDImpute) [9]; replacing missing value by means of nearest neighbours [6]; or replacing the missing values with zeros [19].
Whilst many unsupervised and supervised learning for the analysis of high dimensional metabolomics data require a complete dataset [2,6,7,20]. Consequently, there is a need to analyse and identify correct approaches for the replacement of missing values. Nevertheless, this seemingly important aspect of data pre-processing has not received wide attention in this field. Herein, we investigate this problem using a common set of metabolomics data produced using GC-MS which contained ~15% missing values, where the goal of the study was to analyse cancer cell lines in terms of changes in oxygen level and which metabolite features change during this process [21]. Therefore, we evaluated five relatively simple potential missing values substitutes-zero, mean, median, k-nearest neighbours (kNN) [22] and random forest (RF) [23] replacements-in terms of their influence on unsupervised and supervised learning and thus their impact on the final output(s); these outputs are related to cluster compactness from replicate biological measurements. Moreover, to our knowledge these approaches have not been compared directly for the analysis of GC-MS data.

Cell Culture and Experimental Protocol
The methods used for cell culture have previously been described [21]. Experimental analysis proceeded as follows: MDA-MB-231 cells were seeded and allowed to adhere for 24 h in 95% air and 5% CO 2 . Cells were divided into three groups: one group was placed in a 95% air and 5% CO 2 incubator (normoxia); one group placed in a 1% O 2 , 5% CO 2 balanced with N 2 hypoxybox (hypoxia) and one group was placed in an anoxic chamber (Bactron anaerobic chamber, Sheldon Manufacturing, Cornelius, OR, USA) where 5% CO 2 , 5% H 2 and 90% N 2 (BOC, Manchester, UK) was flowed over a palladium catalyst to remove any remaining oxygen (anoxia) for 24 h. Each of the three groups were split into three sub groups, which were dosed with 0, 0.1 or 1 µM doxorubicin for 16 h whilst remaining in the predefined oxygen condition for a further 24 h.

Methanol Metabolite Extraction
Metabolites were extracted whilst cells remained in the predefined oxygen condition. The hypoxybox was sealed and cycled into the anoxic chamber prior to the extraction of metabolites. Extracellular medium was decanted and cells were washed with three 1 mL aliquots of PBS. A total of 1 mL of methanol (−48 °C) was added to each dish to quench the metabolism [24]. Cell scrappers detached adherent cells from the Petri dish allowing the intracellular metabolites to be extracted into the methanol. The solution was transferred into a pre-weighed Eppendorf tube followed by three freeze thaw cycles (frozen using liquid nitrogen). Cells were centrifuged at 17,000 × g for 15 min to separate the pellet from the supernatant. The supernatant was transferred into a fresh Eppendorf tube and the pellet was lyophilised and weighed. The volume of supernatant lyophilised was normalised according to the dry weight of the pellet for intracellular samples. A QC was produced using 150 µL of remaining supernatant from each sample. Aliquots (1 mL) of the QC solution were placed into Eppendorf tubes. Both the sample and QC were lyophilised.

Metabolic Profiling Using GC-MS and Raw Data Processing
A total of 8 randomised batches were analysed using GC-MS. Of these batches, 7 contained 40 samples and 9 QCs each and 1 batch contained 30 samples and 9 QCs. Samples were analysed on an Agilent 6890 GC (Agilent Technologies, Stockport, UK) coupled to a LECO Pegasus III (Leco Corp., St. Joseph, MO, USA) EI-ToF-MS. The GC-MS instrument setup used has been previously described [5,25]. QCs were used to condition the GC-MS before sample injection. A total of 5 injections of QC samples were performed prior to analysing samples to normalise the chromatographic conditions to the sample matrix. Every 5 sample injections were followed with a single QC injection. In addition, three QC injections were made at the end of each analysis. These QC samples were compared to account for the system stability in each batch analysed. Briefly, the temperature was set at 70 °C for 4 min followed by a 20 °C increase every min until 300 °C was reached and stabilised for 4 min. Samples (2 μL) were injected onto the column. The mass spectrometer was operated from m/z 45-600. The total duration per sample was 25 min.
The processing of raw GC-MS data were performed following the methods described previously, using the LECO ChromaToF v3.25 software package to apply the chromatographic deconvolution algorithm [5,25]. Chromatographic peaks from the QC samples selected from each analytical block were manually integrated and added to a database, which included the peak area, the retention index and mass spectrum to compile a reference database. This reference database was manually inspected to ensure the mass spectra had the required fragmentation patterns and met the criteria for signal-to-noise (S/N) ratio and chromatographic peak width. Metabolites with fewer than 50% of features detected in the QC were removed [26].
Identification of metabolites within the compiled database was performed through searching against an in-house mass spectral and retention index library as described previously, where a high confidence match was assigned when the mass spectral match was greater than 80% and the retention index was within 30 [27]. All other matches that did not meet this specification were not regarded as a high confidence match. The metabolite library was previously developed in-house from authentic standards at the University of Manchester. Metabolites that were not identified in the in-house database were searched against the Golm metabolome database [28] and a similar match score to the in-house database was applied.
In this dataset, 102 peaks were integrated from the GC chromatogram of which 45% were unidentified peaks; these unidentified peaks were removed from the data matrix. The coefficient of variation (CV) was calculated for metabolites detected in QC samples. Two metabolites with CV > 30% were removed from the whole data matrix as these were considered to have poor reproducibility [26]. We followed the metabolite identification protocols outlined by the metabolomics standards initiative [29] where the in-house matches are MSI level 1 (43 metabolites) and all other identification in this manuscript to MSI level 2 (9 metabolites).
Finally the GC-MS data were normalised to the internal standard succinic d 4 acid to account for the technical variability associated with chemical derivatisation and low sample injection volumes. The chromatogram peaks areas were exported as an ASCII file into Microsoft Excel ® and these 52 metabolites were used for all further analyses.

Software Tools
All data analyses used in this work were conducted in R version 2.15.0 and 3.1.0; R is a software environment [30] that contains a variety of packages useful for data analysis and is a free open source program [30]. The following R packages were employed for visualization and statistical analysis; -chemometrics‖ [31] package has been used to determine the number of principal components (PCs) with cross-validation and for hierarchical cluster analysis (HCA) and picturing; -MASS‖ package was used for linear discriminant analysis (LDA) [32]; -rgl‖ package was applied as a 3D visualization device [33]; -mixOmics‖ [34] package was used to perform both PCA and partial least squares-discriminant analysis (PLS-DA); package -impute‖ [35] which is a part of Bioconductor an open source and open development software project to provide tools for the analysis and comprehension of genomic data [36] was applied for substitution of missing values using kNN; and finally -missForest‖ [23,37] for missing value imputation using RF.

Data Preparation
Data pre-treatment is a very important process because of its ability to make the data clearer and suited to analysis [38]. Here, we used autoscaling as it is the probably most reliable method [7,39]. This was conducted after replacing missing values. In autoscaling for each column (input feature) the mean value of that column is subtracted, followed by dividing the row entries for that column by the standard deviation within the same column. In other words, the aim of autoscaling is to break down relation with irregularities in order to produce slighter, well-structured relations, which in practice allow for better data recognition [39,40].

Imputation Methods
In this study five different imputation methods of missing value were applied: zero, mean, median, kNN and RF. Zero imputation was used to replace all missing value with zeros. Mean imputation to replace the missing values with average value of corresponding metabolites (variables). The median imputation corresponds to central value of metabolite that split its values into two equal parts the higher and the lower part.

Imputation of Missing Values Using k-Nearest Neighbours (kNN).
kNN for each metabolite identifies missing values by recognizing the nearest neighbours of that feature using Euclidean distance [22]. The calculation takes into consideration that each candidate neighbours might miss some of the coordinates. As a result, the method uses the average distance from the non-missing locations. As long as the k nearest neighbours for each metabolite is identified, the imputation of missing value is computed by averaging of non-missing values of its neighbours. Caution need to be taken when one want to deal with the cases where all of the neighbours are missing in a particular group. Nevertheless, this obstacle is circumvented by taking the overall column mean for that particular metabolite. Consequently, the technique has the ability to deal with the data set that contains large number of variables that include missing values [22].

Missing Value Imputation Using Random Forest (RF)
RF is a classification and regression technique that can handle both parametric and non-parametric data sets of complex linear and non-linear problems [41]. The approach is based on the estimation of imputation error which is calculated for the bootstrapped samples and therefore no separate cross-validation is required. The process begins with splitting the missing data into training and test sets. The training set is used to generate a group of trees for the observed values in the training data set to predict the missing values. Consequently, the first step is to replace each of the missing data points from training set with a mean of the observed values for that variable. All missing values are replaced and the proximity matrix is calculated for first iteration set. The average weights calculated from the proximity matrix for first iteration are used to replace the missing values in the second iteration and so on until the stopping criterion is satisfied. The stopping criterion is passed based on the difference between the first and second iteration, then between second and third and so on. This is achieved when both differences have become larger once. As soon as this is accomplished the final matrix with replaced values is used to perform further analysis [23].

Unsupervised Learning
Unsupervised learning algorithms attempt to find hidden structure in data, which displays any similarities between data points (e.g., metabolite levels), such as how close they are to each other and whether the data contain any outliers [42]. Additionally, these so-called exploratory methods are mainly based on cluster analysis. Cluster analysis aims to cluster data (i.e., samples) into constituent groups which describe common characteristics that they possess. These clusters should display high internal homogeneity within clusters (or as high as possible), and high external heterogeneity between clusters.

Principal Components Analysis (PCA)
The most frequently used exploratory analysis method is principal component analysis (PCA) which is used for computing linear latent variables (components). PCA is possibly one of the oldest and most well-known methods of multivariate data analysis (MVA) method which was invented in 1901 by Karl Pearson [43]. This technique is one of the simplest multivariate statistical procedures that reduce the dimensionality of data sets comprised of a large number of possibly correlated variables. On finding these correlated observations, Hotelling showed that PCA retains these features allowing any variability present in the data to be explained and visualized [44]. In general, the correlated original variables are transformed into a set of new uncorrelated variables, (principal components (PCs)), ordered according to their decreasing variance, whereby the first PCs describes more variability than any other PC. The absence of correlation is a valuable property because it means that these indices (PCs) are assessing dissimilar correlations in the data [45].
Rather than generate a single PCA model we used k-fold cross-validation in order [46,47]: (i) to generate a generalized scree plot which allows us to assess random sampling of the population; and (ii) to calculate how many times the 3 groups were recovered (for this we calculate Q 2 based on the test sets from the k-folds). Thus for PCA we applied 10-fold cross-validation and we repeated this 100 times. This allowed the explained variance for each model to be assessed. These was performed to determine the optimal number of PCs and from 1 to n PCs were used in this process and the overall percent explain variance represented as a box-whisker plot; further explanation of the PCA procedure can be found in supplementary information (SI).

Hierarchical Cluster Analysis (HCA)
There are numerous clustering approaches such as partitioning clustering which is stepwise process that work by partitioning observations, optimization-partitioning which allows for reassignment of objects to optimize an overall parameter, density or mode-seeking methods where cluster are created based on dense concentration of entities, and clumping with possible overlapping samples. However, in this study we have selected the most commonly used method, namely hierarchical clustering analysis (HCA) [48][49][50]. This approach is a fast and reliable technique, allowing for a huge amount of closely related samples or features to be organized into separated clusters [50].
Joining of clusters is based on similarities and dissimilarities in metabolites that allow for the separation of groups, with the most important factor in the creation of cluster beings the distance between two objects. Distances in clustering can be measured in various ways. One of the most common and frequently used is Euclidean distance. However, this can be only used to measure the distance between samples, but not a group of samples constructing clusters. Therefore, here in this study -Ward‖ linkage was used to link groups of different clusters together and provides a clear visualization of the results [49][50][51][52]; further explanation of the HCA can be found in SI.
In general, unsupervised learning is applied when there are no unlabeled training data (i.e., the modeling is performed on only the X data (SI)). On the other hand, there are supervised techniques that are mainly focused on classification of an input signal, which is to be defined as a member of predefined classes (so called Y data (SI)).

Supervised Learning
Supervised learning might be represented by a variety of methods and a particularly popular method used is discriminant analysis. By way of an example we can consider the scenario where we want to use supervised methods to classify whether cell samples under investigation by GC-MS metabolomics are cancerous or not. In this process a model is first constructed (calibrated) using a data training set of X data (metabolite levels) and Y data (cancer versus control) whose identities are unequivocally known. Obviously, it is essential to establish the training set prior to any analysis and this must be fully representative of the problem domain [53].

Linear Discriminant Analysis (LDA)
Linear discriminant analysis (LDA) (also known as discriminant function analysis (DFA)) [54] is a statistical tool for studying the association between a set of predictors versus a categorical response. This method aims to find the best linear separation between groups by determining the minimum dimensions at which groups can be separated. LDA is used in statistics as a dimensionality reduction technique in many classification problems, as well as being a useful variable selection method. In general, LDA is a procedure to find a linear combination of features, which describe or separate two or more groups of samples. In other words, LDA classifies data by providing more distinguishable classes using multivariate decision boundaries [54]. LDA can be used without prior PCA if the following condition is satisfied [38,[55][56][57]: (1) where N s correspond to number of samples, N g to the number of groups, and N v reflects the number of inputs (features).
However, LDA cannot be used with data which are highly collinear, such as is the case in the present study. Therefore, in this investigation PCA was performed prior to LDA to remove the effect of collinearity and reduce the number of variables in the X data matrix.

Partial Least Squares-Discriminant Analysis (PLS-DA)
Another widely used statistical tool for classification is partial least squares-discriminant analysis (PLS-DA) [58] which can be used either for both classification and variable selection [59]. The approach is an extension of the traditional PLS regression model which was developed for quantitative predictions. The main difference is that in PLS-DA, a -dummy‖ response matrix was used to code multiple classes and the algorithm is calibrated to predict group membership. The maximum separation between classes (groups) is established by simultaneous rotation of the X and Y components [60,61].

Model Validation
To assess the classification rates for PC-LDA and PLS-DA for each of the five missing value imputations we performed bootstrapping (with replacement). This approach proposed by Efron was used to minimize the risk of misrepresentative results as a single training set-test set split may not be sufficient enough to draw valid conclusions. Bootstrapping was performed 100 times and on average training sets will contain ~63.2% and test sets ~36.8% of all samples [62,63]. Therefore, the number of components for both PLS-DA and PC-LDA has been estimated based on percentage of correct classification in validation.

Results and Discussion
As discussed above in order to compare five different missing value substitutes we used data that had been generated from the metabolic profiling of MDA-MB-231 breast cancer cells cultured in three oxygen levels-normoxia, hypoxia and anoxia-that had been treated with 0.1 or 1 μM doxorubicin drug, along with a control where no drug was used. Previously, for this data set the following analyses were performed: two-way ANOVA, PCA, PC-DFA, and correlation analysis [21,64].
In this study for the first approach we used the full data set which included all classes. These were 3 aeration environments times 3 drug doses: (1) normoxia; (2) normoxia + 0.1 μM doxorubicin; (3) normoxia + 1 μM doxorubicin; (4) hypoxia; (5) hypoxia + 0.1 μM doxorubicin; (6) hypoxia + 1 μM doxorubicin; (7) anoxia; (8) anoxia + 0.1 μM doxorubicin; and finally (9) anoxia + 1 μM doxorubicin. Initially we looked at PCA and the discriminant analyses on these 9 groups. However, due to the relatively large number of groups this separation is very difficult to visualize as is illustrated in SI Figure S3. To try and uncouple the environment and the drug dosing we applied multiblock PCA (see SI for details) that should provide us with a graphical overview of separation between these two different projections [65]. However, this approach did not work very well (see SI Figure S4). Therefore, our analysis was performed for three different environments, as the effects of normoxia, hypoxia or anoxia were the dominating phenotype.
Visualization of the scores generated from unsupervised and supervised methods is the main tool used in this study for evaluating the effect of the missing value substitutes (zero, mean, median, kNN and RF). The general idea being that for the three different groups these should group close to each other and away from the other groups. That is to say, (i) more compact clusters are more desirable; and (ii) the greater the distance between clusters is also a preferred outcome.
For the first analysis we performed PCA and evaluated the total explained variance (TEV) captured in the principal components (PCs). Rather than doing this on all of the data in a single experiment we used 10-fold cross-validation, where 9 folds were used to generate each PCA and the 10th fold was left out, this was repeated until each fold had been left out once; this process was repeated 100 times with random splits into the 10 different partitions/folds. Figure 1 is a scree plot which illustrates how each of the different substitution can influence variance across the analysis. Here in this study 80% cut off point has been selected as criterion in estimation of the number of PCs to be used. However, this might vary from study to study and depends on the amount of information that needs to be included in the model. As can be observed in Figure 1 when RF was applied as a substitution of missing value only 1-14 PCs were required to cover over 80% of explained variance. This result is followed by kNN method where 1-18 PCs where needed to cover similar level of explained variance. For zero replacement 1-21 PCs were used to reach comparable level of variance this is followed by median with 1-27 PCs were required to cover 80% of explained variance, and the mean replacement approach required the highest number and the first 37 PCs (out of a possible 52) were needed to cover same level of variance. As ~15% of the data were zeros it is not surprising that such good results for zero were found as the variance will obviously be lower. What is interesting is that the median captures more variance than the mean imputation approach. Indeed this is very evident in just PC1 where for zero and median the first PCs characterize ~30% TEV whilst the mean explains <20%. Such a good results for RF might originate from the fact that the methods do not take into consideration the samples distributions.
Having established the shape of the PCA plots for the five different imputation methods, the scores plots for PC1 vs. PC2 were generated when all the metabolomics data was used (Figure 2). 3D plots of PC1 vs. PC2 vs. PC3 are available in SI ( Figure S5). Visual inspection of these PC scores plots indicates that despite some overall in all five plots, the separations is not quite as good for zero and mean imputation, whereas the median, kNN and RF replacement results in reasonably good separation between three groups.   Figure 3 shows dendrograms based on Euclidean distance and -Wards‖ linkage. Note that these dendrogram can be used to characterise particular clustering and the distance gives an element of multilevel hierarchy [49]. To identify grouping within the dendrogram we have coloured the text according to the atmospheric environment in which the cells were cultivated and used boxes to show where the majority of the same samples fall. From the HCA we can see that different substitutes of missing values have a significant impact on clustering and it is clear than when the mean is used to replace the missing values that the clustering does not reflect the three expected groups. By contrast, all other replacement methods generate dendrograms that more-or-less recover the three environmental condition groups, with obvious exceptions that can be seen in Figure 3. This can also be observed in Table S1 (SI) that represent a purity measure of HCA; the percentage of samples coming from the same class as one cluster.
Although PCA and HCA are useful because they are unsupervised analyses methods the interpretation is somewhat subjective as cluster plots and dendrograms are used. Moreover, the interpretation is after the cluster code is broken (i.e., the labels or colours are put onto the clusters) and this is not really appropriate for providing information about classification rate. Therefore, in the next step we investigated the impact of the different missing value replacements on supervised methods based on discriminant analyses.
As described above bootstrapping was used to construct 100 different PC-LDA or PLS-DA models and inspection of these allows one to predict how accurate classification is. For PC-LDA this is in terms of recovery of test data with appropriate training data in LDA scores space, whilst for PLS-DA the Y predicted matrix is inspected and the binary encoding in this PLS2 model used. That is to say the encoding for normoxia is 100, hypoxia is 010 and for anoxia the output would be 001. Figure 4 highlights how the different replacement approaches can have an impact on the average correct classification rates of the 100 models for PC-LDA ( Figure 4A) and for PLS-DA ( Figure 4B) and also shows the effect of the number of components used on this classification accuracy. To aid interpretation this plot is annotated with the points at which the highest classification occurs. For PC-LDA the preferred model is the one with RF imputation demonstrated an excellent classification rate of 98.02% with 22 latent variables. In addition, this method performed the best along all components in comparison to other replacement techniques. Second best result of 95.63% classification rate was recorded for kNN. Zero imputation as whilst this is slightly lower (92.02%) than the median imputation (92.57%), however fewer latent variables are used, 7 compared with 22 PCs. Whilst for PLS-DA both zero and median imputation give very similar predictive abilities (90.97% vs. 91.76%) and use a very similar number of latent variables (7 vs. 8). In terms of best performance for PLS-DA yet again RF (97.96%) and kNN (96.06%) accomplished the best in comparison to other methods; however the difference between both is smaller as for PC-LDA. What is clear from both PC-LDA and PLS-DA is that mean imputation ( Figure 4) is significantly worse and predictions are only ~78%. However, these differences in classification rates are very small as it could be expected as these both methods should display similar results as reported in the literature [7,58,66,].
An additional output from PC-LDA is the scores plot which shows the best group separation, and Figure 5 shows LD1 vs. LD2 scores. It is clear that the 3 different groups of samples are generally well separated; however it is also clear that median and RF imputation have resulted in the best separation, followed by zero and kNN replacement and lastly mean substitutions for missing values. A very similar pattern can be observed for the PLS-DA where median seems to be the best substitute for missing value (data not shown).

Conclusions
A summary of the class prediction results are shown in Table 1, and it is clear from this that RF performed the best. The RF imputation technique is followed by kNN, median and zero replacement methods that perform relatively well. By contrast, the mean imputation approach does not perform well which is an interesting observation, as this is the default method used by many researchers. An obvious question arises as to why this has been the case? The answer to the question is likely to be based on the concentration distribution profile of the metabolites collected. If a particular metabolite's concentration follows a normal distribution then the mean imputation would be appropriate as the normal distribution has a Gaussian shape. Clearly in our case this is not the case as can be observed in Table 2. Almost half of the samples do not follow normal or near normal distribution as their values are highly positively or negatively skewed. Where for Normoxia samples the highest value (3.01) of positive skewness is recorded for galactose/glucose, followed by glycerol for Hypoxia with the highest value of 4.03 positive skewness, and finally benzoic acid with 4.15 positive skewness for Anoxia. Additional statistics such as percentage of missing values, mean, median, standard deviation, kurtosis and standard error can be found in Tables S2-S4.  Different conclusions were given when mean and median imputations were compared recently for data generated from direct infusion Fourier transform ion cyclotron resonance mass spectrometry (DI-FT-ICR-MS) [8], where mean and median showed nearly same results. This suggests that these DI-FT-ICR-MS data may had been normally distributed with little or no skews or outliers in the metabolite concentrations [6]. However, in our study we have shown that metabolites that display skewness can have a significant impact on the substitute of missing values and therefore, final output of the analysis. Where, skewness is an important feature that can articulate if the data are normally or near normally distributed where normal distribution has skewness of 0 or close to 0. As shown in this study not all metabolomics data follows normal distribution. Hence, we believe that careful analysis of distribution of the data prior to any analysis should be performed. Moreover, the authors recommended kNN, mean and median as imputation methods for missing values in metabolomics field. However, in our study we have shown that much better results can be accomplished when RF is used as a replacement method of missing values.
We have demonstrated that the selection method chosen to replace missing values may have a dramatic impact on GC-MS based metabolomics data, both for unsupervised and supervised learning. Therefore, the handling of missing values is an absolutely vital step in data pre-processing, as has been suggested previously in the literature [7], to which special consideration should be given. Moreover, we have shown that selecting the median metabolite level as a substitute of missing values is a substantially better option than the selection of the mean metabolite response, which is sensitive to metabolite distribution when data are skewed, as shown in this study. As observed in [8] when metabolite levels are normally distributed then both mean and median imputations will provide very similar results, but one cannot be guaranteed that this is always be the case. The median is less vulnerable to outliers than the arithmetic mean and is therefore considered to be more robust [6]. We conclude that before doing a missing value imputation that the distribution of the data be examined; conducting a skewness calculation or more simply observing the ratio of the mean and median of the distribution could identify this. We therefore believe that an awareness of these factors needs to be concomitant with an understanding of the data themselves, prior to any data analysis when missing values are observed [13], else the conclusions of the study could be inaccurate. In summary, in the case of GC-MS metabolomics data studied here our findings recommend that RF should be favored as an imputation of missing value over the other tested methods, as this method provides more robust and trustful results.