Introduction

"Mathematics is biology's next microscope, only better; biology is mathematics' next physics, only better" wrote Joel E. Cohen in 2004. The microscope has allowed many breakthroughs in biology. It has revealed a formerly invisible world, the existence of microorganisms such as microbes, and gave birth to cellular biology. Similarly, mathematics can shine a new light on biological problems. Computational methods can test many candidates to find drugs to block a virus such as HIV. Furthermore, mathematics offer tools to visualize and understand data. Conceiving such tools is an important challenge in today's biology, as the quantity of data is booming. Firstly, scientists' collaboration and internet databases have increased the amount of available data. Secondly, the price to acquire new data has dropped in many fields. For example, in 1990 started the Human Genome Project, which successfully sequenced a human genome after 13 years of work, mobilizing a hundred laboratories and costing three billion dollars. Today, companies consider sequencing someone's genome for less than a thousand dollars, and within thirty hours. Molecular biology is switching from an era of data scarcity to an era of data abundance. One hopes that mathematics will reveal knowledge from large data, as the microscope revealed images of micro-organisms.

For mathematicians, biology raises many challenges. The dynamics of living systems are caused by interactions of many actors, on various scales. In most biological cases, mathematicians have very few data compared to the problem complexity, not to mention the noise surrounding those data. However, challenges in the empirical sciences have fostered discoveries in mathematics, in particular physics. As Joseph Fourier said, "L'étude approfondie de la nature est la source la plus féconde des découvertes mathématiques" 1 . His study of heat propagation lead to modeling by partial derivative equations, Fourier analysis, Fourier series, etc. Similarly, the theory of distributions was inspired by challenges in physics.

The synergy between mathematics, computer science and biology has been highly visible in the last decades. Biological problems have found useful available tools in the two other fields. Yet, there is no doubt that collaboration can go much further, and can increase our understanding of biological mechanisms. In particular, much bioinformatic research is devoted to identify gene regulatory network (GRN). A GRN is a collection of DNA segments that interact indirectly with each other, interact with other substances in the cell and respond to the external environment. For example, in the presence of sugar, a yeast cell will turn on genes to process the sugar to alcohol. Yeast's GRN commanded this process, which is necessary to yeast's living as it made the yeast cell gain energy to multiply -incidentally, this process is necessary to men for wine-making. With the multiplication of data -in particular gene expression data-, scientists hope to discover the GRN's mechanisms. This knowledge would increase our understanding of the living system, and may help creating cure against specific pathology. Considering a living system's malfunction, one could see which genes are involved in that malfunction, then target some genes in the functional pathway to make the living system healthy.

Bioinformaticians have tackled the GRN inference (GRNI) problem, with the relevant mathematical tools. A common approach to GRNI in bioinformatic is reverse-engineering. Given gene expression data, bioinformaticians build a model that mimic the observed genes' behavior. Then, the bioinformaticians know what regulatory interactions happen in his model, and assume that the same interactions happen in the living system. Mathematical tools exist to model an observed system. In particular, many tools are provided by machine learning, a field of computer science and mathematics that sprang up with computational capacities and its applicability to many fields (image recognition, text categorization, spam detection). Nevertheless, gene modeling is a very demanding task: observed data contain much noise, gene have a nonlinear behavior, some actors of the GRN are not observed, the scale of the data demands computationally fast methods, data types are heterogeneous and, compared to the problem complexity, few data are available.

Among the relevant tools, kernel functions provide methods that are robust-tonoise, able to model any nonlinear behavior and computationally fast. Besides, kernel-based models are the result of minimizing a loss function. By adding constraints or modifying the loss function, one can incorporate prior knowledge or external sources of information to reduce the problem's complexity. Nevertheless, kernel methods have scarcely been used for GRN Inference. The reason was mainly their lack of interpretability. In this work, I present two contributions to network inference using interpretable kernel methods. The two approaches originate from the same idea: a model is interpretable if the importance of input features for the output prediction can be weighted. In the first approach, I assume that partial derivatives of a perfect gene model should reflect the GRN. Indeed, if the concentration of a regulator gene changed, it should affect the concentration of the regulated gene. I demonstrate that partial derivatives of kernel-based models can consistently estimate the mean of partial derivatives of the ideal model. Thus I will interpret feature importance in a kernel-based model by its partial derivatives. In the second approach, I turn to multiple kernel-based models with multiple local kernels, each kernel being devoted to one feature. The idea is to find the optimal linear combination of these local kernels for the modeling of a target gene. The linear combination will weight the importance of each local kernel, thus of each feature. To get more stability and to tackle high dimension data in both approaches, ensemble of those models are built using a double scheme of randomization which provides a drastic improvement in terms of performance. Besides, this randomization scheme allow learning from heterogeneous data types. These two ways to interpret kernel models are then used to infer GRNs. On real and realistically simulated datasets, these methods show state-of-the-art performances: on some well-referenced datasets, they perform better than current state-of-the-art methods. They do under-perform on some other datasets. However, combining a kernel GRN inference method with other state-of-the-art GRN inference methods lead to substantial improvement over state-of-the-art. This thesis is organized as follows: in Chapter 1, I introduce informally the GRN Inference problem, as well as the relevant ideas of machine learning and kernel methods, and describe current GRN inference methods. In the Chapter 2, I demonstrate that kernel methods consistently estimate any continuous linear form of the partial derivatives. In Chapter 3, I give several methods to interpret a kernel model, and observe that, on many realistically simulated data, using partial derivatives is the most efficient method for network inference. Then, I use this method on real and real-sized networks, showing complementary results to other GRNI methods. In Chapter 4, I describe prevalent kernel feature selection methods, and develop two methods: one based on multiple kernel learning, another based on kernel alignment. On realistically simulated data, these two methods yield better results than other kernel feature selection methods. In Chapter 5, one of these two feature selection methods is used on real data, and show state-of-the-art performances. I suggest a modification of this method, to take into account prior information and specificities of biological data. Incorporation of reasonable prior knowledge greatly enhances the performances of this method.

Chapter 1 Context 1.1 Introduction

In the twentieth century, many breakthroughs have profoundly transformed the field of biology. Firstly, as early as 1930, new knowledge and technologies have allowed scientists to analyze living systems and phenomena on a molecular scale; most notably, the discovery of DNA. Secondly, starting in the nineties, new tools have been developed to measure gene or protein expression levels, giving scientists another type of data to analyze and understand cells and living systems. The cost of these tools has greatly decreased over the years. Also, scientific collaboration through data publication-in journals or in databases-greatly increased. As a result, the amount of available data augmented, creating new possibilities for scientific research. Along with experimental design, i.e. choosing which experiments to perform, the study of large data is a considerable step toward new knowledge. Mathematics and computer science provide powerful tools for the exploration of large datasets. The inter-disciplinary field of bioinformatics has developed in this context. This is the study of mathematical and computational tools for the analysis of biological systems. Bioinformatics includes storing and visualizing biological data, and assessing relationships between phenomena. In particular, many bioinformatics methods have been developed for the problem of gene regulatory network inference (GRN inference, or GRNI). This thesis is devoted to developing new tools for GRNI. This chapter will explain the matter and current approaches. In the first section, I introduce informally the gene regulatory network inference problem. Secondly, I describe and discuss available data and the technology to produce these data. Thirdly, I introduce a central tool in bioinformatics: machine learning, and, in particular, kernel methods. Fourthly, I present current methodologies and discuss their most efficient version. In the fifth section, I mention key papers comparing GRNI methods, and give their conclusions.

Gene Regulatory Networks -an introduction

All of our cells contain the same genetic information, contained in our DNA. Nevertheless, skin cells are different from liver or kidney cells. These differences come about because different genes are expressed at high levels in different tissues. So, how are genes "expressed" ? The "central dogma of molecular biology" asserts that "DNA make RNA make protein" [START_REF] Leavitt | Deciphering the Genetic Code: Marshall Nirenberg[END_REF], as illustrated in Figure 1.1.

How are genes expressed ?

DNA consists of sequences of nucleobases A, T, G and C1 . These sequences hold all necessary information for the development and functioning of a living system. In particular, these sequences hold the system's genes. Nowadays, many definitions exist for a "gene". The following definition, from [START_REF] Pearson | Genetics: what is a gene?[END_REF], is sufficient for this thesis: "a gene is a locatable region of genomic sequence, corresponding to a unit of inheritance, which is associated with transcribed regions, regulatory regions, and or other functional sequence regions", with the following explanations:

• a transcribed region is a sequence of A, T, G and C that will be "copied" (transcribed) into messenger RNA. Most messenger RNAs are then "read" by ribosomes, which translates such a sequence into a sequence of amino acids, which forms the protein carrying out the function coded in the gene.

• a regulatory region is a segment of DNA capable of increasing or decreasing gene expression. A particular regulatory region, shared by almost all genes, is known as the promoter, which provides a position that is recognized by the transcription machinery when a gene is about to be transcribed and expressed. A gene can have more than one promoter, resulting in different RNAs.

Figure 1.1: Molecular biology's central dogma Only 1.5% of the human genome consists of protein-coding nucleobases. Some of the noncoding DNA is transcribed in microRNAs, which regulate gene expression, but the purpose of most of the noncoding DNA is yet unknown.

How are genes differentially expressed ?

Some of the proteins hold a gene regulatory function. They are called "Transcription Factors" (sometimes called sequence-specific DNA-binding factors). These proteins bind to specific DNA sequences, controlling gene expression, by promoting (activator) or blocking (repressor) the gene's transcription into mRNA. Thus, expression of specialized genes can regulate other genes, leading to a specialized cell.

How is gene expression measured ?

New technologies allow us to learn many aspects of the genome. In particular, gene expression can be evaluated through RNA sequencing, which gives a snapshot of RNA's presence and quantity in a given tissue at a given time. With mRNA concentration, one has a measure of gene expression -even though other post transcriptional gene regulation events exist, such as RNA interference. The main tool for RNA concentration measure consists in DNA microarrays. To measure an mRNA concentration, one uses a probe, which is a complementary sequence of a part of this mRNA. Each spot of a DNA microarray contains millions of copies of a probe. After extraction from a cell, mRNAs are spread over the array, where they bind to their specific complementary sequence. The array is then washed to remove unbound sequences. Then, the array is scanned with a laser. Each probe produces a fluorescent signal, whose intensity is linked to the number of bound mRNAs. This method faces several noise sources, and noise reducing methods have been developed. For example, probes can "cross-hybridize", i.e. bind with the wrong target. To control cross-hybridization, some arrays pair probes that should work (Perfect Match, PM) with probes that should not (Mismatch MM). PM is perfectly complementary to the sequence of interest, and MM is the same as PM for all but one base. Further discussion of DNA microarrays gene expression measure problematics can be found in [START_REF] Baggerly | An introduction to highthroughput bioinformatics data[END_REF]. The reader must be aware of three main sources of noise in gene expression data. First, there is a noise inherent to observing a living system, especially on these scales. Second, probes have different binding affinities to their target mRNA. Thus, it is difficult to tell whether "gene A beats gene B in experiment 1", as opposed to "there is more gene A in experiment 1 than in experiment 2". Microarrays only produce relative measurements of gene expression. Thirdly, the number of genes is well above the number of available experiments. This is known as the "large p small N" framework, where p is the problem dimension (here, the number of genes) and N is the number of data samples (here, the number of available experiments). This is a very challenging framework, with uncertainties on the results.

GRN Inference

The gene regulatory network represents gene interactions at the transcription level, i.e. which gene regulates which gene. Knowledge of a living system's GRN has many potential applications. It would enhance our understanding of the system. As a possible consequence, it could help the development of cures. Considering a pathology resulting from a system's malfunction, one could see which genes are involved in that malfunction, then target some genes in the functional pathway to enhance or inhibit this function, and hopefully cure the pathology. This explains why GRN inference has become a major challenge in biology. GRN can be described by a graph whose nodes are genes and in which a directed edge from node i to node j means that gene i regulates gene j. This is a simplified view, which does not take into account several key players such as microRNAs. In fact, a regulatory interaction involves DNA, mRNAs and proteins. All these elements are merged into one element, the gene representative, see Figure 1.2. Let us call A the adjacency matrix of this graph : a i j = 1 if j regulates i, 0 otherwise. Network inference usually refers to the estimation of this matrix A.

Available data

A p-gene system is observed through mRNA expression levels in steady-state and time-series data. Steady-state data consist in the concentration of each of the p genes in a particular cell. Noting x k the vector containing the p concentrations in experiment k, x i k is the concentration of the i th gene in the k th experiment. The different types of measures that may be available are: steady-state measures on unperturbed individuals, called wild-type data; mRNA concentration of perturbed individuals, on which a gene's mRNA concentration has been increased or diminished; knock-out experiments, in which a gene is knocked out, hence its expression is null. These data are called "perturbed" or "perturbational data". Time-series, measures of genes' mRNA concentrations through time, are also available. The vector x(t, u) groups concentrations of all genes at time t in experiment u. Usually, the individual has received an exogenous perturbation or signal at time t = 0, such as a heat step, presence of a molecule like glucose, or exposure to radiations. A finite number of observations is available; the system is observed at times t 1 , . . . , t k .

Machine Learning

Machine learning aims at extracting information or knowledge from data. GRN inference takes special interest in supervised learning, a branch of machine learning identifying the link between input variables X and a response or output variable Y. Supervised learning is useful in many domains (computer vision, medical imaging, bioinformatics, etc.), and has fostered the development of many performant algorithms. In particular, Gaussian kernel methods have valuable properties in theory (consistency) and in practice (efficiency, robustness-to-noise), thus are good candidates for GRN Inference.

A simplified view

Adjacency matrix

A =           0 0 0 1 0 0 1 1 0           Figure 1.2:
Regulatory interactions are simplified into a network representing only genes. Image modified from [START_REF] Huynh-Thu | Machine learning-based feature ranking: Statistical interpretation and gene network inference[END_REF] In this section, I first present generalities of supervised learning: notations, goals and difficulties. Next, I describe linear methods, and the associated solutions to supervised learning problems. Finally, I present kernel methods, which can be seen as an extension of linear methods.

Generalities on supervised Learning

The objective of supervised learning is to identify the link between some input variables x = x 1 , . . . , x p ′ and an output variable, or response y. The random variables (X, Y) ∈ (X, Y) are assumed to be distributed according to a distribution P. In all the problems of this thesis, X is a finite p-dimensional space, Y is either a subset of R (regression problem), the discrete set {-1, 1} (binary classification) or a discrete set {1, . . . , m} (multi-class classification). A bold upper case letter, e.g X, denotes a random variable, a lower case letter, such as y or x i , a scalar, and a bold lower case denotes a column vector, e.g. x = (x 1 , . . . , x p ) T . Variables (X, Y) are linked through a function f

Y ∼ f (X) + ǫ (1.1)
with ǫ a zero-mean noise. N realizations of (X, Y) are observed and usually assumed independent and identically distributed (i.i.d.). They form the learning set S = (x i , y i ) i=1...N . Given a prediction function g ∈ F (X, Y), let l : Y × Y → R + be a loss function, quantifying the cost for predicting g(x i ) instead of y i (see examples of classic loss functions in Figure 1.3). The risk of a function g, R(g), is defined as the expected loss:

R(g) = E P (l (g(X), Y)) (1.2)
The goal is to find function f * , from H the set of admissible functions, that minimizes the risk:

f * = arg min g∈H R(g) (1.3)
Unfortunately, the function R(.) is not known. However, the training set S is available so that the empirical risk, i.e. the mean of the loss function of the

Binary classification

Hinge loss square loss

l(g(x), y) = max(0, 1 -g(x)y) l(g(x), y) = (g(x) -y) 2 -2 -1 0 1 2 0 1 2 3 g(x)y l(g(x),y) -2 -1 0 1 2 0 1 2 3 4 g(x)-y l(g(x),y)

Regression

Epsilon-insensitive loss square loss l(g(x), y) = max(0, |g(x) -y|ǫ), ǫ = 0.5 l(g(x), y) = (g(x)y) 2

R emp (g) = 1 N N i=1 l(g(x i ), y i ) (1.4) f = arg min g∈H R emp (g) (1.5)

Bias-variance tradeoff

The choice of the functional space H is critical for the performance of the empirical risk minimization approach. If one chooses too "small" a functional space, there could be no function in H that approximate correctly the true function f and the true risk will be high. The model has bias.

On the other hand, the empirical risk minimizer f depends of the learning set S.

If the functional space is too "big", data may be insufficient, and the learned function f may be very different of the true risk minimizer f * . Those two sources of error can be decomposed this way:

Model error = R( f ) -R( f ) , (1.6) = R( f ) -R( f * ) (1) + R( f * ) -R( f ) (2)
.

(

(1) is the error from learning from S and not infinite data. ( 2) is the error from choosing the functional space H. To find the best model, one has to find the right balance between the bias and the variance of the model. This is quantified by the biais-variance tradeoff. Considering that the loss function is the square loss, the expected true risk of f according to the distribution of the learning set S:

E S (R( f )) = E S E P y -f (x) 2 (1.8)
T.Hastie and R.Tibshirani [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF] prove that this can be decomposed in the following manner-proof in appendix, page 125:

E S (R( f )) = E P,S [ǫ 2 ] + E P ( f (x) -f * (x)) 2 + E S,P f * (x) -f (x) 2 (1.9) = σ 2 + (bias) 2 + (variance) (1.10)
with σ 2 be the variance of the noise,

σ 2 = E P [ǫ 2 ]
. This exposes all the sources of error in supervised learning. σ 2 is the best error one can have, when trying to predict y from x. The bias, (bias 2 , expresses the error from choosing functional space H. This term can be high if H is poorly chosen. Finally, if H is too big, the learned function f will greatly vary with the learning set S, resulting in a third source of error.

) 2 = E P ( f (x) -f * (x))

Curse of dimensionality

The expression "curse of dimensionality", coined by Richard [START_REF] Bellman | Adaptive control processes: a guided tour[END_REF], illustrates the problem of the dramatic increase of the volume of data with the increase of dimension. I give an illustrative example in Figure 1.4, where I have drawn 64 points with coordinates in the interval [0, 1]. In one dimension, the whole space is well occupied. In two dimensions, the data become more sparse.

In three dimensions, the space [0, 1] 3 is clearly under-sampled. To have, in the 10-dimensional space [0, 1] 10 , the same space coverage as for 100 points in a 1dimensional space, we would need 10 20 points [START_REF] Bellman | Adaptive control processes: a guided tour[END_REF].

In this context, it becomes clear that one cannot identify the best prediction function in a "big" functional space H by minimizing only the empirical risk.

Overfitting

Related to the "curse of dimensionality" and the bias-variance tradeoff is the danger of overfitting. Figure 1.5 shows a toy classification problem. Noting U(X) the uniform distribution in X, the test and training data have been drawn according to the following distributions:

X ∼ U([0, 1] 2 ) (1.11) Y ∼ sign(-3 + x 2 + 3x 1 + 2 × U([0, 1])) .
(1.12)

The linear model makes mistakes on the training set, but captures the true classification solution. The nonlinear classification function has perfect prediction on the training set, but would give worse results on the test set.To avoid overfitting, an approach is to regularize the prediction function. One minimizes the empirical risk under a constraint on the complexity of the function: with Ω : H → R + a convex penalty on the complexity of the function. This can be equivalently formulated as min g∈H R emp (g) + λΩ(g) .

(

Even with linear functions, one faces the curse of dimensionality as well as the possibility of overfitting if the vector x is high-dimensional, which is typically the case in bioinformatics. I will now show efficient regularized methods.

Linear model

The first functional space to be studied is the space of linear functions:

f lin (x) = p m=1 a m x m (1.16) = < a, x > (1.17)
The choice of the loss function l(., .) and of the complexity penalty Ω(.) will define various algorithms. For regression problems, the most common methods are ridge-regression (or Tikhonov regression) and LASSO. Each method insures smoothness in a different way.

Ridge-regression

In the case of regression, a popular loss function is the quadratic loss l quad (y, f lin (x)) = (yf lin (x)) 2 .

(1.18)

When adding the ℓ 2 -norm as regularization term Ω(.), the problem is known as ridge regression: min

a∈R p i (y i -< a, x i >) 2 + λ a 2 2 . (1.19)
It admits a closed-form solution:

â = X T X + λI p -1 (X T )y , (1.20)
where I p is the identity matrix of order p. The matrix X T X + λI p is always invertible; the problem is well-posed. As can be seen in equation (1.21),

| f lin (x 1 ) -f lin (x 2 )| = |< a, x 1 -x 2 >| ≤ a 2 x 1 -x 2 2 , (1.21)
the norm of a bounds the ratio between the distance of two input points (x 1 , x 2 ) and their image through the function f lin , ( f lin (x 1 ), f lin (x 2 )). The function f lin is smooth in a Lipschitzian way.

LASSO

By choosing the ℓ 1 -norm as a penalty term, the problem is known as LASSO ("Least Absolute Shrinkage and Selection Operator"): min

a∈R p i (y i -< a, x i >) 2 + λ a 1 . (1.22)
With this regularization term, the solution a is sparse, i.e. many of its entries will be zeroes. along them. At their center is a * the minimizer of the quadratic loss function on the learning set. The regularization can be seen as a second term to minimize, but also as a constraint. In the figure, it is interpreted as a constraint. One optimizes a to minimize the empirical risk given that a's ℓ • -norm is below a threshold T . The optimum, noted a ℓ • , will be found at the first intersection between a contour line and the set of admissible values. With the ℓ 1 -norm, the intersection will be found at a vertex of the set of admissible points. The vertices of admissible points with ℓ 1 -norm are points with some zero components. The property of giving sparse solutions makes the LASSO model interpretable. The few non-zero entries are the only relevant ones. This interpretability allowed LASSO to be used as a feature selection method and also to become a popular regression technique. Unfortunately, many regression problems are not solved with sufficient accuracy by linear methods. I present in the following section an extension of those methods to nonlinear models.

Kernel Methods

To extend the linear model methods to non-linear separation, one approach is to map the data in a high dimensional space. In the example in Figure 1.7, taken

: (x 1 , x 2 ) → (z 1 , z 2 ) = (x 2 1 , x 2
2 ), a perfect linear classifier is available.

Clearly better classifiers can be built by mapping data in the proper space, but selecting the map φ for each problem would be a very hard task. This difficulty has been overcome by using kernel functions and the so-called kernel trick, explained below, without calculating the map or the higher dimensional space. A function k : X × X → R is called a kernel if it has the following properties:

• symmetric: for all (x, x ′ ) ∈ X 2 , k(x, x ′ ) = k(x ′ , x) • positive semi-definite: for all N ∈ N, for all (x i ) i=1...N ∈ X N , for all (α i ) i=1...N ∈ R N , N i, j=1 α i α j k(x i , x j ) ≥ 0
Aronszjan showed, in [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], that k(x, x ′ ) is a scalar product in a particular functional space H ⊂ F (X, R). In particular, there exists a mapping φ such as

k(x, x ′ ) =< φ(x), φ(x ′ ) > H .
The functional space H has the property that, for any f ∈ H, for any x ∈ X, < f, k(x, .) > H = f (x); H is called a Reproducing Kernel Hilbert Space (RKHS). The best function in functional space H is searched:

f = arg min g∈H R emp (g) + λΩ(g) (1.23) f (x) = < â, φ(x) > H (1.24)
Using as regularization term an increasing function of g H , the form of f can be deduced from the representer theorem. Theorem 1. Representer theorem. Let X be a nonempty set and k a positivedefinite real-valued kernel on X × X with corresponding reproducing kernel Hilbert space H. Given a training sample (x 1 , y 1 ), . . . , (x N , y N ) ∈ (X × R) N , a strictly monotonically increasing real-valued function Ω : [0, ∞) → R, and an arbitrary empirical risk function L :

(R × R) N → R ∪ {∞}, then for any f * ∈ H satisfying f = arg min f ∈H L ((y 1 , f (x 1 )), . . . , (y N , f (x N ))) + Ω( f ) (1.25)
f admits a representation of the form:

f (.) = N i=1 α i k(., x i ) (1.26) with α i ∈ R for all 1 ≤ i ≤ N
This is the general formulation of the representer theorem, as given in [START_REF] Schölkopf | A generalized representer theorem[END_REF]. Proof is given in appendix, page 126.

The kernel trick [START_REF] Aizerman | Theoretical foundations of the potential function method in pattern recognition learning[END_REF] consists in calculating all scalar products through the kernel k, and never computing the feature map φ

f (x) = N i=1 α i k(x i , x) (1.27) f 2 H = < f, f > H (1.28) = N i, j=1 α i α j < φ(x i ), φ(x j ) > H (1.29) = N i, j=1 α i α j k(x i , x j ) (1.30)
With a convex loss l and a strictly increasing convex function Ω, solving the problem min

f ∈H i l(y i , f (x i )) + λΩ( f H ) (1.31)
reduces to a convex problem in α, which is quickly solved by gradient descent. In the following work, I will mainly use the kernel-ridge regression.

Kernel ridge-regression

We now have tools to learn non-linear functions with regularization enforcing smoothness. I will use them with the quadratic loss and a regularization term

min f ∈H N i=1 (y i -f (x i )) 2 + λ f 2 H . (1.32)
Noting K the Gram matrix (K i j = k(x i , x j )) and using the representer theorem, this amounts to solving:

min α∈R N (y -Kα) T (y -Kα) + λα T Kα , (1.33)
α admits the closed form solution α = (K + λI N ) -1 y .

(1.34)

Note that minimizing the norm of function f give smoother functions in a Lipschtzian way, similarly to the ridge-regression case

| f (x) -f (x ′ )| = | < f, φ(x) -φ(x ′ ) > H | (H RKHS) (1.35) ≤ f H × φ(x) -φ(x ′ ) H (Cauchy-Schwarz) . (1.36)

GRN Inference from gene expression data

Mathematical tools and, in particular, machine learning tools have been used to infer a GRN from gene expression data, without knowledge of existing interactions. Although they may use supervised learning tools, they are unsupervised methods as they try to predict the existence or absence of edges, but they were not trained with examples of existing or non-existing edges. They can be decomposed into two groups: scoring methods, where one evaluates the dependency of one gene to another through a pre-defined metric, and modeling methods, or reverseengineering methods, where one creates a model f with parameters θ such as this model mimics observed expression data. The model is built so f or θ can be interpreted to infer a gene regulatory network. I give below a list of existing methods with short technical details.

Scoring methods

Statistic dependency measures are used in these methods to evaluate the likeliness of a link between gene i and gene j. Adaptations of these measures have also been suggested to be more relevant in the biological context, or to remove spurious interactions.

Correlation

Correlation studies the intensity of a link between two variables. Many forms of correlation measures exist (Spearman's rank correlation [START_REF] Spearman | The proof and measurement of association between two things[END_REF], Kendall's τ [START_REF] Kendall | A new measure of rank correlation[END_REF], Goodman and Kruskal's γ [START_REF] Goodman | Measures of association for cross classifications*[END_REF]); the following methods use Pearson's correlation ρ. With two variables x and y, ρ's value is contained in [-1, 1], and evaluate the linear dependency of those two variables. With x (resp. ȳ) the mean of x (resp. y), (x i , y i ) observed examples of the couple (x, y), ρ is the cosine between the two centered vectors x and y and is defined as:

ρ(x, y) = i (x i -x)(y i -ȳ) i (x i -x) 2 i (y i -ȳ) 2
(1.37)

Figure 1.8 shows some examples, including ρ = 0 for uncorrelated variables and ρ = -1 for y = ax + b if a < 0. If gene j excites (resp. inhibits) gene i, i will have low (resp. high) concentration when j has low concentration, and they will have a highly positive (resp. negative) correlation. If the correlation is superior to a threshold, a link between these two genes is inferred [START_REF] Butte | Relevance networks: A first step toward finding genetic regulatory networks within microarray data[END_REF].

For more plausible biological networks, some researchers have proposed to postprocess the obtained adjacency matrix to produce scale-free networks2 , see the Symmetric-N [START_REF]Extreme self-organization in networks constructed from gene expression data[END_REF] or the Asymmetric-N algorithm [START_REF] Chen | Rank-based edge reconstruction for scale-free genetic regulatory networks[END_REF]. Kuffner et al. used a non-linear correlation coefficient, called η 2 [START_REF] Cohen | Eta-squared and partial eta-squared in fixed factor anova designs[END_REF], for GRN inference. They improved this measure in [START_REF] Küffner | Inferring gene regulatory networks by anova[END_REF] by giving additional weights to perturbed data.

Mutual information

Another metric is the mutual information I between a pair of random variables

I(x i , x j ) = H(x i ) -H(x i |x j ) = x i ,x j p(x i , x j ) log p(x i , x j ) p(x i )p(x j )
Where H(x i ) is the entropy of random variable x i , measuring its unpredictability, and H(x i |x j ) is the entropy of x i once x j is observed, hence is always lower than H(x i ). If gene j regulates gene i, its concentration will be responsible of much of gene i's concentration, thus x i should be much more predictable once x j is observed, and the mutual information will be high. See the work on Relevance Network [START_REF] Butte | Mutual information relevance networks: functional genomic clustering using pairwise entropy measurements[END_REF] for a first use of mutual information. Many network inference algorithm use mutual information, such as ARACNE [START_REF] Margolin | Aracne: An algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context[END_REF]. Faith et al. suggest the CLR algorithm [START_REF] Faith | Large-scale mapping and validation of escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF]. The mutual information between all gene pair (i, j) is computed. Then, authors can compute the p i distribution of mutual information with gene i. For a randomly chosen gene z, the mutual information between gene i and z follows the distribution p i , MI(i, z) ∼ p i . For a gene j, Faith et al. compute Z i j (resp. Z ji ), the unlikeliness of observing MI(i, j) given p i (resp p j ). They score the likelihood of an interaction with the following equation: .38) An edge between genes i and j is inferred if s(i, j) is superior to a threshold.

s(i, j) = Z 2 i j + Z 2 ji . ( 1

Covariance

For the same reason as for correlation, genes interacting together should have relatively high absolute value for their covariance. Opgen-Rhein and Strimmer [START_REF] Opgen-Rhein | Accurate ranking of differentially expressed genes by a distribution-free shrinkage approach[END_REF] suggest an estimation of covariance through an efficient multi-dimensional estimator: the James-Stein estimator [START_REF] Stein | Inadmissibility of the usual estimator for the mean of a multivariate normal distribution[END_REF]. Gaussian graphical models (GGM), also called concentration graph, covariance selection or Markov random field model [START_REF] Whittaker | Graphical models in applied multivariate statistics[END_REF][START_REF] Dempster | Covariance selection[END_REF][START_REF] Kollar | Probabilistic graphical models: principles and techniques[END_REF], estimate the concentration matrix Ω = Σ -1 , i.e. the inverse of the covariance matrix Σ, assuming all gene expression levels are distributed according to a multivariate normal distribution. The partial correlation between gene i and j knowing all other genes is related to the concentration matrix according to:

ρ(i, j|{1, . . . , p}\{i, j}) = -Ω i j Ω ii Ω j j (1.39)
In particular, Ω i j should be equal to 0 if two genes do not interact, even if a third gene z interacts with both of them. Various methods exist to evaluate Ω or to test if the observed value of Ω i j is sufficiently high to infer an edge between gene i and j. See [START_REF] Johnstone | Needles and straw in haystacks: Empirical bayes estimates of possibly sparse sequences[END_REF][START_REF] Wille | Sparse graphical gaussian modeling of the isoprenoid gene network in arabidopsis thaliana[END_REF][START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF][START_REF] Li | Gradient directed regularization for sparse gaussian concentration graphs, with applications to inference of genetic networks[END_REF].

Z-score

Static, dynamic and perturbed data are available. In particular, an experiment with gene j perturbed should highlight a change of behavior in every gene that gene j regulates. Pinna et al. [START_REF] Pinna | From knockouts to networks: Establishing direct cause-effect relationships through graph analysis[END_REF] measure the Z-scores: for each gene i, they calculate from wild type experiments its mean µ i and standard deviation σ i . Then the Zscore is computed as:

Z i j = G j i -µ i σ i (1.40)
with G j i the concentration of gene i in the experiments where gene j has been knocked out. Z i j measures how "unlikely" the state G j i would be reached by chance if j did not regulate i. The authors also used a method to remove indirect edges. This last approach performed very well on synthetic data, where perturbed data is available for all genes.

Modeling

Another approach to gene regulatory network inference is through modeling of the genes. The parameter θ is sought such that:

x i = f i stat (x -i , θ) + ǫ for static data , (1.41) x(t + τ) = f dyn (x(t), θ) + ǫ t for time-series , (1.42)
Where the function f is in a particular functional space. The network is learned from the parameter vector θ. I describe here the tried functional space and the procedure to identify θ.

Linear Methods

A linear dynamical model has the form:

x i (t + τ) = p j=1 θ i j x j (t) (1.43)
Where θ i j gives the importance of gene j for the prediction of gene i. One would usually add constraints to obtain a sparse θ, and consequently have an interpretable model. Thus, the LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] was used in [START_REF] Gustafsson | Constructing and analyzing a large-scale gene-to-gene regulatory network-lasso-constrained inference and biological validation[END_REF]. In biological systems, there may be a delay between the presence of a gene i and its influence on another gene j. Shojaie and Michailidis [START_REF] Shojaie | Adaptive thresholding for reconstructing regulatory networks from time-course gene expression data[END_REF] build linear models depending on the system's state at several previous times (x(t), x(t -1), . . . ). Besides, linear methods offer a framework to learn from heterogeneous data. Given dynamic data, illustrating a dynamics x(t + τ) = f dyn (x(t)), and steadystate data, illustrating x i = f stat (x -i ), one can learn a linear regression for each model θ dyn and θ stat that result from the same network. One wants θ dyn and θ stat to be null on genes not interacting with gene i, thus (θ dyn ) j = 0 if (θ stat ) j = 0, and reciprocally. Marbach et al. [START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF] used group-LASSO [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] and bootstrap samples; J.Chiquet et al. [START_REF] Chiquet | Sparsity with sign-coherent groups of variables via the cooperative-lasso[END_REF] also wanted (θ dyn ) j and (θ stat ) j to share the same sign, which they obtained through cooperative-LASSO. Another approach [START_REF] Wang | Inferring gene regulatory networks from multiple microarray datasets[END_REF] consists in learning θ separately on each dataset, then find θ, the linear regression that would be closer to all the found regressions.

The state-of-the-art method TIGRESS [START_REF] Haury | Tigress: Trustful inference of gene regulation using stability selection[END_REF] uses linear models combined with stability selection [START_REF] Meinshausen | Stability selection (with discussion)[END_REF]. Given a target gene i and a potential regulator j, TIGRESS wants to evaluate the probability that θ i j is non-null. To this purpose, it uses LARS regressors [START_REF] Efron | Least angle regression[END_REF], that iteratively select the L most important regulators for the linear model. TIGRESS runs R LARS regressor on modified data: the expression levels of candidate transcription factors are multiplied by a random number in the interval [r, 1], and the model is trained on a random subsample of the data. TIGRESS now has, for all candidate transcription factor j and all ℓ ∈ [1, L], the frequency F(i, j, ℓ) with which the TF j was selected by LARS in the top ℓ features to predict the expression of gene i. By selecting a too small value for ℓ, many TF would have 0 score; by selecting too large a value for ℓ, several TF may have the same probability 1. Thus, TIGRESS will infer an edge from j to i with an averaged version of F(i, j, ℓ), with the score s(i, j):

s(i, j) = 1 L L ℓ=1 F(i, j, ℓ) . (1.44)
s(i, j) corresponds to the area under the curve of the probability to select transcription factor j for target gene i when selecting ℓ regulators, with ℓ varying from 1 to L. An edge from gene j to gene i is inferred if s(i, j) is superior to a threshold.

Boolean Model

In the Boolean framework, a gene can only be active [START_REF] Leavitt | Deciphering the Genetic Code: Marshall Nirenberg[END_REF] or inactive (0). Let xi be the discrete state of gene i, xi = x i ≥ θ i , x = x1 , . . . , xp ′ . In the dynamical case, given a system state at time t, x(t), there exists a Boolean function giving the state at time t + 1. This function can be represented by a truth table, or a wiring diagram, as shown in Figure 1.9. For a p-gene system, there exists 2 2 p possible Boolean functions for each output. Finding the optimal Boolean function is a combinatorial problem, and its resolution by brute-force algorithms is computationally too expensive. Heuristics are used to simplify the problem. Akutsu et al. [START_REF] Akutsu | Identification of genetic networks from a small number of gene expression patterns under the boolean network model[END_REF] limit the number of regulators that one gene can have to a constant K, thus drastically reducing the search space. In the REVEAL algorithm [START_REF] Liang | Reveal, a general reverse engineering algorithm for inference of genetic network architectures[END_REF], S.Liang et al. start by identifying gene j which has highest mutual information with output gene i. They iteratively add to the list of regulators of gene i the gene which most improve the Boolean function. Boolean networks focus on generic network behavior rather than quantitative biochemicals details. Despite a simplification of the input variables, Boolean networks succeeded in retrieving meaningful biological information [START_REF] Shmulevich | Binary analysis and optimization-based normalization of gene expression data[END_REF][START_REF] Tabus | Normalized Maximum Likelihood Models for Boolean Regression with Application to Prediction and Classification in Genomics[END_REF][START_REF] Huang | Gene expression profiling, genetic networks, and cellular states: an integrating concept for tumorigenesis and drug discovery[END_REF] and are able to model the behavior of biological systems [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF]. Particularly, they allow to study attractor states for the system.

Graph

Truth table

x 1 x 2 x 3 AND Input Output x 1 x 2 x 3 x 1 x 2 x 3 0 0 0 0 0 1 0 0 1 0 0 1 0 1 0 1 0 1 0 1 1 1 0 1 1 0 0 0 0 0 1 0 1 0 1 0 1 1 0 1 0 0 1 1 1 1 1 0

Bayesian methods

Bayesian networks and Bayesian methods offer a probabilistic framework able to use noisy data and prior knowledge to infer dependencies. A causal relationship between variables is described by a directed acyclic graph G (for an example see Figure 1.10). Let Π i be the parents of gene i in graph G, i.e. all genes j such as the edge j → i exists in G. The state of gene i is assumed to depend only on the state of its parents. Given a graph G, conditional probabilities p(X i = x i |X Π i = x Π i ) are estimated. Bayesian methods search for the graph G and probabilities p(.|.) that maximizes the probability of the observed data, with a penalty term representing Directed acyclic graph Adjacency matrix x 1 (t)

x 1 x 2 x 3 x 4 x 5 G =                     0 1 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0                    
x 2 (t)

x 3 (t)

x 1 (t + τ)

x 2 (t + τ)

x 3 (t + τ)
L((x i ) i=1...N , G, η) = prior(G, η) × N i=1                        p j=1 p η (X j = x i j |X Π j = x iΠ j ) likelihood of observing x i                        likelihood of observing all data (1.45)
The prior is usually a penalty on the density of the graph, to obtain a sparse graph. Bayesian networks usually learn simplified dynamics, with discrete variables, or linear dynamics with continuous variables. Optimization of L((x i ) i=1...N , G, η) is an NP-hard problem, and is computationally prohibitive for brute-force algorithms [START_REF] Chickering | Learning bayesian networks: search methods and experimental results[END_REF][START_REF] Chickering | Learning bayesian networks is np-complete[END_REF], thus most algorithms use heuristics or approximated solutions. Friedman et al. [START_REF] Friedman | Using bayesian networks to analyze expression data[END_REF] use the Sparse candidate algorithm [START_REF] Friedman | Learning bayesian network structure from massive datasets: the sparse candidate algorithm[END_REF], in which parents of a node can only be found in a small subset of genes that are highly correlated with target gene. Another approach [START_REF] Hartemink | Using graphical models and genomic expression data to statistically validate models of genetic regulatory networks[END_REF] is the estimation of a score of several networks through variational approximation methods [START_REF] Jaakkola | Variational probabilistic inference and the qmr-dt network[END_REF][START_REF] Attias | Inferring parameters and structure of latent variable models by variational bayes[END_REF].

The assumption of an acyclic graph is false for biological networks, which contain feedback loops. This limitation is overcome by dynamical Bayesian network, see Figure 1.11. These networks can only learn from time-series data. The state of a gene at time t depends only on gene states at previous times [START_REF] Husmeier | Sensitivity and specificity of inferring genetic regulatory interactions from microarray experiments with dynamic bayesian networks[END_REF][START_REF] Zou | A new dynamic bayesian network (dbn) approach for identifying gene regulatory networks from time course microarray data[END_REF]. In the dynamical case, algorithms have been developped for learning nonlinear continuous dynamics [START_REF] Imoto | Bayesian network and nonparametric heteroscedastic regression for nonlinear modeling of genetic network[END_REF][START_REF] Kim | Dynamic bayesian network and nonparametric regression for nonlinear modeling of gene networks from time series gene expression data[END_REF][START_REF] Morrissey | Inferring the time-invariant topology of a nonlinear sparse gene regulatory network using fully bayesian spline autoregression[END_REF]. Gene expression level are usually measured by mRNA concentration, but there are other components influence a gene's dynamic, such as protein concentration or microRNA. Bayesian methods can furthermore consider those unseen variables through latent variables. Bayesian methods alternatively identify the most likely latent variables given the model and optimize the model given the data and current estimated latent variables [START_REF] Perrin | Gene networks inference using dynamic bayesian networks[END_REF][START_REF] Rangel | Modeling t-cell activation using gene expression profiling and state-space models[END_REF].

Other approaches

Other, more specific, models have been developed. Kalman filters model noisy and partially observed systems. They can model nonlinear dynamics, via variants such as "Extended Kalman filters" or "Unscented Kalman filters", used for GRN inference in [START_REF] Wang | An extended kalman filtering approach to modeling nonlinear dynamic gene regulatory networks via short gene expression time series[END_REF][START_REF] Qian | Inference of noisy nonlinear differential equation models for gene regulatory networks using genetic programming and kalman filtering[END_REF][START_REF] Noor | Inferring gene regulatory networks via nonlinear state-space models and exploiting sparsity[END_REF]. S-Systems offers very generic modeling, but are computationally intensive and no method guarantees to find the best solution for the model. They were used in [START_REF] Thomas | A model-based optimization framework for the inference on gene regulatory networks from dna array data[END_REF][START_REF] Noman | Inference of gene regulatory networks using ssystem and differential evolution[END_REF]. Neural network models mimic the way the brain functions, with neurons sending signals to others to calculate an output signal. They suffer from the same weaknesses as S-systems, and have been used in [START_REF] Ao | Ensemble of elman neural networks and support vector machines for reverse engineering of gene regulatory networks[END_REF][START_REF] Grimaldi | Regnann: Reverse engineering gene networks using artificial neural networks[END_REF]. Gaussian processes can model a nonlinear stochastic process such as gene dynamics, as shown by Aijo et al. [START_REF] Äijö | Learning gene regulatory networks from gene expression measurements using non-parametric molecular kinetics[END_REF]. Logistic regression offers a nonlinear interpretable model, useful for GRN inference [START_REF] Bonneau | The inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo[END_REF]. Lim et al. [START_REF] Lim | Okvar-boost: a novel boosting algorithm to infer nonlinear dynamics and interactions in gene regulatory networks[END_REF] used operator-valued kernels to model genes. Using partial derivatives of their models, Lim et al. were able to infer GRNs. Finally, best performances in several GRNI challenges were reached by the GE-NIE3 method [START_REF] Huynh-Thu | Inferring regulatory networks from expression data using tree-based methods[END_REF]. GENIE3 learns a model f i (x -i ) with Regression trees. For a new input variable x, f i will make a binary test, for example x j ≥ 0.5. Following the answer to this binary test, input variable will either go in the left branch or the right branch of the tree. It will land on another node, with a new binary test classifying x -i , until it lands on a leave, giving the value of f i (x -i ). After each binary test, f (x -i ), the prediction for gene i, is more determined. One can interpret the importance of feature j in f i as the amount of variance it reduces. For example, given a node N, with its left branch N ℓ and its right branch N r . Let S (resp. S ℓ , S r ) be all data from training set S that reaches node N (resp. branch N ℓ , N r ). Noting |S | the number of elements in S , Var(S ) the variance of x i for all x ∈ S , the variance reduction I of node N is computed following:

I(N) = |S |Var(S ) -|S ℓ |Var(S ℓ ) -|S r |Var(S r ) (1.46)
Importance of feature j is the sum of the variance reduction I for all nodes N where the binary test was done on x j . A randomization scheme is added. Instead of uses a single regression tree, Huynh et al. uses random forests, i.e. many trees learned on random subsamples of the data.

Assessment of GRNI methods

Each author of a GRNI method tested his method in one or several of the following manners: inferring network from real gene expression data and comparison of the inferred network to the known GRN; inferring from real data and inspection by an expert of the inferred network; inferring from simulated data and comparison to the true regulatory network. Each method has its advantages and its drawbacks. When using real data and the known GRN, the results are questionable as large-scale GRN are, to a large extend, incomplete. The known GRN is called a bronze standard. To overcome errors from our limited knowledge of GRN, some author also justified their inferred network by observing that the topology of the network they inferred was in agreement with other studies. Using simulated data, one knows the true GRN producing the gene expression data. The inferred network is compared with a gold standard, but simulated data may badly reflect the real behavior of genes. In addition to questions on a dataset's quality to evaluate a GRNI method, most methods have some hyper-parameters to select a-priori. The arbitrariness of these parameters may put the study in question. GRNI methods have also been evaluated by comparative studies. GRNI is an active field of research, and reviews quickly become outdated. I will mention two reviews: Emmert-Streib [START_REF] Emmert-Streib | Statistical inference and reverse engineering of gene regulatory networks from observational expression data[END_REF] gives an overview of correlation or mutualinformation-based methods and summarizes papers comparing these methods; Narendra [START_REF] Narendra | A comprehensive assessment of methods for de-novo reverse-engineering of genome-scale regulatory networks[END_REF] compares 32 state-of-the-art methods on 15 real or realistically simulated data. No method clearly stands out. Authors essentially conclude that some methods are always underperforming, and should not be used. Yet, comparative studies are not a perfectly fair manner to compare methods. Indeed, as the study's authors both use and evaluate methods, they might fine-tune the GRNI method or the dataset to produce better results, such as results stressing their method's advantages. This can be overcome by blind-challenges. In these competitions, organizers released gene expression data on several networks. Each contestant sent the networks inferred by his methods with no knowledge of the target network. In particular, Stolovitzky, Monroe and Califano created the Dialogue for Reverse-Engineering Assessments and Methods (DREAM) conferences and challenges [START_REF] Stolovitzky | Dialogue on reverseengineering assessment and methods[END_REF] in 2006. They released data for many bioinformatic challenges, producing a fair evaluation of researchers' methods and many benchmark sets. Currently, there has been eight editions of DREAM challenges. In particular, the issue of GRNI was tackle in challenge four of the third edition and challenge four of the fifth edition. Organizers surveyed strengths and weaknesses of each method. They also reached several conclusions [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF][START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF]. Firstly, on simu-lated data, when knock-out data are available for each gene, Z-score methods give excellent results, with a simple and fast algorithm. Secondly, results on a particular network depend mainly on the chosen model; two mutual-information-based methods would essentially find the same network. Thirdly, some methods give relatively good results: GENIE3, TIGRESS, ANOVA and CLR. Fourthly, and most importantly, the consensus of several methods usually performs better than methods taken individually. In the DREAM5 Challenge, Marbach et al. compared the results of each contestant to the results of averaging all contestants' responses. The consensus was always ranked in the three best methods, and often was the best. Besides, for consensus to perform best, methods using different modeling should be used: the consensus of a Bayesian, a mutual-information-based and a regression-based methods is expected to perform better than the consensus of three regression-based methods.

Problem formulation

In this manuscript, I suggest model-driven, nonparametric, GRN inference methods, using kernel functions and gene expression data. I set this work in the common simplified view, described in Section 1.2.4 and shown in Figure 1.2. Let p be the number of genes and A be the p × p adjacency matrix: a i j = 1 if gene j regulates gene i, a i j = 0 otherwise. To estimate this adjacency matrix A, I decompose this task into p independent tasks: for each gene i, I will estimate the row vector a i• , assuming that I observe data with an additive noise whose covariance is diagonal. As seen in Section 1.2.5, time-series data are available, allowing us to learn a model:

x i (t + τ) = f i dyn,τ (x -i ) + ǫ i dyn,t+τ , (1.47)
with x -i = (x 1 , . . . , x i-1 , x i+1 , . . . , x p ) and ǫ i dyn,t+τ a zero-mean noise. There also exists data with no time dependence, called steady-state data. They allow us to learn a model:

x i = f stat (x -i ) + ǫ i , (1.48)
with ǫ i a zero-mean noise. Those models are learned using dynamic data,

S dyn,τ = {(x -i (t 1 ), x i (t 1 + τ)), . . . , (x -i (t N τ -τ), x i (N τ ))}, or static data, S stat = {(x -i 1 , x i 1 ), . . . , (x -i N , x i N )}.
Both modeling use the same supervised learning tools, that will be described using the general notation S = {(z 1 , y 1 ), . . . , (z N , y N )}. These data are observations of a biological system, hence they are noisy and result from nonlinear functions. Gene can be modeled through parametric approaches, where one assumes that gene dynamics is ruled by particular ordinary differential equations (ODE), and one has to find optimal parameters. These approaches can tackle nonlinear problems, but, in practice, they face prohibitive computational time on large dataset. Not to mention the relevance of the ODE for gene modeling. Kernel functions, nonparametric models seen in Section 1.3.3, are good candidates for gene modeling. Very few works have applied kernel methods to GRN inference, mainly because they are not easily interpretable, so, even if they perfectly modeled genes, extracting the GRN from the kernel model would not be straightforward. This manuscript describes two GRN inference methods using kernels. In the first part, I propose a method to interpret a kernel model, called Jac. Given f i a model for gene i learned on steady-state or time-series data with kernel functions, Jac estimates

a i• through the partial derivatives of f i âij = ∂ f i ∂x j dx (1.49)
In Chapter 2, I demonstrate that, under some assumptions on the distribution of the genes x and the kernel function k, kernel methods consistently estimate the mean of partial derivatives. In Chapter 3, I describe other methods to interpret variable importance in a model f i . I suggest a method to learn from both steady-state and time-series data. I compare Jac to those other model interpretation methods and to state-of-the-art methods for GRN inference. In the second part, I develop an interpretable kernel-based model, called LocKNI:

x i ≈ f i (x -i , w i ) + ǫ , (1.50)
with w i a feature weighting parameter; w i j measures the importance of gene j for the modeling of gene i. w i is learned from data. Then I estimate the adjacency matrix:

a i. = w i (1.51)
In Chapter 4, I compare LocKNI to other kernel feature weighting methods on realistic and widespread datasets. In Chapter 5, I describe a method to incorporate prior knowledge when learning w. I show improvements when adding reasonable prior knowledge. I also compare LocKNI to state-of-the-art methods on real and realistically simulated datasets. LocKNI shows state-of-the-art performances, and also a behavior complementary to other existing methods.

Part I Estimation of Partial Derivatives for Network Inference

Given a model f i for a gene i, the partial derivative ∂ f i /∂x j should reflect the action of gene j on gene i. Indeed, if j regulates gene i, a change in gene j's concentration would have repercussions on gene i, thus ∂ f i /∂x j 0. This approach raises two questions: how should genes be modeled and how the partial derivatives should be estimated.

Concerning the first question, several characteristics of gene behavior are difficult to model; most notably, the nonlinear dynamics of concentration of gene mRNAs. Research in supervised learning developed several methods able to tackle these difficult problems, such as regression trees, neural networks or kernel-based methods. However, regression trees produce piecewise constant functions, thus their model cannot be derived, while computation of neural networks may be prohibitive for real-sized networks. On the other hand, kernel methods possess desirable properties, such as robustness to noise. Kernel methods seem therefore to be good candidates for gene modeling.

With respect to the second question, many works exist for the estimation of derivatives of univariate functions. Unfortunately, gene models will be multivariate functions, and very few works exist for this problem.

In Chapter 2, I tackle the problem of partial derivative estimation from a theoretical point of view. Given some assumptions on the observations x that should be met by genes, I give sufficient conditions on a learning algorithm K, taking training data S to learn a model f , so that, for any continuous linear form g, g(∂ f /∂x j ) consistently estimates g(∂ f /∂x j ). I then show that some kernel methods meet these conditions. Finally, I test this partial derivative estimation method on toy examples.

In Chapter 3, I implement the use of partial derivatives to infer the gene regulatory network. This method is improved by an ensemble method, allowing to learn from both steady-state and dynamic modeling. Estimation of partial derivatives ∂ f i /∂x j can be seen as a way to interpret the importance of feature j in the model f i . Therefore, I compare this network inference method to other feature importance interpretations. I also compare this network inference method to state-of-the-art GRNI methods on realistic datasets. It provides state-of-the-art results.

Chapter 2

Estimation of Partial Derivatives

Introduction

Learning a relationship f between input variables x and output y from past examples is a common task in machine learning. Many domains (computer vision, medical imaging, bioinformatics, etc) present such problems, fostering the development of many efficient algorithms for learning f , an estimation of f , from observational data. Most learning methods provide of theoretical bounds for convergence and performances. One can mention linear, tree-based or kernel-based methods. However, the partial derivatives estimation problem has fewer theoretical guarantees. Recently, several fields have shown interest in learning the partial derivatives ∂ f /∂x i . For the characterization of nanoparticles, quantitative features, such as the diameter of a nanoparticle, can be more accurately estimated by the derivative of a response function [START_REF] Charnigo | Estimating quantitative features of nanoparticles using multiple derivatives of scattering profiles[END_REF][START_REF] Francoeur | Sensitivity analysis for characterization of gold nanoparticles and agglomerates via surface plasmon scattering patterns[END_REF][START_REF] Charnigo | Estimating multiple derivatives simultaneously: What is optimal[END_REF]. Derivatives allow to quantify the progress of a disease [START_REF] Ramsay | Functional data analysis[END_REF] or to infer gene regulatory networks [START_REF] Lim | Okvar-boost: a novel boosting algorithm to infer nonlinear dynamics and interactions in gene regulatory networks[END_REF]. Besides, independently of the domain, partial derivatives may quantify the relevance of an input variable for a function, thus can be used as a feature selection criterion. Several saliency measures are based on partial derivatives, see [START_REF] Steppe | Feature saliency measures[END_REF][START_REF] Leray | Feature selection with neural networks[END_REF] and references therein. Nevertheless, theoretical consistency of derivative estimation must be proved. Technical difficulties arise because differentiation is not continuous, without assumption on the studied functional space. Thus estimating f by f through a consistent learning algorithm does not imply that ∂ f /∂x i consistently estimate ∂ f /∂x i . Consistency of several derivative estimation methods for univariate problems has been proven. The approach calculates empirical derivatives [START_REF] Brabanter | Nonparametric derivative estimation[END_REF]. Given a sorted

y ′ i = y i+1 -y i-1 x i+1 -x i .
For more robustness to noise, Charnigo et al. calculate an averaged empirical derivative [START_REF] Charnigo | A generalized c p criterion for derivative estimation[END_REF]

y ′ i = k w k y i+k -y i-k x i+k -x i-k with w k a decaying weight, k w k = 1, w k-1 ≥ w k ≥ 0 for all k.
The algorithm in [START_REF] Khemchandani | Regularized least squares support vector regression for the simultaneous learning of a function and its derivatives[END_REF] trains f so that f (x i ) fits y i and f ′ (x i ) fits y ′ i for all i. Another approach is to model f with locally polynomial functions, such as splines [START_REF] Rondonotti | Sizer for time series: a new approach to the analysis of trends[END_REF]. Charnigo et al. [START_REF] Charnigo | A generalized c p criterion for derivative estimation[END_REF] learn a model f (x) = ℓ w ℓ (x)µ ℓ (x), with µ ℓ a polynomial and w ℓ a function determining if x is in the region modeled by µ ℓ . With w ℓ (x) = exp(-(x-x l )/h), Charnigo [START_REF] Charnigo | Self-consistent estimation of mean response functions and their derivatives[END_REF] proves consistency of derivatives of order j. For multivariate problems, there is only, to our knowledge, Mosci's new kernel-based approach that has been proven to estimate ∂ f /∂x i 2 consistently [START_REF] Mosci | Nonparametric sparsity and regularization[END_REF]. In this work , I firstly give sufficient conditions on a target function f and a learning algorithm to consistently estimate g(∇ f ), with g any continuous linear form, and ∇ f the gradient of f , ∇ f = (∂ f /∂x1 , . . . , ∂ f /∂x p ). Secondly, I prove that Gaussian kernel methods satisfy these conditions. Thirdly, I suggest to approximate the integrals by sum, for faster and simpler computation. Finally, I experiment this partial derivative estimation method on toy examples.

Consistent estimation of mean partial derivatives

In the following, (X, Y) are random variables from Ω × Y, with Ω a finite pdimensional space, and Y is a bounded subset of R. X is drawn from a distribution µ. I make the following assumptions:

(A-1) The support of µ, noted X o , is bounded and convex. The closure of X o is denoted by X c . (A-2) The random variable Y follows y = f (x) + u (2.1)
with u a zero-mean noise. I further assume that f is continuous and continuously differentiable on X c , i.e. belongs to the space C 1 (X c , Y).

In theory, assumption (A-1) is met by many empirical distributions, notably, the truncated Gaussian distribution. In practice, many random variables are assumed to follow a Gaussian distribution and the loss of generality induced by using a truncated Gaussian distribution is very small. Besides, assuming that gene concentration is bounded seems reasonable; the convexity assumption of reachable states is more debatable. Assumption (A-2) describes the link between Y and X and the additive nature of the noise. It makes mild constraints on f , namely that ∇ f exists and is continuous, criteria met by many biological models. I will use the following norm notations. L 2 (Ω, µ) is the space integrable function, according to distribution µ, with the norms:

f 2 L 2 (Ω,µ) = Ω f (x) 2 µ(x)dx = X c f (x) 2 µ(x)dx , (2.2)
f ∞ = max x∈X c | f (x)| . (2.3)
The Euclidian norm in space Ω is x 2 = p i=1 (x i ) 2 . The last needed notations and assumptions concern the learning algorithm. Let fℓ be f 's estimate by algorithm K with ℓ data pairs, (x i , y i ) i=1...ℓ . I assume that:

(A-3) fℓ is a consistent estimator: for all ǫ, η > 0, there exists an integer ℓ 0 such that, for all integer ℓ ≥ ℓ 0 , the probability that f -fℓ L 2 (Ω,µ) ≤ ǫ is greater than 1η, and

(A-4) there exists a integer ℓ 0 such that, for all integer ℓ ≥ ℓ 0 , the norm of fℓ 's gradient is bounded by a constant M.

Machine Learning research has developed many algorithm satisfying (A-3). Assumption (A-4) is satisfied by some of them, in particular kernel methods, as we shall see later. If all assumptions are met, any continuous linear form of the derivatives of f will be consistently estimated by the one of fℓ . This is proven using the following two theorems: Theorem 2. Let f be a function belonging to the space C 1 (X c , Y). Assume the gradient of f is bounded by M, ∇ f 2 < M. Then, for any ǫ > 0, there exists C 1,ǫ , independent of f , such as:

f ∞ ≤ C 1,ǫ f L 2 (Ω,µ) + ǫ .
(2.4)

Theorem 3. Let f be a C 1 (X c , Y) function, with gradient bounded by M. Let g be a linear continuous form of C 0 (X c , Y), g : C 0 (X c , Y) → R.
Then, for any ǫ > 0, there exists a constant C 2,g,ǫ such that:

|g(∇ f )| ≤ C 2,g,ǫ f ∞ + ǫ (2.5)
These two theorems are proven in the appendix, pages 128 for theorem 2 and page 129 for theorem 3). I show next that kernel-ridge regression and partial leastsquare regression learning algorithms both satisfy assumptions (A-3) and (A-4).

Kernel methods for partial derivative estimation

There exist universal kernels, i.e. kernels whose RKHS H is the whole space of continuous functions from X c to R (or C) [START_REF] Micchelli | Universal kernels[END_REF]. Using consistency of kernelridge [START_REF] Caponnetto | Optimal rates for the regularized leastsquares algorithm[END_REF] or partial least-square [START_REF] Blanchard | Kernel partial least squares is universally consistent[END_REF] regression and a universal kernel, I obtain a learning algorithm K satisfying assumption (A-3). I give sufficient conditions on a universal kernel k ensuring that the associated learning algorithm K satisfies (A-4):

Lemma 4. Let H be the RKHS of universal kernel k. If, for all x ∈ Ω • the kernel is constant on the set of points (x, x), k(x, x) = c,
• at point (x, x), the gradient of the kernel is null ,∇ z k(x, z)| z=x = 0 p , and

• at point (x, x), the Hessian matrix H(x, z) i j = ∂ 2 k(x,z) ∂z i ∂z j has eigenvalues bounded by a constant M, | < u, H(x, x)u > | ≤ M u 2 for all x ∈ X and all u ∈ R p
then, for all f ∈ H and for all x ∈ X c , the gradient of f satisfies:

∇ f (x) 2 ≤ √ M f H . (2.6)
Sketch of proof. The RKHS property and Cauchy-Schwarz inequality give:

| f (x) -f (x + h)| = < f, φ(x) -φ(x + h) > H (2.7) ≤ f H φ(x) -φ(x + h) H , (2.8)
while we also have:

| f (x) -f (x + h)| = | < ∇ f, h > +o( h )| (2.9)
Using the polarization identity, the distance φ(x)φ(x + h) 2 H can be expressed with scalar products. Using the kernel trick, I can express this distance with function k. Using a Taylor expansion and conditions of lemma 4:

φ(x) -φ(x + h) 2 H = k(x, x) + k(x + h, x + h) -2k(x, x + h) (2.10) = h ′ H(x, x)h + o( h 2 ) (2.11) φ(x) -φ(x + h) H ≤ √ M + o( h ) (2.12)
Equations (2.8) and (2.12) gives the result. A detailed proof is given in the appendix, page 130. The conditions of lemma 4 are satisfied by the Gaussian kernel.

Lemma 5. The Gaussian kernel of bandwidth σ satisfies the hypothese of lemma 4.

The universality of the Gaussian kernel is proven in [START_REF] Micchelli | Universal kernels[END_REF]. The other conditions are proven using calculus. This is done in the appendix, page 131. Consistent estimation of continuous linear forms of partial derivatives with Gaussian kernel methods is obtained by lemma 5 then lemma 4, then theorem 2 and, finally, theorem 3 can be applied, giving us the following theorem: Theorem 6. With samples (x i , y i ) i=1...ℓ i.i.d., for any ǫ and η > 0, and for any continuous linear form g, g : C 0 (X c , Y) → R, there exists ℓ 0 such that if ℓ ≥ ℓ 0 , then, with probability greater than 1η,

X c g(∇ fℓ (x)) -g(∇ f (x))dx ≤ ǫ , (2.13)
where fℓ is the estimator of f based on Gaussian kernel ridge regression or Gaussian partial least-square regression.

A detailed proof is given in the appendix, page 132. Now, the term

X c
g(∇ fℓ (x))dx may be hard to compute. I give below a simple and fast approximation of this integral.

Approximation of integrals

An analytical formula for the integral of the estimate f over X may not be available. To overcome this difficulty, I suggest a simple method to approximate this integral by a sum. By the central limit theorem:

1 ℓ ℓ i=1 ∇ f (x i ) -→ ℓ→+∞ X c ∇ f (x)µ(x)dx (2.14)
Moreover, as the gradient ∇ f is bounded, the central limit theorem gives bounds on the difference between the sum and the integral, provided in the appendix, page 133. If one wants to estimate a linear form different from equation (2.14), one can use the following sum, assuming that the distribution µ is known:

1 ℓ ℓ i=1 h(x i ) µ(x i ) ∇ f (x i ) -→ ℓ→+∞ X c ∇ f (x)h(x)dx (2.15)
In most cases, the distribution µ is unknown, but several accurate methods exist to estimate it. In particular, multivariate kernel density estimation [START_REF] Simonoff | Smoothing methods in statistics[END_REF] provides nonparametric and consistent estimation of µ. Noting κ(u) = (2π) -p/2 exp( u 2 2 /2), a ℓ a sequence such that a ℓ -→ ℓ→+∞ 0 and ℓa ℓ -→ ℓ→+∞ +∞ (for example, a ℓ = ℓ -1/2 ), Simonoff defines:

μℓ (x) = 1 ℓa p ℓ ℓ i=1 κ x -x i a ℓ (2.16)
This estimator μ satisfies:

μℓ -µ ∞ -→ ℓ→+∞ 0 (2.17)
Proofs and bounds can be found in [START_REF] Giné | Rates of strong uniform consistency for multivariate kernel density estimators[END_REF]. Thus, one can estimate the integral according to:

X c ∇ f (x)h(x)dx ≈ 1 ℓ ℓ i=1 h(x i ) μℓ (x i ) ∇ f (x i ) (2.18)

Numerical simulations

In this section, I illustrate the consistency of kernel-ridge regression for the estimation of continuous linear forms of partial derivatives. I build toy examples with input variables x drawn independently and uniformly in [-1, 1] p . The output variable y is a function of the following form:

y = f (x) + gǫ = tanh x T Ax + b T x sin c T x + gǫ (2.19)
with ǫ a normal noise, ǫ ∼ N(0, 1), uncorrelated between samples, and A, b and c independently chosen random parameters:

A ∼ U([0, 1] p×p ) (2.20) b ∼ U([0, 1] p ) (2.21) c ∼ U([0, 2π] p ) (2.22)
See Figure 2.1 for examples of functions given by equation (2.19).

In order to measure the error in partial derivative estimation, I define E, the error in partial derivative estimation for one simulation as follows:

E = p m=1 X c ∂ f * ∂x m µ(x)dx - X c ∂ f * ∂x m µ(x)dx 2 (2.23)
For fixed parameters (number of training points N, dimension of input variables p, noise to signal ratio σ, with g = σ var( f (x))), I simulate 100 models following equation (2.19) and plot the mean and variance of the E in box plots. Figure 2.2 shows that the error decreases with the number of training points. In these examples, the problem dimension p is 2. The noise to signal ratio is 0% in the top left plot, 100% in the top right plot. The bottom plots show the evolution of the error with the noise to signal ratio σ (plot at the bottom left), and with the input dimension p (plot at the bottom right). As expected, the error increases with the noise-to-signal ratio. Besides, the estimator is very sensitive to the dimension of input variables. When p increases, the error greatly increases. We also dispose of a partial derivative estimation to estimate other linear forms of the partial derivatives. To simulate this, I draw x following a centered and reduced Gaussian distribution,

X ∼ N(0, I p ) , (2.24)
E = p m=1 [-1,1] p ∂ f * ∂x m dx - [-1,1] p ∂ f * ∂x m dx 2 .
(2.25)

Results are shown in Figure 2.3, with p = 3 and σ = 25% . These results illustrate the property that kernel-ridge regression consistently estimates partial derivatives.

Conclusion

In this chapter, I tackled the partial derivative estimation problem. Observing the scarcity of results for the multi-dimensional case, their specificity to the learning algorithm used, I showed that, under a few assumptions on the learning algorithm, linear forms of partial derivatives were consistently estimated. Furthermore, I showed that kernel methods met these assumptions. In addition, if a kernel method is proved not to require the i.i.d. assumption to learn f , this results and the theorems in this chapter can be used to prove consistency of the estimates of partial derivatives without the i.i.d. assumption. Finally, the proposed method to estimate derivatives was tested on toy examples. We saw that this method gives good quality results in the low-dimensional case, but its performance decreases quickly with the dimension of the input. In the following chapter, I will use partial derivatives to infer gene regulatory networks.

Introduction

As shown in the first chapter, many model-driven approaches have been developed for gene regulatory network inference. Although kernel methods provide efficient nonparametric modeling, they have rarely been used for gene network inference, mainly due to their "black box" behavior. However, we saw, in the previous chapter, that partial derivatives of a regression model based on a universal kernel, typically a Gaussian kernel, can be used to provide consistent estimates of the mean of the partial derivatives of the target function. Therefore, I propose here to use the partial derivatives to interpret a kernel-based regression a posteriori. Calculating the empirical mean of the partial derivative for a given feature on data is seen as an importance measure of this feature. This chapter presents a study of this new GRN inference method. To improve this measure and to deal with large dimensions, this approach is extended to an ensemble of randomized kernel-based models. The whole approach is presented in Section 2. Then, related works about other measures of feature importance are reviewed. Section 4 describes methods for hyper-parameters selection. Section 5 is devoted to the numerical experiments which involve results on small-scale realistic networks and large scale, real-sized biological networks. Jac, the GRN inference method using partial derivative, is compared to other importance measures applied on kernel-based models.

Partial derivatives of kernel-based models

Jacobian matrix estimation for GRN inference

Several works have studied the importance of a feature through partial derivatives, see [START_REF] Steppe | Feature saliency measures[END_REF][START_REF] Leray | Feature selection with neural networks[END_REF][START_REF] Mosci | Nonparametric sparsity and regularization[END_REF] and references therein. Assume that the behavior of a gene regulatory network is governed at steady-state by the following p models f i :

∀i ∈ {1, . . . , p}, ℓ ∈ {1, . . . , N}, x i ℓ = f i (x ℓ ) + ǫ i ℓ , (3.1)
where S = {x 1 , . . . , x N } are the observed dataset and f i is the true model governing the expression level of gene i in function of the expression levels of the other genes. If gene j excites (resp. inhibits) gene i, then its action should be observed in the partial derivatives: ∂ f i /∂x j should be positive (resp. negative). Note that the effect of a gene j on gene i may be visible only on some of the system's possible states X. For example, gene j may influence gene i only if its concentration is above a threshold θ j , so that ∂ f * ,i • /∂x j = 0 if x j ≤ θ j . Similarly, beyond another concentration γ j , gene i's receptors to gene j's action may saturate, and

∂ f * ,i • /∂x j = 0 if x j ≥ γ j . Therefore the mean partial derivative ∂ f * ,i
• should be estimated on all possible states X. Let µ be the distribution of the system state, I want to estimate the average value of the (i, j) th entry of the Jacobian matrices:

Jac( f ) i, j = X ∂ f i (x) ∂x j µ(x)dx, (3.2)
As seen in Chapter 2, the empirical mean of the partial derivatives of the estimated models f i consistently estimate Jac( f ):

Jac( f ) i, j = 1 |S| x∈S ∂ f i ∂x j (x) . (3.3)
Similarly, the temporal behavior of the gene regulatory network may be assumed to be governed by an autoregressive model that decomposes into p models

f i τ : ∀i ∈ {1, . . . , p}, ℓ ∈ {0, . . . , N -1}, x i (t ℓ + τ) = f i (x(t ℓ )) + ǫ i t ℓ , (3.4) with D τ = {x(t 0 ), . . . , x(t N )}, t ℓ = t 0 + ℓτ and ǫ i t ℓ are i.i.d. realization of a gaussian noise. Then Jac( fτ ) i, j = 1 N t∈{t 1 ,...,t N } ∂ f i τ (x(t)) ∂x j (t) . (3.5)
These Jacobian estimates of the learned model can be used in order to get an estimate of the binary and asymmetric target adjacency matrix of the network. One way to provide such an estimate is to threshold the absolute value of the Jacobian's coefficients given a threshold θ:

Âi,j = H(| Jac( f ) i, j | -θ), (3.6)
with H : R → {0, 1}, the indicator function of R + .

In the remaining part of the chapter, I will consider the problem of the estimation of the i th row in adjacency matrix A. For the sake of simplicity, let i, the index of a target gene, be fixed and consider the following base model that satisfies the following equations:

y ℓ = f (z ℓ ) + ǫ ℓ (3.7)
with S = {(z ℓ , y ℓ ), ℓ = 1, . . . N} and ǫ ℓ i.i.d. realization of a Gaussian noise. Input variables have finite dimension z ∈ R p Both these estimates would greatly suffer from high-dimensionality of the input variable z. To overcome this difficulty, I propose to learn functions f that uses few dimensions of the variables. With an ensemble method, many functions f are learned using Gaussian kernel-ridge regression, but each function relies on a small subset of features of z-another approach is changing the base learned to produce a sparse Jacobian matrix, and will be developed in Chapters 4 and 5. Moreover, note that other partial derivatives open many possibilities to infer the GRN: for example, one may consider that, if gene j regulates gene i, there must exists a state z where gene j exert much its influence on gene i, and the partial derivative ∂ f i /∂z j will reach a high value, thus estimating the GRN by the maximum value obtained by the partial derivatives. In this thesis, I was interested in the average Jacobian.

Ensemble of randomized kernel-based models

From a very general point of view, performances have been improved in many domains by using an ensemble method approach. Instead of using a function f learned with one algorithm K on the whole dataset S, one can create B subsets of S: (E 1 , . . . , E B ), and learn a prediction function h b on each dataset E b . Using the mean function 1 B B b=1 h b is usually more accurate and more robust to noise that using the single function h learned using the whole dataset (see [START_REF] Breiman | Random forests[END_REF][START_REF] Meinshausen | Stability selection (with discussion)[END_REF] for some examples). In this Chapter, there are many motivations to use ensemble methods. In order to estimate f in 3.7, a Gaussian kernel-based model is considered:

h(z) = N ℓ=1 α ℓ k(z ℓ , z) (3.8) k(z ℓ , z) = exp - ||z ℓ -z|| 2 2σ 2 (3.9) (3.10)
Then, the partial derivatives are given by :

∂h(z) ∂z j = N ℓ=1 α ℓ ∂k(z ℓ , z) ∂z j (3.11) ∂k(z ℓ , z) ∂z j = z j ℓ -z j σ 2 exp - ||z ℓ -z|| 2 2σ 2 (3.12) (3.13)
As the Gaussian kernel uses each input feature in the same way, the partial derivatives regarding two features j and m only differs by the factors (z j ℓ -z j ) and (z m ℓ -z m ). To increase this difference, the kernel-ridge regression is used as a base learner in an ensemble method, each base learner is trained on a random subspace of R p and a subsample of the data. Contrary to random forests, the base learner does not select further features; however, for a given size of the training data, models trained on smaller subspaces are expected to give more contrasted partial derivatives. Moreover, working in random subspaces allows to tackle a large number of input features. Finally, I suggest an ensemble method approach allowing to learn from heterogeneous datasets, here steady-state data and irregularly-observed time series.

Learning from heterogeneous data types

In this work, my objective is not obtaining the most accurate prediction function f , but identifying the important features of gene i's behavior. I observe gene i under different perspectives (steady-state or dynamic) in dataset (E stat , E dyn ), but I assume that both behaviors are the results of the same GRN. Thus, learning importance of features from both of these perspectives gives complementary information on the same network.

Dataset building

All measures can be used as static data E stat = (x 1 , . . . , x N ) to learn a model x i = f stat (x -i ) + ǫ stat , with ǫ stat a zero-mean noise, as done by several state-of-the-art methods [START_REF] Huynh-Thu | Inferring regulatory networks from expression data using tree-based methods[END_REF][START_REF] Haury | Tigress: Trustful inference of gene regulation using stability selection[END_REF]. In addition, these data include time-series, for which we also have information about time dependencies. We have at our disposal N dyn data from several time-series. For a time-series u, the system is observed at times t 1,u , . . . , t T u ,u . The observation of the system at time t k,u in time-series u is denoted x(t k,u , u). These data can be used to learn a dynamical model of the system:

x i (t + τ, u) = f τ (x -i (t, u)) + ǫ τ ∀t, u (3.14)
for all 1 ≤ u ≤ N u , and with ǫ τ a zero-mean noise. To build a model f τ from time-series, there must be regularly spaced time-points, which is rarely the case in biological time series. Indeed, when the time-series is the response of the system to an exogenous perturbation, most behavior will be observed soon after the perturbation. After a long time, the system has stabilized. Thus, experimentalists take many measurements right after the perturbation, and few when much time has passed. To overcome the problem of irregularly spaced observation time, I build for all possible time steps τ the datasets E dyna,τ = {(x(t 1 , u 1 ), x(t 1 + τ, u 1 )), . . . , (x(t N ττ, u N τ ), x(t N τ , u N τ ))} of all time point pairs (t, t + τ) available in the data, see Figure 3.1 for a visual example. Each of these datasets allows to learn a model f τ , but not every dataset is interesting: some of them contain too little data, others have a time-step τ very long compared to the characteristic time of the system, and only a return to a steady-state is observed in the model f τ . Without expert knowledge on a system, I relied on the experimentalists' competence, and assumed that the time τ * , for which E τ * is the largest dynamical dataset (τ * = arg max τ |E τ |), is below the characteristic time of the system. Thus I accepted dataset E τ if τ ≤ τ * and |E τ | ≥ 30.

Subspaces and subsamples building

From available datasets (E 0 , E 1 , . . . , E N e ). I build B subsamples E ℓ , ℓ = 1 . . . B with subbagging and subsampling of variables, in a similar way to extremerandomization [START_REF] Haury | Feature selection for gene expression data: molecular signatures for breast cancer outcome prediction and gene regulatory network inference[END_REF]. Subsamples are built according to the following procedure: n var is a hyper-parameter to fix. The set of selected variables is called G ℓ .

Let f i ℓ be the model for gene i learned on subsample E ℓ . One may directly use Jac( fℓ ), but there is no unit nor references to assess that the value Jac( fℓ ) i j is significantly high. The p-value statistical test will give more meaning to the observed Jac( fℓ ) value.

p-value

Several ways exist to infer a GRN from an importance measures, as Jac( f ) i j can be seen as a measure of the importance of gene j to predict gene i-related to saliency measures. One can directly use the importance measure: for example, one can assume that j regulates i if the importance measure is above a threshold or if it is above a certain percentage of measured importances. I shall abbreviate this interpretation method by dir. Another approach is to compute the p-value of Jac( f ) i j . This is the probability of observing such an importance value under the null hypothesis that gene i and j are independent. For example, one can assume that gene j regulates i if the pvalue of observed measured importance of gene j for gene i is smaller than 5% or 1%. I obtain the p-values by a permutation test, theorized by Fisher [START_REF] Fisher | The design of experiments[END_REF] and Pitman [START_REF] Pitman | Significance tests which may be applied to samples from any populations[END_REF]. This test, also called randomization test or re-randomization test, is an exact nonparametric test, making no assumptions on underlying distributions, thus it is legitimate to apply it here. Given an importance measure, Impo, a training set S = (x i , y i ) i=1...N , the test is performed as follows:

1. Randomly choose π, a permutation of {1, . . . , N}. x π(i) is now independent of y i 2. Calculate Impo( j) for the dataset (x π(i) , y i ).

3. Repeat step (1-2) M times -in this work, I chose M = 10000. We now have the distribution of Impo( j) with j independent from output, hence p-values.

Among the advantages of the permutation test are its theoretical soundness and the clear interpretability of its results. Its main drawback is the associated computational time. By grouping several calculations, the permutation test may be done by performing M matrix multiplications. Its full computational cost is O(pMN 2 ), which is the same order of magnitude than the O(pN 3 ) cost for solving each kernel-ridge regression.

Voting procedure

The introduction of a permutation test provides a way to calculate p-values and use it to estimate the target matrix. Let L im be the indexes of subsamples allowing to learn importance of gene m for gene i, i.e. L im = {ℓ, m ∈ G ℓ , i G ℓ }. Let V meth be the matrix such as (V meth ) im is the mean vote for a regulation of gene i by gene m for method meth:

V Impo im = 1 |L im | ℓ∈L im Impo(i, m, ℓ) direct use of Impo V pVal im = 1 |L im | ℓ∈L im -pVal(i, m, ℓ) use of Impo's p-value V pVal<5% im = 1 |L im | ℓ∈L im (pVal(i, m, ℓ) < 5%) subsample vote 1 iff pVal(i, m, ℓ) ≤ 5% V pVal≤1% im = 1 |L im | ℓ∈L im (pVal(i, m, ℓ) ≤ 1%) subsample vote 1 iff pVal(i, m, ℓ) ≤ 1%
The adjacency matrix is estimated by method meth by inferring an edge from variable m to variable j if (V meth ) m j exceeds a threshold θ ≥ 0 âmj = 1 θ ((V meth ) m j ) .

(3.15)

Other measures for feature importance in a model

I proposed one way to measure the importance of a feature in a model using the model's partial derivatives. Alternatively, other measures can be used and are compared to the partial derivatives estimate: (a) information loss resulting from ignoring this feature, (b) the effect of this feature on the functional cost. Inspired by related works, I propose other a posteriori measures of feature importance of a (nonparametric) base model in the ensemble approach.

Feature relevance

Feature relevance has been defined in many ways (see [START_REF] Kohavi | Wrappers for feature subset selection[END_REF] for various formulations and their shortcomings). In particular, Kohavi and John [START_REF] Kohavi | Wrappers for feature subset selection[END_REF] have defined the notion of strong relevance in case of supervised classification. A feature j is strongly relevant if removal of this feature alone will result in performance deterioration of a Bayes optimal classifier. They give the formal definition:

Definition 7. (Strong relevance for supervised classification).

A feature Z j is strongly relevant if and only if there exist some z j , y and z -j such that

p(Y = y|Z j = z j , Z -j = z -j ) p(Y = y|Z -j = z -j ) (3.16)
In case of regression, a similar definition would involve the regression function: z → E[Y|z], being optimal according to the mean squared error criterion.

Definition 8. (Strong relevance for regression).

A feature Z j is strongly relevant if and only if there exist some z j , y and z -j such that

E[Y|Z j = z j , Z -j = z -j ] E[Y|Z -j = z -j ]
(3.17)

However, this optimal model is not available. Therefore, I suggest to evaluate the relevance of a feature j by the change in the mean-square error of the kernelridge regression based on all features, and based on all features but j. Let f be the regression function trained on all features, and fj the one trained on all features but j. Let S train be the set of data used for training f and fj . Let S test be the test data. I define a empirical relevance measure based either on training set or on test set.

rel train ( j) = 1 |S train | (z,y)∈S train y -f (z) 2 -y -f -j (z -j ) 2 , (3.18
)

rel test ( j) = 1 |S test | (z,y)∈S test y -f (z) 2 -y -f -j (z -j ) 2 . (3.19) (3.20)
If feature j is relevant, there will be a loss of information; fj should be less accurate than f , and rel train ( j) and rel test ( j) will be positive. Note that, as we are not using the regression E[Y = y|z] optimal function, rel test ( j) may be negative. If feature j is irrelevant, its removal from training data will reduce noise, and allow fj to be a better regression function than f .

Sensitivity analysis

Le Cun et al. [START_REF] Lecun | Optimal brain damage[END_REF] suggest that the sensitivity of a learning algorithm be measured for each feature. In their work, a feature's importance is the change in the functional cost caused by the removal of the feature. Following Guyon et al. [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF], I measure the sensitivity as follows:

• Train a regression function f (z) = z i ∈S train αi k(z i , z) with algorithm K. Calculate the kernel-ridge functional cost

L(0) = 1 |S train | (z i ,y i )∈S train          y i - z j ∈S train α j k(z j , z i )          2 (3.21) +λ (z i ,z j )∈S 2 train α j αi k(z j , z i ) (3.22)
• With the same coefficients α, calculate the functional cost with kernel k -j , the adaptation of kernel k to use all features but feature j

L( j) = 1 |S train | (z i ,y i )∈S train          y i - z j ∈S train α j k -j (z j , z i )          2 (3.23) +λ (z i ,z j )∈S 2 train α j αi k -j (z j , z i ) (3.24)
• Define the sensitivity of algorithm to feature j by:

sens( j) = L(0) -L( j) (3.25)
The value of sens( j) quantifies how feature j affects the functional cost using all features L(0). Note that this measure follows [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF], but one may also consider its normalized variant (L(0) -L( j))/L(0). These two importance measures will also be implemented with the previously described ensemble methods.

Hyperparameter selection

Several parameters need to be fixed. The kernel function used is the Gaussian kernel

k(z, z ′ ) = exp z -z ′ 2 2 2σ 2 (3.26)
since it is a very popular kernel, with good performance in practice and universal consistency in theory. There remains four hyper-parameters to choose:

• For the subsampling parameters, cross-validation would be a biased procedure. If a model was trained on more data, it would have lower test error. Moreover, by using fewer variables, many models would not have any relevant feature as explanatory variable, and would have a high test error. Our experiments have shown that the subsampling parameters n var and p data have little influence on the GRN inference performance, see I consider that B = 100×p n var is enough subsamples, the vote matrix has converged.

Experiments on small-scale networks

Data: DREAM3 Challenge

The Dialogue for Reverse-Engineering Assessment and Methods (DREAM) challenge allows researchers' team to compete on bioinformatic problems. In the 3 rd edition [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF][START_REF] Prill | Towards a rigorous assessment of systems biology models: the dream3 challenges[END_REF], the 4 th challenge focused on gene regulatory network inference from realistically simulated networks. The networks are modules from known networks of E.Coli and Yeast [START_REF] Marbach | Generating realistic in silico gene networks for performance assessment of reverse engineering methods[END_REF]. Data were generated with GeneNetWeaver [START_REF] Schaffter | Genenetweaver: in silico benchmark generation and performance profiling of network inference methods[END_REF].

With realistic simulations, one knows the real network behind the data, thus one can confidently compare the performances of network inference methods.

The challenge provides networks of size 10, 50 and 100 genes. For each size, there are 5 networks, 2 extracted from E.Coli and 3 from Yeast. For each network, there is a knock-out and a knock-down experiment for every gene, and, for size 10 networks (resp. 50, 100), 4 (resp. 23, 46) time-series of 21 observation points. Knock-down and knock-out experiments bring a lot of information by their metainformation (which gene is perturbed). The challenge's organizers showed that best performances could be reached by only using the perturbed data and Z-score (described in section 1.4.1, page 28). As our methods do not integrate this metainformation, they under-perform with respect to methods using the complete information. Thus, I ignore the perturbed data and learn the network using only time-series data. I compare our results to the best teams of the DREAM3 challenge also using only time-series data.

Performance evaluation

Given a GRN inference method, one obtains a score for each interaction and assumes that there is a regulation from gene i on gene j if  j,m is superior to a threshold θ. Noting R the number of predicted regulations, this prediction is compared with gold standard. By counting the number of true positives (T P), false positives (FP = R -T P), false negatives (FN) and true negatives (T N), performance can be calculated as precision, i.e. the percentage of real regulations in the inferred regulations (pre = T P/(T P + FP)), recall, i.e. percentage of all real retrieved regulations (rec = T P/(T P + FN)), or false positive rate, i.e. percentage of non-regulation predicted as regulation (FPR = FP/(FP + T N)). These statistics depend on the threshold θ. To avoid the problem of poorly chosen θ, one draws the "Receiver Operating Characteristic" curve (ROC: (FPR, pre) = f (θ)) and the precision-recall curve (PR: (pre, rec) = f (θ)). The better the method, the higher the area under these curves. For example, with a perfect predictor, the area under the ROC curve will be 1. With a random predictor, it will be 0.5.

Performance without ensemble methods

These results are shown to complete this work, in Figure 3.6. Results are, as expected, lower if no ensemble method is used. Performances of each importance measure on 50-genes networks are shown in Figure 3.6, using static or dynamic modeling. Overall, static modeling exhibits better performance for both AUROC and AUPR. No importance measure has good performance.

Analysis of p-value

The p-value is the probability of observing such an extreme value under an assumption of independence between input variable x and output y. Thus, some false positive should occur. Considering that a regulation from gene j to gene i exists if the p-value of Impo(i, j) is inferior or equal to θ; calculating M importance measures, we expect θM false regulations, noted EFP, expected false positives. Noting P, for positive, the number of edges considered, the ratio P/EFP quantifies how this importance measure produces extreme values compared to values obtained with features independent from the output. If this ratio is below 1, it means that this importance measure does not produce extreme values. This does not mean that the measure does not rank features well. We see in Table 3.1 the values of the mean of the E/EFP ratio on all 5 networks, and on how many networks this ratio was superior to 1 for the two thresholds θ = 5% and θ = 1%. From this perspective, both sensitivity analysis and partial derivatives' integral do not produce extreme values. On the contrary, relevance on the training data seems to clearly separate independent and dependent cases. Surprisingly, relevance on test data is not as extreme. This difference may arise from the difficult nature of the data. We have few noisy and incomplete data; for example, data contains only RNAm concentrations, and no protein concentration or exogenous stimuli. This may cause the poor generalization of our model, thus the irrelevance of measuring mean-square error on test data.

The plot in Figure 3.2 shows the performance of the adjacency matrix estimate V meth as a function of the number of subsamples, on the DREAM3 size 50 networks for two importance measures: relevance on training data and integral of partial derivatives. The behavior is similar for all four importance measures. For both importance measures, we observe that V •,dir and V •,pVal converge faster than V •,pVal<5% and V •,pVal<1% . Also, with the two latter methods, we observe some gaps, which occur sometimes after many subsamples. Thus, the direct and pVal methods appears to perform better, as far as speed of convergence is concerned.

Performances on DREAM3 50-gene and 100-gene Networks

On size 50 datasets, results are shown in Table 3.4 for AUROC and in Table 3.5 for AUPR. Note that, aside from relevance on training data, using the importance measure directly or through p-values makes very little difference. For relevance on training data, p-values seem to enhance performance, but the change is very small. Figure 3.3 compares each directly used importance measure . Partial derivatives's integral clearly stands out of other methods. On four networks, this method has the highest AUROC; on three networks, it has the highest AUPR. Moreover, Jac has a large difference with second best methods on some networks (+0.09, +0.12 AUROC on networks Ecoli1 and 2), and it is never far from the best performing method (in the worst cases, Jac has -0.01 AUROC and -0.015 AUPR with respect to the best performing method).

Those results are confirmed on 100-gene networks. Jac has the best AUROC on four datasets, and best AUPR on all five. Note that these are difficult datasets. I compared Jac with other GRN inference methods: TIGRESS [START_REF] Haury | Tigress: Trustful inference of gene regulation using stability selection[END_REF], GENIE3 [START_REF] Huynh-Thu | Inferring regulatory networks from expression data using tree-based methods[END_REF], a dynamic Bayesian method, called G1DBN [START_REF] Lèbre | Inferring dynamic genetic networks with low order independencies[END_REF], and the best DREAM3 team using only the time-series data. Results are shown in Figure 3.4. Jac is the best of the five methods in four datasets out of five both in AUROC and AUPR. This shows that a naive implementation of estimation of partial derivatives is very competitive with state-of-the-art methods. 3.2: Characteristics of datasets: # TF is the number of potential transcription factors, p is the number of genes of the system, N is the number of data points, N dyn is the number of data from time-series, among the N data

Method

Experiments on real and real-sized networks

Jac importance measure has been evaluated on four well-referenced benchmarks. The scales and sizes of the datasets are summarized in Table 3.2.

Datasets

The datasets were found in the DREAM5 challenge [START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF]. The first dataset is a realistic simulation of genes-see [START_REF] Marbach | Generating realistic in silico gene networks for performance assessment of reverse engineering methods[END_REF] and [START_REF] Schaffter | Genenetweaver: in silico benchmark generation and performance profiling of network inference methods[END_REF] for a complete explanation-while the last two are real data. For each dataset, gene expression data of heterogeneous types is available: steady-state, time-series and perturbed data. We have a list of potential regulators, called "Transcription Factors" (TF). The gold standard is given by the DREAM organizers. I compared Jac with state-of-the-art methods. DREAM5 organizers provided the contestants' inferred networks. I compared my results to those of the best contestants.

Using other GRN Inference methods: A consensus approach

One talks of wisdom of crowds [START_REF] Surowiecki | The wisdom of crowds: Why the many are smarter than the few and how collective wisdom shapes business[END_REF] when a group, using the vote of each individual, makes better decisions than any individual alone. This phenomenon is well illustrated in gene regulatory network inference. Using different methods, the mean of all inferred networks is often closer to reality than each network taken alone, as Marbach et al. have shown with the DREAM5 challenge [START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF].

To average each prediction, Marbach et al. sorted, for each method, which link was considered most likely. They had n J inference methods. Method i gives a score for each interaction, score stored in B i ∈ R p×p the score matrix of method i.

The higher the score, the more likely this interaction is, according to the method. Marbach et al. [START_REF] Surowiecki | The wisdom of crowds: Why the many are smarter than the few and how collective wisdom shapes business[END_REF] ranked each interaction in the matrix R i .

(R i ) j,m = r (3.27) s.t. ∃( j b , m b ) b=1...r-1 (B i ) j b ,m b ≥ (B i ) j,m (3.28) ∃( j b , m b ) b=r+1...p×p (B i ) j b ,m b ≤ (B i ) j,m (3.29)
and then built the mean average over all network inference method,

R = 1 n J n J i=1 R i . (3.30)
To obtain an adjacency matrix, one can threshold the matrix R to obtain only a given percentage of most likely interactions. One hopes that kernel-based models infer well gene regulatory network, but also that this modelling method, as it is very different from others, will contribute well to enhance the GRN inference when used in a consensus method. For real-sized networks, I compared the performances of the consensus with and without the network inferred by the importance measures obtained in the previous subsection.

Performance

Results are shown in Table 3.3. While Jac does reach best performance for AU-ROC for two networks out of three, it shows very poor performances for AUPR. Jac obtains interesting performance on the difficult and interesting problem of GRN inference, but does not seem competitive with state-of-the-art methods as an individual method.

When observing the consensus, we observe a substantial increase in performance by adding Jac to other methods. Besides, Jac accounts for only 1 vote out of 30, so the performance increase is noteworthy. Marbach et al. have shown that consensus worked better if the methods gathered relied on different models, grasping different types of information. It seems that Jac, by its Gaussian kernel and partial derivatives approach, which are both novel, acquire an information that was not captured before.

Conclusion

In this work, I proposed a reverse-engineering GRN method not by learning an interpretable model, but by interpreting a Gaussian kernel-based model. The Jac importance measure has shown that it could deal with real-size networks, and provide good results. It seems that this approach captures some interactions which were overlooked by previously established methods.

For future work, one can use the Jac measure with regression methods that are more specific to gene modeling. One can add constraints on the functional cost to incorporate prior information, or change the regression model, for example operator-valued kernels, that learn structured output. Particular mention should be made of the work of N.Lim et al. [START_REF] Lim | Okvar-boost: a novel boosting algorithm to infer nonlinear dynamics and interactions in gene regulatory networks[END_REF] that shows very good performances on dynamic data, using partial derivatives of an operator-valued kernel model.

Part II Kernel Feature Weighting for Network Inference

In the previous part, we saw a method to interpret a kernel-based model. In the next part, I design an interpretable kernel model, a method that can model nonlinear behavior and can estimate which genes are regulating the target gene as directly as linear model. In addition, this method produce models whose gradient is sparse, thus the kernel model may be better suited to interpretation through partial derivatives than the one obtained with kernel-ridge regression. An additional advantage of this method is its ability to incorporate prior knowledge. Such a method is made possible by the versatility of kernels. Given M different kernels k 1 , . . . , k M , a positive linear combination of these kernels, k * = m w m k m , with w m ≥ 0 for all m = 1 . . . M, is a kernel. In particular, defining local kernel or component kernel as a kernel function that uses one feature of the input vector k m (x, z) = k(x m , z m ), the weights w m measure the importance of gene m for the target gene. Using kernel k * on the learning set S = (x 1 , . . . , x N ), the model will have the following form:

f (x) = N ℓ=1 α ℓ        M m=1 w m k m (x, x ℓ )        = M m=1 w m g m (x m )
The vector w will be learned from the training set, using multiple kernel learning. In particular, this additive model can be seen as an extension of the linear model to additive nonlinear dynamics.

The vector w can be seen as feature weighting. In Chapter 4, I describe several feature weighting or feature selection methods that can take into account nonlinear dynamics. On realistic datasets, the local kernel approach performs better that other feature selection methods. In Chapter 5, I compare my local kernel approach, called LocKNI (LOCal Kernels for Network Inference), to state-of-the-art methods on real and realistic datasets. LocKNI gives state-ofthe-art results. Besides, when used with other methods in a "wisdom of crowds" approach, LocKNI substantially enhances the performance of the consensus, hinting that local kernels grasp information that is not captured by other models. Finally, I provide a simple approach to incorporate prior knowledge to the multiple kernel learning approach, and use two reasonable biological sources of prior knowledge that greatly enhance LocKNI's performance on a realistic dataset.

Chapter 4

Local Kernels and feature selection

Introduction

In our era of data abundance, the scale of supervised learning problems has changed. Whereas few domains used more than 40 variables in 1997, today, many papers explore domains with hundreds to tens of thousands of variables. For example, text categorization [START_REF] Forman | An extensive empirical study of feature selection metrics for text classification[END_REF], image classification [START_REF] Tuia | Learning relevant image features with multiple-kernel classification[END_REF] or cancer prognostic from genomic data [START_REF] Sotiriou | Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis[END_REF][START_REF] Ein-Dor | Outcome signature genes in breast cancer: is there a unique set?[END_REF] deal with thousands of variables for a classification task. Reducing the problem size has many benefits: from a computational point of view, it reduces data storing requirements and computational time; from a human point of view, it helps understand and visualize the problem; from a pragmatic point of view, it helps defy the curse of dimensionality [START_REF] Bellman | Adaptive control processes: a guided tour[END_REF]. In this context, researchers have developed feature selection methods. Among all the variables, many may be redundant or noisy features, see [START_REF] John | Irrelevant features and the subset selection problem[END_REF][START_REF] Kohavi | Wrappers for feature subset selection[END_REF] for a formal definition of relevant features. In many cases, the best subset of features is the one that minimizes the generalization error. Feature selection methods can be decomposed in three categories [START_REF] Guyon | An introduction to variable and feature selection[END_REF]:

• Filtering methods. They constitute a preprocessing step, independent of the choice of the learning algorithm. One evaluates the dependence between the features and the value to predict through a pre-defined metric, such as mutual information or correlation. Generally, one uses this metric to rank features, and then select the most informative features. Filters are the fastest feature selection methods. Two classical filtering examples are Fisher's linear discriminant [START_REF] Bishop | Neural networks for pattern recognition[END_REF] and mutual information [START_REF] Zaffalon | Robust feature selection by mutual information distributions[END_REF] • Embedded methods. These machine learning methods simultaneously select relevant features and build the prediction function. While slower than filters, these methods select the best features for a specific model. LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and decision trees, as viewed in CART [START_REF] Breiman | Classification and regression trees belmont[END_REF], are examples of such methods

• Wrapper methods. These methods take into account the model of prediction functions and are the most general, but also the slowest methods. They consist in selecting a subset of features F, learning a prediction function using F, and then training a prediction function on other subset of features F ′ until a subset F * is found that is optimal with respect to a criterion, such as accuracy of the prediction function. Sequential Forward Selection [START_REF] Kohavi | Wrappers for feature subset selection[END_REF] is a wrapper method.

Feature weighting -where features is given a weight from 0 to 1-can be seen as an extension of feature selection -where feature are either useful (weight is 1) or ignored (weight is 0)-. The field of gene regulatory network (GRN) inference from gene expression data is keen of feature selection or feature weighting methods. Using these methods, one can evaluate which genes are relevant for the prediction of a target gene, and consequently assume that they regulate this target gene. Filtering methods [START_REF] Margolin | Aracne: An algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context[END_REF][START_REF] Faith | Large-scale mapping and validation of escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF] and linear [START_REF] Haury | Tigress: Trustful inference of gene regulation using stability selection[END_REF][START_REF] Gustafsson | Constructing and analyzing a large-scale gene-to-gene regulatory network-lasso-constrained inference and biological validation[END_REF], Boolean [START_REF] Shmulevich | Binary analysis and optimization-based normalization of gene expression data[END_REF][START_REF] Akutsu | Identification of genetic networks from a small number of gene expression patterns under the boolean network model[END_REF], Bayesian [START_REF] Perrin | Gene networks inference using dynamic bayesian networks[END_REF] or tree-based [START_REF] Huynh-Thu | Inferring regulatory networks from expression data using tree-based methods[END_REF]

Feature weighting methods

In order to infer a GRN from gene expression data, given training data S = (x 1 , . . . , x N ), N observations of the expression level of a p-gene system, one can learn a static or dynamic model. If each gene is modeled independently, a feature weighting/selection algorithm may be used to weight the importance of other genes for the prediction of a target gene. One can learn a static model

x i = f i stat (x -i ) + ǫ (4.1) a dynamic model x i (t + τ) = f i dyn,τ (x(t)) + ǫ (4.2)
Using a feature weighting (resp. selection) algorithm in learning f i • , one obtains the vector of weights w i ∈ R p + (resp. w i ∈ {0, 1} p ) that evaluates the importance of each gene for the prediction of gene i. The estimated adjacency matrix is the juxtaposition of these feature weights:

 = w 1 , . . . , w p T (4.3)

Local Kernel Approach

In order to have a nonlinear interpretable model, I suggested to use local kernels or component kernels, k 1 , . . . , k M , kernel k m relying only on feature m of input variable. Then, a weighted combination k * of those kernels is used as a kernel; k * = M m=1 w m k m , with weights w learned from the training set. This produces a model of the form:

f (x) = N ℓ=1 α ℓ        M m=1 w m k m (x, x ℓ )        , (4.4) = M m=1 w m g m (x m ) , (4.5)
with a sparse vector w. This additive model has several advantages: firstly, it is able to model nonlinear dynamics in order to reflect the true behavior of biological systems. Secondly, by controling the sparsity of w, the sparsity of the partial derivatives of f is controlled (∂ f /∂x m = 0 if w m = 0). Thirdly, extraction of regulators is as direct as in linear models, vector w quantifying each feature's importance. This model can be seen as an extension of linear models. The vector w is learned from data, using Multiple kernel learning (MKL). MKL was introduced in [START_REF] Lanckriet | Kernel-based data fusion and its application to protein function prediction in yeast[END_REF] for protein function prediction. MKL has improved stateof-the-art classification performances in various domains, such as bioinformatics [START_REF] Lanckriet | Kernel-based data fusion and its application to protein function prediction in yeast[END_REF][START_REF] Sonnenburg | Learning interpretable svms for biological sequence classification[END_REF], image classification [START_REF] Guo | Customizing kernel functions for svm-based hyperspectral image classification. Image Processing[END_REF][START_REF] Tuia | Learning relevant image features with multiple-kernel classification[END_REF], or sound localization [START_REF] Takashima | Feature selection based on multiple kernel learning for single-channel sound source localization using the acoustic transfer function[END_REF]. For a review of MKL algorithm, I recommend [START_REF] Gonen | Multiple kernel learning algorithms[END_REF]. In this chapter, I will use three MKL algorithms: SimpleMKL [126], ℓ 2 -MKL, developed in this thesis, and SUpport vector Parsimonious ANOVA (SUPANOVA) [START_REF] Gunn | Structural modelling with sparse kernels[END_REF].

SimpleMKL

In this algorithm, one obtains a linear combination of kernels by minimizing a functional cost with an ℓ 1 -constraint on the kernel weights w: w m ≥ 0 for all m and m w m = 1. These constraints produce sparse weights w, thus giving interpretable feature selection.

To define the functional cost, we notice that each kernel k m has a reproducing kernel Hilbert space H m . For any function f ∈ H m , one has < f, k m (x, .) > H m = f (x).

In particular, for any function of the following type

f (x) = i α i k m (x i , x), one has f 2 H m = i, j α i , α j k m (x i , x j ) = α T K m α. A.
Rakotomamonjy [126] considers the following optimization cost:

L( f 1 , . . . , f M , w) = i l(y i , M m=1 f m (x i )) + λ M m=1 1 w m f m 2 H m (4.6) s.t. w m ≥ 0 ∀m, m w m = 1 (4.7)
If l(., .) is the hinge loss or the ǫ-insensitive loss, Rakotomamonjy et al. prove that solution fm has the form fm (x) = i α i w m k m (x i , x). With convention 0 0 = 0, the functional cost is define for all w ∈ [0, 1] M . Rakotomamonjy et al. proves that the functional cost (4.6) is convex, so it has a unique minimum, which is global. They reformulate the functional cost: Existence and computation of derivatives of J(.) are done in [126]. They obtain the minimizer of J(.) by a reduced gradient method, which converges for such function [START_REF] Luenberger | Linear and nonlinear programming[END_REF]. Once ∇J is calculated, they use a descent direction ∇ red J assuring that w+ s∇ red J still satisfies the constraints (4.9). They perform line-search to find the optimal descent-step. These reduced gradient-descent will give ŵ, the minimizer of J(). 2 . I constrain f m to have the form f m (.) = i α i w m k(x i , .). I minimize the following functional cost:

(x i ) = (y i -m f m (x i ))
L( f 1 , . . . , f M , w) = i (y i - M m=1 f m (x i )) 2 + λ M m=1 1 w m f m 2 H m (4.10) s.t. f m (x) = i α i w m k m (x i , x), ∀m (
α =        m K m + λId        -1 y (4.13)
For (b), I do a reduced-gradient descent. The functional cost L( f 1 , . . . , f M , w) decreases at each iteration and is lowerly bounded by 0, thus it converges to a local minimum.

SUPANOVA In this article [START_REF] Gunn | Structural modelling with sparse kernels[END_REF], sparsity of w is not assured by an ℓ 1constraint, but by an ℓ 1 -regularization term. S.Gunn et al. minimize the following functional cost:

L(α, w) = i         y i - j α j        m w m k m (x j , x i )                2 (4.14) +λ 1 α T        m w m K m        α + λ 2 w 1 (4.15) s.t. w m ≥ 0 ∀m (4.16)
They kept the positivity constraint on w, so m w m k m (., .) is still a semi-definite positive function. Functional cost (4.14) is a quadratic problem in α and in w. They alternatively solve it for α with w fixed and for w with α fixed. For optimal α, they have closed-form formula:

α =        m w m K m + λ 1 Id        -1 y (4.17)
For w, they use an optimization algorithm for quadratic functions under constrains [START_REF] Coleman | A reflective newton method for minimizing a quadratic function subject to bounds on some of the variables[END_REF]. S.Gunn and J.S. Kandola [START_REF] Gunn | Structural modelling with sparse kernels[END_REF] proves convergence to a global minimum by this argument : "In the quadratic case the second order partial derivatives with respect to α and w are always positive ensuring that every slice is convex. This fact combined with the knowledge that the solution is finite in α and w should ensure convergence to the global minimum".

Filtering methods

Many filtering methods are based on correlation and mutual information, see [START_REF] Duch | Filter methods[END_REF][START_REF] Sánchez-Maroño | Filter methods for feature selection-a comparative study[END_REF] for comparative study of filtering methods. I will not experiment with correlation or mutual information-based filters, as they have been used and improved a lot for GRN inference, see section 1.4.1, page 26. In the following, I will test kernel target alignment, a feature selection method which tries to find optimal feature weights for the classification or regression task. This method should take into account multivariate dependencies. I also implement RReliefF, a filter capable to identify nonlinear dependency, and that was not tested for GRN inference, to my knowledge.

Kernel target alignment (KA)

This measure was first published in [START_REF] Cristianini | On kernel target alignment[END_REF]. A classification or regression problem would be easy to solve if the input variable x i were only the label to predict y i = x i for all i from 1 to N. In this case, the perfect predictor would simply be the identity f * (x i ) = f * (y i ) = y i . Using a linear kernel, one would have the perfect Gram-matrix K * linear = yy T whose entries are K * linear i j =< y i , y j >. Cristianini et al. look for the best feature weighting to make the Gram-matrix K w close to the perfect Gram-matrix. With feature weights w, the Gram-matrix K w take the values (K w ) i j = k(w • x i , w • x j ), where • denotes the Hadamard product, w • x i = w 1 x 1 i , . . . , w p x p i . Noting < ., . > F the Euclidian scalar product between two matrices of size N × N:

< A, B > F = N i, j=1 A i j B i j , (4.18)
Cristianini defines the alignment between two matrices, K and K * , as:

A(K, K * linear ) = < K, K * linear > F K F K * linear F , (4.19)
which can be interpreted as the cosine between the two matrices. If the employed Gram-matrix is aligned with the perfect Gram-matrix, the prediction function should be easy to find. For the classification case, Cristianini, in [START_REF] Cristianini | On kernel target alignment[END_REF], shows that the generalization error is bounded by 1 -A(K, K * linear ) and that is has the concentration property 1 . Kernel alignment has been used for regression, in [START_REF] Kandola | On the extensions of kernel alignment[END_REF], with centered output y ′ i = y iȳ, with ȳ the mean of y. In the regression case, kernel alignment still has the concentration property. In my experiments, I shall use a weighted Gaussian kernel:

K w (x, z) = exp        - p i=1 w i (x i -z i ) 2        (4.20)
and optimize w in order to aligne K w with K * linear . I also suggest to align K w with K * Gauss , the Gaussian kernel that would be obtained with perfect input:

K * Gauss i j = exp (y i -y j ) 2 2 (4.21)
I follow the optimization procedure in [START_REF] Ramona | Multiclass feature selection with kernel gram-matrix-based criteria[END_REF], called Scaled Alignment method. Starting from w 0 = [1, . . . 1] T , I perform a gradient descent until I have found a local minimum, see Algorithm 1.

RReliefF

This filter is employed to measure the quality of an attribrute for the task of prediction. First developed for binary classification [START_REF] Kira | A practical approach to feature selection[END_REF], "Relief" randomly selects an instance r of the training set S. It looks at the nearest instance of the same class (nearest hit, h) and of the opposite class (nearest miss, m). It evaluates the distance between two instances a and b according to attribute i with the function dist:

Algorithm 1 Scaled Alignment Selection

Starting values w 0 = (1, . . . , 1)

T , A 0 = 0, n = 0 Choose K * = K * linear/Gaussian while |A n -A n-1 | ≥ ǫ do Compute gradient g n+1 = ∇ w n A(K w , K * ) Obtain descent-step by line search η n+1 = arg max η∈[0,1] A(K w n +ηg n+1 , K * ) w n+1 = w n + η n+1 g n+1 A n+1 = A(K w n+1 , K * ) n ← n + 1 end while dist(i, a, b) = 0 if a i = b i 1 otherwise if attribute i is discrete, (4.22) dist(i, a, b) = |a i -b I | max x∈S x i -min x∈S x i if attribute i is continuous. (4.23)
A good attribute will have a large distance between the selected instance r and its nearest miss, and a small distance with its nearest hit. To measure the quality of each attribute in the vector w, Kim et al. initialize w at 0 p , randomly select n instances in the training set and for each selected instance r with nearest hit h and nearest miss m, increment w according to:

w i ← w i + 1 n (dist(i, r, m) -dist(i, r, h)) (4.24)
This algorithm is given in Algorithm 2. An improvement of "Relief" by selecting the k nearest hit and miss, called "ReliefF" [START_REF] Kononenko | Estimating attributes: analysis and extensions of relief[END_REF],is more robust to noise and can be applied to multiclass problems. Furthermore, it can be interpreted probabilistically. Noting y r the label of instance r, if all attributes are discrete, "ReliefF" is an approximation of the difference between probabilities:

w i = P(dist(i, r, m)|y r y m ) -P(dist(i, r, h)|y r = y h ) (4.25)
Noting the probability of label difference P di f f Y = P(y r y m ), the probability of attribute distance P di f f i = P(dist(i, r, m)) and the probability of label difference knowing attribute distance P di f f Y|di f f i = P(y r y m |dist(i, r, m)), equation (4.25) may be extended following Bayes's rule according to:

w i = P di f f Y|di f f i P di f f i P di f f Y - 1 -P di f f Y|di f f i P di f f i 1 -P di f f Y (4.26)
Inspired by equation (4.26), "ReliefF" has been extended to regression problems under the name "RReliefF" [START_REF] Robnik | An adaptation of relief for attribute estimation in regression[END_REF], . In a regression problem, the notions of nearest hit or miss do not exist. Given a selected instance r 0 with label y 0 and its k nearest neighbors r 1 , . . . , r k with labels y 1 , . . . , y k , (a) P di f f Y evaluates how much the labels vary (b) P di f f i how much each attribute varies and (c) P di f f Y|di f f i how much the labels vary with a variation of the attribute. Robnik-Sikonja et al. [START_REF] Robnik | An adaptation of relief for attribute estimation in regression[END_REF] suggest to weight the influence of the j th nearest neighbor with function I:

I( j) = I 1 ( j) k ℓ=0 I 1 (ℓ)
where I 1 is a method to weight the influence, (4.27)

I 1 ( j) = 1 uniform weight, (4.28) I 1 ( j) = exp -j σ exponentially decreasing weight, (4.29)
Similarly to "ReliefF", they randomly select n instances of the training set. They set P di f f Y = 0, and, for all i, P di f f i = 0 and P di f f Y|di f f i = 0. Then, for each selected instance r 0 , they identify the k nearest neighbors r 1 , . . . , r k and perform the following updates:

P di f f Y ← P di f f Y + k j=1 I( j)|y r 0 -y r j | (4.30)
P di f f i ← P di f f i + j I( j)dist(i, r 0 , r j ) (4.31)
P di f f Y|di f f i ← P di f f Y|di f f i + j I( j)|y r 0 -y r j |dist(i, r 0 , r j ) (4.32)
Once these computations have terminated, they evaluate the weights of each attribute according to:

w i = P di f f Y|di f f i P di f f Y - P di f f i -P di f f Y|di f f i n -P di f f Y . (4.33)
This algorithm is described in Algorithm 3. For a study of theoretical and empirical performances of "RReliefF", see [START_REF] Robnik-Šikonja | Theoretical and empirical analysis of relieff and rrelieff[END_REF].

Algorithm 2 Relief algorithm

Given learning set S Starting w = 0 p for ℓ = 1 to n do Randomly select instance r ∈ S Find nearest hit h and nearest miss m in S for i = 1 to p do

w i = w i + dist(i, r, m) -dist(i, r, h) end for end for Output w Algorithm 3 RReliefF algorithm Given learning set S Starting P di f f Y = 0, P di f f i = 0, P di f f Y|di f f i = 0 for ℓ = 1 to n do Randomly select instance r 0 ∈ S Find k nearest instances r 1 , . . . , r k in S P di f f Y = P di f f Y + k j=1 I( j)|y r 0 -y r j | for i = 1 to p do P di f f i = P di f f i + j I( j)dist(i, r 0 , r j ) P di f f Y|di f f i = P di f f Y|di f f i + j I( j)dist(i, r 0 , r j )|y r 0 -y r j | end for end for w i = P di f f Y|di f f i P di f f Y - P di f f i -P di f f Y|di f f i n-P di f f Y Output w

Embedded methods

Recursive Feature Elimination

One feature selection method iteratively eliminates the least informative feature [START_REF] Kohavi | Wrappers for feature subset selection[END_REF]. Starting from F 0 , the set of all features, one creates nested subsets of features

F 0 ⊃ F 1 ⊃ • • • ⊃ F final until F final
is an optimal subset of features. To identify the least informative feature, one can learn, for each feature i, a predictor ignoring this feature i, but this has a high computational cost. To avoid running the learning algorithm many times, Le Cun et al. [START_REF] Lecun | Optimal brain damage[END_REF] suggest the Optimal Brain Damage algorithm (OBD). A prediction function f α, characterized by parameter α, is learned by minimizing a cost function J( α). To learn the importance of each feature i, Le Cun et al. calculate the change in objective function J for a change in parameter α i = αi + h i :

J( α + h) = J( α) + i h i ∂J ∂α i + i h 2 i 2 ∂ 2 J ∂α 2 i + j i h i h j 2 ∂ 2 J ∂α i ∂α j + O( h 3 ), (4.34)
They suggest the following simplifications:

• the "extremal" approximation: they assume α is an optimum of J, therefore ∂J ∂α i = 0 • the "diagonal" approximation: they assume that a change in J by deleting several parameters is close to the sum of changes caused by individually deleting parameters. Therefore, they neglect cross-terms ∂ 2 J ∂α i ∂α j if j i • The "quadratic" approximation: the functional cost is nearly quadratic, so they neglect the last term of the equation Equation (4.34) thus reduces to:

J( α + h) -J( α) ≈ i h 2 i 2 ∂ 2 J ∂α 2 i . (4.35)
A change h i = -α i corresponds to removing feature i. They denote by:

D i J = α2 i 2 ∂ 2 J ∂α 2 i (4.36)
the importance of feature i. Guyon et al. [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] apply this methodology to SVM with SVM-RFE, meaning SVM Recursive Feature Elimination. In particular, for linear SVM, they obtain a prediction function f α (x) =< α, x > +b and D i J = α 2 i . Starting from the set of all feature F, they train a linear SVM. They eliminate the n F features i with smallest parameter α i , n F being a number of features. They continue (a) learning a linear SVM on the remaining features and (b) removing the least informative features until they have found an optimal subset of features. The corresponding algorithm is described in Algorithm 4. I present later in this section ways to determine the optimal subset of features using several criteria. For non-linear kernel methods, I.Guyon et al. compute K F , the Gram-matrix on all features, and they have a prediction function of the form:

f α,F (x) = n i=1 α i k F (x i , x) (4.37)
with k F the kernel function calculated with all features of set F. With α fixed, they calculate the cost function J if feature i was removed, with kernel k F\{i} . If feature i is not important, the functional cost J will have little change. The algorithm is described in Algorithm 5.

In this chapter, I will use RFE with kernel-ridge regression using linear and Gaussian kernel.

i R = [F( f ), R] F = [F(1 : f -1), F( f + 1 : end)]
Find least informative feature f = arg min i J(K F\{i} , w) R = [F( f ), R] F = [F(1 : f -1), F( f + 1 : end)] end while Output: feature ranked list R

KerNel Iterative Feature Extraction -KNIFE

This algorithm performs feature selection by weighting each feature in the kernel. One chooses a kernel function k. Two instances x, z are compared with the weighted kernel k w (x, z) = k(w • x, w • z), where (w • x) i = w i x i . This feature selection idea was described in [START_REF] Weston | Feature selection for svms[END_REF][START_REF] Grandvalet | Adaptive scaling for feature selection in svms[END_REF]. In the KNIFE algorithm [START_REF] Allen | Automatic feature selection via weighted kernels and regularization[END_REF], the following cost is minimized:

min α,w L(α, w) = L(y, K w α) loss function + λ 1 α T K w α α regularization + λ 2 w 1 w regularization (4.38)
subject to 0 ≤ w j ≤ 1, for all j = 1 . . . p .

The ℓ 1 -regularization produces sparsity in the feature-weights w. The functional is minimized by alternatively optimizing α with w fixed and optimizing w with α fixed. Optimal α is found via standard SVM algorithms. To optimize w, Allen et al. [START_REF] Allen | Automatic feature selection via weighted kernels and regularization[END_REF] consider the linear approximation of the Grammatrix around current feature-weight estimation w (t) , Kij = kw (x i , x j ) where kw (x, x ′ ) = k w (t) (x, x ′ )+ < ∇ w k w (t) (x, x ′ ), ww (t) >. They minimize a convex surrogate cost function:

min w L(α, w) = L(y, Kw , α) + λ 1 α T Kw α + λ 2 w 1 1 (4.40) s.t. 0 ≤ w j ≤ 1, for all j = 1 . . . p (4.41)
Allen proves that this algorithm converges to a local minimum for several loss functions: the hinge loss (support vector machine), squared error loss (kernel-ridge regression) and binomial deviance loss (logistic regression). I shall use this algorithm with the squared error loss and the Gaussian kernel. Note that, since the algorithm scales the inputs with w, the kernel bandwidth σ is no longer an hyper-parameter. There remains only two hyper-parameters, the regularization parameters λ 1 and λ 2 .

Algorithm 6 KNIFE algorithm

Input: training examples X = [x 1 , . . . , x n ] with labels y = [y 1 , . . . , y n ] Randomly initialize weights 0 ≤ w (0) < 1, set t = 0 while w (t)w (t-1) 1 ≥ ǫ do Compute, with a SVM algorithm, α (t) = arg min α L(y, K w (t) ) + λ 1 α T K w (t) α t = t + 1 Compute, by gradient descent, w (t) = arg min w L(α, w) under constraints 0 ≤ w ≤ 1 end while Output: weight of all features w

Hyper-parameter selection method

The described feature selection methods often require some hyper-parameters, be it the Gaussian kernel bandwidth σ or a regularization parameter λ 1 , λ 2 , ǫ. I present several methods to determine hyper-parameter values such that the trained prediction function would generalize well. Then, I will not verify that the prediction function generalize well, I check if the feature selection method selected the genes that regulate a target gene. Given a feature selection method with c hyper-parameters, I define a finite set of values Λ i that hyper-parameter i can take, noting µ ∈ (Λ 1 × • • • × Λ c ) a value of hyper-parameters for the feature selection method.

Cross-Validation (CV)

Given a value µ for the hyper-parameters, the goal is not to measure the quality of the prediction function on the training set but on new unlabeled data, usually by means of the generalization error. Cross-validation is one of the most popular methods to evaluate the generalization error [START_REF] Duan | Evaluation of simple performance measures for tuning svm hyperparameters[END_REF]. To estimate the generalization error, I perform k-fold cross-validation. The training data S is partitionned in k equal size sets, the subsamples E 1 , . . . , E k . One subsample E i is kepts as the validation data to test. Denoting f -i to be the prediction function trained on the k -1 remaining subsamples, the mean square error (MS E) is evaluated according to:

MS E = 1 |E i | (x,y)∈E i (y -f -i ) 2 . (4.42)
For each choice µ of hyperparamaters, the MS E is obtained on each subsample, using a function trained with hyper-parameter µ on the other subsamples. The value of µ for which the MS E is the lowest is selected for the hyper-parameter.

Stability

Breiman [START_REF] Breiman | Heuristics of instability and stabilization in model selection[END_REF] points out that, since cross-validation selects a hyper-parameter value on the training data, the MS E may be greatly underestimated by this procedure when the dataset is not large enough. He suggests aiming for a stable model if we have few data -which is typically the case in gene regulation network inference. If a small variation in the training data leads to significant changes in the prediction function, the model is unstable, hence it is not reliable. To evaluate the stability of a method with hyperparameter value µ, I follow the procedure in [START_REF] Ben-Hur | A stability based method for discovering structure in clustered data[END_REF]. I create n s = 50 subsamples of the training data S:E 1 , . . . , E n s . Each subsample uses a fraction f = 80% of the training data. Denoting w i be the feature weights of the model trained on subsample E i . The stability of an algorithm with hyper-parameter µ is evaluated a the mean alignment of feature weights, whose definition is in equation (4.19):

stab = 2 n s (n s -1) n s i=1 n s j=i+1 A(w i , w j ) (4.43)
I select the value of the hyperparameter µ for which the feature selection method has the highest stability.

Block-Stability

Politis et al. [START_REF] Politis | Automatic block-length selection for the dependent bootstrap[END_REF][START_REF] Patton | Correction to automatic block-length selection for the dependent bootstrap by d. politis and h. white[END_REF] observe that this stability evaluation treats data from training set as if they were independent. When the data is a time-series, the data samples are not independent. To get around this limitation, they introduced a block-stability procedure, that takes into account the dependence of data. After having selected a block-size b and randomly chosen an integer r ∈ [1, n], I build subsample E i = (x r , x r+1 , . . . , x (r+b)mod n ) of the training set S = (x 1 , . . . , x n ). Given subsamples E 1 , . . . , E n s built with the block-stability procedure, stability is then calculated as in the previous method.

Bayesian and Akaike Information Criterion

The fewer parameters a model has, the less likely it is to overfit the data. From this observation, Schwarz [START_REF] Schwarz | Estimating the dimension of a model[END_REF] developed a criterion for selecting hyperparameters, consisting of the likelihood function and a penalty term for the number of parameters in the model. Assuming that the model errors are independent and distributed according to a centered normal distribution, and denoting:

σ2 E = 1 n -1 N i=1 (y i -f (x i )) 2 (4.44)
the empirical variance of the model errors, he measures the likelihood of observing y 1 , . . . , y N given the model f (x 1 , . . . , x N ) by:

p(y| f ) = N i=1 1 2π σ2 E exp (y i -f (x i )) 2 2 σ2 E (4.45)
He looks for a model that maximizes the probability of observations y under model f while minimizing the number of free parameters. In my case, the number of free parameters is the number of features used in the model. Denoting k the number of used features,, i.e. features for which w i > 0, I aim to minimize the Bayesian Information Criterion:

BIC = -2 log (p(y| f )) + k log(N) (4.46) = N log( σ2 E ) + k log(N) + Constant (4.47)
The Akaike Information Criterion [START_REF] Akaike | A new look at the statistical model identification[END_REF] relies on the same idea of maximizing the log likelihood of observed data while minimizing the number of free parameters, and is given by this formula:

AIC = 2k -2 log(p(y| f )) (4.48) = 2k + 2N log( σE ) + Constant (4.49)

Experiments

The experiments are carried out using the DREAM3 Challenge 4 data, described in Chapter 3. Performance will be evaluated for AUROC and AUPR, defined in Chapter 3.

Hyperparameter selection

The methods have been tested with hyper-parameters in the following sets:

• For SimpleMKL, regularization tradeoff C ∈ {0.01, 0.1, 1, 10, 100}, ǫinsensitive loss ǫ ∈ {0.01, 0.05, 0.15} and kernel bandwidth σ ∈ {0.5, 1, 2}.

• For ℓ 2 -MKL, regularization tradeoff λ ∈ {0.01, 0.1, 1, 10, 100} and kernel bandwidth σ ∈ {0.5, 1, 2}.

• For SUPANOVA, α regularization tradeoff λ 1 ∈ {0.01, 0.1, 1, 10, 100}, w regularization tradeoff λ 2 ∈ {0.01, 0.1, 1, 10, 100} and kernel bandwidth σ ∈ {0.5, 1, 2}.

• For KNIFE, α regularization tradeoff λ 1 ∈ {0.01, 0.1, 1, 10, 100} and w regularization tradeoff λ 2 ∈ {0.01, 0.1, 1, 10, 100}.

• For kernel target alignment and RReliefF, there is no hyper-parameter to set.

• For RFE, regularization tradeoff λ ∈ {0.01, 0.1, 1, 10, 100} and the number of used features to model target gene k ∈ [1, . . . , p]. For the number of used features, stability is irrelevant, it will always use all the features. Thus the selection of hyper-parameters must be performed through cross-validation, BIC or AIC for RFE methods.

The performance of a feature selection method with a hyper-parameter selection criterion greatly varies with the network studied, see The AUROC appears to depend more on the dataset than the hyperparameters. Second, the shape of the distribution of the points varies greatly from one network to another. Cross-validation does not seem to be a good hyperparameter selection method when SimpleMKL is used for the task of network inference.

between the performance of the feature selection method and the value of the hyper-parameter criterion. The average of this correlation, across all networks, is shown in Tables 4

Performance on size 50 networks

AUROC and AUPR on all networks are displayed in Figure 4.3. Also shown in this figure are the two best kernel feature selection methods against the two best competitors of the DREAM3 challenge. We see that ℓ 2 -MKL is the best performer in AUROC on all five 50-gene networks. In AUPR, the Gaussian kernel alignment prevails. Performance is still poor, and computation time will be a burden for ℓ 2 -MKL on size 100 networks. As we shall see in Chapter 5, ensemble methods will enable the use of ℓ 2 -MKL for networks of size 100 or more.

Conclusion

Gene regulatory network inference from gene expression data is an active research field, that will increasingly benefit from the abundance of data in the coming years. Many approaches to GRNI have been explored, but kernel feature selection methods has so far been neglected. In this work, I tested nine kernel feature selection methods, and five hyper-parameter selection methods on realistic datasets. The results show the limitations of these methods. No hyper-parameter selection method seems efficient, except BIC for RFE methods, which would need further testing to be validated. Some methods are too time-consuming to be directly applicable to real-size networks. This comparative study also gave some positive indications for constructing methods that would give better results. ℓ 2 -MKL and Gaussian kernel alignment appeared to extract more information for inference of GRN. This motivated the enhancement of the ℓ 2 -MKL method, in particular its enhancement by ensemble methods. Learning the network on many subsamples of data usually yields better result than learning one network on all data [START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF]. In the following chapter, I will show that an enhanced ℓ 2 -MKL method yields state-of-the-art results on real networks.

Chapter 5

LocKNI: Local Kernel for Network Inference

Introduction

Local kernels were motivated by several reasons. First, kernel methods deal with several difficulties of gene modeling, since they are robust-to-noise, they provide a nonlinear model, and the versatility of the loss function allows incorporation of prior knowledge. We saw, in the previous chapter, that ℓ 2 -MKL using local kernels provided some of the most interesting results in GRN inference among kernel feature weighting or feature selection methods, therefore, this method will be implemented with the ensemble method described in Chapter 3. Even on small scale networks, ensemble methods will increase the performance of the GRNI method. I shall call LocKNI the alliance of ℓ 2 -MKL, local kernels and the double randomization scheme. Original contributions can be made when using local kernels in the regular ℓ 2 -MKL framework. A simple modification of the loss function introduces prior knowledge in the regularization framework, thus improving LocKNI's results. Besides, since they rely on a few features, local kernel-based models have sparse gradients, thus have more interpretable partial derivatives. The drawback of these local models is that they are no longer universal, hence they may be unable to approximate well the partial derivatives of a perfect model. In this chapter, I first describe ℓ 2 -MKL and its modification to incorporate prior knowledge. Second, as in Chapter 3, I use randomized ℓ 2 -MKL models with local kernel in an ensemble, exploiting a double scheme of randomization both on variables and on data. In the fourth section, I explain the chosen hyperparameter selection method. In the fifth section, I compare, on real and real-sized networks, LocKNI without prior knowledge with state-of-the-art methods. I also evaluate the partial derivatives of LocKNI as a GRNI method. There will also be a discussion of LocKNI's strength and weaknesses: how does it behave compared to other GRNI methods ? What types of error does LocKNI make ? In the fifth section, I propose two reasonable priors that may be incorporated in LocKNI. I will test, on simulated data, the improvement in LocKNI's performance with incorporation of prior knowledge.

ℓ 2 -MKL and Prior Knowledge Incorporation

ℓ 2 -MKL

We saw, in the previous chapters, that two types of training data S = (x 1 , . . . , x p ) may be available, to learn either a static or a dynamic model:

x i = f i stat (x -i ) + ǫ i static model, (5.1) x i (t + τ) = f i dyn,τ (x -i (t)) + ǫ i (t) dynamic model. (5.2)
These models can be estimated with supervised learning tools on a training set.

With the general notation S = ((z 1 , y 1 ), . . . , (z N , y N )), z being input variables, y being an output variable to predict, a function y = f (z) + ǫ is learned on the training set. This notation will be used to describe the learning algorithm.

In the previous chapter, a multiple kernel approach was presented, using local kernels, i.e. kernels that use only one component of an input vector, k m (z 1 , z 2 ) = k(z m 1 , z m 2 ). We saw how a linear combination of these kernels, k * = m w m k m , could be learned, by optimizing the weight w. The function f and weights w are found by minimizing the following functional cost:

L(w, α) = N i=1 (y i - M m=1 f m (z i )) 2 + λ M m=1 1 w m f m 2 H m (5.3) s.t. (C) : M m=1 w m = 1, w m ≥ 0, ∀m ∈ {1, . . . , M} (5.4) f m (z) = i α i w m k m (z i , z) . (5.5)
This function can be optimized by alternately optimizing L(w, α f ixed ) over w when α is fixed and optimizing L(w f ixed , α) over α when w is fixed. With w fixed, the minimum is obtained in closed-form as:

α = (K * + λI N ) -1 y (5.6) with K * = m w m K m (5.7)
and I N is the N × N identity matrix. With α fixed, w is optimized through reduced gradient descent [START_REF] Luenberger | Linear and nonlinear programming[END_REF].

Prior knowledge incorporation

As kernel-based models result from the minimization of a unique global loss function, they benefit from the possibility of taking prior knowledge into account in the definition of the loss function. In order to incorporate an assumption of the existence of a regulator, I propose to relax the sparsity constraint imposed on the weight w m by dividing λ by a coefficient λ m . If prior knowledge hints that gene m should be a regulator gene i, λ m would be set larger than 1; if prior suggests that m should not be a regulator, λ m < 1:

L(w, α) = N ℓ=1 (y ℓ - M m=1 f m (z ℓ )) 2 + M m=1 λ λ m w m f m 2 H m (5.
f m (z) = i α i w m k m (z i , z) (5.10)
As in Section 5.2.1, this function can be optimized by alternately optimizing L(w, α f ixed ) over w with α fixed and optimizing L(w f ixed , α) over α with w fixed. With w fixed, the minimum can be expressed in closed-form as: With α fixed, w is optimized through reduced gradient descent.

α = K 2 * + λK prior * -1 (K * y) (5

Ensemble method

We saw, in Chapter 3, an ensemble method to learn the GRN from several heterogeneous datasets. From available datasets (E 0 , E 1 , . . . , E N e ), I build B subsamples E ℓ , l = 1 . . . B, with subbagging and subsampling of variables, in a way similar to extreme-randomization [START_REF] Haury | Feature selection for gene expression data: molecular signatures for breast cancer outcome prediction and gene regulatory network inference[END_REF]. Subsamples are built according to the following procedure: E u , with p data a fixed hyper-parameter. 3. Randomly select n var variables to be the potential regulators of the system.

The number n var is a hyper-parameter to fix. The set of selected variables is called G ℓ .

To learn the importance of genes in G ℓ for all the other genes i on subsample E ℓ , the mean vote for a regulation of gene i by gene m is computed as the mean weight w i m of feature m when predicting gene i with ℓ 2 -MKL and local kernels. Mean votes are stored in matrix B:

B im = 1 n im         ℓ,m∈G ℓ ,i G ell w i m         , (5.14)
with n im the number of subsamples where gene m was a potential regulator and gene i was not. The estimated adjacency matrix can then by obtained by thresholding matrix B:

Âim = H(B im -θ) (5.15)
where H is the Heaviside step function.

Experimental results

Hyperparameter choice

The chosen kernel function is the Gaussian kernel

k(z, z ′ ) = exp (z -z ′ ) 2 2σ 2 , (5.16)
which has the advantage of being widely used, with good performances in practice and universal consistency. Four hyperparameters remain to be set:

• I select the kernel bandwidth σ = 1, as I have normalized the empirical variance of the data to 1, and the Gram matrix achieves maximum entropy when σ matches the empirical standard deviation.

• Experiments have shown that the subsampling parameters, n var and p data , have little influence on the capacity of LocKNI to recognize edges. Since LocKNI is faster with smaller subsamples, I set n var = 5 and p data = 0.20.

• For the tradeoff λ between regularization and data fitting, I also use a heuristic. With a sample D of size N D , it can be shown that kernel-ridge regression is consistent if λ = 1/ √ N D [START_REF] Caponnetto | Optimal rates for the regularized leastsquares algorithm[END_REF].

We can consider that convergence has been reached with L = 100×p n var subsamples. For each regression, computational complexity is O(|E l | 3 ) for each variable of each subsample. There are L subsamples and p variables, computational complexity is O(Lp(p data N) 3 ).

Datasets

The evaluation of LocKNI will be two-fold. First, it is evaluated without integration of prior knowledge. LocKNI will be compared to state-of-the-art methods on four real and real-sized networks. The datasets are summarized in Table 5.1. The Jac measure will be used on LocKNI's model as a second approach to GRNI. LocKNI's strength and weaknesses will be studied on simulated data. LocKNI's error will be categorized. Then, LocKNI is evaluated with prior information on a simulated dataset. The purpose of this assessment is to demonstrate the benefit of incorporation of prior knowledge to the reverse-engineering approach. For evaluation against state-of-the-art methods, the first three datasets come from the DREAM5 challenge [START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF], as described in Chapter 3 and the fourth from the Many Microbe Microarray Database (M3D) [START_REF] Faith | Many microbe microarrays database: uniformly normalized affymetrix compendia with structured experimental metadata[END_REF] (EColi version 4 build 6). For EColi from M3D, the gold standard is the currently known network available in RegulonDB v8.1 [START_REF] Salgado | Regulondb v8. 0: omics data sets, evolutionary conservation, regulatory phrases, cross-validated gold standards and more[END_REF].

For the DREAM5 datasets the network inferred by LocKNI is compared to the best networks inferred in the challenge. For the EColi dataset from M3D, LocKNI is compared to three of the best methods in DREAM5, namely GENIE3 [START_REF] Huynh-Thu | Inferring regulatory networks from expression data using tree-based methods[END_REF], CLR [START_REF] Faith | Large-scale mapping and validation of escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF] and TIGRESS [START_REF] Haury | Tigress: Trustful inference of gene regulation using stability selection[END_REF], with their MATLAB implementation and their default parameters. Table 5.2 shows the performance of LocKNI compared to those of GENIE3, TIGRESS and CLR, and other methods presented in this thesis-Jac from Chapter 3 and Jac(LocKNI), the GRN inference through the partial derivatives of LocKNI. On network N4 from the DREAM5 challenge all the methods perform equally and just slightly better than a random guess in terms of AUROC. Such a result suggests that this network inference problem is too difficult given the avail-able data. Jac(LocKNI) has better performance on this dataset, but 55.99 AUROC is close to random guess. Overall, Jac(LocKNI) performs poorly. As LocKNI has a performant feature selection, LocKNI's model should give more weight to the partial derivative of relevant features than a Gaussian kernel-ridge regression. We see that, in practice, this is not the case. It seems that much information about partial derivatives is lost by the additive model. Results on the other networks are of greater interest. All methods perform significantly better than a random guess (e.g. an AUROC of 50%). On Network 3 and E.coli, LocKNI achieves the best performance in terms of AUROC, even outpacing GENIE3 in terms of AUPR for E.coli. On network 1, GENIE 3 leads over all other methods for the AUROC, while TIGRESS performs the best according to AUPR. Overall, LocKNI achieves state-of-the-art performance and improves it in two cases. From these comparisons follows a new question: does the prediction of LocKNI differ from the other methods? Can we benefit from that by including LocKNI in a consensus method that will gather predictions of the whole set of four methods?

Network inference by consensus of methods

As in Chapter 3, consensus of GRNI inference methods is used: one consensus (C) consists of all competing methods, a second consensus (C+L) consists of LocKNI and all competing methods, a third one (C+J) consists of Jac and competing methods, and a fourth consensus (C+L+J) consists of all competing methods, LocKNI and Jac together. On the DREAM5 datasets, the competing methods are the 30 best contestants' submissions; on the E.Coli dataset, the competing methods are TIGRESS, GENIE3 and CLR. Results are shown in Table 5.3. Interestingly, AUROC is always improved by including LocKNI in the ensemble. However this improvement is small, approximately 1%, except for Network 3 of the DREAM5 challenges. In this case LocKNI outperforms all the other methods and also provides improvements when combined with the consensus models. This improvement is obtained while being one vote out of thirty, which indicates that LocKNI grasps different information than the other methods. Note that consensus with LocKNI is usually better than consensus with Jac. This is coherent as LocKNI is a better method on these datasets. Nevertheless, on EColi, the network where LocKNI is the best in AUROC and AUPR, consensus with Jac outperforms the one with LocKNI. This shows that complementary methods are necessary for the consensus approach. Besides, using all methods (C+L+J) shows little improvement over adding LocKNI or Jac to the consensus.

Error analysis on Network N1

As LocKNI differs from the other network inference methods, the different types of errors are studied here, following the categorization introduced in [45]:

• a reverse edge means that a regulation from i to j is inferred instead of the regulation j → i,

• siblings: LocKNI infers a link from i to j because those genes are related through a third gene z. Gene i may be a "grand-father" of j (i → z → j), For all DREAM5 networks, one method's prediction is the rank it gives to each interaction, from most likely to least likely. a "grand-son" ( j → z → i), a "brother" (i ← z → j) or i and j may be a couple (i → z ← j),

• great-siblings: genes i and j are connected through two other genes,

• others: all links i → j which do not belong to any of the previous categories.

LocKNI's errors were inspected on the simulated dataset (DREAM5 Network 1), where the class of each error can be determined. The ROC curve of all types of targets is shown in Figure 5.2. For a category, such as reverse edges, I look at precision and false positive rate of LocKNI's prediction if I consider that reverse edges, without true edges, were the links to find. If LocKNI regularly made one type of mistake, then the area under the ROC curve for this category would be high. Figure 5.2 shows that LocKNI's errors are rather well balanced. It has a small bias towards selecting siblings rather than other edges. LocKNI seems as sensitive to great-siblings as it is to others. Note that LocKNI appears to identify directionality well, as reverse edges are rarely selected.

Incorporation of prior knowledge

As shown in section 5.2.2, in addition to being able to retrieve networks from the experimental data, LocKNI inherits from the regularization framework the ability to incorporate prior knowledge. A simple way to include evidence about a regulation j → i consists of setting the hyperparameter λ j in Eq. 5.8 according to the available information. In knock-out or knock-down experiments information about existence of connections between genes are close to being elucidated. One can infer that gene j regulates gene i if gene i is very perturbed in gene j knock-out experiments. This can be measured by Z-scores [START_REF] Pinna | From knockouts to networks: Establishing direct cause-effect relationships through graph analysis[END_REF]. Let µ i wt be the mean of gene i's expression level in wild type experiments, σ i wt be its standard deviation in these experiments and µ i j be its mean concentration in gene j's perturbation experiments. Pinna et al. [START_REF] Pinna | From knockouts to networks: Establishing direct cause-effect relationships through graph analysis[END_REF] define the metric:

W i j = µ i j -µ i wt σ i wt , (5.17)
which has given very good results for network inference on realistically simulated datasets [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF]. For prediction of target gene i, I normalize the sparsity constraint by dividing λ by the weight λ j = |W i j | if the experiment with gene j knocked-out is available, λ j = 1 otherwise. I examine the algorithm behavior with prior knowledge on a limited size DREAM3 network of 50 genes, where knock-out's are systematically provided. The data set consists of 21 time point measurements and of steady-state measurements. In the first experiment, I measure average AUROC and AUPR obtained Another aspect of networks is their decomposition into modules, where modules correspond to groups of genes strongly connected in the regulation graph. Such an assumption is reasonably realistic and may be obtained in some cases from gene ontology devoted to biological processes. In order to integrate this hypothesis, the sparsity constraint may be modulated by differentiating two kinds of regulation weights in the models: the weights that concern intra-module edges and those that concern inter-module edges. For each subtask P i trained on a given subset and a target gene i, the sparsity constraint is relaxed by dividing the hyperparameter λ by a parameter λ j : edges within a module are encouraged and for them, λ j is chosen above 1 (for instance, 2 in the numerical results) to reduce the effect of the sparsity constraint; edges between two modules are not encouraged and sparsity is imposed with a larger strength by setting λ j = 1. The network of DREAM3 Size 50 E.Coli1 and its decomposition in four modules is shown in Figure 5.4. To decompose this network in two modules, I consider modules 1 and 2 to be one module, modules 3 and 4 to be the other module. Table 5.4 shows the drastic improvement provided when prior knowledge about module decomposition become more precise: when using the assumption that the network can be decomposed into four modules, AUROC reaches 79.5% ad the AUPR doubles when a rough knowledge about modules (two modules) is replaced by a more accurate one (four modules).

Conclusion

I have proposed a new model-driven network inference method, LocKNI, that learns sparse nonlinear models and then extracts an estimate of the target regulation graph matrix from estimated models. Interestingly, this kernel-based method shares some features with linear approaches, such as the regularization framework, and some features with ensemble-based methods such as randomization on both variables and individuals. Compared to tree-based algorithms based on the greedy and incremental minimization of a local loss, kernel-based models are derived from the minimization of a global loss function under constraints, whose choice offers some versatility. I show that the sparsity constraint can be modulated according to the prior knowledge about existing edges. From probabilistic and Bayesian points of view, this is similar to choosing different variances in Laplace priors of the model. In practice, evidence given by knock-out data can be used to relax the sparsity constraint on potential edges. Another kind of prior knowledge, about the decomposition of the target network into modules, can be easily incorporated by imposing a lower degree of sparsity within a module and encouraging sparsity between modules. Another feature of LocKNI is the versatility it inherits from the property that kernels can be built using a convex combination of base kernels. Using 1D projection kernels, as presented in Chapter 4, provides a way to select regulators. This can be further extended to joint regulations by considering 2D projection kernels, that involve pairs of regulators.

Conclusion

In this thesis, I studied the gene regulatory network inference problem, from a bioinformatic point of view. I attempted to develop machine learning tools to extract information about the gene regulatory network (GRN), from gene expression data. This work started by noticing that many methods had been tried for this problem, but kernel methods had received little attention, although they have both relevant theoretical advantages-they are nonparametric, robustness-to-noise; they can estimate any function-and practical advantages-they have low computational costs and give good results in practice on many applications-. I have made two contributions: from a theoretical point of view, I have shown that the mean of partial derivatives is estimated consistently by some Gaussian kernel methods; from a practical point of view, I introduced and developed LocKNI, an interpretable kernel method.

On real and realistically simulated datasets, I have obtained interesting results with Jac, a new partial derivative estimation method, and state-of-the-art results with LocKNI. Currently, no method prevails in GRN inference, a notoriously difficult problem. The best way to infer a GRN is to average networks inferred by several methods. In my opinion, in this "large p, small N" framework with noisy data, combining various methods will remain the best methodology. Even if one had the perfect model, it would not be well learned on such a hard dataset. By using various models, each one makes error in its own way. Gathering base learners that make independent mistakes would give a predictor that makes fewer errors. This thesis was initiated by the idea that kernels grasped an information that could not be seen by other methods-for example, linear methods cannot understand nonlinear behavior; tree-based models provide piece-wise constant functions, and not the smooth functions that may be built using kernels. Jac and LocKNI substantially enhanced the network inferred by the consensus of other methods. Besides, even if LocKNI appears more accurate that Jac, their contribution to the consensus are very similar, hinting that the main information gain was not obtained by the way kernels were used, but by the use of kernels.

The methods introduced in this thesis may be applied to other problems. For example, these methods could be useful in biomarker discovery for breast cancer prognosis. Breast cancer prognosis is an important challenge, as better prognosis can save lives. Studies have shown that gene expression is relevant to detect subclasses of breast cancer. Scientists have searched for a signature, i.e. a list of genes that contain prognostic power for breast cancer. Using expression data from healthy and ill patients, LocKNI could identify a signature, and be added to other classification methods used on breast cancer prognosis. Another example is the differentially networking problem. One wants to identify subparts of the regulatory network that change between healthy and disease-affected tissues. One could infer two networks using LocKNI, one on healthy tissues, another on disease-affected ones. The two networks could be compared.

Several improvements may be introduced. For Jac, I used a very general regression method: kernel-ridge regression. One may use a regression better suited to gene modeling. As long as this regression asymptotically finds the true model f and its partial derivatives are bounded asymptotically, the method will consistently estimate any continuous linear form of the partial derivatives of f . For example, N.Lim [START_REF] Lim | Okvar-boost: a novel boosting algorithm to infer nonlinear dynamics and interactions in gene regulatory networks[END_REF] uses partial derivatives of operator-valued kernels. His regression model learns a structured output, better adapted to gene dynamic modeling.

LocKNI can be improved by adding prior information, as shown in Chapter 5, using Z-score to weight potential regulators of a target gene. Knowing that perturbational data are more informative than others, one may change LocKNI to use Z-score weighting if available, or to give additional weight to perturbational data. Because LocKNI has a functional cost and a simple optimization scheme (solving a quadratic problem on α and doing a gradient descent on w), it can easily be modified to take into account extra information.

Besides, this manuscript also contains several failed attempts, such as using p-values to estimate the adjacency matrix, or trying unsuccessfully many hyper-parameters selection methods. These results should not be overlooked, as they may provide information to further work on GRNI. The hyper-parameter selection problem is strategic. Many GRN inference methods choose beforehand their hyper-parameters. Performances can be increased or degraded by choosing the right or wrong hyper-parameters. Yet, the most common hyper-parameter selection methods does not seem relevant for GRN inference. Maybe they need adjustments. Maybe a new metric has to be created. Computational cost is not a limiting issue when using ensemble methods and evaluation on out-of-bag samples, the computational burden is "only" multiplied by the number of hyper-parameters to try. This is affordable with several state-of-the-art GRN inference methods, and may substantially upgrade the methods.

Finally, many mathematical models have been tried for GRN inference. Many tools are available freely online (GENIE3, TIGRESS, CLR, ANOVA, LocKNI, etc), in a ready-to-use format. An important improvement will come: more data will be available, thus increasing the size of the learning sets and methods will be able to achieve a better identification of their optimal model. Maybe we will then reach the limits of the simplified view, modeling the network only with genes, and scientists will have to incorporate all actors of gene regulation (mRNA, protein, microRNA, etc). This seems far away. For the near future, biologists have at hand powerful tools, and the significant improvements should occur by gathering data and using theses tools.

Bias-variance tradeoff

Considering that the loss function is the squared-loss one, the expected true risk of f according to the distribution of the learning set S is given by: E S (R( f )) = E S E P y -f (x) 2 [START_REF] Albert | Scale-free networks in cell biology[END_REF] = E S E P ǫ + f (x) -f (x) 2 [START_REF] Stumpf | Subnets of scale-free networks are not scale-free: sampling properties of networks[END_REF] = E S E P ǫ 2 + E S E P f (x) -f (x) 2 (20)

+2E S E P f (x) -f (x) ǫ (21)
As E P [ǫ] = 0 and ǫ is independent from x, E P f (x) -f (x) ǫ = 0. Let σ 2 be the variance of the noise σ 2 = E P [ǫ 2 ]. f * is the minimizer of the true risk, thus E S [ f (x)] = f * (x).

E S (R( f )) = σ 2 + E S E P f (x) -f * (x) + f * (x) -f (x)
f admits a representation of the form:

f (.) = N i=1 α i k(., x i ) ( 30
)
with α i ∈ R for all 1 ≤ i ≤ N Proof. Given a kernel function k, thus an RKHS H and a feature map φ (not unique). Let E be the linear span of the mappings φ(x i ) in the RKHS H, E = span ({(φ(x i )) i=1...N }) ⊂ H. Let E T be its orthogonal complement, E T E = H. Let f ∈ H:

f = v E + v E T Decomposition of f in v E ∈ E and v E T ∈ E T (31) f = N i=1 α i φ(x i ) + v Definition of E (32) f (x j ) = v E + v E T , φ(x j )
H Because H is the RKHS of k (33)

f (x j ) = v E , φ(x j )
H Because E T is the orthogonal complement [START_REF] Wille | Sparse graphical gaussian modeling of the isoprenoid gene network in arabidopsis thaliana[END_REF] of E, and φ(x j ) ∈ E (35) f (x j ) = v E (x j ) ∀ j ∈ {1, . . . , N}

Thus a function f ∈ H will have the same value on the training set than its projection in E. Then comes equality of the loss functions L((y 1 , f (x 1 )), . . . , (y N , f (x N ))) = L((y 1 , v E (x 1 )), . . . , (y N , v E (x N ))) ( 38)

f 2 H = < f, f > H (39) = < v E + v E T , v E + v E T > H (40) = v E 2 H + v E T 2 H + 2 < v E T , v E > H =0 because E and E T orthogonal (41) = v E 2 H + v E T 2 H (42) ≥ v E 2 H (43)
So function f has same risk function as v E . f ≥ v E , and Ω is a strictly increasing function, so the minimum of the risk function plus regularization term is reached in E, so the function f * has the form:

f * = N i=1 α i φ(x i ) (44) f * (x) = N i=1 α i k(x i , x) (45) (46)
Theorem 2 Let f be a function in space C 1 (X c , Y). Assume f has gradient bounded by M, ∇ f 2 < M. Then, for any ǫ > 0, there exists C 1,ǫ such as: . Clearly this set covers X c . By the Borel-Lebesgue theorem, we can extract a finite set of those balls (B i ) i=1...N that will cover X c . As X c is compact and f continuous on X c , f ∞ is reached in a point x * ∈ X c . Let B i * be a ball of the finite set that contains x * . For all x in B i * , we have xx * 2 ≤ ǫ/M because they both belong to a ball of radius ǫ/(2M).

f ∞ ≤ C 1,ǫ f L 2 (Ω,µ) + ǫ (47
Thus we have:

B i * ∩X c | f (x)|µ(x)dx ≥ B i * ∩X c f ∞ -ǫ µ(x)dx (51)
≥ f ∞ -ǫ µ(B i * ∩ X c ) (52)
Noting µ(B i * ∩ X c ) = B i * ∩X c µ(x)dx. Using a Cauchy-Schwarz theorem, we have:

B i * ∩X c | f (x)|µ(x)dx ≤ B i * ∩X c f (x) 2 µ(x)dx B i * ∩X c 1µ(x)dx (53)
≤ f L 2 (Ω,µ) µ(B i * ∩ X c ) (54)
Using equation ( 52) and ( 54), we have:

f ∞ -ǫ µ(B i * ∩ X c ) ≤ f L 2 (Ω,µ) µ(B i * ∩ X c ) ( 55
)
f ∞ ≤ f L 2 (Ω,µ) µ(B i * ∩ X c ) + ǫ (56)
The serie (µ(B i ∩ X c )) i is finite, so there is a minimum µ(B * ). As X c is the closure of an open set, µ(B i ∩ X c ) is stricly greater than 0 for all i, so µ(B * ) is strictly greater than 0. So

f ∞ ≤ f L 2 (Ω,µ) µ(B * ) + ǫ (57)
Theorem 3 Let f be a C 1 (X c , Y) function, with gradient bounded by M. Let g be a linear continuous form of C 0 (X c , Y), g : C 0 (X c , Y) → R. Then, for any ǫ > 0, there exists C 2,g,ǫ such that:

|g(∇ f )| ≤ C 2,g,ǫ f ∞ + ǫ ( 58
)
Proof. We are going to prove this in one dimension, without loss of generality:

g ∂ f ∂x j ≤ C 2,g,ǫ f ∞ + ǫ (59)
As we are in a finite-dimensional space, this suffices to prove for g(∇ f ). By Fréchet-Riezs'theorem, there exists h ∈ C 0 (X c , Y) such that:

g ∂ f ∂x j = X c h(x) ∂ f ∂x j (x)dx (60)
By Heine's theorem, as X c is compact, h is uniformly continuous on X c . Thus, there exists δ such that, for all (x, x ′ ) with xx ′ 2 < δ, we have

|h(x) -h(x ′ )| ≤ ǫ MV(X c ) (61)
with V(X c ) = X c 1dx. As we did in the proof of theorem (2), we can cover X c with a set of balls (B i ) i=1...N of radius δ. We note

P i = B i ∩ X c . g ∂ f ∂x j = i P i ∂ f ∂x j (x)( h(x) -h i ≤ǫ/(MV(X c ))
+h i )dx [START_REF] Zou | A new dynamic bayesian network (dbn) approach for identifying gene regulatory networks from time course microarray data[END_REF] with f (x max j )f (x min j )dx -j [START_REF] Kim | Dynamic bayesian network and nonparametric regression for nonlinear modeling of gene networks from time series gene expression data[END_REF] With x max j = arg max u∈P i {u j , u -j = x -j } and x min j = arg min u∈P i {u j , u -j = x -j }. Then, for all f ∈ H, for all x ∈ X c :

h i = h(b i ),
g ∂ f ∂x j ≤ ǫ + i 2 f ∞ h i V(P -j i ) (65)
≤ ǫ + 2V(X -j c ) h ∞ f ∞ (66
∇ f (x) 2 ≤ √ M f H (67)
Proof. Let f ∈ H. We have:

| f (x) -f (x ′ )| = | < f, K x -K x ′ > H | (68) ≤ f H K x -K x ′ H (69)
Using the property of reproducing kernel and Cauchy-Schwarz's theorem. Using polarization identity, the kernel trick, Taylor expansion and lemma (4)'s assumptions, we have:

K x -K x+h 2 H = K x 2 H + K x+h 2 H -2 < K x , K x+h > H (70
) = k(x, x) + k(x + h, x + h) -2k(x, x + h) (71) = 2c -2c -2 < ∇ x k(x, x), h > (72) -< h, H(x, x)h > +o( h 2 2 ) (73) = -< h, H(x, x)h > +o( h 2 2 ) (

From equations [START_REF] Qian | Inference of noisy nonlinear differential equation models for gene regulatory networks using genetic programming and kalman filtering[END_REF][START_REF] Grimaldi | Regnann: Reverse engineering gene networks using artificial neural networks[END_REF], we have:

| f (x) -f (x + h)| = | < ∇ x f, h > Ω +o( h 2 )| (75) ≤ f H h √ M + o( h 2 ) (76)
Thus

∇ f 2 ≤ √ M f H (77)
Lemma 5 The Gaussian kernel of bandwidth σ satisfies the hypothese of lemma [START_REF] Huynh-Thu | Machine learning-based feature ranking: Statistical interpretation and gene network inference[END_REF].

Proof. This result is obtained by a few calculations with the Gaussian kernel:

k(x, z) = exp - x -z 2 2 2σ 2 (78) k(x, x) = 1 (79) ∂k ∂z i (x, z) = -(z i -x i ) σ 2 exp - x -z 2 2 2σ 2 (80) = 0 if x = z (81) ∂ 2 k ∂z i ∂z j (x, z) = (z i -x i )(z j -x j ) σ 4 exp - x -z 2 2 2σ 2 (82) = 0 if x = z (83) ∂ 2 k ∂(z i ) 2 (x, z) = (z i -x i ) 2 σ 4 - 1 σ 2 exp - x -z 2 2 2σ 2 (84) = -1 σ 2 if z = x (85)
We have the Hessian matrix H(x, x) = -Id/σ 2 , which has eigenvalues bounded by 1/σ 2 .

Theorem 6 With samples (x i , y i ) i=1...ℓ i.i.d., for any ǫ, η > 0, any continuous linear form g of C 0 (X c , Y), g : C 0 (X c , Y) → R, there exists ℓ 0 such as if ℓ ≥ ℓ 0 , then, with probability greater than 1η:

X c g(∇ fℓ (x)) -g(∇ f (x))dx ≤ ǫ ( 86
)
with fℓ the estimator of f based on a Gaussian kernel ridge regression or Gaussian partial least-square regression.

Proof. From [START_REF] Caponnetto | Optimal rates for the regularized leastsquares algorithm[END_REF], we have consistency of kernel-ridge regression. From [START_REF] Blanchard | Kernel partial least squares is universally consistent[END_REF], we have consistency of partial least-square regression. In both articles, authors prove that the difference between the estimator and the true function is bounded in both the norms . H and . L 2 (Ω,µ) . Using lemma (4) and bound on flf H , we obtain that, for ℓ sufficiently large, we have, with probability greater than 1η:

∇ fℓ (x) -∇ f (x) 2 ≤ 1 for all x ∈ X c [START_REF] Steppe | Feature saliency measures[END_REF] Let g be the linear form. By Fréchet-Riezs'theorem, there exists ν ∈ C 0 (X c , Y) such that::

g : f → A f (x)ν(x)dx (88)
Using theorem (3):

g ∂ fℓ -f ∂x j ≤ C 2,g,ǫ fℓ -f ∞ + ǫ 3 ( 89
)
Using theorem (2):

fℓ -f ∞ ≤ C 1,ǫ fℓ -f L 2 (Ω,µ) + ǫ 3C 2,g,ǫ (90)
Using consistency result, for ℓ sufficiently large, we have, with probability greater than 1η:

fℓ -f L 2 (Ω,µ) ≤ ǫ C 1,ǫ C 2,g,ǫ 3 (91)
Combining the 3 equations [START_REF] Brabanter | Nonparametric derivative estimation[END_REF][START_REF] Charnigo | A generalized c p criterion for derivative estimation[END_REF][START_REF] Khemchandani | Regularized least squares support vector regression for the simultaneous learning of a function and its derivatives[END_REF], we obtain the result:

A ∇ fℓ (x) -∇ f (x)ν(x)dx ≤ ǫ (92)
Bounds for integral estimation Let f be a C 1 (X c , Y) function whose partial derivatives are bounded by M. Let (x i ) i=1...ℓ be independently and identically distributed random variables, drawn from a distribution µ. The difference between the empirical mean of a partial derivative

S l = 1 ℓ ℓ i=1 ∂ f ∂x j (x i )
and its true mean

e = E ∂ f ∂x j = X c
∂ f ∂x j (x)µ(x)dx follows the normal distribution N(0, ℓ -1 σ 2 ), with σ ≤ M. Thus we have

P        1 ℓ ℓ i=1 ∂ f ∂x j (x i ) - X c ∂ f ∂x j µ(x)dx ≤ m        ≤ Φ mM √ ℓ ( 93
)
with Φ(m) the cumulative distribution of N(0, 1)

Proof. Firstly, I show that the variance of ∂ f ∂x j is lower or equal to M 2 . Then I verify that I meet the central limit theorem's assumption, and apply it to find the previously stated bounds.

Var

∂ f

∂x j = E       ∂ f ∂x j -e 2       (94)
≤ M 2 p(2 -2p) 2 + (1 -p)(2p) 2 (98)
≤ 4M 2 (p(1 -p)) (99)
And p(1p) ≤ 1 / 4 , minimum reached in p = 1 / 2 . Thus:

Var ∂ f ∂x j ≤ M 2 (100)
As x i are i.i.d., ∂ f (x i ) ∂x j are i.i.d.. Mean and variance of ∂ f (x i ) ∂x j are defined and finite, thus, I can apply the central limit theorem. Noting S ℓ = ℓ i=1 ∂ f (x i ) ∂x j , and σ the standard deviation of ∂ f (x i ) ∂x j , the theorem implies that:

S ℓ -ℓe σ √ ℓ -→ ℓ→+∞ N(0, 1) (101)
or

P S ℓ ℓ -e ≤ m = Φ mσ √ ℓ (102)

Figure 1 . 3 :

 13 Figure 1.3: Examples of loss functions

3 Figure 1 . 4 : 2 Figure 1 . 5 :

 314215 Figure 1.4: Illustration of the curse of dimensionality. As dimension increases, data become more sparse.

Figure 1 . 2 Figure 1 . 6 :

 1216 Figure1.6: Illustration that ℓ 1 -norm induces sparsity. a * is the minimizer of the empirical risk. a ℓ • is the minimizer of ℓ • -regularized empirical risk. a ℓ • will be the first intersection between a level set of the empirical risk and the level set of the regularization norm. Due to its form, the ℓ 1 -norm level set intersects the empirical risk level set at a vertex, where many features are null.

Figure 1 . 7 :

 17 Figure 1.7: Illustration of a nonlinear problem transformed in a linear problem by the appropriate map φ

Figure 1 . 8 :

 18 Figure 1.8: Examples of data (x, y) and their Pearson correlation. Each observed couple (x, y) is represented by a blue dot. The regression y = ax is drawn in red. Pearson's correlation can miss nonlinear dependencies

Figure 1 . 9 :

 19 Figure 1.9: Different representations of a Boolean function. Left: wiring diagram. →: positive regulation (excitate); ⊣: negative regulation (inhibit). Right: truth table

Figure 1 . 10 :

 110 Figure 1.10: Example of a Bayesian network, represented as a graph or an adjacency matrix

Figure 1 . 11 :

 111 Figure 1.11: graph of a dynamical Bayesian network

Figure 2 . 1 :

 21 Figure 2.1: Examples of functions of the form given by equation (2.19) with dimension p = 2

  of x error on partial derivative estimation

Figure 2 . 2 :

 22 Figure 2.2: Mean and variance of p m=1 (g(∂ f * /∂x m )g(∂ f /∂x m )) 2 . The top plots show the error in function of the number of training samples N for noise to signal ratio 0% (top left) and 100% (top right), dimension p = 2. The bottom plots show the error in function of the signal-to-noise ratio σ (bottom left -with p = 3, N = 250) and the dimension of x p (bottom right σ = 25%, N = 250)

Figure 2 . 3 :

 23 Figure 2.3: Error in the estimation of partial derivatives when input variables x are drawn from a Gaussian distribution, and the integral is estimated on the cube [-1, 1] p

  Figure 3.5. I chose n var = 5 and p data = 80%. • The regularization-data fitting tradeoff λ and kernel bandwidth σ are selected through cross-validation. With subsampling, test error can be calculated by out-of-bag data, i.e. data not selected in the subsample. Performing cross-validation for M values of λ only costs M times the cost of performing the algorithm one time. σ can take the following values {0.1, 0.5, 1, 2, 5}. With E 0 the largest dataset, p data |E 0 | is the average size of a subsample. Let s = 1/ p data |E 0 |. λ can take the following values {0.01s, 0.1s, s, 10s, 100s}.

Figure 3 . 3 :Figure 3 . 4 :

 3334 Figure 3.3: Comparison of importance measures on DREAM3 size 50 and size 100 networks

Figure 3 . 5 :

 35 Figure 3.5: Sensivity of Jac to parameters p data and n var on DREAM3 size 50 networks

  embedded feature selection methods have been successfully used for GRN inference. In this context, I propose a new feature weighting methods. Using the versatility of kernels, I define local kernels, that rely on only one feature. I propose several methods to optimize local kernel-based models. Moreover, I test various popular kernel feature selection or feature weigthing methods along with the local kernel approach on realistically simulated GRN. Firstly, I describe kernel feature selection methods, and how to use them for GRN inference; Secondly, I describe hyper-parameter selection methods. Finally, I show the experiments and experimental results.

1

 1 ,..., f M L( f 1 , . . . , f M , w)

ℓ 2 -

 2 MKL I use the MKL framework of SimpleMKL, but I change the loss function for the square error loss function l(y i , m f

  (a) over α with w fixed and (b) over w with α fixed. For the (a) step, I have a closed-form solution:

Algorithm 4

 4 SVM-RFE for linear SVM Input: training examples X 0 = [x 1 , . . . , x n ] with labels y = [y 1 , . . . , y n ] Initialize remaining features F = [1, . . . , p], feature ranked list R = [] while F [] do Remaining training features X = X 0 (F, :) Train the classifier w = S V M train(X, y) Find the least informative feature f = arg min i w 2

  end while Output: feature ranked list R Algorithm 5 SVM-RFE for non-linear SVM Input: training examples X 0 = [x 1 , . . . , x n ] with labels y = [y 1 , . . . , y n ] Initialize remaining features F = [1, . . . , p], feature ranked list R = [] while F [] do Compute the Gram matrix K F Train the classifier w = S V M nonlinear train(K F , y)

Figure 4 . 1 .

 41 For each feature selection method and for each size 10 network, I compute the correlation

Figure 4 . 1 :

 41 Figure 4.1:For SimpleMKL, values of the points (MS E evaluated by crossvalidation, AUROC) for the different values of the hyper-parameter µ in the predefined set. The AUROC appears to depend more on the dataset than the hyperparameters. Second, the shape of the distribution of the points varies greatly from one network to another. Cross-validation does not seem to be a good hyperparameter selection method when SimpleMKL is used for the task of network inference.

Figure 4 . 2 :

 42 Figure 4.2: Top: Comparison of the AUPR of described methods on the five size 10 networks. Bottom: Comparison of the AUPR of RFE Gaussian and SU-PANOVA to best competing teams in DREAM3 Challenge

Figure 4 . 3 :

 43 Figure 4.3: Comparison of kernel feature selection methods on DREAM 3 Size 50 Networks

= 1 ,

 1 w m ≥ 0, ∀m ∈ {1, . . . , M}(5.9)

Figure 5 . 1 :

 51 Figure 5.1: PCA on the prediction vector of all methods. For all DREAM5 networks, one method's prediction is the rank it gives to each interaction, from most likely to least likely.

Figure 5 . 2 :

 52 Figure 5.2: Types of links inferred by LocKNI for DREAM5 network 1.

Figure 5 . 3 :

 53 Figure 5.3: (a) Average AUROC versus the number of KO experiments, using E. coli 1, in DREAM3 challenge. (b) AUROC versus the number of chosen KO experiments, using E. coli 1, in DREAM3 challenge.

Figure 5 . 4 :

 54 Figure 5.4: Network E.Coli1 size 50 from DREAM3

2 ( 22 )= σ 2 + 2 ( 23 )E 2 ( 27 )= σ 2 + 28 ) 1

 22222232272281 E S,P ( f (x)f * (x)) 2 + E S,P f * (x) -f (x) +2E S,P ( f (x)f * (x))( f * (x) -f (x)) E S [E P [.]] = E P [E S [.]] = E S,P[.]. f and f * do not depend on the learning set S, and E S f * (x) -f (x) = 0, soE S,P ( f (x)f * (x))( f * (x) -f (x)) = E P             ( f (x)f * (x)) E S f * (x) -f (x) S (R( f )) = E P,S [ǫ 2 ] + E P ( f (x)f * (x)) 2 + E S,P f * (x) -f (x) (biais) 2 + (variance) (Theorem Representer theorem.Let X be a nonempty set and k a positive-definite real-valued kernel on X × X with corresponding reproducing kernel Hilbert space H. Given a training sample (x 1 , y 1 ), . . . ,(x N , y N ) ∈ (X × R) N ,a strictly monotonically increasing real-valued function Ω : [0, ∞) → R, and an arbitrary empirical risk function L : (R × R) N → R ∪ {∞}, then for any f * ∈ H satisfying f = arg min f ∈H L ((y 1 , f (x 1 )), . . . , (y N , f (x N ))) + Ω( f )

) and C 1 ,

 1 ǫ is independant of f . Proof. Consider the infinite set of open balls B(x, ǫ 2M ) x∈X c

|

  f (x)f (x * )| = < ∇ f (c), xx * >for some c, by the mean value theorem (48) ≤ M xx * 2 Cauchy-Schwarz and bounded gradient (49) ≤ ǫ

) Lemma 4 •

 4 Let H be the RKHS of universal kernel k. If, for all x ∈ Ω • The kernel is constant k(x, x) = c• On point (x, x), the gradient of the kernel is null ,∇ z k(x, z)| z=x = 0 p On point (x, x), the Hessian matrix H(x, z) i j = ∂ 2 k(x,z) ∂z i ∂z j has eigenvalues bounded by a constant M, | < u, H(x, x)u > | ≤ M u 2 for all x ∈ X and all u ∈ R p

  Let p = P ∂ f ∂x j (x) ≥ e . In the worst-case scenario where the partial derivative only takes the values {-M, M}, e = M(2p -1).Var ∂ f ∂x j ≤ p(M -M(2p -1)) 2 + (1p)(-M -M(2p -1)) 2(97)

  1. Randomly choose from which dataset E u to extract data. The probability of choosing dataset E u is proportional to its size, p(E u selected) = |E u |/( N e ℓ=1 |E ℓ |). Let E u be the selected dataset. 2. Each data vector x i of E u has probability p data,u to be in subsample E ℓ . To obtain subsamples of similar sizes, p data,u is taken inversely proportional to E u 's size. Assuming, without loss of generality, that E 0 is the largest dataset, I have fixed p data,u = p data E 0 E u , with p data a fixed hyper-parameter. To obtain test data,this probability is capped at 95%. 3. Randomly select n var variables to be the potential regulators of our system.

	Available time-series
	[gene 1] after perturbation 1			[gene 1] after perturbation 2
	[gene 1]						[gene 1]
	1						1
	0.75						0.75
	0.50						0.50
	0.25						0.25
	0 15 30 45 60		120	time (s)			0	30	60	90	120	time (s)
	Dataset data collected
	[gene 1]
	1
	0.75
	0.50
	0.25
		0 15 15 30 30 45 45 60
	E dyna,15s
	[gene 1]
	1
	0.75
	0.50
	0.25
		0	30 15	45 30	60	0	30 30	60 60	90 90	120
	E dyna,30s
	. . .	. . .
	Figure 3.1: Example of irregularly spaced time-series, and the building of datasets
	E dyna,.

Table 3 .

 3 1: Using a threshold p-value of Impo < θ%, one expects to obtain false positives, precisely θ% times number of Impo calculated. Noting EFP the number of Expected False Positives, the table shows the mean ratio of number of p-values superior to threshold divided by the EFP and gives the number of networks (out of 5) for which this ratio exceeded 1.

	Edge threshold Method	pVal ≤ 5% P/EFP #{(P/EFP) > 1} P/EFP #{(P/EFP) > 1} pVal ≤ 1%
	Jac	0.53	1	2.27	3
	rel test	1.41	5	5.71	5
	rel train	3.86	5	16.34	5
	sens	0.24	0	1.02	2

Table 3 .

 3 3: Performance measured as AUROC and AUPR in % on DREAM5 Networks 1, 3 and 4. The results in bold font are the best ones. The first row is the consensus of other methods which are all the participating teams on DREAM5 and GENIE3, TIGRESS and CLR onEColi .

				Static Modeling			Dynamic modeling
			AUROC with static modeling without subsampling		AUROC with dynamic modeling without subsampling
								0.65					jac
		0.65 0.6						0.6					sens rel train rel test
	AUROC	0.45 0.5 0.55					AUROC	0.5 0.55
	AUROC	0.35 0.4	ECol1 rel test jac sens rel train	EColi2	Yeast1	Yeast2	Yeast3	0.4 0.45	ECol1	EColi2	Yeast1	Yeast2	Yeast3
			AUPR with static modeling without subsampling		AUPR with dynamic modeling without subsampling
		0.12	jac					0.1	jac
		0.1	sens rel train rel test					0.08	sens rel train rel test
	AUPR	0.06 0.08					AUPR	0.06
	AUPR	0.02 0.04	ECol1	EColi2	Yeast1	Yeast2	Yeast3	0.02 0.04	ECol1	EColi2	Yeast1	Yeast2	Yeast3
	Figure 3.6: Results of importance measures without subsampling, using either
	static or dynamic modeling

Table 4 .

 4 .1 and 4.2. MS E, BIC and AIC should be negatively correlated with AUROC or AUPR; stability and block-stability should be positively correlated. First, we note that too few examples are available to be able to conclude that any hyper-parameter selection criterion is good. We can conclude that most criteria give poor results, no better than what would be obtained with an educated guess. Second, the correlation is highly unstable. For ℓ 2 -MKL and for AUROC, one would rather use stability than cross-validation: stability has lower correlation (0.44 against -0.54) 2 , but the standard deviation of this correlation is much lower (0.16 against 0.34). Third, for RFE, the BIC criterion seems more relevant than other methods, even though the results have a high variance. 1: Correlation between various hyper-parameter selection criteria and the AUROC for six feature selection methods for network inference. The values displayed are the mean ± standard deviation over size 10 networks. The best hyperparameter selection criterion for a feature selection method is shown in bold. For SimpleMKL and KNIFE, the criterion is either wrongly correlated with the feature selection's performance or the correlation is close to 0. We conclude that none of the considered hyper-parameter selection criteria works well in this context.

	FS Method	CV MSE	Stability	Block Stability	BIC	AIC
	SimpleMKL	0.1±0.48 -0.35±0.15	-0.49±0.08	-0.06±0.14 -0.02±0.14
	ℓ 2 -MKL SUPANOVA	-0.54±0.34 0.44±0.16 -0.23±0 -0.06±0.1	0.43±0.1 -0.02±0.07	0.47±0.36 0.45±0.37 0.22±0 -0.22±0
	KNIFE	0.24±0.18 -0.04±0.5	-0.21±0.5	-0.01±0.42 -0.01±0.42
	RFE linear	0.14±0.93	-	-	-0.4±0.35 -0.09±0.81
	RFE Gaus	0.01±0.7	-	-	-0.27±0.51 0.25±0.53

Table 4 .

 4 2: Correlation between various hyper-parameter selection criteria and the AUPR for six feature selection methods. For every method but RFE, the criterion is either wrongly correlated with the feature selection's performance or the correlation is close to 0. For RFE, the BIC criterion seems the most relevant, but is also highly unstable.

	AUROC	EColi1 EColi2 Yeast1 Yeast2 Yeast3 Avg time (sec)
	SimpleMKL	0.66	0.36	0.56	0.54	0.51	2726.17
	ℓ 2 -MKL SUPANOVA	0.65 0.5	0.41 0.5	0.49 0.5	0.57 0.5	0.63 0.65	200.78 104.83
	KNIFE	0.61	0.57	0.65	0.45	0.33	3842.98
	KA Gaus	0.58	0.54	0.43	0.62	0.53	0.45
	KA Linear	0.36	0.44	0.47	0.48	0.34	0.54
	RReliefF	0.5	0.53	0.54	0.49	0.54	0.86
	RFE Linear	0.36	0.58	0.48	0.61	0.57	0.39
	RFE Gaus	0.5	0.5	0.5	0.5	0.5	9.02
	Team 236	0.62	0.65	0.65	0.44	0.49	-
	Team 190	0.57	0.52	0.63	0.58	0.60	-

Table 4 .

 4 3: Results in AUROC of nine eleven feature selection methods and two DREAM competing methods on Size 10 networks. The best performing method on a given network is highlighted in bold. The last column gives the average computational time on one network for the whole testing process (i.e., computing each hyper-parameter selection criterion and calculating model for each selected hyper-parameter)

4.4.2 Performance on size 10 networks

  The performance is assessed by AUROC, AUPR and the computational time. For SimpleMKL, SUPANOVA, KNIFE, hyper-parameters have been selected with cross-validation, as no method stands out. Stability has been used for ℓ 2 -MKL; BIC for RFE. For AUROC, results are shown in Table4.3. No method is competitive on all networks. Many methods underperform random guessing (AUROC= 0.50) on some dataset. For AUPR, results are shown in Figure4.2. The SUPANOVA and RFE Gaussian methods stand out. We also see that these methods are competitive on these few experiments. Computational time is given in Table4.3. Given a network, the time given is the sum of the time to calculate every hyper-parameter selection criterion, and compute the prediction function f for the selected hyper-parameters. Filter methods (KA, RReliefF) are considerably faster, as expected, but the RFE methods are competitive. ℓ 2 -MKL and SUPANOVA are executed rather quickly, and can be applied to size 50 networks. The computational time of SimpleMKL and KNIFE is too large to be able to apply them to size 50 networks in this study.

  1. Randomly choose from which dataset E • to extract data. The probability p(E u ) of choosing dataset E u is proportional to its size, p(E u ) = |E u |/( N e ℓ=1 |E ℓ |). Let E u be the selected dataset. 2. Each data vector x i of E u has probability p data,u to be in subsample E ℓ . To obtain subsamples of similar sizes, p data,u is taken inversely proportional to E u 's size. Assuming, without loss of generality, that E 0 is the largest dataset, I have fixed p data,u = p data E 0

Table 5 .

 5 1: Characteristics of the datasets. # TF is the number of potential regulators, p is the number of target genes, N is the number of data in the training set, N dyn is the number, among these N data, that are part of a time-series.

			Network			# TF	p	N N dyn
	DREAM5 Network 1 (in-silico)		195 1643 805 463
	DREAM5 Network 3 (E.Coli)		334 4511 805 463
	DREAM5 Network 4 (S.cerevisiae) 333 5950 536 298
			E.Coli			169 4297 466 186
	networks	DREAM 5 N1	DREAM 5 N3	DREAM 5 N4	EColi
		AUROC AUPR AUROC AUPR AUROC AUPR AUROC AUPR
	LocKNI	78.7	25.54	65.6	7.68	51.43		1.98	65.11	8.56
	Jac	71.89	8.64	63.19	3.47	52.57		2.21	63.5	5.75
	Jac(LocKNI)	52.22	9.33	60.97	4.16	55.99		2.32	62.77	5.70
	GENIE3	81.5	29.1	61.7	9.3	51.8		2.1	64.04	6.04
	TIGRESS	78.2	30.1	59.5	6.9	51.7		2.0	64.54	3.98
	CLR	77.3	25.5	59.0	7.5	51.6		2.1	63.26	6.59

Table 5 .

 5 2: Performance of GRN inference methods measured by AUROC and AUPR in percentages on DREAM5 Network 1, 3 and 4, and on EColi, for which the data is from M3D and the gold standard is from RegulonDB. The best results are shown in bold font.

Table 5 .

 5 3: AUROC and AUPR in %. Performance of Consensus Method first with competing methods (C)-competing methods are GENIE3, TIGRESS and CLR on EColi, all DREAM5 contestants on the DREAM5 dataset-, second, also including LocKNI (C+L), third with competing methods and Jac on Gaussian kernel-ridge regression (C+J), and, finally, consensus with competing methods, LocKNI and Jac (C+L+J). The best results are shown in bold font. In the "best algorithm" row, "G3" stands for GENIE3, "TIG" for TIGRESS, "L" for LocKNI and Jac(L) for the network inferred by the partial derivatives of the LocKNI model.

	networks	DREAM 5 N1	DREAM 5 N3	DREAM 5 N4	EColi
	Method	AUROC AUPR AUROC AUPR AUROC AUPR AUROC AUPR
	Consensus
	C	80.89	32.65	64.94	8.99	52.02	2.24	66.06	5.94
	C+L	81.31	31.82	73.11	8.68	53.64	2.29	66.47	5.76
	C+J	81.23	31.81	72.78	8.75	53.58	2.30	67.93	5.99
	C+L+J	81.30	31.83	73.03	8.72	53.44	2.25	67.94	5.96
	Best-of-Single	81.5	31.5	65.6	9.3	55.99	2.32	65.11	8.56
	algorithms
	Best algorithm	G3	TIG	L	G3	Jac(L) Jac(L)	L	L

It seems that, individually, LocKNI and Jac bring complementary results to other methods, but little seems to be learned by using both LocKNI and Jac. The contribution of LocKNI and Jac is similar, probably because kernels would capture information that wasn't identified before. Following Marbach et al.

[START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF]

, I perform PCA on the contestants' predictions to see which methods were close to each other (see Figure

5

.1). LocKNI stands out as yielding different results from other approaches. The improvements in AUROC are obtained nearly without degrading AUPR.

Table 5 .

 5 4: AUROC and AUPR when prior knowledge about module existence is given (DREAM3 network E. Coli size 50). by LocKNI trained on available time courses and enhanced by 12, 25, 37 and 50 uniformly drawn single KOs out of the 50 total knock-outs. Incorporation of prior knowledge significantly improves the performance. However this experiment is very different from what a biologist would do. A biologist, expert in the studied biological system, would choose carefully the KO experiments to run by giving preference to hubs. Figure5.3 shows the results for an increasing number of KO experiments chosen according the number of known targets of transcription factors of the system. Contrary to the previous experiment, which provides average results, this study shows how a biologist can improve drastically performance of a network inference algorithm by producing a very few well chosen KOs. Only 5 well chosen single KO experiments are needed to make AUROC reach 78%.

  with b i the center of ball B i

	g	∂ f ∂x j	≤	i	P i	∂ f ∂x j (x)	ǫ MV(X c )	dx +	P -j i	P i j	∂ f ∂x j (x)h i dx j dx -j (63)
			≤ ǫ +	i	h i	p j i P

"Nature's In-depth study is the most prolific source of mathematic discoveries"

Adenine, Thymine, Guanine and Cytosine

A scale-free network is a network whose degree distribution follows a power law, at least asymptotically. This property is observed on many large-scale networks, such as the scientific collaboration network[START_REF] Albert | Statistical mechanics of complex networks[END_REF] or biological networks[START_REF] Farkas | The topology of the transcription regulatory network in the yeast, saccharomyces cerevisiae[END_REF][START_REF] Barabási | Network biology: understanding the cell's functional organization[END_REF][START_REF] Albert | Scale-free networks in cell biology[END_REF]. Note that the scale-free assumption for biological networks is contested in[START_REF] Stumpf | Subnets of scale-free networks are not scale-free: sampling properties of networks[END_REF][START_REF] Khanin | How scale-free are biological networks[END_REF]

this assumption can be made without loss of generality

the probability of the empirical estimate deviating from its mean can be bounded as an exponentially decaying function of that deviation

Note that stability should be positively correlated to performance, and MSE negatively correlated. So a correlation of -0.54 for MSE is a better result than a correlation of 0.44 for stability.