Marie Perrot-Dockes

Celine Levy Leduc

Julien Chiquet

Laure Sansonnet

Margaux Brégère

Marie-Pierre Etienne

Stephane Robin

Grégory Genta-Jouve

M Perrot-Dock

C L Ès

J Évy-Leduc

L Chiquet

M Sansonnet

É G Br

M.-P Ère

S Étienne

G Robin

Genta-Jouve

A variable selection approach in the multivariate linear model: an application to LC-MS metabolomics data

Keywords: Variable selection, high-dimension, multivariate linear model, metabolomics

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Metabolomics aims to provide a global snapshot (quantitative or qualitative) of the metabolism at a given time and by extension, a phenotypic information, see [START_REF] Nicholson | metabonomics': understanding the metabolic responses of living systems to pathophysiological stimuli via multivariate statistical analysis of biological nmr spectroscopic data[END_REF]. To this end, it studies the concentration in small molecules called metabolites that are the end products of the enzymatic machinery of the cell. Indeed, minor variations in gene or protein expression levels that are not observable via high throughput experiments may have an influence on the metabolites and hence on the phenotype of interest. Thus, metabolomics is a promising approach that can advantageously complement usual transcriptomic and proteomic analyses.

In metabolomics, the analysis of the biological samples is often performed using High Resolution Mass Spectrometry (HRMS), Nuclear Magnetic Resonance (NMR) or Liquid Chromatography-Mass Spectrometry (LC-MS) and produces a large number of features (hundreds or thousands) that can explain a difference between two or more populations, see [START_REF] Zhang | Modern analytical techniques in metabolomics analysis[END_REF]. It is well-known that the identification of metabolites discriminating these populations remains a major bottleneck in metabolomics and therefore the selection of relevant features (variables) is a crucial step in the statistical analysis of the metabolomic data, as explained in [START_REF] Verdegem | Improved metabolite identification with midas and magma through ms/ms spectral dataset-driven parameter optimization[END_REF].

Different supervised machine learning approaches have been used in metabolomics during the last few years, see [START_REF] Saccenti | Reflections on univariate and multivariate analysis of metabolomics data[END_REF]; [START_REF] Ren | Computational and statistical analysis of metabolomics data[END_REF]; [START_REF] Boccard | Exploring omics data from designed experiments using analysis of variance multiblock orthogonal partial least squares[END_REF]. Among them the most widely used is the partial least squares-discriminant analysis (PLS-DA) which has recently been extended to sPLS-DA (sparse partial least squares-discriminant analysis) by Lê [START_REF] Lê Cao | Sparse pls discriminant analysis: biologically relevant feature selection and graphical displays for multiclass problems[END_REF] to include a variable selection step. Nevertheless, [START_REF] Grissa | Feature selection methods for early predictive biomarker discovery using untargeted metabolomic data[END_REF] highlight the need for new development in the process of features selection that would take into account the specificity of metabolomics data which is the dependence that may exist between the different metabolites. In this perspective, our paper proposes a novel feature selection methodology which consists in a variable selection approach based on the Lasso criterion in a multivariate setting taking into account the dependence that may exist between the different metabolites.

More precisely, let us consider a classical metabolomics experiment where n samples have been collected and analyzed. This results in an n × q data matrix where q stands for the number of metabolites. When the n samples have been obtained under various conditions, we are typically interested in understanding the effect of each condition on each metabolite. In the case where C experimental conditions are compared, n c denotes the number of replicates under condition c, where c ∈ {1, . . . C} and C c=1 n c = n. We further denote Y (j) c,r the centered LC-MS signal obtained for the jth metabolite (j ∈ {1, . . . q}) under Condition c for Replicate r (r ∈ {1, . . . n c }). In the following, the set of conditions will be called the "factor", each specific condition being a "level" of this factor. The most popular model to analyze quantitative observations Y as a function of a qualitative variable, that is a factor, is the analysis of variance (ANOVA) model, which we write here as follows:

Y (j) c,r = µ (j) c + E (j) c,r , (1)
where the observations {Y (j) c,r } are assumed to be centered, so that µ

(j)
c can be interpreted as the effect of Condition c (Level c) on Metabolite j and where the residual terms {E (j) c,r } are assumed to be independent and identically distributed (i.i.d.) zero-mean Gaussian random variables. The goal of such a modeling is to highlight which effects among the µ

(j) 1 , µ (j) 2 , . . . , µ (j) C
are the most significant for the metabolite j since the {Y (j) c,r } are assumed to be centered. When the whole n × q data matrix is considered instead of a single column j, the model can be summarized in the following matrix form:

Y = XB + E, (2)
where Y = (Y i,j ) 1≤i≤n, 1≤j≤q is the n × q observation matrix, X is the n × p design matrix, B is the p × q coefficient matrix and E = (E i,j ) 1≤i≤n, 1≤j≤q is the n × q matrix of residual errors. Observe that p corresponds to the number of explicative variables, which is simply C in Model (1). For notational simplicity, the indices c, r in Y (j) c,r are summarized in a unique index i in {1, . . . n}.

In this paper, we pay a special attention to the potential dependence that may exist among the columns of Y , namely the different metabolites. To this aim, we shall assume that for each i in {1, . . . , n},

(E i,1 , . . . , E i,q ) ∼ N (0, Σ q ), (3)
where Σ q denotes the covariance matrix of the ith row of the residual error matrix. Note that the model defined by ( 2) and ( 3) is usually called a general linear model or a multivariate linear model which has been extensively studied in [START_REF] Mardia | Multivariate analysis. Probability and mathematical statistics[END_REF].

The simplest assumption regarding the dependence structure of the noise is Σ q = σ 2 I q , where I q denotes the q×q identity matrix. In this case the different columns of Y are assumed to be independent. In more general cases, the matrix Σ q models the dependence between the different columns of Y , namely the dependence between the metabolites. In the following, we shall moreover assume that (E i,1 , . . . , E i,q ) and (E k,1 , . . . , E k,q ) are independent, when i = k, which means that the individuals are assumed to be independent.

The problem of finding which parameters are significant among the (µ

(j)
c ) 1≤c≤C,1≤j≤q in Model (1) boils down to finding the non null coefficients in the matrix B in Model (2) and hence can be seen as a variable selection problem in the general linear model. Several approaches can be considered for solving this issue: either a posteriori methods such as classical statistical tests in ANOVA models, see [START_REF] Mardia | Multivariate analysis. Probability and mathematical statistics[END_REF]; [START_REF] Faraway | Linear Models with R[END_REF] or methods embedding the variable selection such as Lasso-type methodologies initially proposed by [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]. However, a raw application of such approaches does not take into account the potential dependence between the different columns of Y . This drawback will be illustrated in Section 3.

The goal of our paper is twofold: First, to remedy the limitations of these approaches by proposing a method for estimating the dependence between the columns of Y and second, to deal with the potentially high number of variables by using a Lasso-type approach taking into account this dependence. For this purpose, we shall propose a three-step inference strategy further detailed hereafter.

The paper is organized as follows. Our method is described in Section 2. To support our methodology, some numerical experiments on synthetic data are provided in Section 3. Finally, an application to a metabolomics data set produced by the analysis of African copals samples is given in Section 4.

Statistical inference

The strategy that we propose can be summarized as follows.

• First step: Fitting a one-way ANOVA to each column of the observation matrix Y in order to have access to an estimation E of the residual matrix E. • Second step: Estimating the matrix Σ q by using the methods described in Sections 2.1.1 and 2.1.2. Then, choosing the most convenient estimator Σ q thanks to a statistical test described in Section 2.1.3. • Third step: Thanks to the estimator Σ q , transforming the data in order to remove the dependence between the columns of Y . Such a transformation will be called "whitening" hereafter. Then, applying to these transformed observations the Lasso approach described in Section 2.2.

The first step provides a first estimate B of B. An estimate of the residual matrix is then defined as E = Y -X B. In the following, we shall focus on the last two steps.

2.1. Estimation of the dependence structure of E. We propose hereafter to model each row of E as realizations of a stationary process and hence we shall use time-series models in order to describe the dependence structure of E. We refer the reader to [START_REF] Brockwell | Time Series: Theory and Methods[END_REF] for further details on time series modeling.

We shall consider hereafter a large variety of models ranging from the simplest parametric to the most general nonparametric dependence modeling. In each case we focus on the estimation of Σ -1/2 q since using the following transformation:

Y Σ -1/2 q = XB Σ -1/2 q + E Σ -1/2 q (4)
removes the dependence between the columns of Y . Such a procedure will be called "whitening" hereafter.

2.1.1. Parametric dependence. The simplest model among the parametric models is the autoregressive process of order 1 denoted AR(1). More precisely, for each i in {1, . . . , n}, E i,t satisfies the following equation:

E i,t -φ 1 E i,t-1 = W i,t , with W i,t ∼ W N (0, σ 2 ), (5)
where φ 1 is a real number and W N (0, σ 2 ) denotes a zero-mean white noise process of variance

σ 2 , namely, if Z t ∼ W N (0, σ 2 ), then E(Z t ) = 0, E(Z t Z t ) = 0 if t = t and E(Z 2 t ) = σ 2 .
In this case, the inverse of the square root of the covariance matrix Σ q of (E i,1 , . . . , E i,q ) has a simple closed-form expression given by

Σ -1/2 q =         1 -φ 2 1 -φ 1 0 • • • 0 0 1 -φ 1 • • • 0 0 0 . . . . . . . . . . . . . . . . . . . . . -φ 1 0 0 • • • 0 1         . ( 6
)
Hence, to obtain the expression of Σ -1/2 q , it is enough to have an estimation of the parameter φ 1 and to replace it in (6). For this, we use the estimator E of the residual errors matrix obtained by fitting a standard ANOVA model to the observations, which corresponds to the first step of our method. Then φ 1 is estimated by φ 1 defined by

φ 1 = 1 n n i=1 φ 1,i ,
where φ 1,i denotes the estimator of φ 1 obtained by the classical Yule-Walker equations from ( E i,1 , . . . , E i,q ), see [START_REF] Brockwell | Time Series: Theory and Methods[END_REF] for more details. More generally, it is also possible to have access to Σ -1/2 q for more complex processes such as the ARMA(p, q) process defined as follows: For each i in {1, . . . , n},

E i,t -φ 1 E i,t-1 -• • • -φ p E i,t-p = W i,t + θ 1 W i,t-1 + . . . θ q W i,t-q , (7)
where W i,t ∼ W N (0, σ 2 ), the φ i 's and the θ i 's are real numbers.

2.1.2. Nonparametric dependence case. In the situation where a parametric modeling is not relevant for Σ q , it can be estimated by

Σ q =      γ(0) γ(1) • • • γ(q -1) γ(1) γ(0) • • • γ(q -2) . . . γ(q -1) γ(q -2) • • • γ(0)      , (8)
with

γ(h) = 1 n n i=1 γ i (h),
where γ i (h) is the standard autocovariance estimator of γ i (h) = E(E i,t E i,t+h ), for all t. Usually, γ i (h) is referred to as the empirical autocovariance of the E i,t 's at lag h (i.e. the empirical covariance between ( E i,1 , . . . , E i,n-h ) and ( E i,h+1 , . . . , E i,n )). For a definition of the standard autocovariance estimator we refer the reader to Chapter 7 of [START_REF] Brockwell | Time Series: Theory and Methods[END_REF].

The matrix Σ -1/2 q is then obtained by inverting the Cholesky factor of Σ q .

2.1.3. Choice of the whitening modeling. In order to decide which dependence modeling is the most adapted to the data at hand we propose hereafter a statistical test. If the whitening modeling used is well chosen then each row of

E = E Σ -1/2 q
should be a white noise. One of the most popular approach for testing whether a random process is a white noise is the Portmanteau test which is based on the Bartlett theorem, for further details we refer the reader to (Brockwell and Davis, 1991, Theorem 7.2.2). By this theorem, we get that under the null hypothesis (H 0 ): "For each i in {1, . . . , n}, ( E i,1 , . . . , E i,q ) is a white noise",

q H h=1 ρ i (h) 2 ≈ χ 2 (H), as q → ∞, (9)
for each i in {1, . . . , n}, where ρ i (h) denotes the empirical autocorrelation of ( E i,t ) t at lag h and χ 2 (H) denotes the chi-squared distribution with H degrees of freedom. Thus, by ( 9), we have at our disposal a p-value for each i in {1, . . . , n} that we denote by Pval i . In order to have a single p-value instead of n, we shall consider

q n i=1 H h=1 ρ i (h) 2 ≈ χ 2 (nH), as q → ∞, (10)
where the approximation comes from the fact that the rows of E are assumed to be independent. Equation ( 10) thus provides a p-value: Pval. Hence, if Pval < α, the null hypothesis (H 0 ) is rejected at the level α, where α is usually equal to 5%. In such a situation taking the dependence into account and estimating the dependence by one of the previous methods should highly improve the modeling and the variable selection step.

Estimation of B.

2.2.1. Lasso based approach. Let us first explain briefly the usual framework in which the Lasso approach is used. We consider a high-dimensional linear model of the following form

Y = X B + E, (11)
where Y, B and E are vectors. Note that, in high-dimensional linear models, the matrix X has usually more columns than rows which means that the number of variables is larger than the number of observations but B is usually a sparse vector, namely it contains a lot of null components.

In such models a very popular approach initially proposed by [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] consists in using the Least Absolute Shrinkage eStimatOr (LASSO) criterion for estimating B defined as follows for a positive λ:

B(λ) = Argmin B Y -X B 2 2 + λ B 1 , (12)
where, for u = (u 1 , . . . , u n ), u 2 2 = n i=1 u 2 i and u 1 = n i=1 |u i |, which is usually called the 1 -norm of the vector u. Observe that the first term of ( 12) is the classical least-squares criterion and that λ B 1 can be seen as a penalty term. The interest of such a criterion is the sparsity enforcing property of the 1 -norm ensuring that the number of non-zero components of the estimator B of B is small for large enough values of λ. Such a criterion is very relevant in our framework since the problem of finding the significant variables boils down to finding the non null coefficients in the matrix B.

This methodology cannot be directly applied to our model since we have to deal with matrices and not with vectors. However, as explained in Appendix A, Model (2) can be rewritten as in (11) where Y, B and E are vectors of size nq, pq and nq, respectively. Hence, retrieving the positions of the non null components in B is a first approach for finding relevant variables.

However, this approach does not take into account the dependence between the columns of Y . Hence, we propose hereafter a modified version of the standard Lasso criterion ( 12) taking into account this potential dependence.

As explained previously, our contribution consists first in "whitening" the observations, namely removing the dependence that may exist within the observations matrix, by multiplying (2) on the right by

Σ -1/2 q , see (4) where Σ -1/2 q is replaced by Σ -1/2 q
. By using the same vectorization trick that allows us to transform Model (2) into Model ( 11), the Lasso criterion can be applied to the vectorized version of Model (4) where

Σ -1/2 q is replaced by Σ -1/2 q
. The specific expressions of Y, X , B and E are given in Appendix B.

Note that this idea of "whitening" the observations has also been proposed by [START_REF] Rothman | Sparse multivariate regression with covariance estimation[END_REF] in which the estimation of Σ q and B is performed simultaneously. An implementation is available in the R package MRCE. In our approach, Σ q is estimated first and then its estimator is used in (4) instead of Σ q before applying the Lasso criterion. Hence, our method can be seen as a variant of the MRCE method in which Σ q is estimated beforehand. Moreover, after some numerical experiments, we observed that for the values of n and q that we aim at using, the computational burden of the approach designed by [START_REF] Rothman | Sparse multivariate regression with covariance estimation[END_REF] is too high for addressing our datasets, contrary to ours. Hence, in the following, we shall not consider the method of [START_REF] Rothman | Sparse multivariate regression with covariance estimation[END_REF] anymore.

2.2.2. Model selection issue. We can see that the estimator defined in (12) depends on a parameter λ which tunes the sparsity level in B. We propose to mix two standard approaches to estimate the positions of the non null components in B.

We first use the 10-fold cross-validation method to choose the λ denoted λ CV minimizing the cross-validation criterion. This λ is then used in the stability selection approach of [START_REF] Meinshausen | Stability selection[END_REF] which guarantees the robustness of the selected variables. This latter approach can be described as follows. The vector of observations Y is randomly split into several subsamples of size nq/2 which is possible thanks to the whitening step. For each subsample, the LASSO criterion is applied with λ = λ CV and the indices i of the non null B i are stored. Then, for a given threshold, we keep in the final set of selected variables only the variables appearing a number of times larger than this threshold. In practice, we generated 5000 subsamples of Y.

Concerning the choice of the final threshold: we propose either to take the one leading to the largest p-value of the whitening test described in (10) or the threshold 1. As we shall see in Section 3, with the first choice, mostly all the positions of the non null variables in B are retrieved with some false positive. With the second choice, all the true positions are not recovered but there are no false positive. Moreover, the second choice guarantees a stability of the selected variables since only the variables which are chosen at each of the 5000 subsamplings of the data are finally kept.

Simulation study

The goal of this section is to assess the statistical performance of our methodology. In order to emphasize the benefits of using a whitening approach from the variable selection point of view, we shall first compare our approach to standard methodologies. Then, we shall analyze the performance of our statistical test for choosing the best dependence modeling. Finally, we shall investigate the performance of our model selection criterion.

To assess the performance of these different methodologies, we generated observations Y according to Model (2) with q = 1000, p = 3, n = 30 and different dependence modelings, namely different matrices Σ q corresponding to the AR(1) model described in (5) with σ = 1 and φ 1 = 0.7 or 0.9.

Note that we have chosen the values of the parameters p, q and n in order to be as close as possible to the real data that we plan to analyze in Section 4.

We shall also investigate the effect of the sparsity and of the signal to noise ratio. In the following, the sparsity level corresponds to the number of non null elements in B divided by the total number nq of elements of B. Different signal to noise ratios are obtained by multiplying B in (2) by a coefficient κ.

3.1.

Variable selection performance. The goal of this section is to compare the performance of our different whitening strategies presented above to standard existing methodologies. More precisely, we shall compare our approaches to the classical ANOVA method (denoted ANOVA), the standard Lasso (denoted Lasso), namely the Lasso approach without the whitening step and to sPLSDA devised by Lê Cao et al. (2011) and implemented in the mixOmics R package, which is widely used in the metabolomics field. By ANOVA, we mean the classical one-way ANOVA applied to each column of the observations matrix Y without taking the dependence into account.

In the following, the different whitening approaches that we propose will be denoted by AR1 and Nonparam. They are described in Sections 2.1.1 and 2.1.2, respectively. These methods will also be compared to the Oracle approach which assumes that the matrix Σ q is known, which is never the case in practical situations.

For comparing these different methods, we shall use two classical criteria: ROC curves and AUC (Area Under the ROC Curve). ROC curves display the true positive rates as a function of the false positive rates and the closer to one the AUC the better the methodology. Since sPLSDA only selects relevant metabolites but does not assign them to a level of the factor, we shall consider that as soon as a relevant metabolite is selected it is a true positive which gives an advantage to sPLSDA.

We can see from Figure 1 that in the case of an AR(1) dependence, taking into account this dependence provides better results than sPLSDA and than approaches that consider the columns of the residual matrix as independent. Moreover, we observe that the performance of the non parametric modeling are on a par with those of the parametric and the oracle κ = 2 FPR Figure 1. Means of the ROC curves obtained from 200 replications for the different methodologies in the AR(1) dependence modeling; κ is linked to the signal to noise ratio (first row: κ = 1, second row κ = 2); φ 1 is the correlation level in the AR(1) and s the sparsity level (i.e. the fraction of nonzero elements) in the vector of true parameters. ones. We also note that the larger the sparsity level the smaller the difference of performance between the different approaches. However, the larger the signal to noise ratio the better the performance of the different methodologies.

3.2. Choice of the dependence modeling. The goal of this section is to assess the performance of the dependence modeling strategy that we proposed in Section 2.1.3. We generated observations Y with the parameters described at the beginning of Section 3 in the case of an AR(1) dependence, for a sparsity level of 0.01 and when κ = 1. The corresponding results are displayed in Figure 2.

We observe from this figure that our test provides p-values close to zero in the case where no whitening strategy is used (Lasso) and that when one of the proposed whitening approaches is used the p-values are larger than 0.7.

Choice of the model selection criterion.

We investigate hereafter the performance of our model selection criterion described in Section 2.2.2.

Figure 3 displays the means of the p-values of the test described in 2.1.3 obtained from 5000 replications of the observations Y generated with the parameters described at the beginning of Section 3 in the case of an AR(1) dependence with φ 1 = 0.9 and κ = 1. q q q q q q q q AR1 nonparam no_whitening We observe from this figure that the p-values are all the more high that the thresholds are large.

Figure 4 displays with bullets ('•') the positions of the variables selected by our three-step approach for the two possible choices of thresholds from 500 replications of Y obtained with the parameters described at the beginning of Section 3 in the case of an AR(1) dependence with φ 1 = 0.9 and κ = 10.

We observe from this figure that mostly all the positions of the non null variables in B are retrieved with some false positive when the threshold is obtained by maximizing the p-value. When the threshold is equal to 1, there are no false positive but all the true positions are not recovered.

3.4. Numerical performance. In order to investigate the computational burden of our approach, we generated matrices Y satisfying Model (2) with n = 30 and q ∈ {100, 200, . . . , 1000}. q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 2 3 250 500 750

Position

Level qq q q q qq q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q qqqq q q q q q q q q q qqqq q qq q q q q q q q q q q q q q q q q q q q q qq q q q q q qqq q q q q q q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q qqqq q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq qqq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q qq q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q qqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqq q q q q qq q q q q q q q qqqq q q q q q q q q q q q q q q q q q q qqqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q qqqqqq qqq qqqqq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q qqqq q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqq q q q q q q q q q q qqq q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q qq q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qq q q q qq q q q q q q q qq q q q q q q q qq q q q q q q q q q q q q q qqqq q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q qqq q q q q qqqqqqqq q q q qqq q q q q q q q q q q q qqq q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q qqq q q q q qq q q q q q qqqq q q q q q q q q qqq q q q q q q q q q q q q q q q q qqq q q q q q q q q q q qqq q q q q q q q q qq q q q qqq q q q q qqqq qqq q q q q q q q q q q q q q q q q qqqq qqqq q q q qqqq q q q q q q q q q q qqq qqq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q qqq q q q q q q qqq q q q q q q q q q q q q q q q q qqqqqqqq q q qq q q q q q q q qqq q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q qq q q q q qqqq q q q q q q q q q q q q q q q q q q q qqqq q q q q q q q q q q q qq q q q q q q q qq q q q q q q q q q q q q q q q qqq qqq q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q qq q q q q q q qqq q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q qqqqq q q q q q q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q q q qqq q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqq q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqq qqq q q q q qqq q q q q q q qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q qq q q qq q q q q q qq q q q q qq q q qqq q q q q q q q q q q q q q q q q qqqq q q q q qq q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q qqqq qq q qqqq q q q q q q q q q q q q q q q q q q q q q qqqqq qqq q q q q q qqqq q q q q q q q q q q qqqq q q q q q qqq q q q q q q q q qqq q q q qqqqqq qq q q q q qqqqqq q q q q q q qqq q q qq q q q q q q q q q qq q qqq qqq q q q q q q q q qq q q q q qqq qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q qq q q q qq q q q q q q q q qqq q q q q q qqqqqqqqqq qqqqq q q q q qq q q q q qqq q q q q q q q q q q qqq q q q q qqqqq q q q q q q q q q qq q q q q qqq q qq q q q q q q qqq qqq q q q q q q q q q q q q q q q q q q q q qqqq qqqqqqqqqq q qqq q q q q q q q q q q q q q q qq q q q q qq q qqq qqq q q q q q q qq q qqq q q q q q qqqq q q q qqq q qqq q q q qqqqqqqq q q q q q qqq q q q q q q qqq q qqqqqq qqqqqq q q q q q q q q q q q q q q q q qqq q qqq qqqqqq q q q q q q q q q q q q q qqq q q q q q q q q q q q q q qqq q q q q qqqq qqqqq q q qq qqqq q q q qq qqqqq q q q q q q q q q qqqqq q q qq q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q qqqq q q q q q qq qq q q q q q q q q q qqqqqqqq q q q q q q q q q q q q q q qqq q qq q q q q q q q q qqq q q q qq q q q qqq q qqq q q q q q q q q q q q q q q q qqqq qq q q q q q q q q qqqq q q q q q q q qq q q q q q q qqq q q q qqqqqq q q q q qqqq qqq qqqq q q q q qqq q q q q q q q q qq q q q q q q qq q qq qq qq q q q q q qqqqq q q q q q q q qq q q q q q q q q q q qq qqqqqq q q q q qqq q q q q q qq q q q q q q q q q q q q qqq qqq qqq q q q q q qqqqqqq q q q q qqq q q q q q q q qq q q q qq q q q qq q q q q q q q q qqq q q q q q qq q q q q q q q q qqq q q q q q q qq q q qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq 1 2 3 0 250 500 750 1000

Position

Level Frequency q q q q 0.25 0.50 0.75 1.00 q Estimate True Figure 4. Positions of the variables selected by our approach ('•') when κ = 10. Values on the y-axis correspond to the 3 levels of the factor. The results obtained when the threshold is equal to 1 are on the left and the results when the threshold is obtained by maximizing the p-value are on the right. The size of the bullets are all the more large that the number of times where a variable has been selected is large. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 2 3 Here, the rows of the matrix E are generated as realizations of an AR(1) process and the level of sparsity s of B is equal to 0.01.

Figure 5 displays the computational times of MultiVarSel obtained from a computer having the following configuration: RAM 16 GB, CPU 8×3.6 GHz for different number of replications in the stability selection stage. We can see from this figure that the computational burden of MultiVarSel is very low and that it takes only a few seconds to analyze matrices having 1000 columns.

Application to a LC-MS dataset

In this section, our three-step methodology implemented in the R package MultiVarSel and available from the CRAN, is applied to a LC-MS (Liquid Chromatography-Mass Spectrometry) data set made of African copals samples. The samples correspond to ethanolic extracts of copals produced by trees belonging to two genera Copaifera (C) and Trachylobium (T) with a second level of classification coming from the geographical provenance of the Copaifera samples (West (W) or East (E) Africa). Since all the Trachylobium samples come from East Africa, we have a single factor having three levels: CE, CW and TE such that n CE = 9, n CW = 8 and n TE = 13.

In this section, we also compare the performance of our method with those of other techniques which are widely used in metabolomics.

4.1. Data pre-processing. LC-MS chromatograms were aligned using the R package XCMS proposed by [START_REF] Smith | XCMS: Processing mass spectrometry data for metabolite profiling using Nonlinear peak alignment, matching, and identification[END_REF] with the following parameters: a signal to noise ratio threshold of 10:1 for peak selection, a step size of 0.2 min and a minimum difference in m/z for peaks with overlapping retention times of 0.05 amu. Sample filtering was also performed: To be considered as informative, as suggested by [START_REF] Kirwan | Characterising and correcting batch variation in an automated direct infusion mass spectrometry (dims) metabolomics workflow[END_REF], a peak was required to be present in at least 80% of the samples. Missing values imputation was realized using the KNN algorithm described in [START_REF] Hrydziuszko | Missing values in mass spectrometry based metabolomics: an undervalued step in the data processing pipeline[END_REF]. Subsequently, the spectra were normalized to equalize signal intensities to the median profile in order to reduce any variance arising from differing dilutions of the biological extracts and probabilistic quotient normalization (PQN) was used, see [START_REF] Dieterle | Probabilistic quotient normalization as robust method to account for dilution of complex biological mixtures. application in 1h nmr metabonomics[END_REF] for further details. In order to reduce the size of the data matrix, selection of the adducts of interest [M+H] + was then performed using the CAMERA package of [START_REF] Kuhl | Camera: an integrated strategy for compound spectra extraction and annotation of liquid chromatography/mass spectrometry data sets[END_REF]. A n × q matrix Y was then obtained and submitted to the statistical analyses. 4.2. Application of our three-step approach. The observations matrix Y is first centered and scaled in order to ensure that the empirical mean in each column is 0 and that the empirical variance is 1. 4.2.1. First step. A one-way ANOVA is fitted to each column of the observation matrix Y in order to have access to an estimation E of the residual matrix E. Then, the test proposed in Section 2.1.3 is applied. We found a p-value equal to zero which indicates that the columns of E cannot be considered as independent and hence that applying the whitening strategy should improve the results. 4.2.2. Second step. The different whitening strategies described in Section 2.1 were applied and the highest p-value for the test described in Section 2.1.3 is obtained for the nonparametric whitening. More precisely, the p-values obtained for the AR(1) and the nonparametric dependence modeling are equal to 1.5 × 10 -4 and 0.5107, respectively. Hence, in the following we shall use the nonparametric modeling. 4.2.3. Third step. The Lasso approach described in Section 2.2 was then applied to the whitened observations where Σ q is obtained by using the nonparametric modeling. The stability selection is then used with 5000 replications and a threshold equal to 1 in order to avoid false positive.

The Venn diagram of Figure 6 displays the repartition of the selected metabolites among the different classes CE, TE and CW. We can see from this figure that at least one metabolite is selected as a marker for each class (20 for TE, 22 for CW and 1 for CE) for a total of 39 unique metabolites. More precisely, our methodology leads to a list of metabolites that mainly characterize a single class. As recommended by Lê Cao et al. ( 2011), we used two components for PLS-DA and sPLS-DA. Moreover, in order to make sPLS-DA comparable with our approach, 20 variables are kept for each component in the sPLS-DA methodology. The corresponding results are displayed in Figure 7. We can see from this figure that sPLS-DA exhibits better classification performance than the standard PLS-DA.

Since PLS-DA does not include a variable selection step we shall compare our approach only to sPLS-DA in the following. For comparing these methodologies Figure 8 displays the PCA obtained when all the metabolites are kept on the one hand and when the metabolites are those selected by sPLS-DA or by our methodology on the other hand. We can see from this figure that, one the hand, the approaches containing a variable selection step exhibit better classification performance and that, on the other hand, sPLS-DA and our method show similar performance from the classification point of view even if our approach is not designed for this purpose.

Figure 9 displays the positions of the metabolites selected by our approach and sPLS-DA. We can see from this figure that out of the 39 selected metabolites, 6 metabolites are selected by both sPLS-DA and our methodology. The major difference between these two variable selection techniques is that our method selects metabolites having a ratio m/z smaller than 300 whereas the metabolites chosen by sPLS-DA lie within the range 300-400 m/z.

In order to further compare our methodology with sPLS-DA, we first propose to assess the stability of the selected variables (or metabolites). For this purpose, we performed 10 bootstrap resamplings of our original data and we compared the variables selected by both approaches. The results are displayed in Figure 10 and in Table 1. Figure 10 displays the frequencies at which each metabolite has been selected by the two methods. We can see from this figure that the highest selection frequency of sPLS-DA is around 0.8 and that a lot of variables have a selection frequency smaller than 0.5. Moreover, we can see from Table 1 which provides the number of metabolites which have been selected once (first row), twice X-variate 1: 26% expl. var X-variate 2: 23% expl. var Legend q q q q q q q q q CE CW TE (second row)..., that our approach selects 4 metabolites with a frequency equal to 1 which does not occur for sPLS-DA. Hence, from this point of view, our approach is more stable than sPLS-DA. PC1: 43% expl. var PC2: 15% expl. var Legend q q q q q q q q q CE CW TE Finally, we compare hereafter the set of variables provided by sPLS-DA and our approach from the classification error point of view. Since our method is not designed for yielding a classification we give the selected variables to PLS-DA in order to obtain such a classification. The estimation of the classification error rates are then obtained by using a 10-fold crossvalidation. The corresponding results are displayed in Tables 2 and3. We observe that the classification error rates of our approach are on a par with those of sPLS-DA.

Nb of selection

q qq qq q qqq q q q qq q q qq q q q q q q qq q q q q q q q q q q q q q q MultiVarSel sPLSDA 200 300 400 500 600 700 m/z TE CE CW TE 0.92 0.33 0.00 CE 0.08 0.67 0.00 CW 0.00 0.00 1.00 Table 2. Classification error rates for sPLS-DA. TE CE CW TE 0.92 0.22 0.00 CE 0.08 0.67 0.00 CW 0.00 0.11 1.00 Table 3. Classification error rates for our approach coupled with PLS-DA.

We observe from these different investigations that our approach provides similar results as sPLS-DA in terms of classification even if our approach was not designed for this purpose and that it yields more stable variables (metabolites) than sPLS-DA for characterizing the different classes.

Conclusion

In this paper, we proposed a novel approach for analyzing LC-MS metabolomics data by introducing a new Lasso-type approach taking into account the dependence that may exist between the columns of the data matrix. Our approach is implemented in the R package Mul-tiVarSel which is available from the The Comprehensive R Archive Network (CRAN). In the course of this study, we have shown that our method has two main features. Firstly, it is very efficient from a statistical point of view for selecting a restricted number of stable metabolites characterizing each level of the factor of interest. Secondly, its very low computational burden makes its use possible on very large LC-MS metabolomics data.

Figure 2 .Figure 3 .

 23 Figure2. Means and standard deviations of the p-values of the test described in Section 2.1.3 for the different approaches in the AR(1) dependence modeling when φ 1 = 0.7 (left) and φ 1 = 0.9 (right).

Figure 5 .

 5 Figure 5. Computational times (in seconds) of MultiVarSel.

Figure 6 .

 6 Figure 6. Venn diagram of the metabolites selected for each class by Multi-VarSel using a threshold equal to 1 in the stability selection stage.

  1: 22% expl. var X-variate 2: 21% expl. var

Figure 7 .

 7 Figure 7. 2D scores plot of the PLS-DA and the sPLS-DA.

Figure 8 .

 8 Figure8. PCA with all the metabolites and with the metabolites selected by sPLS-DA and our approach MultiVarSel.

Figure 9 .Figure 10 .

 910 Figure 9. Comparison of the metabolites selected by our approach Multi-VarSel and by sPLS-DA.

Table 1 .

 1 Number of times the different metabolites have been selected by sPLS-DA and our approach MultiVarSel.

		Nb of selected metabolites Nb of selected metabolites
		by sPLS-DA	by MultiVarSel
	1	117	143
	2	41	54
	3	26	24
	4	4	20
	5	8	15
	6	5	6
	7	3	6
	8	2	3
	9	0	3
	10	0	4

Acknowledgements

This project has been funded by La mission pour l'interdisciplinarité du CNRS in the frame of the DEFI ENVIROMICS (project AREA). The authors thank the Musée François Tillequin for providing the samples from the Guibourt Collection.

where we used that (Mardia et al., 1979, Appendix A.2.5)

In this equation, B denotes the transpose of the matrix B. Thus,

where X = I q ⊗ X and Y, B and E are vectors of size nq, pq and nq, respectively.

Appendix B. Let us apply the vec operator to Model (4) where

UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005, Paris, France E-mail address: marie.perrot-dockes@agroparistech.fr UMR CNRS 8638 Comète -Université Paris-Descartes, CNRS, 75006 Paris France.