Camilo Broc
email: camilo.broc@univ-pau.fr

Marina Evangelou
email: m.evangelou@imperial.ac.uk

Pascal Guenel
email: pascal.guenel@inserm.fr

Therese Truong
email: therese.truong@inserm.fr.

Benoit Liquet
email: benoit.liquet@univ-pau.fr

Borja Calvo

Investigating Gene-and Pathway-environment Interaction analysis approaches

Keywords: Gene-environment interactions, Generalized Linear Models, Pathway analysis, Resampling methods Interaction gène-environnement, Modèles linéaire généralisés, Analyse par pathway, Méthodes de rééchantillonage AMS 2000 subject classifications: 62F03, 62F40, 62P10 Batch Effect, High dimensional data, Partial Least Square, Sparse methods

Durant les dernières décennies, la quantité de données disponibles en génétique a considérablement augmenté. D'une part, une amélioration des technologies de séquençage de molécules a permis de réduire fortement le coût d'extraction du génome humain. D'autre part, des consortiums internationaux d'institutions ont permis la mise en commun de la collecte de données sur de larges populations. Cette quantité de données nous permet d'espérer mieux comprendre les mécanismes régissant le fonctionnement de nos cellules. Dans ce contexte, l'épidémiologie génétique est un domaine cherchant à déterminer la relation entre des caractéristiques génétiques et l'apparition d'une maladie. Des méthodes statistiques spéci ques à ce domaine ont dû être développées, en particulier à cause des dimensions que les données présentent : en génétique, l'information est contenue dans un nombre de variables grand par rapport au nombre d'observations. Dans cette dissertation, deux contributions sont présentées. Le premier projet appelé PIGE (Pathway-Interaction Gene Environment) développe une méthode pour déterminer des interactions gène-environnement. Le second projet vise à développer une méthode de sélection de variables adaptée à l'analyse de données provenant de di érentes études et présentant une structure de groupe de variables.

Le document est divisé en six parties. Le premier chapitre met en relief le contexte, d'un point de vue à la fois biologique et mathématique. Le deuxième chapitre présente les motivations de ce travail et la mise en oeuvre d'études en épidémiologie génétique. Le troisième chapitre aborde les questions relatives à l'analyse d'interactions gène-environnement et la première contribution de la thèse y est présentée. Le quatrième chapitre traite des problématiques de méta-analyses. Le développement d'une nouvelle méthode de réduction de dimension répondant à ces questions y est présenté. Le cinquième chapitre met en avant la pertinence de la méthode dans des cas de pleiotropie. En n, le sixième et dernier chapitre dresse un bilan du travail présenté et dresse des perspectives pour le futur.

Remerciements/Acknowledgements

It is with great pleasure that I can dedicate the following lines to acknowledge all the people that have supported and helped me during my work during those three years.

First, I acknowledge my supervisors Benoit Liquet and Borja Calvo for directing my thesis. I felt grateful and fortunate to have you directing my work. Thank you Benoit for the constant interest you put in me. Thank you Borja for your insight and your advices that always came at the right moment.

I want to thank a lot Stéphane Robin and Christophe Ambroise for reviewing my manuscript and all the interest you showed for my work. I would like to thank also Hélène Jacqmin-Gadda for supervising the defense of this thesis and Astrid Jourdan for participating to it.

I thank the LMAP (Laboratoire de Mathématiques et leurs Applications de Pau) for hosting me. I would like to thank all the people from the center of Anglet for all the daily moments at work. I would like to thank also the Intelligent System Group in Donostia for their warm welcome. I acknowledge the UPPA and the UPV/EHU for funding my thesis and hosting me. I would like also to thank the french association Ligue contre le Cancer for its nancial support.

Molecular biology

The material in this dissertation deals with the statistical analysis of biological data. Before going into the details of the contributions of this thesis we need to review some fundamental concepts. In the rst part of this chapter basic molecular biology concepts are presented. Then, computation biology and its vocabulary is introduced. Finally data used for analyses are presented.

Biological background 1.Behaviour of a cell

The cell is the structural and functional unit of all known living organisms (cell comes from the Latin cellula, meaning, a small room). Some organisms have only one cell (they are called unicellular), whereas, other organisms are composed of several ones (they are called multi-cellular). For a human, the number of cells is estimated to be 10 14 and a typical cell size is 10 µm (10 -5 m) with an average mass of 1 ng (10 -9 g). Each cell is a small zone of water and chemicals wrapped by a membrane, which can interact with its close environment, and perform specialized functions. Two types of cell exist: eukaryotic and prokaryotic. Eukaryotic cells are often found in multi-cellular organisms. Prokaryotic cells are usually independent and they show the simplest structure. Eukaryotic cells present in general additional internal structures compared to prokaryotic cells (Figure 1.1). Notably, the major di erence between those is that eukaryotic cells contains a nucleus in opposition to prokariotic ones: a membrane-delineated compartment that houses the eukaryotic cell's DNA (Deoxyribonucleic Acid) as presented in Figure 1.2) . Eukariotic cells have the following notable internal entities that ensure its behaviour:

• The mitochondrion provides the energy needed by the cell through enzymes exchanging electrons. In a cell, mitochondria extract the energy contained in the nutrients used by the cell, as well as doing many other specialized tasks. Each human mitochondrion has a chromosome composed of mitochondrial DNA. This DNA is distinct from the DNA that is located in the cell's nucleus.

• The ribosomes are a large complex of RNA (Ribonucleic Acid) and protein molecules, and they are essential in the production of proteins. Those proteins will determine the behaviour of the cell.

1. M Figure 1.1 -Illustration of a the structure of a cell including: the nucleolus which is inside the nucleus, the ribosomes, a mitochondrion, the cilia, the lysosome, the centrioles, the microtubules, the golgi smooth endoplasmic reticulum, rough endoplamsic reticulum (cell membrane), the cytoplasm, the chromatin. Image taken from [1].

• The cell nucleus is the most noticeable structure in an eukaryotic cell. The nucleus has a spherical shape and is separated from the rest of the cell by a double membrane called the nuclear envelope. It is the location of the chromosomes (the structures that contain the genome), and it is responsible for maintaining them and controlling the activities of the cell by regulating the gene expression (i.e. the protein production). The nucleolus is a specialized region within the nucleus where ribosome are assembled. It is where the transcription process takes place (Section 1.1.3).

The nucleus of human cells (except for notable exceptions such as sperm cells) contains two sets of 23 chromosomes (22 regular ones and X or Y sex chromosome), each set inherited from each parent. Each chromosome is a very long DNA molecule that is wrapped tightly around proteins. Thus, chromosomes are structures of signi cant size (it can be seen under a light microscope). Di erences in size and composition allow the 24 di erent chromosomes to be distinguished from each other through an analysis called a karyotype. Near the center of a chromosome is the centromere which is a narrow region that divides the chromosome into a long arm and a short arm. Chromosomes contain the information necessary to build a human being which is called the genome. It is often put in opposition with the phenotype which represents the external features of an organism.

DNA

DNA is made of a pair of molecules forming an helix structure. The molecules are held together by weak hydrogen bonds between base pairs of nucleotides. The double helix can be seen as an extremely long ladder twisted into a helix. The sides of the ladder are formed by a backbone of sugar and phosphate molecules, and the rungs consist of nucleotide bases weakly joined in the middle by the hydrogen bonds. The DNA chain is 2.2 to 2.6 nm wide, and one nucleotide unit is 0.33 nm long. Although each individual base is very small, the The chromosome is a long molecule with a double helix structure (right). Image taken from [1] full molecule is a combination of millions of base units and is much larger. For instance, the largest human chromosome has around 300 million base pairs (bp).

DNA is coded in an alphabet with 4 basic units, the nucleotides. A DNA nucleotide is made of a ve carbon sugar (a deoxyribose), a molecule of phosphoric acid, and a molecule called a base. The code letters for nucleotides are A, C, G, and T, which stand respectively for adenine, cytosine, guanine, and thymine. In DNA base pairing, adenine always pairs with thymine, and guanine always pairs with cytosine. Based on this, the information contained in one strand of the DNA can be inferred looking at the other strand. This information relies on the order of the nucleotides. DNA sequencing aims at retrieving this order. The genome of an organism is in fact made of the sequence of all the nucleotide bases for all the existing chromosomes.

Raw information in the genome must go through a certain chemical process before being able to be exploited by the cell. This genetic information must rst be translated into proteins (amino acid sequences) by living cells. The canonical genetic code de nes a mapping between DNA and the formation of the proteins.

From DNA to proteins

DNA is transcripted in RNA (Ribonucleic acid), a nucleic acid molecule similar to DNA.

While DNA is only one molecule, the RNA is replicated in several copies in the cell which are then interpreted for creating proteins. The amount of copies of RNA molecules changes from one cell to another. This explains that even though every cell has the same DNA, their behaviour will be di erent depending on which part of the DNA is transcripted.

Contrarily to DNA, RNA is usually single-stranded and the sugar molecule of its nucleotides is a ribose rather than a deoxyribose. Furthermore, the nucleotide base thymine (T), that 1. M is present in DNA, is replaced by uracil (U) in RNA. RNA play crucial roles in protein synthesis and other cell activities.

RNA is formed following DNA information. Synthesis of RNA is usually catalysed by the RNA polymerase which is an enzyme (a protein) that assembles the RNA from ribonucleotides. The transcription begins with the binding of the enzyme to a promoter sequence in the DNA. The DNA double helix is unwound by the enzyme and, then, the enzyme progresses along the template strand synthesizing a complementary RNA molecule. A DNA sequence indicates where RNA synthesis stops. There are more than thirty classes of RNA molecules. Among the most important ones we can cite the following three:

• Messenger RNA (mRNA) carries information from DNA to the ribosome where proteins are synthesized. The amino acid sequence in the protein is determined by the coding sequence of the mRNA.

• Transfer RNA (tRNA) is a short-chain type of RNA present in cells that carry a speci c amino acid through the cell.

• Ribosomal RNA (rRNA) forms a site for the synthesis of the proteins.

The identi cation of the role of each part of the DNA (called genes) in the coding of speci c proteins is then interesting. All humans have in common most of the genome, but the relatively few genomic di erences between individuals account for most of the di erences among individuals. There are estimated between 20 000 and 25 000 human protein-coding genes, although this number could drop. Human genes are distributed unevenly across the chromosomes: each chromosome contains various gene-rich and gene-poor regions. A gene contains both coding sequences that determine what the gene does, and non-coding sequences that determine when the gene is expressed (i.e. active) which have a crucial role in the control of the expression.

The most basic unit of genetic variation is the single nucleotide polymorphism or SNP. SNPs are single base-pair changes in the DNA sequence that occur with high frequency in the human genome [2]. Despite the high number of existing SNP, only some of them have an impact on the behaviour of cells through the transcription that leads to the production of proteins [3]. Genetic variants can also be called alleles.

SNPs exist in a variable proportion in human populations. The rst information that is taken for a SNP is it's Minor Allele Frequency (MAF) that indicates whether the alleles know for that SNP are common or not. For instance, a SNP with a MAF of 0.30 implies that 30% of a population has the allele versus the more common allele (the major allele) which is found in 70% of the population (assuming that it has only 2 variants, which is the most common situation). An allele is called "rare" below 0.5% MAF, "low frequency" between 0.5% and 5% and "common" above 5% [4]. Depending on the categorization, the allele can be speci c to small or large populations and the statistical tools used to assess an association can di er.

When a SNP is scarse in the population (i.e. rare variant), the SNP can be called a mutation. Cystic brosis is an example involving those rare variants. Results have shown that certain mutations show an association with a disease [5]. In those cases, the SNPs have a strong role in the biological process. Unfortunately, data on common variants can be easier to obtain and analyze, but SNP-disease association can be hard to nd [6].

Computational biology disciplines

Studies on rare variants deal with small portions of the population being a ected by a genetic trait. The idea of working with common variants aims at making the genetic studies bene cial for a larger portion of the population. It relies on the idea that common diseases have a di erent underlying genetic architecture than rare disorders which led to the development of the common disease/common variant (CD/CV) hypothesis [7].

The concept of linkage disequilibrium (LD) is commonly used for the correlation between SNPs. " Equilibrium " stands for the fact that, when a mutation appears in a population, corresponding data will be correlated in this population in opposition with other populations.

Knowing the LD gives a priori information on the structure of the data.

Many measures of LD have been proposed, though all are ultimately related to the di erence between the observed frequency of co-occurrence for two alleles and the frequency expected if the two markers were independent. The choice of measures for LD has been discussed in the litterature [8,9].

Computational biology disciplines

With the increase in the computational capacity of computers and the amount of available biological data "bioinformatics" and "biostatistics" became key-words for various elds. Computational biology can be de ned as a new interdisciplinary eld that applies the techniques of computer science, applied mathematics and statistics to address biological problems. Being a very general de nition, it encompasses many di erent elds. Here is a list of di erent approaches that have been taken.

• Computational biomodeling builds computational models of biological systems. These biomodels try to mathematically emulate the biological mechanisms involved in a particular system.

• Computational genomics is a eld that studies the genomes of the organisms. It aims at understanding the genome and, particularly, the principles of how DNA controls the biology of any species at the molecular level. This is the axis studied in this document.

• Molecular modeling similarly to the computational biomodeling, models the behaviour of molecules (hence in a smaller scale). The elds of application can range from small chemical systems to large biological molecules. Typical examples are potential or energy functions that simulate biological behaviours by mathematics.

• Protein structure prediction aims at discovering the shape of the proteins, how it folds in di erent environments. The observation of the 3D structure of the proteins being extremely expensive, mathematical models are built to guess the possible shapes.

• Computational biochemistry makes extensive use of structural modeling and simulation methods in an attempt to elucidate the kinetics and thermodynamics of protein functions. Di erential equations without a close solution and computational optimization are used to adjust a kinetic model corresponding to the real measures obtained by experimental methods.

1. M

Genetic epidemiology

Among the previously listed domain of computational biology, genomics focuses on DNA analysis. After decades of research, several biological mechanisms have been highlighted and di erent statistical approaches have been developed. Genetic epidemiology aims at understanding the role of genetic factors in the development of diseases. The issue consists in knowing which genetic factors are associated with the onset of a particular illness. For this purpose it is necessary to compare the genotype of people with the phenotype traits in question.

Since past decades, data analysis applied to high dimensional data has become essential, especially in biostatistics. Extracting information from ever larger data has become a trend in numerous elds. The key word "Big Data" has arisen to describe this kind of analysis which covers a variety of elds such as experimental data in physics (plasma data, astronomy data), or internet data (data from Facebook, emails). Medicine and biology are also involved in the increasing amount of data generated [10]. While the existence of the genome was discovered early in the 20th century, the ability to extract all the genome information on consequent amount of population have been quite recent. In the 70's extraction of sequences of the genome was performed with extremely expensive methods and huge means. Nowadays, a whole genome can be sequenced within the day at an a ordable cost (less than 1 000 dollars); technological leaps have permitted to gather data at a much higher pace.

This amount of data have permited to explain more disease/trait heritability and have led to advances in genetic epidemiology [11].

Because of this need for managing ever-growing data the development of methods adapted for high-dimensional data is rising. The problem concerns not only the volume of data but also the variety of data.

Statistical focuses

From a statistical point of view, data are collected in a data set. Observations are performed on variables. Some variables are related to the genome and are predictors. Other variables (also called "traits") are related to the phenotype which is related to a disease in genetic epidemiology, they are outcome that must be explained from the predictors. Additional variables can exist and help building models and they are called covariates.

Variables can be either quantitative or qualitative. Qualitative variables represent types and are not countable, where as quantitative variables are the results of a measure. Models are inferred on the data by means of a given method or algorithm. They can either explain the data or being able to make predictions about the outcome when predictors are observed for a new case. When outcome data is qualitative, the prediction is called classi cation.

Association Studies

A rst area of investigation consists in identifying the association between variants of SNPs. It permits to characterize a population by establishing correlation between variants. The study of this correlation leads to establishing of linkage disequilibrium information which represents the correlation between SNPs.

Genetic epidemiology

At the larger scale of association studies, Genome Wide Association Studies (GWAS) analyse DNA sequence variations at the genome level.

The rst sequencing of a whole genome was achieved in 2003 [12]. Since then, the rst Genome Wide Association Study, was conducted in 2005. 96 patients with age-related macular degeneration (ARMD) were compared with 50 healthy control individuals [13].

It was the rst time the whole genome of patients was gathered on a population and the study identi ed two SNPs with signi cantly altered allele frequency from one group to the other. The establishing of a link from a statistical point of view was followed by the establishing of a link from a biological point of view.

A major motivation on GWAS is to highlight DNA sequences associated to a phenotype with the hope of revealing new biological mechanisms. Another motivation is providing models that help for the medication of patient after a genetic test. For instance, after the discovery of the in uence of DNA sequences on warfarin dosing for blood treatments [14], a genetic test was proposed. Using genetic information from the patients to adapt and optimize their treatment is one of the main goals of personalized medicine.

Gene-Environment Interaction

Another approach consists in considering the gene-environment interaction: the exposition to a certain environment can act as a catalyst for the expression of a given genetic feature.

Identifying such an interaction allow to highlight patient that would be more likely to have a certain disease under the exposition to a given environment. From a statistical point of view, the identi cation of interaction e ect requires speci c statistic tools. The problem is discussed in further details in Chapter 3.

Multi-level data analysis

As it will be explained in Section 1.3.7, data can be related to either the genome, the transcription, or the proteins. Most of the analyses in genetic focus on one type of data, but some approaches build models taking into account several of those type of data. In the context of this manuscript, only one type of data (genomic) is considered.

Meta-analysis

Meta-analysis is the study of data from a given omic data type but coming from di erent datasets. In this case biases between them cannot be quanti ed. In this case, separate models can be applied on each source, but the comparison of the results and the formulation of global model need a particular work. Some issues can appear like batch e ects due to biases. The matter is explained in more details in Chapter 4.

Pleiotropy

Pleiotropy is the eld that studies the case where one or more genetic factor have an e ect on one or more phenotypes [15]. The problem can be seen either as a multivariate analysis or a meta-analysis problem and, thus, it can be tackled with several statistical approaches such as multivariate models, set based models or meta-analysis models. The subject is discussed in further detail in Chapter 5.

1. M

Data sources

A recent English neologism, omics, has been used to address biology elds that study cell molecules in their entirety. The rst example is genomics, which stands for the study of the whole genome. In general, the term omics focuses on large scale and holistic data research to understand information contained in cells. Genomic data are the focus of this document and a list of the current database sources is given below. The study of genomic data is called genomics.

There are other types of omic worth of mentioning. Due to the success of high-throughput biological devices and new analytical tools, the su x -om-has also been picked up by a wide array of other large-scale quantitative biology elds. Some of them are very recent terms that are not fully accepted by all the research community. Some of the most established ones are presented in the following.

Genomics

At rst DNA data was sequenced pair by pair [16]. Since then, the methods of acquisition have been improved. Nowadays, two ways of acquisition of DNA data are mainly used: Microarrays and deep sequencing technologies (NGS). Microarrays are based on the hybridization of DNA [17]. Genetic information is retrieved through probes that correspond to prede ned location in the genome. NGS (next-generation sequencing) relies on the fragmentation of DNA into pieces that are then sequenced and assembled to recontruct the complete sequence [18]. Microarrays can retrieve information faster to process but also need an a priori knowledge of locations in the genome while NGS is longer to process but can nd new sequences.

The collection at genome wide level of DNA sequencing have been lead by consortia of institutions and laboratories. A list of the main ones is:

• The international HapMap Project established a catalog of SNPs mostly in European populations but also in Africa (the Yoruba population), in China (the Han population) and in Japan (Tokyo) [8], [19]. To increase the number of SNPs new projects have been lead and the project has been reaching 11 populations on its last version [4]. This collection of data has been made possible thanks to the edi cation of a nomenclature for the human genome project [20].

• The 1000 Genomes Project was concluded in 2015, with 2504 genomes sequenced from 26 populations [21].

• The UK 10K project propose an even larger number of observations [22].

• The 100k Genomes Project is one of the latest project [23].

• Precision Medicine Initiative is another of the latest project [24].

These databases allow to de ne the genetic variation in large populations. In genetic epidemiology, the genome of individuals showing given diseases are the main focus. The problem is more speci c and hence, databases have less observations, and can show a reduced number of variables. Despite not being necessarily genome-wide studies, those genome studies are still subject to the same issues in the analysis.

Conclusion

Transcriptomic data

Transcriptomcis deals with the recollection of sequences of mRNA. Data on this kind of information can come from either microarrays or RNA sequencing and are analyzed through di erent steps. In microarrays, sequences of interest in the mRNA are a priori known. They are isolated through an array of probes, each probe corresponds to a feature that have been determined a priori. Data are compiled recollecting information from each probe. In RNA sequencing the full sequence of ribonucleic acid is obtained and is analysed.

From one side, RNA sequencing allow to discover new RNA markers. From the other side, microarray can recollect data at a higher pace. Both approaches are used in a back and forth process where RNA sequencing help building a priori knowledge for making RNA sequencing.

We can note that genomic data and transcriptomic data have similarities: they have a large number of variables, each variable being related to a gene. This is the reason why statistical methods used for one of the type of data is also used for the other. They give an insight into the role of a gene in the behavior of a cell. Especially, it can give an apriori knowledge of the parts of the DNA that are most likely to have an impact on a cell (remember most of the DNA is not transcripted [25]). eQTLs (Expression quantitative trait loci) identify loci related to gene expressions. A large amount of loci has already been determined [26].

Proteomic data

Proteomic data are related to the proteins in the cell. The sequencing of the proteins on a basis of 20 amino acids is determined. Then, the shape of the molecule, i.e. how the molecule is folded, is inferred. This shape determines the function of the protein [27]. The comparison of proteins present within each cell permits to understand the dynamics in the cell. This information can be used for diagnostic but also for designing drugs with speci c targets [27]. The eld has its own databases like [28] and protein-protein interaction can be studied similarly to SNP data [29,30].

Conclusion

Genetic studies aim at understanding the mechanisms ruling the behaviour of the cell. Genetic epidemiology raises hope for new medical approaches for treating diseases with preventive medicine or risk assessments in the application of drugs. The increase of results in the elds is due to the recent arrival of a deluge of new data. In terms of statistics presented, new challenges are raised by this data, advancing in these challenges is the motivation of this thesis. Especially, in the following chapters, studied data will be genomic data and the aim is to deal with genetic epidemiology.

C

2

Statistical methods for genetic epidemiology

Genetic epidemiology aims at explaining the in uence of the genome on the onset of a disease. It leads to studying the association between a genotype and a phenotype, i.e. predictors and an outcome. In this chapter statistical methods aiming at analyzing such associations are presented. First the mathematical notations for biological data are introduced. Secondly the challenges of the analysis are developed. Finally, statistical methods for nding simple e ects in genomics are presented.

Mathematical notations

Data are represented by X ∈ R n×p and Y ∈ R n×q , two matrices, representing n observations of p predictors and q independent variables. For any matrix A of size (a, b), for i ∈ {1, . . . , a}, its rows are noted A i,• and, for j ∈ {1, . . . , b}, its columns are noted A •,j .

For subsets ã ⊂ {1, . . . , a} and b ⊂ {1, . . . , b} resp. row and column sub-matrices are noted A ã,• and A •, b. For any vector v of size a, for i ∈ {1, . . . , a}, its elements are noted v i and for subsets ã ⊂ {1, . . . , a}, v ã represents the elements of the vector corresponding to the positions in the subset.

The Frobenius norm on matrices is denoted F . We note X T the transpose matrix of X and the cardinal of a set S is noted #S. The positive value of a real number x is noted and is equal to the number if it is positive and equal to zero otherwise.

(x) + = |x|+x
The notations for submatrices permit to represent a data matrix considering sets of observations and/or group of variables. Let us consider M di erent sets of observations in the data. Noting, for m ∈ N, M m a subset of {1, . . . , n}, let M = (M m ) m=1..M be a partition of {1, ..., n} corresponding to the observations sets. We note #M m = n m . Row blocks are de ned by this partition and a representation is given in Figure 2.1 (observations are assumed to be ordered by observation set).

Let us consider that the variables are gathered in K groups. Let P = (P k ) k=1..K be a partition of {1, ..., p} corresponding to this variable group structure. We note #P k = p k . We then we have K k=1 p k = p. Column blocks are de ned by this partitions as presented in Figure 2.2 (observations are assumed to be ordered by observation set). In general, observation sets can represent the fact that di erent sets of observations come from di erent sources and must be analyzed accordingly. For instance, data coming from 2.2. Statistical challenges for genomics analysis Figure 2.3 -Illustration of data structured by groups of variables and groups observations. Observations and variables are assumed to be ordered by resp. observations sets and variable groups. p represents the number of variables of matrix X, q the number of variables of matrix Y and n is the number of observations. n 1 , • • • , n M are resp. the number of observations of each observation set, p 1 , • • • , p K , are the resp. the number of variables in each group of variables. di erent studies may present biases. Group of variables can represent either a set of variables that is known to be quite correlated or a group of variables that must be treated together. For instance, in genetics, a locus de nes a group of SNP variables.

Statistical challenges for genomics analysis

Due to the nature of the biological process and the method of acquisition of data, SNP data present particular statistical challenges. The main challenges are discussed in this section.

Ill-posed problems

In genomic data, the number of variables is often much bigger than the number of observations. Indeed, the number of variables is extremely large due to the length of the studied molecule while the number of observations remains smaller due to the cost of their acquisition. A large portion of the existing data analysis tools are not suited for this kind of data which is the reason why, in mathematics, such a case is called an "ill-posed" problem. As an example, the limitation of linear regression in this case is presented.

Let us consider the following linear model :

Y = Xβ + ε (2.1)
where Y is a variable to explain, X are the predictors β are the coe cent of the linear regression and ε is the noise term which is assumed to follow a normal distribution.

S

The estimated coe cients of the linear regression β are chosen in order to minimize the residuals.

β = argmin β Y -Xβ 2 (2.2)
The solution is

β = (X T X) -1 X T Y (2.3)
and then the predictive outcome is

Ŷ = X(X T X) -1 X T Y (2.4)
However such a model requires the inversion of the matrix X T X which cannot be calculated when the number of variables is larger than the number of observations.

In the following, each method used copes with the problem of p n.

Correlated variables

In genomics, a large number of variables can be correlated. Notably, it is known that genes associated to a same gene will be highly correlated. In statistics, correlation can have a huge impact in the accuracy of the methods. Having correlated data imply that similar information belongs in several variables. Hence, it a ects the performance of statistical methods for which the lesser variables exist, the better they behave. Furthermore, the addition of variables makes the results harder to interpret as several variables will have synonymous roles in the model. After years of analysis of the genome, sets of SNPs or genes have been highlighted. The location of a gene or a SNP of interest is often called locus (loci for plural). A set of genes that are known to be involved in a same biological process is called a pathway. This common representation leads to consider group of SNPs in general, as SNPs from a same gene are often highly correlated. We can gain power and reduce the number of tests by combining weak signals from SNP-level analysis. This is why in genomics, making inference at a group of variables is preferred. Over the recent years, numerous pathway analysis for GWAS data have been proposed in the literature for nding pathways associated with the studied disease [12,31].

Two approaches compete for pathway analyses. From one side, self-contained GSA methods test whether genes within a pathway are collectively associated to the phenotype of interest [32]. Most widely used self-contained GSA methods combine either SNP-level or gene-level p-values into pathway test statistics, where their signi cance is typically assessed through permutation testing. From the other side pathway analysis has also been shown in enriching GWAS results by identifying SNPs that were missed by single-SNP analyses [START_REF] Evangelou | A method for genebased pathway analysis using genomewide association study summary statistics reveals nine new type 1 diabetes associations[END_REF][START_REF] Carbonetto | Integrated enrichment analysis of variants and pathways in genome-wide association studies indicates central role for il-2 signaling genes in type 1 diabetes, and cytokine signaling genes in crohn's disease[END_REF][START_REF] Eleftherohorinou | Pathway analysis of gwas provides new insights into genetic susceptibility to 3 in ammatory diseases[END_REF].

Variety of type of variables

While the rst two challenges where related to the nature of the predictors, another challenge concerns the de nition of the outcome. In genomic studies, phenotype of interest can either be quantitative or qualitative and statistical method exist for each case. When 2.3. Methods for genetic epidemiology a physical feature can be measured in an organism, the outcome is quantitative. When an expert performs a diagnostic on an organism, a trial is set and the variable is then qualitative. However, statistical methods must be adapted to the nature of the outcome.

Some tricks can allow to pass from one case to another. For instance, threshold can be set on a quantitative outcome to create categories, but a biological incentive is needed. For the contrary case, a qualitative variable with r di erent possible values can be seen as r dummy variables which are equal to 0 or 1 and that are then considered as quantitative variable between 0 and 1. However, a loss of power in the method can arise due to this transformation.

When a method is developed for genomics, it is important to be aware to which case it can be applied.

Methods for genetic epidemiology

In this section, a list of methods that have been developed for genetic epidemiology is presented.

Group of variables structure

As we have seen in Section 1.1.3, certain SNPs tend to be highly correlated due to their being transmitted together. Moreover, a single gene typically contains several SNPs and several genes take part in the same pathway. These facts can be used to reduce the complexity of the resulting models, incorporating the information in the methods. One way to achieve this is by considering a group structure for the variables, as shown in Figure 2.2.

The approaches dealing with this structure can be classi ed into two groups. A rst group of approaches consist in summarizing the SNP data to a gene-level data or a pathway level data in a pre-process step. The obtained data are supposed to overcome the challenges listed bellow and a second step analysis can be performed on data showing less intrinsic correlation and the number of variables is reduced. A second group of approaches uses this group structure and integrates into the framework of a statistic tool.

In the rst set of approaches requires, there is a loss of information in the pre-process and a well-established a priori biological knowledge is required. Conversely, the second set of approaches, keeps all information and is then less stringent on the a priori assumptions and less in uenced by a pre-preocessing step. Its drawback is that integrating the a priori information to the framework needs investment in the mathematical analysis and in the coding of the implementation.

Pre-processing methods

As a rst simple approach, a pre-processing can be performed. A rst pre-processing method consists in summurizing closely related variables in a weighted sum. The idea is that single variables may have a low e ect on a disease but the overall e ect may be more noticeable. The use of the weighted sum principle has been proposed on linear model [START_REF] Li | Atom: a powerful gene-based association test by combining optimally weighted markers[END_REF] and on Fourier transform analysis [START_REF] Wang | Improved power by use of a weighted score test for linkage disequilibrium mapping[END_REF]. Other methods called screening methods, perform low time cost tests on variablse to identify variables with higher chances of being relevant for association. A pre-selection of variables is made which can be followed by a second 2. S step analysis on data with a reduced number of variables [START_REF] Murcray | Gene-environment interaction in genome-wide association studies[END_REF][START_REF] Kooperberg | Increasing the power of identifying gene× gene interactions in genome-wide association studies[END_REF][START_REF] Murcray | Sample size requirements to detect gene-environment interactions in genomewide association studies[END_REF][START_REF] Hsu | Powerful cocktail methods for detecting genome-wide gene-environment interaction[END_REF]. The methods presented in the following integrate the group structure to their framework.

Dimension reduction methods

Dimension reduction techniques is a set of popular methods that have been used in genetic studies, being the most common ones Canonical Correlation Analysis (CCA) [START_REF] Jose A Seoane | Canonical correlation analysis for gene-based pleiotropy discovery[END_REF], Principal Component Analysis (PCA) [START_REF] Alkes L Price | Principal components analysis corrects for strati cation in genome-wide association studies[END_REF] and Partial Least Square (PLS) [START_REF] Wang | A partial least-square approach for modeling gene-gene and gene-environment interactions when multiple markers are genotyped[END_REF]. All these methods rely on the projection of the data into a subspace of lower dimension which represents most of the variation of the data, and they are often posed as an eigen value problem [START_REF] De Bie | Eigenproblems in pattern recognition[END_REF]. PLS and CCA are both supervised models and di er from the norm used whereas PCA is unsupervised. PCA was developed in the early days of statistics [START_REF] Pearson | Liii. on lines and planes of closest t to systems of points in space[END_REF] and have been raised interest in last decades since the number of available variables has become larger and larger. It has been one of the rst methods to be used to study GWAS . PLS was introduced later by H. Wold [START_REF] Wold | Path models with latent variables: The nipals approach[END_REF]. It has had succes in domains with a large number of variables to analyze, like Chemometrics [START_REF] Lldiko | A statistical view of some chemometrics regression tools[END_REF], and GWAS [START_REF] Wang | A partial least-square approach for modeling gene-gene and gene-environment interactions when multiple markers are genotyped[END_REF]. CCA was proposed by Hotelling [START_REF] Hotelling | Relations between two sets of variates[END_REF] and its application to GWAS came even later than the two previous methods [START_REF] Jose A Seoane | Canonical correlation analysis for gene-based pleiotropy discovery[END_REF]. All methods called "dimension reduction methods" derive from those ones.

Principal Component Analysis

In PCA we search for two linear regressions, u, of the variables of X that explains most of the covariance

u = argmax u 2 =1
Xu 2 (2.5)

Partial Least Square

For PLS, we search for a linear combinations, u and v, of the variables of X and Y that explain most of the covariance of X and Y :

(u, v) = argmax u 2 = v 2 =1 u T X T Y v 2 (2.6)

Canonical Correlation Analysis

The CCA is similar to PLS, we search for linear combinations, u and v of the variable of X and Y that explain most of the correlation, instead of the covariance for the PLS, between X and Y .

(u, v) = argmax u 2 = v 2 =1 u T X T Y v 2 2 Xu 2 Y v 2 (2.7)
Both PLS and CCA are closely related and are generally applied in elds with similar properties but they have di erence statistical meaning. The study of the correlation (CCA) impose a normalization of the analyzed data. In practice, for most of the analysis, data are already standardized and hence PLS and CCA give similar results.

Sparse formulations

In these dimension reduction techniques, results are formulated with new variables that are linear combinations of the original ones. These combinations can be hard to interpret due to the huge number of coe cient they represent. To cope with this problem, Lasso methods have been used [START_REF] Zou | Sparse principal component analysis[END_REF]. Lasso's penalization shrinks to zero the participation to the model of the least relevant variables. Results highlight a smaller number of variables that are easier to explain. In addition, the noise in the signal is reduced and the power of the methods is boosted. In particular sparse Partial Least Square (sPLS), has shown encouraging results [START_REF] Cao | A sparse pls for variable selection when integrating omics data[END_REF]. In the rest of the document the Lasso penalization is applied only to the weights associated to X but the Lasso penalization could also be applied to the weights associated to Y .

For PLS the sparse formulation with a Lasso penalization on u is as follows :

min ||u||=1,||v||=1
Zuv T 2 F + λP (u)

P (u) = p i=1 |u i | 1 Z = XY T (2.8)
where λ is a parameter controlling the sparsity of the model. The larger the parameters, the less the number of non-null elements in the linear combination.

Group of variables

The Lasso penalization can be adapted when groups of variables are known a priori. In genetics, incorporation of this grouping structure is becoming increasingly common due to the success of gene set enrichment analysis approaches [START_REF] Subramanian | Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression pro les[END_REF]. Using a model that takes into account this variable group structure allow to improve the performance and the readability of the results. To this purpose, the group Lasso has been proposed for the linear regression and used in genetics [START_REF] Simon | A sparse-group lasso[END_REF]. The sparse group Partial Least Square (sgPLS) adapts the penalization to the PLS [START_REF] Liquet | Group and sparse group partial least square approaches applied in genomics context[END_REF] where two overlaid Lasso penalizations translate the group structure in the Partial Least Square formulation. A structure with group and sub-groups can also be handled by this generalization with three overlaid Lasso penalizations which is called sparse group sparse PLS (sgsPLS) [START_REF] Sutton | Sparse partial least squares with group and subgroup structure[END_REF].

The application of a group Lasso for the PLS is :

{u opt , v opt } = argmin u,v Z -uv T 2 F + λP group (u) with P group (u) = K k=1 √ p k u P k 2 and Z = X T Y.
where P k refers to the notation for group of variables from section 2.1.

A group Lasso level penalization and a classical lasso penalization can be added up in a same model.

{u opt , v opt } = argmin u,v Z -uv T 2 F + λ (1 -α) P group (u) + λαP variable (u) with P group (u) = K k=1 √ p k u P k 2 , P variable (u) = p i=1 u i 2
and Z = X T Y.

S

Where α is a parameter between 0 and 1 controlling weather the group penalization or the individual variable penalization have more impact, the lower the parameter, the more group penalization is prevalent.

Case of the overlapping groups

The solution of the optimization problems presented above can be obtained through an SVD (Singular Value Decomposition) algorithm or a NIPALS algorithm [START_REF] Wold | Soft modelling by latent variables: the non-linear iterative partial least squares (nipals) approach[END_REF]. However if the groups of variables are overlapped the solution is hard to obtain. Some methods proposed to integrate the overlap in the algorithm of the solution whereas Jacob [START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF] proposed to duplicate the variables related to the overlap. Pathways group Lasso with adaptive weights accounts for the presence of overlapping pathways [START_REF] Silver | Fast identi cation of biological pathways associated with a quantitative trait using group lasso with overlaps[END_REF]. It has been used to identify pathways associated with a trait of interest and its application [START_REF] Silver | Pathways-driven sparse regression identi es pathways and genes associated with high-density lipoprotein cholesterol in two asian cohorts[END_REF]. Several other methods were also proposed for handling overlapping group structure and it is still a eld of interest [START_REF] Obozinski | Group lasso with overlaps: the latent group lasso approach[END_REF][START_REF] Yuan | E cient methods for overlapping group lasso[END_REF]. In this dissertation, the method from Jacob is preferred as duplicating do not force a modi cation of the framework of an algorithm. It is the easiest approach to implement the possible increment of the number of variables.

P-value combination methods

Methods presented above have the common feature of trying to make an overall model of the data based on an underlying structure assumption. Other methods only aim at testing associations. It is the case of the methods presented in this section. While former approaches can provide a larger interpretability of the model and formulate predictions, latter need less stringent hypotheses on the data.

The idea behind p-value combination methods is to perform statistical tests for single variables and combine them to a group of variables result. In epidemiology genetics, groups of variables are typically genes (group of SNPs) or pathways (group of genes).

Fisher's method

Fisher's method is a well-established p-value combination method that combines the pvalues from multiple statistical tests. The FM test statistic equals

FM = -2 L =1 log (p ) = -2 log L =1 p , (2.9)
where L is the number of tests combined (for example the number of SNPs within a pathway). Under null hypothesis (that the p-values are independent) FM test statistic follows a χ 2 distribution with 2L degrees of freedom.

Adaptive rank truncated product (ARTP)

The idea behind the ARTP is to truncate the highest p-values in the FM method. The only p-values left are the most signi cant ones. We denote the ordered statistics of those p-values p (1) ≤ . . . ≤ p (L) , with p ( ) being the -th smallest p-value. The original RTP statistic given by

W K = K k=1 log (p (k) ) = log K k=1 p (k) (2.10)
combines the K smallest p-values statistics of the tested pathway [START_REF] Dudbridge | Rank truncated product of p-values, with application to genomewide association scans[END_REF]. In the adaptive RTP method J di erent truncation K 1 ≤ . . . ≤ K J are investigated. Let ŝ(K j ) be the estimated p-value for W K j , (1 ≤ j ≤ J). The ARTP statistic is then de ned using the minimum p-value procedure

M inP = min 1≤j≤J ŝ(K j ). (2.11) Let p (0) 1 , . . . , p (0)
L be the p-values for each interaction test on the null hypothesis based on the observed data. B datasets under the complete null hypothesis H 0 = H 0,1 ∩ . . . ∩ H 0,L using appropriate resampling procedure are generated. Let p

j = K j i=1 log (p (b) (i) ), 0 ≤ b ≤ B, 1 ≤ j ≤ J (2.12)
is calculated for each truncation point, for both the observed data set and each of the B simulated datasets. A p-value is computed

ŝ(b) j = B b * =0 I W (b * ) j ≤ W (b) j B + 1 , 0 ≤ b ≤ B, 1 ≤ j ≤ J (2.13)
for each W j . The p-value for the ARTP statistic M inP (0) of the group of variable is estimated as

p ART P = B b * =0 I M inP (b) ≤ M inP (0) B + 1 , (2.14)
where

M inP (b) = min 1≤j≤J ŝ(b) j , 0 ≤ b ≤ B, 1 ≤ j ≤ J. (2.15)
In order to compute this minimum, several p-values are needed for the same variable. A resampling of the data can provide those p-values, it can either be, permutations, bootstraps or parametric bootstraps.

Resampling methods

Both ARTP and FM rely on appropriate resampling strategy to generate data set under the null hypothesis considered.

Yu [START_REF] Yu | Pathway analysis by adaptive combination of p-values[END_REF] proposed a permutation procedure to evaluate the signi cance level of the ARTP statistic. An alternative to the permutation procedure is to used a parametric bootstrap procedure [START_REF] Efron | An Introduction to the Bootstrap[END_REF][START_REF] Liquet | Correction of the signi cance level when attempting multiple transformations of an explanatory variable in generalized linear models[END_REF]. For permutation procedures, the number of valid permutation can be limited while their is not such a problem for parametric bootstrap [START_REF] Good | Permutation Tests: Permutation Tests: A Practical Guide to Resampling Methods for Testing Hypotheses[END_REF].

Methods that perform a gene-based pathway analysis rst combine the p-values of single-SNP analysis into gene-level test statistics (or p-values) that are subsequently summarized into pathway-level associations. The SNP-level p-values can also lead directly to a pathwaylevel analysis.

2. S

Variance component methods

Another group of methods has been popular in GWAS: the variance component methods.

A joint test is propose at a group of variables level. The originality is to test the variance instead of means or weighting to capture information even with variables in di erent direction. This family of methods use commonly a variance component test [START_REF] Lin | Variance component testing in generalised linear models with random e ects[END_REF] which has been used in genetic studies [START_REF] Michael C Wu | Rarevariant association testing for sequencing data with the sequence kernel association test[END_REF]. The burden test is a similar test that have also been used in genetics [START_REF] Li | Atom: a powerful gene-based association test by combining optimally weighted markers[END_REF].

The most popular application of these types of methods to genetics is called SKAT [START_REF] Michael C Wu | Rarevariant association testing for sequencing data with the sequence kernel association test[END_REF]. It is especially suited for dealing with the rare variants case, but it can also be used for common variants. Indeed, a formulation combining rare and common variants has been formulated [START_REF] Ionita-Laza | Sequence kernel association tests for the combined e ect of rare and common variants[END_REF]. Improvements to its power and robustness have been also developed [START_REF] Lee | Optimal uni ed approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies[END_REF][START_REF] Derkach | Robust and powerful tests for rare variants using sher's method to combine evidence of association from two or more complementary tests[END_REF].

A genotype can generally be related to the phenotype by a regression model. Let us consider a dichotomous phenotype y and a set of variables i 1 , • • • , i r . The following logistic regression model can be proposed:

logit(π i ) = α 0 + r l=1 β i l X i l (2.16)
where logit is the function x → log p 1-p and π i denotes the disease probability. X i is the predictor and β i , the regression parameter.

The regression model can be rewritten as

logit(π i ) = α 0 + β C r l=1 w i l X i l (2.17)
where β C iks a general coe cient parameter and w i l are weights that can depend on the minor allele frequency.

The burden score for testing the null hypothesis that β C = 0 is:

Q B =   n i=1 (y i -π i )   m j=1 w j X i,     2 (2.18)
where y j is the phenotype j.

The SKAT score is then:

Q S = (y -π) K(y -π) (2.19)
where K = XW W X is and n × n kernel matrix de ned by the matrix of predictors X and a diagonal matrix of weights W .

The SKAT-O score combines burden and SKAT statistics in a weighted sum.

Q ρ = ρQ B + (1 -ρ)Q S (2.20)
where ρ is a parameter that determine the participation of the burden score and the SKAT score.

The optimal Q ρ is considered for SKAT-O

Q optimal = min 0≤ρ≤1 Q ρ (2.21)

Other methods

The methods presented so far are the basis for the work in this dissertation, but are not the only approaches to the problems posed in genetic epidemiology. For frequentist side, methods are for instance c-alpha test [START_REF] Benjamin M Neale | Testing for an unusual distribution of rare variants[END_REF], the group additive regression model [START_REF] Luan | Group additive regression models for genomic data analysis[END_REF], Tukey's model [START_REF] Chatterjee | Powerful multilocus tests of genetic association in the presence of gene-gene and gene-environment interactions[END_REF] and entropy-based methods [START_REF] Zhao | An entropy-based statistic for genomewide association studies[END_REF]. Methods can also be Bayesian [START_REF] Shahbaba | A pathway analysis method for genome-wide association studies[END_REF]77] or been take from the machine learning domain [77]. We can also cite pairwise similarity based model [START_REF] Mukhopadhyay | Association tests using kernel-based measures of multi-locus genotype similarity between individuals[END_REF] and U-statistic models [START_REF] Tzeng | On the identi cation of disease mutations by the analysis of haplotype similarity and goodness of t[END_REF][START_REF] Daniel | Nonparametric tests of association of multiple genes with human disease[END_REF][START_REF] Wessel | Generalized genomic distance-based regression methodology for multilocus association analysis[END_REF].

Conclusion

Data analysis on omics data sets particular challenges in terms of statistics. Data have a low number of observations compared to the number of variables. Models must be able to extract information from this too wide variety of information and summarize it as readable results. A large number of variables are known to be correlated, which is a wide-spread problem in statistic in general and models must be resilient to it. Furthermore, variables are often gathered into groups of variables (like when di erent SNP data refer to the same gene). Models need to be able to give interpretations at group of variables results. In this manuscript, a special focus is given to dimension reduction methods and p-value combining statistics, and extensions are presented in both cases. In the rest of the document, the contributions of this dissertation to cope with these problems will be presented. The Chapter 3 is devoted to gene-environment interaction methods, while Chapter 4 and Chapter 5 deals with meta-analysis and pleiotropy. In both cases, after a brief summary of the work is presented and, then, the original papers where the contribution was published are reproduced.

Novel statistical methods for genetic epidemiology C 3

Gene-environment interaction methods

This part covers the study of gene-environment interactions (G×E) which accounts for situations in genetic epidemiology when the exposure to a certain environment catalyst the appearance of a disease. First, G×E interactions are presented. Then the statistical challenges raised by this kind of analyses are explained. After, current approaches used in the literature are presented. Finally, a novel method is presented. It is an extension of a combination test framework for a set based G×E interactions analysis, and it is the contribution of this thesis to G×E interactions studies. The methods have been highlighted by a piblication in the Journal of the French Society of Statistics [START_REF] Broc | Investigating gene-and pathway-environment interaction analysis approaches[END_REF]. The article is presented right after the chapter.

De nition of gene-environment interaction

A rst de nition of the phenomenom was given two decades ago [START_REF] Ottman | Gene-environment interaction: de nitions and study design[END_REF]. In the case of dichotomous environments and phenotypes for instance, the analysis relies on having a population with individuals i) either exposed or not to an environment ii)with di erent disease diagnostic iii) with di erent genotype information. Observations need to span each possible combination. Then the question is to know if a population exposed to the environment and with speci c genetic variants has disease appearance larger than the one that can be expected from a simple environment e ect or a simple genetic e ect. Formal genetic evidence for G×E interaction can consist in the observation that a certain exposure has di erent e ects in di erent populations or ethnic groups or in people with di erent genetically determined phenotypes.

A rst motivation, for these studies, is the idea of personalized recommendations in a public health perspective. Identifying G × E interaction helps identifying high-risk individuals if their environmental exposure is known. A second motivation is pharmacogenetics. The idea is to determine conditions enhancing the response of a patient to a treatment [START_REF] Taye | Genome-wide gene-environment study identi es glutamate receptor gene grin2a as a parkinson's disease modi er gene via interaction with co ee[END_REF].

Interactions have been highlighted for a large variety of etiologies. Cigarette smoking has been identi ed as interacting with genetic markers for cancers [START_REF] García-Closas | Nat2 slow acetylation, gstm1 null genotype, and risk of bladder cancer: results from the spanish bladder cancer study and meta-analyses[END_REF][START_REF] Rothman | A multi-stage genome-wide association study of bladder cancer identi es multiple susceptibility loci[END_REF] and alcohol consumption interacts with genes in upper aerodigestive cancers [START_REF] Hashibe | Multiple adh genes are associated with upper aerodigestive cancers[END_REF]. Environment in uence has also been noted for Parkinson disease [START_REF] Allen | The e ects of an exercise program on fall risk factors in people with parkinson's disease: a randomized controlled trial[END_REF], for example with co ee consumption [START_REF] Taye | Genome-wide gene-environment study identi es glutamate receptor gene grin2a as a parkinson's disease modi er gene via interaction with co ee[END_REF].

3. G

Statistical challenges

From a statistical point of view interactions are often related to an interaction test. However, interaction e ects are known to be harder to highlight. Especially, is has been settled that a four times larger database in terms of observations is needed for highlighting an interaction e ect compared to a simple disease/trait e ect [START_REF] Smith | The design of case-control studies: the in uence of confounding and interaction e ects[END_REF]. This means that G×E questions are expected to need more extensive studies than usual.

The design of the recollection of clinical data takes a major role in this problem. For instance, in a rst step, case-only studies were developed for G×E [START_REF] Albert | Limitations of the case-only design for identifying gene-environment interactions[END_REF]. It appears that such studies rely heavily on strong hypothesis and the performances can be hindered when the hypothesis are violated. Then, methods using a case-control study were also used [START_REF] Chatterjee | Semiparametric maximum likelihood estimation exploiting gene-environment independence in case-control studies[END_REF].

Unfortunately, case-control analysis su ers from low statistical power. In contrast, the case-only studies can be powerful in certain scenarios, although violation of the assumption of independence between the genetic and environmental factors can greatly bias the results. Some propositions to deal with the hypothesis reliance has been developed. A bayesian method trying to cope with this reliance was proposed [START_REF] Mukherjee | Exploiting gene-environment independence for analysis of case-control studies: An empirical bayes-type shrinkage estimator to trade-o between bias and e ciency[END_REF] with a trade o between the hypothesis reliance and the method e ciency. Another model combines both case-only and case-control approaches in one method [START_REF] Li | Detecting gene-environment interactions using a combined case-only and case-control approach[END_REF]. Nowadays the case-control approach is more common due to its freedom to assumptions.

Current methods

In general the global framework of statistical methods for detecting interaction can be quite similar to the simple e ect case. If the methods rely on summary statistics at single variable levels, we can just take a summary statistic suited for interaction detection instead of simple e ect detection, although interaction e ects are harder to highlight.

iSKAT and GESAT

The SKAT method presented in Section 2.3.5 has been adapted to G×E interaction studies. The Gene-Environment Set Association Test (GESAT) [START_REF] Lin | Test for interactions between a genetic marker set and environment in generalized linear models[END_REF] method adapts the SKAT framework for the analysis of common variants and the iSKAT method [START_REF] Peng-Lin Lin | Pathway analysis incorporating proteinprotein interaction networks identi ed candidate pathways for the seven common diseases[END_REF] adapt it for rare variants. The main di erence with SKAT relies on the di erence of regression models considered. For instance, for a dichotomous phenotype variable, a simple e ects can be modelled by a regression model as

logit[P (Y = 1)] = α i + β i X i (3.1)
Where logit is the function x → log p 1-p , X i refers to the i-th SNP variable and Y is the outcome.

Under the null hypothesis, the model would be

logit[P (Y = 1)] = α i (3.2)
For G × E studies, the regression model would be

logit[P (Y = 1|X i , E)] = α i + βE + β i X i + γX i × E (3.3)

Contribution to gene-environment interaction methods

and under the null hypothesis, the model would be

logit[P (Y = 1|X i , E)] = α i + βE + β i X i (3.4)

P-value combination methods

Methods for combining p-values can be used for G×E interaction. In these methods, the p-value is computed for a G × E interaction like in iSKAT, as it has been described in Section 3.3.1.

A change in the resampling framework is also needed especially for permutations. Yu [START_REF] Yu | Pathway analysis by adaptive combination of p-values[END_REF] use a permutation procedure to evaluate the signi cance level of the ARTP statistic in the context of disease-pathway association.

For the interaction model, the permutation procedure to generate data set under the complete null hypothesis xes SN P and E and permutes Y for generating data. Although the procedure seems easy to implement, the number of a valid permutations can be restricted and alter the quality of the results. The alternative is a parametric bootstrap procedure [START_REF] Efron | An Introduction to the Bootstrap[END_REF] which implies less stringent assumptions [START_REF] Good | Permutation Tests: Permutation Tests: A Practical Guide to Resampling Methods for Testing Hypotheses[END_REF][START_REF] Liquet | Correction of the signi cance level when attempting multiple transformations of an explanatory variable in generalized linear models[END_REF]. In our context, this procedure is chosen in the following.

Contribution to gene-environment interaction methods

In our contribution, we have decided to compare the Fisher's product test statistic and ARTP approach because they have been discussed in the literature to be the most powerful pathway analyses among combination tests [START_REF] Evangelou | Comparison of methods for competitive tests of pathway analysis[END_REF][START_REF] Su | Adaptive set-based methods for association testing[END_REF].

We choose to use an ARTP (Adaptative Rank Truncated Product) method. For SNPenvironment interactions, p-values are computed and a result at the gene-environment and pathway environment is inferred. The chosen method uses a null hypothesis corresponding to a self contained approach [START_REF] Evangelou | Comparison of methods for competitive tests of pathway analysis[END_REF]. In most of the p-value combination methods used in genetic epidemiology the null hypothesis is chosen to highlight simple e ect. Instead, an interaction test is performed in the presented method. Also, the framework of the ARTP requires a resampling procedure. The permutation is widely used for this purpose. However, promising results with parametric bootstrap have been shown which is interesting as the method is more adaptable to cases where the number of observations is low [START_REF] Buzkova | Permutation and parametric bootstrap tests for gene-gene and gene-environment interactions[END_REF]. The method still hasn't been used in a case of interactions. The proposed model aims at lling this gap.

The contribution has been used on real data.

The analysis is conducted in a population-based case-control study from France including 1 126 breast cancer cases and 1 174 controls. Data are composed of 23 genes gathered in one pathway: the circadian pathway. This is a case-control study conducted in Côte d'Or and Ille-et-Vilaine (France) [START_REF] Truong | Breast cancer risk, nightwork, and circadian clock gene polymorphisms[END_REF]. Eligible cases were women aged between 25 and 75 years, resident in one of these two areas and diagnosed breast cancer between 2005 and 2007. A total of 1 232 incident breast cancer cases were included in the study. Controls were selected among women living in the same areas with no history of breast cancer. In total 1 317 controls were enrolled in the study. Genotyped data from a microarray for oncology (called oncoarray) are available on 1019 cases and 999 controls. The OncoArray chip targets 3. G up to 500 000 variants with a genome-wide backbone of 250 000 tag SNPs. A preselection of SNP of interest can then be performed by epidemiologists.

The developed methods have been implemented in an R package [START_REF] Liquet | PIGE: Self Contained Gene Set Analysis for Gene-And Pathway-Environment Interaction Analysis[END_REF].

Conclusion

This part presents the G × E interaction questions in general. We can note that there are not so many statistical methods dedicated to this problem and we can wonder how to use properties of existing statistical methods in genomics and apply them to G × E interaction. This approach have been followed in this chapter where the parametric bootstrap as a resampling is used while it has only been highlighted in the literature for simple e ects but no G × E interaction e ects.

Investigation Gene-and Pathway-environment Interaction analysis approaches

Camilo Broc, Marina Evangelou, Pascal Guenel, Therese Truong and Benoit Liquet

Published in Journal de la Société Française de Statistiques

Introduction

During the last decade, genome-wide association studies (GWAS) have been successful in identifying several hundred single nucleotide polymorphisms (SNPs) associated with multiple cancer types (http://www.genome.gov/gwastudies/). However, such findings are not enough for explaining the genetic heritability of these cancers. Several reasons have been discussed that possibly explain the "missing heritability" in complex diseases such as the fact that most of these genetic associations were identified through single-SNP analyses (each SNP tested independently). It has been raised that polygenic effects, gene-gene and gene-environment (GxE) interactions are not fully explored in traditional methods [START_REF] Manolio | Functional and genomic context in pathway analysis of gwas data[END_REF]. Several approaches were developed in order to complete the agnostic GWAS in the discovery of additional genetic risk factors or to provide additional insights into the mechanisms involved in the studied disease. One such approach is pathway analysis that consists of aggregating signals from SNPs (and/or genes) to pathways. Pathways are sets of genes that work together for the production of a specific biological outcome. Pathway analysis therefore incorporates the available biological knowledge of genes and SNPs for a better understanding of the genetic and biological mechanisms of the studied disease [START_REF] Mooney | Functional and genomic context in pathway analysis of gwas data[END_REF][START_REF] Pers | Gene set analysis for interpreting genetic studies[END_REF]). One of the main thrust of the statistical analyses will be to gain power and reduce the number of tests by combining weak signals from SNP-level analysis. Over the recent years, numerous pathway analysis for GWAS data have been proposed in the literature for finding pathways associated with the studied disease (a non exclusive list includes the methods proposed by Wang These methods can be divided by the null hypothesis they test, namely the competitive (enrichment) or self-contained (association) null hypotheses. The self-contained null hypothesis states that no pathway genes are associated with the phenotype. On the other hand, the competitive hypothesis states that the statistics of the pathway genes are no more associated with the phenotype than the statistics of the genes outside of the pathway. A pathway where the competitive null hypothesis is rejected, is said to be an enriched one. The self-contained null hypothesis can be tested in both GWAS and candidate gene analysis, since only the statistics from a selection of genes is required. By contrast, competitive methods are usually used in GWAS data as all pathways are tested simultaneously. As discussed in the literature self-contained methods are generally more powerful than competitive methods [START_REF] Evangelou | Comparison of methods for competitive tests of pathway analysis[END_REF].

Self-contained methods could also be classified into marginal approaches, which are based on the combination of p-values of individual SNPs (such as for instance Fisher's Method (FM), Adaptative Rank Truncated Product (ARTP)), or joint approaches that jointly model and test the effect of all the SNPs in the set (such as random and mixed effect models, Sequence Kernel association Test (SKAT) proposed by Ionita-Laza et al., 2013). p-value combination test statistics are usually combined with phenotype permutations for estimating their significance. Phenotype permutations avoid making any assumptions about the distribution of the effect of the genetic variants on the disease.

Another distinction factor is whether the pathway analysis is considered at the gene-level or at the SNP-level. Methods that perform a gene-based pathway analysis they first combine the p-values of single-SNP analysis into gene-level test statistics (or p-values) that are subsequently summarised into pathway-level associations. On the other hand, SNP level pathway methods skip the intermediate gene level and map SNPs directly to pathways.

Although a number of statistical approaches have been proposed to test for pathway association with disease, the literature has not been greatly extended for testing for GxE interactions at the pathway level. [START_REF] Lin | Test for interactions between a genetic marker set and environment in generalized linear models[END_REF] proposed a computationally efficient GxE set association test (GESAT), a variance component score test statistic is proposed that extends the SNP-set Kernel association test for GxE testing. The proposed method tests each set of SNPs independently from the other sets and it is a SNP-based pathway analysis approach.

In addition to this, [START_REF] Jiao | Sberia: set-based gene-environment interaction test for rare and common variants in complex diseases[END_REF] proposed the set based gene environment interaction (SBERIA) method and two more extensions that overcome the limitations of SBERIA [START_REF] Jiao | Powerful set-based gene-environment interaction testing framework for complex diseases[END_REF] for both rare and common variants. SBERIA firstly computes the correlation between the environmental factors and all SNPs in the set, where a z-score of correlation is obtained. These scores are translated into weights based on a preselected threshold that are included in a regression model that tests whether they are needed or not in the model. The first extension, named enhanced set-based GxE testing (eSBERIA), is composed of two steps: the first one tests the null hypothesis that the gene-environment weights are not associated with the response. The second step implements the SKAT statistic that accounts for any residual effects that might have been missed by the logistic regression model with the gene-environment interactions. As the two tests are independent, their p-values are combined using Fisher's product statistic. The third proposed approach is coSBERIA which combines SBERIA and SKAT tests for the case-only test. The case-only GxE test for a single SNP has been found to improve the power for testing for GxE under the assumption that G and E are independent [START_REF] Albert | Limitations of the case-only design for identifying geneenvironment interactions[END_REF].

In this work we were interested in extension of combination tests to the analysis of set based GxE interactions for which we have looked at replacing the phenotype permutation procedure for testing the significance of each pathway (and/or gene) by the bootstrap approach proposed by [START_REF] Buzkova | Permutation and parametric bootstrap tests for gene-gene and geneenvironment interactions[END_REF]. We have decided to compare the Fisher's product test statistic and ARTP approach as these two approaches have been discussed in the literature to be the most powerful pathway analyses among combination tests [START_REF] Evangelou | Comparison of methods for competitive tests of pathway analysis[END_REF][START_REF] Su | Adaptive set-based methods for association testing[END_REF]. In contrast to the other proposed approaches GESAT and SBERIA, we are considering the case of genebased pathway analyses over SNP-based ones. Further, [START_REF] Su | Adaptive set-based methods for association testing[END_REF] discussed the need for a fast algorithm to test for GxE interactions through pathway analysis and in this conducted work we aim to fill this gap.

A brief description of the context which has first motivated the development of our R package PIGE (Pathway Interaction Gene Environment, Liquet et al., 2017) is presented in Section 2. In Section 3, Fisher's Method (FM) and ARTP approaches are presented in context of gene and pathway-environment interaction. Both permutation and parametric bootstrap resampling methods are presented. In Section 4, a simulation study is presented to analyse the performances of FM and ARTP methods combined with both permutation and parametric bootstrap approaches. The methods are applied on genotype data from the CECILE case-control study in Section 5. Concluding remarks are presented in Section 6.

Motivation

In a first step, we propose to test the performance of the proposed methods in simulated datasets that mimic the application dataset. In a second step, the methods will be illustrated on genotype data from the CECILE case-control study, in which we are interested in the interaction between nightwork, a binary environmental factor (defined as ever worked at night more than two years: yes/no) and polymorphisms from genes in the circadian rhythm pathway [START_REF] Truong | Breast cancer risk, nightwork and circadian clock gene polymorphisms[END_REF].

In the case of a binary response Y , the null hypothesis that there is no association between the response and the interaction term between SNP and environment is evaluated through the following logistic model:

logit[P(Y = 1|SNP , E)] = α + β SNP + β E, E + γ E × SNP . ( 1
)
where E presents the environmental factor. The likelihood ratio test (LRT) could be used to test the evidence of the interaction term (H 0, : γ = 0 versus H 1, : γ = 0), resulting to the -th p-value (p ). As discussed earlier, there are usually multiple SNPs within each gene and multiple genes within each pathway. The questions that we will answer through our conducted work are: how to combine these results to get (i) association evidence between gene-environment interaction and the outcome, (ii) association evidence between pathway-environment interaction and the outcome?

In the context of a gene-based pathway analysis, a two-step procedure is needed. At the first level the association evidence between a gene and the response is found and at the second level these gene-level p-values are combined into a test statistic for the disease-pathway association.

Phenotype permutations are usually implemented for computing the null distribution of the test statistic that can be used for obtaining a p-value for the global null hypothesis of no association between the gene with the response. In this work, we are investigating the performance of two alternative resampling approaches one based on phenotype permutations and a second one on the bootstrap approach proposed by [START_REF] Buzkova | Permutation and parametric bootstrap tests for gene-gene and geneenvironment interactions[END_REF] that has been proposed for interaction models. Both these resampling approaches are presented in Section 3.3.

Methods

In this section, we first present two frequentist approaches for combining p-values under investigation FM, and ARTP methods. We subsequently present the two resampling approaches. Finally, we shortly present an alternative frequentist approach iSKAT.

Fisher's method

Fisher's method is a well established association method that combines the results from multiple statistical tests. The FM test statistic equals

FM = -2 L ∑ =1 log (p ) = -2 log L ∏ =1 p , (2)
where L is for example the number of SNPs within a pathway. Under null hypothesis the FM test statistic follows a χ 2 distribution with 2L degrees of freedom when the p-values are independent.

In the presence of linkage disequilibrium, the correlation between SNPs leads to dependent teststatistics. We have used the resampling approaches (presented in Section 3.3) to approximate the empirical distribution of the FM test statistics.

Adaptive rank truncated product (ARTP)

The idea behind the ARTP is to truncate the highest p-values in the FM method. The only pvalues left are the most significant ones. To simplify the presentation of the ARTP proposed by Yu et al. ( 2009), we consider a pathway consisting of L SNPs and we want to test the null hypothesis that there is no pathway-environment interaction associated to the disease phenotype. Using model ( 1), we can perform a LRT test on individual interaction E×SNPs within the considered pathway. We denote the ordered statistics of those p-values p (1) ≤ . . . ≤ p (L) , with p ( ) being the -th smallest p-value. The original RTP statistic given by

W K = K ∑ k=1 log (p (k) ) = log K ∏ k=1 p (k) (3)
combines the K smallest p-values E×SNP statistics of the tested pathway (Dudbridge and Koeleman, 2003). In the adaptive RTP method J different truncation K 1 ≤ . . . ≤ K J are investigated. Let ŝ(K j ) be the estimated p-value for W K j , (1 ≤ j ≤ J). The ARTP statistic is then defined using the minimum p-value procedure

MinP = min 1≤ j≤J ŝ(K j ). (4)
Note that for a single truncation point (J = 1), the ARTP method is the RTP method and the RTP statistic simplifies to the FM test statistic when the truncation point K is fixed to L. Two levels of resampling approach are required to get the adjusted p-value for MinP: (1) for estimating ŝ(K j ), (2) for the adjustment for multiple testing over different truncation points. To avoid this computational issue specially when the number of test L is large,

L be the p-values for each interaction test on the null hypothesis based on the observed data. We generate B datasets under the complete null hypothesis H 0 = H 0,1 ∩ . . . H 0,L using appropriate resampling procedure (see section 3.3). Let p

W (b) j = K j ∑ i=1 log (p (b) (i) ), 0 ≤ b ≤ B, 1 ≤ j ≤ J (5)
is calculated for each truncation point, for both the observed data-set and each of the B simulated datasets. Then Ge's algorithm is used to estimate the p-value

ŝ(b) j = ∑ B b * =0 I W (b * ) j ≤ W (b) j B + 1 , 0 ≤ b ≤ B, 1 ≤ j ≤ J (6)
for each W j . The p-value for the ARTP statistic MinP (0) of the pathway is estimated as

p ART P = ∑ B b * =0 I MinP (b) ≤ MinP (0) B + 1 , (7)
where

MinP (b) = min 1≤ j≤J ŝ(b) j , 0 ≤ b ≤ B, 1 ≤ j ≤ J. (8)
Remark. The adjusted p-value for MinP (b) , the ARTP statistic from the b-th dataset, can also be estimated similarly using

∑ B b * =0 I(MinP (b * ) ≤MinP (b) ) B+1 .
Thus this procedure can give an evidence of association between a pathway-environment interaction and the disease outcome. It is called a SNP-based strategy. We describe in the following the gene-based strategy consisting to used the ARTP method for both derive the gene-environment interaction level summary and to combine gene-environment interaction level p-values across all genes within a pathway. This procedure adapted for interaction investigation is the one described in [START_REF] Yu | Pathway analysis by adaptive combination of p-values[END_REF].

Consider a pathway composed of L genes, with the -th composed of n SNPs, 1

≤ ≤ L. Let p (0)
,i be the p-value for the association test on the i-th interaction SNP×environment of the -th gene based on the observed dataset. We then generate using resampling approach B datasets under the null hypothesis, and define p (b) , j the p-value for the test on the ith interaction SNP×environment of the -th gene based on the b-th generated dataset, 1 ≤ b ≤ B. The ARTP is then applied (with a predetermined set of candidate truncation points, which could be varied from gene to gene) to combine interaction SNPs×environment-level evidence of association within a gene. For the -th gene, we apply the minimum p-value procedure (MinP) given earlier on, 1 ≤ i ≤ n , 0 ≤ b ≤ B, to obtain p * (0) , the interaction gene×environment-level pvalue for the observed data, and p * (b) , the interaction GxE level p-value for the b-th permuted dataset. Finally in order to get a evidence of interaction pathway×environment the ARTP statistic is used to combine the gene×environment-level p-values for the observed and the resampling "null" data sets. We use the MinP procedure one more time to obtain the adjusted p-value for the pathway×environment-level ARTP statistic. Note that the same set of generated "null" datasets are exploited each time for the MinP procedure to derive interaction gene×environment-level and interaction pathway×environment-level evidence. Thus the full procedure overcomes the expensive computational multi-layer resampling issue. The same procedure is used for FM method.

Resampling methods

Both ARTP and FM rely on appropriate resampling strategy to generate data set under the null hypothesis considered. For gene-and pathways-environment interaction, we consider the global null hypothesis:

H 0 = H 0,1 ∩ . . . H 0, . . . ∩ H 0,L , with H 0, : γ = 0 (see equation (1)) ( 9
)
where L is the number of SNPs within a considered pathway.

Yu et al. ( 2009) use a permutation procedure to evaluate the significance level of the ARTP statistic in the context of disease-pathway association which for example corresponds to the situation of the simplified model:

logit[P(Y = 1|SNP , E)] = α + β SNP (10)
with H 0, : β = 0. In this situation, there is no difficulty to define the permutation procedure for the complete null hypothesis. One just need to permute the phenotype Y . However, for the interaction model (1), a valid permutation procedure to generate data set under the complete null hypothesis ( 9) is complex to define. As noted by [START_REF] Buzkova | Permutation and parametric bootstrap tests for gene-gene and geneenvironment interactions[END_REF], fixing SNP and E and permuting Y generates data in which the generated phenotype Y * is independent of SNP and E. This procedure fails to generate data set under the null hypothesis since in model ( 1) the phenotype Y is not independent of SNP and E. Indeed, for the logistic model there is no permutation procedure which can be used to generate data set for the complete null hypothesis [START_REF] Edgington | Randomization tests[END_REF]). An alternative to the permutation procedure is to used a parametric bootstrap procedure (Efron and Tibshirani, 1994; Liquet and Riou, 2013) which implies less stringent assumptions (Good, 2000). In our context, the procedure could be defined in the following.

For each SNP ( = 1, . . . , L):

1. Fit the model under the null hypothesis H 0, , using the observed data, and obtain α , β , βE, , the maximum likelihood estimate (MLE) of respectively α , β and β E,

2.

Generate a new outcome Y * i, for each subject from the probability measure defined under H 0, . For example, for model ( 1), we generate Y * i, according to:

P(Y * i, = 1|SNP , E) = exp ( α + β SNP + βE, E) 1 + exp ( α + β SNP + βE, E) .
Repeat this for all the subjects to obtain a sample noted s * = {Y * i, , SNP i, , E i } which is related to the -th SNP. Remark: In case of marginal association of both SNP and environmental factor, step 2 might generate unbalanced data which could affect the statistical power of the resampling methods. A screening investigation on the marginal association might be used before using the booststrap method.

iSKAT

The other frequentist approach is iSKAT proposed by [START_REF] Lin | Test for rare variants by environment interactions in sequencing association studies[END_REF]. The method uses the spirit of SKAT-O methods (Wu et al. ( 2011)) and apply it to an interaction test context. From one side burden tests are know to be an efficient test in many cases but they struggle when rare variants are involved in the data. From the other side kernel test can handle those rare variants. The idea behind the algorithm is to separate from the data the rare variants from the rest and to take advantage of both burden tests and kernel tests. Furthermore, iSKAT offer the possibility of weighting the covariates to take into account extra information. However, no weight have been added in the use of iSKAT in this article. The method GESAT is a particular case of the iSKAT method.

Simulation Study

In this section, we compare the FM and ARTP methods through a simulation study investigating their control of Type-I error and FWER and their power performance. Two resampling approaches (permutation and bootstrap) are compared for a range of sample sizes (n = 200, 500, 1000). The combination methods FM and ARTP are compared to iSKAT and to the popular frequentist approach MinP which combines p-value by considering only the most signifactive p-value: MinP = min ∈{1,...,L} p .

Let's note that MinP method doesn't have the same meaning than the quantity MinP used in the intermediary steps of ARTP (see equation ( 4)).

Data Generated

We work on generated data which are supposed to mimic experimental data. The parameters of the generation are inspired from [START_REF] Buzkova | Permutation and parametric bootstrap tests for gene-gene and geneenvironment interactions[END_REF]. The genetic structure simulated is composed by one pathway containing I genes (genes are called G 1 , G 2 , . . . , G I ). Each gene contains several SNPs. The SNPs are binary variables. In order to generate the i-th gene, composed by k i SNPs, SNP i 1 . . . SNP i k i , we use the following procedure:

S i ∼ Bern(0.2) logit(p j ) = logit(0.2) + S i for j ∈ {1, . . . , k i } SNP i j |S i ∼ Bern(p j ) for j ∈ {1, . . . , k i }.
Hence, conditionally on the latent polymorphism S i , for a given gene i, the individual SNP i j are independent and identically distributed, but they are marginally dependent.

A binary environment variable is also simulated that is marginally dependent with one gene i E and generated with the following procedure:

logit(p E ) = a + bS i E P[E = 1] = p E
Finally, a binary outcome variable is simulated. It is generated from a logistic model using SNPs from gene i Y . Among those s SNPs only s SNPs are associated to the response variable Y as specified in the following equations:

λ 1 , . . . , λ s ∈ {0, 1} and s ∑ l=1 λ l = s (11) logit[P(Y = 1|G i Y , E)] = α + β E E + k i Y ∑ j=1 λ j β SNP i Y j SNP i Y j + k i Y ∑ j=1 λ j γ SNP i Y j SNP i Y j × E (12)
The parameters λ l control the choice of the SNP involved in the generation of Y and the parameter s controls the number of those SNPs. Different choices for parameters β j and γ j are chosen in order to highlight different results. For Type I error results the parameters γ j are set to 0 whereas for power results parameters β j and γ j are chosen in order to evidence the different rejection of the null hypothesis for different parameters and methods used.

Simulation design

We present nine different simulation models. The first four ones are used to investigate the Type-I error and the Family Wise Error Rate (FWER) results while the others are used to explore the power of the different methods.

As discussed above, different resampling methods are used: (i) permutation that permutes the outcome Y and (ii) parametric bootstrap. We set to 1,000 the number of permutations and bootstrap resampling. The sample size n of the simulation datasets are 200, 500 and 1000. The results we look for are the p-values at gene level that are based on SNP-level information. Then a pathway-level p-values is computed from this information.

We use cases 1, 2, 3 and 4 (defined in the following) to investigate the control of the Type-I error rate for each gene. We also investigate the control of the FWER at the gene level and at the pathway level. An empirical FWER for gene-environment interaction is defined by the number of times the procedure detects wrongly at least one significant gene-environment interaction (from all the genes within the Pathway) over the N = 500 simulation replications. We also add the empirical FWER computed using a Bonferroni correction (i.e., each gene-environment interaction p-value is divided by the number of investigated genes). Indeed, the output of each method is a set of p-values (one for each investigating gene-environment interaction). This set of p-values is associated to a set of null hypotheses which define our family of hypotheses. Then it is important to control an overall error for these hypotheses. The empirical FWER for pathway-environment interaction is defined by the number of times the procedure detects wrongly a significant pathwayenvironment interaction over the 500 replications. As only one pathway is considered in our simulation study, the control of the FWER at the pathway level is similar to the Type-I error rate control of the pathway investigated.

Simulation cases for investigating Type-I error rate and FWER

-Case 1: Data is composed of 5 genes with 10 SNPs each. True model is based on the main effect of E and the main effect of 5 randomly selected SNP from the first gene. The environment is marginally correlated with the first gene but not with the other genes. The outcome is therefore generated from the following model:

λ 1 , . . . , λ 10 ∈ {0, 1} and 10 ∑ j=1 λ j = 5 logit[P(Y = 1|G 1 , E)] = α + β E E + 10 ∑ j=1 λ j β SNP 1 j SNP
λ j = 2 logit[P(Y = 1|G 1 , E)] = α + β E E + 5 ∑ j=1 λ j β SNP 1 j SNP 1 j
This model enables us to see how the methods perform when the pathway gene members have different sizes. -Case 4: Data is composed of 20 genes with 20 SNPS for the first gene and 10 SNPs for the others. True model is based on: the main effect of E; the main effect of 10 randomly selected SNPs from the first gene and 5 from the second genes. The environment is marginally correlated with first gene but not with other genes. The outcome is therefore generated from the following model:

λ 1 1 , . . . ,
∑ j=1 λ 2 j = 5 logit[P(Y = 1|G 1 ,G 2 , E)] = α + β E E + 20 ∑ j=1 λ 1 j β SNP 1 j SNP 1 j + 10 ∑ j=1 λ 2 j β SNP 2 j SNP 2 j
Simulation cases for power performance -Case 5: Data is composed of 5 genes with 20, 10, 10, 10, 10 SNPs. True model is based on the main effect of E and the main effect of 10 randomly selected SNPs from the first gene and the interaction between the environment with each of the selected SNPs. The environment is marginally correlated with first gene but not with other genes. The outcome is therefore generated from the following model:

λ 1 , . . . , λ 20 ∈ {0, 1} and 20 ∑ j=1 λ j = 10 logit[P(Y = 1|G 1 , E)] = α + β E E + 20 ∑ j=1 λ j β SNP 1 j SNP 1 j + 20 ∑ j=1 λ j γ SNP 1 j SNP 1 j × E
-Case 6: Data is composed of 5 genes with 20, 10, 10, 10, 10 SNPs. True model is based on the main effect of E and the main effect of 2 randomly selected SNPs from the first gene and the interaction between the environment with each of the selected SNPs. The environment is marginally correlated with first gene but not with other genes. The outcome is therefore generated from the following model:

λ 1 , . . . , λ 20 ∈ {0, 1} and 20 ∑ j=1 λ j = 2 logit[P(Y = 1|G 1 , E)] = α + β E E + 20 ∑ j=1 λ j β SNP 1 j SNP 1 j + 20 ∑ j=1 λ j γ SNP 1 j SNP 1 j × E -Case 7:
Data is composed of 20 genes with 20 SNPS for the first gene and 10 SNPs for the others. True model is based on: the main effect of E; the main effect of 10 randomly selected SNPs from the first gene and 5 from the second genes; the interactions between the environment with each of the selected SNPs. The environment is marginally correlated with first gene but not with other genes. The outcome is therefore generated from the following model:

λ 1 1 , . . . , λ 1 20 ∈ {0, 1} with 20 ∑ j=1 λ 1 j = 10 and λ 2 1 , . . . , λ 2 10 ∈ {0, 1} with 10 ∑ j=1 λ 2 j = 5 logit[P(Y = 1|G 1 ,G 2 , E)] = α + β E E + 20 ∑ j=1 λ 1 j β SNP 1 j SNP 1 j + 10 ∑ j=1 λ 2 j β SNP 2 j SNP 2 j + 20 ∑ j=1 λ 1 j γ SNP 1 j SNP 1 j × E + 10 ∑ j=1 λ 2 j γ SNP 2 j SNP 2 j × E
-Case 8: Data is composed of 20 genes with 20 SNPS for the first gene and 10 SNPs for the others. True model is based on: the main effect of E; the main effect of 2 randomly selected SNPs from the first gene and 2 from the second genes; the interactions between the environment with each of the selected SNPs. The environment is marginally correlated with first gene but not with other genes. The outcome is therefore generated from the following model:

λ 1 1 , . . . , λ 1 20 ∈ {0, 1} with 20 ∑ j=1 λ 1 j = 2 and λ 2 1 , . . . , λ 2 10 ∈ {0, 1} with 10 ∑ j=1 λ 2 j = 2 logit[P(Y = 1|G 1 , E)] = α + β E E + 20 ∑ j=1 λ 1 j β SNP 1 j SNP 1 j + 10 ∑ j=1 λ 2 j β SNP 2 j SNP 2 j + 20 ∑ j=1 λ 1 j γ SNP 1 j SNP 1 j × E + 10 ∑ j=1 λ 2 j γ SNP 2 j SNP 2 j × E -Case 9:
Data is composed of 2 pathways with 10 genes in each pathways. Each genes includes 10 SNPS. True model is based on: the main effect of E; the main effect of 2 randomly selected SNPs from the first gene and 2 from the second genes of each pathways;

the interactions between the environment with each of the selected SNPs. The outcome is therefore generated from the following model:

λ 1 1 , . . . ,
∑ j=1 λ 12 j = 2 logit[P(Y = 1|G 1 , E)] = α + β E E + ∑ k∈K 10 ∑ j=1 λ k j β SNP k j SNP k j + ∑ k∈K 10 ∑ j=1 λ k j γ SNP k j SNP k j × E
where K = {1, 2, 11, 12}.

Simulation parameters

The different coefficient used in our cases are gathered in table 1. The notation used refer to part 4.1. For each simulation case, FM, and ARTP methods are applied for the 9 simulation cases to investigate the presence of interaction effects of gene-and pathway-environment based on a gene-based strategy (see end of Section 3.2).

For all of the 9 cases the truncation points of the ARTP parameters are optimized like in previous ARTP results [START_REF] Yu | Pathway analysis by adaptive combination of p-values[END_REF]. Let m be an integer; k G and k SNP be real numbers; k i be the number of SNP in the i-th gene; I the number of genes. Let k G × I , 2 × k G × I , . . . , m × k G × I be a set of truncation points for genes and, for each gene i, let k SNP × k i , 2 × k SNP × k i , . . . , m × k SNP × k i be a set of truncation points for the SNPs of this gene. The notation x gives the largest integer that does not exceed x (if x = 0 we set the value to 1) . We define p k SNP ,k G ,m the p-value of the ARTP computed with this set of truncation points. The optimal pvalue of the ARTP is defined as:

min k SNP ∈A ,k G ∈A p k SNP ,k G ,m with A = {2%, 4%, . . . , 20%}.
The parameter m is fixed to 5 in our study.

Type I error rate and FWER

Type I error rate and FWER of the methods are computed in cases 1, 2, 3 and 4. The data are generated under the null hypothesis (i.e. no interaction). The SNP-level tests are performed under the interaction assumption with a significance level of 0.05. Hence, the expected value of all pvalues at gene-level and pathway level are 0.05. In the section, we study the behavior of the different methods for this case. Tables 2, 3, 4 , 5 , 6 and 7 present the results. The expected value of the average p-values is 0.05. Computing a binomial model, the average of the p-values on the 500 iterations should be between 3% and 7%.

p E a b α β E β SNP i Y j i Y ∈{1,...,k i Y } γ SNP i Y j i Y ∈{1,...,k i Y } case 1 0.2 logit(2) log(2) -2 2 ∈ {3, 2, 1} = 0 case 2 0.2 logit(2) log(2) -2 2 ∈ {1.5, 1.0, 0.5} = 0 case 3 0.2 logit(2) log(2) -2 2 ∈ {3, 2, 1} = 0 case 4 0.5 logit(2) 2 -1 1 ∈ {0.3, 0.2, 0.1} = 0 case 5 0.2 logit(2) log(2) -2 0.4 ∈ {0.06, 0.04, 0.02} = 0.5 case 6 0.2 logit(2) log(2) -2 2 ∈ {0.3, 0.2, 0.1} = 1.5 case 7 0.5 logit(2) 2 -1 0.1 ∈ {0.03, 0.02, 0.01} = 0.3 case 8 0.5 logit(2) 2 -1 0.1 ∈ {0.075, 0.050, 0.025} = 0.7 case 9 0.5 logit(2) 2 -1 0 = 0 = 0.5
The permutation approach obtained very low error rate for both approaches (FM, MinP, ARTP). The bootstrap approach gives good results for controlling the Type-I error rate for both FM and ARTP methods. As expected the FWER at the gene level is not controlled. The FWER computed can then be corrected using the Bonferroni method which is known to be conservative and more trustable. Finally, a pathway p-value is given by each combining method (FM, MinP, ARTP). When the number of genes is low (cases 1 to 3), the Type-I error rate at the pathway level is well controlled using the bootstrap approach for both FM, MinP and ARTP methods while permutation approach give conservative results. When the number of genes is higher (case 4), the type-I error rate at the pathway level of the ARTP and iSKAT is slightly inflated.

Power performance

Tables 8, 9, 10, 11 and 12 present the results for the power of the methods. As expected for all methods power performances increase with larger sample size.

In those results the proportion of significant SNP in the true model have a huge importance on the performances. In cases 5 and 7 the proportions of significant SNPs are high whereas in cases 6, 8 and 9 they are low. For higher proportions bootstrap is slightly but consistently more powerful than the permutation. For lower proportions permutation and bootstrap results are equivalent. We can see that when the proportion is high FM and ARTP have equivalent results and MinP has abysmal results. This is due to the fact that MinP truncates too much of the information contained in the data. When the proportion are lower, ARTP and then MinP have a good performance but FM have lower ones. This is due to the fact that FM results take into account too much irrelevant SNPs in its combination. The ARTP have the merit of having good power whatever is the proportion of significant SNP in the true model. FM and MinP seems to detect different patterns but ARTP can detect both. When we compare the combining methods with iSKAT we can see that the level of performance of ARTP and iSKAT is similar. We notice that in general, ARTP is more powerful on small sample sizes (200 and 500).

Application: Breast cancer and night work

Circadian rhythm is a roughly 24 hours cycle of biological processes that are synchronized by external cues such as light or temperature, and regulated endogenously by periodic transcription of a set of genes that form a network of self-regulated feedback loop. The circadian rhythm pathway plays a key role in the maintenance of various endocrine, physiological factors and behavioral functions including cell cycle regulation, hormone secretion, body temperature and sleep/wake cycle. Shift work that involves circadian disruption was classified as probably carcinogenic to humans (group 2A) by the International Agency for Research on Cancer in 2007 [START_REF] Straif | Carcinogenicity of shift-work, painting, and fire-fighting[END_REF]). An increased risk of breast cancer was reported in women working at night by several studies (Hansen and Lassen, 2012; [START_REF] Menegaux | Night work and breast cancer: a population-based case control study in france (the cecile study)[END_REF] and it was hypothesized that this association could be modulated by polymorphisms in the circadian pathway genes. As the circadian pacemaker requires multiple molecular interactions to generate the circadian rhythms, single-SNP analyses may not be sufficient to analyze the association between circadian genes and breast cancer. Therefore, we have investigated the role of circadian clock gene polymorphisms and their interaction with nightwork in breast cancer risk using a pathway analysis. This work was previously described in more details using only the ARTP method with a modified permutation procedure that permutes the outcome, the environmental factor and the adjustment variables together [START_REF] Truong | Breast cancer risk, nightwork and circadian clock gene polymorphisms[END_REF]. Here, we present the results using FM, MinP and ARTP methods using permutation and Bootstrap resampling procedures as well as iSKAT method for comparison.

Briefly, the analyses are conducted in a population-based case-control study from France including 1126 breast cancer cases and 1174 controls.

We considered the circadian pathway as defined in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database that included 23 genes (CLOCK, ARNTL, NPAS2, CRY1, CRY2, PER1, PER2, PER3, RORA, RORB, RORC, BHLHE40, BHLHE41, SKP1, FBXW11, CUL1, TIMELESS, FBXL3, NR1D1, CSNK1D, CSNK1E, RBX1, and BTRC). These genes constitute a complex regulatory network with multiple negative and positive feedback loops. A selection of tag SNPs from these genes were selected in order to capture SNPs within 5 kb of each genes (pairwise approach with a square of correlation coefficient r 2 > 0.8) with a minimum minor allele frequency of 0.05 in the CEU population from HapMap project. After quality controls, we have included 577 SNPs from the 23 genes. The circadian pathway was additionnally divided into two subpathways: the core circadian genes which are involved in the same transcriptional feedback loop (CLOCK, ARNTL, NPAS2, CRY1, CRY2, PER1, PER2, PER3,CSNK1E) and the other genes that are involved in other auxiliary loops.

Odds ratios (OR) and corresponding 95% confidence intervals (CI) were calculated using unconditional logistic regression models adjusted for the matching factors (age, area of residence) and for established risk factors of breast cancer (age at menarche, age at first full-term pregnancy, parity, current use of menopausal hormone therapy, body mass index, alcohol consumption and tobacco consumption). An OR of 1.42 (95% CI: 1.08-1.88) (p=0.01) was observed in women that have a lifetime duration of nightwork greater than 2 years compared to less. The interaction between the polymorphisms in circadian genes and nightwork were first analysed using a SNP by SNP approach and no interaction term was statistically significant after correction for multiple testing (results not shown). Gene-level and pathway-level interaction p-values obtained by the FM, MinP, ARTP and iSKAT are shown in Table 13 for 1000 resampling. The parameters of the ARTP are calibrated in the same way than in the simulation part (see section 4.2). At the gene level, we observed that both methods FM and ARTP highlight the same two genes PER1 and ARNTL in the interaction analysis with nightwork, while only PER1 is significant with the MinP method and only ARNTL is significant with the iSKAT method. Bootstrap resampling method tends to give lower p-values than permutations for these two genes in particular. This is in accordance with the simulation section in which we shown that the parametric bootstrap method is more powerful for large sample size.

At the pathway level, a significant interaction p-value (see Table 13) was observed for the overall circadian pathway for both FM and ARTP methods when parametric bootstrap is used while no association is observed using permutation resampling approach. This association is observed only for the core circadian genes subpathway that includes the genes PER1 and ARNTL. No significant association was observed while using the methods MinP and iSKAT.

To summarize, FM and ARTP gave similar results in our data. Significant interaction p-values were observed at the gene and pathway levels using the boostrap resampling method, while only significant results at the gene level were observed using the permutation resampling method. MinP and iSKAT methods highlighted only part of the genes that were found significant by FM and ARTP methods and reported non-significant interaction at the pathway level.

PER1 and ARNTL which are highlighted in the gene level analysis, are important components of the circadian system which is regulated by molecular feedback loops. Heterodimers composed of ARNTL and either of the two related proteins CLOCK or NPAS2 are transcriptional factors that induce the expression of PER and CRY genes by binding to their promoters, which in turn will act on the ARNTL-CLOCK/NPAS2 complex to repress their own transcription.

Variants in both genes has been previously associated to breast cancer risk (Hansen and Lassen, 2012; [START_REF] Zienolddiny | Analysis of polymorphisms in the circadian-related genes and breast cancer risk in norwegian nurses working night shifts[END_REF]. The finding with PER1 from the interaction analysis may be of particular interest, as a variant in PER1 (rs2735611) was previously associated with an extreme morning preference [START_REF] Carpen | A silent polymorphism in the per1 gene associates with extreme diurnal preference in humans[END_REF], a condition that was associated with an increased breast cancer risk among Danish military women working in night shifts (Hansen and Lassen, 2012).

Running time performance

The most demanding part of the p-value algorithms in terms of time computation is the resampling part. All p-value combining methods have been ran with the same resampling samples. We focus on the mesure of the running time related to this part of the algorithm. Tables 14 and15 presents the running time performances of permutation and bootstrap approaches. The results given are computed on one standard core, and results are running times on the application data (see table 14) and on 500 iterations of simulation case 5 (see table 15). We can see that bootstrap and permutation have similar running times. The running time of iSKAT is added in comparison. p-value combining methods have a much higher running time than iSKAT. Hopefully it can be computed in parallel, whereas iSKAT cannot. ). Permutation and parametric bootstrap approaches have been implemented. Our simulation study suggests sightly better results from bootstrap compared to permutation, especially when the number of significant SNP is high. Furthermore we have shown that our proposed methods can be competitive and even slightly more powerful then the cutting edge methods like iSKAT.

The cornerstone of the implemented approaches are the running time of the resampling approaches which could be problematic in presence of large data set (i.e., large sample size and large number of genetic information). To overcome this issue, PIGE offers a parallel implementation of these approaches. As an example, our application on interaction between circadian genes and night work in breast cancer risk which includes n = 2300 subjects and p = 577 SNPs took 45 minutes with the permutation procedure and 1 hours 5 minutes using 4 cores and 1000 resampling.

In this application study, using ARTP method with the parametric bootstrap approach, we highlighted significant interactions at the pathway-level which were missed when using the permutation procedures. Our results suggest that polymorphisms in the circadian rhythms pathway could modulate the association between night work and breast cancer. This association seems to be driven mostly by the genes PER1 and ARNTL.

Note that our approaches can deal in the context of p > n as the methods are based on combining individual p-values. Finally, our proposed approaches are not restricted to a binary casecontrol outcome. In this study, we focus the presentation on an binary environment variable which was motivated by binary environment data of our application. The method is not restricted to binary environment variable and has been extended and implemented in our R package PIGE [START_REF] Liquet | PIGE: Self Contained Gene Set Analysis for Gene-And Pathway-Environment Interaction Analysis[END_REF] to any quantitative environment variable. Further, our package also offers the possibility to deal with survival outcome variable or quantitative outcome in general. It is also possible to investigate gene-and pathway-environment interaction for more than one pathway during the same analysis.

Meta-analysis methods for genomics

This chapter deals with meta-analysis in the context of genetic epidemiology. Meta-analysis consists in comparing results between di erent studies. In genetic epidemiology, gathering data coming from di erent studies is an e cient way to increase the amount of data available which can enhance the power of statistical methods. In a rst section, the batch e ect problem due to the meta-analysis is explained. In a second section, methods used for meta-analysis in pleiotropy context are presented. Finally a novel method for tackling this issue is presented. The method has been presented in an article in the Arabian Journal of Mathematics which is reproduced after this chapter.

The batch e ect

For a meta analysis, data are gathered di erence sources, i.e., from di erent institutions on di erent populations and with similar but potentially di erent technologies and protocols.

In multi study data sets the accuracy of the results is impacted by the di erent experimental bias between each study. Those di erences in the recollection of data can lead to involuntary bias, known as batch e ect. In general, a proper step of harmonization of the di erent chunks of data is necessary. Several methods have been proposed to "standardize" (or normalize) data according to the sources [START_REF] Johnson | Adjusting batch e ects in microarray expression data using empirical bayes methods[END_REF][START_REF] Hornung | Combining location-and-scale batch e ect adjustment with data cleaning by latent factor adjustment[END_REF]. In general, before the application of a method each column of X and Y is centered according to its mean and scaled according to its standard deviation. As a basic approach for meta-analysis, the normalization has to be changed to take this diversity into account. For instance, the normalization can be performed study by study [START_REF] Rohart | Mint: A multivariate integrative method to identify reproducible molecular signatures across independent experiments and platforms[END_REF]. This decomposition can even be used to retrieve a study by study model from dimension reduction models like in MINT method [START_REF] Rohart | Mint: A multivariate integrative method to identify reproducible molecular signatures across independent experiments and platforms[END_REF].

Another approach consists in integrating the knowledge in the framework of the method [START_REF] Listgarten | Correction for hidden confounders in the genetic analysis of gene expression[END_REF] but this approach is less used than standardization methods. In general a meancentering method is commonly used. It is easy to implement and need few hypothesis on data [START_REF] Andrew | The removal of multiplicative, systematic bias allows integration of breast cancer gene expression datasets-improving meta-analysis and prediction of prognosis[END_REF]. More advanced methods have been discussed extensively in notable articles [START_REF] Gagnon-Bartsch | Using control genes to correct for unwanted variation in microarray data[END_REF].

In the following, we consider that study by study standardization is applied to data.

M

Methods for meta-analysis in genetic epidemiology

In this section three approaches to the meta-analysis question are presented. The rst one is ASSET which is one of the most used methods in this eld for meta-analysis. The method relies on the exploration of subsets of studies for the presence of true association signals that are in either the same direction or possibly opposite directions. CPBayes is a bayesian equivalent of the method. The second one is an extension of SKAT (see Section 2.3.5) named meta-SKAT. Finally, extensions of the Lasso methods are discussed.

A large number of other methods also exists, but are not developed in this dissertation. ASSET and meta-SKAT are presented because they are the most established methods in the eld, while the last approach sets the background for the introduction of the novel method "sparse group Partial Least Square for structured data" which is one the contribution of this thesis.

Note that the mathematical notation are introduced in Section 2.1. Especially, for this chapter, the structure of data presented in Figure 2.3 is considered.

ASSET and CPBayes

ASSET and CPBayes are methods outputting statistics for meta-analysis. Both frameworks rely on common ideas, except that ASSET is a frequentist method whereas CPBayes is Bayesian.

ASSET

ASSET is a method suited for meta-analysis providing a p-value across studies [START_REF] Bhattacharjee | A subset-based approach improves power and interpretation for the combined analysis of genetic association studies of heterogeneous traits[END_REF]. The input of the method are single variables summary statistics which are combined by the method. ASSET exhaustively explores subsets of studies for the presence of true association signals that are in either the same direction or possibly opposite directions.

For a given variable i, the considered summary statistics are β i,m and s i,m resp. the regression parameter and standard error of a model for variable i and study m. The standard xed e ect of the variable i is calculated as a weighted sum of Z statistics:

Z i = M m=1 w i,m Z i,m (4.1)
with

Z i,m = β i,m
s i,m and where w i,m is a weight calculated as

w i,m = 1/s i,k M m=1 1/w 2 i,m (4.2)
With those weights, the statistic is known to be asymptotically equivalent to a pooled analysis of the studies:

Z i = M m=1 Z i,m /s i,m 1/s i,k (4.3)

ASSET and CPBayes

In this framework, the hypothesis that some studies do not contribute to the pleiotropy e ect is considered. Hence for a given subset of study S ⊂ {1, • • • , M }, the pooled Z statistic associated is

Z i (S) = m∈S π m (S)Z i,m (4.4)
in which π m (S) = n m / m∈S n m denotes the sample size of study m compared to the subset S.

Z meta-analysis i = max possible S |Z i (S)| (4.5)
The number of possible subsets S is high: 2 K-1 but all the combinations are based on the same statistics Z i,m which reduces the computational cost.

In order to address the possibility of e ect with opposite directions a "two sided test" is proposed (whereas the statistic Z meta-analysis is called "one sided test"). Positive and negative e ect are detected separately through two statistics Z meta-analysis,+ and Z meta-analysis,-.

CPBayes

CPBayes [START_REF] Majumdar | An e cient bayesian meta-analysis approach for studying cross-phenotype genetic associations[END_REF] uses summary statistics combinations similar to ASSET but is based on a Bayesian framework with and an approach called spike and slab priors and a MCMC Gibbs sampling. The evidence of pleiotropy is measured by the local false discovery rate (locFDR) and with Bayes factors (BF). CPBayes also estimates the posterior probability of association (PPA) and coe cient for each phenotype.

We can note that neither ASSET nor CPBayes take into account a group structure of the variables. And, thus, they do not take advantage of any SNP-gene-pathway structure knowledge.

Meta-SKAT

SKAT is a method to detect association between rare variants in a region and a phenotype (continuous or binary). It is a supervised test for joint e ects of multiple variants in a region on a phenotype. Meta-SKAT can do the same but aggregating several studies. This method outputs a p-value corresponding to a set of variables, for instance a gene or a pathway. The method is based on a weighted sum of SKAT statistics of the di erent studies [109].

Lasso penalization for meta analysis on dimension reduction methods

Di erent methods have been proposed for using dimension reduction methods in the case of di erent sets of observations. They rely on simultatenous selections of variables for each study. The fuse Lasso proposed by Tibshirani [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] encourages the selected variables to be similar from one study to another. A joint group Lasso that selects the same variables for each study has also been proposed [START_REF] Obozinski | Joint covariate selection and joint subspace selection for multiple classi cation problems[END_REF] .

For each study m, u m and v m optimal weight vectors are computed to maximize the covariance cov(X m u m , Y m v m ). Let us gather in two matrices those weight vectors, U a p × M matrix and V a q × M matrix in which columns are respectively the vectors u m and v m . The rows correspond to all the weights related to a same variable. Then, we want 4. M to set to zero all the weights related to a same variable at the same time. This is why a L 2 penalization is introduced on the columns (U •,i ) of U .

Application of those penalization to the sparse PLS are presented and they can be generalized to sparse group PLS. Data have the structured presented in Figure 2.3 from Section 2.1. The fuse Lasso for a sparse PLS (penalizing X) is:

{U opt , U opt } = argmin U,V M m=1 Z (m) -U •,m V T •,m 2 F + λ 1 P Lasso (U ) + λ 2 P Fused (U ) with P Lasso (U ) = p i=1 U i,• 1 with P Fused (U ) = p-1 i=1 U i,• -U i+1,• 1 and Z (m) = X T M,• Y M,• (4.6)
where λ 1 and λ 2 are parameters driving the fuse Lasso penalization P . P Lasso is the classical Lasso penalization whereas P Fused is the fused term.

The common Lasso introduced by Obozinski is:

{U opt , U opt } = argmin U,V M m=1 Z (m) -U T m,• V m,• 2
F + λP (U ) with P (U ) = p i=1 U i,• 2 and Z (m) = X T M,• Y M,• (4.7)
where λ a parameter driving the penalization P .

The fuse Lasso forces the coe cients corresponding to di erent observation sets to be similar whereas the second penalization forces variables shrunk to zero to be shrunk to zero for all observation sets at the same time.

This second penalization can then be seen as a variable selection method. This is why we are interested in studying this kind of penalization.

Sparse group Partial Least Square for structured data

This part introduces the contributions developed during this PhD: the sgPLS for structured data. The theoretical background of the method has been presented in [START_REF] Broc | Penalized partial least square applied to structured data[END_REF]. The method has also been presented in two conferences at the French Society of Statistics [113] [114]. The content of those presentations is not included in this section as they do not add additionnal material to this dissertation and the presented journal articles. The contribution is a Partial Least Square to which the Lasso penalization presented by Obozinski [START_REF] Obozinski | Joint covariate selection and joint subspace selection for multiple classi cation problems[END_REF] is adapted.

In this part, we consider that data that are composed of independent observation sets and group of variables and this structure is knwon a priori. The presented methods allow 4.4. Sparse group Partial Least Square for structured data us to use the information about the edi cation of the data set in order to improve the performance of the analysis. Although this theory has been developed with the aim to answer a problem occurring in genomic public data sets, it can be applied to any eld where a certain observation set structure exists. The method called "penalized PLS for structured data" is de ned where separate PLS model are linked together with a common-Lasso penalization similar to the one developed in [START_REF] Obozinski | Joint covariate selection and joint subspace selection for multiple classi cation problems[END_REF]. Variables selected by the model are the same for all observation sets but the underlying model computes separated models for each observation set, giving both readability and exibility to the model. The theoretical background for this method is presented.

The implementation have been performed through an R package [START_REF] Liquet | Sparse Group Partial Least Square Methods[END_REF].

Framework of proposed method

The method is based on the following key results. In the same way that for Section 2.3.3 Lasso penalization is applied to X in the formulas but can be also applied to Y . Taking the notations from Section 4.3.2, we can introduce the following sparse group PLS formulation:

min U,V M m=1 Z (m) -U •,m V •,m T 2 F + P λ (U ) with P λ (U ) = λ K k=1 (p k ) U P k ,• F with Z (m) = X T Mm,• Y Mm,• (4.8)
The biconvex resolution use the following formulas. Fixing each V •,m 2 = 1 the optimal

U (k,m) is U (•,m) =     1 - λ 2 M m=1 Z (m) P k ,• V •,m 2
F     + Z (m) P k ,• V •,m (4.9)
and xing U 2 = 1 the optimal V is

V (•,m) , = Z (m)T V (•,m) with Z (m) = X T Mm,• Y Mm,• . (4.10)
where u m implies another thresholding function. The thresholding function sets to zero all the weights of a same gene at the same time. Gene selection is performed across all the study in this way.

Group Lasso and single variable Lasso can be combined in a Sparse Group Sparse PLS for 4. M structured data. The minimization problem becomes:

min U,V M m=1 Z (m) -U •,m V T •,m 2 F + (α) P (single) λ (U ) + (1 -α) P (group) λ (U ) with P (single) λ (U ) = λ p i=1 U i,• 2 P (group) λ (U ) = λ K k=1 (p k ) U P k ,• F with Z (m) = X T M ,• Y M ,• (4.11)
The biconvex resolution use the following formulas. Fixing each V (•,m) 2 = 1 the optimal U (k,m) is

U P k ,• =   1 - λ (1 -α) 2 ŨP k ,• 2
F    + ŨP k ,• with Ũi,• = 1 - λα 2 Z i,• V 2 + Z i,• V i,• (4.12)
and xing U = 1 the optimal V is

V •,m , = Z (m)T V •,m with Z (m) = X T Mm,• Y Mm,• (4.13)
The Joint-Lasso needs to set two penalization parameters. The larger the rst parameter is the lesser variables are selected in the model. The parameters can be optimized for instance by minimizing the error prediction of the model under a cross-validation procedure. The minimum zone being quite at we allow the parameters to be in a neighborhood of the minimum.

Remark: If the number of observations is di erent from one observation set to another, the observation sets with the largest number of observation can prevail in the model. In order to cope with this, we propose to divide data from a an observation set m by √ n m .

Conclusions

This chapter presents the problems speci c to meta-analysis. An extension of the PLS is proposed which has the merit of taking into account both the group of variables structure and the study structure. We will see in next chapter that the method can be used for a particular case of meta-analysis where few method are well established: pleiotropy studies. After this chapter, the original work where the method was proposed is reproduced.

Introduction

Since past years data analysis applied to high dimension in all domains has arisen [1]. Extracting information from ever larger data has become a trend in numerous fields and a large number of observation need to be gathered in order to evaluate statistical models. When data are hard to retrieve, gathering existing data sets is an efficient way for assembling data of high dimension. However this technique have its drawbacks : existing independent data sets can present intrinsic bias which can decrease the performance of the models used.

Those biases imply an unwanted underlying structure that will interfere with the signal we want to find. Bias can come from a difference in the source of information or the process used during the recollection of the data. This set structure has to be taken into account in order to improve the efficiency of the models. For instance, in genomics, data can be gathered from different studies because of the cost of the experimentation. Each clinical study may have been performed with its own chemistry protocol, with its own experimental material and on its specific populations, and bias can arise among the different data sets obtained. This "batch effect" is known and can significantly decrease the power of the analysis [2]. Another bias can occur in particular analysis where different "dynamics" exist between the studies : a predictor can be highly correlated with independent variable, but the direction of the correlation depends on the study. For instance, pleiotropy [3] is a field of genetics where a gene (predictor) can have a particular effect on different phenotypes (independent variables). Data can be gathered from different studies where the nature of the phenotype differs. Therefore, a gene can be highly correlated with each phenotype but an overall model struggles to catch the particularity of those effects.

In the article we tackle the problem of "batch effect" for dimension reduction such as Partial Least Square (PLS) method introduced by Wold [4]. Common dimension reduction techniques are Canonical Correlation Analysis (CCA) [5], Principal Component Analysis (PCA) [6] and PLS [7]. All these methods rely on the projection of the data into a subspace of lower dimension which represents most of the variation of the data. They are often posed as an eigen value problem [8]. PLS and CCA are both analysis two blocks of data and differ from the norm used whereas PCA analyses one block. Aiming to apply our method to supervised analysis, PLS approach is considered in this article.

In these dimension reduction techniques, results are formulated with new variables that are linear combination of the original ones. These combination can be hard to interpret due to the huge number of coefficient they represent. To answer this problem, Lasso methods have been used. Introducing this penalization shrink to zero the participation to the model of the least relevant variables. Results highlight a smaller number of variable that are easier to explain. In addition, noise of the signal is reduced and the power of the methods is boosted. These are called sparse method and have been developed for linear regression [9] [10], CCA [11], PCA [12] and Partial Least Square (sPLS). The sparse PLS (sPLS) has shown encouraging results [13] [14] and is the object of analysis in this article. The PLS and sPLS methods have also been used to control "the batch effect" when related studies are combining to increase sample size combining independent but related studies ( [15], [16]). In particular combining sPLS separating models and linking them can be an option like in the Multivariate INTegrative method (MINT) proposed in [16]. However, this approach cannot identify the true signal in presence of different dynamics.

For high dimensional regression problems, using problem-specific prior information improves the accuracy of the prediction and the interpretability of the model [17]. For example, in genomics, genes within the same pathway have similar functions and act together in regulating a biological system. Incorporation of this grouping structure is becoming increasingly common due to the success of gene set enrichment analysis approaches [18]. Using a model taking into account this variable group structure allow to improve the performance and the readability of the results. To this end sparse group Partial Least Square (sgPLS) have been developed [19] where two overlaid Lasso penalizations translate the group structure in the Partial Least Square formulation. A structure with group and sub-groups can also be handle by its generalization with three overlaid Lasso penalizations (sgsPLS [20]). Methods such as MINT don't take into account this kind of group structure.

In this article we consider data that are composed of independent observation sets. The observation sets are assumed to be known and are expected to introduce bias in the data. The presented methods allow us to use the information about the edification of the data set in order to improve the performance of the analysis. Although this theory have been developed with the aim to answer a problem occurring in genomic public data sets, it can be applied to any field where a certain observation set structure exists. Different method using Lasso penalization on data structured toward observation sets are discussed. In particular a "penalized PLS for structured data" is defined where separate PLS model are linked together with a common-Lasso penalization. In the end variables selected by the model are the same for all observation sets but the underlying model computes separated models for each observation set, giving both readability and flexibility to the model. We present the theoretical background for this method. Especially, we can show that the common-Lasso constraint that is used (i.e. a penalization across studies) can be be written as a standard Lasso with an overlaid group structure in an equivalent formulation of the PLS problems. We extend also this idea of common-Lasso constraint to a case where an a priori structure is known, where the variables are gathered into groups.

Notations

Before going into further details, the notation used in this article are introduced. Data are represented by X ∈ R p×n and Y ∈ R qn , two matrices, representing n observations of p predictors and q independent variables. Then X is a (n, p) matrix and Y a (n, q) matrix. For any matrix A of size (a, b) , for i ∈ {1, . . . , a} its rows are noted A (i,•) and for j ∈ {1, . . . , b} its columns are noted A (•,j) and for subsets ã ⊂ {1, . . . , a} and b ⊂ {1, . . . , b} resp. row and column sub-matrices are noted A (ã,•) and A (•, b) . For any vector ω of size a , for i ∈ {1, . . . , a} its elements are noted ω (i) and for subsets ã ⊂ {1, . . . , a} ω (ã) represents the elements of the vector corresponding to the positions in the subset. Matrices will always be in uppercase letters and vectors in lowercase letters to avoid any confusion.

The Frobenius norm on matrices is denoted F . We note X T the transpose matrix of X. The cardinal of a set S is noted #S. The positive value of a real number x is noted (x) + = |x|+x 2 .

Data with observation sets

Some data may present a structure among the observations gathered around groups of observations. For instance data can be composed of different studies, each one presenting its own mechanisms and bias. Let us consider M different sets in the data. Noting, for m ∈ N, M m a subset of {1, . . . , n}, let M = (M m ) m=1..M be a partition of {1, ..., n} corresponding to the observation sets. We note #M m = n m . Row blocks are defined by this partitions in Figure 1 (observations are assumed to be ordered by observation set).

Data with group of variables

Some data may present a structure among the variable gathered around groups. Let us consider that the variables are gather in K groups. Let P = (P k ) k=1..K be a partition of {1, ..., p} corresponding to this variable group structure. We note #P k = p k . We then have K k=1 p k = p. This partition can define column blocks among the variables if variables are assumed to be ordered by variable group. Both observation set structure and variable group structure can be defined at the same time like in Figure 2.

Formulation of the sparse Partial Least Square

In the literature, two formulation of the Partial Least Square exist, some extension of the PLS follow a first one usually called PLS1 [21] and other extensions follow a second one called "PLS2" [14]. In the context of the article we study exclusively the first one.

PLS and sPLS

Let X be a predictor matrix of size (n, p) and Y a matrix of independent variables of size (n, q). PLS finds successively couples of vector {u 1 , v 1 }, . . . , {u r , v r } for r < min(p, q) where the couples are composed of vectors of length resp. p and q, maximizing Cov(Xu i , Y v i ) for any i ∈ {1, . . . , r}, under the constraint that the family of vectors u 1 , . . . , u r and v 1 , . . . , . . . , v r are both of them orthogonal families [4]. It can be solved considering successive maximization problems [22], for h ∈ {1, . . . , r}

max ||u h ||2=||v h ||2=1 Cov(X h-1 u h , Y h-1 v h ), where X 0 = X, Y 0 = Y and X h-1 , Y h-1 are deflated matrices computed from u h-1 ,v h-1 , X h-2 , Y h-2 for h ∈ {2, .
. . , r}. The deflation depends on the PLS mode that is chosen ( [23], [4]). In this article we focus on the enhancement of the optimization problem and its Lasso formulation in its h-th step. According to [22] this step can be written as

{u opt , v opt } = argmin ||u||2=||v||2=1 X T Y -uv T 2 F + λP (u)
Lasso Penalty term for sparse PLS .

(

) 1
where the notation h is removed in order to simplify the formulation because we are interested in only one of the r steps of the PLS.

The sparse PLS introduces a penalization in this formulation of the problem. The penalty P (•) forces lowest values of u to be set to zero. The parameter controlling the degree of sparsity in the model is λ. In the presented formula the sparsity is applied only to the vector u, but a similar penalization can be define for v. In the context of this article we treat only the penalization of u but all the results stand also for a v penalization. The following sections compare different ways of writing the sPLS optimization problem presented in Equation ( 1) taking into account an observation or/and variable set structure.

Remark: Before analysis, the X and Y matrices are transformed by subtracting their column averages. Scaling each column by their mean and standard deviation is also often recommended [24]. Thus, the cross-product matrix X T Y is proportional to the empirical covariances between X-and Y -variables when the columns of X and Y are centered. When the columns are standardized, X T Y is proportional to the empirical correlations between X-and Y -variables. In this article the standardization is an important step to overcome the issue of the "batch effect" or to aggregate observations from different studies.

Formulation of the penalized PLS

Six different formulations of the sPLS are presented in this article. The first four correspond to data presenting an observation set structure like in Figure 1. The two last correspond to data presenting an observation set structure and a variable group structure like in Figure 2 which correspond to sgPLS models (see [19]). We can note that Problem 5 is a particular case of Figure 2 where there is only one observation set (M = 1). Loading vectors introduced in those figures refer to vectors formulated in the following problems. The study of Problems 4 and 6 are the main contribution of the article.

-Problem 1 (standard sPLS): This approach consists in simply considering all the observation set as one set. Data are standardized across all the sets, i.e. X and Y are standardized. The formulation is a standard sPLS problem

{u opt , v opt } = argmin ||u||2=||v||2=1 X T Y -uv T 2 F + λP (u). ( 2
)
In the model the loading u is composed of p elements and the loading v is composed of q elements. The sparsity of u is controlled by the parameter λ : for a given λ, s λ elements of u will be non-zero. -Problem 2 (MINT): Introduced in [16], this approach consists in considering M different sPLS problems corresponding to each of the M observation sets. Data are standardized within each observation set, i.e. for every m ∈ {1, . . . , M }, X (Mm,•) and Y (Mm,•) are standardized instead of X and Y . The sPLS problem is the same than in previous problem in Equation ( 2).

In the model, the loading u is composed of p elements and the loading v is composed of q elements. The sparsity of u is controlled by the parameter λ : for a given λ, s λ elements of u will be non-zero. -Problem 3 (multiple sPLS): This approach consists in considering all the observation set as one set, i.e. X (Mm,•) and Y (Mm,•) are standardized. Data are standardized within each observation set, i.e. for every m ∈ {1, . . . , M }, X (Mm,•) and Y (Mm,•) are standardized instead of X and Y . Formulation is a classic sPLS problem

{u m,opt , v m,opt } = argmin ||um||2=||vm||2=1 X (Mm,•)T Y (Mm,•) -u m v T m 2 F + λ m P (u m ). ( 3
)
In the model the set of loading {u m } m∈{1,...,M } is composed of p × m elements (p elements per u m ). The set of loading {u m } m∈{1,...,M } is composed of q × m elements (q elements per v m ). The sparsity of u m is controlled by the parameter λ m : for a given λ m , s m,λm elements of u m will be non-zero. Therefore, variables concerned by the shrinkage to zero will depend on the observation set m. -Problem 4 (" sparse PLS for structured data"): This approach consists in considering M different sPLS problems corresponding to each of the M observation sets. Data are standardized within each observation set, i.e. for every m ∈ {1, . . . , M }, X (Mm,•) and Y (Mm,•) are standardized instead of X and Y . All problems are solved at the same time with a common-Lasso. The formulation of the problem is

{U opt , V opt } = argmin U,V M m=1 Z m -U (•,m) V (•,m) T 2 F + λP (U ) with P (U ) = p i=1 U (i,•) 2 and Z m = X (Mm,•)T Y (Mm,•) . ( 4
)
In the model the set of loading U is composed of p × m elements (p elements per U (•,m) ). The set of loading V is composed of q × m elements (q elements per V (•,m) ).

The sparsity of all U (•,m) is controlled by the parameter λ : for a given λ, the same s λ elements of each U (•,m) will be non-zero. -Problem 5 (classical sgPLS): When variables can be gathered in groups (Figure 2), the sgPLS propose to add a group-Lasso penalization to the classical PLS. Data are standardized within each observation set, i.e. for every m ∈ {1, . . . , M }, X (Mm,•) and Y (Mm,•) are standardized instead of X and Y . The formulation of the problem is

{u opt , v opt } = argmin u,v Z -uv T 2 F + λ (1 -α) P group (u) + λαP variable (u) with P group (u) = K k=1 √ p k u (P k ) 2 , P variable (u) = p i=1 u (i) 2
and Z = X T Y.

(

) 5
In the model the loading vectors u and v is composed of resp. p and q elements. The penalization P variable forces single variables to be set to zero whereas the penalization P group forces sets of variables to be set to zero. The degree of sparsity in general in the model is λ whereas the parameter controlling the balance between both kind of sparsity is α. In this model elements of u corresponding to least relevant variables and least relevant group of variables are set to zero. -Problem 6 (" sgPLS for structured data"): In the same spirit of adapting problem 2 into problem 4, problem 5 can be adapted with a common-Lasso penalization. Data are standardized within each observation set, i.e. for every m ∈ {1, . . . , M }, X (Mm,•) and Y (Mm,•) are standardized instead of X and Y . The formulation of the problem is

{U opt , V opt } = argmin U,V Z m -U (•,m) V (•,m) T 2 F + λ (1 -α) P group (U ) + λαP variable (U ) with P group (U ) = K k=1 √ p k U (P k ,•) F , P variable (U ) = p i=1 U (i,•) 2 and Z m = X (Mm,•)T Y (Mm,•) . ( 6
) In the model the set of loading U is composed of p × m elements (p elements per U (•,m) ). The set of loading V is composed of q×m elements (q elements per V (•,m) ). In this model elements of U corresponding to least relevant variables and least relevant group of variables are set to zero. In this model the same variables and variable groups corresponding to least significant variables are set to zero for all U (•,m) , m ∈ {1, . . . , M }.

Solutions of the penalized PLS

The classical sPLS can be seen as a biconvex optimization problem. It can be solved by successively optimizing the loading u and v [22]. For a given v an optimized ũ is computed and the value of u is updated. Then the same is performed permuting the roles of u and v. This optimization process relies on solving the problems

u opt = argmin ||u||2=1 X T Y -uv T 2 F + λP (u) v opt = argmin ||v||2=1 X T Y -uv T 2 F . (7)
The solution of problems 1 to 3 (composed of standard sPLS methods) is given by the following theorem :

Theorem 1 The marginal optima in ũ and ṽ in the sPLS (Equation ( 1)) are : Fixing v, the optimal u opt for ( 7) is

u (i) opt = u (i) 0   1 - λ 2 u (i) 0 2   + , u 0 = X T Y v . ( 8
)
Fixing u, the optimal v opt for ( 7) is

v opt = Y T Xu. (9)
In this formula a soft thresholding sets down to zero loadings corresponding to variables whose scores are too low. Setting λ equal to zero we find the formulation of the PLS problem without Lasso constraint. A proof can be find in [13]. For problems 4, 5 and 6 the solution is more complex. Problem 4 introduces a common-Lasso penalization, problem 5 introduces a variable group structure and the problem 6 introduces both common-Lasso penalization and variable group structure. We can note that problem 4 is a particular case of problem 6 where there is no group penalty, i.e. α = 1. Problem 5 is a particular case of problem 6 where there is only one observation set, i.e. M = 1. The solution of problem 6 is given in theorem 2 (presented in the following) whereas solutions of problem 4 and 5 are corollaries of this theorem and can be found after the proof (corollary 1 and 2).

Theorem 2 The marginal optima in U and V in the "sparse group PLS for structured data" (Equation ( 6)) are :

Fixing V , the optimal U opt for ( 6) is :

U (P k ,•) opt = U (P k ,•) 1     1 - λ (1 -α) 2 i∈P k U (i,•) 1 2 2     + = U (P k ,•) 1   1 - λ (1 -α) 2 U (P k ,•) 1 F   + With U (i,•) 1 = U (i,•) 0   1 - λα 2 U (i,•) 0 2   + , U (•,m) 0 = Z m V (•,m) and Z m = X (Mm,•)T Y (Mm,•)
(10) Fixing U , the optimal V opt for ( 6) is :

V (•,m) opt = Z T m U (•,m) (11)
Proof. The proof is composed of three steps. In step 1 we settle the sub-gradient equation corresponding to the minimization problem. In step 2, we make the sPLS thresholding emerge in the equation. In step 3, we make emerge the group thresholding and prove the theorem. Let's settle the sub-gradient equation. The optimal U for a given V is min

U M m=1 Z m -U (•,m) V (•,m)T 2 F + λ (1 -α) K k=1 √ p k U (P k ,•) F + λ p i=1 U (i,•) 2 .
We note that the problem can be formulated making appearing the column blocks corresponding to the variable groups. A second formulation of the problem would be min

U K k=1 M m=1 Z (P k ,•) m -U (P k ,m) V (•,m) T 2 F + λ (1 -α) K k=1 √ p k U (P k ,•) F + λ K k∈1 i∈P k U (i,•) 2 .
We can see that the problem can be separated in K distinct problems for every k ∈ {1, . . . , K} min

U (P k ,•) M m=1 Z (P k ,•) m -U (P k ,m) V (•,m)) T 2 F + λ (1 -α) √ p k U (P k ,•) F + λ i∈P k U (i,•) 2 .
In order to solve this problem, let's consider the k-th problem developing the Frobenius norm min

U (P k ,•) M m=1 Trace Z (P k ,•) m Z (P k ,•)T m -2Trace Z (P k ,•) m V (•,m) U (P k ,m)T + Trace U (P k ,m) U (P k ,m)T +λ (1 -α) √ p k U (P k ,•) F + λα i∈P k U (i,•) 2 .
Taking the sub-gradient, the optimal U opt verify for m ∈ {1, . . . , M } -U

(P k ,m) opt + U (P k ,m) 0 = λ (1 -α) √ p k 2 Θ (P k ,m) g + λα 2 Θ (P k ,m) v with the (p × M ) matrix U 0 such that U (P k ,m) 0 = Z (P k ,•) m V (•,m) with the (p × M ) matrix Θ g such that Θ (P k ,m) g =    U (P k ,•) opt U (P k ,•) opt F if U (P k ,•) = 0 Θ g ∈ {Θ g , Θ g F ≤ 1} if U (P k ,•) opt = 0 and with the (p × m) matrix Θ v Θ (i,•) v =    U (i,•) opt U (i,•) opt 2 if U (i,•) opt = 0 Θ g ∈ {Θ v , Θ v 2 ≤ 1} if U (i,•) opt = 0 . ( 12
)
We can note that when there is no penalty (i.e. λ = 0,) U opt = U 0 is the solution of the non sparse problem.

The sub-gradient equation is settle (step 1). Let's now make emerge the thresholding of sPLS. We investigate in which case U (P k ,•) opt = 0, i.e. when loading corresponding to a group of variables is set to zero. If U (P k ,•) opt = 0 then U (i,•) opt = 0 for every i ∈ P k . Hence we have

U (P k ,m) 0 = λ (1 -α) √ p k 2 Θ (P k ,m) g + λαΘ (P k ,m) v 2 with Θ (P k ,•) g 2 2 ≤ 1
and with

Θ (i,•) v 2 2 ≤ 1. (13)
and for i ∈ P k we have also

U (i,•) 0 - λαΘ (i,•) v 2 = λ (1 -α) √ p k 2 Θ (i,•) g with Θ (i,•) g 2 2 ≤ 1
and with

Θ (i,•) v 2 2 ≤ 1.
Let's define

U (i,•) 1 =   1 - λα 2 U (i,•) 0 2   + U (i,•) 0 . (14)
We can establish in the following lemma, which makes emerge the variable thresholding term of sPLS like in (1) in Equation ( 12). Lemma 1

U (i,•) 1 2 ≤ U (i,•) 0 - λαΘ (i,•) v 2 2 ( 15
)
and there is a Θ v such that

U (i,•) 1 = U (i,•) 0 - λαΘ (i,•) v 2 . ( 16
)
Proof. if U (i,•) 0 2 ≤ λα 2 2 then 1 - λα 2 U (i,•) 0 2 + = 0 and then U (i,•) 0 = 0. The in-
equality is then true. Furthermore, there is a

Θ (i,•) v = U (i,•) 0
reach the equality ( 16). Otherwise,

U (i,•) 1 = 0 and U (i,•) 0 - λαΘ (i,•) v 2 = U (i,•) 0 - λαU (i,•) 0 2 U (i,•) 0 2 = U (i,•) 1
The inequality ( 15) is true because the equality ( 16) is reached. In any case the lemma 1 is proved.

From lemma 1 we can infer that in ( 12)

U (i,•) 1 2 ≤ λ (1 -α) √ p k 2 Θ (i,•) g 2
and the inequality can be reached as an equality.

We have

i∈P k Θ (i,•) g 2 2 = Θ (P k ,•) g 2 2
and we have also

i∈P k U (i,•) 1 2 2 ≤ Θ (P k ,•) g 2 2
and the inequality can be reached.

Therefore

Θ (P k ,•) g 2 2 ≤ 1 stand if and only if i∈P k U (i,•) 1 2 2 ≤ 1.
In the end we have

U (i,•) = 0 if i∈P k U (i,•) 1 2 2 ≤ 1.
Let's now consider that U

(P k ,•) opt = 0 or U (i,•) opt = 0 for at least one i ∈ P k then U (i,•) opt = U (i,•) 0 - λα 2 Θ (i,•) v - λ (1 -α) 2 U (i,•) opt U (P k ,•) opt F . If U (i,•) 0 2 ≤ λα 2 then we can set Θ (i,•) v such that U (i,•) 0 -λα 2 Θ (i,•) v = 0. Otherwise U (i,•) 0 -λα 2 Θ (i,•) v = U (i,•) 0 -λα 2
U (i,•) U (i,•) 2
.

In both cases we can consider that U

(i,•) 0 - λα 2 Θ (i,•) v = U (i,•) 1
.

From this point we find successively

U (i,•) opt = U (i,•) 1 - λ (1 -α) U (i,•) opt 2 U (P k ,•) opt F , U (i,•) opt   1 + λ (1 -α) 2 U (i,•) opt 2   = U (i,•) 1
and

U (i,•) opt   1 + λ (1 -α) 2 U (i,•) opt 2   = U (i,•) 1
.

Summing the square for every element of P k we have

i∈P k U (i,•) opt 2 2 1 + λ(1-α) 2 U (i,•) opt 2 2 = i∈P k U (i,•) 1 2 2 .
and hence

U (i,•) opt 2 F 1 + λ(1-α) 2 U (i,•) opt 2 2 = i∈P k U (i,•) 1 2 F = U (P k ,•) 1 2 F .
After extracting the value of U (P k ,•) F from this equation we finally find that

U (P k ,•) opt = U (P k ,•) 1   1 - λ (1 -α) 2 U (P k ,•) 1 F   .
Corollary 1 The solution to Equation ( 4) can be seen as a biconvex optimization problem. Fixing V , the optimal U opt for (4) is :

U (i,•) opt = U (i,•) 0   1 - λα 2 U (i,•) 0 2   + , U (•,m) 0 = Z m V (•,m) and Z m = X (Mm,•)T Y (Mm,•) (17)
Fixing U , the optimal V opt for (4) is :

V (•,m) opt = Z T m U (•,m) (18)
Corollary 2 The solution to Equation ( 5) can be seen as a biconvex optimization problem is : Fixing v, the optimal u opt for (5) is

u (P k ) = u (P k ) 1     1 - λ (1 -α) 2 i∈P k u (i) 1 2 2     + = u (P k ) 1   1 - λ (1 -α) 2 u (P k ) 1 F   + With u (i) 1 = u (i) 0   1 - λα 2 u (i) 0 2   + , u 0 = Zv and Z = XY T . (19)
Fixing u, the optimal v opt for ( 5) is

v opt = Z T u. (20) 5 Discussion
Problems 1 to 4 are discussed in this part, but we consider that the following remarks can be transposed to Problems 5 and 6, Problem 5 and 6 being resp. the equivalents of Problems 2 and 4 for data with a variable groups structure.

Size of the data

The larger data (in term of number of observations) are, the better models are supposed to perform. We can see that Problem 1 and 2 have the merit of performing an sPLS on data containing n observations whereas Problem 3 and 4 performs M different sPLS methods on data with resp. M m observations for m ∈ {1, . . . , M }. For some observation set the number of observations can be significantly smaller than the size of the hole data which can have a negative impact on the result.

Number of loading elements in the model

Number of loading elements is an important parameter to control. From one side, the bigger the number is, the more information can be stored by the model, from the other side having too much loading elements can give results harder to interpret and there is higher risk of over-fitting. Problems 1 and 2 have only p loadings for u whereas Problems 3 and 4 have M × p ones. For Problem 4 the number of loadings is important but the number of non-zero variables will vary between 1 and p in the same way than in Problem 1 and Problem 2. The Problem 4 gives readable results while keeping the flexiblity of a model with higher number of loading elements. For Problem 3 the non-zero variables can be different from one study to another, we cannot control weather a variable will be null for all studies and the number of non-zero variable will be significantly higher than for other problems.

Sensibility to batch effect

Batch effect can arise when data provided from different source present a bias. This effect can happen when observation sets are expected to introduce its intrinsic error. Therefore the cross-product matrices could be represented by a model like :

Z m = Z + E m for m ∈ {1, . . . , M }
Where Z follows a given law and E m are Gaussian noise with parameters depending on m. Under this hypothesis the standardization within studies can bypass this bias. Therefore Problem 3 to 4 can correct this kind of batch effect.

However, more complex bias can exist. For instance what happens if different observation sets have different dynamics? Let us consider a variable that is positively correlated in some observation sets and negatively correlated in others. In Problems 1 and 2 and their overall sPLS, the variable will have a small corresponding loading because positive and negative effect compensate each other and the variable will be cut because of the sparsity heuristic. In Problem 3, the distinction between all dynamics will be highlighted by the model. Finally, Problem 4 will select the same variable because it has a significant loading on every observation set. In the end, Problem 4 can handle more cases where the observation sets introduce bias.

Relation between Problem 4 and classic sgPLS

We can establish also that Problem 4 (sPLS method with a common-Lasso penalization) applied to matrices X and Y of size resp. (n, p)and (n, q) can be equivalent to a classical sgPLS without a standard Lasso on well chosen matrices X and Ỹ of size resp. (n, p × M ) and (n, q × M ). Those matrices are constructed by shifting the row blocks of X and Y : they are diagonal bloc matrices whose blocs are resp. X (Mm,•) and Y (Mm,•) for m ∈ {1, . . . , M )}. The corresponding loading vectors of size resp. p × M and q × M are called here resp. u e and v e . The representation of those objects is shown in Figure 3. On those basis the formulation of the sgPLS problem searching for optimal u e,opt and u e,opt would be :

In this formulation loading vectors u e and v e can be seen as the concatenation of the rows of resp. U and V in a unique uni-dimensional vector. This notation is interesting from a theoretical point of view because it ensure that Problem 4 can inherit properties from sPLS. However this notation is not wise for computational efficiency because the matrices X and Ỹ are M times bigger than X and Y , where M is the number of observation set. For implementation, computing directly the solution from equations ( 10) and (11) seems wiser.

Application on simulated data

Presented methods are illustrated on simulated data. A first simulation case presents data where a batch effect exists and a second simulation case presents data where different dynamics exist among different observation sets. For each case different noise levels are considered. Every simulation is performed 50 times. The code can be found at https://github.com/camilobroc/sgPLS_for_structured_data.

Design of the simulated data

In the following, a training data set of 900 observations gathered in 3 observation sets of 300 observations are generated and then a test data set of 300 observations gathered in 3 observation sets of 100 observations (for the training data M = 3 and n 1 = n 2 = n 3 = 300, for the test data n 1 = n 2 = n 3 = 100).

Batch effect cases

In first the simulation case, applications of the methods with data presenting a batch effect are performed. The simulation are performed with different noise levels. Data have an observation set structure and a group of variable structure as shown in figure 2. A batch effect implies that one same physical process is observed but the methods of measurement vary among the different group of observation. We represent this difference of measurements by a bias in depending on the observation set.

A matrix X with 1000 variables gathered in 50 groups of 20 variables (K = 50 and p 1 = • • • = p K = 20) and a matrix Y with 3 variables (q = 3) are generated. In order to mimic a batch effect, the generation procedure of the matrices resp. X and Y has different parameters depending on different sub-matrices of resp. X and Y . Those matrices are composed of a signal and a noise. The signal corresponds to a PLS model with one latent variable. For m ∈ {1, . . . , M }:

X (Mm,•) = H (Mm) C T Signal λ X,B m + µ X,B m 1 nm×p Batch + E (Mm,•) X Noise Y (Mm,•) = H (Mm) D T Signal λ Y,B m + µ Y,B m 1 nm×q Batch + E (Mm,•) Y Noise (22)
Signal The latent variable H is a (n × 1) column vector where each element follows a normal distribution of mean 0 and standard deviation 1. The loadings associated to this latent variable are resp. C a p × 1 column vector and D a q × 1 column vector corresponding resp. to X and Y .

Batch effect

The signal is blurred by a batch effect. The parameters λ

(X,B) m , µ (X,B) m , µ (Y,B) m and λ (Y,B) m
are real numbers depending on the observation set m. They control the shape of the batch effect. The notations 1 nm,p and 1 nm,q correspond to the matrices which elements are equal to 1 and of respective size n m × p and n m × q.

Noise

The noise is represented by E X , a (n × p) matrix, and E Y , a (n × q) matrix. The matrix E X is constructed by group of variables : for k ∈ {1, . . . , K}, the rows of

E (•,P k ) X follow a multivariate normal distribution N p k (O p k , λ X,E Σ p k ,ρ
) where ρ and λ X,E are real parameters and Σ p k ,ρ is a (p k × p k ) matrix which diagonal elements are equal to 1ρ and non-diagonal elements are equal to ρ. The notations p k stands fro the vector of size p k and which elements are all equal to 0. The rows of the matrix E Y follow a multivariate normal distribution N q (1 q × 0, λ Y,E Σ q,ρ ) where λ Y,E is a real parameter and Σ q,ρ is a (q × q) matrix which diagonal elements are equal 1ρ and non-diagonal elements are equal to ρ. The notations q to the vector of size q and which elements are all equal to 1. The parameter ρ represent a correlation between variables of a same group and λ Y,E and λ Y,E represent the noise levels.

The non null parameters of C are the 15 first variables of the first 4 group of variables. Among those elements resp. 15, 30 and 15 are equal to resp. 1, -1, 1.5 and the values are randomly distributed. Other parameters are given in table 1 and the noise levels are indicated in table 2.

Effects of different magnitudes among group of observations

This simulation case mimic data presenting different dynamics among observation sets. The generation process follow the same formulas as the previous one but with different parameters. The main difference with the previous cases is that the parameters λ X,B m for m ∈ {1, . . . , M } can have opposite signs. While in the first case a batch could be represented by a difference of magnitude, the effects can have here opposite directions. In this simulation case we are not interested in a bias concerning µ

(X,B) m or µ (Y,B) m
, the parameters are set to zero. The non null parameters of C are the 15 first variables of the first 4 group of variables. Among those elements resp. 15, 30 and 15 are equal to resp. 1, -1, 1.5 and the values are randomly distributed. Other parameters are given in table 1 and the noise levels are indicated in table 2.

X,B 1 2 0 µ X,B 2 -1 0 µ X,B 3 -1 0 λ X,B 1 1 1 λ X,B 2 0.8 -0.8 λ X,B 3 1.5 1.05 µ Y,B 1 2 0 µ Y,B 2 0 0 µ Y,B 3 -2 0 λ Y,B 1 0.6 0.6 λ Y,B 2 1.4 1.4 λ Y,B 3 1 1 D {1, -1, 1.5} {1, -1, 1.5}

Compared methods

In the first simulation case , methods corresponding to problems 1, 2 and 5 are compared whereas in the second simulation case, methods corresponding to problems 1,2, 4 and 6 are compared. For the methods corresponding to problems 1, 2 and 4 the penalization parameter λ is set such that the number of variables is equal to the true number of variables having an effect (in this case 60). For the method corresponding to problems 5 and 6, the penalization parameter λ is set such that the number of groups of variable is equal to the true one ( in this case 4) while α is chosen from the set {0.1, . . . , 0.9} by cross validation : the value giving the best Mean Square Error Prediction is kept.

results

The performances of the method are measured through the True Positive Rate (TPR), the Total Discordance (TD) and the Mean Square Error Prediction. The TPR is defined as :

T P R =
True Positives True Positive + False Negatives and TN is defined as :

T D = False Postives + False Negatives
Results of first and second simulation are given in table 2. In the first simulation case, data present a bias that depends on the observation set. Different noise levels are generated (λ X,E = 2, 20, 30). We can see that when noise is small (λ X,E = 2), MINT and "sparse group PLS for structured data" can retrieve the true variables whereas a classic sPLS cannot. We can also see that the MSEP is better for the "sparse group PLS for structured data" than for MINT. We can note that "sparse group PLS for structured data" misses a few true variables which gives a non null TD. This is due to the fact that the calibration of the method doesn't seek for a selection of the true number of variables, hence, a small number of true variables can be missed. When noise is greater (λ X,E = 20), a difference in terms of detection of true variables is observed. The classical PLS have a much worse TPR and TD while "sparse group PLS for structured data" is above MINT. When noise is even greater (λ X,E = 30), "sparse group PLS for structured data" clearly outperforms MINT.

X,E 1 = λ X,E 2 = 2 λ X,E 1 = λ X,E 2 = 20 λ X,E 1 = λ X,E 2 =
X,E 1 = λ X,E 2 = 2 λ X,E 1 = λ X,E 2 = 10 λ X,E 1 = λ X,E 2 =
In the second simulation case, data present a magnitude in the latent variables that depends on the observation set. Different noise levels are generated (λ X,E = 2, 10, 20). We can see that when noise is small (λ X,E = 2) we can see that in order to retrieve the true variables "s(g)PLS for structured data" is better than MINT which is better than sPLS. In the same way that in the first simulation case, "sgPLS for structured data" miss a few number of the true variables whereas "sPLS for structured data" do not because of the specificity of the calibration. We can see at noise level λ X,E = 10, that only the methods calibrated "for structure data" are able to retrieve the true variables. At the highest noise level (λ X,E = 20) we can see that sgPLS stands clearly above "sPLS for structured data", and those two methods outperform the existing ones.

Conclusion

In the end different ways of formulating an sPLS problem on data presenting an observation set structure have been discussed. The MINT formulation have the merit of being easy to implement and correct the batch effect. The novel method "sparse PLS for structured data" can also correct it. Further more it allows to take into account a lot of different bias, especially when the different observation set don't have the same dynamics. Despite its high number of parameters, the common-Lasso penalization ensures that the result is readable with a small number of selected variables in the overall analysis.

This article proved it's ability to inherit properties of sPLS. It's adaptation to variable groups developed in this article, called "sparse group PLS for structured data", is a notable example of " sparse PLS for structured data" benefiting from an extension of the sPLS. We can note also that it can be applied on either quantitative or qualitative variables as any sPLS can.

A simulation shows that the new methods can outperform existing methods for detecting a small signal in a large noise. Because its requirements on the nature of the data are very general we are confident that the method can be applied to the wide area of domains where sPLS is competitive.

Meta-analysis methods applied to pleiotropy

This chapter tackles the problem of pleiotropy in genetic epidemiology. At rst, an overview of the de nitions covered by the name pleiotropy are presented. Especially, its formulation as a meta-analysis method is explained. Then the application of the method introduced in Section 4.4 to pleiotropy studies is presented. The method was applied to a thyroid and breast cancer study and results have been submitted to BMC bioinformatics journal. The current version of this paper (draft) is included after this chapter.

Pleiotropy de ntions

Pleiotropy de nes cases in genetic epidemiology with complex relations between genome and phenotype, especially when one or more genetic traits have an e ect on several phenotype traits. It occurs in genetics when a genetic factor in uences di erent phenotypes. The de nition is quite large and encompasses di erent statistical approaches [15].

• It can be seen as a multivariate analysis where one genetic trait have an in uence on several phenotype traits.

• If phenotypes come from di erent sources, it can lead to a meta-analysis problem where several data sets have common genetic predictors but di erent phenotypes. The goal is to nd genetic markers related to several phenotype across di erent studies.

Pleiotropy enriches our insights into the mechanisms involved in the appearance of the phenotypes. Gene set analysis with pleiotropy is expected to improve the chances of revealing the underlying genetic architecture of complex phenotypes. Highlighting pleiotropy provides opportunities for understanding the shared genetic underpinnings among associated diseases.

The scope of application in genetic epidemiology is wide as it can cover a large number of di erent traits at the same time. For instance, 42 traits or diseases have been studied in a single analysis approach recently [START_REF] Joseph K Pickrell | Detection and interpretation of shared genetic in uences on 42 human traits[END_REF]. In this case, the goal is to search new relations between diseases. On the opposite side, other pleiotropy analysis can take into account a fewer number of disease that are already expected to have common biological mechanisms.

For instance, the large-scale Collaborative Oncological Gene-environment Study (COGS)

5. M studies breast, ovarian and prostate cancers and several cancer showed peliotropy e ects [START_REF] Sakoda | Turning of cogs moves forward ndings for hormonally mediated cancers[END_REF]. Common mechanisms in psychotic disorder disease has also been investigated [START_REF] Hong Lee | Genetic relationship between ve psychiatric disorders estimated from genome-wide snps[END_REF].

In this dissertation, we are interested in cases with phenotypes coming from di erent sources. The motivation comes from oncology where the onset of a disease is rare and nding individuals presenting di erent diseases (type of cancers) can be very hard. Hence, di erent data sets with di erent phenotypes have te be used.

For the rest of the chapter we are interested in pleiotropy tackled as a problem of metaanalysis where the goal is the identi cation of predictors common to two or more data sets.

Challenges speci c to pleiotropy

Meta-analysis is often necessary Genetic epidemiology aims at studying populations presenting given disease, which reduces the number of individuals available. For pleiotropy, the population of interest must present several illness traits which constraints even more the population sample and, thus, nding patients presenting targeted symptoms of disease may be impossible. In this case, an alternate is to use patients coming from di erent studies and each patient having only one disease.

An example of such issues exist in oncology which is the eld studying cancer. In oncology study, patients rarely present several cancer. In general, data related to di erent cancer are gather in distinct data bases that are usual analysed separately. To overcome this limitation, a meta-analysis is needed. This is the main motivation for studying pleiotropy as a meta-analysis in this dissertation.

Opposite e ects

Sometimes a variable can have a positive e ect for one phenotype and a negative one for another. In this case, some methods can have problems identifying an e ect because a mean procedure will compensate the positive and the negative contribution of the model. This e ect can be problematic especially in pleiotropy. Indeed, antagnostic pleiotropy has been observed in some studies, and hypothesis regarding their origins from an evolutionary point of view has been debated [START_REF] George | Pleiotropy, natural selection, and the evolution of senescence[END_REF]. Hence, statistical methods need to cope with this aspect.

Interest of sparse group PLS for structured data Regarding pleiotropy, many meta-analysis statistical tools can be used. So far, this method was used to investigate candidate genes across di erent cancer types and allowed to identify novel regions with pleiotropic e ects [START_REF] Rayjean | Cross cancer genomic investigation of in ammation pathway for ve common cancers: lung, ovary, prostate, breast, and colorectal cancer[END_REF]. Although a number of statistical approaches have been proposed in the literature for identifying the existence of pleiotropy at SNP-level, few work covers pleiotropy at gene-and pathway-levels. For instance, ASSET does not take into account a group structure of variables while SKAT only considers homogeneous e ect across studies which does not take into account e ect of opposite directions. The method proposed aims at lling this gap.

Challenges speci c to pleiotropy

Application to thyroid and breast cancer

In this contribution, the method presented in the previous chapter/paper is applied to the study of thyroid and breast cancer data. The results of the study are included in a paper that has been submitted to BMC Bioinfomatics.

The developed statistical approaches will be applied for enriching our insights about the genetic mechanisms of thyroid and breast cancer types. We are interested in exploring gene level as well as pathway level associations between the individual cancer types as well as at both cancer types together, as breast and thyroid cancers share common hormonal risk factors [START_REF] Sarah M Nielsen | The breast-thyroid cancer link: a systematic review and meta-analysis[END_REF]. Both cancers are hormonally mediated and occur more frequently in women than in men and have been found to be in uenced by reproductive and hormonal factors. Breast cancer is the most common form of cancer in women and thyroid cancer is a less common cancer and was associated with early age at menarche, late menopause, late age at rst full-term pregnancy, or low parity. Hormonal treatments, particularly those combining estrogens and synthetic progestagens can be a cause. Obesity is also a recognized risk factors for breast cancer in post-menopausal women [START_REF] Lof | Dietary fat intake and gestational weight gain in relation to estradiol and progesterone plasma levels during pregnancy: a longitudinal study in swedish women[END_REF]. On the other hand, thyroid cancer risk has been associated with late age at menarche and high parity. Conversely, breastfeeding and oral contraceptive use were associated with a lower risk of thyroid cancer [START_REF] Truong | Breast cancer risk, nightwork, and circadian clock gene polymorphisms[END_REF][START_REF] Xhaard | Menstrual and reproductive factors in the risk of di erentiated thyroid carcinoma in young women in france: a population-based case-control study[END_REF]. Anthropometric factors such as tall height and large body size have been also consistently associated with risk of thyroid cancer [START_REF] Fachal | From candidate gene studies to gwas and post-gwas analyses in breast cancer[END_REF]. Regarding breast cancer risk, established genetic variants currently explain ∼49% of heredity risk, with the 90 variants identi ed in GWAS accounting for about 14% [START_REF] Fachal | From candidate gene studies to gwas and post-gwas analyses in breast cancer[END_REF]. Thyroid cancer has been shown to be the only cancer for which the contribution of inherited genetic factors exceeds that of environmental factors [START_REF] Hemminki | Familial risk of cancer by site and histopathology[END_REF]. However, recent GWAS have identi ed a few susceptibility loci explaining about 4% of familial heredity [START_REF] Figlioli | Novel genetic variants in di erentiated thyroid cancer and assessment of the cumulative risk[END_REF].

Presentation of the data

The data used for this contribution and its expected structure is the one depicted in Chapter 1 and 2, that is, n observations are available with p variables for the genotype (SNP data) and q variables for the phenotype (traits), and the number of observations, n is inferior to p (see Figure 2.3). Data has been provided by the INSERM institute. SNP variables can be organized in groups corresponding either to genes or pathways. We suppose that one SNP is in exactly one gene and one pathway. Observations can also be gathered by groups corresponding to di erent studies. We use Beluhca data set (see Table 5.1) which includes data from CECILE, a french casecontrol study on breast cancer (1 125 cases, 1 172 controls) and from CATHY and Young-thyr, two french case-control studies on thyroid cancer (463 female cases and 482 female controls). All these individuals were genotyped using a customized microarray including 8716 genetic variants from 28 candidate pathways (648 genes) selected from KEGG database and from a litterature review. After quality controls, we retained 6 677 SNPs available for both type of cancers.

We went through the following test to pre-process the datasets:

• missing values were imputed using median among case/control and data centered to µ = 0. Then, we studied the correlation between SNPs to remove those with a coe cient equal to 1. This step is necessary because some methods don't allow correlated data. • Couple of extremely correlated SNPs belonging to same genes were eliminated.

• As some methods do not deal with overlapping groups, 10 non-overlapping pathways were selected and only the 3766 SNPs related to those groups were kept in the nal database. After all this steps, the new dataset is composed of 3766 SNPs, grouped in 337 genes and 10 non-overlapping pathways.

After, this pre-processing, data analyzed in this contribution are composed of 3766 SNP data corresponding to 337 genes that participate in 10 pathways.

Conclusion

As we have seen, pleiotropy is deals with meta-analysis in genetics and raises statistical challenges. From our knowledge, statistic methods for meta-analysis tackling both the opposite e ect possibility and the group of variable structure hasn't been quite investigated.

The proposed method joint-sgPLS aims at lling this gap. It leads to new perspectives for modeling pleiotropy in genetic epidemiology.

Penalized Partial Least Square for pleiotropy

Camilo Broc, Thérèse Truong and Benoit Liquet

Submitted to BMC Bioinformatics (2020)

Background

The increasing number of genome-wide association studies (GWAS) have revealed several loci that are associated to multiple distinct phenotypes, suggesting the wide existence of pleiotropic effects. Highlighting these cross-phenotype genetic associations could help to identify and understand common biological mechanisms underlying some diseases. Common approaches test the association between genetic variants and multiple traits at the SNP-level. In this paper, we propose a novel gene-and a pathway-level approach in the case where several independent GWAS on independent traits are available. The method is based on a generalization of the sparse group Partial Least Squares (sgPLS) to take into account variable groups, and a Lasso penalization that links all independent data sets. This method, called joint-sgPLS, is able to convincingly detect signal at the variable level and at the group level.

Results

Our method has the merit of proposing a global readable model while coping with the structure of data. It can outperform traditional methods and provides a wider insight in terms of a priori information. Results are also provided on synthetic data and on an application to real data. Genetic markers common to breast cancer and thyroid cancer are investigated in a real data application.

Conclusion

The joint-sgPLS shows interesting properties for detecting a signal. As an extension of the PLS, the method is suited for data with a large number of variables. The choice of Lasso penalization copes with structures of variables groups and observations sets. Furthermore, although the method has been applied to a genetic study, its formulation is adapted to any data with high number of variables and an exposed a priori structure in other application fields.

Keywords: Genetic epidemiology; High dimensional data; Lasso Penalization; Meta-analysis; Oncology; Partial Least Square; Pathway analysis; Pleiotropy; Sparse methods; Variable selection

Background

Genome-wide association studies (GWAS) have identified numerous genetic markers linked to multiple phenotypes, suggesting the existence of pleiotropy that occurs when a single variant or gene can influence several phenotype traits [1,2,3,4].

Highlighting pleiotropy provides opportunities for understanding the shared genetic underpinnings among associated diseases. As the genetic effects for most complex traits are small, several methods were proposed to combine results across studies of different phenotypes in order to improve the power of detecting cross-phenotype or pleiotropic associations.

We propose to focus on the case where several GWAS on different traits are available from different independent sources. One way to analyze pleiotropy is to analyze each trait separately in each independent datasets. In order to gain statistical power to detect pleiotropy, we propose a method that can combine the different traits into a meta-analysis taking into account that effect in opposite directions may exist across the different traits. We also propose a gene-based approach that combines the association signals from the single nucleotide polymorphisms (SNP) into a signal at the gene level or at the pathway-level. Genes will be defined by their coordinates and pathways are group of genes involved in a common biological mechanism. Independent data sets give an observation set structure to the data while the gene and the pathways give a group of variables structure. The challenge of pleiotropy is then to take advantage of this structure. In addition, possible biases between observation sets can be introduced, which is a common problem for meta-analyses, and methods for pleiotropy must take it into account [5]. Furthermore, such methods must cope with the case where a genetic variable have a positive effect on one disease and a negative effect on an other disease. Those opposite effects, may be hard to highlight statistically.

In this article, an extension of the Partial Least Square (PLS) method suited for meta-analysis for pleiotropy is proposed. It deals with observation sets and group of variables information while taking into account the possibility of opposite effects.

PLS is a dimension reduction developed by Wold [6] and that have been widely used for the analysis of data with large number of variables [7]. Unlike, its cousin method, the Principal Component Analysis (PCA) [8], the PLS deals with two blocks of data and is then widely used for genotype-phenotype analyses. Moreover its sparse extension using Lasso penalization have been successful for providing readable models [9]. Especially sparse group Partial Least Square can take into account group of variables as a priori information [10] [11]. For different group of studies an alternative Lasso penalization has been proposed by Obozinski [12] for a linear regression to deal with data made of different sets of observations. An adaptation of the Lasso penalization, the joint-sgPLS, has recently been proposed for the PLS [13], answering the specific structure of both groups of variables and sets of observations. In this article, we exploit the same idea to leverage pleiotropy effects, especially because the method copes with the challenge of detecting small possible opposite effects. The method is compared to two well established statistical methods in genetic studies. The first one, ASSET [14] extends standard fixed-effects meta-analysis to allow for potential opposite directions of the same SNP on different traits. However the method does not take into account the group level information such as pathways. The second one metaSKAT [15] permits to carry out gene-based meta-analysis, but all effects are in the same direction.

The developed statistical approaches will be applied to real dataset for enriching our insights about the genetic mechanisms of thyroid and breast cancer types. We are interested into exploring gene level as well as pathway level associations between the individual cancer types as well as at both cancer types together.

Notations

Data are represented by X ∈ R n×p and Y ∈ R n×q , two matrices, representing n observations of p predictors and q independent variables. The Frobenius norm on matrices is denoted F . We note X T the transpose matrix of X and the cardinal of a set S is noted #S. The positive value of a real number x is noted (x) + = |x|+x 2 and is equal to the number if the number is positive and equal to zero otherwise. In general, observation sets can represent the fact that different sets of observations come from different sources and must be analyzed accordingly. For instance, data coming from different studies may present biases. Variables groups can represent either a set of variables that is known to be quite correlated or a group of variables that must be treated together. For instance, in genetics a gene defines an established group of SNP variables and pathway define established group of genes. From one side, let us consider M different sets of observations in the data. Noting, for m ∈ N, M m a subset of {1, . . . , n}, let M = (M m ) m=1..M be a partition of {1, ..., n} corresponding to the observation sets. We note #M m = n m . Row blocks are defined by this partition. From the other side, let us consider that variables are gathered in K groups. Let P = (P k ) k=1..K be a partition of {1, ..., p} corresponding to this variable group structure. We note #P k = p k . We then we have

K k=1 p k = p.
Column blocks are defined by this partitions. Both observation set structure and variable group structure can be defined at the same time as shown in Figure 1. For matrices, the notation • is used to refer to blocks of matrices. For instance X •,P k is the block of matrix of X corresponding the columns of the k-th group of variables and X Mm,• is the block of matrix of X corresponding the columns of the m-th set of observations. Figure 1: Illustration of data structured by groups of variables and groups observations. Observations and variables are assumed to be ordered by resp. observations sets and variable groups. The notation p represents the number of variables of matrix X, q the number of variables of matrix Y , n is the number of observations. n 1 , • • • , n M are the resp. number of observations of each observation set. p 1 , • • • , p K are the resp. number of variables in each group of variables.

Sparse Partial Least Square for structured data In the literature, several formulations of the PLS exist [16]. In the scope of this article, only one of them, the called "PLS1" is considered [17].

PLS finds successively couples of vector {u 1 , v 1 }, . . . , {u r , v r } for r < min(p, q) where the couples are composed of vectors of length resp. p and q, maximizing Cov(Xu i , Y v i ) for any i ∈ {1, . . . , r}, under the constraint that the family of vectors u 1 , . . . , u r and v 1 , . . . , . . . , v r are both of them orthogonal families [6]. It can be solved considering successive maximization problems [18], for h ∈ {1, . . . , r}

max ||u h ||2=||v h ||2=1 Cov(X (h-1) u h , Y (h-1) v h ), (1)
where 2) for h ∈ {2, . . . , r}. The deflation depends on the PLS mode that is chosen ( [6,19]). In the following, the notation h is removed in order to simplify the formulation because we are interested in only one of the r steps of the PLS. The sparse PLS (sPLS) propose to add a penalization to the loading vectors u and v. The following equivalence is used:

X 0 = X, Y 0 = Y and X (h-1) , Y (h-1) are deflated matrices computed from u (h-1) ,v (h-1) , X (h-2) , Y (h-
argmax ||u||2=||v||2=1,u∈R p ,v∈R q Cov(Xu, Y v) = argmin ||u||2=||v||2=1,u∈R p ,v∈R q X T Y -uv T 2 F (2)
and the proof can be found in [10].

The sPLS [18] can be written as

{u opt , v opt } = argmin ||u||2=||v||2=1,u∈R p ,v∈R q X T Y -uv T 2 F + λP (u)
Lasso Penalty term for sparse PLS .

The sparse PLS introduces a penalization in this formulation of the problem. The penalty P (•) forces lowest values of u to be set to zero. The parameter controlling the degree of sparsity in the model is λ. In the presented formula the sparsity is applied only to the vector u, but a similar penalization can be define for v. In the context of this article we treat only the penalization of u but all the results stand also for a v penalization.

Remark : Before analysis, the X and Y matrices are transformed by subtracting their column averages. Scaling each column by their mean and standard deviation is also often recommended [20]. Thus, the cross-product matrix X T Y is proportional to the empirical covariances between X-and Y-variables when the columns of X and Y are centered. When the columns are standardized, X T Y is proportional to the empirical correlations between X-and Y-variables. In this article the standardization is an important step to overcome the issue of the "batch effect" or to aggregate observations from different studies. The subject has been discussed in [13].

Extensions of the sparse Partial Least Square

In the following, extension of the sPLS take into account an observation or/and variable set structure are presented. The last method has been recently developed [13] and deals with both kind of structures.

In order to cope with the group structure, sgPLS have been proposed [10]:

{u opt , v opt } = argmin ||u||2=||v||2=1,u∈R p ,v∈R q Z -uv T 2 F + λ (1 -α) P group (u)+ λαP variable (u) with P group (u) = K k=1 √ p k u P k 2 , P variable (u) = p i=1 u i 2 and Z = X T Y. (4)
where the loading vectors u and v is composed of resp. p and q elements. The penalization P variable forces single variables to be set to zero whereas the penalization P group forces sets of variables to be set to zero. The degree of sparsity in general in the model is λ whereas the parameter controlling the balance between both kind of sparsity is α. In this model elements of u corresponding to least relevant variables and least relevant group of variables are set to zero.

An extension using the joint Lasso penalization from Obozinski [12] has been proposed. Its formulation for the sgPLS is:

{U opt , V opt } = argmin U ∈R p×M and V ∈R q×M ||U•,m|| 2 =||V•,m|| 2 =1 for m∈{1,••• ,M } M m=1 Z (m) -U •,m V •,m T 2 F +λ (1 -α) P group (U ) + λαP variable (U ) with P group (U ) = K k=1 √ p k U P k ,• F , P variable (U ) = p i=1 U i,• 2 and Z (m) = X T Mm,• Y Mm,• , (5)
where the set of loading U is composed of p × m elements (p elements per U •,m ). The set of loading V is composed of q × m elements (q elements per V •,m ). In this model elements of U corresponding to least relevant variables and least relevant group of variables are set to zero. In this model the same variables and variable groups corresponding to least significant variables are set to zero for all U •,m , m ∈ {1, . . . , M }.

For all sparse methods, parameters driving the penalization (λ and α) must be chosen. In general, the choice is made through a K-fold cross validation, take an error in prediction as criteria. An example of the procedure can be found in the implementation of many extension of the sPLS [10] [21] [22].

Remark : The proposed joint penalization is biconvex, and thus multiple local minima may exist. In this work we implement an algorithm based on the alternating convex search [23]. It is known to be an effective method for biconvex optimisation and recent work has shown alternating methods for minimisation of biconvex functions converge with high probability to estimates with similar statistical accuracy to a global optimum [24].

Benchmark methods

Both ASSET and metaSKAT are considered as benchmark methods.

ASSET is a method suited for meta-analysis providing a p-value across studies [14]. The input of the method are single variables summary statistics which are combined by the method. ASSET exhaustively explores subsets of studies for the presence of true association signals that are in either the same direction or possibly opposite directions. However ASSET does not exploit the group information such as genes or pathways. Further, the current version of ASSET provides pleiotropy result for each variant which should be corrected using a FDR correction in order ro control possible false positive pleiotropy effect.

SKAT is a method to detect association between rare variants in a region and a phenotype (continuous or binary). It is a supervised test for joint effects of multiple variants in a region on a phenotype. The metaSKAT method can do the same but aggregating several studies. This method outputs a p-value corresponding to a set of variables, for instance a gene or a pathway. The method is based on a weighted sum of SKAT statistics of the different studies [15].

Results

The code used for running the methods is available on github (https://github. com/camilobroc/BMC_joint_sgPLS).

Simulation

Presented methods are illustrated on simulated data presenting the structure in Figure 1. From one side, SNP genotypes are coded as minor allele counting {0, 1, 2} and a certain correlation is expected within a group of SNP from the same linkage disequilibrium block. From the other side, phenotype data are binary and have a true effect from one or more genetic markers. In order to simulate the correlation between SNPs, for a group of variables P k , a multivariate normal distribution with n observations x

(continuous) k ∼ N p k (µ k , Σ k ) is
simulated where µ k is a null vector of size p k and Σ k is a p k × p k matrix with 1 on the diagonal and ρ k , coefficients controlling the correlation between SNPs within a group, outside of the diagonal. A simulation of this variable gives a matrix which represents simulated observations for group of variables k. Those blocks are concatenated in a n×p matrix, X continuous that represents the whole data.

In order to have {0, 1, 2} genotype data, a discretization is performed. For a given variable j ∈ P k , we aim at simulating a SNP variable with a Minor Allele Frequency (MAF), which we note MAF j . This MAF means that:

P (x j = 0) = (1 -MAF j ) 2 P (x j = 1) = 2MAF j (1 -MAF j ) P (x j = 2) = MAF 2 j .
To this aim, for a given MAF j , quantiles q (j)

1 and q (j)

2 are chosen such as P (x j ≤ q 1 ) = (1 -M AF j ) 2 and P (x j ≤ q

2 ) = (1 -M AF j ) 2 + 2M AF j (1 -M AF j ) A discrete genotype, X (discrete) , is computed such that X (discrete) i,j =      0 if X (discrete) i,j ≤ q (j) 1 1 if X (discrete) i,j ≤ q (j) 2 2 if X (discrete) i,j > q (j) 2 ,
where i ∈ {1, • • • , n} are simulated observations and j ∈ P k is a variable of the k-th group of variables.

For each observation i, a binary phenotype y i is simulated with a logit model

logit(π i ) = log( π i 1 -π i ) = α + p j=1 X (discrete) i,j β j ,
where π i = P (y i = 1|data), β j for j ∈ {1, • • • , p} is a regression parameter.

Then different simulations of the process can be performed successively in order to simulated several studies.

In this article, 8 combinations of parameters are considered. Data composed of two studies. The simulated genotype has 25 groups of 20 variables. The number of observations, n can be either 200 (cases 1,3,5,7) or 400 (cases 2,4,6,8). The intragroup correlation parameters ρ k are equal to 0.5 and the MAF is equal to 0.3 for each variable. The first 5 groups have an effect in the model of the simulation. Either all the SNPs of the group have an effect (cases 1,2,5,6) or half of the SNPs have an effect (cases 3,4,7,8). For each group, half of the non-null regression parameters are positives (taken at random) while the other half is negative. In cases 1,2,5 and 6, the absolute value of those parameters is set to exp(0.1) whereas in cases 3,4,7 and 8, the absolute value of those parameters is set to exp(0.5). Effect are either in the same direction from one study to another (cases 1,2,3,4) or in opposite direction (cases 5,6,7,8).

For all methods 50 replications of the data are performed. For the implementation of the sgPLS and joint-sgPLS, penalisation parameters must be chosen similarly to [10]. The penalization parameter λ and α are optimized through a K-fold penalization procedure with an error of prediction as criteria. Choosing a parameter λ is equivalent to set a number of selected groups [10]. In this simulation the grid of number of selected groups {1, • • • , 25} is used and the grid for α is {0.1, 0.5, 0.9}. Figures 2 and3 show the error of prediction performances through a cross-validation procedure of the sgPLS and joint-sgPLS in a simulation of case 1, for different levels for α and different levels of group selection. The observed mean and the variance of the error rate over 50 replication are presented. In the framework of the method the set of parameters corresponding to the lowest error of prediction rate is kept for the model.

For ASSET, sgPLS and joint-sgPLS, the variables selected by the models are compared to the variable having an effect on the true model. For metaSKAT, sgPLS We can see that at the SNP level sgPLS and joint-sgPLS have a higher TP than ASSET and have a better TP than metaSKAT at the group level (tables 1, 2, 3 and 4). The number of true positives is high for sgPLS, joint-sgPLS and ASSET however. We can remark that when effects in the true model are in opposite directions the joint-sgPLS stands out as its TP is clearly higher than other methods (tables 3 and 4). We can also see that the number of FP for sgPLS and joint-sgPLS are inflated when only 10 SNPs per true-groups are involved in the model of simulation (tables 2 and 4). This is probably due to the fact that the method selects too many variables with a selected group.

Overall, we can see that sgPLS and joint-sgPLS performs better for detecting effect in the same direction while joint-sgPLS is the method with the best performance for detecting opposite effects. MetaSKAT has a lower true positive rate than the two previous methods while having less total discordance which means it is more conservative on this simulation. A good compromise would be considering joint-sgPLS as both opposite and same direction effects are well selected.

Pleiotropy investigation on breast and thyroid cancer

The developed statistical approaches will be applied for enriching our insights about the genetic mechanisms of thyroid and breast cancer types. Thyroid and breast cancers share a lot of similarities in their biology: both are more frequent in women, are influenced by hormonal and reproductive factors and are hormonally-mediated. Moreover, individuals diagnosed with breast cancer are more likely to develop thyroid cancer as a secondary malignancy than patient diagnosed with other cancer types, and vice-versa [25]. Genetic factor contributing to the incidence of breast cancer have been extensively studied, and it is known that genetic variants explain approximately 49 percent of the familial risk to develop this disease. Using GWAS, 313 risk variants were identified for breast cancer [26]. On the other hand, GWAS studies on thyroid cancer have been scarce, due to the lesser incidence of this disease as well as the lack of data. However, we know that thyroid cancer is the only cancer for which genetic factors contribute more than environmental factors [27]. Only 4 loci have been associated with thyroid cancer risk and have been replicated in other studies [28]. One of them, 2q35, was also previously reported to increase risk of breast cancer [29]. To date, no study has been conducted to identify common genetic factors between breast and thyroid cancer. Exploring the genetic relationship between the two cancers would help to elucidate the common mechanisms between both disease and could permit to improve their diagnostic and therapeutic management.

We propose to illustrate the methods on real datasets, by investigating the pleiotropic effect of genetic variants from candidate pathways in breast and thyroid cancers.

Beluhca dataset includes data from CECILE, a french case-control study on breast cancer (1 125 cases, 1 172 controls) and from CATHY a french case-control study on thyroid cancer (463 female cases and 482 female controls). All these individuals were genotyped using a customized microarray including 8 716 genetic variants from 28 candidate pathways (648 genes) selected from KEGG database and from a literature review. After quality controls, we retained 6 677 SNPs available for both type of cancers. Missing values were imputed using the median among cases or controls and data were centered to µ = 0. When 2 SNPs were correlated at r 2 = 1, one of the SNP was removed and couple of extremely correlated (r 2 > 0.98) SNPs belonging to same genes were eliminated.

As dimension reduction methods such as sgPLS and joint-sgPLS need to be extended in case of overlapping groups of variables [30], 10 non-overlapping pathways were selected and only the 3766 SNPs related to those groups were kept in the final database. After all this steps, the new dataset is composed of 3766 SNPs, grouped in 337 non-overlapping genes and 10 non-overlapping pathways. The list of the pathways and there genes is displayed in Tables 6 and7 in Appendix A.

The methods implemented in this article are : ASSET, metaSKAT, sgPLS and joint-sgPLS. For metaSKAT, sgPLS and joint-PLS, SNP-level, gene-level and pathway-level results are given by the methods whereas ASSET gives only SNPlevel results. Hence, in the case of ASSET, genes corresponding to selected SNPs are considered. For each SNP i, an univariate logistic model for gene-disease association can be considered separately for thyroid data and breast data (thyroid and breast cancer, Figure 4).

As it has been presented before, for sgPLS and joint-sgPLS, a calibration of the parameters is generally performed through a cross-validation process. This process relies on the definition of a measure of performance: the error of prediction of the model. However, in genetic studies, the effects are small and the prediction performances based on genetic units are usually very low. The prediction performance results of sgPLS and joint-sgPLS shows no significant difference from different sets of penalization parameters (i.e., numbers of genes and pathways). In order to facilitate the interpretation, we present the results for calibration parameters set to 20 genes and 3 pathways and α = 0.5. We explore the stability of the methods using a bootstrap strategy.

The method sgPLS and joint-sgPLS are computed with those parameters on the data. A 100 bootstraps procedure is performed on the data. The methods sgPLS and joint-sgPLS are then implemented on each bootstrap. We are interested in knowing whether genes selected in the computation on original data are still selected in the bootstrap. Genes and pathways selected by the methods applied on the original data are considered as preselected. Then, the rates of selection of selected genes and pathways and non selected genes and pathways over the 100 bootstraps are compared. Figures 5 and6 present this rate for preselected and non preselected features. A gene and resp. a pathway is kept in the final selection if and only if it is preselected and its rate of selection among the bootstraps is higher than any other gene (resp. pathway) that is not preselected. We can see that for joint-sgPLS less genes are selected than for other methods (4 against resp. 20 and 18 for metaSKAT and sgPLS on both data).

Results of the selection are presented in Table 5 where the name of genes and pathways is presented. "sgPLS single" stand for the use of the sgPLS on thyroid and breast data separately while "sgPLS both" stands for the use of the method on a concatenation of both data set standardizing by study.

From one hand we can see that for non meta-analysis methods (SKAT and sgPLS), except for INSR, genes selected for thyroid cancer are distinct from genes selected for breast cancer. From the other hand, meta-analysis methods, which are metaSKAT, sgPLS, joint-sgPLS, select gene and pathways for both data sets.

We can note that LRRN6C is selected for thyroid data with the univariate model, but the gene is selected only by ASSET and not by the others. A similar remark can be said about RORA which is selected for breast cancer with the univariate model and ASSET and not with the other methods. The fact that those genes are selected only by SNP based methods would mean that only one SNP of the gene has an effect on the phenotype, but the rest of the gene do not have much effect.Longer genes has more chance to be selected by these SNP based methods.

The gene MAP2K2 is selected by all methods for thyroid and by joint-sgPLS as well. The gene seems to be associated with thyroid cancer, but the fact that it is selected also by sgPLS means it could be a candidate for pleiotropy. The same can be said for PLA2G6, MTHFD2, ERCC3 that is initally associated to breast cancer but could be pleiotropic.

We can see that joint-sgPLS selects the pathways of inflammatory response and obesity and obesity-related phenotypes that are not selected by other methods.These results are interesting as obesity and inflammation are known to play a role in both cancers. We can wonder if there is an association at the pathway-level that is not detected by other methods.

Remark : Results based on different choice of calibration parameters for sgPLS and joint-sgPLS (50, 100 genes and 5 pathways) showed similar patterns.

Discussion

In this article, the properties of the joint-sgPLS is presented and is compared to the classical sgPLS, the ASSET method and metaSKAT. From one hand, joint-sgPLS and ASSET are design to retrieve effect of opposite direction. From the other hand, joint-sgPLS, sgPLS and metaSKAT can give results at a gene-level instead of a SNPlevel whereas ASSET can not. We have seen through a simulation that joint-sgPLS can be competitive against any of the other methods considered. Hence, joint-sgPLS seems perfectly suited for meta-analysis where effects in opposite directions can exist which invite us to pursue further investigation with it in complex studies for genetic epidemiology such as pleiotropy.

Conclusion

We do believe that further investigation can be done on the same subject. In this article, sgPLS and joint-sgPLS have been applied with one component, but several components could be considered. This could lead to the selection of variables that are orthogonal to the selection of the first component but that have still a large participation to the covariance matrix. We acknowledge that on the application the stability of the method is an important point due to the fact that the cross-validation procedure is not satisfying for choosing the parameters of penalization. One improvement could consist in exploiting different the criteria of the procedure (the error prediction) with, for instance, stability measures [31]. Another improvement could consist in adaptating the adaptative Lasso [32] for our method which could bypass the stability questions.

In order to advance on the application, this study should be replicated on a larger data base. Particularly, thyroid cancer has been less studied than breast cancer, and data for thyroid are still scarce in this application. Other cases of pleiotropy could be investigated, for instance for the case where the phenotype is multivariate for each subject. The joint-sgPLS is suitable for any kind of phenotype, continuous or qualitative. R code is available from the author to reproduce the results and is available on github (https://github.com/camilobroc/BMC_joint_sgPLS).

Appendix A: Tables of genes and pathways

Conclusions and futur work

Genetic epidemiology aims at understanding better the in uence of the genome on the onset of a disease. It has promising perspectives it terms of personalized medicine and for pharmacology. The rise of high dimensional data in the eld, which are called genomics, has raised the stakes even further in our expectations in terms of ndings. The use of the term "big data" for this trend must be very cautious. Indeed, genomic data are massive but are structure in an atypical way compared to other data analysis domains because the number of variables is extremely high while the number of observations is moderate. Furthermore, the domain needs to cope carefully with the biological meaning of its models. The treatment in an appropriate way of groups of variables, groups of observations and covariates must translate medical knowledge.

Contributions

This dissertation deals with the development of statistical methods for genomic data where two genetic problems has been considered.

Both contributions answer a general framework speci c to genetic studies. The idea consists in nding algorithms with mathematical properties which are adapted to the eld of application. It leads to the development of methods suited for the study of "ill-posed" problems. It also leads to treat a group of variables structure as an a priori information. For this purpose, truncation, sparsity, weighting are di erent terms referring to distinct statistical method that follow a general philosophy: concentrate the models on sets of variables that have the highest contribution to the signal. However it leads generally to the choice of hyper-parameters to which results are heavily dependent. A proper calibration needs generally repeated measure and resampling techniques are a common answer to this. The present dissertation deals with all those aspects at the di erent levels of the contributions.

The rst one, gene-environment interaction, studies the in uence of a gene on the appearance of a disease. It aims at target speci c environmental factors for predicting risks in a population regarding their genome or at determining conditions for the application of a medical treatment. This work was published in the Journal of the French Society the implementation is available in the PIGE package for R (in CRAN).

The second one, pleiotropy, studies complex links when several genetic traits and several 5. M phenotypes are involved in a same biological mechanisms. The goal is, for instance, to highlight processes implying several diseases. In this thesis, the method has been applied to oncology data. This method has published in the Arabian Journal of Mathematics and its application to oncology will be submitted to BMC Bioinformatics. the implementation is available in the sgPLS package for R (in CRAN).

Discussion

Some aspects presented in this dissertation leave some open question and provide perspectives for further investigations.

From a statistical point of view, we do believe that the development of statistical methods that have been done in this thesis can be extended to an even larger number of frameworks.

• The rst contribution introduces the parametric bootstrap as a resampling method in the case a framework for interaction studies. The choice of resampling methods can be overlooked when developing a new method and putting an emphasis in the choice of the resampling methods, in the same spirit that in this contribution, is an interesting path of investigation.

• The second contribution introduces an extension of the sgPLS in order to take into account e ect in opposite directions in a meta-analysis. Especially, its use with a group-Lasso is proposed. From our knowledge tackling both opposite e ects and a group of variable structure in a same method hasn't been quite investigated. In the same spirit, other perpendicular statistical techniques for PLS could be integrated to the presented framework. For instance, an elastic net can deal with stability of the results and adaptive Lasso and its oracle properties can solve the choice of hyperparameters. Altough, this material hasn't been discussed in this thesis, it can be a promising development for future works.

From an epidemiologic point of view, the statistical methods could be applied to a larger number of biological problems.

• The rst contribution proposes to investigate G × E interaction while the second one propose to investigate pleiotropy. The next step would be to study G × E interaction with pleiotropy. The method could be applied to the data provided by the french health institute INSERM. From a statistical point of view, the method of the rst contribution could be adapted using a meta-analysis framework like presented in this dissertation. This could open new biological interpretations.

• Diseases studied in this dissertation are mainly cancers. However, pleiotropy can be related to other kind of diseases. The developed method could be applied to other elds than oncology.

• Other elds of research have data with features similar to genetics. Mainly, chemometrics or neuroscience handle data with a large number of variables compared to the number of observations and results at a group of variables level can be bene cial.

For instance in neuroscience, data are collected through probes that captures the 5.5. Discussion neurons signal in given areas of the brain overtime. Observations corresponding to related time frames and related brain areas can be correlated and can be used as an a priori structure of the variables.

List of Symbols, Notations and concepts

Biology and bio-statistics vocabulary D DNA Molecule in the kernel of human cells. The genetic information is the same in each cell of a human.

Gene

Identi ed section of the DNA. A gene can exist in di erent variations in a population and has an e ect on the functioning of a cell.

SNP

The single nucleotide polymorphism or SNP is the smallest measured genetic variations : single base-pair changes in the DNA sequence that occur with high frequency in the human genome. Several SNP can generally be found in one gene.

Variants, Alleles

The di erent versions of a SNP are called variants or alleles.

MAF

The Minor Allele Frequency is the proportion of the most rare variant of a SNP in the population. For instance, a SNP with a MAF of 0.30 implies that 30% of a population has the allele versus the more common allele (the major allele) which is found in 70% of the population. Rare variant Below 0.5% MAF, the variant is called rare. It can also be called a mutation. Low-frequency variant between 0.5% and 5% the allele is called low frequency. Common variant and above 5% the variant is called common variant.

Locus (pl. loci)

A position on a chromosome, either a gene or a SNP. Pathway A group of gene involved in a same biological process. Data A set of variables and observations. For each observation and each variable a value is collected. Data can be represented in a matrix where the rows represent the observations and the columns represent the variables. Predictors,explanatory variables Variables that are assumed to have an e ect on an outcome L S Outcome,dependent variable, response

The variable that we want to explain.

Covariate

Variables that gives an additional information and can be used.

Statistic

A statistic is a value computed with a method that characterize a data set. p-value A p-value is a value between 0 and 1 associated to a statistic. It indicates weather a statistic is abnormal under a certain hypothesis. The p-value has classically been at the center of discussions in statistics. It is considered by many statisticians as the corner stone of statistics, eventhough some opposite positions are arising. Biomarker Measurable indicator on an organisms. They are the variable of interest we look for in bio-statistics. Notably, Genes, SNPs, loci can be biomarkers.

1. 1 .Figure 1 . 2 -

 112 Figure 1.2 -Represenation of the DNA and its location in the cell. The cell (top left) contains the nucleus (bottom left) where chromosome are stored (middle).The chromosome is a long molecule with a double helix structure (right). Image taken from[1]

2

 2

Figure 2 . 1 -

 21 Figure 2.1 -Illustration of data structured by groups of observations. Observations are assumed to be ordered by observation set. p represents the number of variables of matrix X, q the number of variables of matrix Y and n is the number of observations. n 1 , • • • , n M are the resp. the number of observations of each observation set.

Figure 2 . 2 -

 22 Figure 2.2 -Illustration of data structured by groups of variables. Variables are assumed to be ordered by observation set. p represents the number of variables of matrix X, q the number of variables of matrix Y and n is the number of observations. p 1 , • • • , p K are the resp. the number of variables in each group of variables.

L

  be the p-values for each interaction test on the null hypothesis based on the b-th generated data set, 1 ≤ b ≤ B. The RTP statistic W (b)

  et al., 2007; Yu et al., 2009; Holmans et al., 2009; O'Dushlaine et al., 2009; Shahbaba et al., 2012; Carbonetto and Stephens, 2013; Evangelou et al., 2014a,b; Su et al., 2016). The challenges, properties and statistical methods for conducting pathway (and gene-set) analysis for GWAS data have been discussed and reviewed by Wang et al. (2011); Fridley and Biernacka (2011); de Leeuw et al. (2016).

L

  be the p-values for each interaction test on the null hypothesis based on the b-th generated dataset, 1 ≤ b ≤ B. The RTP statistic

3 .

 3 Generate B new datasets s * b, , b = 1, . . . , B by repeating B times the steps 1, 2 and 3.

Fig. 1

 1 Fig. 1 Illustration of data structured by group of observation. Observations are assumed to be ordered by observation set.

Fig. 2

 2 Fig. 2 Illustration of data structured by group of variables and group observation. Variables are assumed to be ordered by variable group.

2 F√ p k u e (P k ) 2 ,

 22 {u e,opt , v e,opt } = argmin ue,ve Zu e v e T + λ (1α) P group (u e ) + λαP variable (u e ) with P group (u e ) = K k=1

Fig. 3

 3 Fig.3Notation of X (grey rectangle) and Ỹ (red rectangle) to write the sPLS for structured data as a sgPLS.

Figure 2 :

 2 Figure 2: Mean and variance of the error of prediction in cross-validation of sgPLS, for one simulation of case 1 of the simulations. The cross-validation is performed for α ∈ {0.1, 0.5, 0.9} and for levels of group selection corresponding to {1, • • • , 25}.

Figure 3 :

 3 Figure 3: Mean and variance of the error of prediction in cross-validation of joint-sgPLS, for one simulation of case 1 of the simulations. The cross-validation is performed for α ∈ {0.1, 0.5, 0.9} and for levels of group selection corresponding to {1, • • • , 25}.

  (a) Percent of selection of genes for sgPLS on thyroid data, 100 bootstraps (b) Percent of selection of genes for sgPLS on breast data, 100 bootstraps (c) Percent of selection of genes for sgPLS on all data, 100 bootstraps (d) Percent of selection of genes for joint-sgPLS on all data, 100 bootstraps

Figure 5 :

 5 Figure 5: Percent of selection of genes for sgPLS and joints sgPLS on 100 bootstraps. (a) sgPLS on thyroid data (b) sgPLS on breast data (c) sgPLS on both data (d) joint-sgPLS. Genes selected on original data (preselected ones) are in blue while other genes (non preselected ones) are in red.

  (a) Percent of selection of pathways for sgPLS on thyroid data, 100 bootstraps (b) Percent of selection of pathways for sgPLS on breast data, 100 bootstraps (c) Percent of selection of pathways for sgPLS on all data, 100 bootstraps (d) Percent of selection of pathways for joint-sgPLS on all data, 100 bootstraps

Figure 6 :

 6 Figure 6: Percent of selection of pathways for sgPLS and joints sgPLS on 100 bootstraps. (a) sgPLS on thyroid data (b) sgPLS on breast data (c) sgPLS on both data (d) joint-sgPLS. Pathways selected on original data (preselected ones) are in blue while other pathways (non preselected ones) are in red. The pathways are noted: 1) Cell cycle 2) Circadian rhythm 3) Folate metabolism 4) Other glycan degradation 5) Obesity and obesity-related phenotypes 6) DNA repair 7) Metabolism of xenobiotics 9) Precocious or delayed puberty 10) Inflammatory response.

  Data is composed of 5 genes with 10 SNPs The true model is based on the main effect of E and the main effect of all SNPs from the second gene. The environment is marginally correlated with first gene but not with other genes. The outcome is therefore generated from the following model:logit[P(Y = 1|G 2 , E)] = α + β E E +Data is composed of 5 genes, 4 of them with 5 SNPs and one with 50 SNPs (the last one).True model is based on the main effect of E and 2 randomly selected SNPs from the first gene. The environment is marginally correlated with the first gene but not with other genes. The outcome is therefore generated from the following model:

	10 ∑ j=1	β SNP 2 j SNP 2 j
	-Case 3:
		1

j Journal de la Société Française de Statistique, Vol. 159 No. 2 56-83 http://www.sfds.asso.fr/journal © Société Française de Statistique et Société Mathématique de France (2018) ISSN: 2102-6238 -Case 2: λ 1 , . . . , λ 5 ∈ {0, 1} and 5 ∑ j=1

TABLE 1 .

 1 Generating parameters for cases 1 to 9. The notation used refer to part 4.1.

TABLE 2 .

 2 Simulation of case 1 with 1000 permutations and 1000 bootstrap resampling with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0	0.03	0.036	0.04	0.024	0.038	0.024
	Gene2	10	0.016	0.08	0.04	0.04	0.038	0.076	0.052
	Gene3	10	0.022	0.074	0.048	0.064	0.042	0.078	0.054
	Gene4	10	0.008	0.056	0.04	0.052	0.028	0.05	0.05
	Gene5	10	0.02	0.07	0.04	0.05	0.042	0.066	0.038
	FW ER BF	0.014	0.066	0.046	0.052	0.042	0.068	0.038
		FW ER	0.066	0.282	0.188	0.23	0.166	0.29	0.2
	Type-I Error: Pathway	0.014	0.066	0.048	0.052	0.028	0.078	0.052
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.002	0.016	0.008	0.01	0.01	0.02	0.03
	Gene2	10	0.028	0.062	0.052	0.062	0.044	0.064	0.044
	Gene3	10	0.028	0.046	0.028	0.044	0.028	0.05	0.038
	Gene4	10	0.018	0.034	0.03	0.044	0.028	0.038	0.03
	Gene5	10	0.028	0.054	0.036	0.042	0.042	0.054	0.066
	FW ER BF	0.028	0.052	0.038	0.042	0.042	0.064	0.042
		FW ER	0.098	0.188	0.148	0.188	0.142	0.208	0.186
	Type-I Error: Pathway	0.03	0.052	0.038	0.042	0.042	0.078	0.04
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.008	0.032	0.026	0.028	0.018	0.028	0.032
	Gene2	10	0.032	0.04	0.03	0.048	0.036	0.052	0.042
	Gene3	10	0.022	0.032	0.036	0.044	0.04	0.05	0.038
	Gene4	10	0.044	0.054	0.046	0.064	0.06	0.076	0.054
	Gene5	10	0.036	0.036	0.032	0.046	0.036	0.042	0.042
	FW ER BF	0.026	0.04	0.036	0.04	0.04	0.048	0.03
		FW ER	0.138	0.178	0.16	0.21	0.178	0.228	0.196
	Type-I Error: Pathway	0.026	0.04	0.036	0.04	0.042	0.048	0.044
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 3 .

 3 Simulation of case 2 with 1000 permutations and 1000 bootstrap resampling with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.018	0.09	0.032	0.046	0.03	0.094	0.046
	Gene2	10	0.002	0.028	0.054	0.044	0.036	0.032	0.04
	Gene3	10	0.008	0.062	0.03	0.046	0.02	0.054	0.038
	Gene4	10	0.012	0.084	0.042	0.056	0.04	0.104	0.062
	Gene5	10	0.016	0.076	0.03	0.038	0.032	0.07	0.038
	FW ER BF	0.008	0.084	0.03	0.034	0.04	0.076	0.038
		FW ER	0.056	0.304	0.172	0.212	0.152	0.31	0.204
	Type-I Error: Pathway	0.008	0.086	0.03	0.036	0.03	0.09	0.028
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.02	0.046	0.028	0.05	0.022	0.05	0.032
	Gene2	10	0.002	0.032	0.036	0.038	0.03	0.03	0.058
	Gene3	10	0.028	0.07	0.038	0.046	0.034	0.066	0.054
	Gene4	10	0.028	0.062	0.042	0.054	0.044	0.076	0.044
	Gene5	10	0.018	0.05	0.036	0.052	0.026	0.046	0.026
	FW ER BF	0.016	0.046	0.042	0.042	0.042	0.06	0.036
		FW ER	0.09	0.232	0.164	0.214	0.146	0.238	0.196
	Type-I Error: Pathway	0.016	0.046	0.048	0.044	0.032	0.066	0.036
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.056	0.07	0.044	0.064	0.056	0.078	0.046
	Gene2	10	0.014	0.03	0.022	0.03	0.024	0.026	0.04
	Gene3	10	0.038	0.054	0.03	0.042	0.042	0.06	0.038
	Gene4	10	0.032	0.042	0.028	0.036	0.036	0.046	0.062
	Gene5	10	0.042	0.048	0.028	0.038	0.038	0.046	0.038
	FW ER BF	0.032	0.048	0.032	0.038	0.034	0.07	0.038
		FW ER	0.164	0.216	0.144	0.194	0.178	0.226	0.204
	Type-I Error: Pathway	0.032	0.048	0.032	0.04	0.042	0.076	0.044
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 4 .

 4 Simulation of case 3 with 1000 permutations and 1000 bootstrap resampling with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	5	0.02	0.064	0.042	0.05	0.034	0.054	0.044
	Gene2	5	0.026	0.07	0.048	0.068	0.038	0.07	0.07
	Gene3	5	0.03	0.068	0.04	0.06	0.044	0.074	0.054
	Gene4	5	0.018	0.04	0.026	0.036	0.02	0.034	0.04
	Gene5	50	0.01	0.144	0.066	0.066	0.062	0.16	0.192
	FW ER BF	0.02	0.144	0.044	0.046	0.052	0.11	0.072
		FW ER	0.096	0.32	0.19	0.242	0.176	0.326	0.348
	Type-I Error: Pathway	0.022	0.146	0.046	0.048	0.058	0.132	0.008
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	5	0.026	0.048	0.036	0.05	0.032	0.048	0.058
	Gene2	5	0.028	0.042	0.038	0.052	0.036	0.056	0.042
	Gene3	5	0.038	0.056	0.048	0.056	0.044	0.068	0.06
	Gene4	5	0.042	0.062	0.032	0.062	0.04	0.06	0.052
	Gene5	50	0.032	0.098	0.032	0.038	0.06	0.11	0.088
	FW ER BF	0.038	0.08	0.042	0.05	0.044	0.092	0.088
		FW ER	0.156	0.278	0.172	0.232	0.198	0.304	0.278
	Type-I Error: Pathway	0.038	0.084	0.044	0.052	0.034	0.096	0.02
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	5	0.016	0.03	0.034	0.042	0.022	0.028	0.04
	Gene2	5	0.046	0.05	0.04	0.054	0.046	0.054	0.046
	Gene3	5	0.04	0.05	0.042	0.056	0.046	0.054	0.052
	Gene4	5	0.03	0.042	0.038	0.046	0.032	0.04	0.04
	Gene5	50	0.042	0.05	0.036	0.05	0.06	0.098	0.05
	FW ER BF	0.024	0.044	0.044	0.054	0.044	0.056	0.038
		FW ER	0.16	0.208	0.176	0.222	0.188	0.248	0.21
	Type-I Error: Pathway	0.024	0.044	0.046	0.056	0.04	0.06	0.012
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 5 .

 5 Simulation of case 4 for sample size of 200 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1		0.046	0.062	0.036	0.05	0.058	0.076	0.028
	Gene2		0.05	0.062	0.044	0.056	0.05	0.064	0.034
	Gene3		0.038	0.042	0.034	0.036	0.046	0.058	0.024
	Gene4		0.034	0.044	0.046	0.052	0.04	0.046	0.032
	Gene5		0.048	0.066	0.038	0.042	0.048	0.062	0.042
	Gene6		0.058	0.07	0.062	0.064	0.06	0.068	0.056
	Gene7		0.044	0.042	0.026	0.034	0.036	0.038	0.036
	Gene8		0.054	0.068	0.044	0.04	0.044	0.06	0.04
	Gene9		0.046	0.064	0.052	0.062	0.06	0.066	0.044
	Gene10		0.032	0.036	0.026	0.042	0.026	0.04	0.032
	Gene11		0.042	0.056	0.038	0.044	0.046	0.058	0.034
	Gene12		0.054	0.072	0.054	0.062	0.06	0.072	0.046
	Gene13		0.038	0.052	0.036	0.052	0.048	0.058	0.026
	Gene14		0.062	0.078	0.044	0.048	0.066	0.076	0.056
	Gene15		0.032	0.044	0.05	0.06	0.04	0.048	0.038
	Gene16		0.038	0.044	0.042	0.05	0.038	0.048	0.034
	Gene17		0.062	0.076	0.028	0.038	0.048	0.07	0.056
	Gene18		0.058	0.07	0.032	0.042	0.044	0.054	0.05
	Gene19		0.038	0.042	0.032	0.048	0.042	0.046	0.028
	Gene20		0.06	0.076	0.04	0.05	0.072	0.086	0.048
	FW ER BF	0.02	0.056	0.024	0.026	0.054	0.08	0.018
	FW ER	0.612	0.704	0.564	0.642	0.646	0.732	0.556
	Type-I Error: Pathway	0.028	0.07	0.028	0.04	0.084	0.152	0.004
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 6 .

 6 Simulation of case 4 for sample size of 500 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1		0.054	0.060	0.046	0.056	0.068	0.073	0.038
	Gene2		0.048	0.053	0.048	0.053	0.064	0.070	0.056
	Gene3		0.044	0.043	0.046	0.049	0.052	0.051	0.032
	Gene4		0.056	0.062	0.042	0.043	0.052	0.060	0.052
	Gene5		0.048	0.043	0.052	0.049	0.052	0.053	0.038
	Gene6		0.046	0.043	0.046	0.049	0.04	0.047	0.028
	Gene7		0.054	0.058	0.042	0.053	0.05	0.049	0.042
	Gene8		0.038	0.030	0.044	0.041	0.052	0.041	0.02
	Gene9		0.066	0.073	0.052	0.053	0.064	0.071	0.054
	Gene10		0.038	0.039	0.052	0.053	0.052	0.058	0.036
	Gene11		0.036	0.045	0.038	0.043	0.054	0.058	0.04
	Gene12		0.054	0.062	0.054	0.064	0.064	0.066	0.056
	Gene13		0.05	0.058	0.04	0.047	0.06	0.075	0.042
	Gene14		0.038	0.036	0.05	0.049	0.04	0.051	0.026
	Gene15		0.048	0.053	0.06	0.062	0.056	0.058	0.044
	Gene16		0.06	0.053	0.058	0.062	0.06	0.066	0.048
	Gene17		0.056	0.060	0.074	0.068	0.058	0.056	0.044
	Gene18		0.052	0.056	0.042	0.045	0.056	0.068	0.048
	Gene19		0.046	0.043	0.07	0.071	0.07	0.079	0.036
	Gene20		0.032	0.032	0.03	0.041	0.04	0.047	0.028
	FW ER BF	0.046	0.056	0.042	0.058	0.078	0.111	0.032
	FW ER	0.626	0.637	0.628	0.658	0.676	0.712	0.544
	Type-I Error: Pathway	0.056	0.068	0.056	0.071	0.122	0.156	0.02
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 7 .

 7 Simulation of case 4 for sample size of 1000 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1		0.038	0.038	0.05	0.052	0.052	0.054	0.04
	Gene2		0.066	0.074	0.058	0.054	0.076	0.07	0.054
	Gene3		0.052	0.054	0.042	0.044	0.046	0.056	0.04
	Gene4		0.044	0.054	0.038	0.038	0.042	0.044	0.03
	Gene5		0.058	0.06	0.038	0.038	0.044	0.052	0.046
	Gene6		0.052	0.058	0.066	0.068	0.06	0.06	0.036
	Gene7		0.058	0.06	0.058	0.054	0.064	0.066	0.046
	Gene8		0.048	0.052	0.044	0.044	0.046	0.046	0.034
	Gene9		0.036	0.036	0.038	0.038	0.032	0.036	0.028
	Gene10		0.062	0.06	0.064	0.064	0.07	0.072	0.052
	Gene11		0.052	0.054	0.042	0.042	0.046	0.056	0.046
	Gene12		0.044	0.048	0.032	0.036	0.038	0.03	0.034
	Gene13		0.046	0.044	0.062	0.07	0.052	0.058	0.036
	Gene14		0.05	0.048	0.054	0.046	0.05	0.058	0.036
	Gene15		0.06	0.054	0.054	0.056	0.062	0.06	0.052
	Gene16		0.034	0.036	0.052	0.048	0.054	0.05	0.028
	Gene17		0.034	0.032	0.048	0.05	0.042	0.044	0.028
	Gene18		0.056	0.06	0.038	0.032	0.052	0.048	0.048
	Gene19		0.048	0.048	0.052	0.046	0.044	0.046	0.038
	Gene20		0.044	0.048	0.06	0.066	0.056	0.06	0.036
	FW ER BF	0.044	0.056	0.044	0.05	0.076	0.092	0.046
	FW ER	0.614	0.61	0.63	0.622	0.634	0.67	0.536
	Type-I Error: Pathway	0.058	0.058	0.054	0.066	0.12	0.13	0.104
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 8 .

 8 Simulation of case 5 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.1	0.122	0.054	0.052	0.1	0.122	0.21
	Power: Pathway	0.062	0.096	0.044	0.042	0.084	0.098	0.012
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.512	0.536	0.162	0.17	0.496	0.526	0.532
	Power: Pathway	0.308	0.326	0.078	0.078	0.328	0.35	0.356
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene1	20	0.876	0.884	0.432	0.42	0.86	0.872	0.872
	Power: Pathway	0.722	0.744	0.184	0.196	0.722	0.726	0.836
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 9 .

 9 Simulation of case 6 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.076	0.09	0.046	0.052	0.094	0.086	0.186
	Power: Pathway	0.056	0.068	0.054	0.04	0.1	0.098	0.048
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.264	0.254	0.222	0.19	0.364	0.348	0.382
	Power: Pathway	0.162	0.152	0.11	0.094	0.226	0.192	0.212
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.62	0.616	0.73	0.7	0.798	0.782	0.752
	Power: Pathway	0.406	0.426	0.448	0.408	0.616	0.598	0.684
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 10 .

 10 Simulation of case 7 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.108	0.118	0.060	0.080	0.108	0.134	0.096
	Gene2	10	0.072	0.084	0.054	0.072	0.062	0.088	0.072
	Power: Pathway	0.060	0.062	0.040	0.034	0.086	0.128	0.004
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.272	0.278	0.144	0.148	0.270	0.272	0.268
	Gene2	10	0.158	0.170	0.100	0.100	0.138	0.144	0.132
	Power: Pathway	0.112	0.132	0.074	0.070	0.214	0.240	0.168
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.572	0.580	0.262	0.268	0.562	0.578	0.592
	Gene2	10	0.318	0.324	0.180	0.186	0.290	0.288	0.3
	Power: Pathway	0.298 FW ER BF stands for FWER results using Bonferroni correction 0.334 0.080 0.090 0.416	0.444	0.464

TABLE 11 .

 11 Simulation of case 8 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.102	0.114	0.108	0.110	0.134	0.153	0.124
	Gene2	10	0.092	0.094	0.104	0.102	0.128	0.125	0.120
	Power: Pathway	0.044	0.060	0.050	0.048	0.142	0.157	0.016
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.224	0.230	0.160	0.160	0.256	0.250	0.302
	Gene2	10	0.310	0.318	0.274	0.272	0.340	0.332	0.398
	Power: Pathway	0.108	0.116	0.096	0.100	0.262	0.280	0.168
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	20	0.376	0.372	0.436	0.438	0.476	0.467	0.550
	Gene2	10	0.538	0.548	0.544	0.544	0.602	0.618	0.700
	Power: Pathway	0.266 FW ER BF stands for FWER results using Bonferroni correction 0.284 0.292 0.288 0.516	0.526	0.568

TABLE 12 .

 12 Simulation of case 9 with 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP. 500 replications of the synthetic data are performed.

					sample size 200
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.070	0.080	0.040	0.045	0.070	0.065	0.090
	Gene2	10	0.095	0.105	0.070	0.075	0.095	0.095	0.080
	Power: Pathway1	0.055	0.065	0.04	0.030	0.110	0.130	0.012
	Gene3	10	0.110	0.125	0.095	0.090	0.100	0.120	0.120
	Gene4	10	0.085	0.100	0.085	0.100	0.110	0.115	0.118
	Power: Pathway2	0.065	0.095	0.04	0.050	0.105	0.105	0.016
	Power: Pathway all	0.055	0.070	0.05	0.025	0.115	0.155	0.012
					sample size 500
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.150	0.145	0.130	0.135	0.150	0.170	0.194
	Gene2	10	0.155	0.160	0.125	0.125	0.160	0.160	0.190
	Power: Pathway1	0.085	0.10	0.075	0.075	0.145	0.150	0.132
	Gene3	10	0.145	0.150	0.100	0.105	0.150	0.145	0.204
	Gene4	10	0.135	0.145	0.125	0.110	0.140	0.145	0.166
	Power: Pathway2	0.100	0.09	0.080	0.085	0.155	0.165	0.084
	Pathway all	0.065	0.10	0.075	0.070	0.260	0.230	0.148
					sample size 1000
	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size	Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap iSKAT
	Gene1	10	0.245	0.245	0.260	0.260	0.295	0.290	0.366
	Gene2	10	0.235	0.260	0.225	0.220	0.275	0.270	0.364
	Power: Pathway1	0.145	0.150	0.135	0.12	0.240	0.255	0.324
	Gene3	10	0.230	0.235	0.240	0.240	0.270	0.285	0.364
	Gene4	10	0.295	0.305	0.270	0.265	0.345	0.335	0.392
	Power: Pathway2	0.160	0.195	0.150	0.14	0.250	0.265	0.304
	Pathway all	0.160	0.220	0.170	0.18	0.425	0.465	0.484
			FW ER BF stands for FWER results using Bonferroni correction

TABLE 13 .

 13 Results of the investigation of gene-environment interaction of Circadian Pathway using 1000 permutations and 1000 bootstrap resamplings with FM, MinP and ARTP.

	Circadian Pathway	FM		MinP		ARTP		iSKAT
	Gene	Size Permutation Bootstrap Permutation Bootstrap Permutation Bootstrap
	ARNTL	24	0.04	0.001	0.1359	0.1578	0.0539	0.001	0.0078
	PER1	5	0.0849	0.035	0.032	0.04	0.049	0.0619	0.1132
	NPAS2	62	0.2647	0.1249	0.0879	0.1309	0.2957	0.1608	0.4338
	CSNK1E	9	0.3337	0.3427	0.4166	0.4965	0.4655	0.5265	0.6011
	CRY1	7	0.4935	0.5235	0.4985	0.5504	0.5295	0.5694	0.6081
	CRY2	9	0.5425	0.7722	0.6603	0.9111	0.6374	0.8761	0.3475
	PER2	11	0.9071	0.967	0.9401	0.983	0.8721	0.9491	0.734
	PER3	15	0.8901	0.984	0.8571	0.976	0.8941	0.995	0.5695
	CLOCK	11	0.9181	0.99	0.982	1	0.964	0.999	0.8104
	subpathway	0.2967	0.007	0.2537	0.2997	0.2957	0.0020	0.2095
	FBXL3	7	0.1658	0.0769	0.1239	0.1848	0.1628	0.1918	0.2986
	SKP1	4	0.4046	0.4236	0.2128	0.2318	0.2827	0.3057	0.4056
	CSNK1D	3	0.3706	0.3906	0.4256	0.4655	0.3736	0.4056	0.4412
	RBX1	2	0.4146	0.3876	0.5005	0.4825	0.4505	0.4296	0.4977
	BHLHE40	9	0.3277	0.3237	0.5864	0.7113	0.4126	0.4396	0.7204
	RORA	288	0.3836	0.3027	0.6184	0.7612	0.4515	0.4456	0.4576
	NR1D1	8	0.4476	0.4985	0.3906	0.4945	0.4206	0.4825	0.5439
	RORC	14	0.1958	0.1139	0.4226	0.6434	0.2687	0.5055	0.4815
	CUL1	23	0.4585	0.5814	0.1009	0.1578	0.3916	0.5105	0.1035
	TIMELESS	7	0.6643	0.7632	0.4406	0.5504	0.5395	0.6513	0.7471
	BTRC	13	0.977	0.998	0.4236	0.6064	0.7153	0.7013	0.8507
	RORB	34	0.5974	0.7163	0.5994	0.7972	0.5614	0.7293	0.7152
	FBXW11	8	0.8741	0.952	0.8711	0.9361	0.9121	0.962	0.7715
	BHLHE41	4	0.959	0.983	0.8891	0.957	0.9211	0.97	0.8227
	subpathway	0.9231	0.6474	0.7682	0.9101	0.8042	0.8951	0.5552
	circadien		0.6054	0.009	0.5085	0.6114	0.6374	0.02	0.4166

TABLE 14 .

 14 Running time: permutation and bootstrap performances using 1000 resampling related to the application. Results are in seconds.

		Running time
	permutation bootstrap iSKAT
	20187.5	20264.8	486.3

TABLE 15 .

 15 Running time: permutation and bootstrap performances using 1000 resampling related to the simulation case 5. Results are the average time in seconds over 500. replications.Based on the work of[START_REF] Yu | Pathway analysis by adaptive combination of p-values[END_REF], we have proposed an efficient practical tool for investigating gene-and pathway-environment interaction. Both FM and ARTP methods are extended in this context and available through our R package PIGE[START_REF] Liquet | PIGE: Self Contained Gene Set Analysis for Gene-And Pathway-Environment Interaction Analysis[END_REF]

					Running time
		size 200			size 500			size 1000
	Permutation Bootstrap iSKAT Permutation Bootstrap iSKAT Permutation Bootstrap iSKAT
	488	491	2.37	611	709	2.81	1170	1376	3.58

Table 1

 1 Table of the parameters used in first and second simulation cases.

		First simulation case Second simulation case
	ρ	0.05	0.05
	µ

Table 2

 2 Results for first and second simulation case. Results in terms of MSEP, TPR and TD are presented for each noise level.

		Simulation case 1
	Noise level	λ

Table 5 .

 5 1 -Data from Beluhca data set after pre-processing.

	5. M
	Cancer	SNPs Observations Missing values Pairs of SNPs removed
	Breast	6 677 2 297	34 869	42
	Thyroid 6 677 945	13 215	20

Table 1 :

 1 Perfomances in terms of mean number of TP, FP, FN and TN over 50 replications for methods sgPLS, joint-sgPLS, ASSET and metaSKAT in simulations cases 1 and 2 (when all SNPs of the 5 groups are involved in the true model and the effects are in the same directions from one study to another). The number of observations per study is n.

					n=200
	Variable level performances		Group level performances
		TP	FP	FN	TN		TP	FP	FN	TN
	sgPLS	46.78	29.24	53.22	370.76	sgPLS	3.54	3.82	1.46	16.18
	joint-sgPLS	40.84	27.16	59.16	372.84	joint-sgPLS	3.44	3.78	1.56	16.22
	ASSET	29.98	22.14	70.02	377.86	metaSKAT	2.22	1.08	2.78	18.92
					n=400
	Variable level performances		Group level performances
		TP	FP	FN	TN		TP	FP	FN	TN
	sgPLS	75.76	139.34	24.24	260.66	sgPLS	4.74	11.9	0.26	8.1
	joint-sgPLS	66.12	76.44	33.88	323.56	joint-sgPLS	4.48	7.62	0.52	12.38
	ASSET	47.74	25.4	52.26	374.6	metaSKAT	3.14	1.26	1.86	18.74

Table 2 :

 2 Perfomances in terms of mean number of TP, FP, FN and TN over 50 replications for methods sgPLS, joint-sgPLS, ASSET and metaSKAT in simulations cases 3 and 4 (when 10 SNPs of the 5 groups are involved in the true model and the effects are in the same directions from one study to another). The number of observations per study is n.

					n=200
	Variable level performances		Group level performances
		TP	FP	FN	TN		TP	FP	FN	TN
	sgPLS	17.96	36.28	82.04	363.72	sgPLS	1.46	4.32	3.54	15.68
	joint-sgPLS	43.44	25.2	56.56	374.8	joint-sgPLS	3.46	3.5	1.54	16.5
	ASSET	30.58	22.48	69.42	377.52	metaSKAT	2.2	0.98	2.8	19.02
					n=400
	Variable level performances		Group level performances
		TP	FP	FN	TN		TP	FP	FN	TN
	sgPLS	75.76	139.34	24.24	260.66	sgPLS	4.74	11.9	0.26	8.1
	joint-sgPLS	66.12	76.44	33.88	323.56	joint-sgPLS	4.48	7.62	0.52	12.38
	ASSET	47.74	25.4	52.26	374.6	metaSKAT	3.14	1.26	1.86	18.74

Table 3 :

 3 Perfomances in terms of mean number of TP, FP, FN and TN over 50 replications for methods sgPLS, joint-sgPLS, ASSET and metaSKAT in simulations cases 5 and 6(when all SNPs of the 5 groups are involved in the true model and the effects are in opposite directions from one study to another). The number of observations per study is n.

Table 4 :

 4 Perfomances in terms of mean number of TP, FP, FN and TN over 50 replications for methods sgPLS, joint-sgPLS, ASSET and metaSKAT in simulations cases 7 and 8 (when 10 SNPs of the 5 groups are involved in the true model and the effects are in the opposite directions from one study to another). The number of obervations per study is n.

Table 6 :

 6 First pathways and their corresponding genes. The number of SNPs for each gene is presented. The pathways are 1) Nicotinate and nicotinamide metabolism 2) Circadian rhythm 3) Folate metabolism 4) Other glycan degradation 5) Obesity and obesity-related phenotypes 6) DNA repair

	6)		6)		7)		7)
	gene	count	gene	count	gene	count	gene	count
	GTF2H5	5	POLL	3	UGT1A8	58	CYP1B1
	LIG3	5	RPA2	3	UGT2A1	29	UGT2B4
	MAPK8	5	SSBP1	3	MGST2	22	ADH4
	NTHL1	5	CCNH	2	CYP2C8	17	ADH5
	OGG1	5	ERCC4	2	AKR1C2	16	AKR1C1
	RAD50	5	GTF2H3	2	CYP2C9	15	ALDH3B1
	RPA4	5	HMGB1	2	CYP2B6	14	CYP2E1
	APEX1	4			EPHX1	14	NQO1
	CDKN1A	4			MGST3	14	ALDH1A3
	ERCC2	4			AKR1C4	13	CYP2F1
	FEN1	4			COMT	12	DHDH
	GTF2H4	4			CYP2C19	12	SOD2
	MBD4	4			CYP2S1	12	GSTA2
	POLE3	4			CYP2C18	11	GSTM3
	RAD51	4			ADH1B	10	GSTO2
	RFC2	4			ADH7	10	GSTP1
	UNG	4			GSTA4	10	CYP1A1
	CDKN2D	3			GSTZ1	10	CYP1A2
	CETN2	3			MGST1	10	CYP3A43
	CUL4A	3			NAT2	10	GSTA3
	ERCC1	3			ADH1C	9	GSTM4
	MLH1	3			AHR	9	UGT2B7
	MPG	3			NAT1	9	AKR1A1
	POLB	3			ADH6	8	ALDH3B2
	POLE4	3			AKR1C3	8	CYP2A6
	7)		8)		9)		10)
	gene	count	gene	count	gene	count	gene	count
	CYP3A4	3	TGFBR3	77	TGFB2	22	EGFR	103
	CYP3A7	3	EBF2	60	IL18RAP	13	CCND3	29
	GSTM2	3	BCAT1	51	CYP4F11	12	MAPT	19
	GSTM5	3	VDR	34	EPHX2	12	MAP2K4	15
	UGT2A3	3	KAL1	26	IL7	11	EGF	11
	UGT2B11	3	TM7SF3	19	IGFBP1	9	MAP2K2	11
	ADH1A	2	CASC1	17	IGFBP3	9	GSK3B
	CYP2D6	2	FGFR1OP2	12	IL17A	9	PLA2G6
			FGFR1	11	IL10	8	PTEN
			KRAS2	10	IGFBP4	7	MAP2K1
			CCR3	8	IL15	7	MYBL2
			KISS1	8	MMP25	7	AKT1
			PROK2	7	IL16	6	MAP2K3
			LIF	5	IL12A	5	FGFR3
			PROKR2	5	IL13	5	MAP2K7
			PTH1R	4	IL18	5	MAPK12
			NKX2-1	2	IL19	5	TP53I3
					IL2	5	CCNA2
					IL9	5	MAPK7
					PLA2G4B	5
					IL3	4
					TGFB1	4
					IL1B	3
					IL4	3
					IL1A	2
					IL23A	1

Table 7 :

 7 Last pathways and their corresponding genes. The number of SNPs for each gene is presented. The pathways are 6) DNA repair 7) Metabolism of xenobiotics 8) Precocious or delayed puberty 9) Inflammatory response 10) Cell cycle.

Journal de la Société Française de Statistique, Vol. 159 No. 2 56-83 http://www.sfds.asso.fr/journal © Société Française de Statistique et Société Mathématique de France (2018) ISSN: 2102-6238

Journal de la Société Française de Statistique, Vol. 159 No. 2 56-83 http://www.sfds.asso.fr/journal © Société Française de Statistique et Société Mathématique de France (2018) ISSN: 2102-6238 Gene-and Pathway-environment Interaction analysis

Acknowledgement

Most of the computations presented in this paper were performed thanks to the computing facilities the Direction du Num�rique of the Universit� de Pau et des Pays de l'Adour provided us.

Acknowledgements

The authors acknowledge Pascal Guénel for providing the breast and thyroid cancer data. The authors acknowledge also the calculus center MCIA (Mésocentre de Calcul Intensif Aquitain) for providing its facilities. The "Ligue contre le Cancer" is acknowledged as well for its support for Cross Cancer Genomic Investigation of Pleiotropy project.

Funding

This study was supported by the "Ligue contre le Cancer" for it's Cross Cancer Genomic Investigation of Pleiotropy project.

Declarations

Ethics approval and consent to participate Written informed consent for the present study was obtained from all participants. Study protocols were approved by the French ethic committees (CNIL, CCPPRB) (reference numbers 05-3144 for CATHY study and 04-53 for CECILE study).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions CB and BL designed the novel method. CB implemented the analysis. TT perfomed interpretation on real data. CB, TT and BL wrote the manuscript. All authors read and approved the final manuscript.

Author details

1 Laboratory for Data Sciences and Decision (Digiteo), LIST, CEA, , Gif-sur-Yvette, France. 2 Laboratoire De Mathématiques et de leurs Applications de PAU E2S UPPA, CNRS, , Pau, France. 3 Université Paris-Saclay, UVSQ, Inserm, CESP, 94807, , Villejuif, France. 4 Gustave Roussy, 94805, , Villejuif, France. 5 Centre of Excellence for Mathematical and Statistical Frontiers and School of Mathematical Sciences at Queensland University of Technology, , Brisbane, Australia.

Mathematical notations

Each notation refer to the example introduced above it.

R, N represent resp. the set of the real numbers and the set of the positive integers. Element of a set i, j, n, p, q ∈ N ∈ indicates that i,j,n, p, q are elements of a set. Subset S, S i , S j ⊂ R ⊂ indicates a subset of a set. Vector v ∈ R p this notation indicates that v is a p-vector of R, i.e., a nite ordered list of p elements of R.

Element of a vectors

is the vector made of the elements of a vector v ∈ R p with indices corresponding to a subset S ∈ R Matrix M ∈ R p×q M is a matrix of size p × q, i.e. with p rows and q columns. In general, in the rest of the document, Matrices are noted in capital letters. Element of a matrix M i,j M i,j is the element of M ∈ R p×q corresponding to the i-th row and j-th column. Submatrix M S i ,S j M S i ,S j is a matrix corresponding to the rows of M indicated by the subset S i ⊂ {1,

List of Figures

1.1 Illustration of a the structure of a cell including: the nucleolus which is inside the nucleus, the ribosomes, a mitochondrion, the cilia, the lysosome, the centrioles, the microtubules, the golgi smooth endoplasmic reticulum, rough endoplamsic reticulum (cell membrane), the cytoplasm, the chromatin. Image taken from [1]. 1.2 Represenation of the DNA and its location in the cell. The cell (top left) contains the nucleus (bottom left) where chromosome are stored (middle). The chromosome is a long molecule with a double helix structure (right). Image taken from [1] Observations and variables are assumed to be ordered by resp. observations sets and variable groups. p represents the number of variables of matrix X, q the number of variables of matrix Y and n is the number of observations. n 1 , • • • , n M are resp. the number of observations of each observation set, p 1 , • • • , p K , are the resp. the number of variables in each group of variables. .