Stéphane Canu

Catherine Matias D'avoir Accepté De Rapporter

La Thèse Ainsi Qu'à

Edouard Duchesnay

Mélina Gallopin

Vincent Segura D'

Merci À Ludivine

Odélia

Ayoub, Perrine, Kyllian, Claire Salim

Arnaud

Elisabetta Liudmila

Antoine Kamari

Dominique, Eugénie Maurice Valérie

Merci

Vincent R Guillem

Merci Notamment À Mariem

Fazia Merci À Elisabeth

Théophile, Camélia Svetlana Guillemette

Jean Paul

Edith Max D'avoir Contribué

and sometimes assume a known or unknown underlying structure. These priors enable summarizing the interactions among variables in a data set.

To estimate graphs at several levels of granularity and uncover a hierarchical clustering structure on variables, we developed a novel graphical inference approach. This approach relies on neigh-borhood selection and convex hierarchical clustering techniques.

The resulting pseudo-likelihood function is optimized via a continuation method using Nesterov smoothing.

We also analyzed omics data from natural populations of poplars, jointly studying dierent types of omics data to highlight the mechanisms of regulation between epigenetic and genetic data. To address the heterogeneity of data sources, we transformed variables to return to the Gaussian framework.

-Introduction

Inferring networks from biological data is important to gain a more complete understanding of the biological mechanisms underlying a particular phenomenon. Networks are a natural way to integrate and describe how biological variables interact, allowing the identication of complex interactions. With the emergence of highthroughput sequencing techniques, it is possible to generate large amounts of omics data related to the genome, transcriptome, genetic variations and metabolome.

The highly dimensional nature of these data is the main diculty from a statistical and interpretation point of view. The objective of the research is to propose and study a network inference method that takes into account a group structure or hierarchy between variables. Identifying groups of densely connected nodes in the network may correspond to biological variables with related functions and oers the possibility of constructing multiscale structures to synthesize the information retrieved by the groups and improve interpretability. This clustering task can be considered before or in conjunction with the network inference task. Probabilistic graphical models represent a well-suited model for inferring relationships between variables. The research will focus on a general subclass of graphical models where the conditional distribution at the nodes is Gaussian and examine methods for estimating clusters with convex criteria.

The research is motivated by the EPITREE (Evolutionary and functional impact of epigenetic variations in forest trees) project, which seeks to understand how epigenetics, in this case DNA methylation, gene expression and allelic variation, inuence mechanisms of adaptation and phenotypic plasticity in forest trees.

Epigenetics is the study of heritable changes that aect gene expression without altering the DNA, while phenotypic plasticity is the ability of an individual genotype to express dierent values of a given phenotypic trait under dierent environmental conditions [START_REF] Rey | Adaptation to global change: a transposable elementepigenetics perspective[END_REF]. Trees are remarkable organisms that are longlived, have complex life cycles, and produce wood, while providing a wide range of ecosystem services. In recent decades, widespread forest dieback due to drought and heat stress has been observed worldwide [START_REF] William Rl Anderegg | Meta-analysis reveals that hydraulic traits explain cross-species patterns of drought-induced tree mortality across the globe[END_REF]. These events highlight the vulnerability of forest ecosystems to environmental change and the urgent need to understand how trees respond to environmental stress. Indeed, climate change is the factor that will have the greatest impact on biodiversity by 2100, after land use [START_REF] Stuart | Consequences of changing biodiversity[END_REF]. It is known that forest trees have complex mechanisms that allow them to adapt to environmental stressors [START_REF] Toby | Stressful memories of plants: evidence and possible mechanisms[END_REF]. Understanding the molecular mechanisms underlying their adaptation is therefore essential for developing conservation and management strategies for forest ecosystems. The EPITREE project studies the molecular mechanisms underlying tree adaptation, focusing on two tree models, poplar and oak. These species were chosen for their genetic diversity and their potential to adapt to changing environments.

The project consists of several work packages, which include screening of candidate regions for epigenomic analysis and genome sequencing using modern omics technologies. These technologies have generated data on methylated single polymorphisms, dierentially methylated regions, gene expression and single nucleotide polymorphisms. The PhD research is primarily motivated by the fourth work package, which aims to perform an integrative analysis to model the multi-scale relationships between quantitative traits and their molecular determinants. Specically, this work package aims to quantify the contribution of genetic and epigenetic diversity to phenotypic variation, to study the impact of oak and poplar evolution on epigenomic plasticity, and to determine whether dierentially methylated regions at the gene level are conserved between the two species. In addition, this project aims to improve models for predicting quantitative trait variation by combining genetic and epigenetic information. Graphical models are a useful tool for inferring interactions between genetic information and methylation patterns, and may provide answers to some of the questions raised in the project.

Probabilistic graphical models (PGMs, Lauritzen (1996); [START_REF] Koller | Probabilistic Graphical Models: Principles[END_REF]) have become a popular tool for analyzing high-dimensional data and capturing the interactions between variables. They are widely used in various applications, such as genomics and image analysis, to reduce the number of parameters by selecting the most relevant interactions between variables. One class of PGMs that is particularly useful in Gaussian settings is the undirected Gaussian graphical models (GGMs). In high-dimensional statistics, Gaussian graphical models are often assumed to be sparse, meaning that only a small number of variables interact compared to the total number of possible interactions. This sparsity assumption oers both statistical and computational advantages by simplifying the dependence structure between variables [START_REF] Dempster | Covariance Selection[END_REF] and enabling the development of ecient algorithms. To support this approach, many researchers have developed methods to estimate sparse GGMs from data. These methods include neighborhood selection and penalized maximum likelihood estimation.

In undirected Gaussian graphical models, inferring the conditional independence graph (CIG) involves identifying the support of the precision matrix Ω (the inverse of the variance-covariance matrix). To learn the CIG of GGMs, several ℓ 1 -penalized methods have been proposed in the literature. One popular method is the neighborhood selection (MB, [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF]) approach based on nodewise regression using the least absolute shrinkage and selection operator (LASSO). This method focuses on learning only the structure of the network.

The MB method has spawned a long line of work in nodewise regression methods, including extensions with various forms of sparsity-inducing penalties such as the Dantzig selector [START_REF] Yuan | High dimensional inverse covariance matrix estimation via linear programming[END_REF] and the Clime estimator [START_REF] Cai | A constrained l1 minimization approach to sparse precision matrix estimation[END_REF].

Another family of sparse CIG inference methods directly estimates Ω via the minimization of the ℓ 1 -penalized negative log-likelihood [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF]. This method, called the graphical LASSO (GLASSO, [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]), benets from many optimization algorithms [START_REF] Yuan | Model selection and estimation in the Gaussian graphical model[END_REF]Rothman et al., 2008a;[START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF][START_REF] Hsieh | Quic: Quadratic approximation for sparse inverse covariance estimation[END_REF].

LASSO-type regularization methods are widely used for conditional independence graph estimation and have been shown to be robust to high-dimensional problems. However, they have limitations in the presence of strongly correlated variables, which are well-known and have been discussed in the literature [START_REF] Bühlmann | Correlated variables in regression: clustering and sparse estimation[END_REF]Park et al., 2006). To overcome these limitations and improve the estimation procedure, previous works have attempted to integrate clustering structures among the variables. Several studies have proposed dierent methods for incorporating clustering structures, including [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF], [START_REF] Ambroise | Inferring sparse gaussian graphical models with latent structure[END_REF], Mazumder and Hastie (2012a), Tan et al. (2015), [START_REF] Devijver | Block-Diagonal Covariance Selection for High-Dimensional Gaussian Graphical Models[END_REF], and [START_REF] Yao | Clustered gaussian graphical model via symmetric convex clustering[END_REF].

The methods discussed earlier utilize the group structure to simplify the graph inference problem and infer the conditional independence graph between single variables. However, the inference of the CIG between groups of variables or representative variables of the groups has received less attention. Although some works have addressed this problem, they have mostly focused on two-level estimations, i.e., at the level of single variables and provided known groups (see, e.g., [START_REF] Cheng | Multilevel Gaussian graphical model for multilevel networks[END_REF]). The research problem addressed in this work aims to dene an inference method that allows for more than two levels of granularity estimations with unknown groups. This problem is mainly motivated by applications in biological data analysis where data from multiple sources, typically multi-omic data, need to be analyzed. In such cases, it might be necessary to group variables sharing the same characteristics and simultaneously take that into account in the network inference procedure, using a unique cost function instead of alternating clustering and network inference tasks.

Our research aims to address the inference of hierarchical clustering structures that are more intuitive for interpretation and focus on learning the inferred network structure rather than estimating the coecients of the precision matrix.

Although the applications are mainly motivated by biological questions from the EPITREE project, answering these questions may require using other dedicated machine learning tools. From a mathematical background, our research is located at the crossroads of probabilistic graphical inference, clustering, and convex optimization. From a statistical biological background, we address various biological questions, primarily concerning the impact of epigenetics on poplar local adaptation, which requires tools such as dierential gene analysis, enrichment analysis of gene sets, transformation of count data, and gene selection methods. Omics applications include applications on transcriptomic (gene expression), epigenetic (especially DNA methylation), genetic (especially SNPs), and an illustration on metagenomic data (microbial abundance).

The proposed methodology for graphical inference is called Multiscale Graphi-cal LASSO (MGLASSO), which is a pseudo-likelihood-based method for estimating hierarchical clustering structures and graphical models that depict the conditional independence structure between clusters of variables at each level of the hierarchy.

MGLASSO combines neighborhood selection with a fused-LASSO type penalty for clustering [START_REF] Pelckmans | Convex clustering shrinkage[END_REF][START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF][START_REF] Lindsten | Clustering using sum-of-norms regularization: With application to particle lter output computation[END_REF].

While the use of fusion penalties in Gaussian graphical model inference has been widely studied, previous works have mainly focused on penalized likelihood and investigated fusion penalties for enforcing local constancy in the nodes of the inferred network [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF][START_REF] Yao | Clustered gaussian graphical model via symmetric convex clustering[END_REF][START_REF] Lin | Estimation of sparse gaussian graphical models with hidden clustering structure[END_REF]. In contrast, MGLASSO employs a pseudo-likelihood criterion that is more computationally ecient and establishes a link with multiscale graphical models. Although the criterion used is similar to that used in supervised convex clustering [START_REF] Hallac | Network lasso: Clustering and optimization in large graphs[END_REF][START_REF] Chu | Adaptive convex clustering of generalized linear models with application in purchase likelihood prediction[END_REF], MGLASSO behaves more like a multitask learning problem [START_REF] Chiquet | Inferring multiple graphical structures[END_REF] due to its coupling with Gaussian graphical inference. In biological applications, MGLASSO relies on various data transformations adapted to the nature of the data, including the center-log ratio [START_REF] Aitchison | The statistical analysis of compositional data[END_REF] for compositional data.

MGLASSO, like [START_REF] Yao | Clustered gaussian graphical model via symmetric convex clustering[END_REF], is a method that combines Gaussian graphical models and convex clustering. However, unlike their work, we focused on the neighborhood selection framework and proposed adding a sparsity-inducing penalty (LASSO) to produce sparser results. We have also made available a beta version R package on the CRAN that implements the approach [START_REF] Sanou | mglasso: Multiscale Graphical Lasso[END_REF].

The algorithm uses a basic path algorithm to highlight the estimated multiscale structures. Our approach can also be seen as an extension of the SpiecEasi (Kurtz et al., 2015) method to multiscale networks when applied to compositional data with the centered log ratio transformation. Our biological applications, within the framework of the EPITREE project, demonstrate that the DNA methylation epigenetic mark for poplars can be used as markers of the genetic structure of the studied populations [START_REF] Mamadou | Epigenetic variation in tree evolution: a case study in black poplar (populus nigra)[END_REF].

The remaining of the manuscript is structured into three chapters. Chapter 2 provides a foundational understanding of graphical modeling, convex clustering, and convex optimization, which are necessary for comprehending the rest of the manuscript. In Chapter 3, the Multiscale Graphical Lasso (MGLASSO) is introduced as a novel approach to Gaussian graphical model inference that combines Vectors and matrices X = (X 1 , . . . , X p ) p-dimensional random vector or data matrix if ∈ R n×p X A subset of X with variables indices taken in A V al (X) set of values that X can take X k k-th column in X X i i-th row in X X \(i,j) X deprived of columns i-th and j-th columns β i ∈ R p-1 regression vector β i k multiple regression coecient of i-th variable on k-th variable V ec (.) convert a matrix into a column vector. X ⊥ ⊥ Y|Z X is conditionally independent of Y given Z. This chapter which can be read independently establishes the necessary mathematical background for the rest of the manuscript and reviews existing research of the eld.

Gaussian distribution

Σ covariance matrix S empirical covariance matrix Ω = Σ -1 ,
The manuscript focuses on a subclass of graphical models with no directional arrows between the edges, as described in Section 2.1. The foundations of the Gaussian graphical model, a particular case of the undirected graph model, are introduced. The model and the central notion of conditional independence are dened. Section 2.1.3 presents an overview of state-of-the-art inference approaches for Gaussian graphical models. Section 2.2 presents convex clustering, a method used to group data points based on a convex criterion and will be used in conjunction with the Gaussian graphical model inference problem later in the manuscript. Some essential properties of the convex approach and the links between convex clustering and clustering approaches, such as hierarchical agglomerative clustering and the k-means method, are discussed.

In Section 2.3, non-smooth convex optimization techniques are reviewed, including subgradient methods, proximal methods such as the alternating direction multiplier method, and smoothing methods, with emphasis on the Nesterov smoothing technique.

Undirected Graphical Models

Probabilistic graphical models [START_REF] Steen | Graphical models[END_REF][START_REF] Koller | Probabilistic Graphical Models: Principles[END_REF], often seen as a marriage between probability and graph theory [START_REF] Barber | Bayesian reasoning and machine learning[END_REF], are widely used in high-dimensional data analysis to synthesize the interactions between variables. In many applications, such as statistical physics, genomics, image analysis, or social network analysis, graphical models can reduce the number of parameters by selecting the most relevant interactions through parsimony. A graph consists of a set of nodes and edges between the nodes. There are two prominent traditional families of graphical models depending on the nature of the edge that connects the nodes:

Markov random elds (MRF), also called undirected graphical models: These models were rst applied in 1902 by Gibbs (Bryan, 1902) to describe the behavior of a system of interacting particles. MRFs are based on undirected graphs with only undirected edges, hence useful to describe soft constraints between variables. The genotype corresponding to the blood type is the pair (X F , X M ) which takes values in {A, B, O} 2 . These types of models can be used for genetic counseling, for example, to predict the genotypes of expected children.

Genotype

C Genotype F Blood type F Genotype M Blood type F Blood type C Figure 2.2: Directed acyclic graph
Another family of graphs, less used, derived from the previous ones is the mixed graphs [START_REF] Sadeghi | Markov properties for mixed graphs[END_REF]. They contain more than one type of edge: directed, undirected and bidirected. Regardless of the type of graph, graphical models gravitate around the central concept of conditional independence. It occurs when two sets of variables are independent given an additional set.

The following sections focus on the undirected Gaussian graphical model and its inference problem. The three Markov properties are recalled as well as their links with conditional independence. Then, an overview of existing sparse inference approaches is given. Finally, inference methods that assume the existence of an underlying clustering structure on the variables are presented.

Markov properties and factorization

This section contains some theoretical results related to undirected graphical models.

Markov properties for Undirected Graphical Models

As mentioned, a graph is undirected if all edges have no directional arrows.

We usually denote undirected graphs as G. Let V = {1, . . . , p} and let P 2 (V ) be the subsets of V of size 2.

Denition 2.1.

For E ⊆ P 2 (V ), an undirected graph is a pair G = (V, E)
where V is the set of vertices and E the set of edges. The graph is complete if

E = P 2 (V ). For a subset C ⊂ V , G C = (C, E C ) is the subgraph of G induced by C, with E C the set of edges in C × C. Whenever G C is complete, C is said to be a clique of G.
For example, the graph in Figure 2.3 consists of four nodes forming a complete graph with six undirected edges. The subgraphs induced by the sets of nodes

{X 1 , X 2 , X 3 }, {X 2 , X 3 , X 4 } or {X 1 , X 2 , X 4 } are cliques of size 3 of the graph G. X 1 X 2 X 3 X 4 Figure 2.3: The set C 1 = {X 1 , X 2 , X 3 } is an example of clique. The set C 2 = {X 1 , X 2 , X 3 , X 4 } is a complete graph.
Undirected Graphical Models can be translated to conditional independence constraints through a set of rules called Markov properties.

Denition 2.2 (Conditional independence). Let X A , X B , X C be sets of random variables. X A is said to be independent of X B given X C in a joint probability distribution P and writes

X A ⊥ ⊥ X B |X C i. P (X A = x A , X B = x B |X C = x C ) = P (X A = x A |X C = x C )P (X B = x B |X C = x C ) for all values x = (x A , x B , x C ) ∈ V al(X)
, where V al(X) denote the set of values that X can take.

The concept of conditional independence provides a precise and formal meaning to the idea of information irrelevance. The notation X A ⊥ ⊥ X B |X C can be interpreted as follows: given the information contained in X C , knowledge of X B is irrelevant for understanding X A .

Let us consider conditional independence in the context where X = (X v ), where v ∈ V , is a p-dimensional random vector taking values in some space X p = ⊗ p s=1 X s , with a joint probability distribution P . Depending on the application, the space X p can be continuous or discrete. When the collection of random variables is associated with an undirected graph G, three Markov properties can be established.

Let us rst introduce the notion of separation, which is based on the idea that specic subsets of nodes in a graph can block the ow of information between other subsets. Denition 2.3 (Separation). Let A, B, S be subsets of V. The set S is said to separate set A from set B if any path from any element of A to any element of B passes through S.

Denition 2.4 (Global Markov property). Let A, B, S ⊂ V be three disjoint subsets such that S separates A from B in the graph G. X satises the global Markov property with respect to the graph G i.

X A ⊥ ⊥ X B |X S .
(2.1)

where

X A = (X k )\k ∈ A.
The global Markov property is closely related to the other two Markov properties: the local and pairwise Markov properties.

Denote ne(i) = {u ∈ V : {u, i} ∈ E} the neighbourhood or Markov blanket of node i and cl(i) = ne(i) ∪ {i}, the closure of node i. The local Markov property states that each node is conditionally independent of its non-neighbors, given its neighbors.

Denition 2.5 (Local Markov property). X satises the local Markov Property with respect to the graph G i.

X i ⊥ ⊥ X V \ cl(i) |X ne(i) (2.2)
for any node i ∈ V.

In contrast, the pairwise Markov property states that two non-adjacent nodes are conditionally independent given their common neighbors.

Denition 2.6 (Pairwise Markov property). X satises the pairwise Markov property with respect to the graph G i.

X i ⊥ ⊥ X j |X V \{i,j} (2.3)
whenever there is no edge between nodes i and j ie (i, j) / ∈ E.

The global Markov property (2.1) implies the local property (2.2), which in turn implies the pairwise property (2.3) [START_REF] Steen | Graphical models[END_REF]. When X has a continuous and strictly positive joint density with respect to the Lebesgue measure, the three properties are equivalent.

Thanks to the Markov properties, particularly the global Markov property, the graph in Figure 2.4 can be interpreted as follows: X 3 ⊥ ⊥ X 4 |X 2 . This is because variable X 2 separates the variables X 1 and X 4 .

X 1 X 2 X 3 X 4

Factorization

Another alternative to connect the graphical and probabilistic structure is through factorization, which is used as a basis for many inference algorithms.

Denition 2.7. The density function f of probability distribution P with respect to product measure ν is said to factorize with respect to the graph G if it can be represented as follows

f (x) = 1 Z C∈C ϕ C (x C ),
(2.4)

where C denotes the set of maximal cliques, ϕ = (ϕ C , c ∈ C) is a collection of positive potential functions and Z is a normalization constant.

The factorization properties of graphical models make it possible to perform tractable computations with multivariate distributions. By representing a joint probability distribution as a product of factors, where each factor depends only on a subset of variables corresponding to a clique in the graph, computations can be performed locally on the cliques and then combined using factorization. This enables the development of ecient inference and learning algorithms for large and complex models.

It is also worth noting that factorization implies the global Markov property for any probability distribution P [START_REF] Steen | Graphical models[END_REF]. In the special case of strictly positive distributions, there is an equivalence between factorization and the three Markov properties, through the following theorem:

Theorem 2.1 (Hammersley-Cliord). A strictly positive and continuous probability distribution P with respect to a product measure ν factorizes with respect to a graph G = (V, E) if and only if P satises the pairwise Markov property with respect to G.

In other words, a probability distribution thus satises the Markov properties with respect to an undirected graph if and only if it can be expressed as a positive product of potential functions, where each potential function depends only on a subset of variables corresponding to a clique in the graph. An undirected graphical model can then be considererd as a pair (G, P ) where G is a graph with undirected edges and P is a distribution that factorizes with respect to G.

The choice of the potential functions ϕ C determines the structure of the dependencies between variables in the graphical model.

Pairwise Markov networks

In a pairwise Markov network, the joint distribution over a set of random variables is represented as a product of potential functions, each of which involves at most two variables. Log-linear models are a way of parameterizing the potential functions in a pairwise Markov network. Instead of directly specifying the potential functions, they model the logarithm of these functions as linear combinations of features. Each feature represents a particular conguration of the variables, and its weight determines the contribution of that conguration to the overall potential function. The potential functions are dened as:

ϕ C (x C ) = exp (-ϵ C (x C )) , where ϵ C (x C ) = -log ϕ C (x C
) is called an energy function [START_REF] Koller | Probabilistic Graphical Models: Principles[END_REF]. Generally, ϵ C is chosen such as

ϵ C (x C ) = -w C f C (x C )
where w C is the clique weight and f C a feature function over the clique.

Denition 2.8. Pairwise Markov Networks can be dened as a subclass of undirected graphical models for which the factorizing density function can be written as

f (x) = 1 Z exp   p i=1 w i f i (x i ) + i<j w ij f ij (x i , x j )   (2.5)
where w i , w ij , ∀i ̸ = j are weights, f i , f ij feature functions and Z is a normalization constant.

This family of Markov networks is well-known for its simplicity in parameteri- Gaussian graphical models, which use Markov random elds, are a popular type of pairwise Markov network.

Gaussian graphical models

Gaussian graphical models (GGMs) or covariance selection models [START_REF] Steen | Graphical models[END_REF] are a special class of undirected graphical models used in gaussian settings.

Let X = (X 1 , . . . , X p ) T ∈ R p be a p-dimensional Gaussian random vector, with zero mean and covariance matrix Σ ∈ S p ≻0 , where S p ≻0 denote the set of real symmetric p × p positive denite matrices. Some properties of graphical models specic to Gaussian distributions are given below.

Proposition 2.1. The conditional independence structure of X ∼ N p (0, Σ) is characterized by the graph G which is uniquely determined by the support of the precision or concentration matrix Ω = Σ -1 .

Proof. The probability density function f of the multivariate normal distribution is dened as

f (X) = (det(Ω)) 1/2 (2π) p/2 exp - 1 2 X T ΩX (2.6)
The quadratic term can be rewritten as

exp - 1 2 X T ΩX = exp   - 1 2 (i,j)∈E Ω ij X i X j   = (i,j)∈E ϕ (i,j)∈E (X i , X j ) with ϕ (i,j)∈E (X i , X j ) = exp -1 2 Ω ij X i X j .
Hence, the distribution can be factorized in terms of potentials over cliques composed of at most two nodes. The Hammersley-Cliord theorem (2.1) thus ensures that, for any two vertices i, j / ∈ V , Ω ij = 0 if and only if the i-th and j-th variables are conditionally independent given the others i.e. X i ⊥ ⊥ X j |X \(i,j) .

Proposition 2.2. The entries of the precision matrix are proportional to partial correlation coecients.

Indeed, the partial correlation between X i and X j given X\{X i , X j } is equal to

-Ω ij √ Ω ii Ω jj
. The reader may refer to [START_REF] Steen | Graphical models[END_REF] for proof.

The following corollary can be derived from partial correlations coecients being directly related to regression coecients.

Proposition 2.3. Given the regression problem X i = X \i β i + ϵ i , where ϵ i is the normal residuals vector, the regression coecient is given by

β i k = -Ω ik /Ω ii .
Proposition 6.1 suggests that the GGMs can be estimated by a series of regressions as outlined by [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF]. The next section will introduce some GGMs estimation methods.

Inference of Gaussian Graphical Models

Let P be an unknown probability distribution that factorizes over a graph G. Given a set of independent and identically distributed (iid) samples from P , the task of learning a Gaussian graphical model is to estimate the potential functions that best t the distribution [START_REF] Maathuis | Handbook of graphical models[END_REF]. In other words, the goal is to infer the edges of the graph and the parameters of the distribution. However, we also include inference methods in this review that consistently recover the graph structure without necessarily providing consistent parameter estimation.

We dierentiate learning approaches into three main classes: those without sparsity constraints, those which include sparsity-based approaches, and those with additional constraints on the node structure.

Maximum Likelihood Estimation

A natural way to estimate a zero-mean GGM is by using the maximum likelihood estimator (MLE). The goal is to maximize the following strictly concave [START_REF] Hastie | Statistical learning with sparsity[END_REF] log-likelihood function given x a realization of X:

l(Ω) = n i=1 log f (x i |Ω) ∝ log det(Ω) -tr(ΩS)
(2.7) where S = 1 n X T X is the empirical covariance matrix and log det the logarithm of the determinant of the matrix. The solution to this problem is unique and given by ΩMLE = S -1

(2.8) whenever S is non singular [START_REF] Maathuis | Handbook of graphical models[END_REF]. However, MLE is not computable when n is lower than p. Even if calculable, it often performs poorly (Mazumder and Hastie, 2012a). Indeed the estimator can result in complete graphs. In the context of high-dimensional statistics, GGMs are generally assumed to be sparse, meaning that a small number of variables interact compared to the total number of possible interactions. [START_REF] Dempster | Covariance Selection[END_REF] introduced the idea of estimating the network structure by setting some elements of Ω to zero.

This assumption has been shown to simplify the structure of dependencies between variables.

Edge recovery with sparsity constraints

Following the idea of [START_REF] Dempster | Covariance Selection[END_REF], several authors proposed approaches to recover the precision matrix support using sparsity.

Lasso penalized pseudo-likelihood estimation [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] originally introduced a nodewise regression approach for neighbourhood selection based on the least absolute shrinkage and selection operator (Lasso, [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]). They regress each variable X i , i ∈ V on the predictors X \{i} := {X k |k ∈ V \{i}}, taking advantage of the link between regression coecients and precision matrix entries (see proposition 6.1).

Denote βi ∈ R p-1 the regression vector when X i is considered as response variable. ∀i ∈ [1, p], the Lasso-based neighborhood regression is solved by optimizing the following problem: βi := βi (λ) = argmin

β i ∈R p-1 1 n X i -X \i β i 2 2 + λ β i 1 .
(2.9)

where λ is a non negative regularization parameter and X \i denotes the matrix X deprived of column i. The neighbourhood is thus all the vertices corresponding to the non-zero regression coecients ie ne(i) = {u ∈ V \{i}| βi j ̸ = 0}. The dierent regressions are then combined applying an `AND' or `OR' rules to infer the conditional independence graph. `AND' add an edge between X i and X j if β i j ̸ = 0 and β j i ̸ = 0. While in `OR', an edge is added when β i j ̸ = 0 or β j i ̸ = 0.

Later, authors like [START_REF] Guilherme | A path following algorithm for sparse pseudo-likelihood inverse covariance estimation[END_REF], [START_REF] Ambroise | Inferring sparse gaussian graphical models with latent structure[END_REF] showed that the neighborhood selection could be seen as a pseudo-likelihood [START_REF] Besag | Statistical analysis of non-lattice data[END_REF] approximation of the global likelihood.

Proposition 2.4. The neighborhood selection with Lasso applied to each node of a gaussian graphical model is equivalent to maximizing a penalized pseudolikelihood estimation.

A proof is given in the appendix 7.1. Penalized pseudo-likelihood approximations are less time-consuming to optimize than penalized maximum likelihood methods that will be introduced later. They might not necessarily yield better parameter estimates, but in the context of GGMs edge recovery, they lead to consistent graph structure estimation under `AND' or `OR' rules [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF].

Other penalized pseudo-likelihood methods include the work of Rocha et al.

( 2008) who merged all the linear regressions in [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] into a single one and added symmetry constraints. A symmetrized version have also been proposed by [START_REF] Friedman | Applications of the lasso and grouped lasso to the estimation of sparse graphical models[END_REF] and [START_REF] Peng | Partial correlation estimation by joint sparse regression models[END_REF].

Lasso penalized global likelihood estimation

Several authors have subsequently taken the idea of Lasso penalization of [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] and applied it to the overall likelihood of the model.

In contrast to neighborhood selection approaches, which are mainly concerned with estimating the structure of the graph, likelihood-based approaches allow learning the graph structure and consistently estimating its parameters simultaneously.

They seek to optimize the following convex function [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF])

Ω = arg max Ω∈S p ≻0 log det(Ω) -tr(ΩS) -λ ∥Ω∥ 1 (2.10)
where S is the empirical covariance matrix, tr(A) the trace of a matrix A, λ a positive penalty parameter and ∥Ω∥ 1 = i,j |Ω| ij . Computation is one of the main challenges of problem (2.10). A variety of optimization algorithms have been proposed over the last decades. Some of them are highlighted here. [START_REF] Yuan | Model selection and estimation in the Gaussian graphical model[END_REF] solved the problem (2.10) by applying the interior points methods. They exploited the link between (2.10) and the determinant maximization problem [START_REF] Vandenberghe | Determinant maximization with linear matrix inequality constraints[END_REF]. However, interior point algorithms are not ecient for large-scale problems [START_REF] Hastie | Statistical learning with sparsity[END_REF]. [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF] introduced a dual block coordinate descent approach interpretable as recursive ℓ 1 -norm penalized regressions. The approach has been rened by [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] and called graphical Lasso. First order optimization methods include d' Aspremont et al. (2008); [START_REF] Lu | Adaptive rst-order methods for general sparse inverse covariance selection[END_REF]. Later other approaches emerged, which allowed improving the estimation process signicantly. Mazumder and Hastie (2012b) and [START_REF] Daniela M Witten | New insights and faster computations for the graphical lasso[END_REF] used the block diagonal screening rule to speed up the graphical Lasso algorithm. They showed that the sample covariance matrix's thresholding leads to the conditional independence graph connected components. Indeed, the inverse of a block diagonal matrix has the same blocks as the matrix.

The inference problem can then be decomposed into low-dimensional subproblems within each connected component. [START_REF] Hsieh | Big & quic: Sparse inverse covariance estimation for a million variables[END_REF] developed a proximal Newton approach based on the work of QUIC algorithm [START_REF] Hsieh | Quic: Quadratic approximation for sparse inverse covariance estimation[END_REF] that can scale up to problems with a million variables. It is a block coordinate descent algorithm combined with METIS graph clustering. Their algorithm converges faster than the methods mentioned above.

The Graphical Lasso properties have been largely studied in [START_REF] Yuan | Model selection and estimation in the Gaussian graphical model[END_REF]; [START_REF] Ravikumar | Highdimensional covariance estimation by minimizing l1-penalized log-determinant divergence[END_REF]; Rothman et al. (2008b). Hence the inference of GGMs with Lasso penalty has then been extended to dierent contexts to address specic applications/issues: for multiple sources data integration [START_REF] Chiquet | A multiattribute gaussian graphical model for inferring multiscale regulatory networks: an application in breast cancer[END_REF], gene regulation networks [START_REF] Charbonnier | Weighted-lasso for structured network inference from time course data[END_REF], clustered samples [START_REF] Chiquet | Inferring multiple graphical structures[END_REF][START_REF] Danaher | The joint graphical lasso for inverse covariance estimation across multiple classes[END_REF], multiple response variables (Chiquet et al., 2017b) or even accounting missing data (Robin et al., 2019).

Other sparsity penalties

Remark that in the literature of GGMs, sparsity penalties dierent from Lasso have also been explored. The ℓ 2 penalty has been studied by [START_REF] Mo Kuismin | Precision matrix estimation with rope[END_REF].

They called the approach ROPE, i.e., the Ridge type operator for precision matrix estimation. Using the ridge penalty yields a closed form solution to the inference problem contrarily to Graphical Lasso methods. However, the estimated precision matrix might not be sparse.

Nonconcave penalties as SCAD (Smoothly Clipped Absolute Deviation, [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]) and adaptive Lasso [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] have been used in place of Lasso by [START_REF] Fan | Network exploration via the adaptive lasso and scad penalties[END_REF] to address the ℓ 1 penalty bias that can occur because of the linear increase of the penalty on regression coecients. However, adaptive Lassolike penalties are sensitive to the choice of the initial estimate of the precision matrix, and the nonconvexity of SCAD can make the computation dicult.

Other penalties that belongs to the successive regressions family include the Dantzig selector [START_REF] Candes | The dantzig selector: Statistical estimation when p is much larger than n[END_REF], the Sqrt-Lasso [START_REF] Belloni | Square-root lasso: pivotal recovery of sparse signals via conic programming[END_REF], the scaled-Lasso [START_REF] Sun | Scaled sparse linear regression[END_REF], the grouped-Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] and the constrained ℓ 1 minimizaton introduced by Yuan (2010), [START_REF] Liu | Tiger: A tuning-insensitive approach for optimally estimating gaussian graphical models[END_REF], [START_REF] Sun | Sparse matrix inversion with scaled lasso[END_REF], [START_REF] Friedman | Applications of the lasso and grouped lasso to the estimation of sparse graphical models[END_REF] and [START_REF] Cai | A constrained l1 minimization approach to sparse precision matrix estimation[END_REF] respectively. [START_REF] Kovács | Graphical elastic net and target matrices: Fast algorithms and software for sparse precision matrix estimation[END_REF] proposed the graphical Elastic Net. The combination of ℓ 1 and ℓ 2 penalties in the inference problem leads to stable estimations when dealing with highly correlated variables. The approach includes in the estimation procedure the addition of a target matrix which is prior knowledge that is provided beforehand. Indeed, in strongly correlated variables presence, Lasso performance is known to be impaired [START_REF] Bühlmann | Correlated variables in regression: clustering and sparse estimation[END_REF]Park et al., 2006;[START_REF] Vigneau | Clustering of variables for enhanced interpretability of predictive models[END_REF][START_REF] Grimonprez | Mlgl: An r package implementing correlated variable selection by hierarchical clustering and group-lasso[END_REF]. In the following, we focus on approaches that tackle this Lasso bias.

Inference while taking into account underlying structure

In a penalized regression problem, the Lasso selects only one feature from a group of correlated features [START_REF] Bühlmann | Correlated variables in regression: clustering and sparse estimation[END_REF]. A variety of solutions have been proposed using dierent sparsity penalties. Among them, The Elastic-net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] applies a linear combination of Lasso and ridge penalties which encourages a grouping eect and can select groups of variables. OSCAR [START_REF] Howard | Simultaneous regression shrinkage, variable selection, and supervised clustering of predictors with oscar[END_REF] achieved that by mixing Lasso and ℓ ∞ penalizations.

For its part, the clustered Lasso [START_REF] She | Sparse regression with exact clustering[END_REF] dened a sort of generalized Fused Lasso [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] criterion where there is no order on the variables.

Another family of approaches associates a preliminary clustering step with the model tting problem. Variables from the same cluster are then averaged to form new representative variables. Among them, we have [START_REF] Bühlmann | Correlated variables in regression: clustering and sparse estimation[END_REF] for the Cluster Representative Lasso and Park et al. (2006).

In the graph inference problem, to overcome this, in addition to sparsity, several previous works attempt to estimate CIG by integrating clustering structures among variables for the sake of both statistical sanity and interpretability. A nonexhaustive list of works that integrate a clustering structure to speed up or improve the estimation procedure includes [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF]; [START_REF] Ambroise | Inferring sparse gaussian graphical models with latent structure[END_REF]; Mazumder and Hastie (2012a); Tan et al. (2015); [START_REF] Yao | Clustered gaussian graphical model via symmetric convex clustering[END_REF]; [START_REF] Devijver | Block-Diagonal Covariance Selection for High-Dimensional Gaussian Graphical Models[END_REF]. Duchi et al. (2012) proposed to penalize prior known groups of variables together via block ℓ 1 penalties. In [START_REF] Ambroise | Inferring sparse gaussian graphical models with latent structure[END_REF], the clustering structure is uncovered via a mixture model. Dierent penalty levels are then used for each group. The penalization is lowered when the variables belong to the same cluster and increased in the opposite case. Tan et al. (2015) derives from the block diagonal screening rule a link between graphical Lasso and single linkage hierarchical clustering. They then propose to adjust the penalization parameters to each cluster of connected components. A similar approach is presented in [START_REF] Devijver | Block-Diagonal Covariance Selection for High-Dimensional Gaussian Graphical Models[END_REF] where clusters are selected in a non-asymptotic fashion. Other two-stage methods include Marlin and [START_REF] Benjamin | Sparse gaussian graphical models with unknown block structure[END_REF]. [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF] infer GGMs graph with a prior knowledge of local neighborhood called local constancy. They suppose that when X i is neighbor to X j , then X k spatial neighbor of X i is likely to be neighbor to X j . The same goes for non-neighbors nodes. They optimize the following cost function:

max Ω∈S p ≻0 log(det(Ω)) -tr(SΩ) -λ 1 ||Ω|| 1 -λ 2 ||D ⊗ Ω|| 1 , (2.11)
where λ 1 and λ 2 non negative penalty weights, D ∈ R q×p the matrix of prior neighbourhood relationships. The problem (2.11) is solved using a coordinate descent-like algorithm. They showed that their method outperforms the graphical Lasso of [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], covariance selection of [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF] and neighborhood selection of [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] in easy and hard settings. [START_REF] Ganguly | Local neighborhood fusion in locally constant gaussian graphical models[END_REF] extends the work of [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF] to the neighborhood selection case. Their neighborhood-fused Lasso consists in estimating the regression coecients through βi (λ 1 , λ 2 ) = arg min

β i 1 n ||X i -X -i β i || 2 + λ 1 ||β i || + λ 2 ||D i β i || 1
(2.12)

In a recent work, [START_REF] Lin | Estimation of sparse gaussian graphical models with hidden clustering structure[END_REF] proposed to simultaneously estimate a precision matrix and uncover a clustering structure for variables. They maximize the following function:

max Ω∈S p ≻0    log det(Ω) -tr(SΩ) -λ 1 i<j |Ω ij | -λ 2 i<j s<t |Ω ij -Ω st |    (2.13)
The term s<t |Ω ij -Ω st | do the pairwise dierences between all the concentration coecients.

In the next section, we introduce convex clustering. This is a tool that can be coupled with the Gaussian graphical model inference problem in order also to recover clustering and graph structures simultaneously.

Convex Clustering

As mentioned earlier, LASSO-type regularization methods commonly used for estimating conditional independence graphs are known to struggle in cases with highly correlated variables [START_REF] Bühlmann | Correlated variables in regression: clustering and sparse estimation[END_REF]Park et al., 2006). To address this limitation and improve the estimation process, some researchers have proposed incorporating clustering structures into the analysis. Various methods have been proposed for this purpose, including those described by [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF], [START_REF] Ambroise | Inferring sparse gaussian graphical models with latent structure[END_REF], Mazumder andHastie (2012a), Tan et al. (2015), [START_REF] Devijver | Block-Diagonal Covariance Selection for High-Dimensional Gaussian Graphical Models[END_REF], and [START_REF] Yao | Clustered gaussian graphical model via symmetric convex clustering[END_REF].

Clustering is an approach that aims to detect group patterns in data, and various algorithms exist for this purpose. A comprehensive review can be found in books such as [START_REF] John A Hartigan | Clustering algorithms[END_REF]. Here, we will focus on hierarchical agglomerative clustering (HAC), k-means, and their convex relaxation.

Hierarchical Agglomerative Clustering

Hierarchical Agglomerative (or bottom up) Clustering [START_REF] Stephen | Hierarchical clustering schemes[END_REF]) is a partitioning algorithm of a data set into successively big clusters. Clusters at each level of the hierarchy are obtained by the merging of lower levels clusters. The process starts with n observations or clusters and ends with one cluster. This can be achieved by dening a dissimilarity between clusters also known as linkage function.

Contrarily to k-means or k-medoids algorithms, there is no need to specify rst the number of clusters. Indeed, the cluster structure is recovered at multiple levels of granularity.

Some common linkages functions are given below. Let A and B be two clusters. Denote d(A, B) their dissimilarity and d(x i , x j ) the distance between observations x i ∈ A and x j ∈ B.

In single linkage,

d(A, B) = min {d(x i , x j )} , ∀x i ∈ A, x j ∈ B (2.14) In complete linkage, d(A, B) = max {d(x i , x j )} , ∀x i ∈ A, x j ∈ B (2.15) In average linkage, d(A, B) = 1 n A n B n A i=1 n B i=1 d(x i , x j ) (2.

16)

where n A , n B are the number of observations in clusters A and B respectively. The Ward's distance is given by

d(A, B) = n A n B n A + n B d(g A , g B ) 2
(2.17)

where g A and g B are the barycenters of clusters A and B respectively.

In order to reduce the computational cost of the HAC algorithm, some dissimilarity update formulas between clustering tree levels have been proposed. Among them, the Lance-Williams formula (Lance and Williams, 1967) which allows to compute distance between a newly formed cluster and the other observations as a function of previous level distances.

The result of HAC can be represented in the form of a tree called dendrogram whose leaves are observations and nodes intermediate clusters.

HAC is an attractive algorithm in term of ease of interpretation of the hierarchical clustering tree. However, the clustering errors made at lower levels remain, as one's go up in the hierarchy. Moreover, the algorithm is highly sensitive to perturbations in the input dataset [START_REF] Eric | Recovering trees with convex clustering[END_REF] and is dened in an iterative way without cost function.

K-means clusetring

Let X = {x 1 , . . . , x n } be a dataset of n observations of a random variable x. In k-means [START_REF] Lloyd | Least squares quantization in pcm[END_REF][START_REF] Macqueen | Classication and analysis of multivariate observations[END_REF], the objective is to nd the assignment of the unlabeled n observations to the desired number of clusters K, and recover a set of centroids vectors {c k , k = 1, . . . , K}. Let r ik be the binary variable which indicates if point x i belongs to cluster k. We seek to minimize the following cost function known as the inertia, with respect to r = {r ik } and c = {c k }:

J(r, c) = n i=1 K k=1 r ik ∥x i -c k ∥ 2 2 .
(2.18)

This criterion ensures that each data point x i is assigned to its closest centroid c k . Remark that dissimilarity measures dierent from the Euclidean distance can be used.

The problem (2.18) can be solved using alternating minimization. First, for each centroid, the set of closest observations is determined. Secondly, the centroid is updated by the mean of assigned observations as follows

c k = 1 n k i∈A k x i ,
where n k = n j=1 r jk and A k is the set of observations indices composing cluster k. The operation is repeated until the algorithm converges. Note that the previous algorithm is closely related to the Expectation-Maximization algorithm [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] for Gaussian Mixture Models (GMMs). As a reminder, a GMM assign observations to overlapping clusters while dening a probabilistic model and has the following form:

p(x|π 1 , . . . , π K , µ 1 , . . . , µ K , Σ 1 , . . . , Σ K ) = K k=1 π k p k (x) (2.19)
where the k-th mixture component p k is the normal distribution with mean µ k and covariance matrix Σ k , π k are weights with k π k = 1. K-means can be seen as a special case of GMM [START_REF] Murphy | Probabilistic machine learning: an introduction[END_REF] in which Σ k = I the identity matrix and π k = 1/K.

When c k is dened as the observation in X which average dissimilarity to the cluster's observations is minimal, it's called a medoid. A partitioning around medoids algorithm [START_REF] Lkpj Rdusseeun | Clustering by means of medoids[END_REF]) also known as k-medoids can thus be derived and is more robust to outliers.

Despite being a popular clustering algorithm, K-means algorithm does not necessarily provide global optima results and is sensible to the algorithm initialization values. Moreover, it might be subject to instabilities due to its non-convex objective functions. Some theoretical results of the method have been surveyed in [START_REF] Steinley | K-means clustering: a half-century synthesis[END_REF].

Convex relaxation of k-means and HAC

The convex relaxation of k-means and HAC [START_REF] Pelckmans | Convex clustering shrinkage[END_REF][START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF][START_REF] Lindsten | Clustering using sum-of-norms regularization: With application to particle lter output computation[END_REF] also known as sum of norms clustering can be formulated as follows. Given X = {x 1 , . . . , x n } ∈ R n×p , we look to minimize with respect to the centroids matrix α ∈ R n×p the criterion

1 2 n i=1 ∥x i -α i ∥ 2 2 + λ i<j w ij ∥α i -α j ∥ q (2.20)
where λ is a sparsity penalization parameter, {w ij } are symmetric positive weights, α i ∈ R p is the centroid to which observation x i is assigned to, and ∥.∥ q is the ℓ q -norm on R p with q ≥ 1. The reader is refered to [START_REF] Lindsten | Clustering using sum-of-norms regularization: With application to particle lter output computation[END_REF] for the proof of the link between the formulation (6.2) and k-means. The relation between HAC and (6.2) is proved in [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF]. Indeed clustering methods like k-means and HAC as outlined by [START_REF] Lindsten | Clustering using sum-of-norms regularization: With application to particle lter output computation[END_REF]; Chi and Lange (2015); [START_REF] Radchenko | Convex clustering via l1 fusion penalization[END_REF] and [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF] have a greedy strategy.

A simple approach like convex clustering solves the clustering problem for a grid of penalization parameters λ independently and thus provide a global optimum.

The optimization problem being convex, the nal solution is independent from the initialization.

Unlike sparsity in GGMs used for features selection, in convex clustering sparsity allows to determine a clustering structure. The data attachment term (rst term) allows centroids to be kept close to the observations that make up their cluster.

The penalization term ∥α i -α j ∥ q also knowned as fusion term is a fused-group Lasso like penalty [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF][START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] when q > 1 and fused-Lasso term when q = 1. It encourages the centroids to be sparse in their dierences. The common used norms are the ℓ 1 and the ℓ 2 norms. When q = 1, the regularization term tends to produce sparse entries in the vectors α i . In the case q = 2, sparsity is observed in the whole vector instead of its entries.

The penalty λ controls the tradeo between the model t and the number of clusters. Two observations x i and x j belong to the same cluster when their estimated centroids are identical ie αi = αj . Usually, strict equality is not required for the assignement to clusters. A fusion threshold can be dened instead. When λ increases, the fusion strength increases too and centroids tend to fuse together.

For a value of λ large enough, all the clusters merge into a single one. Now let's consider the regularization path of solutions also so called clusterpath in [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF] obtained while convex clustering (6.2). Like for HAC, a

dendrogram can be recovered under some conditions. Recovering tree structures is closely related to the type of norm and weights used. Some theoretical properties and algorithms derived for hierarchical convex clustering problems will be addressed in the following sections.

Convex hierarchical clustering

While solving problem (6.2) for a grid of λ values, it may happen that observations that were once in the same clusters for a given grid value, split for others.

Some authors have studied conditions under which a tree structure without splits is obtained. These conditions are highly dependent on the weights {w ij }.

When weights are identity ie w ij = 1, [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF] proved that a tree structure is recovered in the ℓ 1 -norm space . Chiquet et al. (2017a) extended the conclusions of [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF]. They showed no split occurs for arbitrary norms ℓ q , q = {1, . . . , ∞}.

Theorem 2.2. (Chiquet et al., 2017a) The solutions path S(λ, w) of the problem (6.2) contains no splits when q ∈ {1, . . . , ∞} and w ij = 1.

Theoritical guarantees on cluster recovery with identity weigths have been derived by [START_REF] Panahi | Clustering by sum of norms: Stochastic incremental algorithm, convergence and cluster recovery[END_REF]. Tan and Witten (2015) studied some statistical properties and thus proved a link between single linkage clustering and convex Figure 2.5: Eect of weights w ij choice on the clustering path of 10 data points using distance decreasing weights of the form exp(-γ ∥x i -x j ∥ 2 2 ) [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF]. Solutions paths while varying λ the fusion penalty, are observed for ℓ 1 , ℓ 2 and ℓ ∞ for two dierent weights choice. For norms dierent from ℓ 1 , fusions are abrupt when γ = 0 i.e. w ij = 1 (rst-row graphs). Adding weights improves the tree structure (second-row graphs).

clustering. They also derived upperbounds of penalty parameter value λ for which all the clusters merge.

The use of distance decreasing weights have also been treated in Chiquet et al. (2017a); [START_REF] Eric | Recovering trees with convex clustering[END_REF]. Figure 2.5 from [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF] illustrates the eect of a distance decreasing weight on the solutions path of a 10 observations convex clustering problem.

Theorem 2.3. (Chiquet et al., 2017a) The solutions path S(λ, w) of the problem (6.2) contains no splits when q = 1 and w ij = f (|x i -x j |) with f , a decreasing positive function. [START_REF] Sun | Convex clustering: Model, theoretical guarantee and ecient algorithm[END_REF] derived conditions for perfect recovery in the general weighted convex clustering model.

A variety of algorithms can be used to solve the convex clustering problem. [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF] proposed a homotopy algorithm to recover the solutions path of problem (6.2) for ℓ 1 -norm with identity weights. This is inspired from the general fused Lasso path algorithm proposed by [START_REF] Hoeing | A path algorithm for the fused lasso signal approximator[END_REF]. A more general subgradient descent algorithm has been proposed by the same authors for ℓ 1 , ℓ 2 norms spaces with arbitrary weights. A Franck-Wolf algorithm is used for ℓ ∞ -norm space. Chi and Lange (2015) used two methods which are Alternating direction method of multipliers (ADMM) and Alternating minimization algorithm (AMA).

In a recent work, [START_REF] Sun | Convex clustering: Model, theoretical guarantee and ecient algorithm[END_REF] used a semismooth Newton based augmented Lagrangian method. This approach performs better than the previously mentioned methods.

Convex clustering can thus allow, to achieve clustering and estimation simultaneously. This can motivate the association with graphical modelling approaches to recover simultaneously networks and clustering partition.

We thus made a review of the litterature on inference of GGM with structure constraints. In the following of the manuscript, we contribute to bring some light on the eld by connecting graphs inference with convex clustering.

Non-smooth convex optimization

Convex clustering and Gaussian graphical model inference, especially neighborhood selection, rely on optimizing a convex function that involves the sum of a smooth term (e.g., squared error loss) and non-smooth penalty functions (e.g., LASSO, group-fused LASSO, or Elastic-net). The problem of optimizing nonsmooth cost functions with structured parsimony has attracted some attention over the last few years. We refer to [START_REF] Bach | Convex optimization with sparsity-inducing norms[END_REF] for an exhaustive survey on the topic.

Dene the following optimization problem

minimize x ∈ X J (x)
(2.21)

where J : X → R is a cost function with X a generic nite-dimensional continuous parameter space. In non-smooth optimization, the cost function is not continuously dierentiable everywhere. Indeed, there exist points where the gradient of J is not well dened.

Denition 2.9. The gradient of a function f : R n → R at a point x is the vector of its partial derivatives:

∇f = ∂f ∂x 1 , . . . , ∂f ∂x n ⊤ .
(2.22)

We focus on optimization problems where the non-smooth cost function can be partitioned as the sum of smooth terms J s (i.e., dierentiable) and non-smooth terms J ns :

J (x) = J s (x) + J ns (x).
(2.23)

Moreover, we require the function J to be convex.

Denition 2.10. Let C be a convex set. For any x, x ′ ∈ S, we have

γx + (1 -γ)x ′ ∈ C, for all γ ∈ [0, 1].
(2.24)

Denition 2.11. A function f (x) is considered to be convex when it is dened on a convex set and if, for any x, y ∈ C, and for any 0 ≤ γ ≤ 1, we have:

f (γx + (1 -γ)y) ≤ γf (x) + (1 -γ)f (y).
(2.25)

The function f is strictly convex when the inequality (2.25) is strict.

Denition 2.12. The function f is strongly convex with parameter c > 0 if for all x, y in R n , the following holds:

(∇f (x) -∇f (y)) ⊤ (x -y) ≥ c ∥x -y∥ 2 2 .
(2.26)

Strong convexity implies strict convexity. In the next sections we review some optimization methods that can tackle the problem dened in equation (2.21).

In the following, we consider three class of approaches that can be used for unconstrained non-smooth optimization problems: subgradient methods, proximal algorithms and smoothing techniques.

Subgradient methods

Historically, the subgradient approaches were the rst to be used for nonsmooth optimization problems [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]. The usual gradient based methods for smooth optimization become ineective.

Denition 2.13. Let f : R n → R be a convex function. The subgradient of f at x is any vector v that satises the inequality

f (y) ≥ f (x) + v ⊤ (y -x) for all y.
The function f is said to be subdierentiable at x if there exist at least one subgradient at x. The set of the subgradients is denoted ∂f (x). To minimize f, at iteration k, the subgradient method takes a step of size η k > 0 in the opposite direction from v (k) :

x

(k+1) = x (k) -η k v (k) .
(2.27)

When f is dierentiable, v (k) = ∇f (x (k) ).
Let f ⋆ denote the optimal value of the minimization problem. Denote {η k } the sequence of step sizes. Several methods exist for its choice [START_REF] Boyd | Subgradient methods[END_REF].

Some of them are given below with their convergence results.

Constant step size: η k = η with η > 0. The algorithm is guaranteed to converge within some range of the optimal value. When f is dierentiable, the algorithm converges to the optimal value.

Constant step length:

η k = d v (k) 2 -1
with d > 0. The algorithm is guaranteed to converge within a range of the optimal value. Square summable but not summable:

∞ k=1 η 2 k < ∞ and ∞ k=1 η k = ∞.
The algorithm is guaranteed to converge to the optimal value.

Nonsummable diminishing: lim k→∞ η k = 0 and ∞ k=1 η k = ∞. The algorithm is guaranteed to converge to the optimal value.

Proximal methods

The proximal algorithms are a class of algorithms generally used for nonsmooth, constrained, or distributed problems [START_REF] Parikh | Proximal algorithms[END_REF], typically when the objective function can be decomposed as a sum of a smooth and a non-smooth term. These approaches solve the convex optimization problem while using proximal operators to approximate the non-dierentiable components of the objective function.

Denition 2.14. Let f : R n → R be a closed proper convex function. The proximal operator prox of f is dened as

prox λf (ν) = arg min x f (x) + 1 2λ ∥x -ν∥ 2 2
(2.28)

The parameter λ > 0 controls to which extent the proximal operator maps a given point towards the minimum of f. Proximal operators can be evaluated for various functions for which an analytical solution exists to the convex optimization problem. When a closed-form solution is unavailable for the proximal operator, generic optimization algorithms like the subgradient or the gradient method can be used. In the following, we recall the proximity operators for some common norms functions.

Let f be the ℓ 2 norm in R n . Its proximal operator is given by the block-wise soft-thresholding:

prox λf (v) = 1 - λ ∥ν∥ 2 + ν.
The proximal operator of the ℓ 1 norm writes as follows for all j = 1, . . . , n:

(prox λf (v)) j = 1 - λ ∥ν j ∥ 1 + ν j .
This is also known as element-wise soft thresholding.

The sum of norms function f = g∈G ∥x g ∥ 2 where G is a partition of {1, . . . , p} admits the following proximal operator:

prox λf (v) = 1 - λ ∥ν g ∥ 2 + ν g .

The alternating direction method of multipliers

The ADMM is an approach dedicated to composite objective function optimization. It solves problems for which the objective function can be decomposed into multiple smaller subproblems, simpler to solve, and coupled through some constraint. Consider the problem minimize f (x) + g(y)

subject to x = y.

(2.29)

where f, g : R n → R are non-smooth (and smooth) convex functions. The ADMM combines an alternate minimization of the augmented Lagrangian with respect to variables x and y with an update of the Lagrange multiplier.

The augmented Lagrangian is dened as the Lagrangian of problem (2.29) augmented with a penalty term for the constraint. It takes the form

L ρ (x, y, u) = f (x) + g(y) + u ⊤ (x -y) + ρ/2 ∥x -y∥ 2 2 ,
(2.30)

where u is the Lagrange multiplier, x and y ∈ R n the primal variables.

The updates of the ADMM can be written as follows:

x

(k+1) = arg min x L ρ (x, y (k) , u (k) ) y (k+1) = arg min y L ρ (x (k+1) , y, u (k) ) u (k+1) = y (k) + ρ(x (k+1) -y (k+1) ).
(2.31) [START_REF] Parikh | Proximal algorithms[END_REF] show the close connexion between the formulation of the ADMM in terms of augmented Lagrangian and the proximal operators. By dening the scaled dual variable z (k) = (1/ρ)u (k) , and λ = 1/ρ it can be shown that the updates (2.31) are equivalent to:

x (k+1) = prox λf (y (k) -z (k) ) y (k+1) = prox λg (x (k+1) -z (k) ) z (k+1) = z (k) + x (k+1) -y (k+1) .
(2.32)

In practice, the algorithm can be stopped when the primal and dual residuals are below a xed tolerance. Following [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], the primal residuals p (k) and dual residuals d (k) are dened as:

p (k+1) = x (k+1) + y (k+1)
(2.33)

d (k+1) = ρ(y (k+1) -y (k) )
(2.34)

Other stopping criteria based on combined residuals can also be used. The combined residual measures both the primal and dual residuals simultaneously and have been shown to lead to more robust solutions (see, e.g., [START_REF] Stanley H Chan | Plug-and-play admm for image restoration: Fixed-point convergence and applications[END_REF]).

Under mild conditions on functions f , g and the unaugmented Lagrangian (see [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]), the ADMM algorithm satises for each iteration, ∀ρ > 0:

Residual convergence: p (k) → 0 as k → ∞. Objective convergence: f (x (k) ) + g(y (k) ) → p ⋆ as k → ∞ where p ⋆ is the optimal value of the problem (2.29). Dual variable convergence: u (k) → u ⋆ as k → ∞ with u ⋆ the dual optimal point.
In some cases, the algorithm can be slow to converge with a xed parameter ρ. [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] proposed the residual balancing technique, which consists of adjusting ρ to keep primal and dual residuals norms within a certain factor of one another. However, proceeding so does not ensure convergence. Other variations of the ADMM have been proposed for accelerated performances (see, e.g., Buccini et al. ( 2020)) for a general framework of ADMM acceleration with a guarantee of convergence. In the case of using the proximal version of ADMM, the algorithm speed also depends on how eciently the proximity operators are solved.

Smoothing methods

Smoothing methods are a set of algorithms that transform non-smooth optimization problems into dierentiable problems. They create a smooth approximation of the original problem, which can then be solved via traditional optimization methods. Some popular smoothing techniques include the Moreau-Yosida regularization, Nesterov smoothing, and Tikhonov regularization. In the following, we will focus on Nesterov's technique.

Nesterov smoothing

Let f be a non-dierentiable convex function f : R n → R of the form:

f (x) = max u∈U u ⊤ Ax -ϕ(u)
(2.35) with ϕ : U → R and U a convex, compact set. Let µ be a positive smoothing parameter, the following smooth approximation is proposed by [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]:

f µ (x) = max u∈U u ⊤ Ax -ϕ(u) -µd(u) (2.

36)

where A is matrix in R m×n and d(u) is a called a prox function. The Nesterov smoothing method consists in adding a regularization term to the objective function f , which is often chosen to be the quadratic function d(u) = 1 2 ∥u∥ 2 2 . Denote u ⋆ µ (x) the optimal solution of problem (2.36), the solution is unique since d(u) is strongly convex with convexity parameter σ > 0 [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF].

Proposition 2.5. The function f µ (x) is well dened, continuously dierentiable, and convex. The gradient

∇f µ (x) = A ⊤ u ⋆ µ (x) (2.37)
is Lipschitz continuous with constant

L µ = 1 µ σ ∥A∥ 2 2 .
(2.38)

with ∥A∥ 2 2 the spectral norm of A.

Proposition 2.6. Denote M = max u d(u), the following inequality holds for

any x ∈ R n f µ (x) ≤ f (x) ≤ f µ (x) + µM.
(2.39)

CONESTA

The CONESTA algorithm is dedicated to a general class of composite functions formulated as the sum of a dierentiable loss function and two regularization functions. The algorithm is based on a Nesterov smoothing approach and dynamically adapts the smoothing parameter according to a duality gap calculation.

Consider the problem

minimize J (x) = f (x) + g(x) + h(x)
(2.40)

where f : R n → R is a smooth convex function and g, h : R n → R are nonsmooth convex functions. An analytical expression of function g proximal operator is assumed to be available. The function h proximal operator is assumed to be expensive to compute or not known and has the form of an l 1,2 group-norm:

h(x) = ϕ∈Φ ∥D Φ x Φ ∥ 2
with Φ the set of possibly overlapping groups of indices, and D Φ a linear operator on the group. Denote h µ (x) the smooth approximation of function h by the Nesterov approach and ∇h µ (x) its gradient. The new smoothed optimization function can be written as:

J µ (x) = f (x) + g(x) + h µ (x).
(2.41)

where

(f + h µ )(x) is the smooth part. Assuming the smooth function f has the form f (x) = 1 2 ∥Ax -y∥ 2 2 where
A ∈ R m×n is a matrix of features and y the response variable, it can be rewritten as a function of Ax with the constraint z = Ax. The Lagrangian function of (2.40)

is thus:

L µ (x, ν) = J µ (x) + ν ⊤ (Ax -z) (2.42)
where ν ∈ R m is the dual variable composed of the multipliers of the equality constraint. The Lagrangian dual function of problem (2.40) is obtained by minimizing the Lagrangian function with respect to x :

H µ (ν) = inf x L µ (x, ν).
(2.43)

For x primal feasible i.e., equality constraint satised, and ν dual feasible i.e. solution to the dual problem, δ µ (x, ν) = J µ (x) -H µ (ν) is called the duality gap between x and ν.

Denition 2.15. The Fenchel conjugate of a function l : R n → R is the

function l ⋆ : R n → R dened by l ⋆ (ω) = sup x∈R n ω ⊤ x -l(x) .
(2.44)

Exploiting the link between the Fenchel conjugate function and the Lagrangian dual function [START_REF] Borwein | Convex Analysis[END_REF], the duality gap can be written as:

δ µ (x, ν) = J µ (x) + f ⋆ (ν) + g ⋆ (ν) + h ⋆ µ (ν).
(2.45)

where f ⋆ , g ⋆ and h ⋆ are the Fenchel conjugates of f, g and h respectively.

Moreover, the duality gap upperbounds the distance to the minimum. The following inequality holds:

δ µ (x, ν) ≥ J µ (x) -J µ (x ⋆ ) ≥ 0.
(2.46)

A zero duality gap is thus equivalent to optimality.

In CONESTA, the duality gap is used as a stopping criterion as it measures the algorithm's progress. The duality gap is calculated at each iteration and compared to a threshold value ϵ. At convergence for a xed smoothing parameter, the obtained approximation satises:

δ µ (x, ν) ≥ ϵ µ (2.47)
The precision is updated dynamically throughout the algorithm and is chosen to be a linear function of the smoothing parameter. The optimization algorithm can be guided toward a better solution by dynamically updating the tolerance level. The smoothing parameter µ is selected to minimize the number of iterations needed to reach the desired precision.

Proposition 2.7. For any tolerance level ϵ > 0, the optimal smoothing parameter which minimizes the worst case bound on the number of iterations is given by

µ opt (ϵ) = -M ∥D∥ 2 2 + (M ∥D∥ 2 2 ) 2 + M L(∇(f )) ∥D∥ 2 2 )ϵ M L(∇(f )) (2.48)
where M is the maximum of the Nesterov smoothing prox function, L(∇(f )) the Lipschitz constant of the gradient of f and D the vertical concatenation of the matrices D ϕ .

The updates of the CONESTA algorithm write:

ϵ (k) µ = ϵ (k) -µ (k) M x (k+1) = arg min x J µ (k) (x) ϵ (k) = δ µ (k) (x (k+1) , ν (k+1) ) + µ (k) M ϵ (k+1) = τ ϵ (k) µ (k+1) = µ opt (ϵ (k+1) )
(2.49) where τ ∈]0, 1[ is the factor by which the distance to the minimum is geometrically decreased. The smoothed minimization problem can be solved using classical proximal algorithms. Thanks to the continuation of Nesterov smoothing in a shrinkage-thresholding algorithm [START_REF] Hadj-Selem | Continuation of Nesterov's Smoothing for Regression with Structured Sparsity in High-Dimensional Neuroimaging[END_REF], the method achieves faster convergence speed than other competing approaches such as alternating direction methods of multipliers (ADMM, [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. A quick empirical convergence analysis is presented for that purpose. We then show how we implemented the approach in practice. The graphs inferred by the MGLASSO are characterized by their multiscale structure, following a grid of fusion regularization parameters. The proposed model selection method based on the stability approach to regularization selection (StARS, [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF]) focuses on selecting the LASSO regularization parameter. We discuss some other selection methods that have been proposed for the LASSO model selection problem. The performances are evaluated in dierent simulation settings, from graph models to phylogenetic tree-based models. The application of MGLASSO on real data is available in the next chapter.

-Multiscale

In Section 3.1 we present the basic framework of MGLASSO. In Section 3.2, we present how the structure of the graphical model is learned. Section 3.3 describes the model selection procedure. In Section 3.4, MGLASSO is compared to some network inference and clustering approaches on synthetic data.

Model presentation

Problem formulation

The Multiscale Graphical LASSO method aims at inferring a graphical Gaussian model while hierarchically grouping variables. It infers conditional independence between dierent groups of variables. The approach is based on neighborhood selection [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] and considers an additional fused-LASSO type penalty for clustering. In the spirit of hierarchical convex clustering, the hierarchical structure is recovered by spanning the regularization path.

Let X = (X 1 , . . . , X p ) ⊤ be a p-dimensional Gaussian random vector, with mean vector µ ∈ R p and positive denite covariance matrix Σ ∈ R p×p . Let G = (V, E) be a graph encoding the conditional independence structure of the normal distribution N (µ, Σ), where V = {1, . . . p} is the set of vertices and E the set of edges. The graph G is uniquely determined by the support of the precision matrix Ω = Σ -1 [START_REF] Dempster | Covariance Selection[END_REF]. Specically, for any two vertices i ̸ = j ∈ V , the edge (i, j) does not belong to the set E if and only if Ω ij = 0. The variables X i and X j are said to be independent conditionally to the remaining variables X \(i,j) . We note,

X i ⊥ ⊥ X j |X \(i,j) ⇔ Ω ij = 0. Let X = (X ⊤ 1 , . . . , X ⊤ n )
⊤ be the n × p-dimensional data matrix composed of n i.i.d samples of the Gaussian random vector X. To perform graphical model inference, [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] consider p separate linear regressions of the form: βi (λ) = argmin

β i ∈R p-1 1 n X i -X \i β i 2 2 + λ β i 1 , (3.1)
where λ is a non-negative regularization parameter, X \i denotes the matrix X deprived of column i, β i = (β i j ) j∈{1,...,p}\i is a vector of p -1 regression coecients and ∥.∥ 1 is the ℓ 1 -norm. These LASSO regularized problems estimate the neighborhoods, one variable at a time. The nal edge set estimates Ê can be deduced from the union of the estimated neighborhoods using an AND or OR rule [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF]. The Meinshausen-Bûhlmann (MB) approach is based on the central relationship between simple linear regression and precision matrix coecients. It can be shown that β i j = - , 1996).

Ω ij Ω ii (Lauritzen
Consider the data matrix X ∈ R n×p without any underlying distribution and the clustering analysis of the p points in R n . The convex clustering problem [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF][START_REF] Lindsten | Clustering using sum-of-norms regularization: With application to particle lter output computation[END_REF][START_REF] Pelckmans | Convex clustering shrinkage[END_REF] is the minimization of the quantity

1 2 p i=1 X i -α i 2 2 + λ i<j w ij α i -α j q (3.2)
with respect to the matrix α ∈ R p×n , where λ is a sparsity penalization parameter, {w ij } are symmetric positive weights, α i ∈ R n is the centroid to which X i is assigned to, and ∥.∥ q is the ℓ q -norm on R p with q ≥ 1. Points X i and X j are assigned to the same cluster if αi ≈ αj . We propose merging the p independent LASSO regressions of the MB approach into a single optimization criterion where a convex clustering fusion penalty in ℓ 2 is applied on the regression vectors considered as cluster centers. Namely, the Multiscale Graphical LASSO (MGLASSO) pseudo-likelihood problem minimizes in a Gaussian framework the following quantity:

J λ 1 ,λ 2 (β; X) = 1 2 p i=1 X i -X \i β i 2 2 + λ 1 p i=1 β i 1 + λ 2 i<j β i -τ ij β j 2 , (3.3) with respect to β := [β 1 , . . . , β p ] ∈ R (p-1)×p , where X i ∈ R n denotes the i-th column of X, λ 1 and λ 2 are penalization parameters, τ ij ∈ R (p-1)×(p-1
) is a permutation matrix, which permutes the coecients in the regression vector β j such as

β i -τ ij β j 2 = k∈{1,...,p}\{i,j} (β i k -β j k ) 2 + (β i j -β j i ) 2 ,
as illustrated in Figure 3.1. The coecient β i k is to be read as the multiple regression coecients of X i on X k .

Grouping eect

The MGLASSO criterion can be seen as a multitask regression problem where the set of responses is identical to the set of predictors. The LASSO penalty term encourages sparsity in the estimated coecients. The group-fused term encourages fusion in the regression coecients β i and β j . An ℓ 1 -norm of the dierences would aect each variable individually [START_REF] Degras | Sparse group fused lasso for model segmentation: a hybrid approach[END_REF]. Moreover, some simulations studies by Tan and Witten (2015) showed that the convex clustering in ℓ 2 outperforms the ℓ 1 case.

The clustering eect in MGLASSO occurs at two scales: rst, by using the norm penalty ℓ 2 to select groups of correlated variables, and second by considering the parsimony of the dierences of the regression vectors. The natural tendency of this type of penalty would be to produce groups of nearly equal regression vectors [START_REF] Zeng | A link-free sparse group variable selection method for single-index model[END_REF]. Using this form of penalization is helpful when no prior information is available on the groups. Intuitively, one can notice that in a framework where the data X has a block-diagonal correlation structure with similar correlation levels, the model excels when these correlation levels are high enough (Park et al., 2006).

Let us illustrate by an example the eect of the fusion term in the proposed approach. Two variables i and j are in the same group when ∥β iτ ij β j ∥ 2 ≈ 0.

Considering a cluster C of q variables, it is straightforward to show that ∀(i, j) ∈ C 2 , we have βi j = β C , where β C is a scalar. Thus the algorithm is likely to produce precision matrices with blocks of constant entries for a given value of λ 2 , each block corresponding to a cluster. In the same vein as Park et al. (2006), a cluster composed of variables that share the same coecients can be summarized by a representative variable. A component-wise dierence between two regression vectors without reordering the coecients would not necesarily cluster variables which share the same neighborhood. The permutation τ ij reoders coecients in such a way that dierences are taken between symmetric coeecients and those corresponding to the same set of predictors. The model is thus likely to cluster together variables that share the same neighboring structure and encourages symmetric graph structures.

(β i , τ ij β j ) =     β i 1 β i 2 . . . β i j . . . . . . β i k . . . β i p β j 1 β j 2 . . . β j k . . . . . . β j i . . . β j p    
(3.4) In practice, when external information about the clustering structure is available, the problem can be generalized into:

min β p i=1 1 2 X i -X \i β i 2 2 + λ 1 p i=1 β i 1 + λ 2 i<j w ij β i -τ ij β j 2 ,
where w ij is a positive weight. In the remainder of the manuscript, we will assume that w ij = 1.

Local constancy

We succintly show the link between the MGLasso model and the local constancy 2009) enforced that by using the following penalty on the precision matrix in addition to the Lasso penalty:

P λ 2 (Ω) = λ 2 ∥D ⊘ Ω∥ 1 ,
where D is a m × p matrix, with m the number of local neigbors and ⊘ the diagonal excluded matrix product [START_REF] Ganguly | Local neighborhood fusion in locally constant gaussian graphical models[END_REF]. The k-th row D k , . = e i -e j where (X i , X j ) are local neighbors and e t a canonical basis vector in R p with the t-th element set to 1.

If the node X 1 is independent (or dependent) of node X 2 , a local neighbor X 1 ′ of X 1 is more likely to be independent or (dependent) of X 2 . Denote G local the domain or prior knowledge graph [START_REF] Ganguly | Local neighborhood fusion in locally constant gaussian graphical models[END_REF]) and E local the associated set of edges. As a remark, the local edges don't necessarily belong to E, the set of edges of graph G. An illustration of the local constancy is given in Figure

3.2. X 1 X 2 X 3 X 4 Figure 3
.2: Illustration of locally constant interactions. The prior knowledge on the neighborhood structure is represented by dashed lines. X 1 and X 2 are local neighbors and mutually connected to X 4 in the true graph (with solid lines). The interaction between {X 1 , X 2 } and X 4 is locally constant. The edge (X 2 , X 3 ) is not locally constant as the local neighbor X 1 is not connected to X 3 .

The dierence matrix D according to the penalty of [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF],

corresponding to Figure 3.2 is D =   1 -1 0 0 0 1 -1 0 0 0 1 -1  
In the MGLasso model with unitary weights, i.e. all w ij = 1, implicitly all nodes are supposed to be local neighbors. Indeed, the fusion penalty includes all the pairwise dierences between the p regression vectors. The local graph G local is hence the complete graph where all the nodes are connected. This assumption might be relevant when no locality information is available. For a model with 4

variables, the local graph is given in Figure 3.3.

X 1 X 2 X 3 X 4
P λ 2 (β 1 , . . . , β p ) = λ 2 i<j β i -τ ij β j 2 .
Let ne(X 1 ) be the markov blanket or neighborhood of node X 1 . A local neighbor X 1 ′ of X 1 is more likely to have the same markov blanket as node X 1 in the MGLasso.

One might notice that when local neighbors are neighbors in the estimated graph Ĝ, the MGLasso model is likely to produce clustering structures where all the variables belonging to the same cluster are mutually connected in term of conditional dependence. The permutation τ ij encourages a symmetry structure.

Model learning

This section introduces a complete numerical scheme to apply MGLASSO in practice using convex optimization algorithms. The objective function in 6 is the sum of three convex components: a smooth function (squared loss) and two nonsmooth penalty functions: the LASSO function, which is separable i.e. can be broken down into independent components involving single entries corresponding to predictors coecients, and the group-fused LASSO penalty, which is nonseparable. Several approaches can be used to tackle its minimization problem.

Among them, we have subgradients methods [START_REF] Zuselevich | Minimization methods for non-dierentiable functions[END_REF]

Optimization algorithms 3.2.1.1 Optimization via subgradient descent

The rst natural approach that can be used to solve the non-smooth optimization problem 6 is the subgradient descent algorithm [START_REF] Boyd | Subgradient methods[END_REF]. The objective in the equation 6 is a non-smooth convex function in β and hence admits a non-unique solution for λ 1 ≥ 0 and λ 2 ≥ 0. That solution is characterized by the following Karush-Kuhn-Tucker conditions, which are the MGLASSO equations for the subgradient: for all i = 1, . . . , p :

∂J λ 1 ,λ 2 (β i ) = -X \i ⊤ X i -X \i β i + λ 1 Sign(β i ) + λ 2 j̸ =i β i -τ ij β j ∥β i -τ ij β j ∥ 2 1 β i ̸ =τ ij β j + γ ij 1 β i =τ ij β j = 0 p , (3.5) with γ ij ∈ R p / ∥γ ij ∥ 2 ≤ 1 and Sign(x) = (Sign(x i )), (1 ≤ i ≤ p) Sign(x i ) =    1 if x i > 0, ∈ [-1, 1] if x i = 0, -1 if x i < 0,
for all x ∈ R p . These equations are sucient and necessary for the solution.

The subgradient method solves problem 6 iteratively. The update at iteration

k is given in algorithm 1 where ν (k) > 0 ∈ ∂J λ 1 ,λ 2 (β (k) ) and η k = 0.1/ √ k is the non-summable diminishing step size i.e. lim k→∞ η k = 0 and ∞ k=1 η k = ∞.
Algorithm 1: Subgradient method update

The algorithm is stopped when ∥β (k+1) -β (k) ∥ 2 ≤ ϵ with ϵ the stopping tolerance. The results obtained through the subgradient method are not always sparse, and the convergence is slow [START_REF] Bach | Convex optimization with sparsity-inducing norms[END_REF].

Proposition 3.1. For diminishing step sizes η k , with lim k→∞ η k = 0 and ∞ k=1 η k = ∞, the subgradient algorithm is guaranteed to converge to an optimal value J ⋆ . Namely,

lim k→∞ J λ 1 ,λ 2 β (k) = J ⋆ .
Proof. Following [START_REF] Boyd | Subgradient methods[END_REF], the proof is based on the fact that we assume J to be a Lipschitz continuous function and hence have the subgradients bounded by the Lipschitz constant.

Optimization via ADMM

A second strategy for solving the optimization problem (6) for xed λ 1 and λ 2 is the Alternating Direction Method of Multipliers (ADMM, [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]).

The approach is adapted to cost functions with a separable structure. We introduce additional variables to account for the dierence between two regression vectors: for (i, j) ∈ {1, . . . , p} 2 , i < j, let ∆ i,j ∈ R (p-1) . We convert the unconstrained optimization (6) to the following constrained problem:

arg min p i=1 1 2 ||X i -X \i β i || 2 2 + λ 1 p i=1 ||β i || 1 + λ 2 i<j ||∆ i,j || 2 such that ∆ i,j = β i -τ ij β j for all (i, j) ∈ {1, . . . , p} 2 .
(3.6)

The augmented Lagrangian function writes, for β ∈ R p(p-1) , ∆ = (∆ i,j ) ∈ R (p-1) p(p-1)/2 , Γ = (Γ i,j ) ∈ R (p-1) p(p-1)/2 and ρ > 0:

L ρ (β, ∆, Γ) = 1 2 p i=1 ||X i -X \i β i || 2 2 + λ 1 p i=1 ||β i || 1 + λ 2 i<j ||∆ i,j || 2 + i<j Γ i,j , ∆ i,j -β i + τ ij β j + i<j ρ/2∥∆ i,j -β i + τ ij β j ∥ 2 2 , (3.7)
where ⟨u, v⟩ denotes the standard scalar product in R (p-1) . ADMM consists of separate iterative updates of the primal variables β and ∆, and of the dual variable Γ. At iteration k + 1, we apply the following update rules described in algorithm 2.

Algorithm 2: ADMM updates

1 while not converged do 2 (i) β (k+1) = arg min β∈R p(p-1) L ρ (β, ∆ (k) , Γ (k) ); 3 (ii) ∆ (k+1) = arg min β∈R p(p-1) L ρ (β (k+1) , ∆, Γ (k) ); 4 (iii) Γ (k+1) i,j = Γ (k) i,j + ρ(∆ (k+1) i,j
β i (k+1) + β j (k+1) );

end 50

Standard optimization procedures can be applied to carry out each one of these updates.

Update (i) boils down to a weighted LASSO problem. Denote b = Vec(β) ∈ R p(p-1) , obtained by stacking the rows of β above each other. Denote also, for i ∈ {1, . . . , p}, J i ∈ R (p-1)×p(p-1) , the block matrix containing p sub-matrices of size (p -1) × (p -1), each null except for the i-th block which is equal to the identity matrix of order (p -1). Up to the terms which do not depend on β, minimizing the augmented Lagrangian with respect to b is equivalent to

b (k+1) ∈ arg min b∈R p(p-1) 1 2 p i=1 X i -X \i J i b 2 2 + i<j ρ 2 ∆ (k) i,j + ρ -1 Γ (k) i,j -(J i -J j )b 2 2 + λ 1 ∥b∥ 1 , (3.8)
which is a standard ℓ 1 penalized quadratic minimization problem. The main issues here are the problem's dimensionality and the non-separability of the quadratic term with respect to the β i , i ∈ {1, . . . , p}. On the other hand, update (ii) corresponds to a classical ℓ 2 penalized problem. Denoting, for (i, j) ∈ {1, . . . , p} 2 , v

(k) i,j = -ρ -1 Γ (k) i,j + (J i -J j )b (k+1) , we solve ∆ (k+1) ∈ arg min ∆ i<j ρ 2 ∆ i,j -v (k) i,j 2 2 + λ 2 i<j ∥∆ i,j ∥ 2 , (3.9)
which is separable across (i, j) pairs and solved in closed form. We obtain the following update rule for all (i, j) ∈ {1, . . . , p} 2 :

∆ (k+1) i,j = 1 -λ 2 ρ -1 ∥v (k) i,j ∥ -1 2 + v (k) i,j .
(3.10)

In contrast to the general residual balancing technique used to stop the ADMM algorithm [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], the criterion we apply is inspired by the work of [START_REF] Stanley H Chan | Plug-and-play admm for image restoration: Fixed-point convergence and applications[END_REF], which sums up the primal and dual residuals. We then stop the algorithm when the combined residual ζ (k+1) ≤ ϵ where ϵ > 0 is the tolerance and

ζ (k+1) = 1 p(p -1) b (k+1) -b (k) 2 + 1 p(p -1) 2 /2 ∆ (k+1) -∆ (k) 2 + Γ (k+1) -Γ (k) 2 . (3.11)
The factors

1 √ p(p-1) and 1 √ p(p-1) 2 /2
account for the fact that the ℓ 2 norms are in R p(p-1) and R p(p-1) 2 /2 respectively.

The MGLASSO problem is convex. It can be shown that the ADMM algorithm converges, i.e., in primal and dual residuals, when the penalty parameter ρ is held constant [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF].

Proposition 3.2. The iterates β (k) in the ADMM algorithm converge to the global minimizer β ⋆ of the multiscale graphical LASSO criterion J λ 1 ,λ 2 (β).

Proof. The proposition is a direct corollary of the convergence of the ADMM algorithm for convex problems' optimization.

Although the ADMM algorithm converges for the MGLASSO optimization problem, it results in a slow convergence in its vanilla version presented afore. The following section introduces a faster algorithm for the MGLASSO optimization problem, which is based on a continuation with Nesterov's smoothing [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]

f ( β) = 1 2 ||Y -X β|| 2 2 + λ 1 || β|| 1 + λ 2 i<j ||D ij β|| 2 , (3.12) where Y = Vec(X) ∈ R np , β = Vec(β) ∈ R p(p-1) , X is a R [np]×[p×(p-1)] block- diagonal matrix with X \i on the i-th block. The matrix D ij is a (p -1) × p(p -1) matrix chosen so that D ij β = β i -τ ij β j .
Notice that the above equation is a multivariate linear regression problem. Denote s( β) = i<j ||D ij β|| 2 the ℓ 1,2 group-norm. A smooth function can approximate this non-smooth penalty with a known gradient computed using Nesterov's smoothing [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]. Given a smoothness parameter µ > 0, the smooth approximation is dened as follows:

s µ ( β) = max α∈K α ⊤ D β - µ 2 ∥α∥ 2 2 , (3.13)
where K is the cartesian product of ℓ 2 -unit balls, D is the vertical concatenation of the matrices D ij and α is an auxiliary variable resulting from the dual reformulation of s( β). Note that lim µ→0 s µ ( β) = s( β).

A Fast Iterative Shrinkage-Thresholding Algorithm (FISTA, [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]) step can then be applied after computing the gradient of the smooth part.

Denote g( β) = ||Y -X β|| 2 2 , the gradient is:

∇ β (g + λ 2 s µ ) = -X ⊤ (Y -X β) + D ⊤ α ⋆ µ ( β), (3.14)
where α ⋆ µ ( β) is the solution to problem (3.13).

The main ingredient of CONESTA remains in determining the optimal smoothness parameter using the duality gap, which minimizes the number of FISTA iterations for a given tolerance ϵ. Let the smooth approximation of the function f ( β)

be:

f µ ( β) = l( X β) + Q µ ( β) (3.15) where l( X β) = ∥ X β -Y∥ 2 2 and Q µ ( β = λ 1 || β|| 1 + α ⋆ µ ( β) ⊤ D β -µ 2 ∥α ⋆ µ ∥ 2 2 .
The duality gap is given by

GAP µ ( β (k) ) = f µ ( β (k) ) + l ⋆ X β (k) -Y + Q ⋆ µ,k -X ⊤ X β (k) -Y (3.16)
where l ⋆ and Q ⋆ µ,k are the Fenchel conjugates.

The specication of µ is subject to dynamic update. A sequence of decreasing optimal smoothness parameters is generated to adapt the FISTA algorithm stepsize towards ϵ dynamically. Namely, µ (k) = µ opt (ϵ (k) ). The smoothness parameter decreases as one gets closer to β ⋆ , a solution of problem (3.12). Since β ⋆ is unknown, the approximation of the distance to the minimum is achieved via the duality gap. Indeed

GAP( β (k) ) ≥ f ( β (k) ) -f ( β ⋆ ) ≥ 0.
The CONESTA routine is spelled out in the algorithm 3 where L(g + λ 2 s µ ) is the Lipschitz constant of ∇(g + λ 2 s µ ), k is the iteration counter for the inner FISTA updates and i is the iteration counter for CONESTA updates.

Algorithm 3: CONESTA solver Input :

f ( β) = g( β) + h( β) + s( β), ϵ > 0, λ 1 > 0, λ 2 > 0, τ ∈ (0, 1) Output: β 1 Initialize β (0) , ϵ (0) = τ GAP µ=10 -8 ( β (0) ), µ (0) = µ opt (ϵ (0) ) ; 2 repeat 3 ϵ (i) µ = ϵ (i) -µ (i) λ 2 d 2 ; 4 /* FISTA */ 5 k = 2 /* new iterator */ 6 β (1) FISTA = β (0) FISTA = β (k) /* Initial parameters value */ 7 t µ = 1 L(g+λ 2 sµ) /* Compute stepsize with L(g + λ 2 s µ ) the Lipschitz constant of ∇(g + λ 2 s µ ) */ 8 repeat 9 z = β (k-1) FISTA + k-2 k+1 ( β (k-1) FISTA -β (k-2) FISTA ) β (k) FISTA = prox λ 1 h (z -t µ ∇(g + λ 2 s µ )(z)) 10 until GAP µ ( β (k) FISTA ) ≤ ϵ (i) µ ; 11 β (k+1) = θ k FISTA ; 12 ϵ (i) = GAP µ=µ i β (k+1) + µ (i) λ 2 d 2 ; 13 ϵ (i+1) = τ ϵ (i) ; 14 µ (i+1) = µ opt (ϵ (i+1) ) 15 until ϵ (i) ≤ ϵ;
Proposition 3.3. The iterates β (k) in the CONESTA algorithm converge to the global minimizer β ⋆ of the multiscale graphical LASSO criterion J λ 1 ,λ 2 (β).

Proof. The proposition is straightforwardly deduced from the convergence theorem of the CONESTA algorithm in the continuation scheme (see Theorem 3 Hadj-Selem et al. ( 2018)).

Empirical convergence analysis

In this section, we present the numerical results for convergence analysis. We use the objective function criterion at β (k) to assess the convergence of algorithms for dierent levels of tolerance for simulated data sets. We show the computation times results for dierent simulation models and the scalability of the algorithms according to the stopping tolerance and the diculty of the problem (number of variables in the data set). i.e f (β (k) ) -f (β ⋆ ) is used to assess convergence. Indeed, the authors proposed in [START_REF] Löfstedt | Simulated data for linear regression with structured and sparse penalties: Introducing pylearn-simulate[END_REF] an approach to simulate a response variable Y and predictors X from a regression vector β ⋆ for linear regression with control over the penalization parameters. This simulation framework does not t the specic case of reformulated Gaussian graphical models in which the equality Y = vec(X) is required.

Details on the simulation models are provided in Section 3.4.1.1. We used an Erdös-Rényi simulation framework parameterized by α = 0.1. The stochastic block model contains 5 blocks of equal size. The scale-free model contains as many expected edges as the number of variables. We set λ 1 = 0.5 and λ 2 = 10 to run the three optimization algorithms. The number of observations is xed at n = 20. The number of variables p is chosen in {10, 20, 40}. The stopping tolerance levels xed are 10 0 , 10 -1 and 10 -2 . The algorithms have been implemented in Python and launched via an R program. The timing results are shown in minutes. We choose the ADMM penalty parameter ρ xed and proportional to the parameter λ 2 . The algorithms are stopped when the stopping tolerance is not reached after a maximum number of iterations equals 10000. algorithm for all stopping tolerance values. However, the algorithm does not converge and does not yield sparse solutions. ADMM is faster than CONESTA for higher stopping tolerance (ϵ = 1 and ϵ = 0.1). However, for lower stopping tolerance (ϵ = 0.01) CONESTA is preferable to ADMM.

The results of our convergence analysis corroborate the more expanded analysis done in [START_REF] Hadj-Selem | Continuation of Nesterov's Smoothing for Regression with Structured Sparsity in High-Dimensional Neuroimaging[END_REF]. The CONESTA is a superior approach for lower stopping tolerances and is the one chosen in practice for the MGLASSO problem optimization.

Practical implementation 3.2.2.1 Path algorithm

In practice, the MGLASSO is implemented as a path-algorithm in which the LASSO penalty is xed. At each evaluation of the criterion 6 for a pair (λ 1 , λ 2 ), the CONESTA solver is applied. In order to determine clustering assignments as the fusion penalty increases, the scheme described in algorithm 4 has been proposed.

If d(i, j) = βiτ ij βj 2 ≤ ϵ f use then the elements of the pair (i, j) are assigned to the same cluster.

The MGLASSO function

The CONESTA solver is loaded from a library that provides optimization methods for regression models with structured and sparse penalties. In order to solve Input : p ∈ N, k ∈ {1, . . . , p}, Output: The MGLASSO routine is called via the function mglasso() to estimate graphs and clusters for a xed LASSO penalty λ 1 and a grid of group fused LASSO regularization parameters λ 2 s. R> mglasso_data <-scale(mglasso_data) R> fit_mglasso <-mglasso(mglasso_data, lambda1 = 0.2*n, lambda2_start = 0.1, fuse_thresh = 1e-3, verbose = FALSE) mglasso() returns a list with the following entries: fit_mglasso$out returns a list for which each element corresponds to a λ 2 value and a clustering level. An element fit_mglasso$out$levelk contains the regression vectors' matrix beta beta and the computed cluster partition clusters with k clusters; fit_mglasso$lambda1 stores the LASSO parameter.

A k 1 l = 1 2 for i ← 1 to p -1 do 3 for j ← i + 1 to p do 4 if i = k then 5 A k [l, (i -1) × p + k] = 1 ; 6 A k [l, (j -1) × p + j] = -1 ; 7 else if j = k then 8 A k [l, (i -1) × p + j] = 1 ; 9 A k [l, (j -1) × p + i] = -1 ; 10 else 11 A k [l, (i -1) × p + k] = 1 ; 12 A k [l, (j -1) × p + k] = -

Clustering path

We explore the clustering path using the function plot_clusterpath(). It returns the algorithm regularization path representation when λ 1 is kept xed, and 59 λ 2 varies. We plot the estimated X at each level onto the two principal components of the input data in Figure 3.7.

R> library(ggplot2) R> library(ggrepel) R> mglasso:::plot_clusterpath(as.matrix(mglasso_data), fit_mglasso)

Figure 3.7: Clustering path

Graphs path

The function plot_mglasso provides the adjacency matrices for each clustering level in Figure 3.8.

R> plot_mglasso(res)

Figure 3.8: Estimated adjacency matrices for dierent levels of λ 2 values 60

Model selection

A crucial question for practical applications consists in dening a rule to select the penalty parameters (λ 1 , λ 2 ). This selection problem operates at two levels: λ 1 controls the sparsity of the graphical model, and λ 2 controls the number of clusters in the optimal clustering partition. First, these two parameters can be dealt with separately: the sparsity parameter λ 1 , chosen via model selection, while the clustering parameter λ 2 varies across a grid of values to obtain graphs with dierent levels of granularity. Secondly, considering λ 1 xed, a model selection rule can be applied to tune λ 2 . Finally, the parameters can be selected simultaneously. The approach privileged in practical applications of MGLASSO is the rst scheme that allows highlighting the multiscale structure in a Gaussian data set. Nonetheless, the selection criterion used for that scheme can be easily extended to others.

Selection of LASSO regularizer

The problem of model selection in graphical models is dicult in the high dimensional case where the number of samples is small compared to the number of variables, as classical Akaike information criterion (AIC, Akaike (1998)) and

Bayesian information criterion (BIC, [START_REF] Schwarz | Estimating the dimension of a model[END_REF]) tend to perform poorly [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF].

In the following, we focus on the StARS stability selection approach proposed by [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF] as suggested by some preliminary tests where we compared the Extended BIC (EBIC, [START_REF] Foygel | Extended bayesian information criteria for gaussian graphical models[END_REF]), the BIC calibrated with slope heuristics [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF], the Rotation invariant criterion implemented in the Huge package [START_REF] Zhao | The huge package for high-dimensional undirected graph estimation in r[END_REF], the GGMSelect procedure [START_REF] Giraud | Graph selection with ggmselect[END_REF], cross-validation [START_REF] Bien | Sparse estimation of a covariance matrix[END_REF] and StARS. The results of the comparative analysis are shown in Section 3.3.2.

The StARS method uses k subsamples of data to estimate the associated graphs for a given range of λ 1 values. For each value, a global instability of the graph edges is computed. The optimal value of λ 1 is chosen so as to minimize the instability, as follows. Let λ

(1) 1 , . . . , λ (K) 1 be a grid of sparsity regularization parameters, and S 1 , . . . , S N be the N bootstrap samples obtained by sampling the rows of the data set X. For each k ∈ {1, . . . , K} and for each j ∈ {1, . . . , N }, we denote by A k,j (X) the adjacency matrix of the estimated graph obtained by applying the inference algorithm to S n with regularization parameter λ (k) 1 . For each possible edge (s, t) ∈ {1, . . . , p} 2 , the probability of edge appearance is estimated empirically by

θ(k) st = 1 N N j=1 A k,j st . Dene ξst (λ (k) 1 ) = 2 θ(k) st 1 - θ(k)
st the empirical instability of edge (s, t) (that is, twice the variance of the Bernoulli indicator of edge (s, t)). The instability level associated with λ (k) 1 is given by D(λ

(k) 1 ) = s<t ξst (λ (k) 1 ) p 2
.

StARS selects the optimal penalty parameter as follows

λ = max k λ (k) 1 : D(λ (k) 1 ) ≤ υ, k ∈ {1, . . . , K} ,
where υ is the threshold chosen for the instability level.

Other selection approaches

The j ∈ {1, . . . , p|Z j = Z i }, ω ij is given by :

ω ii := 1 + ρ(p Z i -2) 1 + ρ(p Z i -2) -ρ 2 (p Z i -1)
;

ω ij := -ρ 1 + ρ(p Z i -2) -ρ 2 (p Z i -1)
.

If α ll were to be equal to one, this construction of Ω would make it possible to control the level of correlation between the variables in each block to ρ. Introducing a more realistic scheme with α ll = 0.75 allows only to have an approximate control.

Multilevel or hierarchically structured models

The motivation of these models is to simulate hierarchically-structured groups of variables to illustrate how the performances of the MGLASSO compare to other standard clustering methods. To do so, we will construct variance-covariance matrices with underlying hierarchical clustering structure, using stochastic phylogeny trees. We rst recall briey the principles of stochastic phylogeny trees, before describing the detailed model. Phylogeny trees were initially used to describe relationships between species and their evolution through time. The tree model describes when speciation occurs, and the evolutionary model describes how the species' quantitative traits evolve. The rst can be chosen in a family available in the ape R package [START_REF] Paradis | Ape: analyses of phylogenetics and evolution in r language[END_REF], and the second is usually a stochastic process, e.g., the Brownian motion or the Ornstein-Uhlenbeck process. The tree model and the evolutionary model are intimately related. When speciation occurs from a specie, the children carry the parent trait and evolve independently. Hence, a correlation or covariance can be dened between the children's traits, and this quantity is known to be proportional to the shared evolution time [START_REF] Bastide | Detection of adaptive shifts on phylogenies by using shifted stochastic processes on a tree[END_REF]. Here we are interested in the covariance structure between the leaves, i.e., at the end of the evolution process, to derive the hierarchically structured covariance matrix.

Let us introduce some terminology about trees. Let G = (V, E) be an undirected graph with edge set E and vertices set V. Let T = (V, F ) be any rooted ultrametric tree whose leaves correspond with V. Ultrametric trees are a class of phylogenetic trees where leaves are equidistant to the root. Each edge f ∈ F has an associated branch length l f . For leaf i, t i = l f = h denotes the distance from the root to the leaf, and pa(i) is the leaf 's parent. For leaves i, j ∈ V , mrca(i, j) denotes their most recent common ancestor and t mrca(i,j) = t i,j is the distance from the root to mrca(i, j). The phylogenetic distance d i,j between the leaves is d i,j = t i + t j -2t i,j . Note that T and G do not share the same edge set. Denote X i1≤i≤Card(V ) a sequence of continuous random variables describing a given trait at each leaf. Let us assume that the branch l f has child leaf i and parent node pa(i). A quantitative trait evolves on this branch according to a stochastic process (W f t , 0 ≤ t ≤ l f ) with distribution P(ω e ), independently from other species, conditionally on W e 0 = X pa(i) . At leaf i, we have X i = W f l f .

A possible choice of the tree model is the coalescent model. We refer the reader to Degnan and Salter (2005) for details on coalescent models. We assume the trait of interest evolves along the coalescent tree under the Ornstein-Uhkenbeck (OU) process. A univariate OU process W t , 0 ≤ t ≤ h is characterized by an optimal value θ. Its stochastic dierential equation can be written as :

dW t = -α(W t -θ)dt + dϵ t
where α is the selection strength. The Brownian motion ϵ ∼ N (0, σ 2 ) where σ 2 is the variance term. Note that when α = 0, the Ornstein-Uhlenbeck process is equivalent to the Brownian motion.

If W 0 is known and xed, it can be shown that

E(W t ) = W 0 e -αt + θ(1 -e -αt )
and

Cov(W t , W s ) = 1 2α
σ 2 e -α|t-s| -e -α(t+s) .

At leaf i in the ultrametric case, we have X i denote the OU process with constant length branch l f = h. The distribution of X i conditionally to the parent node X pa(i) can be written as:

X i |X pa(i) ∼ N X pa(i) e -αh + θ i (1 -e -αh ), 1 2α σ 2 (1 -e -2αh
) .

By integrating the distance constraint on the tree, the covariance between two leaves can be written as:

Σ i,j = σ 2 2α (e -αd ij -e -2αh
). The simulation procedure follows: we rst simulate a coalescent tree. Then with an OU process on its branches, we compute the covariance at leaves. This covariance matrix is nally used to simulate a multivariate Gaussian dataset if its inverse exists. Figure 3.13 shows the covariance and inverse covariance matrices derived with the OU process for a coalescent tree model. The benet of this simulation scheme for the performance evaluation is that the ground truth tree can be cut into any number of clusters to match the number of estimated clusters by the compared methods.

Support recovery

In a small simulation study, we compare our approach to the GLASSO in its neighborhood selection version. The goal is to measure the accuracy of the compared methods for edge detection in sparse graphs using ROC curves. There are 3 scenarios with a number of variables xed to 40 and a varying ratio of sample size over number of variables n p ∈ {0.5, 1, 2}. The data were generated from Gaussian distributions using the models described in Section 3.4.1.1. The stopping tolerance for the MGLASSO method was chosen to be 0.001. As mentioned earlier in Section 3.2.2, the MGLASSO is coded in both R and Python, and the GLASSO is performed using the huge package in R. For the approaches to be comparable, LASSO penalties used for MGLASSO, need to be rescaled by a factor of n in huge implementation. Timings results are not computed. However, GLASSO is much faster than MGLASSO in all the described settings.

In GLASSO and MGLASSO, the sparsity is controlled by a regularization parameter λ 1 , and they are selected through the StARS approach. Moreover, MGLASSO admits an additional regularization parameter, λ 2 , which controls the strength of the convex clustering. To compare the two methods, in each ROC curve, we vary the parameter λ 1 while the parameter λ 2 (for MGLASSO) is kept constant. We have chosen to compare ROC curves for dierent fusion penalty parameters instead of the results for a particular value of λ 2 to highlight the fusion penalty eect slightly. The GLASSO method does not aim at proposing clusters of variables and focuses on the graph inference task. For this reason, we did not derive clusters from its estimated graphs, as the approach chosen for clustering might be subject to extensive discussion. We computed ROC curves for 3 dierent penalty levels for the λ 2 parameter. Since GLASSO does not depend on λ 2 , the GLASSO ROC curves are replicated.

In a decision rule associated with a sparsity penalty level λ 1 , we recall the denition of the two following functions.

The true positive rate is given by T P (λ 1 )

T P (λ 1 )+F N (λ 1 ) . The false positive rate is dened as follows 1 -

T N (λ 1 ) T N (λ 1 )+F P (λ 1 ) , where T P is the number of true positives, T N the number of true negatives, F N the number of false negatives and F P the number of false positives. The ROC curve represents the true positive rate as a function of the false positive rate. For a given level of true positive rate, the best method minimizes the false positive rate. Within each panel, the ROC curve shows the True positive rate (y-axis) vs. the False positive rate (x-axis) for both MGLASSO (blue) and GLASSO (brown). Since GLASSO does not have a fusion penalty, its ROC curves were replicated for panels belonging to the same row. We also plot the random classier (dotted grey line). The results have been averaged over 50 simulated datasets and suggest that MGLASSO performs no worse than GLASSO. For λ 2 = 0, the MGLASSO approach is equivalent to GLASSO in its neighborhood selection version. Figures 3.14a,3.14b,3.14c show the average values of ROC curves for MGLASSO and GLASSO for dierent congurations as averaged over 50 simulations. Based on these empirical results, we rst observe that, in all the considered simulation models, MGLASSO improves over GLASSO in support recovery in the high-dimensional setting where p < n. In addition, in the absence of a fusion penalty, i.e., λ 2 = 0, MGLASSO performs no worse than GLASSO in each one of the 3 models. This conguration corresponds to the baseline and thus suggests that the CONESTA algorithm achieves similar results as the huge algorithm [START_REF] Zhao | The huge package for high-dimensional undirected graph estimation in r[END_REF] for GLASSO neighborhood selection with no added fusion term. The entire regularization path of the MGLASSO according to the fusion penalty is not presented, and dierences between the non-zeros fusion penalty parameter for MGLASSO are minimal. Indeed, in a multiple conguration analysis, the grid of penalty values may not induces the same amount of fusion according to the size of the simulated dataset. A ner grid can be used for each conguration taken independently to handle this problem; however, it would be prohibitively expensive to compute. Note cluster proportions is xed at π = (1/5, . . . , 1/5). Hundred Gaussian data sets were then simulated for each conguration (ρ, n/p xed).The optimal sparsity penalty for MGLASSO was chosen by the Stability Approach to Regularization Selection (StARS) method [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF]. The parameter λ 2 has been varied. The results shown in Figures 3.15b and 3.15a suggest that, particularly at the lower to medium levels of the hierarchy (between 20 and 10 clusters), the hierarchical clustering structure recovered by MGLASSO is comparable to popular clustering methods used in practice. For the higher levels (5 clusters), the performances of MGLASSO deteriorate. As expected, the three compared methods also underperform as the level of correlation inside clusters decreases. Note that the MGLASSO performances may be sensible to fusion threshold ϵ f use set to 10 -6 for this simulation study. Using non-trivial weights could also help improve the overall performance.

The expected empirical evidence according to which the MGLASSO would work fairly well for correlated variables is somehow highlighted with Figures 3.15b -3.15a. However, the result suggests that controling the correlation between groups of predictors used for each node wise regression in the simulation model, is not an easy task in a multitask learning framework where the set of predictors is the same as the set of responses.

As also noticed in [START_REF] Wang | Sparse convex clustering[END_REF], some preliminary tests suggest that the MGLASSO clusters recovery performance is robust to the Lasso penalty parameter used. Higher values for the penalty contributes, however, to increase the strength of the fusion, i.e., the number of variables merging simultaneously. Nevertheless, the performances of the MGLASSO get signicantly improved in the hierarchicallystructured model described in the next paragraph. The signicant dierence with the previous simulation framework is that the oracle partition number is changed via cutting the tree used for simulation, according to the number of estimated clusters by the compared methods.

The convex approaches (MGLASSO, classical convex clustering) do not share the same number of boxplot points as HAC, k-means, and spectral clustering. In a multiple simulation framework, it appears dicult to control or select a desired number of clusters for these approaches. Each boxplot represents the distribution of 130 simulations for HAC, k-means, and spectral clustering, and at least 1 matching simulation for MGLASSO and convex clustering if the results are available. The MGLASSO performance decreases with the number of clusters. So it is easier to nd pertinent clusters at higher levels of the tree than at lower levels. The higherlevel tree performances for MGLASSO are relatively better than the other clustering approaches results. For the classical convex clustering with unitary weights, clusters tend to fuse abruptly, comparatively to MGLASSO. Intermediate partitions with a number of clusters between 1 and 40 are harder to be found.

Another GGM inference method that is based on the theory of convex clustering is the Clustered Gaussian graphical model developed by [START_REF] Yao | Clustered gaussian graphical model via symmetric convex clustering[END_REF]. This approach would be a potential candidate for the evaluation of the performance simultaneously in terms of support recovery and clustering results. However, the approach does not lead to sparse graphs after some preliminary tests. So it would be unfair to compare it to the MGLASSO. An alternative comparison method with a graph and clustering method would be to apply a two-step approach where the graph is rst estimated, and community detection or classical clustering algorithms are then used to infer clusters. However, proceeding in such a scheme diverts from the initial objective of the MGLASSO, which is the simultaneous clustering and graph inference. Nonetheless, a not yet published approach that attempts to achieve the simultaneous task of graph inference and cluster estimation appears to be the work of [START_REF] Lin | Estimation of sparse gaussian graphical models with hidden clustering structure[END_REF]. However, an explicit rule for cluster deduction is not proposed in the paper, nor an explicit link with the theory of convex clustering is established, which leaves it impractical to make an eventual comparison with the MGLASSO.

Inference of microbial networks via MGLASSO

This section presents a second application of the MGLASSO model on real data, particularly microbial abundance data that was collected as part of the American gut project [START_REF] Mcdonald | American gut: an open platform for citizen science microbiome research[END_REF]. The use of neighborhood selection in sparse microbial network inference has already been reported in works like Kurtz et al. (2015). We show here how the MGLASSO model, also belonging to neighborhood selection approaches, can be used for the problem while drawing a parallel with the cited work. Human gut microbes play a critical role in homeostasis, disease, and digestion. Knowledge of their interactions can serve as an initial community-level description of the underlying microbial ecosystem [START_REF] Yoon | Microbial networks in spring-semi-parametric rank-based correlation and partial correlation estimation for quantitative microbiome data[END_REF]. We perform multi-scale network clustering and inference analysis in which we focus on deriving representative variables for microbial species, unlike the previous application of MGLASSO, in which networks were drawn across the whole set of variables in the model. The material of the analysis is described in Section 3.5.1. We briey reintroduce the inference approach proposed by Kurtz et al. (2015) in Section 3.5.2. MGLASSO learning is presented in Section 3.5.3. We show some results in Section 3.5.4.

Material

We used the microbial abundance dataset called amgut1.filt included in the SpiecEasi package (Kurtz et al., 2015). It describes the microbial community structure of n = 289 samples via an abundance table for p = 127 types of microbes (operational taxonomic units). The data has been pre-processed by Kurtz et al. (2015) and is a subset acquired from the rst round of the American gut project (AGP, McDonald et al. (2018)). Samples of feces, skin and other body regions were taken from thousands of participants and proled using 16S rRNA sequencing.

The project aimed to study the associations between the human microbiome and factors such as diet.

SpiecEasi method

We briey recall the SpiecEasi method for the inference of sparse microbial association networks and refer to the seminal article for more details. The technique essentially combines transformation for relative abundance data with a neighborhood selection approach for Gaussian graphical inference [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF]) coupled with the StARS scheme for model selection [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF]. The transformation applied to the data is the centered log ratio and is presented below.

For j = 1, . . . , n, let Y j = (Y j1 , . . . , Y jp ) ∈ N p be the p-dimensional rowvector of number of OTUs observed from the j-th sample with an associated relative abundance vector Y j =

Y j p i=1 Y ij
. Note that { Y j } j=1,...,p belongs to the p-dimensional unit simplex :

S p = x = (x 1 , . . . , x p ) ∈ R p p i=1 x i = 1, x ≥ 0 .
(3.17) Kurtz et al. (2015) then applied the centered log ratio transformation (clr, Aitchison (1982)) to the relative abundances:

clr Y j = log    Y j p i=1 Y ij 1/p    , for all j = 1, . . . , p. (3.18) where p i=1 Y ij 1/p
is the geometric mean of the relative abundance vector. A pseudo-count is added to {Y j } j=1,...,n to avoid numerical problems due to zero counts.

MGLASSO learning

We applied the MGLASSO method introduced in Chapter 3, on the OTUs previously transformed by the clr approach. The MGLASSO with a value of λ 2 = 0 is therefore equivalent to the SpiecEasi method.

For the model tting, the LASSO penalty parameter λ 1 is selected via the StARS approach with a variability threshold set at 0.05 in a grid of values dened according to the same rules as in Section 4.4.3. The fusion penalty λ 2 is chosen in the interval [0, 20] with irregular steps. This was done in an attempt to control abrupt merges in the clustering path and to reduce the number of values to be evaluated due to computation time. We chose 20 equidistant values in each of the following intervals: [0.20], [1.20] and [0.4]. The CONESTA solver precision is set to ϵ = 0.01.

To represent the graphs estimated on the grid of penalization parameters λ 2 , we rst selected the representative variables by clusters to contrast with the representation used in the section 4.4 where the matrices of adjacency on the whole set of nodes have been presented. The representative variable is one of the cluster's variablesthe fusion principle. Indeed, theoretically, variables that belong to the same cluster share the same neighborhood. To compute clustering partitions, the fusion threshold was set to ϵ f use = 0.001.

Results

Using the parameters above, we compute a clustering path of MGLASSO solutions and then display the estimated networks for specic values of the λ 2 fusion penalty.

Figure 3.17 shows how the predicted data computed from the estimated regression vectors scales on the λ 2 fusion penalty parameter grid. The variables were not clearly separated on the rst two principal components. The results are therefore displayed on the components 3 and 4.

Discussion

The analysis presented above was exploratory and aimed not only to illustrate the MGLASSO model but also to show its connection to the SpiecEasi approach when applied in a microbial network inference framework.

Elements of omics data

Omics is a eld of study in biology that involves the comprehensive characterization and quantication of molecules, such as genes, proteins, and metabolites, within a biological system. This can include genomics (study of genes and genomes), proteomics (study of proteins), metabolomics (study of metabolites), and other subelds. The goal of omics is to understand the interactions and relationships between these molecules in order to gain a better understanding of the overall function and regulation of biological systems. The omics data are generated from high-throughput sequencing technologies, also known as next-generation sequencing (NGS). NGS opposed to rst-generation technologies, are new sequencing approaches that read genomes at a higher speed and a relatively cheaper cost.

In this section, we review some concepts of biology and some omics data type that will help understand the remainder of the chapter.

Genomics

The deoxyribonucleic acid of an organism (DNA) is a helix molecule of twopaired strands that contains information about its functioning, development, reproduction, and growth. It is structured into chromosomes found in the cell's nucleus. The genome is the complete DNA sequence of an organism. Its size diers between species. For example, the poplar (Populus trichocarpa) genome has 19

chromosomes for about 0.485 billion base-pairs while the human (Homo sapiens) genome has 42 chromosomes for approximately 3 billion base-pairs. The genomic sequence can be segmented into genes and intergenic regions. Genes can, in turn, be segmented into units called exons and introns. The exons correspond to coding regions in the gene, while introns are the non-coding regions. The expression of genes can be measured as the number of RNA transcripts produced by the gene.

Transcriptomics

Transcriptomics studies the ribonucleic acid (RNA) available in a cell or a set of cells. The transcriptome is the complete set of RNA transcripts, some coding (messenger) and others not (e.g., ribosomal, transfer, or small nuclear).

De novo RNA-sequencing (RNA-Seq) is a breakthrough technology in biology that easily captures any RNA and, therefore, measures genes' expression. It refers to a method of transcriptome assembly where a reference genome is not used to align the RNA-seq reads. Instead, the reads are assembled into contigs and then into longer transcript sequences (isoforms) without the use of a reference genome. De novo assembly is useful for organisms for which a reference genome

is not yet available or for organisms that have a high degree of genetic diversity.

It can also be used to identify novel transcripts and isoforms, which can provide important information about the regulation of gene expression. The de novo RNAseq pipeline includes several steps: RNA isolation, quality control, RNA sequencing, assembly and annotation. This approach is computationally intensive and requires large amounts of computational resources and memory, specially when dealing with large or complex transcriptomes.

In the RNA-Seq experiment, the total numbers of mapped reads, i.e., sequencing depths, are dierent for samples (see Figure 4.2), and more reads would map to longer genomic regions. In order to have comparable counts between these samples, RNA-Seq data are usually normalized via diverse approaches, which might take into account other within and between samples non-uniformities [START_REF] Soneson | A comparison of methods for dierential expression analysis of rna-seq data[END_REF].

Genetics

Population genetics is the study of the variations in DNA sequence between individuals in a population, while quantitative genetics focus on individuals' traits and their heritability [START_REF] Harmon | Phylogenetic comparative methods: learning from trees[END_REF]. The genotype of an organism is its heritable genetic identity.

Multiple variants exist in the genomes of a population of a given species and can be regrouped into 2 signicant categories, which are the single nucleotides variations and the structural variations which extend to multiple nucleotides [START_REF] Manzoni | Genome, transcriptome and proteome: the rise of omics data and their integration in biomedical sciences[END_REF]. The SNVs are also known as singles nucleotide polymorphisms (SNPs, Figure 4.3) when common in a population (frequency higher than 1%).

For most SNPs, two variants are only observed and called alleles [START_REF] Neuvial | Statistical analysis of Single Nucleotide Polymorphism microarrays in cancer studies[END_REF]. Since organisms inherit one copy of each SNP position from each parent, the organism's genotype at an SNP position is either AA, AB or BB with A and B the two alleles [START_REF] Laframboise | Single nucleotide polymorphism arrays: a decade of biological, computational and technological advances[END_REF]. The minor allele frequency (MAF) is the proportion of the less common allele in a population (with a proportion lower than 0.5). Thus, an SNP is characterized by its genomic position, alleles, and MAF. The SNPs' genetic variants can be captured by whole genome sequencing based on NGS technologies.

Epigenetics

Epigenetics studies heritable changes that aect gene expression without changing the DNA sequence [START_REF] Plomion | Forest tree genomics: 10 achievements from the past 10 years and future prospects[END_REF]. The epigenome is the pattern of these changes that mark the genome. For an organism, the genome is usually stable across the cells. However, the epigenome is highly dynamic across cell types and in time (Robinson et al., 2014a). Its proling is thus more complex than proling the genome.

DNA methylation is a commonly studied epigenetic mark. It consists of the The gold-standard method for methylation analysis is whole genome bisulte sequencing (WGBS). The WGBS uses the sodium bisulte treatment to convert unmethylated cytosines to thymine (via uracil) while the methylated cytosines remain untouched [START_REF] Wreczycka | Strategies for analyzing bisulte sequencing data[END_REF]. This allows quantifying the proportion of methylated cytosines over a number of reads (read coverage). When interested in the sequencing of a specic region, approaches such as methyl capture sequencing (MC-Seq, Teh et al. ( 2016)) can be used. They are cheaper than WGBS and need the prior denition of probes to survey the genomic loci.

Microbiomics

Microbiomics studies the microbiome, which is the collection of microorganisms (such as viruses and bacteria) that live in a particular environment 1 . Therefore, the metagenome is the total genomic DNA of the microorganisms, and the metatranscriptome is their total transcribed RNA (Morgan and Huttenhower, 2012).

Microbial communities dier from multicellular organisms such as poplars or humans, which have a unique genomic signature. They do not necessarily carry identical genomes, and it can be impractical to sequence each cell's genome entirely. Some molecular markers (DNA sequences), such as the 16S ribosomal RNA gene, have thus been dened in microbial ecology, which can uniquely identify distinct genomes.

85

16S rRNA gene, which is a conserved gene found in all bacteria and archaea. The 16S rRNA gene is used as a marker to identify and classify dierent bacterial and archaeal species. This method is typically used to study the composition and diversity of bacterial and archaeal communities in a sample, but it does not provide information about the functional genes or metabolic pathways present in the sample. In summary, shotgun metagenomics is a more comprehensive approach that allows for the identication of both known and unknown species and functional genes, while 16S metagenomics is more focused and is used specically to identify and classify dierent bacterial and archaeal species.

The operational taxonomic unit (OTU) is a cluster of microorganisms with a similar marker gene sequence beyond a xed threshold. It replaces species because named species genomes are not always available for specic marker genes. A typical microbiome dataset is the count of sequences per OTUs for a given number of samples. In order to have comparable OTUs between samples, relative counts are usually computed, which gives rise to compositional data.

When interested in identifying SNPs in the microbial sequences, whole DNA sequencing, also known as whole metagenome shotgun (WMS), can be conducted.

Then, the sequenced genome can be compared to reference genomes.

The EPITREE Project

The evolutionary and functional impact of epigenetic variation in forest trees In the framework of the PhD research, we contribute to highlight poplar populations' genetic and epigenetic structure through SNP data and DNA methylation data, respectively. We also show how methylation is used as a marker of population dierentiation by proposing genes whose methylation pattern closely follows the genetic structure of populations. We briey study the link between methylation and expression patterns for a specic class of genes whose expression proles are stable between trees belonging to the same meta-populations. These chosen areas of the analysis result from discussions with experts from the EPITREE project. We also laid the groundwork for further analysis by briey exploring the structure and interplay between transcriptomic and methylation patterns in dierent methylation contexts using the MGLASSO model. The study detailed in Section is summarized in the application section on real data of the MGLASSO article [START_REF] Sanou | mglasso: Multiscale Graphical Lasso[END_REF].

Data and material

Three datasets of dierent omic natures were used for the analysis. These are SNP data, DNA methylation data, and gene expression data. We present them in this section and briey describe some applied preprocessing and transformations.

Note that details of plant material, experimental designs, and bioinformatics preprocessing can be found in Sow ( 2019

SNP data

The SNP data contained 9.7 million SNPs across the 20 genotypes and was collected via the whole genome sequencing (WGS) method. An illustration of the dataset is given in Table 4.1. Missing values are ignored and the SNPs are ltered out, so those with a minor allele frequency (MAF) less than 5% are discarded (see Figure 4.6). This step made it possible to select 5.1 million SNPs. Methylation proportions per cytosine are computed by taking the ratio between the number of methylated reads and the total number of reads in the sample. Filters based on standard deviation [START_REF] Akalin | Computational Genomics with r[END_REF]

Methylation data

Expression data

The raw genes expression data are count data. They were collected via the RNA sequencing method (RNA-seq) on the genotypes of interest. After being normalized with the trimmed mean of M -values (TMM, [START_REF] Mark | A scaling normalization method for dierential expression analysis of rna-seq data[END_REF]), the counts per millions (CPM) are computed. A mixed linear model is then tted on the gene expressions to correct for the eects of cofactors such as date, block, or time of sampling. The best linear unbiased predictors of between and within population random eects are summed to obtain the best linear unbiased prediction for each gene. The process is described in detail in Chateigner et al. where π ij = π(x ij ). When N ij = 1, we return to classical logistic regression for binary variables. We introduce dummy variables to encode the contribution of the categorical variable as follows: dummy(

x ij ) = [1(x ij = 1), . . . , 1(x ij = K)].
Consider the generalized linear model with the link function h

(π) = log π 1-π , with π ∈]0, 1[: logit(π ij ) = log π ij 1 -π ij = K k=1 β jk 1(x ij = k) and k β jk = 0.
The log-likelihood is dened as follows:

log L(β) = n i=1 {ỹ ij log(π ij ) + (N ij -y ij ) log(1 -π ij )} (4.1)
The model without intercept is tted via the following negative convex log-likelihood function, except for one constant:

arg min β j n i=1 -ỹ ij K k=1 β jk 1(x ij = k) + N ij log 1 + exp K k=1 β jk 1(x ij = k) subject to k β jk = 0.
(4.2)

In order to establish whether there is an association between the groups and the response, the deviance of the null model is compared to the deviance of the tted model M via:

G = 2(log Lnull -log L M ),
where Ls and L M are respectively the maximized likelihoods of the null model and the tted model M.

The alternative hypothesis stating H 1j : β jk ̸ = 0 for at least one coecient k was tested against the null hypothesis H 0j : β jk = 0 for all k. The statistic G follows a chi-square probability distribution χ 2 with K degrees of freedom under the null hypothesis. In practice, with a modality of the categorical variable chosen as a reference class, G ∼ χ 2 (K -1). Bonferroni corrections are applied to correct for multiple tests. The dierential analysis procedure described is available in the methylKit package and named weighted fractional logistic regression.

ANOVA model for expression proles

The problem addressed here is the denition of a method to identify dierentially expressed genes from multiple samples. Various tools and methods can be used depending on the nature of the gene expression data [START_REF] Soneson | A comparison of methods for dierential expression analysis of rna-seq data[END_REF]. Regarding counts, methods based on Poisson and Negative Binomial models are the most common. In our case, the data has already been transformed from counts to continuous data, which eases standard methods based on hypothesis testing.

Dierential analysis is conducted for gene expression data by one-way univariate analysis of variance (ANOVA). Let y j = (y 1j , . . . , y nj ) ⊤ ∈ R n be the vector encoding the normalized gene expression values in the n samples. Let X j ∈ {1, . . . , K} n be the n-dimensional categorical covariate with K ≥ 2 categories corresponding to the groups to which the samples belong.

Consider the generalized linear model with the identity link function h(y) = y, where y ∈ R :

y ij = K k=1 β jk 1(x ij = k),
and k β jk = 0 for identiability constraints.

The least squares estimation procedure can be used to estimate the model coecients, which are the solution of:

arg min β j n i=1 y ij - K k=1 β jk 1(x ij = k) 2 subject to k β jk = 0. (4.3)
An ANOVA table for hypothesis testing can be constructed by partitioning the sum of squares and using a Fisher statistic. Denote ŷij the tted values for gene j and observation i and ȳj the mean expression in the sample. The regression residuals can be broken down as follows:

n i=1 (y ij -ŷij ) 2 = n i=1 (y ij -ȳj ) 2 - n i=1 (ŷ ij -ȳj ) 2 .
(4.4)

The statistic

f = n i=1 (ŷ ij -ȳj ) 2 /(K -1) n i=1 (y ij -ŷij ) 2 /(N -K) ∼ F (K -1, N -K), (4.5)
where F (df 1 , df 2 ) is the Fisher distribution with df 1 and df 2 degrees of freedom. The null hypothesis is formulated as follows: H 0j : The gene j has the same expression level in all groups. Bonferroni corrections are applied for the correction of multiple tests. The ANOVA model is tted using the stats package.

Clustering of SNP and methylation data

Identifying the genetic structure of a population of species is common in genetics (Robin and Ambroise, 2019). As part of our analysis, the genetic structure is used to nd epigenetic markers that mimic the same population structure. We briey recall the clustering procedure developed by a co-author of the paper Sow et al. ( 2023) to nd the genetic structure of populations. Hierarchical ascendant clustering with Ward's method is applied to the VanRaden Genomic Relationship Matrix (GRM, VanRaden ( 2008)) converted into a dissimilarity matrix. The R package stats is used.

The epigenetic structure is also uncovered using hierarchical clustering with Ward's method on the correlation similarity matrix between populations. The approach is implemented in the methylKit package.

Gene sets testing

We perform enrichment analysis to provide a biological interpretation of differentially methylated cytosines. We seek signicantly enriched gene cluster in gene ontology terms, recovered from clustering on the cytosines. This approach is slightly dierent from classical gene set enrichment analysis (GSEA), which considers both a set of dierentially expressed or methylated genes and the complementary set.

Cytosine groups are found using k-means with euclidean distances on the methylation proles. The centered log-ratio previously transforms these proles.

Note that the clustering step can make it possible to nd co-methylated proles and help characterize genes whose biological function is not yet known [START_REF] Rau | Statistical methods and software for the analysis of transcriptomic data[END_REF].

The number of clusters is selected using the slope heuristic [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF] with a number of clusters varying between 1 and 30. Clustering is done using the coseq package [START_REF] Rau | coseq[END_REF] initially developed for clustering RNA-Seq data.

In the second step, a map is established between the positions of the cytosines and the known genes of poplar using the genome of Populus trichocarpa v3.1. The Gene Ontology database (GO, [START_REF] Ashburner | Gene ontology: tool for the unication of biology[END_REF]) is used for functional gene annotations. We limit ourselves to the GO terms of biological processes. Molecular functions and cellular components are not taken into account.

Finally, the enrichment analysis is performed using the topGo package after the clustering step. The tests are based on the number of genes per ontology belonging to the cluster. For cluster k and GO term Y , the following question is answered: Is there an association between the genes of the cluster k and the annotations for the GO term Y ? The contingency table 4.5 can therefore be built.

The p-values are calculated using Fisher's exact (hypergeometric) test. hypothesis H 0 states: the genes of the cluster k are not overrepresented in GO Y.

Only ontologies with a p-value less than 0.05 after Benjamini-Hochberg correction are considered signicant.

Complementary GO enrichment analysis is also performed using Metascape [START_REF] Zhou | Metascape provides a biologist-oriented resource for the analysis of systems-level datasets[END_REF] bioinformatics software with default settings. This approach is based on the sets of dierentially and non-dierentially methylated genes. Additionally, it performs clustering of GO terms and returns the GO with the lowest p-value per cluster. DNA methylation data's epigenetic analysis must be considered a triple analysis according to the methylation contexts (CpG, CHG, and CHH). To make the reading pleasant and relevant, we focus only on the results obtained in the CpG context.

Note that the analysis is readily reproducible in other methylation contexts. The clustering tree obtained on the methylation data in a CpG context is given in Figure 4.9. The structure is dierent from the genetic structure of the population.

However, the two genotypes belonging to the same natural population merged all the same before joining other geographical sites.

The three meta-populations identied earlier are used for the dierential analysis step to identify epigenetic markers of the population's genetic structure and, therefore, implicitly of the geographical structure.

Epigenetic markers of the genetic population structure

To characterize which cytosines dierentiate the three meta-populations based on the genetic structure (Section 4.3.4.1), we performed a dierential analysis for DNA methylation data in the CpG context. Figure 4.10 illustrates the histogram of unadjusted p-values for cytosines. After Bonferroni's correction, 69189 cytosine positions out of 260078, or about 26%, were signicantly associated with the genetic structure. The overall signicance level was set at 0.01.

For the enrichment analysis, we selected the top 1500 of the most dierenti- We conducted an enrichment analysis on epigenetic markers, identifying 8

clusters obtained using the clustering approach described in Section 4.3.3. The mapping step identied 746 genes. Figure 4.12 shows the enriched GO for a selected cluster. Compared to the other clusters, this displayed cluster is enriched in genes involved in biological regulation and response to stimuli, among others.

A gene set enrichment analysis is also done using the Metascape approach in Figure 4.13 and also highlights relevant biological processes.

Methylation and expression proles

The objective is to verify if there are patterns between the expression proles and the methylation proles of the housekeeping and dierentially methylated genes within the meta-populations.

Housekeeping genes are, by denition, genes whose expression level is stable between populations. In order to identify them, we performed a dierential analysis between genotypes belonging to the same meta-populations for the RNA-Seq data and selected the non-dierentially expressed genes. This dierential analysis within For expression data (34229 genes), we identied 34227, 33332, and 32660 nondierentially expressed genes in the rst meta-population (Paglia, Basento), the second (Dranse, Kuhkopf, Rhin, Ticino) and the third (Adour, Val d'Allier, Loire, Ramières), respectively. For DNA methylation data, we obtained 8558, 6154, and

Discussion

The population structure was investigated for epigenetic data (DNA methylation in a CpG context) and genetic data (SNPs). As mentioned earlier in the introductory part of the review, studies exploring genetic and epigenetic variation in trees are less well explored than in annual plants.

The results obtained on the poplar forest tree suggest that, as well as the SNPs data, there is a relationship between DNA methylation and the local adaptation of trees. This result has also been noted by [START_REF] Gina F Lamka | Epigenetics in ecology, evolution, and conservation[END_REF]. Indeed, the hierarchical clustering trees have shown that it is always possible to merge rst the nucleotides adenine, cytosine, guanine, and thymine, respectively. As a reminder, methylation consists of adding a methyl group to a cytosine of the genome and occurs in three contexts (CpG, CHG, and CHH, where H ∈ {A, C, T }). Methylation can be measured on promoter and gene-body regions. In promoter regions, methylation is related to gene silencing, while it is related to tissue-specic expression or alternative splicing in gene-body regions (Sow, 2019).

For the genes methylation data, we used data that have been normalized via the reads by density approach and then passed to a logarithm function log 2 (x + 1) with x ∈ R. The reads by density is the product between the number of methylated reads and the proportion of methylated reads. The approach was dened in Sow et al. ( 2023). RNA-Seq data were normalized via trimmed mean of M -values followed by a correction for experimental design eects via a mixed linear model [START_REF] Chateigner | Gene expression predictions and networks in natural populations supports the omnigenic theory[END_REF].

An average value is calculated for each pair of genotypes per population. Integrating the 10 natural populations on the 7 omic entities gives a set of 70 proles. The analysis was restricted to a group of 151 target genes selected from 24926, which explains the most variability in the data. The gene selection process is described in Section 4.4.2.

Gene selection through sparse PCA

We perform a sparse principal component analysis for variable selection based on singular value decomposition and the LASSO penalty on loading vectors. The most signicant genes, i.e., those for which the loading vector entries are non-zero, are selected. The approach is implemented in the R package mixomics [START_REF] Rohart | mixomics: An r package for `omics feature selection and multiple data integration[END_REF], which is based on a low-rank matrix approximation [START_REF] Shen | Sparse principal component analysis via regularized low rank matrix approximation[END_REF].

Consider a dataset X ∈ R n×p with n observations and p variables. Suppose the columns of X are centered. Vanilla principal component analysis can be performed via the singular value decomposition (SVD) of X. Denote r = rank(X), the SVD writes:

X = UDV ⊤ , (4.6)
where U ∈ R n×r is the matrix of left singular vectors, V ∈ R p×r the matrix of right singular vectors, and D ∈ R r×r is a diagonal matrix whose diagonal elements are ordered singular values. The columns of UD are the projections of the data i.e. the principal components. The columns of V are the corresponding loadings or direction vectors, which are used to interpret the principal component axes. Interpretation is improved when loadings are sparse (see [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF] Section 14.5.5).

The intuition of the sparse PCA approach introduced by [START_REF] Shen | Sparse principal component analysis via regularized low rank matrix approximation[END_REF] is to exploit the link between the SVD and the low-rank approximation method [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] while solving a penalized least-squares regression problem. The matrix X can be approximated by another matrix X (l) of lower rank In the convex clustering problem introduced by [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF], the classic illustration used for the method results is the path of convex clustering solutions over the fusion penalty, projected onto the principal components of the original dataset. In the MGLASSO model, we adopted the same idea with the LASSO penalty xed. However, the solutions, i.e., the estimated regression vectors, do not belong to the same space as the data. Consequently, the path is drawn on the predicted data.

Graphs computed using MGLASSO are symmetrized using the AND rule (Meinshausen and [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] i.e., an edge (i, j) is considered to be present when the estimated coecient of variable i on variable j and the estimated coecient of variable j on variable i are both non-zero. The desired precision for stopping the CONESTA-based optimization algorithm [START_REF] Hadj-Selem | Continuation of Nesterov's Smoothing for Regression with Structured Sparsity in High-Dimensional Neuroimaging[END_REF] is xed to ϵ = 0.01.

Results

We present the clustering path and the estimated adjacency matrices generated via the MGLASSO model. We observe some observations' splits after their fusion. However, there is a clear separation between the proles of natural populations corresponding to RNA-Seq data, gene-body methylation data, and promoter methylation data. In The graphs also suggest some relationships between expression and methylation for some natural populations. As the merging penalty increases, the blocks corresponding to the three methylation contexts merge rst, then follow the upper left block corresponding to the expression data. For λ 2 = 30.94, all natural populations merge into a single cluster and complete graph.

Figure 4.17a shows the adjacency matrix with xed Lasso penalty and λ 2 = 0.

The graphs suggest a higher number of interactions between natural populations belonging to the same omic entity. There are fewer interactions between natural populations from dierent omics entities. Regarding these interactions between proles from dierent omic entities, we observe a clear interaction trend between proles belonging to the same methylation context.

As the fusion penalty increases, in Figure 4.17b we observe a grouping effect in the edges, which suggests the same 3 clusters identied in the clustering path. The interactions inside the 3 groups are irrelevant because almost all the nodes have merged. However, the relationships between the three groups may suggest potentially natural populations whose expression, gene-body methylation, and methylation in the promoter region are correlated.

As the fusion penalty continues to increase, all methylation proles merge in For a maximum fusion penalty value, the clustering eect extends to all nodes in the graph resulting in a complete graph.

Discussion

We conducted an integrative analysis on multiple omics data to model the multi-scale relationships between natural populations of poplar and also illustrate the performance of the MGLASSO method in a context of real data where the results can be proof checked according to the biological context.

We accounted for conditional dependencies within natural populations measuring the same omic entity and conditional dependencies between natural populations belonging to dierent omic entities. It is well known that methylation regulates the expression of genes [START_REF] Neuvial | Contributions to statistical inference from genomic data[END_REF], and some of the results suggest potential poplar natural populations for which there is a strong link between their methylation and expression proles. When calculating the MGLASSO clustering path, the path was not agglomerative, but the method highlighted three pure clusters that correspond to strictly dierent omics entities. These results also suggest that gene- 107

-Conclusion and perspectives

This thesis proposes a novel approach to enhance the inference of Gaussian graphical models by incorporating a structure of the variables. This structure is assumed to be unknown. Multilevel Gaussian graphical models have gained increasing attention in recent years, yet their application has been limited to the inference of networks for single variables and a priori known groups of variables.

There is a growing interest in exploring their potential for analyzing heterogeneous data, especially omics data consisting of variables from dierent probability distributions. In omics projects, such as EPITREE, various types of data, including genetic, genomic, and epigenetic data, have been produced to gain a better understanding of the molecular mechanisms involved in the adaptation of forest trees throughout their lifespan. However, one question that has received less attention is the quantied impact of DNA methylation, an epigenetic mark, on the trees' local adaptation. This question can be addressed by exploring the links that can be established between the dierent omics data. By proposing an approach that can account for the structure of variables in Gaussian graphical models, this thesis contributes to the advancement of statistical inference methods. Furthermore, by demonstrating the potential of multilevel Gaussian graphical models for analyzing heterogeneous omics data, this research has the potential to inform on the molecular mechanisms of forest tree adaptation and improve our understanding of the role of DNA methylation in this process.

In this research, we proposed a novel inference model that combines the theories of convex clustering and neighborhood selection. Convex clustering is a group-fused LASSO penalized convex relaxation of hierarchical clustering that can produce robust results. Meanwhile, neighborhood selection is a LASSO penalized pseudo-likelihood approximation that reduces the computational burden in estimating graphical models and simplies extensions to multiple probability distributions.

By linking the inference of Gaussian graphical models with the convex clustering theory, we can construct networks with nested structures, similar to the hierarchical tree structure, under certain conditions. While some previous research has combined both theories, the idea of multi-scaling structures has not been clearly demonstrated. Our proposed model overcomes this limitation and provides a more comprehensive framework for analyzing complex networks.

Our proposed model, the Multiscale Graphical LASSO (MGLASSO), aims to improve the interpretability of networks by providing graphs at multiple levels of granularity. This feature is particularly useful in biological data analysis where it's necessary to zoom in and out of inferred networks. The procedure starts by estimating a network based on single variables, followed by the addition of a fuse-group LASSO penalty to encourage similarity in the estimated vectors of parameters and recover networks on groups of variables. The approach uses convex optimization and minimizes the neighborhood selection objective, penalized by a hybrid regularization that combines a sparsity-inducing norm and a convex clustering penalty.

To apply MGLASSO in practice, we developed a complete numerical scheme that includes an optimization algorithm based on CONESTA [START_REF] Hadj-Selem | Continuation of Nesterov's Smoothing for Regression with Structured Sparsity in High-Dimensional Neuroimaging[END_REF] and a model selection procedure based on StARS [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF]. We implemented the approach as an R package based on Python libraries [START_REF] Sanou | mglasso: Multiscale Graphical Lasso[END_REF].

Although the MGLASSO was designed for graph inference, it is also a convex clustering technique. In the literature of convex clustering, this is the rst attempt to use a continuation method based on Nesterov smoothing.

Our simulation results on synthetic and real datasets demonstrated that MGLASSO outperforms GLASSO [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] in network support recovery in the presence of groups of correlated variables. We further illustrated the method's eectiveness through its application on heterogeneous data, which combined transcriptomic and DNA methylation data, as well as on non-heterogeneous metagenomic data. Our approach accounted for heterogeneous data and non-Gaussian data by applying variable transformations to return to the Gaussian framework. In an eort to better understand the impact of DNA methylation on the local adaptation of forest trees, specically the poplar, the results of a joint study suggest that DNA methylation can be utilized as a marker of trees population dierentiation based on genetic structure [START_REF] Mamadou | Epigenetic variation in tree evolution: a case study in black poplar (populus nigra)[END_REF].

The research work has some limitations that could be addressed in future work. For example, the MGLASSO method struggles to produce tree structures without splits, and the calibration of the threshold for merging variables can be challenging. Data transformation can also introduce biases, and there is currently no selection criterion for the fusion penalty parameter. Additionally, the method is most eective with a relatively low number of variables, which may limit its utility in certain applications. In the framewrk of the EPITREE project, it would be valuable to explore the inference of gene regulatory networks from omics data, in order to identify interactions between epigenetic markers of the genetic structure of natural poplar populations. Furthermore, while the proposed R software is functional, it could benet from more documentation to help users fully understand and utilize its capabilities. Overall, this research lays the groundwork for further advancements in this area, and future work can build upon these ndings to address the existing limitations and extend the method's capabilities. In the theory of convex clustering, weights w ij are often used in the fusion penalty to measure the similarity between variables X i and X j in order to improve clustering results. However, as seen in the case of MGLASSO, the convex clustering solution path may display abrupt fusions or splits of clusters when using identity weights. It should be noted that while convex clustering is a convex relaxation of hierarchical clustering, splits are uncommon in standard hierarchical clustering.

Perspectives

Several rules have been proposed for selecting weights that result in a clustering solution path that reects a tree structure without splits. For example, [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF] conjectures that for decreasing weights w ij = exp(-X i -X j 2

2 ) in the case of ℓ 2 fusions, the clustering does not display splits. In addition, Chiquet et al. (2017a) proposed weights with theoretical guarantees that prevent splits in the clustering solution path. Another study by [START_REF] Eric | Recovering trees with convex clustering[END_REF] suggests that if the weights w ij reect a tree structure among the variables, then the expected solution path will exactly reconstruct the tree. However, how to select weights in the MGLASSO framework to ensure paths with no splits is a topic for future research.

Selection of clusters' fusion threshold

As mentioned throughout the work, in practice, exact equality between estimated vectors is not required before assigning variables to the same cluster as done in most research work about convex clustering. For example, in [START_REF] Sun | Convex clustering: Model, theoretical guarantee and ecient algorithm[END_REF], the fusion threshold is arbitrarily set to ϵ = 10 -5 . Dening the fusion threshold based on a rigorous rule can improve the results. While we assumed a transitivity relation by using an approximate test to derive the clusters, in practice, we have no strong guarantee that β iτ ij β j 2 ≤ ϵ and β jτ jk β k 2 ≤ ϵ ensure that β iτ ik β k 2 ≤ ϵ. As suggested by [START_REF] Jiang | On identifying clusters from sum-of-norms clustering computation[END_REF] for the case of vanilla convex clustering, it may be worth exploring how the MGLASSO is sensitive to threshold selection and how it can be selected more rigorously.

Bounds on the regularization parameters

The MGLASSO algorithm is currently implemented as a basic path algorithm in which the LASSO penalty is kept constant and the clustering penalty parameter is varied in a grid of values. Following iterations can benet from the previous iterations' estimations as warm starts. However, this approach could be further improved by considering theoretical guarantees and incorporating them into the algorithm. In Hoeing (2010), a path following algorithm is proposed for the Fused LASSO signal approximator. This algorithm considers the computation of the next hitting times and violation times as fusion parameter values for which the clusters will fuse or split, respectively. By adopting a similar approach for the MGLASSO, the burden of dening the right grid of values for the fusion penalty parameter could be removed, and the computation burden related to a path algorithm could be optimized. In the worst case scenario, if the hitting and violation

Estimation of precision matrix coecients

The MGLASSO algorithm is not based on the full Gaussian likelihood of the model, unlike the graphical LASSO (GLASSO). Instead, it focuses on the neighborhood selection scheme of inference of Gaussian graphical models using a pseudolikelihood approximation. While it would be interesting to aim for an exact estimation of the precision matrix, it is not guaranteed to lead to faster optimization procedures. However, it may overcome uncertainties related to edge estimations.

In neighborhood selection approaches, OR or AND rules are typically applied, and the consistency of the approach for the MGLASSO case has not been proven with theoretical properties. To prove the consistency of MGLASSO, the necessary assumptions need to be explicitly addressed, as done in [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF].

The estimation of the precision matrix can be approached in two ways. The rst method involves directly maximizing the full likelihood criterion, which is given by:

log det Ω -tr(SΩ) -λ 1 j,k |Ω jk | -λ 2 i<j ∥Ω i. -τ ij Ω j. ∥ 2 , (5.1)
where S is the empirical covariance matrix, Ω is the precision matrix, λ 1 and λ 2 are the regularization parameters for the LASSO and clustering penalties, respectively.

The second method is to rst estimate the precision matrix using the vanilla MGLASSO and then estimate the diagonal elements of the precision matrix. Ideas about this approach can be found in [START_REF] Balmand | On estimation of the diagonal elements of a sparse precision matrix[END_REF]. An additional two-step approach involves retting a full likelihood model with the edge constraints found via MGLASSO estimations, as done in Section 17.3.1 of [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF].

In practice, a selection criterion has not been dened for the MGLASSO since our focus was on exploratory applications. However, it may be worth considering how a selection criterion can be dened for the fusion penalty parameter. This could improve the algorithm's performance and make it more ecient in real-world applications.

Mixed graphical models

The initial aim of this research was to develop a class of graphical models that can handle variables following dierent probability distributions. As mentioned earlier, one of the advantages of the pseudo-likelihood formulation of the model is that it can be extended to the class of probability distributions belonging to the exponential family. Some introductory work on this topic can be found in [START_REF] Yang | Graphical models via generalized linear models[END_REF]. For handling heterogeneous data, mixed graphical models can be considered [START_REF] Yang | Mixed graphical models via exponential families[END_REF]. However, when coupled with convex clustering, dening the fusion penalty term may not be straightforward.

Conclusion notes

The discussion sections of the thesis directly mention perspectives on biological data analysis. Although some quality control type data analyses for the EPITREE project were carried out as part of the thesis, they were not reported in this manuscript. Moving forward, the MGLASSO package could be signicantly improved by enhancing the graphic rendering functions and by implementing objectoriented programming in an S4 framework.

-Résumé en français

This chapter is a summary of the present work, written in French.

Chapitre 1 Cette hypothèse de parcimonie ore des avantages à la fois statistiques et computationnels. En simpliant la structure de dépendance entre les variables [START_REF] Dempster | Covariance Selection[END_REF], elle permet de développer des algorithmes ecaces.

L'inférence du graphe d'indépendance conditionnelle dans les modèles graphiques gaussiens non dirigés, implique l'identication du support de la matrice de précision Ω (l'inverse de la matrice de variance-covariance). Pour apprendre ce graphe, plusieurs méthodes ℓ 1 -pénalisées ont été proposées dans la littérature. Une des approches populaires est la méthode de sélection de voisinage (MB, [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF]) fondée sur des régressions successives par variable en utilisant l'opérateur LASSO. Cette méthode se focalise sur l'apprentissage de la structure du graphe. La méthode MB a donné lieu à une longue série de travaux sur les méthodes de régressions successives par variable, y compris des extensions avec diverses formes de pénalités induisant la parcimonie, comme le sélecteur de Dantzig [START_REF] Yuan | High dimensional inverse covariance matrix estimation via linear programming[END_REF] et l'estimateur Clime [START_REF] Cai | A constrained l1 minimization approach to sparse precision matrix estimation[END_REF]. Une autre famille de méthodes d'inférence de graphes d'indépendance conditionnelle parcimonieux estime directement Ω au travers de la minimisation de la log-vraisemblance négative ℓ 1 -pénalisée [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF]. Cette méthode, appelée LASSO graphique (GLASSO, [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]), bénécie de nombreux algorithmes d'optimisation [START_REF] Yuan | Model selection and estimation in the Gaussian graphical model[END_REF]Rothman et al., 2008a;[START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF][START_REF] Hsieh | Quic: Quadratic approximation for sparse inverse covariance estimation[END_REF].

Les [START_REF] Pelckmans | Convex clustering shrinkage[END_REF][START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF][START_REF] Lindsten | Clustering using sum-of-norms regularization: With application to particle lter output computation[END_REF].

Bien que l'utilisation de pénalités de fusion dans l'inférence de modèles graphiques gaussiens ait été largement étudiée, les travaux précédents se sont principalement concentrés sur la vraisemblance complète pénalisée et ont étudié les pénalités de fusion pour renforcer la constance locale des n÷uds du réseau inféré [START_REF] Honorio | Sparse and locally constant gaussian graphical models[END_REF][START_REF] Yao | Clustered gaussian graphical model via symmetric convex clustering[END_REF][START_REF] Lin | Estimation of sparse gaussian graphical models with hidden clustering structure[END_REF]. Cependant le MGLASSO utilise un critère de pseudo-vraisemblance qui est plus ecace d'un point de vue computationnel et permet de proposer des structures multi-échelles. Bien que le critère utilisé soit similaire à celui utilisé dans le clustering convexe supervisé [START_REF] Hallac | Network lasso: Clustering and optimization in large graphs[END_REF][START_REF] Chu | Adaptive convex clustering of generalized linear models with application in purchase likelihood prediction[END_REF], l'approche proposée peut-être vue comme un problème d'apprentissage multitâche [START_REF] Chiquet | Inferring multiple graphical structures[END_REF] Propriétés de Markov et Factorisation Denition (Indépendance Conditionnelle). Soient X A , X B , X C des ensembles de variables aléatoires. X A est indépendant de X B étant donné X C dans une distribution de probabilité jointe P , si et seulement si :

P (X A = x A , X B = x B |X C = x C ) = P (X A = x A |X C = x C )P (X B = x B |X C = x C ) pour toutes les valeurs x = (x A , x B , x C ) ∈ V al(X). On note X A ⊥ ⊥ X B |X C .
Soit G, un graphe non dirigé, déni par un ensemble de noeuds V et un ensemble d'arêtes E, et dont les arêtes ne portent pas de èches directionnelles.

La structure du graphe est liée à l'indépendance conditionnelle par la notion de séparation.

Denition (Séparation). Un sous-ensemble S est dit séparer les ensembles A et B dans un graphe si tout chemin de A à B doit passer par S.

Cela permet d'introduire les trois propriétés de Markov associées aux modèles graphiques non dirigés : Denition (Propriété de Markov Globale). Dans un graphe où S sépare les ensembles A et B, X satisfait la propriété de Markov globale par rapport au graphe G si : Denition (Factorisation). La fonction de densité f de la distribution de probabilité P par rapport à une mesure produit ν est dite se factoriser par rapport au graphe G si elle peut être représentée comme suit :

X A ⊥ ⊥ X B |X S . où X A = (X k ) \ k ∈ A.
f (x) = 1 Z C∈C ϕ C (x C ),

Modèles graphiques gaussiens

Les modèles graphiques gaussiens (GGM), également connus sous le nom de modèles de sélection de covariance [START_REF] Steen | Graphical models[END_REF], sont une catégorie spéciale de modèles graphiques non dirigés utilisés dans des cadres gaussiens. Soit X = (X 1 , . . . , X p ) T ∈ R p un vecteur aléatoire gaussien de dimension p, de moyenne nulle et de matrice de covariance Σ ∈ S p ≻0 , où S p ≻0 désigne l'ensemble des matrices réelles symétriques dénies positives de taille p×p. Certaines propriétés des modèles graphiques spéciques aux distributions gaussiennes sont présentées ci-dessous.

Proposition. La structure d'indépendance conditionnelle de X ∼ N p (0, Σ)

est caractérisée par le graphe G, qui est uniquement déterminé par le support de la matrice de précision ou de concentration Ω = Σ -1 .

Proposition. Les entrées de la matrice de précision sont proportionnelles aux coecients de corrélation partielle.

En eet, la corrélation partielle entre

X i et X j étant donné X\{X i , X j } est égale à -Ω ij √ Ω ii Ω jj
. Le lecteur peut se référer à Lauritzen (1996) pour la preuve.

Le corollaire suivant peut être dérivé du fait que les coecients de corrélation partielle sont directement liés aux coecients de régression.

Proposition 6.1. Dans le cadre du problème de régression X i = X \i β i + ϵ i , où ϵ i est le vecteur résiduel gaussien, le coecient de régression est donné par

β i k = -Ω ik /Ω ii .
La Proposition 6.1 suggère que les modèles graphiques gaussiens peuvent être estimés par une série de régressions, comme cela a été exposé par [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF]. La prochaine section présentera quelques approches d'inférence de modèles graphiques gaussiens.

Inférence des Modèles Graphiques Gaussiens

Soit P une distribution de probabilité inconnue qui se factorise par rapport au graphe G. En utilisant un ensemble d'échantillons indépendants et identiquement distribués (i.i.d) de P , l'objectif de l'apprentissage d'un modèle graphique gaussien est d'estimer les fonctions potentielles qui correspondent le mieux à la distribution [START_REF] Maathuis | Handbook of graphical models[END_REF]. En d'autres termes, on cherche à estimer les arêtes du graphe et les paramètres de la distribution. Nous regroupons les approches d'apprentissage en trois catégories principales : celles sans contraintes de parcimonie, celles basées sur la parcimonie et celles avec des contraintes supplémentaires sur la structure des noeuds.

Estimation du Maximum de Vraisemblance (MLE)

Une approche courante pour estimer un modèle graphique gaussien est d'utiliser l'Estimateur du Maximum de Vraisemblance. L'objectif est de maximiser la fonction de log-vraisemblance concave basée sur les données observées x provenant du vecteur aléatoire X : l(Ω) = n X T X est la matrice de covariance empirique. La solution à ce problème est unique et donnée par ΩMLE = S -1 pour toute matrice S non singulière. Cependant, le MLE n'est pas calculable lorsque le nombre d'échantillons n est inférieur à la dimension p. Même s'il est calculable, il peut être peu performant et donner lieu à des graphes complets. En général, les modèles graphiques gaussiens sont souvent supposés être parcimonieux, ce qui signie qu'un petit nombre de variables interagissent par rapport à toutes les interactions possibles. Une approche, introduite par [START_REF] Dempster | Covariance Selection[END_REF], consiste à estimer la structure du réseau en xant certains éléments de la matrice de précision Ω à zéro, simpliant ainsi la structure des dépendances entre les variables.

Inférence avec des contraintes de parcimonie Suivant l'idée de Dempster (1972), plusieurs auteurs ont proposé des approches pour retrouver la structure de la matrice de précision en utilisant la parcimonie.

Estimation de la pseudo-vraisemblance pénalisée par Lasso [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] Cette approche communément connu sous le nom de Graphical Lasso [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] a été largement étudiée [START_REF] Yuan | Model selection and estimation in the Gaussian graphical model[END_REF][START_REF] Ravikumar | Highdimensional covariance estimation by minimizing l1-penalized log-determinant divergence[END_REF]Rothman et al., 2008b) et étendue à diverses applications [START_REF] Chiquet | A multiattribute gaussian graphical model for inferring multiscale regulatory networks: an application in breast cancer[END_REF][START_REF] Charbonnier | Weighted-lasso for structured network inference from time course data[END_REF][START_REF] Danaher | The joint graphical lasso for inverse covariance estimation across multiple classes[END_REF]Robin et al., 2019). Elle reste un outil puissant pour l'inférence des modèles graphiques gaussiens.

Autres pénalités de parcimonie

Dans la littérature sur les modèles graphiques gaussiens, des pénalités de parcimonie alternatives au Lasso ont été étudiées. En présence de variables fortement corrélées, le Lasso peut être moins performant. Ainsi, diérentes pénalités, telles que la pénalité ℓ 2 , le SCAD (Smoothly Clipped Absolute Deviation), le Lasso adaptatif [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF], le Lasso groupé [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] Inference tout en prenant en compte la structure sous-jacente Dans un problème de régression pénalisée, le Lasso a tendance à sélectionner une seule variable parmi un groupe de variables corrélées [START_REF] Bühlmann | Correlated variables in regression: clustering and sparse estimation[END_REF].

Diverses solutions ont été proposées en utilisant diérentes pénalités de parcimonie.

Parmi elles, l'Elastic-net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]

Clustering convexe

Étant donné un ensemble de données X = {x 1 , . . . , x n } ∈ R n×p et une matrice de centroïdes α ∈ R n×p , l'objectif du clustering convexe est de minimiser le critère suivant :

1 2 n i=1 ∥x i -α i ∥ 2
2 + λ i<j w ij ∥α i -α j ∥ q (6.2) où λ est un paramètre de pénalisation, {w ij } sont des poids symétriques positifs, α i ∈ R p est le centroïde attribué à l'observation x i , et ∥ • ∥ q est la norme ℓ q sur R p avec q ≥ 1. Contrairement à la parcimonie dans les modèles graphiques gaussiens, utilisée pour la sélection de variables, dans le clustering convexe, la parcimonie est utilisée pour déterminer une structure de clustering. Le terme d'attache aux données (le premier terme) garantit que les centroïdes restent proches des observations qui composent leur groupe. Le terme de régularisation ∥α i -α j ∥ q , également connu sous le nom de "terme de fusion", est une pénalité de type fused-group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF][START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] lorsque q > 1 et un terme fused-Lasso lorsque q = 1. Cela encourage la parcimonie dans les diérences entre centroïdes.

Le paramètre de pénalisation λ contrôle le compromis entre l'ajustement du modèle et le nombre de clusters. Deux observations x i et x j sont considérées comme appartenant au même cluster lorsque leurs centroïdes estimés sont presque identiques, c'est-à-dire αi ≈ αj . À mesure que λ augmente, la force de fusion augmente également, ce qui provoque la fusion des centroïdes. Pour une valeur de λ susamment grande, tous les clusters fusionnent en un seul.

Le chemin de régularisation des solutions, également connu sous le nom de clusterpath [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF]

Modèle

Les outils ayant servi à la construction du modèle sont inspirés de la sélection de voisinage via le Lasso introduite par [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] et de la théorie du clustering convexe [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF][START_REF] Lindsten | Clustering using sum-of-norms regularization: With application to particle lter output computation[END_REF][START_REF] Pelckmans | Convex clustering shrinkage[END_REF].

Nous proposons de fusionner les p régressions LASSO indépendantes de l'approche de [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF] en un seul critère d'optimisation avec ajout d'une pénalité de fusion dans ℓ 2 sur les vecteurs de régression considérés comme centres de cluster. Plus précisément, le problème de pseudo-vraisemblance du MGLASSO minimise, dans le cadre gaussien, la quantité suivante :

J λ 1 ,λ 2 (β; X) = 1 2 p i=1 X i -X \i β i 2 2 + λ 1 p i=1 β i 1 + λ 2 i<j β i -τ ij β j 2 ,
par rapport à β := [β 1 , . . . , β p ] ∈ R (p-1)×p , où X i ∈ R n désigne la i-ème colonne de X, λ 1 et λ 2 sont des paramètres de pénalisation, τ ij ∈ R (p-1)×(p-1) est une matrice de permutation, qui permute les coecients dans le vecteur de régression β j comme suit :

β i -τ ij β j 2 =
k∈{1,...,p}\{i,j}

(β i k -β j k ) 2 + (β i j -β j i ) 2 .
Le critère MGLASSO peut être vu comme un problème de régression multi-tâches où l'ensemble des réponses est identique à l'ensemble des prédicteurs. Le terme La sélection du paramètre de régularisation LASSO est faite via l'approche StARS [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF] qui est une méthode de sélection basée sur une technique de ré-échantillonnage, en gardant le paramètre de pénalité group-fused LASSO nul.

Nous avons également discuté d'autres méthodes de sélection qui ont été proposées pour le problème de sélection de modèle de modèles graphiques gaussiens telles que le BIC étendu (EBIC, Foygel and Drton (2010)), la validation croisée [START_REF] Bien | Sparse estimation of a covariance matrix[END_REF] et la méthode GGMSelect [START_REF] Giraud | Graph selection with ggmselect[END_REF].

Performances

Nous avons mené une étude de simulation pour évaluer les performances de (7.1)

This function is concave with respect to the parameters µ and Σ. Now let's derive the normal conditionals.

Given the set of random variables X \j = Z, the conditional distribution of node X j = Y is gaussian [START_REF] Steen | Graphical models[END_REF]. We can partitioned Σ in four blocks

Σ = Σ Y Y Σ Y Z Σ ZY Σ ZZ
and thus have,

Y |Z = z ∼ N (µ Y |Z , Σ Y |Z ) ∼ N (µ Y + (z -µ Z ) T Σ -1 ZZ Σ ZY , Σ Y Y -Σ T ZY Σ -1 ZZ Σ ZY ) (7.2)
Using Schur complement for block inverse matrix and partitionning Ω as we did for Σ, we have the following

Y |Z = z ∼ N (µ Y + (z -µ Z ) T (-Ω Y Z /Ω Y Y ), Ω -1 Y Y )
Let µ = (µ Y , µ Z ) = 0 p for the sake of simplicity. The log-likelihood of the univariate-conditional normal distribution is given by

log p(Y |Z) = log 1 (2πΣ Y |Z ) 1/2 exp - 1 2 (y -µ Y |Z ) 2 Σ Y |Z = 1 2 log(Σ -1 Y |Z ) - 1 2 Σ -1 Y |Z (y -µ Y |Z ) 2 - 1 2 log(2π) = 1 2 log(Ω Y Y ) - 1 2 Ω Y Y y + Ω Y Z Ω Y Y z 2 + const
Let's x Ω Y Y to a constant and denote β = -Ω Y Z /Ω Y Y , we have:

log p(Y |Z) = - 1 2 Ω Y Y (y -Zβ) 2 + const
Penalizing the pseudo-likelihood with a sparsity term and coming back to our original notations amounts to:

pl(Ω; X) -λ 1 ||Ω|| 1 = n i=1 p j=1 -Ω jj 2 (X j i -X -j i β j ) 2 + λ 1 ||β j || 1 = p j=1
-Ω jj 2 ||X j -X -j β j || 2 2 + λ 1 ||β j || 1 (7.3) with diagonal entries of Ω jj not estimated.

  Titre: Intégration de données omiques complexes par inférence de modèle graphique gaussien multiéchelles Mots clés: Sélection de voisinage Classication hiérarchique convexe Modèles graphiques gaussiens Données omiques Optimisation non lisse Résumé: Cette thèse se concentre sur l'inférence de modèles graphiques gaussiens multi-échelles appliqués à des données omiques. Les nombreuses méthodes statistiques existantes pour l'inférence de réseaux supposent généralement que le réseau est parcimonieux (peu d'interactions réelles parmi les interactions possibles) et font parfois l'hypothèse de l'existence d'une structure sousjacente, qu'elle soit connue ou non. Ces a priori permettent d'obtenir un résumé synthétique des interactions présentes entre les variables d'un ensemble de données. Dans un premier temps, nous avons développé une nouvelle approche d'inférence de graphes permettant d'estimer des graphes à plusieurs niveaux de granularité tout en recouvrant une structure de classication hiérarchique sur les variables. Pour cela, nous nous sommes basés sur les techniques de sélection de voisinage et de classication hiérarchique convexe. La fonction de pseudovraisemblance dérivée a été optimisée grâce à une méthode de continuation utilisant le lissage de Nesterov. Dans un second temps, nous avons eectué des analyses de données omiques provenant de populations naturelles de peupliers. Ces analyses ont consisté à étudier conjointement des données omiques de diérentes natures an de mettre en lumière en particulier les mécanismes de régulation entre données épigénétiques et données génétiques. Nous avons également pris en compte le problème de la nature hétérogène des sources de données grâce à des transformations de variables permettant de revenir au cadre gaussien. Title: Integration of complex omics data through Multiscale Gaussian Graphical Models Keywords: Neighborhood selection Convex hierarchical clustering Gaussian graphical models Omic data -Nonsmooth optimization Abstract: This thesis addresses the inference of multiscale Gaussian graphical models with applications to omic data. Most existing statistical methods for inferring networks assume sparsity (few fundamental interactions among possible interactions)

  Figure 2.1 illustrates a 3-states Potts model. This model (Potts, 1952) is a generalization of the Ising model which arose in statistical physics to model interactions between spins of atoms. In the K-states Potts model, each node takes values in a discrete space = {0, 1, . . . , K -1} where K > 2 is an integer. Nodes can only interact with their nearest neighbors, and the model encourages neighboring nodes on the square lattice to be in the same state.

Figure 2

 2 Figure 2.1: Potts model

Figure 2 .

 2 Figure 2.2 shows a DAG that models the transmission of blood type from parents F , M to child C according to their respective genotypes. Human genetic material is encoded in DNA strings stored in 23 pairs of chromosomes. For each pair, one chromosome comes from the father and the other from the mother. Chromosomes can be divided into regions called loci,

Figure

  Figure 2.4: An undirected graph

  zation. Pairwise Markov networks restrict the potential functions to be over single or pairs of variables, which signicantly reduces the computational burden and improves the interpretability of the model during the learning step. Some examples of pairwise Markov networks for which conditional probability distributions of nodes belong to the exponential family include the Poisson model for count data, Ising or Potts models for categorical data, and the Gaussian model for continuous variables.

  presentation . . . . . . . . . . . . . . . 44 3.1.1 Problem formulation . . . . . . . . . . . . . . 44 3.1.2 Grouping eect . . . . . . . . . . . . . . . . . 46 3.1.3 Local constancy . . . . . . . . . . . . . . . . . 47 3.2 Model learning . . . . . . . . . . . . . . . . . . 48 3.2.1 Optimization algorithms . . . . . . . . . . . . 49 3.2.2 Practical implementation . . . . . . . . . . . 56 3.3 Model selection . . . . . . . . . . . . . . . . . . 61 3.3.1 Selection of LASSO regularizer . . . . . . . . 61 3.3.2 Other selection approaches . . . . . . . . . . 62 3.4 Model performance . . . . . . . . . . . . . . . 65 3.4.1 Synthetic data models . . . . . . . . . . . . . 65 3.4.2 Support recovery . . . . . . . . . . . . . . . . 69 3.4.3 Clustering . . . . . . . . . . . . . . . . . . . . 73 3.5 Inference of microbial networks via MGLASSO 76 3.5.1 Material . . . . . . . . . . . . . . . . . . . . . 77 3.5.2 SpiecEasi method . . . . . . . . . . . . . . . . 77 3.5.3 MGLASSO learning . . . . . . . . . . . . . . 77 3.5.4 Results . . . . . . . . . . . . . . . . . . . . . 78 3.5.5 Discussion . . . . . . . . . . . . . . . . . . . . 80 In this chapter, we introduce the Multiscale Graphical Least Absolute Shrinkage Operator (MGLASSO), a novel method to estimate simultaneously a hierarchical clustering structure, and graphical models depicting the conditional independence structure between variables at multiple levels of granularity. Some of the chapter's material has been published in Sanou et al.. We focus on the neighborhood selection framework and the convex clustering theory presented in the previous chapter to propose a convex optimization problem with a hybrid penalty term combining graphical LASSO and group-fused LASSO penalties. The method allows for highlighting common patterns in the data through clustering while computing interactions between those clusters. In the estimated graphs, the variables belonging to the same clusters are likely to share the same neighborhood.

  The regularization path of solutions to problem in (3.2) can be represented as a dendrogram. The path properties have been studied in Chi and Lange (2015) and Chiquet et al. (2017a), among others.

Figure 3 . 1 :

 31 Figure 3.1: Illustration of the permutation between regression coecients in the MGLASSO model.

  notion introduced by Honorio et al. (2009) and derived a sort of local constancy denition in the sense of the MGLasso. Honorio et al. (2009) introcuded a Gaussian graphical model in which the locality information i.e., known interactions in a dataset are taken in account as a prior for learning the graph structure G. Local constancy encourages the search of dependencies between clusters of variables, instead of variables taken individually. Honorio et al. (

Figure 3 . 3 :

 33 Figure 3.3: Illustration of locality graph in the MGLasso model with 4 variables. All the nodes are expected to be local neighbors. The pairwises dierences includes all the existing pairs of variables.

Figure 3 . 4 :

 34 Figure 3.4: Scalability of the algorithms on three types of graph structures. The mean running time in minutes is plotted against the dimension of the data set for dierent tolerance levels.

Figure 3 . 5 :

 35 Figure 3.5: Comparison of the eect of dierent tolerance levels. We plot the convergence as a function of the running time.

Figure 3 .

 3 Figure 3.5 and Figure 3.4 indicate that the subgradient method is the fastest

  the MGLASSO problem, we construct a linear operator A = {A k } k∈{1,...,p} which encodes the structure of the fusion penalty. Instead of constructing the matrix D (see the equation (3.13)) with p(p-1) 2 components, each corresponding to the pairwise dierence structure, the component A k in A encodes the neighborhood structure of the k-th variable. In other words, A = {A k } k∈{1,...,p} is a reorganization of the elements of D = {D ij } 1≤i<j≤p so that the structure of the dierences involving the variable k is stored in the component A k (see Algorithm 5). This step, not directly accessible, is done internally via the following python routine: Python> A_=linear_operator_from_num_variables(p, type_, W_) Algorithm 5: Component k of linear operator A.

  Figure 3.6: Correlation matrix for a 3-block diagonal model with 3 variables per clusters.

  goal is to outline how the model selection methods compare in dierent simulation settings. The performances are evaluated in Erdös-Rényi, scale-free, and stochastic block models. The selection methods compared are cross-validation (CV), CAPUSHE for selection based on slope heuristic, extended Bayesian information criterion (EBIC), GGMSelect, rotation invariant criterion (RIC), and the stability-based approach StARS. We used an Erdös-Rényi simulation framework parameterized by α = 0.1. The stochastic block model contains 5 blocks with equal sizes. The scale-free model contains as many expected edges as the number of variables. Details on the simulation models are provided in Section 3.4.1.1. The number of observations is xed to n = 20. The number of variables p is chosen in {10, 20, 40}. For each simulation setting, data are replicated 200 times. Performances are assessed via the false positive rate (FPR), true positive rate (TPR), and structural Hamming distance. The parameter λ 1 takes values in [0, 1]. The parameter λ 2 is xed to 0. In Figures, 3.

  Figure 3.9: False positive rate, True positive rate, and structural Hamming distance for dierent model selection criteria in a neighborhood selection problem. Measures are taken on 200 data sets simulated with p chosen in {20, 40, 80} and n = 40.

Figure 3 .

 3 Figure 3.10: Adjacency matrix of a stochastic block model dened by K = 5 classes with identical prior probabilities set to π = 1/K, inter-classes connection probability α kl = 0.75, k ̸ = l, intra-classes connection probability α ll = 0.01 and p = 40 vertices.

Figure 3 .

 3 Figure 3.11: Adjacency matrix of an Erdös-Renyi model with probability of connection α = 0.1 and p = 40 vertices.

Figure 3 .

 3 Figure 3.12: Adjacency matrix of a scale-free model with 40 edges and p = 40 nodes.

  We refer the reader to[START_REF] Bichat | Hierarchical correction of p-values via a tree running ornstein-uhlenbeck process[END_REF] and[START_REF] Bastide | Detection of adaptive shifts on phylogenies by using shifted stochastic processes on a tree[END_REF] for details and proofs of the derivations.

Figure 3 .

 3 Figure 3.13: Tree, Covariance and precision matrices of the phylogeny based hierarchical model

  (c) Mean ROC curves for MGLASSO and GLASSO graph inference in the stochastic block model.

Figure 3 .

 3 Figure 3.14: We varied the fusion penalty parameter of MGLASSO λ 2 ∈ {0, 3.33, 10} alongside the ratio n p ∈ {0.5, 1, 2}.Within each panel, the ROC curve shows the True positive rate (y-axis) vs. the False positive rate (x-axis) for both MGLASSO (blue) and GLASSO (brown). Since GLASSO does not have a fusion penalty, its ROC curves were replicated for panels belonging to the same row. We also plot the random classier (dotted grey line). The results have been averaged over 50 simulated datasets and suggest that MGLASSO performs no worse than GLASSO. For λ 2 = 0, the MGLASSO approach is equivalent to GLASSO in its neighborhood selection version.

  (b) Correlation level ρ = 0.3.

Figure 3 .

 3 Figure 3.15: Boxplots of Rand index for HAC, k-means and MGLASSO in the stochastic-block model. The number of estimated clusters {5, 10, 15, 20} vary alongside the ratio n p ∈ {0.5, 1, 2}. Within each panel, the boxplots of ARI between true partition (with 5 classes) and estimated clustering partitions on 100 simulated datasets for k-means (blue), hierarchical agglomerative clustering (yellow), and MGLASSO (brown) methods are plotted against the ratio n p . The cluster assignments of MGLASSO are computed from a distance between estimated regression vectors for a given value of λ 2 . Missing boxplots for MGLASSO thus mean computed partitions in the grid of values of λ 2 do not yield the xed number of clusters. The higher the ARI values, the better the estimated clustering partition is.

Figure 3 .

 3 Figure 3.16: Comparisons of adjusted Rand indices for MGLASSO and some clustering methods in the hierarchically structured model. The Rand index is computed between the estimated partition in k clusters and the partition obtained after cutting the tree used for simulation in k clusters. The boxplots of the indices distribution are provided for dierent ratios n/p. The compared approaches are the HAC, k-means, spectral clustering and vanilla convex clustering. Missing boxplots for MGLASSO or convex clustering mean computed partitions do not yield the xed number of clusters.

Figure 3 .

 3 Figure 3.16 compares the performances of MGLASSO, HAC, k-means, spectral clustering, and convex clustering under dierent n/p ratios for dierent numbers of clusters using the hierarchically structured covariance matrix simulation framework.

20

  Figure 3.17: Clustering path of MGLasso solutions on human microbiome data composed of 127 operational taxonomic units. OTUs are colored according to their phylum classication. The path displays abrupt merges. The pure cluster on the graph's left side (down) corresponds to the phylum Bacteroidetes.

  The DNA sequence is an oriented sequence made up of four character states, also known as bases: adenine (A), thymine (T), cytosine (C), and guanine (G). In the paired structure, C is linked to G, and A is linked to T to form a base pair (bp). Some DNA regions (coding regions) can be transcribed into messenger ribonucleic acids (mRNA). The mRNAs can, in turn, be translated into proteins (chains of amino acids) which are involved in various cellular processes the central dogma of molecular biology (see Figure4.1).

Figure 4 . 1 :

 41 Figure 4.1: Central dogma of molecular biology (source: National Human Genome Research Institute).

Figure 4 . 2 :

 42 Figure 4.2: Genes in two samples with dierent numbers of reads (source: (Rau, 2017))

Figure 4 . 3 :

 43 Figure 4.3: Single nucleotide polymorphism (source: David Hall)

  addition of a methyl group to a cytosine position. Methylation occurs predominantly on CpG(CG) dinucleotides. Non-CpG methylation sites include CHG and CHH contexts where H is a DNA base dierent from the guanine.

Figure 4 . 4 :

 44 Figure 4.4: Dierent methylation contexts (image readapted from Colette Picard and Robert Erdmann)

  Shotgun metagenomics and 16S metagenomics are two dierent methods used to study the microbial diversity in a sample. Shotgun metagenomics involves sequencing all the DNA in a sample, regardless of its origin. The DNA is then broken down into smaller fragments, which are sequenced in parallel. The resulting sequences are then assembled into contigs and then into larger sequences called scaolds. Shotgun metagenomics provides a comprehensive overview of the microbial diversity in a sample, including both known and unknown species. It also allows for the identication of functional genes and metabolic pathways present in the sample. 16S metagenomics, on the other hand, focuses specically on the 1 https://www.genome.gov/genetics-glossary/Microbiome

(

  EPITREE, ANR-17-CE32) is a forest research project that focuses on how genetic and epigenetic variations contribute to phenotypic plasticity and adaptation to the local environment. EPITREE was born from a need to understand the mechanisms underlying forest tree adaptation to manage genetic resources better. Indeed, over the last few years, widespread forest die-o has been observed due to drought constraints. These trees play an essential role in providing ecosystem balance on earth. Two tree species are studied in the project, namely the poplar and oak.The project is one of some, interested in the genetic bases of trees' local adaptation. However, the existing studies usually focus on the contribution of SNPs. Epigenetic mechanisms are not studied in depth. The dynamic nature of the epigenome makes it an interesting subject to study in these long-lived organisms. EPITREE focuses on studying the variations of DNA methylation, gene expression, and genome structural variations for a better understanding of the contribution of epigenetic variations in local adaptation and phenotypic plasticity. The project is subsetted in multiple work packages which center around the identication of epigenomic candidate regions, the characterization of epigenomic variation in natural populations and its functional consequences, the characterization of the epigenomic plasticity and its functional consequences in response to environmental constraints, data generation, and integrative multi-omics analysis.Understanding the evolutionary pattern of species is essential in biology, espe-cially for long-lived organisms such as trees. The study of genetic and epigenetic variations of these trees' natural populations can help shed light on their mechanisms of adaptation to their local environment, which is relevant in the current context of climate change. Simultaneous interest in genetic and epigenetic variation is a recently explored area of research focused on annual plants, which are short-lived organisms[START_REF] Mamadou | Narrow-sense heritability and pst estimates of dna methylation in three populus nigra l. populations under contrasting water availability[END_REF]. Characterizing these variations for trees remains an open question that can bring added value in forest management[START_REF] Amaral | Advances and promises of epigenetics for forest trees[END_REF]. One of the objective of the EPITREE project is to provide answers to this concern by producing diversied data on populations of trees located on dierent geographical sites, focusing on oak and poplar.

  ),[START_REF] Mamadou | Epigenetic variation in tree evolution: a case study in black poplar (populus nigra)[END_REF][START_REF] Chateigner | Gene expression predictions and networks in natural populations supports the omnigenic theory[END_REF]. Twenty genotypes corresponding to ten natural populations (two genotypes per population) of poplars (Populus nigra), representative of the diversity of the black poplar collection in Western Europe (see Figure4.5), were investigated.

Figure 4 . 5 :

 45 Figure 4.5: Map of Western Europe illustrating natural populations of black poplars. The ten natural populations are partitioned into six ancestral populations whose colors are shown on the map. Black dots indicate sites where populations of Populus nigra were present, including sites in France, Germany, and Italy.

  are applied to select the cytosines that vary the most between the samples. The histogram of the distribution of standard deviations in the CpG context is given in Figure 4.7. We chose an arbitrary threshold of at least 20. After the ltering step and the selection of complete cases, we end up with 260078 SMPs in the CpG context, 92690 SMPs in the CHG context, and 82436 SMPs in the CHH context.

Figure 4

 4 Figure 4.6: Histogram of minor allele frequencies. The SNPs with a minor allele frequency (MAF) less than 5% are discarded from the study.chr start end strand coverage1 numCs1 numTs1 10 Chr01 387 387 + 13 2 11 14 Chr01 416 416 + 21 2 19 18 Chr01 466 466 + 37 6 31 20 Chr01 504 504 + 44 9 35 21 Chr01 506 506 -25 3 22

Figure 4

 4 Figure 4.7: Histogram of the distribution of standard deviations of methylation levels per genomic position in the CpG context. Positions with a standard deviation greater than 20 are selected. geneid ratio1.Adour.gbM.CG ratio2.Adour.gbM.CG 1 Potri.001G000200 0.80 0.74 2 Potri.001G000300 0.70 0.64 3 Potri.001G000400 0.29 0.33 4 Potri.001G000500 0.01 0.01 5 Potri.001G000600 0.01 0.01

  ' genetic and epigenetic structure The recovered groups from SNP data and the natural groups from the geographical structure are concordant in dierentiating the 10 natural populations of poplar genotypes (see gures 4.8 and 4.5). Indeed, each pair of poplars per geographical site rst merged before joining other sites. One genotype from the Rhine population was discarded due to some ltering constraints on the methylation data and removed from the SNPs and expression datasets. The hierarchical tree in Figure 4.8 summarizes the pattern of evolutionary relatedness among the group of poplars' species. With a partition into 3 clusters, we note that the rst group, composed of the populations of Basento and Paglia, corresponds to geographical sites in the south of the map (Figure 4.5). The second group, composed of Dranse, Kuhkopf, Rhin, and Ticino, corresponds to the geographical sites of the center. The last group, composed of populations from Adour, Val d'Allier, Loire, and Ramières, corresponds to the geographical sites of the West.

Figure 4

 4 Figure 4.8: Genetic structure of natural poplar populations from SNP data. The clustering dendrogram is obtained by hierarchical ascendant clustering based on the similarity matrix of genomic relationships between genotypes. A 3 cluster partition is indicated by dierent colors for the branches and leaves of the tree.

Figure 4

 4 Figure 4.9: Epigenetic structure of natural poplar populations from DNA methylation data in the CpG context. The clustering dendrogram is obtained by hierarchical ascendant clustering based on the similarity matrix of correlation between genotypes.

Figure 4 .

 4 Figure 4.10: Histogram of unadjusted p-values

Figure 4 .

 4 Figure 4.11: Epigenetic structure of natural poplars based on 1500 markers from DNA methylation data in the CpG context. The hierarchical tree on DNA methylation data is plotted against the one computed from the whole SNP data.

Figure 4 .

 4 Figure 4.12: Top 8 enriched gene ontology terms for the biological process ontology in one selected cluster. The bar plot displays the -10 log of the adjusted p-value for the Fisher's exact test.

Figure 4 .

 4 Figure 4.13: Encriched gene ontology with Metascape.

Figure 4 .

 4 Figure 4.14: Methylation proles plotted against expression proles for housekeeping and dierentially methylated genes. The samples are the mean expression and methylation levels of the genotypes belonging to the natural population of Paglia in the CpG context. We do not observe a signicant linear trend between the two genomic entities.

Figure 4 .

 4 Figure 4.16 shows the resulting clustering paths on the entire grid of fusion penalty parameters (Figure 4.16a) and a truncated grid (Figure 4.16b).

Figure 4 .

 4 Figure 4.16: Clustering path of solutions on DNA methylation and transcriptomic samples. The gure shows 3 distinct clusters which correspond to omics data of dierent natures: transcriptomic (right), methylation on the promoter (bottom), and methylation on gene-body (top left).

  the group composed of gene-body proles, we can see a split between the other methylation contexts and the CpG context. These proles in the CpG context can be subdivided into two main groups consisting of the natural populations of the Rhine, Kuhkopf, Paglia, and Ticino on the one hand and the natural populations of Basento, Ramières, Loire Adour, Dranse, and Val d'Allier on the other hand. The expression data suggest three groups of natural populations: Tessin and Ardour; Rhine and Kuhkopf; Val d'Allier, Loire, Ramières, Paglia, Basento, and Dranse.

Figure 4 .

 4 Figure4.17 shows some graphs (adjacency matrices) of the MGLASSO network inference path. As λ 2 increases, the resulting multi-scale graphs are block-diagonal graphs with the number of blocks in {7, 3, 2, 1}.

  Figure4.17: Adjacency matrices for dierent fusion penalty parameters. The rst graph shows the inferred network when no fusion penalty is added to the model. In that graph, the rst block of size 10 × 10 variables corresponds to RNA-Seq samples. The second sparser block of size 30 × 30 corresponds to gene-body DNA methylation data in the three methylation contexts. The last sparse block of the same size corresponds to promoter methylation. The edge bands suggest a relationship between DNA methylation measurements that belong to the same context. For example, the Loire methylation sample in the CpG context is likely related to the Loire samples in the CHG and CHH contexts. The graphs also suggest some relationships between expression and methylation for some natural populations. As the merging penalty increases, the blocks corresponding to the three methylation contexts merge rst, then follow the upper left block corresponding to the expression data. For λ 2 = 30.94, all natural populations merge into a single cluster and complete graph.

Figure 4 .

 4 Figure 4.17d. The interactions highlighted may suggest potential natural populations with a stronger relationship between their transcriptomic and methylation proles.

  body methylation signals for CpG are clearly dierentiated from other methylation contexts in an integrative analysis framework. Note, Sow et al. (2023) indicated that methylation in the CpG context has more adaptive power than methylation in the other contexts.The analysis presented is preliminary; however, it lays the groundwork for further research. It is based on natural populations and not genes and contrasts with integrative analysis approaches performed on genes (Section 5.1 in Rau (2017);[START_REF] Chiquet | A multiattribute gaussian graphical model for inferring multiscale regulatory networks: an application in breast cancer[END_REF]).

  times cannot be dened, one could consider deriving an upper-bound on the fusion penalty parameter and avoid the blind denition of the grid of values. Tan and Witten (2015)'s work on convex clustering and Tibshirani and Taylor (2011)'s work on the generalized LASSO could provide some insights on how the upper-bound can be derived. Overall, incorporating theoretical guarantees into the MGLASSO algorithm could lead to more ecient and accurate clustering results. 5.1.2 Avenues in probabilistic graphical models inference 5.1.2.

  où C désigne l'ensemble des cliques maximales, ϕ = (ϕ C , c ∈ C) est une collection de fonctions potentielles positives, et Z est une constante de normalisation. Les propriétés de factorisation des modèles graphiques permettent d'eectuer des calculs tractables sur des distributions multivariées. En représentant une distribution de probabilité jointe comme un produit de facteurs, où chaque facteur ne dépend que d'un sous-ensemble de variables correspondant à une clique dans le graphe, les calculs peuvent être eectués localement sur les cliques et ensuite combinés à l'aide de la factorisation. Cela permet de développer des algorithmes d'inférence et d'apprentissage ecaces pour des modèles larges et complexes.

  n i=1 log f (x i |Ω) ∝ log det(Ω) -tr(ΩS) (6.1) 121 où S = 1

  ont introduit une approche basée sur la régression pour la sélection de voisinage. Ils utilisent la régularisation Lasso et régressent chaque variable sur les autres, en tirant parti du lien entre les coecients de régression et la matrice de précision. Les diérentes régressions sont ensuite combinées pour inférer le graphe d'indépendance conditionnelle. Plus tard, des auteurs comme Rocha et al. (2008), Ambroise et al. (2009) ont montré que la sélection de voisinage pouvait être considérée comme une approximation de la vraisemblance globale via une pseudo-vraisemblance. Estimation de la vraisemblance globale pénalisée par Lasso D'autres auteurs ont étendu la pénalisation Lasso aux vraisemblances globales du modèle. Contrairement aux approches de sélection de voisinage, qui s'intéressent principalement à l'estimation de la structure du graphe, les approches basées sur la vraisemblance globale permettent d'apprendre la structure du graphe et d'estimer simultanément ses paramètres de manière cohérente. Elles optimisent une fonction convexe en utilisant divers algorithmes.

  et l'Elastic-net Kovács et al. (2021), ont été explorées selon la nature du problème d'inférence.

  applique une combinaison linéaire des pénalités Lasso et Ridge, favorisant un eet de regroupement et permettant la sélection de groupes de variables. OSCAR[START_REF] Howard | Simultaneous regression shrinkage, variable selection, and supervised clustering of predictors with oscar[END_REF] y parvient en mélangeant des pénalités Lasso et ℓ ∞ . Le clustered Lasso[START_REF] She | Sparse regression with exact clustering[END_REF] dénit un critère généralisé de Fused Lasso[START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] où il n'y a pas d'ordre spécique sur les variables. Dans le problème de l'inférence des graphes, pour surmonter ce problème, en plus de la sparsité, plusieurs travaux antérieurs tentent d'estimer la CIG en intégrant des structures de regroupement parmi les variables, dans un souci d'équilibre statistique et d'interprétabilité. Une liste non exhaustive de travaux qui intègrent une structure de clustering pour accélérer ou améliorer la procédure d'estimation comprend Honorio et al. (2009); Ambroise et al. (2009); Mazumder and Hastie (2012a); Tan et al. (2015); Yao and Allen (2019); Devijver and Gallopin (2018).

  La relation entre cette formulation (6.2) et le clustering k-means est détaillée dans Lindsten et al. (2011); et la relation avec la classication ascendante hiérarchique est établie dans Hocking et al. (2011). Le clustering convexe résout le problème de clustering de manière indépendante pour une plage de paramètres de pénalisation λ. Étant donné que le problème d'optimisation est convexe, la solution nale est indépendante de l'initialisation des centroides.

  et similaire au dendrogramme de la classication ascendante hiérarchique, peut être obtenu après une tâche de clustering convexe. Chapitre 3 Dans ce chapitre, nous introduisons la méthode d'inférence de graphe dénommée Lasso Graphique multi-échelle (MGLASSO, Sanou et al.), une nouvelle approche permettant d'estimer simultanément une structure de clustering hiérarchique et des graphes d'indépendance conditionnelle entre des variables, à plusieurs niveaux de granularité.

  de pénalisation LASSO encourage la parcimonie dans les coecients estimés. Le terme fused-group LASSO encourage la fusion dans les coecients de régression β i et β j . Le MGLASSO permet ainsi de mettre en évidence de façon simultanée les interactions entre variables tout en proposant des groupes. Les variables qui appartiennent au même groupe sont suceptibles d'avoir la même structure de voisinage.OptimisationLa fonction objective introduite est la somme de trois composantes convexes : une fonction diérentiable (perte quadratique) et deux fonctions de pénalité non diérentiables : la fonction de pénalité LASSO, qui est séparable car décomposable en composantes indépendantes correspondant chacune aux coecients des prédicteurs, et la pénalité group-fused LASSO, qui n'est pas séparable. Plusieurs approches peuvent être utilisées pour résoudre ce type de problème de minimisation. Parmi celles-ci, se trouvent entre autres la méthode du sous-gradient (Shor, 2012), la méthode ADMM Boyd et al. (2011) ou encore les méthodes de continuation combinées avec des techniques de lissage pour les parties non diérentiables du critère telle que CONESTA (Hadj-Selem et al., 2018). Nous avons comparé les performances de convergence empirique de ces 3 approches d'optimisation. Sauf pour le cas de la méthode du sous-gradient, les algorithmes sont généralement appliqués à des versions reformulées du critère initial MGLASSO, en un problème de régression unique. Avec l'algorithme d'optimisation CONESTA (Hadj-Selem et al., 2018), le MGLASSO atteint une vitesse de convergence plus rapide que les autres approches concurrentes. Une brève analyse empirique de la convergence a été présentée à cette n. En pratique, le MGLASSO est implémenté dans le package R mglasso (version 0.1.2, Sanou (2022)). La pénalité LASSO est gardée xe et le paramètre de pénalité group-fused LASSO varie dans une grille de valeurs. Pour un certain seuil ϵ f use , les variables i et j sont assignées au même groupe si d(i, j) = βiτ ij βj 2 ≤ ϵ f use .

  la méthode MGLASSO, tant en terme de clustering que de recouvrement de support. Les courbes ROC sont utilisées pour évaluer l'adéquation des graphes inférés avec les graphes de référence, pour le MGLASSO et le GLASSO dans sa version 125 de sélection de voisinage, dans les cadres des modèles Erdös-Rényi[START_REF] Erd®s | On the evolution of random graphs[END_REF], scale-free[START_REF] Mark Ej Newman | Random graphs with arbitrary degree distributions and their applications[END_REF], et des modèles blocs stochastiques(SBM, Fienberg and Wasserman (1981)). Les indices de Rand ajustés sont utilisés pour comparer les partitions obtenues avec MGLASSO, la classication hiérarchique ascendante, et la classication k-means pour des modèles bloc stochastique et hiérarchique. Le MGLASSO ne fait pas pire que le GLASSO en terme de performance. En présence de l'ajout d'une pénalité de fusion, l'approche MGLASSO donne lieu à des performances meilleures en terme de courbe ROC. Chapitre 4 Ce chapitre se concentre sur les applications du modèle MGLASSO et d'autres modèles statistiques sur les données omiques. Diverses questions biologiques sont abordées, en particulier dans le cadre du projet EPITREE, qui vise à étudier l'impact évolutif et fonctionnel de la variation épigénétique chez les arbres forestiers. La méthode MGLASSO est illustrée dans une analyse intégrative des données transcriptomiques et de méthylation d'ADN provenant du projet EPITREE. L'impact évolutif et fonctionnel des variations épigénétiques chez les arbres forestiers (EPITREE, ANR-17-CE32) est un projet de recherche forestière qui se concentre sur la façon dont les variations génétiques et épigénétiques contribuent à la plasticité phénotypique et à l'adaptation à l'environnement local. EPITREE est né de la nécessité de comprendre les mécanismes sous-jacents à l'adaptation des arbres forestiers an de mieux gérer les ressources génétiques. En eet, au cours des dernières années, un dépérissement généralisé des forêts a été observé en raison des contraintes liées à la sécheresse. Ces arbres jouent un rôle essentiel dans l'équilibre des écosystèmes de la planète. Le projet s'intéresse aux bases génétiques de l'adaptation locale des arbres. Cependant, les études existantes se concentrent généralement sur la contribution des SNP. Les mécanismes épigénétiques ne sont pas étudiés en profondeur. La nature dynamique de l'épigénome en fait un sujet intéressant à étudier chez ces organismes à longue durée de vie. EPITREE se concentre sur l'étude des variations de la méthylation de l'ADN, de l'expression des gènes et des variations structurelles du génome an de mieux comprendre la contribution des variations épigénétiques à l'adaptation locale et à la plasticité phénotypique. Le projet est subdivisé en plusieurs modules qui s'articulent autour de l'identication de régions épigénomiques candidates, de la caractérisation des variations épigénomiques dans les populations naturelles et de leurs conséquences fonctionnelles, de la caractérisation de la plasticité épigénomique et de ses conséquences fonctionnelles en réponse aux contraintes environnementales, de la génération de données et de l'analyse multiomique intégrative. Comprendre le schéma évolutif des espèces est essentiel en biologie, en particulier pour les organismes à longue durée de vie tels que les arbres. L'étude des variations génétiques et épigénétiques des populations naturelles de ces arbres peut contribuer à mettre en lumière leurs mécanismes d'adaptation à leur environnement local, ce qui est pertinent dans le contexte actuel du changement climatique. L'intérêt simultané pour les variations génétiques et épigénétiques est un domaine de recherche récemment exploré, centré sur les plantes annuelles, qui sont des organismes à courte durée de vie (Sow et al., 2018). La caractérisation de ces variations pour les arbres reste une question ouverte qui peut apporter une valeur ajoutée dans la gestion forestière (Amaral et al., 2020). L'un des objectifs du projet EPITREE est d'apporter des réponses à cette préoccupation en produisant des données diversiées sur des populations d'arbres situées sur diérents sites géographiques, en se concentrant sur le chêne et le peuplier. Dans le cadre de la recherche doctorale, nous avons contribué à mettre en évidence la structure génétique et épigénétique des populations de peupliers grâce aux données SNP et aux données de méthylation de l'ADN, respectivement. Nous montrons également comment la méthylation est utilisée comme marqueur de la diérenciation des populations en proposant des gènes dont le prol de méthylation suit de près la structure génétique des populations. Nous étudions brièvement le lien entre la méthylation et les prols d'expression pour une classe spécique de gènes dont les prols d'expression sont stables entre les arbres appartenant aux mêmes métapopulations. Ces domaines d'analyse choisis résultent de discussions avec des experts du projet EPITREE. Nous avons également exploré brièvement la structure et l'interaction entre les prols transcriptomiques et de méthylation dans diérents contextes de méthylation à l'aide du modèle MGLASSO. between neighbourhood selection and pseudo-likelihood optimizationProof. The pseudo log-likelihood pl of a multivariate normal random vector X ∼ N p (µ, S) is given by

  LASSO and fuse-group-LASSO penalties. Finally, in Chapter 4, the results of applying MGLASSO and other data analysis tools in the context of the EPITREE project are presented. . Sow et al., Epigenetic Variation in Tree Evolution: a case study in black poplar (Populus nigra)Mamadou, bioRxiv 2023.07.16.549253.

	Notations
	E. Sanou, mglasso: Multiscale Graphical Lasso, CRAN package, 2022.
	https://CRAN.R-project.org/package=mglasso In this section we describe the notations used in the remaining of the manuscript.
	Specic sets
	X	discrete or continuous space
	S p ≻0	set of real symmetric p × p positive denite matrices

Publications E. Sanou, C. Ambroise, G. Robin, Inference of Multiscale Gaussian Graphical Model, Computo, 2023.

M

  Undirected Graphical Models . . . . . . . . . 18 2.1.1 Markov properties and factorization . . . . . 19 2.1.2 Gaussian graphical models . . . . . . . . . . . 23 2.1.3 Inference of Gaussian Graphical Models . . . 24 2.2 Convex Clustering . . . . . . . . . . . . . . . . 29 2.2.1 Hierarchical Agglomerative Clustering . . . . 29 2.2.2 K-means clusetring . . . . . . . . . . . . . . . 30 2.2.3 Convex relaxation of k-means and HAC . . . 31 2.2.4 Convex hierarchical clustering . . . . . . . . . 32 2.3 Non-smooth convex optimization . . . . . . . 34 2.3.1 Subgradient methods . . . . . . . . . . . . . . 35 2.3.2 Proximal methods . . . . . . . . . . . . . . . 36 2.3.3 Smoothing methods . . . . . . . . . . . . . . 38

	2 -Mathematical background
	Contents
	2.1
	precision matrix

  Genomics . . . . . . . . . . . . . . . . . . . . 82 4.1.2 Transcriptomics . . . . . . . . . . . . . . . . . 83 4.1.3 Genetics . . . . . . . . . . . . . . . . . . . . . 84 4.1.4 Epigenetics . . . . . . . . . . . . . . . . . . . 84 4.1.5 Microbiomics . . . . . . . . . . . . . . . . . . 85 4.2 The EPITREE Project . . . . . . . . . . . . . 86 4.2.1 Data and material . . . . . . . . . . . . . . . 87 4.3 Dierential analysis and selection of markers 91

	Contents
	4.1 Elements of omics data . . . . . . . . . . . . . 82
	4.1.1 This chapter focuses on applications of the MGLASSO model and other statisti-
	It has the advantage of cal models on omics data. In Section 4.3, various biological questions are answered,
	proposing multi-scale networks with clustering of OTUs, which makes it possible in particular within the framework of the EPITREE project, which aims to study
	to nd interesting links between OTUs that do not necessarily belong to known the evolutionary and functional impact of epigenetic variation in forest trees. Note
	taxonomic classications. The partition into 2 groups found in data from the that we have carried out other analyses within the project that are not reported in
	American gut project suggests a relationship between a group composed exclusively this manuscript. In Section 4.4, MGLASSO is illustrated in an integrative analysis
	of OTUs belonging to the phylum Proteobacteria and the other phyla represented of transcriptomic and DNA methylation data from the EPITREE project. Section
	by the phylum Firmicutes. Setting weights for the fusion penalty can improve 3.5 presents another example of applying MGLASSO to the inference of sparse
	MGLASSO's clustering path with less abrupt fusions. microbial networks.

4.3.1 Dierential analysis . . . . . . . . . . . . . . 91 4.3.2 Clustering of SNP and methylation data . . . 93 4.3.3 Gene sets testing . . . . . . . . . . . . . . . . 94 4.3.4 Results . . . . . . . . . . . . . . . . . . . . . 95 4.3.5 Discussion . . . . . . . . . . . . . . . . . . . . 99 4.4 Integrative analysis of methylation and transcriptomic data through MGLASSO . . . . . 101 4.4.1 Data and pretreatments . . . . . . . . . . . . 4.4.2 Gene selection through sparse PCA . . . . . . 4.4.3 MGLASSO learning . . . . . . . . . . . . . . 4.4.4 Results . . . . . . . . . . . . . . . . . . . . . 4.4.5 Discussion . . . . . . . . . . . . . . . . . . . .

Table 4 .

 4 Datasets are available for the three methylation contexts, with 3.53 million SMPs in the CpG context, 10.41 million in the CHG context, and 53.72 million in the CHH context. The DNA methylation data provided are count data, and a sample is shown in Table 4.2.

1: Sample of SNP data: allele count proles for ve genomic positions over six genotypes. The geographical location of the genotypes is as follows: 1-A26, 1-J31: Ramière; ALL-014, ALL-019: Val d'Allier; AST-005, BDX-003: Adour. DNA methylation levels per cytosine (single methylated polymorphism, SMP) are measured by the whole genome bisulte sequencing (WGBS) method for all 20 samples.

Table 4 .

 4 2: Sample of DNA methylation data in CHG context for ve genomic positions: chr, start, end and strand locate the genomic position while coverage1 , numCs1 and numTs1 are the number of reads, the number of methylated reads and the number of unmethylated reads, respectively for sample 1 (Ramière). proportions of DNA methylation per gene in the gene body and promoter regions are also available. They result from the correspondence between SMPs and genes. A sample of the data is given inTable 4.3. After discarding missing values, 40238, 40318 and 40407 genes methylation levels are available for the CpG, CHG and CHH contexts, respectively.

	The

Table 4 .

 4 3: Sample of DNA methylation data in the CpG context for ve genes (gene body region). The methylation proportions are given for two genotypes located in the Adour area.

(2020)

. A sample of the RNA-seq normalized data is given in Table

4

.4. The transcriptomic dataset used contains 34229 genes.

Table 4 .

 4 5: Contingency table between cluster's genes and ontology term.

  L'inférence de réseaux à partir de données biologiques est utile pour obtenir une compréhension plus complète des mécanismes biologiques sous-jacents à un phénomène particulier. Les réseaux sont un moyen naturel d'intégrer et de décrire la manière dont les variables biologiques interagissent, ce qui permet d'identier des interactions complexes. Avec l'émergence des techniques de séquençage à haut débit, il est possible de générer de grandes quantités de données omiques liées au génome, au transcriptome, aux variations génétiques et au métabolome.

	aux changements environnementaux et le besoin urgent de comprendre comment
	les arbres réagissent au stress environnemental. En eet, le changement clima-
	tique est le facteur qui aura le plus d'impact sur la biodiversité d'ici 2100, après
	l'utilisation des terres (Chapin Iii et al., 2000). On sait que les arbres forestiers
	possèdent des mécanismes complexes qui leur permettent de s'adapter aux facteurs
	de stress environnementaux (Bruce et al., 2007). Il est donc essentiel de compren-
	dre les mécanismes moléculaires qui sous-tendent leur adaptation pour élaborer
	des stratégies de conservation et de gestion des écosystèmes forestiers. Le projet
	EPITREE étudie les mécanismes moléculaires qui sous-tendent l'adaptation des
	arbres, en se concentrant sur deux modèles d'arbres, le peuplier et le chêne. Ces
	espèces ont été choisies pour leur diversité génétique et leur potentiel d'adaptation
	à des environnements changeants.
	La recherche est motivée par le projet EPITREE (Impacts évolutif et fonction-
	nel de variations épigénétiques chez des arbres forestiers), qui cherche à compren-Les modèles graphiques probabilistes (MGP, Lauritzen (1996); Koller and Fried-
	dre comment l'épigénétique, en l'occurrence la méthylation de l'ADN, l'expression man (2009)) sont un outil populaire d'analyse de données en grande dimension et
	des gènes et la variation allélique, inuence les mécanismes d'adaptation et la de capture d'interactions entre des variables. Ils sont largement utilisés dans di-
	plasticité phénotypique chez les arbres forestiers. L'épigénétique est l'étude des verses applications, telles que la génomique et l'analyse d'images, pour réduire le
	changements héréditaires qui aectent l'expression des gènes sans modier l'ADN, nombre de paramètres en sélectionnant les interactions les plus pertinentes entre
	tandis que la plasticité phénotypique est la capacité d'un génotype individuel à les variables. Une classe de MGP particulièrement utile dans un contexte gaussien
	exprimer diérentes valeurs d'un trait phénotypique donné dans des conditions en-est celle des modèles graphiques gaussiens non dirigés. En grande dimension, ces
	vironnementales diérentes (Rey et al., 2016). Les arbres sont des organismes modèles sont souvent supposés parcimonieux. On suppose que seul un petit nom-
	remarquables qui vivent longtemps, ont des cycles de vie complexes et produisent bre de variables interagissent par rapport au nombre total d'interactions possibles.
	du bois, tout en fournissant un large éventail de services écosystémiques. Au cours
	des dernières décennies, un dépérissement généralisé des forêts dû à la sécheresse
	et au stress thermique a été observé dans le monde entier (Anderegg et al., 2016).
	Ces événements mettent en évidence la vulnérabilité des écosystèmes forestiers

La grande dimension de ces données constitue la principale diculté d'un point de vue statistique et d'interprétation. L'objectif de la recherche est de proposer et d'étudier une méthode d'inférence de réseau qui prend en compte une structure de groupe ou une hiérarchie entre les variables. L'identication de groupes de n÷uds densément connectés dans le réseau peut correspondre à des variables biologiques ayant des fonctions apparentées et ore la possibilité de construire des structures multi-échelles pour synthétiser l'information récupérée par les groupes et améliorer l'interprétabilité. Cette tâche de regroupement peut être envisagée avant ou en même temps que la tâche d'inférence du réseau. Les modèles graphiques probabilistes représentent un modèle bien adapté pour déduire les relations entre les variables. La recherche se concentrera sur une sous-classe générale de modèles graphiques où la distribution conditionnelle aux n÷uds est gaussienne et examinera les méthodes d'estimation de groupes basées sur un critère convexe. Le projet se compose de plusieurs modules de travail, qui comprennent entre autres le criblage de régions candidates, l'analyse épigénomique et le séquençage du génome à l'aide de technologies omiques modernes. Ces technologies ont permis de générer des données sur les polymorphismes simples méthylés, les régions diérentiellement méthylées, l'expression des gènes et les polymorphismes de nucléotides simples. La recherche doctorale est principalement motivée par le quatrième module de travail, qui vise à eectuer une analyse intégrative pour modéliser des relations multi-échelles entre les caractères quantitatifs et leurs déterminants moléculaires. Plus précisément, ce module vise à quantier la contribution de la diversité génétique et épigénétique à la variation phénotypique, à étudier l'impact de l'évolution du chêne et du peuplier sur la plasticité épigénomique, et à déterminer si les régions diérentiellement méthylées au niveau des gènes sont conservées entre les deux espèces. En outre, ce module vise à améliorer les modèles de prédiction de la variation des caractères quantitatifs en combinant les informations génétiques et épigénétiques. Les modèles graphiques sont un outil utile pour déduire les interactions entre l'information génétique et les schémas de méthylation, et peuvent apporter des réponses à certaines des questions soulevées dans le projet.

  d'ensembles de gènes, la transformation des données de comptage et les approches de sélection des gènes. Les applications omiques compren-

	nent des applications sur les données transcriptomiques (expression des gènes),
	épigénétiques (notamment la méthylation de l'ADN), génétiques (notamment les
	SNP) et un exemple d'illustration sur des données métagénomiques (abondance
	microbienne).
	La méthodologie proposée pour l'inférence de modèle graphique, dénommée
	LASSO graphique multi-échelle (MGLASSO), est une méthode basée sur la pseudo-
	vraisemblance pour estimer des structures de clustering hiérarchique et des modèles
	graphiques qui dépeignent la structure d'indépendance conditionnelle entre des
	groupes de variables à chaque niveau de la hiérarchie. La méthode MGLASSO
	combine la sélection de voisinage et le clustering via une pénalité de type fused-
	LASSO
	Nos recherches se focalisent sur l'inférence de structures de clustering hiérar-
	chique, plus intuitives pour l'interprétation, et sur l'apprentissage de la structure
	de réseau inférée plutôt que sur l'estimation des coecients de la matrice de pré-
	cision. Bien que certaines des applications soient principalement motivées par des
	questions biologiques du projet EPITREE, répondre à ces questions peut nécessiter
	l'utilisation d'autres outils d'apprentissage automatique dédiés autre que l'inférence
	de graphes. D'un point de vue mathématique, nos recherches se situent à la croisée
	de l'inférence graphique probabiliste, du clustering et de l'optimisation convexe. Du
	point de vue de la biologie statistique, nous abordons diverses questions biologiques,

méthodes discutées précédemment utilisent la structure de groupe pour simplier le problème d'inférence de graphe et inférer le graphe d'indépendance conditionnelle entre des variables non groupées. Cependant, l'inférence du graphe entre des groupes de variables ou des variables représentatives de ces groupes a reçu peu d'attention. Bien que certains travaux aient abordé ce problème, ils se sont surtout concentrés sur des estimations à deux niveaux, c'est-à-dire au niveau des variables prises individuellement et des groupes connus a priori (voir, par exemple, Cheng et al. (

2017

)). Le problème de recherche abordé dans ce travail vise à dénir une méthode d'inférence qui permet des estimations à plus de deux niveaux de granularité avec des groupes inconnus. Ce problème est principalement motivé par des applications en analyse de données biologiques où des données provenant de sources multiples, typiquement des données multi-omiques, sont analysées. Dans de tels cas, il peut être nécessaire de regrouper des variables partageant les mêmes caractéristiques et de prendre simultanément en compte la structure de regroupement dans la procédure d'inférence de réseau, en utilisant une fonction de coût unique au lieu d'alterner les tâches de clustering et d'inférence de réseau. concernant principalement l'impact de l'épigénétique sur l'adaptation locale du peuplier, ce qui nécessite des outils tels que l'analyse diérentielle des gènes, l'analyse d'enrichissement

  La propriété de Markov globale est liée à deux autres propriétés de Markov : celles de Markov locale et de Markov par paires. La propriété de Markov locale stipule que chaque noeud est conditionnellement indépendant des noeuds non-voisins, étant donné ses noeuds voisins. La propriété de Markov par paires quant à elle stipule que deux noeuds non adjacents sont conditionnellement indépendants, étant donné leurs voisins communs. La factorisation, un autre concept important, joue un rôle signicatif dans la compréhension des modèles graphiques. Un graphe est considéré comme complet lorsque chaque arête possible existe entre les noeuds. On note E = P 2 (V ). (Clique). Dans un graphe, un sous-ensemble C est appelé clique si le sous-graphe G C induit par C est complet.

	Denition

Some popular methods based on proximal operators include, among other approaches, the proximal gradient method, the fast iterative shrinkage-thresholding algorithm (FISTA,[START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]), and the alternating direction method of multipliers[START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. In the following, we focus on the specic case of the alternating direction method of multipliers (ADMM).

β (k+1) = β (k) -η k ν (k) ;

(a) Mean ROC curves for MGLASSO and GLASSO graph inference in the Erdos-Rényi model. (b) Mean ROC curves for MGLASSO and GLASSO graph inference in the scale-free model.

dierentially methylated genes in the rst, second, and third metapopulations, respectively.A simple linear regression model is then tted between the expression and

Acknowledgments

Algorithm 4: MGLASSO algorithm Input : X = {X 1 , . . . , X p } ∈ R n×p , λ 1 ≥ 0, λ 2initial > 0, η > 1, ϵ f use ≥ 0

Output: β(λ 1 , λ 2 ) ∈ R p×(p-1) , C(λ 1 , λ 2 ) for all (λ 1 , λ 2 )

The mglasso package

The algorithm 4 is available as an R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] package mglasso (version 0.1.2, [START_REF] Sanou | mglasso: Multiscale Graphical Lasso[END_REF]) available on the comprehensive R archive network (CRAN) and Github (https://github.com/desanou/mglasso). It is based on the Python (vanRossum, 1995) library pylearn-parsimony [START_REF] Hadj-Selem | Continuation of Nesterov's Smoothing for Regression with Structured Sparsity in High-Dimensional Neuroimaging[END_REF] and the R package reticulate [START_REF] Ushey | reticulate: Interface to Python[END_REF] which provides an R interface to Python code. The mglasso package provides both Python and R functions and can be installed at https://cran.r-project.org/package=mglasso. R> install.packages("mglasso") R> library(mglasso)

In the following, we show how to use in practice the principal functions and highlight a hidden function that builds the linear operator A in the reformulated criterion 3.12.

Installing Python libraries from Github into R

After installing mglasso, the required python libraries are delay loaded. That allows the package to be loaded even when the Python engine is not correctly set up.

The next step is then installing the Python dependencies. For libraries not available on Conda distribution or PyPi, reticulate does not oer an easy option to install them directly. The function install_pylearn_parsimony() in mglasso needs to be run in order to install all the libraries not automatically installed by reticulate.

The following code installs the libraries in the Conda environment rmglasso. R> install_pylearn_parsimony(envname = "rmglasso", method = "conda") R> reticulate::use_condaenv("rmglasso", required = TRUE) R> reticulate::py_config() ## Initialize the Python engine

Model performance

In this Section, we conduct a simulation study to evaluate the performance of the MGLASSO method, both in terms of clustering and support recovery. Receiver Operating Characteristic (ROC) curves are used to evaluate the adequacy of the inferred graphs with the ground truth, for the MGLASSO and GLASSO in its neighborhood selection version, in the Erdös-Rényi [START_REF] Erd®s | On the evolution of random graphs[END_REF] , 1981), the Erdös-Renyi model [START_REF] Erd®s | On the evolution of random graphs[END_REF] and the Scale-Free model [START_REF] Mark Ej Newman | Random graphs with arbitrary degree distributions and their applications[END_REF]. In each case, Gaussian data is generated by drawing n independent realizations of a multivariate Gaussian distribution N (0, Σ) where Σ ∈ R p×p and Ω = Σ -1 . The support of Ω, equivalent to the network adjacency matrix, is generated from the three dierent models.

The diculty level of the problem is controlled by varying the ratio n p with p xed at 40: n p ∈ {0.5, 1, 2}.

Stochastic Block-Model

We construct a block-diagonal precision matrix Ω as follows. First, we generate the support of Ω as shown in Figure 3.10, denoted by A ∈ {0, 1} p×p . To do this, the variables are rst partitioned into K = 5 hidden groups, noted C 1 , . . . , C K described by a latent random variable Z i , such that

where α kl is the probability of inter-cluster connectivity, with α kl = 0.01 if k ̸ = l and α ll = 0, 75. For k ∈ {1, . . . , K}, we dene p k = p i=1 1 {Z i =k} . The precision matrix Ω of the graph is then calculated as follows. We dene Ω ij = 0 if Z i ̸ = Z j ; otherwise, we dene Ω ij = A ij ω ij where, for all i ∈ {1, . . . , p} and for all that a shrinkage eect might be observed in the regression vectors' estimates as the λ 2 penalty parameter increases. This shrinkage eect of group-fused penalty terms was also observed in [START_REF] Chu | Adaptive convex clustering of generalized linear models with application in purchase likelihood prediction[END_REF]. Moreover, the MGLASSO may fail to recover the edges in a stochastic block model framework where the inter-clusters edge connection probability is high.

Clustering

In this section we evaluate the performance of MGLASSO and compare with existing methods in the stochastic-block model and the hierarchically structured covariance matrix model.

Dierential analysis 4.3.1.1 Binomial regression model for methylation data

The problem addressed here is the denition of a method to identify positions of the genome for which the methylation proportions are signicantly dierent between sample groups. A treatment variable or clustering usually determines the groups. The issue of dierential analysis for methylation data has been tackled in Robinson et al. (2014b); [START_REF] Wreczycka | Strategies for analyzing bisulte sequencing data[END_REF]. Common methods are based on standard statistical tests, regression models, or even hidden Markov Models.

Here, the focus is put on a regression-based method. These approaches are the most suited when dealing with multiple groups in the samples and allow adding additional covariates in the model if needed.

Denote y j = (y 1j , . . . , y nj ) ⊤ ∈]0, 1[ n the vector encoding the proportion of methylated cytosines at a position j of the DNA, in the n samples. Denote N j ∈ N n the vector of reads coverages. The proportion of methylated cytosines equals the ratio between the number of cytosines and the read coverage. Let X j ∈ {1, . . . , K} n be the categorical n-dimensional covariate with K ≥ 2 categories corresponding to the groups to which the samples belong.

Let's introduce Ỹj , the response variable deduced from y j corresponding to the number of successes in the N j trials. We have Ỹij

genotypes from the same geographical sites. Our results also suggest that DNA methylation can be used as a marker of population dierentiation according to the genetic structure. To this end, we have proposed a list of dierentially methylated genes markers of the genetic structure of the population.

When examining methylation within meta-populations, we did not nd a clear pattern between methylation and expression levels for housekeeping genes. The analysis could be renewed by focusing on a smaller sample of genes and by dening, together with the biologists, much more stringent lters. We refer to Sow et al.

(2023) for more details on the biological ndings of this study. 100

Integrative analysis of methylation and transcriptomic data through MGLASSO

This section presents an application of our MGLASSO approach (Chapter 3), to transcriptomic and DNA methylation data. They were collected as part of the EPITREE project [START_REF] Maury | Epigenetics in forest trees: role in plasticity, adaptation and potential implications for breeding in a context of climate change[END_REF] and have already been presented in the previous section (Section 4.3). We performed an integrative analysis in which we modeled the multi-scale relationships between natural poplar populations for a set of selected genes. Unlike Section 4.3, where clustering is performed on omics data of dierent natures separately, in what follows, we incorporate several types of omic data to generate a clustering path of the natural populations of poplars and multi-scale interaction graphs. We laid the groundwork for further analysis by briey exploring the structure and interplay between transcriptomic and methylation patterns in dierent methylation contexts. The study detailed here is summarized in the application section on real data of the MGLASSO article [START_REF] Sanou | mglasso: Multiscale Graphical Lasso[END_REF]. We present the data in Section 4.4.1. Section 4.4.2 introduces the approach used for gene selection. Section 4.4.3 gives details on how the model is learned. We present some results in Section 4.4.4.

(4.7)

For l = 1, we can write X (1) = ũṽ ⊤ where ũ is a norm-1 n-vector and ṽ is a p-vector. The rank one matrix X ⋆ = ũ⋆ ṽ⋆ is estimated as follows, with parsimony constraints on the loading vector:

where ∥.∥ F is the Frobenius norm and λ ≥ 0 the LASSO penalty parameter. The authors propose an iterative procedure to obtain the rst sparse loading vector ṽ/∥ṽ∥ 2 . The other sparse loading vectors are also computed via a rank one approximation.

In practice, the selection of the parameter λ is solved in [START_REF] Rohart | mixomics: An r package for `omics feature selection and multiple data integration[END_REF] by dening the degree of parsimony as the desired number of signicant genes, i.e., the number of non-zero entries in the vectors loading. We applied a sparse PCA on the omics datasets taken separately; then, on the rst 3 principal components, we selected 15 genes. This yielded a set of 315 candidate genes. Removing duplicates leaves us with a nal set of 151 genes on which we integrate omics data.

MGLASSO learning

We applied the MGLASSO method introduced in Chapter 3, on the material of the study. Natural populations are considered as variables, and genes as observations. The sample size is, therefore, n = 151, and the number of variables p = 70.

The bounds on the penalization parameters introduced in Section 3.3 were written after the study. They have, accordingly, not been tested in practice with MGLASSO applications. Following the discussion on model selection in the section mentioned above, we selected only the LASSO penalty parameter, and this via the StARS approach [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF] to provide maximal network stability. In proceeding so, MGLASSO is equivalent to the neighborhood selection method [START_REF] Meinshausen | High-dimensional graphs and variable selection with the Lasso[END_REF]. The LASSO penalty parameter is chosen from a grid of 50 equidistant values on a logarithmic scale [λ 1min , . . . , λ 1max ] where λ 1min = λ 1max /1000 and λ 1max is the maximum of the smallest values of the LASSO tuning parameters for which the node-wise regression coecients are all zero [START_REF] El Ghaoui | Safe feature elimination in sparse supervised learning technical report no[END_REF]. In practice, we used the huge package [START_REF] Zhao | The huge package for high-dimensional undirected graph estimation in r[END_REF] to compute the StARS selection with a grid of penalization parameters values provided. The StARS cut-point value (variability threshold) was set to 0.05.

We drew 50 subsamples from the data with a subsampling block size set to 0.8n.

The sequence of fusion penalty parameters λ 2 is a grid of 20 equally-spaced values in the interval [0, 30.94].