M Perrot-Dock
email: marie.perrot-dockes@agroparistech.fr

C L Ès

L Évy-Leduc

Rajjou

ESTIMATION OF LARGE BLOCK STRUCTURED COVARIANCE MATRICES: APPLICATION TO "MULTI-OMIC" APPROACHES TO STUDY SEED QUALITY

Estimation of large block structured covariance matrices: Application to "multi-omic" approaches to study seed quality

Introduction

Plant functional genomics refers to the description of the biological function of a single or a group of genes and both the dynamics and the plasticity of genome expression to shape the phenotype. Combining multi-omics such as transcriptomic, proteomic or metabolomic approaches allows us to address in a new light the dimension and the complexity of the different levels of gene expression control and the delicacy of the metabolic regulation of plants under fluctuation environments. Thus, our era marks a real conceptual shift in plant biology where the individual is no longer considered as a simple sum of components but rather as a system with a set of interacting components to maximize its growth, its reproduction and its adaptation. Plant systems biology is therefore defined by multidisciplinary and multi-scale approaches based on the acquisition of a wide range of data as exhaustive as possible.

In this context, it is crucial to propose new methodologies for integrating heterogeneous data explaining the co-regulations/co-accumulations of products of gene expression (mRNA, proteins) and metabolites. In order to better understand these phenomena, our goal will thus be to propose a new approach for estimating block structured covariance matrix in a highdimensional framework where the dimension of the covariance matrix is much larger than the sample size. In this setting, it is well known that the commonly used sample covariance matrix performs poorly. In recent years, researchers have proposed various regularization techniques to consistently estimate large covariance matrices or the inverse of such matrices, namely precision matrices. To estimate such matrices, one of the key assumptions made in the literature is that the matrix of interest is sparse, namely many entries are equal to zero. A number of regularization approaches including banding, tapering, thresholding and 1 minimization, have been developed to estimate large covariance matrices or their inverse such as, for instance, [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], [START_REF] Bickel | Covariance regularization by thresholding[END_REF], [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF], [START_REF] Bien | Sparse estimation of a covariance matrix[END_REF] and [START_REF] Rothman | Positive definite estimators of large covariance matrices[END_REF] among many others. For further references, we refer the reader to [START_REF] Cai | Adaptive covariance matrix estimation through block thresholding[END_REF] and to the review of [START_REF] Fan | An overview of the estimation of large covariance and precision matrices[END_REF].

In this paper, we shall consider the following framework. Let E 1 , E 2 , • • • , E n , n zero-mean i.i.d. q-dimensional random vectors having a covariance matrix Σ such that the number q of its rows and columns is much larger than n. The goal of the paper is to propose a new estimator of Σ and of the square root of its inverse, Σ -1/2 , in the particular case where Σ is assumed to have a block structure without limiting ourselves to diagonal blocks. An accurate estimator of Σ can indeed be very useful to better understand the links between the columns of the observation matrix and may highlight some biological processes. Moreover, an estimator of Σ -1/2 can be very useful in the general linear model in order to remove the dependence that may exist between the columns of the observation matrix. For further details on this point, we refer the reader to Perrot-Dockès et al. (2018), Perrot-Dockès et al. (2018) and to the R package MultiVarSel in which such an approach is proposed and implemented for performing variable selection in the multivariate linear model in the presence of dependence between the columns of the observation matrix.

More precisely, in this paper, we shall assume that

(1)

Σ = ZZ + D,
where Z is a q × k sparse matrix with k q, Z denotes the transpose of the matrix Z and D is a diagonal matrix such that the diagonal terms of Σ are equal to one. Two examples of such matrices Z and Σ are given in Figure 1 in the case where k = 5 and q = 50 and in the case where the columns of Σ do not need to be permuted in order to see the block structure. Based on (1), our model could seem to be close to factor models described in Johnson and Wichern (1988) and [START_REF] Fan | An overview of the estimation of large covariance and precision matrices[END_REF]. However, in Johnson and Wichern (1988), the high-dimensional aspects are not considered and in [START_REF] Fan | An overview of the estimation of large covariance and precision matrices[END_REF] the sparsity constraint is not studied. [START_REF] Blum | Sparse factor model for co-expression networks with an application using prior biological knowledge[END_REF] proposed a methodology which is based on the factor model but with a sparsity constraint on the coefficients of Z which leads to a sparse covariance matrix. Note also that the block diagonal assumption has already been recently considered by [START_REF] Devijver | Block-diagonal covariance selection for high-dimensional gaussian graphical models[END_REF] for estimating the inverse of large covariance matrices in high-dimensional Gaussian Graphical Models (GGM).

We also propose a methodology to estimate Σ in the case where the block structure is latent; that is, permuting the columns and rows of Σ renders visible its block structure. An example of such a matrix Σ is given in Figure 2 in the case where k = 5 and q = 50.

Our approach is fully data-driven and consists in providing a low rank matrix approximation of the ZZ part of Σ and then in using a 1 regularization to obtain a sparse estimator of Σ. When the block structure is latent, a hierarchical clustering step must be applied first. With this estimator of Σ, we explain how to obtain an estimator of Σ -1/2 .

Our methodology is described in Section 2. Some numerical experiments on synthetic data are provided in Section 3. An application to the analysis of "-omic" data to study seed quality is performed in Section 4.

Statistical inference

The strategy that we propose for estimating Σ and Σ -1/2 can be summarized as follows.

• First step: Low rank approximation. In this step, we propose to approximate the part ZZ of Σ by a low rank matrix using a Singular Value Decomposition (SVD). • Second step: Detecting the position of the non null values. In this step, we use a Lasso criterion to yield a sparse estimator Σ of Σ. • Third step: Positive definiteness. We apply the methodology of [START_REF] Higham | Computing the nearest correlation matrix -a problem from finance[END_REF] to Σ to ensure that the final estimator Σ of Σ is positive definite. • Fourth step: Estimation of Σ -1/2 . In this step, Σ -1/2 is estimated from the spectral decomposition of Σ obtained in the previous step.

2.1. Low rank approximation. By definition of Z in (1), ZZ is a q × q low rank matrix having its rank smaller or equal to k q. In the first step, our goal is thus to propose a low rank approximation of an estimator of ZZ .

Let S be the sample q × q covariance matrix defined by

S = 1 n -1 n i=1 E i -E E i -E , with E = 1 n n i=1 E i ,
where E i = (E i,1 , . . . , E i,q ) . The corresponding q × q sample correlation matrix R = (R i,j ) is defined by:

(2)

R i,j = S i,j σ i σ j , ∀1 ≤ i, j ≤ q,
where

σ 2 i = 1 n -1 n =1 (E ,i -E i ) 2 , with E i = 1 n n =1 E ,i , ∀1 ≤ i ≤ q.
Let us also consider the (q -1) × (q -1) matrix Γ defined by:

Γ i,j = R i,j+1 , ∀1 ≤ i ≤ j ≤ q -1, (3) Γ i,j = Γ j,i , ∀1 ≤ j < i ≤ q -1.
If S was the real matrix Σ, the corresponding matrix Γ would have a rank less than or equal to k. Since S is an estimator of Σ, we shall use a rank r approximation Γ r of Γ. This will be performed by considering in its singular value decomposition only the r largest singular values and by replacing the other ones by 0. By [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF], this corresponds to the best rank r approximation of Γ. The choice of r will be discussed in Section 2.5.

2.2.

Detecting the position of the non null values. Let us first explain the usual framework in which the Lasso approach is used. We consider a linear model of the following form (4)

Y = X B + E,
where Y, B and E are vectors and B is sparse meaning that it has 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] is the Least Absolute Shrinkage eStimatOr (Lasso), which is defined as follows for a positive λ:

(5)

B(λ) = Argmin B Y -X B 2 2 + λ B 1 , where, for u = (u 1 , . . . , u n ), u 2 2 = n i=1 u 2 i and u 1 = n i=1 |u i |, i.e.
the 1 -norm of the vector u. Observe that the first term of ( 5) 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 λ.

Let (6) Y = vec H (Γ r ),
where vec H defined in Section 16.4 of [START_REF] Harville | Matrix Algebra: Exercises and Solutions: Exercises and Solutions[END_REF] is such that for a n × n matrix A,

vec H (A) =      a 1 * a 2 * . . . a n *     
, where a i * is the sub-vector of the column i of A obtained by striking out the i -1 first elements. In order to estimate the sparse matrix ZZ , we need to propose a sparse estimator of Γ r . To do this we apply the Lasso criterion described in (5), where X is the identity matrix.

In the case where X is an orthogonal matrix it has been shown in [START_REF] Giraud | Introduction to High-Dimensional Statistics[END_REF] that the solution of ( 5) is:

B(λ) j = X j Y(1 -λ 2|X j Y| ), if |X j Y| > λ 2 0,
otherwise, where X j denotes the jth column of X . Using the fact that X is the identity matrix we get

(7) B(λ) j = Y j (1 -λ 2|Y j | ), if |Y j | > λ 2 0, otherwise.
We then reestimate the non null coefficients using the least-squares criterion and get:

(8) B(λ) j = Y j , if |Y j | > λ 2 0, otherwise,
where Y is defined in (6). It has to be noticed that Γ r obtained in (7) satisfies the following criterion:

Γ r = Argmin Θ { Γ r -Θ F + λ|Θ| 1 } ,
where • F denotes the Frobenius norm defined for a matrix A by A 2 F = Trace(A A), |M | 1 = vec(M ) 1 denotes the 1 -norm of the vector formed by stacking the columns of M . It is thus closely related to the generalized thresholding estimator defined in [START_REF] Wen | Positive definite estimation of large covariance matrix using generalized nonconvex penalties[END_REF] and to the one defined in [START_REF] Rothman | Positive definite estimators of large covariance matrices[END_REF] with τ = 0 except that in our case |Θ -| 1 is replaced by |Θ| 1 where Θ -corresponds to the matrix Θ in which the diagonal terms are replaced by 0. The diagonal terms of Σ were indeed already removed in Γ r . Hence, we get Γ r by elementwise soft-thresholding that is by putting to zero the value of Γ r that are under a given threshold and by multiplying the non null values by a coefficient containing this threshold.

Here, we choose to estimate Γ r by Γ r (λ) defined through B(λ) in ( 8) which corresponds to a hard-thresholding and we set the upper triangular part of the estimator Σ(λ) of Σ to be equal to Γ r (λ). Since the diagonal terms of Σ are assumed to be equal to 1, we take the diagonal terms of Σ(λ) equal to 1. The lower triangular part of Σ(λ) is then obtained by symmetry.

The choice of the best parameter λ denoted λ final in the following will be discussed in Section 3.2.

2.3.

Positive definiteness. To ensure the positive definiteness of our estimator Σ of Σ, we consider the nearest correlation matrix to Σ(λ final ) which is computed by using the methodology proposed by [START_REF] Higham | Computing the nearest correlation matrix -a problem from finance[END_REF] and which is implemented in the function nearPD of the R package Matrix, see [START_REF] Bates | Matrix: Sparse and Dense Matrix Classes and Methods[END_REF].

2.4. Estimation of Σ -1/2 . Even if providing an estimator of a large covariance matrix can be very useful in practice, it may also be interesting to efficiently estimate Σ -1/2 . Such an estimator can indeed be used in the general linear model in order to remove the dependence that may exist between the columns of the observations matrix. For further details on this point, we refer the reader to Perrot-Dockès et al. (2018), Perrot-Dockès et al. (2018) and to the R package MultiVarSel in which such an approach is proposed and implemented for performing variable selection in the multivariate linear model in the presence of dependence between the columns of the observation matrix.

Since Σ is a symmetric matrix, it can be rewritten as U DU , where D is a diagonal matrix and U is an orthogonal matrix. The matrix Σ -1/2 can thus be estimated by U D -1/2 U where D -1/2 is a diagonal matrix having its diagonal terms equal to the square root of the inverse of the singular values of Σ. However, inverting the square root of too small eigenvalues may lead to poor estimators of Σ -1/2 . This is the reason why we propose to estimate Σ -1/2 by ( 9)

Σ -1/2 t = U D -1/2 t U ,
where D

-1/2 t is a diagonal matrix such that its diagonal entries are equal to the square root of the inverse of the diagonal entries of D except for those which are smaller than a given threshold t which are replaced by 0 in D -1/2 t . The choice of t will be further discussed in Section 3.7.

2.5. Choice of the parameters. Our methodology for estimating Σ depends on two parameters: The number r of singular values kept for defining Γ r and the parameter λ which controls the sparsity level namely the number of zero values in B(λ) defined in (8).

For choosing r, we shall compare two strategies in Section 3.1:

• The Cattell criterion based on the Cattell's scree plot described in [START_REF] Cattell | The scree test for the number of factors[END_REF] and

• the PA permutation method proposed by [START_REF] Horn | A rationale and test for the number of factors in factor analysis[END_REF] and recently studied from a theoretical point of view byDobriban (2018).

To choose the parameter λ in (8), we shall compare two strategies in Section 3.2:

• The BL approach proposed in [START_REF] Bickel | Covariance regularization by thresholding[END_REF] based on cross-validation and • the Elbow method which consists in computing for different values of λ the Frobenius norm R -Σ(λ) F , where R and Σ(λ) are defined in (2) and at the end of Section 2.2, respectively. Then, it fits two simple linear regressions and chooses the value of λ achieving the best fit.

Numerical experiments

Our methodology described in the previous section is implemented in the R package BlockCov and is available from the CRAN (Comprehensive R Achive Network) and from GitHub.

We propose hereafter to investigate the performance of our approach for different types of matrices Σ defined in (1) and for different values of n and q. The four following cases considered correspond to different types of matrices Z, the matrices D being chosen accordingly to ensure that the matrix Σ has its diagonal terms equal to 1.

• Diagonal-Equal case. In this situation, Z has the structure displayed in the left part of Figure 1, namely it has 5 columns such that the numbers of the non values in the five columns are equal to 0.1 × q, 0.2 × q, 0.3 × q, 0.2 × q and 0.2 × q, respectively and the non null values are equal to • Extra-Diagonal-Unequal case. Z has the same structure as in the Extra-Diagonal-Equal case except that the values are randomly chosen as in the Diagonal-Unequal case except for the fourth column where the additional non values are still equal to -0.5 in the range [0.35 × q, 0.45 × q]. For n ∈ {10, 30, 50} and q ∈ {100, 500}, 100 n × q matrices E were generated such that its rows E 1 , E 2 , • • • , E n are i.i.d. q-dimensional zero-mean Gaussian vectors having a covariance matrix Σ chosen according to the four previous cases: Diagonal-Equal, Diagonal-Unequal, Extra-Diagonal-Equal or Extra-Diagonal-Unequal.

√ 0.7, √ 0.75, √ 0.65, √ 0
3.1. Low rank approximation. The approaches for choosing r described in Section 2.5 are illustrated in Figure 3 in the Extra-Diagonal-Unequal case. We can see from this figure that both methodologies find the right value of r which is here equal to 5.

To go further, we investigate the behavior of our methodologies from 100 replications of the matrix E for the four different types of Σ. Figure 4 displays the barplots associated to the estimation of r made in the different replications by the two approaches for the different scenarii. We can see from this figure that the PA criterion seems to be slightly more stable than the Cattell criterion when n ≥ 30. However, in the case where n = 10, the PA criterion underestimates the value of r. Moreover, in terms of computational time, the performance of Cattell is much better, see Figure 5.

Positions of the non null values.

For the four scenarios, the performance of the two approaches: BL and Elbow described in Section 2.5 for choosing λ and hence the number of non null values in Σ(λ) is illustrated in Figure 6. This figure displays the True Positive Rate (TPR) and the False Positive Rate (FPR) of the methodologies from 100 replications of the matrix E for the four different types of Σ and for different values of n and q. We can see from this figure that the performance of Elbow is on a par with the one of BL except for the case where n = 10 for which the performance of Elbow is slightly better in terms of True Positive Rate. Moreover, in terms of computational time, the performance of Elbow is much better, see Figure 7.

3.3.

Comparison with other methodologies. The goal of this section is to compare the statistical performance of our approach with other methodologies.

Since our goal is to estimate a covariance matrix containing blocks, we shall compare our approach with clustering techniques. Once the groups or blocks have be obtained, Σ is estimated by assuming that the corresponding matrix estimator is block-wise constant except for the diagonal blocks for which the diagonal entries are equal to 1 and the extra-diagonal terms are assumed to be equal. This gives a great advantage to these methodologies in the Diagonal-Equal and in the Extra-Diagonal-Equal scenarii. More precisely, let ρ i,j denote the value of the entries in the block having its rows corresponding to Group (or Cluster) i and its columns to Group (or Cluster) j. Then, for a given clustering C:

(10) ρ i,j =            1 #C(i)#C(j) k∈C(i), ∈C(j) R k, , if C(i) = C(j) 1 #C(i)(#C(i)-1) k∈C(i), ∈C(i),k = R k, , if C(i) = C(j)
,

where C(i) denotes the cluster i, #C(i) denotes the number of elements in the cluster C(i) and R k, is the (k, ) entry of the matrix R defined in Equation ( 2).

For the matrices Σ corresponding to the four scenarios previously described, we shall compare the statistical performance of the following methods:

• empirical which estimates Σ by R defined in (2),

• blocks which estimates Σ using the methodology described in this article with the criteria PA and BL for choosing r and λ, respectively, • blocks fast which estimates Σ using the methodology described in this article with the criteria Cattell and Elbow for choosing r and λ, respectively, • blocks real which estimates Σ using the methodology described in this article when r and the number of non null values are assumed to be known which gives access to the best value of λ, • hclust which estimates Σ by determining clusters using a hierarchical clustering with the "complete" agglomeration method described in [START_REF] Hastie | The Elements of Statistical Learning[END_REF] and then uses Equation ( 10) to estimate Σ, • Specc which estimates Σ by determining clusters using spectral clustering described in von Luxburg ( 2007) and estimates Σ with Equation ( 10), • kmeans which estimates Σ by determining clusters from a k-means clustering approach described in [START_REF] Hastie | The Elements of Statistical Learning[END_REF] and then uses Equation (10) to estimate Σ.

In order to improve the performance of the clustering approaches: hclust, Specc and kmeans, the real number of clusters has been provided to these methods. The performance of the different approaches is assessed using the Frobenius norm of the difference between Σ and its estimator.

Figure 8 displays the mean and standard deviations of the Frobenius norm of the difference between Σ and its estimator for different values of n and q in the four different cases: Diagonal-Equal, Diagonal-Unequal, Extra-Diagonal-Equal and Extra-Diagonal-Unequal. We can see from this figure that in the case where n = 10, the performance of blocks fast is on a par with the one of blocks real and is better than the one of blocks. In the case where n = 50, the performance of blocks is slightly better than the one of blocks fast and is similar to the one of blocks real. Moreover, in all cases, either blocks fast or blocks outperforms the other approaches.

Then, the estimators of Σ derived from blocks, blocks fast and blocks real were compared to the PDSCE estimator proposed by [START_REF] Rothman | Positive definite estimators of large covariance matrices[END_REF] and implemented in the R package PDSCE and to the estimator proposed by [START_REF] Blum | Sparse factor model for co-expression networks with an application using prior biological knowledge[END_REF] and implemented in the FANet package Blum et al. (2016a). Since the computational burden of PDSCE is high for large values of q, we limit ourselves to the Extra-Diagonal-Equal case when n = 30 and q = 100 for the comparison. Figure 9 displays the results. We can see from this figure that blocks, blocks fast and blocks real provide better results than PDSCE and FANet. However, it has to be noticed that PDSCE is not designed for dealing with block structured covariance matrices but just for providing sparse estimators of large covariance matrices. 3.4. Columns permutation. In practice, it may occur that the columns of E consisting of the rows E 1 , E 2 , . . . , E n are not ordered in a way which makes blocks appear in the matrix Σ. To address this issue, we propose to perform a hierarchical clustering on E beforehand and use the obtained permutation of the observations which guarantees that a cluster plot using this ordering will not have crossings of the branches. Let us denote E ord the matrix E in which the columns have been permuted according to this ordering and Σ ord the covariance matrix of each row of E ord . Then, we apply our methodology to E ord which should provide an efficient estimator of Σ ord . In order to get an estimator of Σ the columns and rows are permuted according to the ordering coming from the hierarchical clustering.

To assess the corresponding loss of performance, we generated for each matrix E used for making Figure 8 a matrix E perm in which the columns of E were randomly permuted. The associated covariance matrix is denoted Σ perm . Then, we applied the methodology described in the previous paragraph denoted blocks samp and blocks fast samp in Figure 10 thus providing Σ perm . The performance of this new methodology was compared to the methodology that we proposed in the previous sections (denoted blocks and blocks fast in Figure 10) when the columns of E were not permuted. The results are displayed in Figure 10. We can see from this figure that the performance of our approach does not seem to be altered by the permutation of the columns.

3.5. Numerical performance. Figure 11 displays the computational times for estimating Σ with the methods blocks and blocks fast for different values of q ranging from 100 to 3000 and n = 30. The timings were obtained on a workstation with 16 GB of RAM and Intel Core i7 (3.66GHz) CPU. Our methodology is implemented in the R package BlockCov which uses the R language (R Core Team, 2017) and relies on the R package Matrix. We can see from this figure that it takes around 3 minutes to estimate a 1000 × 1000 correlation matrix.

3.6. Choice of the threshold t for estimating Σ -1/2 . Since we are interested in assessing the ability of Σ -1/2 t defined in (9) to remove the dependence that may exist between the columns of E, we shall consider the Frobenius norm of Σ -1/2 t Σ Σ -1/2 t -Id q which should be close to zero, where Id q denotes the identity matrix of R q . Figure 12 displays the Frobenius

norm of Σ -1/2 t Σ Σ -1/2 t
-Id q for different threshold t. A threshold of 0.1 seems to provide a small error in terms of Frobenius norm. Hence, in the following, t will be equal to 0.1 and

Σ -1/2 0.1 will be referred as Σ -1/2
. This technique was applied to all of the estimators of Σ discussed in Section 3.3 to get different estimators of Σ -1/2 . The Frobenius norm of the error

Σ -1/2 Σ Σ -1/2
-Id q is used to compare the different estimators obtained by considering the different estimators of Σ. The results are displayed in Figure 13. We observe from this figure that in the case where n = 10 the estimators of Σ -1/2 derived from the empirical, the blocks fast and the blocks real estimators of Σ perform similarly and seem to be more adapted than the others to remove the dependence among the columns of E. However, when n = 50, the behavior is completely different. Firstly, in the Diagonal-Equal case, the estimator of Σ -1/2 derived from the hclust estimator of Σ seems to perform better than the others. Secondly, in the Diagonal-Unequal case, the estimator derived from blocks, blocks fast and blocks real perform similarly than the one obtained from hclust. Thirdly, in the Extra-Diagonal case, the estimators derived from blocks, blocks fast and blocks real methodology perform better than the other estimators.

Then, the estimators of Σ -1/2 derived from blocks, blocks fast and blocks real were compared to the GRAB estimator proposed by [START_REF] Hosseini | Learning sparse gaussian graphical models with overlapping blocks[END_REF]. Since the computational burden of GRAB is high for large values of q, we limit ourselves to the Extra-Diagonal-Equal case when n = 30 and q = 100 for the comparison. Figure 14 displays the results. We can see that blocks and blocks real provide better results than GRAB. However, it has to be noticed that the latter approach depends on a lot of parameters that were difficult to choose, thus we used the default ones. 3.7. Use of Σ -1/2 to remove the dependence in multivariate linear models. Eventually, we assess the performance of the BlockCov methodology to remove the dependence in the columns of an observation matrix in order to be used for variable selection in the multivariate linear model as it is performed in the MultiVarSel R package:

(11) Y = XB + E,
where Y is a n × q response matrix, X is a n × p design matrix, B is a coefficients matrix and E is an error matrix. Here, E 1 , E 2 , • • • , E n are n zero-mean i.i.d. q-dimensional Gaussian

-1/2 t Σ Σ -1/2 t -Id q , where Σ -1/2 t is computed for different thresholds t.
random vectors having a covariance matrix Σ. To achieve this goal, we generate observations Y according to this multivariate linear model. We choose q = 100, p = 3, n = 30 and bX is the design matrix of a one-way ANOVA model.We compared our methodology with the one proposed by [START_REF] Perthame | Stability of feature selection in classification issues for high-dimensional correlated data[END_REF] and implemented in the FADA R package [START_REF] Perthame | FADA: Variable Selection for Supervised Classification in High Dimension[END_REF]. We shall investigate the effect of the sparsity of B and of the signal to noise ratio (SNR) for the four scenarii defining Σ on the selection of the non null values of B in (11). Different signal to noise ratios are obtained by multiplying B in (11) by a coefficient κ.

Since the results are barely influenced by the scenario chosen for Σ, only the Extradiagonal-Equal case is displayed in Figure 15, the other scenarii are available in Annexe 6.1. We can -Id q in the Extra-Diagonal-Equal case. see from this figure that when the signal to noise ratio is low and the value of s is high, meaning that there is a lot of non-zero values, the FADA methodology performs better than the BlockCov methodology. Nevertheless, in the three other cases the performance of BlockCov is either better or on a par with the one of FADA methodology. obtained from 100 replications comparing the variables selected by the Mul-tiVarSel strategy using either Σ -1/2 obtained by BlockCov to remove the dependence or the methodology proposed by FADA methodology. κ is linked to the signal to noise ratio and s denotes the sparsity levels i.e the fraction of non-zero elements in B.

Application to "multi-omic" approaches to study seed quality

Climate change could lead to major crop failures in world. In the present study, we addressed the impact of mother plant environment on seed composition. Indeed, seed quality is of paramount ecological and agronomical importance. They are the most efficient form of dispersal of flowering plants in the environment. Seeds are remarkably adapted to harsh environmental conditions as long as they are in a quiescent state. Dry mature seeds (so called "orthodox seeds") are an appropriate resource for preservation of plant genetic diversity in seedbanks. It has been reported that the temperature regime during seed production affects agronomical traits such as seed germination potential, see [START_REF] Huang | The effect of temperature on reproduction in the summer and winter annual arabidopsis thaliana ecotypes bur and cvi[END_REF][START_REF] Macgregor | Seed production temperature regulation of primary dormancy occurs through control of seed coat phenylpropanoid metabolism[END_REF] and [START_REF] Kerdaffrec | The maternal environment interacts with genetic variation in regulating seed dormancy in swedish arabidopsis thaliana[END_REF]. In order to highlight biomarkers of seed quality according to thermal environment of the mother plant, Arabidopsis seeds were produces under three temperature regimes (14-16 o C, 18-22 o C or 25-28 o C under a long-day photoperiod). Dry mature seeds were analysed by shotgun proteomic and GC/MS-based metabolomics [START_REF] Durand | Combined proteomic and metabolomic profiling of the arabidopsis thaliana vps29 mutant reveals pleiotropic functions of the retromer in seed development[END_REF]. The choice to use the model plant, Arabidopsis, was motivated by the colossal effort of the international scientific community for its genome annotation. This plant remains at the forefront of modern genetics, genomics, plant modelling and system biology, see [START_REF] Provart | 50 years of arabidopsis research: highlights and future directions[END_REF]. Arabidopsis provides a very useful basis to study gene regulatory networks, and develop modelling and systems biology approaches for translational research towards agricultural applications.

In this section, we apply our R packages BlockCov and MultiVarSel Perrot-Dockès et al. ( 2019) to metabolomic and proteomic data to better understand the impact of the temperature on the seed quality. More precisely, we use the following modeling for our observations:

(12) Y = XB + E,
where Y is a n × q matrix containing the responses of the q metabolites (resp. the q proteins) for the n samples with n = 9, q = 199 (resp. q = 724) for the metabolomic (resp. proteomic) dataset, X is a n × 3 design matrix of a one-way ANOVA model, such that its first (resp. second, resp. third) column is a vector which is equal to 1 if the corresponding sample grows under low (resp medium, resp. elevated) temperatures and 0 otherwise. B is a coefficient matrix and E is such that its n rows E 1 , E 2 , • • • , E n are n zero-mean i.i.d. q-dimensional random vectors having a covariance matrix Σ. We used our R package BlockCov to estimate Σ and Σ -1/2 assuming that there exists a latent block structure in the covariance matrix of the rows of E. More precisely, we assume that there exists some groups of metabolites (resp. proteins) having the same behavior since they belong to the same biological process. Then, we plugged this estimator into our R package MultiVarSel to obtain a sparse estimation of B. Thanks to this estimator of B, we could identify the metabolites (resp. proteins) having a higher (resp. lower) concentration when the temperature is high or low.

4.1.

Results obtained for the metabolomic data. We first estimated the matrices Σ and Σ -1/2 associated to E defined in Equation ( 12) by using the methodology developed in this paper, namely the BlockCov package. By the results of Section 3, we know that the PA and BL approaches performed poorly when n = 10. Since here n = 9, we used the The estimation of Σ obtained with our methodology is displayed in Figure 17 once the rows and the columns have been permuted according to the ordering provided by the hierarchical clustering to make visible the latent block structure.

Using the estimator of Σ -1/2 provided by the BlockCov package in the R package MultiVarSel provides the sparse estimator of the matrix B defined in Model 12 and displayed in Figure 18. We can see from this figure that for the metabolite X5MTP the coefficient of the matrix B is positive when the temperature is high which means that the production of the metabolite X5MTP is larger in high temperature conditions than in low temperature conditions.

In order to go further in the biological interpretation, we wanted to better understand the underlying block structure of the estimator of the correlation matrix of the residuals based on metabolite abundances Σ. Thus, we applied a hierarchical clustering with 8 groups to this matrix in order to split it into blocks. The corresponding dendogram is on the left part of Figure 17. The matrix containing the correlation means within and between the blocks or groups of metabolites is displayed in Figure 19. The composition of the metabolites groups is available in Appendix 6.2.

Interestingly, we could observe that X5MTP belongs to Group 6 which displays an high correlations mean equal to 0.8 between the 14 metabolites that make it up. At least 6 metabolites of this group belong to the same family, namely glucosinolates (i.e. X4MTB, 4methylthiobutyl glucosinolate; X5MTP, 5-methylthiopentyl glucosinolate; X6MTH, 6-methylthiohexyl glucosinolate; X7MTH, 7-methylthiohexyl glucosinolate; X8MTO, 8-methylthiooctyl glucosinolate; UGlucosinolate140.1, unidentified glucosinolate). Glucosinolates (GLS) are specialized metabolites found in Brassicaceae and related families (e.g. Capparaceae), containing a β-thioglucose moiety, a sulfonated oxime moiety, and a variable aglycone side chain derived from a α-amino acid. These compounds contribute to the plant's overall defense mechanism, see [START_REF] Wittstock | Glucosinolate research in the arabidopsis era[END_REF]. Methylthio-GLS are derivated from methionine. Methionine is elongated through condensation with acetyl CoA and then, are converted to aldoximes through the action of individual members of the cytochrome P450 enzymes belonging to CYP79 family, see [START_REF] Field | Glucosinolate and amino acid biosynthesis in arabidopsis[END_REF]. The aldoxime undergoes condensation with a sulfur donor, and stepwise converted to GLS, followed by the side chain modification. The present results suggest that the accumulation of methionine-derived glucosinolate family is strongly coordinated in Arabidopsis seed. Moreover, we can see that they are influenced by the effect of the mother plant thermal environment. The matrix containing the correlation means within and between the blocks or groups of proteins is displayed in Figure 21. We can see from this figure that Group 8 has the highest correlation mean equal to 0.47. It consists of 34 proteins which are given in Appendix 6.3.

A basic gene ontology analysis (http://geneontology.org/) showed that proteins involved in response to stress (biotic and abiotic), in nitrogen and phosphorus metabolic processes, in photosynthesis and carbohydrate metabolic process and in oxidation-reduction process are overrepresented in this group, see Figure 22. Thus, the correlation estimated within Group 8 seems to reflect a functional coherence of the proteins of this group.

The variable selection in the multivariate linear model using the R package MultiVarSel provided 31 proteins differentially accumulated in seeds produced under low, medium or elevated temperature. An aspartyl protease (AT3G54400), belongs to both, the Group 8 and to the proteins selected by MultiVarSel. This cell wall associated protein was up-acccumulated in dry seeds produced under low temperature. The gene encoding for this protease was described as a cold responsive gene assigned to the C-repeat binding factor (CBF) regulatory pathway, see [START_REF] Vogel | The cbf cold response pathways of arabidopsis and tomato[END_REF]. This pathway is requested for regulation of dormancy induced by low temperatures, see [START_REF] Kendall | Induction of dormancy in arabidopsis summer annuals requires parallel regulation of dog1 and hormone metabolism by low temperature and cbf transcription factors[END_REF]. Consistently, in Figure 23, two other proteins related to cell wall organization, a beta-glucosidase (BGLC1, AT5G20950) and a translation elongation factor (eEF-1Bβ1, AT1G30230) were differentially accumulated in seeds produced under contrasted temperature. eEF-1Bβ1 is associated to plant development and is involved in cell wall formation, see [START_REF] Hossain | The translation elongation factor eef-1bβ1 is involved in cell wall biosynthesis and plant development in arabidopsis thaliana[END_REF]. These results suggest that cell wall rearrangements occur under temperature effect during seed maturation. As displayed in Figure 23, 6 other proteins involved in mRNA translation: AT1G02780, AT1G04170, AT1G18070, AT1G72370, AT2G04390 and AT3G04840 were selected. The absolute failure of seed germination in the presence of protein synthesis inhibitors underlines the essential role of translation for achieving this developmental process, see [START_REF] Rajjou | The effect of α-amanitin on the arabidopsis seed proteome highlights the distinct roles of stored and neosynthesized mrnas during germination[END_REF]. Previous studies highlighted the importance of selective and sequential mRNA translation during seed germination and seed dormancy, see [START_REF] Galland | Dynamic proteomics emphasizes the importance of selective mrna translation and protein turnover during arabidopsis seed germination[END_REF], [START_REF] Bai | Extensive translational regulation during seed germination revealed by polysomal profiling[END_REF] and [START_REF] Bai | Combined transcriptome and translatome analyses reveal a role for tryptophan-dependent auxin biosynthesis in the control of DOG1-dependent seed dormancy[END_REF]. Thus, exploring translational regulation during seed maturation and germination through the dynamic of mRNA recruitment on polysomes or either neosynthesized proteome are emerging fields in seed research.

Conclusion

In this paper, we propose a fully data-driven methodology for estimating large block structured sparse covariance matrices in the case where the number of variables is much larger than the number of samples without limiting ourselves to block diagonal matrices. Our methodology can also deal with matrices for which the block structure only appears if the columns and rows are permuted according to an unknown permutation. Our technique is implemented in the R package BlockCov which is available from the Comprehensive R Archive Network and from GitHub. In the course of this study, we have shown that BlockCov is a very efficient approach both from the statistical and numerical point of view. Moreover, its very low computational load makes its use possible even for very large covariance matrices having several thousands of rows and columns. and the Plant Observatory-Chemistry/Metabolism platform (IJPB, Versailles; Clement G) for the analysis of GC/MS-based metabolome analyses. 6.2. Groups of metabolites.

6.3. Groups of proteins.

Figure 1 .

 1 Figure 1. Examples of matrices Σ generated from different matrices Z leading to a block diagonal or to a more general block structure (extra-diagonal blocks).

Figure 2 .

 2 Figure 2. Examples of matrices Σ of Figure 1 in which the columns and rows are randomly permuted.

  respectively. • Diagonal-Unequal case. In this scenario, Z has the same structure as for the Diagonal-Equal case except that the non null values in the five columns are not fixed but randomly chosen in [ Diagonal-Equal case. Here, Z has the structure displayed in the right part of Figure 1. The values of the columns of Z are the same as those of the Diagonal-Equal case except for the fourth column which is assumed to contain additional non values equal to -0.5 in the range [0.35 × q, 0.45 × q].

Figure 3 .

 3 Figure3. Illustration of PA and Cattell criteria for choosing r when q = 500 and n = 30 in the Extra-Diagonal-Unequal case. The value of r found by both methodologies is displayed with a dotted line, the straight lines obtained for the Cattell criterion and the eigenvalues of the permuted matrices in the PA methodology are displayed in grey.

Figure 4 .Figure 5 .

 45 Figure 4. Barplots corresponding to the number of times where each value of r is chosen in the low-rank approximation from 100 replications for the two methodologies in the different scenarii for the different values of n et q.

Figure 6 .Figure 7 .

 67 Figure 6. Boxplots comparing the TPR (True Positive Rate) and the FPR (False positive Rate) of the two methodologies proposed to select the parameter λ from 100 replications in the different scenarii.

Figure 8 .

 8 Figure 8. Comparison of the Frobenius norm of Σ-Σ for different estimators Σ of Σ and for different Σ.

Figure 9 .

 9 Figure 9. Comparison of the Frobenius norm of Σ -Σ in the Extra-Diagonal-Equal case for n = 30 and q = 100.

Figure 10 .Figure 11 .

 1011 Figure 10. Comparison of the Frobenius norm of Σ -Σ, and Σ perm -Σ perm .

Figure 12 .

 12 Figure 12. Frobenius norm of Σ

Figure 13 .-Figure 14 .

 1314 Figure 13. Comparison of the Frobenius norm of the error Σ -1/2 Σ Σ -1/2 -Id q , for different estimators Σ of Σ.

Figure 15 .

 15 Figure15. Means of the ROC curves (left) and Precision Recall curves (right) obtained from 100 replications comparing the variables selected by the Mul-tiVarSel strategy using either Σ -1/2 obtained by BlockCov to remove the dependence or the methodology proposed by FADA methodology. κ is linked to the signal to noise ratio and s denotes the sparsity levels i.e the fraction of non-zero elements in B.

Figure 16 .

 16 Figure 16. Illustration of the Cattell and Elbow criteria.

Figure 17 .

 17 Figure 17. Estimator of the correlation matrix Σ of the rows of E once the rows and the columns have been permuted according to the ordering provided by the hierarchical clustering.

Figure 18 .Figure 19 .

 1819 Figure 18. Sparse estimator of the coefficients matrix B obtained thanks to the package MultiVarSel with a threshold of 0.95.

4. 2 .

 2 Results obtained for the proteomic data. The same study was conducted on the proteomic data. The estimator of the correlation matrix of the residuals based on proteine abundances Σ obtained with our methodology is displayed in Figure20once the rows and the columns have been permuted according to the ordering provided by the hierarchical clustering to make visible the latent block structure. To better understand the underlying block structure of Σ, we applied a hierarchical clustering with 9 groups to this matrix in order to split it into blocks. The corresponding dendogram is on the left part of Figure20.

Figure 20 .

 20 Figure 20. Estimator of the correlation matrix of the residuals of the protein accumulation measures once the rows and the columns of the residual matrix have been permuted according to the ordering provided by the hierarchical clustering.

Figure 21 .

 21 Figure 21. Means of the correlations between the groups of proteins.

Figure 23 .

 23 Figure 23. Values of the coefficients obtained using the package MultiVarSel with a threshold of 0.95 on the proteomic dataset.

Figure 24 .Figure 25 .

 2425 Figure24. Means of the ROC curves obtained from 100 replications comparing the variables selected by the MultiVarSel strategy using either Σ -1/2 obtained by BlockCov to remove the dependence or the methodology proposed by FADA methodology. κ is linked to the signal to noise ratio and s denotes the sparsity levels i.e the fraction of non-zero elements in B.

Photosynthesis Carbohydrate metabolic process Oxidation-reduction process Phosphorus metabolic process Response to biotic stimulus Response to abiotic stimulus Nitrogen compound metabolic process

	Response to stress
	0	5	10	15
		Expected from Arabidopsis genome		Proteins, group 8(34)
	Figure 22. Gene ontology (GO) term enrichment analysis of the 34 pro-
	teins belonging to Group 8. Data from PANTHER overrepresentation test
	(http://www.geneontology.org); One uploaded id (i.e. AT5G50370) mapped
	to two genes. Thus, GO term enrichment was performed on 35 elements. Blue
	bars: observed proteins in Group 8; Orange bars: expected result from the
	reference Arabidopsis genome.

Acknowledgments

We thank the members of the EcoSeed European project (FP7 Environment, Grant/Award Number: 311840 EcoSeed, Coord. I. Kranner). IJPB was supported by the Saclay Plant Sciences LABEX (ANR-10-LABX-0040-SPS). We also thank the people who produced the biological material and the proteomic and metabolomic analysis. In particular, we would like to thank the Warwick University (UWAR, Finch-Savage WE and Awan S) for the production of seeds, the Plant Observatory-Biochemistry platform (IJPB, Versailles; Bailly M, Cueff G) for having prepared the samples for the proteomics and metabolomics, the PAPPSO Proteomic Plateform (GQE-Moulon; Balliau T, Zivy M) for mass spectrometry-based proteome analysis

Allantoin U1530.2.314 U2513.2.296 Arabitol Erythritol U2053.6.321.1 U2513.2.296.1 myo.Inositol Ethanolamine U2109.3.305 U2692.9.361 Glycerol.3.P Sorbitol U2197.2.494 U2798.377 myo.Inositol.1.P Threitol U2315.2.245 U2942.2.556 Phosphate Xylitol U3898.1.204 U3063.0.361 U2206.2.299 Ethylphosphate U3415.9.498 Fructose Glucose Glucopyranose..H2O. Mannose U1154.3.156 Raffinose U1393.172 Ribose U1541.8.263 U1127.5.140 U1647.2.403 U1172.9.281 U1705.2.319.pentitol. U1559.4.217 U1729.0.273 U1628.9.233 U1816.2.228 U1849.2.285 U1859.2.246 U1927.0.204 U2076.9.204 U1939.1.210 U2170.6.361 U1983.0.217 U2184.1.299 U1983.0.217.1 U2251.5.361 U2012.7.361 U2278.6.361 U2282.4.349 U2550.7.149 U2400.1.179 U2731.2.160 U2779.9.361 U2857.8.342 U2929.1.297 U3041.1.361 U3080.7.361 U3100.8.361 Group 5 Group 6 Group 7 Group 8 Quercitrin X4MTB BenzoylGlucosinolate.2Breakdown Maleate Dehydroascorbate X5MTP Hexanenitrile.6methylthio Pentonate.4 Fumarate X6MTH Sinapinate.trans U1408.4.298 Sinapinate.cis X7MTH Anhydroglucose U1617.8.146 Arabinose X8MTO U1125.1.140 U1767.3.243 Galactose UGlucosinolate140.1 U1290.198 U1904.9.204 U1127.4.169 U1129.9.184 U1371.5.151 U2828.8.361 U1718.0.157 U1270.1.240 U1549.7.130 U2839.3.312 U1931.5.202 U1897.2.327 U1568.5.313 U2882.5.297 U2261.0.218 U2473.361 U1592.8.217 U3008.3.457 U2412.1.157 U2529.8.361 U1700.6.288 U3168.2.290 U2588.9.535 U2756.4.271 U1759.4.331 U3218.5.297 U2688.5.333 U2924.3.361 U1852.0.217 U3910.6.597.Trigalactosylglycerol. U3213.1.400 U3279.7.361 U1872.1.204.methyl.hexopyranoside. U2443.7.217 U1380.5.184 U1958.217 U2053.6.321 U2087.6.321 U2150.9.279 U2271.6.249 U3188.1.361 U3701.368 U4132.5.575