Max Min Average

Max Min

a class representative model for pph-ud

A B S T R A C T

The interdisciplinary field of systems biology has evolved rapidly over the last few years. Different disciplines have contributed to the development of both its experimental and theoretical branches. Although computational biology has been an increasing activity in computer science for more than a two decades, it has been only in the past few years that optimization models have been increasingly developed and analyzed by researchers whose primary background is Operations Research (OR). This dissertation aims at contributing to the field of computational biology by applying mathematical programming to certain problems in molecular biology. Specifically, we address three I am very grateful to Professor Martine Labbé, who introduced me to the world of research during my master thesis and who encouraged me to pursue a Ph.D. degree in Brussels. I also want to thank the remaining members of my thesis committee for the time and effort they devoted to my work. In particular, I want to convey my gratitude to Professor Tom Lenaerts, who has been part of my comité d'accompagnement at ULB from the beginning of this thesis. I always appreciated his sincerity and moral support during all these years. At this point, I would like to thank Professor Bjarni V. Haldòrsson, for inviting me to visit his group for a research internship. During this stay I enjoyed immensely his collaboration, his motivation, focus on his research and openness to collaborative work are to this day a great inspiration for me. xiii I would be remiss if I did not mention Professor Daniele Catanzaro, who introduced me to the topic of computational biology in the very beginning. But mostly, I would like to thank him for his support and friendship that have been invaluable to me and for being the person I've always looked up to.

Many thanks go out to all member of the Department of Computer Science that has offered me a contract as Teaching Assistant throughout my PhD. This experience made me passionate about teaching and taught me how to be a better teacher every day.

Then I would like my friends and colleagues both from G.O.M. and the Computer Science Department (former and current) for always being there to discuss science but also, and probably most importantly, for a coffee break or a small chat to recharge batteries. They provided a very pleasant working environment, which I believe played an important part in my achievements, and some of them are now very good friends.

And finally, I wish to convey special thanks to my family and my friends for helping and supporting me in so many ways throughout these years.

Finding a maximum node-weighted clique in a

Max-Point Tolerance Graph (MPTG) . . . . . . . . .

3.2.3

Initializing the columns of the RMP . . . . . . . An example of a trio, i. e., a set of two parents and their offspring, and their genotypes. SNPs are highlighted in red. Specifically, the first highlighted column represents a SNP inconsistent with a loss of heterozygosity; the second highlighted column represents a SNP consistent with a loss of heterozygosity; the third highlighted column represents a SNP

showing an evidence of a loss of heterozygosity. .

Figure

The ALU Graph (ALUG) is based on the reads showed in Figure 12b. In particular, Figure 13 shows the mapping of improper reads to a segment of the reference genome, while Figure 13b shows Performances of the CRM for PPH-UD when considering input data having an error ratio of 10%.

Table 3 Performances of the CRM for PPH-UD when considering input data having an error ratio of 15%. et al. [2004], Bła żewicz et al. [2005] and [START_REF] Banga | Optimization in computational systems biology[END_REF].

Optimization and optimality are certainly not new concepts in biology. The structures, movements, and behaviours of animals, and their life histories, have been shaped by the optimizing processes of evolution or of learning by trial and error [START_REF] Alexander | Optima for Animals[END_REF][START_REF] Sutherland | The best solution[END_REF]. Moreover, optimization theory not only explains current adaptations of biological systems, but it also helps to predict new designs that may yet evolve [START_REF] Alexander | Optima for Animals[END_REF][START_REF] Sutherland | The best solution[END_REF]. Some reviews on the use of optimization in the close fields of computational biology and bioinformatics can be found in [START_REF] Greenberg | Opportunities for Combinatorial Optimization in Computational Biology[END_REF] and [START_REF] Larrañaga | Machine learning in bioinformatics[END_REF].

In the following section, we focus our attention on a specific methodology used in systems biology for the studies of common diseases and complex traits in human beings.

genome-wide association studies

The combinations of genetic and environmental factors are the causes of common diseases such as cancer, obesity, diabetes, cardiovascular, and inflammatory diseases [START_REF] Consortium | The international hapmap project[END_REF]]. Discovering these genetic factors would provide fundamental new insights into the diagnosis and treatment of human diseases. The recent completion of the sequencing phase of the Human Genome Project [START_REF] Frazer | A second generation human haplotype map of over 3.1 million SNPs[END_REF] has shown that any two copies of the human genome differ from one another by approximately 0.1% of nucleotide sites i.e., one variant or polymorphism per 1000 nucleotides on average [START_REF] Li | Low nucleotide diversity in man[END_REF][START_REF] Wang | Largescale identification, mapping, and genotyping of single-nucleotide polymorphisms in the human genome[END_REF][START_REF] Cargill | Characterization of single-nucleotide polymorphisms in coding regions of human genes[END_REF][START_REF] Halushka | Patterns of single-nucleotide polymorphisms in candidate genes for bloodpressure homeostasis[END_REF]. This result has suggested as a possible approach to identifying genetic risk factors the search for an association between the variant sites of a specific chromosome region and a disease [START_REF] Risch | The future of genetic studies of complex human diseases[END_REF].

A number of association studies which focused on candidate genes, have already led to the discovery of genetic risk factors for several diseases. Examples include type 1 diabetes [START_REF] Risch | The future of genetic studies of complex human diseases[END_REF][START_REF] Dorman | Worldwide differences in the incidence of type I diabetes are associated with amino acid variation at position 57 of the HLA-DQ beta chain[END_REF][START_REF] Nistico | The CTLA-4 gene region of chromosome 2q33 is linked to, and associated with, type 1 diabetes. Belgian Diabetes Registry[END_REF], type 2 diabetes [START_REF] Deeb | A Pro12Ala substitution in PPAR γ 2 associated with decreased receptor activity, lower body mass index and improved insulin sensitivity[END_REF][START_REF] Altshuler | The common PPARgamma Pro12Ala polymorphism is associated with decreased risk of type 2 diabetes[END_REF]], Alzheimer's disease [START_REF] Strittmatter | Apolipoprotein E and Alzheimer's disease[END_REF], deep vein thrombosis [START_REF] Dahlbäck | Resistance to activated protein C caused by the factor V R506Q mutation is a common risk factor for venous thrombosis[END_REF]], inflammatory bowel disease [START_REF] Hugot | Association of NOD2 leucine-rich repeat variants with susceptibility to Crohn's disease[END_REF][START_REF] Ogura | A frameshift mutation in NOD2 associated with susceptibility to Crohn's disease[END_REF][START_REF] Rioux | Genetic variation in the 5q31 cytokine gene cluster confers susceptibility to Crohn disease[END_REF], hypertriglyceridemia [START_REF] Pennacchio | An apolipoprotein influencing triglycerides in humans and mice revealed by comparative sequencing[END_REF], schizophrenia [START_REF] Stefansson | Neuregulin 1 and susceptibility to schizophrenia[END_REF], asthma [START_REF] Van Eerdewegh | Association of the ADAM33 gene with asthma and bronchial hyperresponsiveness[END_REF], stroke [START_REF] Gretarsdottir | The gene encoding phosphodiesterase 4D confers risk of ischemic stroke[END_REF], and myocardial infarction [START_REF] Ozaki | Functional SNPs in the lymphotoxin-alpha gene that are associated with susceptibility to myocardial infarction[END_REF].

One approach to perform association studies consists of testing each putative variant site for correlation with the disease (the direct or brute-force approach [START_REF] Consortium | The international hapmap project[END_REF]. At present, this approach is limited to sequencing the functional parts of suspected genes, selected on the basis of a previous functional or genetic hypothesis [START_REF] Consortium | The international hapmap project[END_REF]]. Unfortunately, the direct approach entails the sequencing of numerous patient samples to identify the responsible variant sites and hence it is prohibitively expensive.

An alternative approach, called the indirect approach [START_REF] Consortium | The international hapmap project[END_REF]], consists of exploiting human sequence variation as genetic markers. In fact, over 90% of sequence variation among individuals is due to common variant sites [START_REF] Li | Low nucleotide diversity in man[END_REF], most of which arose from single historical mutation events on the ancestral chromosome [Consortium 2005]. Hence, in a group of people affected by a disease, the variant sites causing the disease will be enriched in frequency compared with its frequency in a group of unaffected ones. This observation proved of considerable value, for example, in the identification of the genes responsible for cystic fibrosis and diastrophic dysplasia [START_REF] Consortium | The international hapmap project[END_REF][START_REF] Frazer | A second generation human haplotype map of over 3.1 million SNPs[END_REF]].

The indirect approach is generally preferred to the direct one because it requires neither sequencing multiple patient samples nor prior knowledge of putative functional variant sites. However, in order to be applicable, the indirect approach requires determination of the common patterns of DNA sequence variation in the human genome , by characterizing sequence variants, their frequencies, and correlations between them [START_REF] Consortium | The international hapmap project[END_REF]]. In general, this is not an easy task because the current molecular sequencing methods only provide information about the combination (or conflation) of the paternal and the maternal of an individual (also called genotype) [Halldórs-son et al. 2003].

Genome-Wide Association Studies (GWAS) is one of the most promising approach in the literature [START_REF] Stranger | Progress and promise of genome-wide association studies for human complex trait genetics[END_REF], and has been able

to successfully identify a very large number of polymorphism associated with disease (e.g. [START_REF] Styrkarsdottir | Multiple genetic loci for bone mineral density and fractures[END_REF], Gudbjartsson et al. [2008], and [START_REF] Rivadeneira | Twenty loci associated with bone mineral density identified by large-scale metaanalysis of genome-wide association datasets[END_REF]). In this dissertation, our purpose is to illustrate the capabilities, opportunities and benefits that mathematical optimization can bring to research in systems biology.

First, we introduce several basic concepts to help the readers unfamiliar with mathematical optimization.

combinatorial optimization

The integer programming books by [START_REF] Wolsey | Integer Programming[END_REF] and by [START_REF] Conforti | Integer Programming[END_REF] are the inspirational references on which the following section is based.

The key elements of mathematical optimization problems are the The set of feasible solutions to (P),

S := {(x, y) ∈ Z n × R p + : Ax + Gy b}
is called a mixed-integer linear set. While a mixed 0,1 linear set is set of the form

S := {(x, y) ∈ {0, 1} n × R p + : Ax + Gy b}
in which the integer variables are restricted to take value 0 or 1. A mixed 0,1 linear program is an mixed integer program whose set feasible solution is a mixed linear 0,1 set.

Solving a MILP is a hard problem, in general. A commonly used solution approach consists of (i) solving a relaxation (that is easier to solve numerically) and (ii) of giving a good approximation to the optimal objective value.

Given a mixed integer set S ⊆ Z n × R p , a linear relaxation of S is a set of the form

P := {(x, y) ∈ R n × R p : A x + G y b } that contains S.
A Linear Programming (LP) relaxation of (P) is a linear program max x,y {c T x + h T y : (x, y) ∈ P }. For the mixed integer linear set S, there is a natural linear relaxation given by

P 0 := {(x, y) ∈ R n + × R p + : Ax + Gy b}
obtained from S by discarding the integrality requirement on the vector x. The natural LP relaxation of (P) is, thus max

x,y {c T x + h T y : (x, y) ∈ P 0 } LP relaxations are used to provide bounds to the optimal solution because (i) it exists a very efficient algorithms to solve linear programs; (ii) one can generate a sequence of linear relaxations of S that provide increasingly tighter approximations of the set S.

The tightest formulation for the set S is given by its convex hull, which is the minimum convex set that contains S.

conv(S)

:= {(x, y) ∈ R n + × R p + : (x, y) is a convex combination of points in S}
Figure 1 shows a mixed integer set S and two formulations P and P for S. The convex hull of S is represented with a dashed line. The vector (c, h) indicates the direction in which one optimizes according is an approximation of the optimum value, evaluated at u. One can see that P provides a better approximation than P in the proximity of the optimum. The convex hull provides the best approximation to the mixed integer set S through linear inequalities. Unfortunately, for a given mixed integer set S, determining the linear inequalities that define the convex hull of the set can be a very difficult problem. Therefore, one strives to construct a model whose natural linear relaxation is a good approximation to the convex hull of feasible solutions and solves the natural LP relaxation of such a model to obtain a bound on the optimal solution. Once one has bounded on the solution, one of the following are standard solving schemes for a MILP can be used.

Branch-and-Bound

The Branch-and-Bound (B&B) method searches for an optimal solution to the MILP by branching, i.e. partitioning the set S into subsets, and attempting to prune the enumeration by bounding the subproblems generated by the partition. In particular, at any point during the solution process, the status of the solution is described by a pool of yet unexplored subsets and the best solution found so far. Initially only one subset exists and the best solution found so far has an infinite objective value. The unexplored subspaces are represented as nodes in a dynamically generated search tree, which initially only contains the root, and each iteration of a classical B&B algorithm processes one such node. The iteration has three main components: (i) selection of the node to process, (ii) bound calculation,(iii) and branching.

The sequence of these may vary according to the strategy chosen for selecting the next node to process. If the selection of next subproblem is based on the bound value of the subproblems, then the first operation of an iteration after choosing the node is branching, i.e.

subdivision of the solution space of the node into two or more subspaces to be investigated in a subsequent iteration. For each of these, it is checked whether the subspace contains an integer solution, in which case it is compared to the current best solution keeping the best of these. Otherwise the bounding function for the subspace is calculated and compared to the current best solution. If it can be established that the subspace cannot contain the optimal solution, the whole subspace is discarded, else it is stored in the pool of unexplored nodes together with it's bound. This is called the eager strategy by [START_REF] Clausen | On the best search strategy in parallel branch-and-bound:Best-First Search versus Lazy Depth-First Search[END_REF] for node evaluation, since bounds are calculated as soon as nodes are available.

The alternative is to start by calculating the bound of the selected node and then branch on the node if necessary. The nodes created are then stored together with the bound of the processed node. This strategy is called lazy by [START_REF] Clausen | On the best search strategy in parallel branch-and-bound:Best-First Search versus Lazy Depth-First Search[END_REF] and is often used when the next node to be processed is chosen to be a live node of maximal depth in the search tree.

The search terminates when there is no unexplored nodes left, and the optimal solution is then the one recorded as "current best".

Branch-and-Cut

In a B&B approach the tightness of the bound is crucial to prune the explored search tree. Tighter bounds can be obtained by applying cut generation at each of the subproblems in the search tree prior to branching. This leads to a Branch-and-Cut (B&C) approach. In some cases, the problem structure is such that the natural linear relaxation of (P) coincides with the convex hull of the mixed integer set S. In these cases, one refers to a perfect or ideal formulation. The main take-away from this section is the need to develop formulations for MILP problems, whose linear relaxations are as tight as possible, in the sense that they are good approximations for the convex hull of feasible solutions, at least in the vicinity of the optimum. The tighter a LP relaxation is for a MILP, the better the bound on the optimal solution. This can lead to exploring smaller search trees in the B&B solving procedure. To improve the tightness of the bound provided by the relaxation, one uses valid cuts to separate fractional points from the mixed integer set S.

Column Generation

Column generation is a general technique to solve linear programming problems with a very large number of variables. At each iteration, the method keeps a manageable subset of variables, solves the linear programming problem restricted to these variables, and either concludes that the optimal solution of the restricted problem corresponds to an optimal solution of the whole problem, or finds one or more "candidate variables" to improve the current solution.

Let us call the following linear program the Master Problem (MP).

z MP := max x,y j∈J c T λ s.t. j∈J a j λ b (MP) λ j 0, j ∈ J
with |J| = n variables and m constraints. In many applications n is exponential in m and working with the MP explicitly is not an option because of its sheer size. Instead, consider the Reduced Master Problem (RMP) which contains only a subset J ⊆ J of variables. An optimal solution λ to the RMP needs not be optimal for the MP, of course. Denote by π an optimal dual solution to the RMP. In the pricing step of the simplex method we look for a non-basic variable of negative reduced cost to enter the basis. To accomplish this in column generation, one solves the pricing problem (or subproblem).

z PP := max{c j -π a j | j ∈ J} (PP)
When z PP < 0, the variable λ j and its coefficient column (c j , a j ) corresponding to a maximizer j are added to the RMP; this is solved to optimality to obtain optimal dual variable values, and the process iterates until no further improving variable is found. In this case, λ optimally solves the MP as well. In particular, column generation inherits finiteness and correctness from the simplex method, when cycling is taken care of.

Branch-and-Price

To solve the MP, one can apply the B&B procedure, solving the MP relaxation by column generation to compute a bound at each node.

This approach is known as Branch-and-Price (B&P). B&P is similar to B&P, except that procedure focuses on column generation rather than row generation. In fact, Pricing and Cutting are complementary procedures for tightening a LP relaxation. In particular, B&P integrates B&B and column generation methods for solving large-scale Integer Programmings (IPs) . At each node of the B&B, columns may be generated to tighten (improve) the LP relaxation. In B&P, sets of columns are left out of the LP relaxation of large IPs because there are too many columns to handle efficiently and most of them will have their associated variables equal to zero in an optimal solution anyway. Then to check optimality, a sub-problem, also called the "pricing problem"

is solved to identify columns to enter the basis. If such columns are found, the LP is re-optimized. Branching occurs when no columns "price" out to enter the basis and the LP solution does not satisfy integrality conditions.

outline of the thesis

The thesis is organized as follows.

Chapter We embody this algorithm in a Branch-and-Price (B&P) algorithm for the PLOHP, whose computational performance is 10-30x faster than the previous approach described in the literature. Finally, Chapter 5 summarize the main contributions of this thesis and outline directions for future research.

A C L A S S R E P R E S E N TAT I V E M O D E L F O R P U R E PA R S I M O N Y H A P L O T Y P I N G U N D E R U N C E R TA I N D ATA
Diploid organisms, such as humans, are characterized by having the DNA organized in pairs of chromosomes, of which one copy is inherited from the father and the other from the mother. The recent completion of the sequencing phase of the Human Genome Project [START_REF] Venter | The Sequence of the Human Genome[END_REF] showed that such copies are extremely similar and that the genomes of two different individuals are identical in more than 99% of the overall number of nucleotides. Nevertheless, differences at the genomic level (also known as polymorphisms) occur, on average, every 1,000 bases [START_REF] Chakravarti | It's raining SNPs, hallelujah?[END_REF]] and are (excluding the recombination process) the predominant form of human variation as well as of genetic diseases [START_REF] Hoehe | Sequence variability and candidate gene analysis in complex disease: association of µ opioid receptor gene variation with substance dependence[END_REF][START_REF] Terwilliger | Linkage disequilibrium mapping of complex disease: fantasy or reality?[END_REF] Haplotypes represent a fundamental source of information for disease association studies. In fact, over 90% of sequence variation among individuals is due to common variant sites, most of which arose from single historical mutation events on the ancestral chromosome [START_REF] Li | Low nucleotide diversity in man[END_REF]. Hence, in a group of people affected by a disease, the SNPs causing or associated with the disease will be enriched in frequency compared with the corresponding frequencies in a group of unaffected individuals. This observation was of considerable assistance, for example, in the identification of the genes responsible for type 1 diabetes [START_REF] Bell | A polymorphic locus near the human insulin gene is associated with insulin-dependent diabetes mellitus[END_REF][START_REF] Dorman | Worldwide differences in the incidence of type I diabetes are associated with amino acid variation at position 57 of the HLA-DQ beta chain[END_REF][START_REF] Nistico | The CTLA-4 gene region of chromosome 2q33 is linked to, and associated with, type 1 diabetes. Belgian Diabetes Registry[END_REF], type 2 diabetes [START_REF] Altshuler | The common PPARgamma Pro12Ala polymorphism is associated with decreased risk of type 2 diabetes[END_REF][START_REF] Deeb | A Pro12Ala substitution in PPAR γ 2 associated with decreased receptor activity, lower body mass index and improved insulin sensitivity[END_REF]], Alzheimer's disease [START_REF] Strittmatter | Apolipoprotein E and Alzheimer's disease[END_REF], deep vein thrombosis [START_REF] Dahlbäck | Resistance to activated protein C caused by the factor V R506Q mutation is a common risk factor for venous thrombosis[END_REF]],

inflammatory bowel disease [START_REF] Hugot | Association of NOD2 leucine-rich repeat variants with susceptibility to Crohn's disease[END_REF][START_REF] Ogura | A frameshift mutation in NOD2 associated with susceptibility to Crohn's disease[END_REF][START_REF] Rioux | Genetic variation in the 5q31 cytokine gene cluster confers susceptibility to Crohn disease[END_REF], hypertriglyceridaemia [START_REF] Pennacchio | An apolipoprotein influencing triglycerides in humans and mice revealed by comparative sequencing[END_REF],

schizophrenia [START_REF] Stefansson | Neuregulin 1 and susceptibility to schizophrenia[END_REF], asthma [START_REF] Van Eerdewegh | Association of the ADAM33 gene with asthma and bronchial hyperresponsiveness[END_REF], stroke [START_REF] Gretarsdottir | The gene encoding phosphodiesterase 4D confers risk of ischemic stroke[END_REF]], myocardial infarction [START_REF] Ozaki | Functional SNPs in the lymphotoxin-alpha gene that are associated with susceptibility to myocardial infarction[END_REF], cystic fibrosis and diastrophic dysplasia [START_REF] Consortium | The international hapmap project[END_REF][Consortium , 2005]].

Haplotyping an individual therefore consists of determining, for each copy of a given chromosome region, a pair of haplotypes. Haplotyping populations of individuals has attracted an increasing amount of attention in recent years [START_REF] Helmuth | Map of the Human Genome 3.0[END_REF][START_REF] Marshall | Drug Firms to Create Public Database of Genetic Mutations[END_REF]] because of its importance in the analysis of fine-scale genetic data [START_REF] Clark | Haplotype structure and population genetic inferences from nucleotide-sequence variation in human lipoprotein lipase[END_REF][START_REF] Schwartz | Methods for inferring block-wise ancestral history from haploid sequences[END_REF]. For example, haplotypes are necessary in evolutionary studies to extract the information needed to detect diseases and to reduce the number of tests to be carried out. In functional genomics, haplotypes are used to discover a functional gene or to study an altered response of an organism to a particular therapy. In human pharmacogenetics, haplotypes explain why people react differently to different types or amounts of drugs. In fact, because SNPs affect the structure and function of proteins and enzymes, they may influence the way in which a drug is absorbed and metabolized.

Direct sequencing of haplotypes via experimental methods is both time-consuming and expensive, and therefore current molecular sequencing methods generally provide more general genotype information. Specifically, genotype data provide information about the multiplicity of each SNP allele of a given individual, i.e., knowledge about its homozygous or heterozygous nature. Unfortunately, a drawback of using genotype data is that information regarding which heterozy-gous site SNP variants came from the same chromosome copy remains unknown [START_REF] Wang | Haplotype inference by maximum parsimony[END_REF]]. Hence, in silico haplotyping methods become attractive alternatives and in some cases the only viable way for haplotyping populations [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF]].

The simplest way for haplotyping a population is described in [START_REF] Bonizzoni | The haplotyping problem: A view of computational models and solutions[END_REF], [START_REF] Bafna | Haplotyping as perfect phylogeny: A direct approach[END_REF], and [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF] and can be summarized as follows: first, to experimentally obtain genotype data, and subsequently, for each individual, to retrieve the haplotypes computationally -i.e., to find a set of haplotypes such that, if they are assumed to be the corresponding set of chromosome copies, then computing the multiplicity of each SNP allele one can obtain exactly the genotypes given. However, this approach requires the presence of some haplotyping criterion.

Several criteria have been proposed for haplotyping populations, each of them based on biological motivations (see, for example, [START_REF] Bafna | Haplotyping as perfect phylogeny: A direct approach[END_REF], [START_REF] Clark | Haplotype structure and population genetic inferences from nucleotide-sequence variation in human lipoprotein lipase[END_REF], [START_REF] Eskin | Efficient reconstruction of haplotype structure via perfect phylogeny[END_REF], [START_REF] Excoffier | Maximum Likelihood estimation of molecular haplotype frequencies in a diploid population[END_REF], [START_REF] Fallin | Accuracy of haplotype frequency estimation for biallelic loci via the expectation maximization algorithm for unphased diploid genotype data[END_REF], [START_REF] Lancia | A polynomial case of the parsimony haplotyping problem[END_REF], Niu et al. [2002a], Niu et al. [2002b], [START_REF] Stephens | A comparison of bayesian methods for haplotype reconstruction from population genotype data[END_REF], and [START_REF] Stephens | A new statistical method for haplotype reconstruction from population data[END_REF]).

In this chapter we consider the parsimony criterion [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF].

The idea at the core of the parsimony criterion is that under many plausible explanations of an observed phenomenon, the one requiring the fewest assumptions should be preferred [START_REF] Semple | Phylogenetics[END_REF]].

Hence, because the number of distinct haplotypes observed in a population is much smaller than the number of possible haplotypes, the parsimony-based approaches aim to determine the minimum number of different haplotypes that, combined in pairs over time, have given rise to a set of observed genotypes.

The problem of haplotyping populations under parsimony (hereafter denoted as the Pure Parsimony Haplotyping (PPH) problem) is known to be APX-hard [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF]]. This result justifies the development of several enumerative optimization algorithms that aim to solve exactly instances of PPH.

Specifically, [START_REF] Bafna | Haplotyping as perfect phylogeny: A direct approach[END_REF] first proposed an integer-programming model to tackle instances of PPH. He described a model, exponential in size, characterized by two kinds of variables -one for haplotypes and the other one for haplotype pairs -and by the exhaustive generation of the set of all haplotypes compatible with some genotype in the input. Similar integer programming models were also used by [START_REF] Brown | A New Integer Programming Formulation for the Pure Parsimony Problem in Haplotype Analysis[END_REF] and [START_REF] Bertolazzi | Solving haplotyping inference parsimony problem using a new basic polynomial formulation[END_REF]. To minimize the number of distinct haplotypes, [START_REF] Brown | A New Integer Programming Formulation for the Pure Parsimony Problem in Haplotype Analysis[END_REF] proposed constructing haplotype vectors by associating a variable to each site; they subsequently used constraints to establish the exact haplotype structures. On the other hand, [START_REF] Bertolazzi | Solving haplotyping inference parsimony problem using a new basic polynomial formulation[END_REF] first formulated PPH as a minimization problem characterized by a polynomial number of variables and constraints. Then the authors turned the problem into a maximization problem and strengthened the model by using clique inequalities, symmetry breaking, inequalities, and dominance relations. Whereas [START_REF] Bafna | Haplotyping as perfect phylogeny: A direct approach[END_REF] and [START_REF] Bertolazzi | Solving haplotyping inference parsimony problem using a new basic polynomial formulation[END_REF] used commercial mixedinteger programing solvers (CPLEX and Xpress-MP, respectively) to get solutions to their models, [START_REF] Brown | A New Integer Programming Formulation for the Pure Parsimony Problem in Haplotype Analysis[END_REF] used a Branch-and-Cut (B&C) algorithm to solve their polynomial model. A comparison of their results shows that the [START_REF] Brown | A New Integer Programming Formulation for the Pure Parsimony Problem in Haplotype Analysis[END_REF] polynomial model is well suited for big dimension samples, whereas the [START_REF] Bafna | Haplotyping as perfect phylogeny: A direct approach[END_REF] and [START_REF] Bertolazzi | Solving haplotyping inference parsimony problem using a new basic polynomial formulation[END_REF] models are more efficient for medium dimensions, and specifically, when the recombination level (i.e., the parameter that affects the structure of the haplotypes) increases. So far, to our best knowledge, data sets containing 68 genotypes and 75

SNPs represent the limit size instances of PPH that can be exactly ana-lyzed [Brown et al. 2006].

In the next section, we introduce some notations that will prove useful to state the PPH problem. In particular, we state the problem of detecting haplotypes from genotype data in terms of an optimization problem. Then, we extend the concepts to the general case in which SNPs are missing or affected by errors.

notation and problem statement

The human genome is divided in 23 pairs of chromosomes thereof, one copy is inherited from the father and the other from the mother.

When a nucleotide site of a specific chromosome region shows a variability within a population of individuals then it is called SNP. Specifically, a site is considered a SNP if for a minority of the population a certain nucleotide is observed (called the least frequent allele) while for the rest of the population a different nucleotide is observed (the most frequent allele) [START_REF] Aringhieri | Optimal solutions for the balanced minimum evolution problem[END_REF]. The least frequent allele, or mutant type, is generally encoded as '1', as opposed to the most frequent allele, or wild type, generally encoded as '0' [START_REF] Zhang | Models and algorithms for haplotyping problem[END_REF].

A haplotype is a set of alleles, or more formally, a string of length p over an alphabet Σ = {0, 1} [START_REF] Catanzaro | The pure parsimony haplotyping problem: overview and computational advances[END_REF].

Extracting haplotypes from a population of individuals is not an easy task. In fact, the current molecular sequencing techniques only provide information about the conflation of the paternal and maternal haplotypes of an individual (also called genotype) rather than haplotypes themselves [START_REF] Halldórsson | Combinatorial Problems Arising in SNP and Haplotype Analysis[END_REF]]. When the family-based genetic information of a population is available, haplotypes can be retrieved experimentally [START_REF] Lu | Haplotype Information and Linkage Disequilibrium Mapping for Single Nucleotide Polymorphisms[END_REF]]. However, the experimental approach is generally laborious, cost-prohibitive, requires advanced molecular isolation strategies [START_REF] Clark | Statistical estimation and pedigree analysis of CCR2-CCR5 haplotypes[END_REF]], and sometimes not even possible [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF]. In absence of a family-based genetic information, a valid alternative to the experimental approach is provided by computational methods which estimate, by means of specific criteria, haplotypes from the set of genotypes extracted from a population of individuals. denotes heterozygous sites. As an example, the sequence 0, 2, 1 encodes a genotype in which: the first SNP is homozygous of wild type;

the second SNP is homozygous of mutant type; and finally the third SNP is heterozygous. A genotype is said to be degenerate if it does not contain '1's. A genotype g k is said to be resolved from a pair of haplo-

types {h i , h j }, in symbols g k = h i ⊕ h j , if the p-th entry of g k , denoted
as g kp , is equal to the sum of the p-th entries of h i and h j , denoted as h ip and h jp , respectively. For example, the genotype 1, 1, 0 andh j = 1, 1, 0, 0 Any feasible solution to PPH is constituted by (i) a set H of (at most 2m) haplotypes and (ii) for each genotype g k ∈ G, a specification of a pair of haplotypes, say {h i , h j }, resolving g k , i.e., such that g k = h i ⊕ h j . For example, the alternative (minimum) solutions relative to the instance of PPH shown in Figure 3 satisfy both conditions (i) and

g k = 1, 2, 1, 0 is resolved from h i = 0,
(ii). However, note that, although having the same set of haplotypes, Solutions 1 and 2 of Figure 3 are different, as they satisfy condition (ii) in a different way.

Hence, a feasible solution to PPH can be represented by means of a bipartite graph in which each vertex g k ∈ G is of degree 2 and the two other vertices, say h i and h j , adjacent to g k satisfy g k = h i ⊕ h j . As an example, the bipartite graphs corresponding to Solutions 1 and 2 of (iii) every pair of subsets intersects in at most one genotype. As an example, the haplotypes of Solution 1 in Figure 3 induce the family of subsets of Figure 5a satisfying properties (i)-(iii). Specifically, the subsets are induced by the following four haplotypes:

h 1 = 1, 1, 0, 0 inducing the subset S 1 = { 1, 1, 1, 1 , 2, 1, 0, 1 }, h 2 = 0, 0, 1, 1 induc- ing the subset S 2 = { 1, 1, 1, 1 , 0, 1, 2, 1 }, h3 = 1, 0, 0, 1 inducing the subset S 3 = { 2, 1, 0, 1 }, and finally h 4 = 0, 1, 1, 0 inducing the sub- set S 4 = { 0, 1, 2, 1 }.
Similarly, the haplotypes of Solution 2 in Figure 4 induce the family of subsets of Figure 5b satisfying properties (i)-(iii).

Specifically, the subsets are induced by the following four haplotypes: 1, 1, 1 , 2, 1, 0, 1 }, andfinally h 4 = 0, 1, 1, 0 inducing the subset S 4 = { 1, 1, 1, 1 , 0, 1, 2, 1 given a genotype and denoted n as the number of its heterozygous sites, there exist 2 n-1 possible haplotypes that may resolve it [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF]. This insight entails the use of a criterion to select pairs of haplotypes among plausible alternatives. [START_REF] Gusfield | Inference of haplotypes from samples of diploid populations: Complexity and algorithms[END_REF] and [START_REF] Wang | Haplotype inference by maximum parsimony[END_REF] observed that the number of distinct haplotypes existing in a large population of individuals is generally much smaller than the overall number of distinct genotypes observed in that population.

h 1 = 1, 1, 0, 0 inducing the subset S 1 = { 2, 1, 0, 1 }, h 2 = 0, 0, 1, 1 inducing the subset S 2 = { 0, 1, 2, 1 }, h 3 = 1, 0, 0, 1 inducing the subset S 3 = { 1,
This insight has suggested that, for low-rate recombination genes at least, the criterion of minimizing the overall number of haplotypes necessary to resolve a set of genotypes may have good chances to recover the biological haplotype set. This criterion, formally introduced by [START_REF] Gusfield | Inference of haplotypes from samples of diploid populations: Complexity and algorithms[END_REF], is known as the pure parsimony criterion of haplotype estimation and was of considerable assistance, for example, in the identification of the genes responsible for psoriasis and severe alopecia areata [START_REF] Aringhieri | Optimal solutions for the balanced minimum evolution problem[END_REF]. Haplotyping a set of genotypes under the parsimony criterion involves solving an optimization problem, called pure parsimonious haplotyping problem, that can be stated as follows:

Problem 1 (The Pure Parsimony Haplotyping (PPH) problem). Given a set G of m non-degenerating genotypes, having s SNPs each, find the minimum set H of haplotypes such that for each genotype g k ∈ G there exists a pair of

haplotypes {h i , h j } ∈ H resolving g k .
PPH is known to be polynomially solvable when each genotype has at most two heterozygous sites [START_REF] Lancia | A polynomial case of the parsimony haplotyping problem[END_REF], and NP-hard when each genotype has at least three heterozygous sites [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF].

Recently, Brown et al. [2006] introduced an interesting version of PPH called the Pure Parsimony Haplotyping under Uncertain Data (PPH-UD).

This version mainly arises when the input genotype set G is not accurate, i.e., when some SNPs are missing or affected by errors, a situation that often occurs in practice. In this case, the input of the problem may include also a binary matrix B, called the error mask matrix, whose generic entry b kp is equal to 1 if the p-th SNP of genotype g k is uncertain (i.e., missing or affected by an error), and 0 otherwise. When a given SNP is uncertain its actual value could sensibly deviate from its true value. For example, the true value of a wild type homozygous SNP affected by uncertainty could be homozygous of mutant type or even heterozygous. Similarly, the true value of a heterozygous SNP affected by uncertainty could be homozygous of wild or mutant type.

The presence of uncertain data modifies the standard definition of resolution for a genotype. Specifically, Brown et al. [2006] stated that when uncertainty occurs in the input data a genotype g k is resolved by a pair of haplotypes {h i ,

h j } if g kp = h ip + h jp + b kp e kp ,

for all

SNPs p, being e kp being an integer variable taking values in the set Brown et al. [2006] described an integer programming model able to tackle instances of PPH affected by uncertain data.

E = {-2, -1, 0, 1, 2}.
Unfortunately, the authors did not offer experimental evidence of the performances of their model due to its unpractical runtimes.

In the next section we address this critical issue by introducing a possible integer linear programming model to solve exactly instances of PPH-UD. The model is based on an extension of Catanzaro et al. [2010b] Class Representative Model (CRM), currently one of the best integer programming model for PPH described in the literature. We propose several reduction rules to improve the performances of our model and show how to embody perfect phylogeny constrains in order to represent known familial relationships among members of the input population.

design and implementation

As shown in Catanzaro et al. [2010b], any feasible solution of PPH induces a family of subsets of genotype such that: (i) each subset represents one unique haplotype with elements in the subset being genotypes carrying the haplotype, (ii) each genotype belongs to exactly two subsets, and (iii) every pair of subsets intersects in at most one genotype. This principle can be exploited also when dealing with PPH-UD. Specifically, let associate an index to each subset S of genotypes induced by a haplotype h. If i is the smallest index of a genotype belonging to S, then i is the index associated to S and the subset will be denoted as S i . Since each genotype k belongs to exactly two subsets (as it must be explained by exactly two haplotypes) it may happen that k is itself the genotype with smallest index in both subsets. In this case a dummy genotype k is added, and the subset S k is created. As an example, one can imagine that the haplotype h 1 induces the subset S i = {g i , g j , g k , . . .}, h 2 induces the subset

S i = {g i , g l , g r , g s , . . .}, h 3 induces the subset S k = {g k , g l , g s , g t , . . .},
and so on. We remark that the index k can be considered only if k was previously used, i.e., if the subset S k already exists.

Since at most 2m haplotypes are necessary to resolve m genotypes [START_REF] Lancia | Haplotyping Populations by Pure Parsimony: Complexity of Exact and Approximate Algorithms[END_REF], then the indices i of the subsets S i can vary inside the index set Q = K ∪ K , where K = {1, . . . , m} and K = {1 , . . . , m }.

Assume that an order is defined on Q in such a way that 1 < 1 < 2 < 2 < . . . < m < m . Define x i , ∀ i ∈ Q, as a decision variable equal to 1 if, in the solution, there exists a haplotype inducing a subset S i of genotypes whose smallest index genotype is g i , and 0 otherwise.

Denote y k ij , ∀ k ∈ K, ∀ i, j ∈ Q, as a decision variable equal to 1 if genotype k belongs to the subsets S i and S j , and 0 otherwise. Denote P as the set of the input SNPs and z ip , ∀ i ∈ Q, p ∈ P, as a decision variable equal to 1 if the haplotype inducing the subset S i of genotypes has such a value at p-th site, and 0 otherwise. Variables z ip describe explicitly the haplotypes of the solution.

For every non-null entry of the error mask matrix B denote e c kp as a decision variable accounting for the difference between the value of g kp and the true underlying value. Specifically, e c kp is equal to 1 if the p-th entry of genotype g k is corrected with a value c ∈ C = {-2, -1, 1, 2}, and 0 otherwise. Finally, let E LB and E UB be a lower and an upper bounds on the overall number of errors in G. Then, the following model is a valid formulation of PPH-UD: translate the sum operation among haplotypes when uncertainty occurs in the input data. Specifically, constraints (1.7) account for the correction imposed on the p-th SNP of haplotypes h i and h i when g i has its p-th SNP equal to 0. In this case, two situations may occur at the p-th SNP: either no correction is performed, or a correction is performed by setting to 1 one between e 1 kp and e 2 kp . If a correction is performed and e 1 kp is set to 1 at the p-th SNP then one haplotype will be homozygous of wild type and the other of mutant type. On the contrary, if e 2 kp is set to 1 then both haplotypes will be homozygous of mutant type. Constraints (1.8) account for the correction imposed on the p-th SNP of haplotypes h i and h i when g i has its p-th SNP equal to 2. Similarly to constraints (1.7), also in this case two situations may occur: either no correction is performed, or a correction is performed by setting to 1 one between e -1 kp and e -2 kp . If a correction is performed and e -1 kp is set to 1 at the p-th SNP then one haplotype will homozygous of wild type and the other of mutant type. On the contrary, if e -2 kp is set to 1 then both haplotypes will be homozygous of wild type. Finally, constraints (1.9) account for the correction imposed on the p-th SNP of haplotypes h i and h i when g i has its p-th SNP equal to 1. If a correction is performed and e 1 kp is set to 1 at the pth SNP then both haplotypes will be homozygous of mutant type. On the contrary, if e -1 kp is set to 1 then both haplotypes are homozygous of wild type. Constraints (1.10) establish the relations between variables z is and y k ij in absence of uncertainty. Specifically, they force the p-th site of the haplotype h i to be equal to 0 when at least one genotype g k , whose p-th entry equal to 0, belongs to the induced subset S i . By analogy, constraints (1.11) force the p-th site of the haplotype h i to be equal to 1 when at least one genotype g k , whose p-th entry equal to 2, belongs to the induced subset S i . Constraints (1.12-1.13) force one of the two p-th sites of haplotypes h i and h j to be equal to 1 when the p-th entry of genotype g k is equal to 1. Constraints (1.14-1.17) are the analogous version of constraints (1.10-1.13) in presence of uncertainty in the input data. Constraints (1.18-1.20) impose that at most one variable e c kp can be equal to one in presence of uncertainty in the input data. Finally, constraints (1.21-1.22) impose the upper and lower bounds on the error variables e c kp .

Formulation 1. Class Representative Model (CRM) for PPH-UD min i∈Q x i (1.1) s.t. x i x i ∀i ∈ Q (1.2) i,j∈Q y k ij 1 ∀k ∈ K (1.3) j∈Q y k ij x i ∀k ∈ K, ∀i ∈ Q (1.4) y k kk x k ∀k ∈ K (1.5) z kp + z k p = g kp ∀k ∈ K, ∀p ∈ P : b kp = 0 (1.6) z kp + z k p = e 1 k,p + 2e 2 k,p ∀k ∈ K, ∀p ∈ P : g kp = 0, b kp = 1 (1.7) z kp + z k p = 2 + (-e -1 k,p -2e -2 k,p )∀k ∈ K, ∀p ∈ P : g kp = 2, b kp = 1 (1.8) z kp + z k p = 1 + (e 1 k,p -e -1 k,p ) ∀k ∈ K, ∀p ∈ P : g kp = 1, b kp = 1 (1.9) z ip 1 - j∈Q y k ij ∀k ∈ K, ∀p ∈ P : g kp = 0, b kp = 0, ∀i ∈ Q : i = k, k (1.10) z ip j∈Q y k ij ∀k ∈ K, ∀p ∈ P : g kp = 2, b kp = 0, ∀i ∈ Q : i = k, k (1.11) z ip + z jp y k ij ∀k ∈ K, ∀p ∈ P : g kp = 1, b kp = 0, ∀i, j ∈ Q (1.12) z ip + z jp 2 -y k ij ∀k ∈ K, ∀p ∈ P : g kp = 1, b kp = 0, ∀i, j ∈ Q (1.13) z ip 1 - j∈Q y k ij + (e 1 k,p + 2e 2 k,p ) ∀k ∈ K, ∀p ∈ P : g kp = 0, b kp = 1, ∀i ∈ Q, i = k, k (1.14) z ip j∈Q y k ij + (-e -1 k,p -2e -2 k,p ) ∀k ∈ K, ∀p ∈ P : g kp = 2, b kp = 1, ∀i ∈ Q, i = k, k (1.15) z ip + z jp y k ij + (e 1 k,p -e -1 k,p ) ∀k ∈ K, ∀p ∈ P : g kp = 1, b kp = 1, ∀i, j ∈ Q (1.16) z ip + z jp 2 -y k ij + (e 1 k,p -e -1 k,p ) ∀k ∈ K, ∀p ∈ P : g kp = 1, b kp = 1, ∀i, j ∈ Q (1.17) e 1 k,p + e 2 k,p 1 ∀k ∈ K, ∀p ∈ P : g kp = 0, b kp = 1 (1.18) e -1 k,p + e -2 k,p 1 ∀k ∈ K, ∀p ∈ P : g kp = 2, b kp = 1 (1.19) e 1 k,p + e -1 k,p 1 ∀k ∈ K, ∀p ∈ P : g kp = 1, b kp = 1 (1.20) k∈K p∈P:b kp =1 c∈C e c kp E LB (1.21) k∈K p∈P:b kp =1 c∈C e c kp E UB (1.22) x i , z ip , y k ij , e c kp ∈ {0, 1} (

Reducing model size

The particular nature of the set of indices Q can be exploited to reduce the size of the CRM for PPH-UD. As observed in Catanzaro et al. [2010b], this operation proves fundamental to improve the efficiency of whole model.

Given that y k ij = 1 if and only if k belongs to two subsets having g i and g j for smallest index genotype, we need to define variables y k ij only when i < j k or i = k and j = k . For example, variable y 2 1,1 does not need to be defined as well as all variables y k ii for all i, k ∈ K, k = i. Similarly, variables y k ik or y k k i (depending on whether k is smaller or bigger than i) do not need to be defined for i ∈ Q with i = k and i = k . In fact, if y k ik = 1, then k belongs to two subsets, one represented by i and the other one by k , which contradicts the assumption that the dummy genotype k can be considered only if k is already used. By extending this analysis to all the possible cases in which variables y k ij are redundant and assuming that variable y 1 11 = 1, it is easy to see that variables y k ij do not need to be defined whenever they belongs to one of the following sets:

R 1 = {y k ij : k ∈ K, i, j ∈ K ∪ K , j < i < k}, (2)
R 2 = {y k ik : k ∈ K, i ∈ K ∪ K , i (k -1) }, (3)
R 3 = {y k ii : k ∈ K, i ∈ K ∪ K , 2 i k -1}. ( 4
)
The sets of redundant variables can be further expanded by observing that for each triplet of genotypes {g i , g j , g k } such that the respective p-th SNP is g ip = 0, g jp = 0, g kp = 1, and b ip = b jp = b kp = 0, variable y k ij is necessarily equal to 0 since the containment of genotype g k to the subsets S i and S j would violate the sum operator among haplotypes at least at p-th SNP. Extending this argument to all the possible combinations of triplets of genotypes that violate the haplotype sum operator, we see that the following proposition holds:

Proposition 1. The set of variables R 4 = {y k ij : i, j, k ∈ K, p ∈ P, g kp = 1, (5)
g ip = g jp = 1, b kp = b ip = b jp = 0} is redundant.
Proof. By contradiction, if the set R 4 is not redundant, then for some p ∈ P, genotype k may belong to two subsets S i and S j . Let assume g kp = 1 and g ip = g jp = 0 and that uncertainty does not affects the p-th SNP. Since haplotype h i is associated with the subset S i and explains genotype i, it must have the p-th SNP equal to 0 otherwise the sum operator would be violated. At the same time, since haplotype h j is associated with the subset S j and explains genotype j, it must have the p-th SNP equal to 0 otherwise the sum operator would be violated. In turn, this implies that genotype k cannot be resolved by h i and h i since a necessary condition for its resolution is that the p-th SNP of h i or h j be equal to 1. In the second case, let assume g kp = 1

and g ip = g jp = 2 and that uncertainty does not affects the p-th SNP. Since haplotype h i is associated with the subset S i and explains genotype i, it must have the p-th SNP equal to 1 otherwise the sum operator would be violated. At the same time, since haplotype h j is associated with the subset S j and explains genotype j, it must have the p-th SNP equal to 1 otherwise the sum operator would be violated.

In turn, this implies that genotype k cannot be resolved by h i and h i since a necessary condition for its resolution is that the p-th SNP of h i or h j be equal to 0.Hence, y k ij = 0 in any feasible solution of PPH.

Similar arguments can be used to prove the following proposition:

Proposition 2. The sets of variables

R 5 = {y k ij : i, j, k ∈ K, p ∈ P, g kp = 0, (6)
g ip = 2 or g jp = 2, b kp = b ip = b jp = 0} R 6 = {y k ij : i, j, k ∈ K, p ∈ P, g kp = 2, (7)
g ip = 0 or g jp = 0, b kp = b ip = b jp = 0} are redundant.
Proof. By contradiction, if the set R 5 is not redundant, then for some p ∈ P, genotype k may belong to two subsets S i and S j . Without loss of generality, assume g kp = 0 and g ip = 2 and that uncertainty does not affects the p-th SNP. Since haplotype h i is associated with the subset S i and explains genotype i, it must have the p-th SNP equal to 1 otherwise the sum operator would be violated. In turn, this implies that genotype k cannot be resolved by h i since a necessary condition for its resolution is that the p-th SNP of h i be equal to 0. Hence, independent of j, y k ij = 0 in any feasible solution of PPH. A similar approach can be used to prove that R 6 is redundant.

Note that, removing the redundant variables y k ij of Propositions 1 and 2 can be performed in O(m 3 s). Finally, a similar process of reduction can be applied to variables z ip both by removing those whose value is fixed by constraints (1.6) (e.g., when g kp = 0 or when g kp = 2). In this way, only variables z ip involved in constraints (1.6) when g kp = 1 and in constraints (1.7-1.9) need to be defined.

Accounting for perfect phylogeny constraints

In some practical situations it may turn out useful to know not only the haplotype set that resolves a given set of genotypes, but also how the familial relationships among members of a population of individuals reflect on the estimated haplotypes. This version of PPH is known in the literature as the Minimum Perfect Phylogeny Haplotyping problem (MPPH) problem and was first introduced by [START_REF] Bafna | A note on efficient computation of haplotypes via perfect phylogeny[END_REF].

MPPH requires to find a minimum set of haplotypes resolving the input genotypes and forming a perfect phylogeny, i.e., a haplotype set H that does not contain the following matrix [Brown et al. 2006;Pe'er et al. 2004]:

F =           0 0 0 1 1 0 1 1           .
Excluding matrix F from the haplotype set H is equivalent to impose the following constraints on z ip variables:

α 1 z ip 1 + α 2 z ip 2 + α 3 z jp 1 + α 4 z jp 2 + α 5 z kp 1 + α 6 z kp 2 + α 7 z qp 1 + α 8 z qp 2 3 ∀ i = j = k = q ∈ Q, p 1 = p 2 ∈ P (8) 8 s=1 α s = 0 (9) α s ∈ {-1, 1} (10)
Adding ( 8), ( 9) and ( 10) to the class representative model could drastically slow down any mixed integer linear programming solver due to the large number of constraints to be added. For this reason, these constraints should be imposed dynamically (e.g., by means of a B&C algorithm). Note also that not all instances of PPH-UD admit a perfect phylogeny. For example, if the genotype set G contains a sub-matrix

Ω =           0 0 0 2 2 0 2 2          
whose entries are not affected by uncertainty, then we can conclude that G does not admit solution. The same situation persists if any row of Ω (possibly all) contains a heterozygous site. Thus, when a perfect phylogeny is required, preprocessing techniques based on sufficient conditions (such as those described in [START_REF] Bafna | Haplotyping as perfect phylogeny: A direct approach[END_REF] able to exclude a-priori the existence of a solution satisfying a perfect phylogeny could be used to speed up the computation. Testing the performances of the CRM for PPH-UD under the perfect phylogeny requirement is out of the scope of the present thesis and warrants additional analysis.

experimental results

In this section we analyze the performances of the class representative model to solve instances of the pure parsimony haplotyping problem under uncertain data. Similar to Brown et al. [2006] and Catanzaro et al. [2010b], we emphasize that our experiments aim simply to evaluate the runtime performance of our model for solving PPH. We neither attempt to study the efficiency of PPH for haplotype inference nor compare the accuracy of our algorithm to haplotype inference solvers that do not use the parsimony criterion. This analysis has been already performed by [START_REF] Bafna | Haplotyping as perfect phylogeny: A direct approach[END_REF], [START_REF] Wang | Haplotype inference by maximum parsimony[END_REF], and [START_REF] Marchini | A comparison of phasing algorithms for trios and unrelated individuals[END_REF], and we refer the interested reader to their respective articles.

Implementation

The results have been obtained by implementing the CRM for PPH-UD in Mosel 2.0 of Xpress-MP, Optimizer version 18, running on a Pentium 4, 3.2 GHz, equipped with 2 GByte RAM and operating system Gentoo release 7 (kernel linux 2.6.17). In our experiments we activated the Xpress-MP Optimizer automatic cuts, the Xpress-MP presolving strategy, and used the Xpress-MP primal heuristic to generate the first upper bound.

Benchmark instances

As in Catanzaro et al. [2010b], we used the standard Brown et al. In Tables 1-3 we show the performances of the CRM for PPH-UD under different error ratios. Specifically, the columns of Tables 1-3 evidence the average, the maximum, and the minimum of: the solution time, the gap (i.e., the difference between the optimal value found and the value of linear relaxation at the root node of the search tree, divided by the optimal value), and the number of nodes expanded in each group of instances belonging to a given dataset.

Computational performances

In order to obtain a qualitative measure of the running time perfor- In terms of accuracy, the performances of CRM2 in the uniform datasets are very good. Specifically, the average accuracy is over 90%

in the majority of the analyzed datasets and independently from the error ratio. However, it is worth noting that in some instances the accuracy may decrease sensibly (see, e.g., datasets 30x50 and 30x75) and proportionally to the increment of the error ratio, by suggesting, as general trend, the fact that the higher the error ratio the more difficult is to recover the correct haplotype set.

Nonuniform Datasets

The general trend observed in the uniform datasets persists also in the nonuniform datasets. Specifically, as shown in In terms of accuracy, the performances of CRM2 in the nonuniform datasets are still good, but slightly poorer than in the uniform datasets. Specifically, the average accuracy is over 80% in the majority of the analyzed datasets and independently from the error ratio. Similarly to the uniform datasets, in some instances the accuracy may decrease sensibly (see, e.g., datasets 30x75). However, in the worst case, the decrease turns out to be much smaller than the corresponding one in the uniform datasets.

Biological Datasets

To complete the performance analysis on Brown et al. [2006]'s datasets, we tested CRM2 on the biological datasets. Once again, the general trend observed in the uniform and nonuniform datasets persists also in the biological datasets: CRM2 is sensibly slower than CRM1, a part from datasets CHR10-CEU and CHR21-CEU in which the trend is inverted due to the peculiar nature of both datasets. While the average gap of CRM1 never overcomes 2.6%, the average gap of CRM2 overcomes 10% confirming the hardness of the biological datasets.

However, we stress once again the fact that PPH and PPH-UD are de facto two different problems, hence intrinsic values such as the gap cannot be directly compared. Our analysis just aims at offering experimental evidence of the tightness of the class representative model in tackling instances of the pure parsimony haplotyping problem under uncertain data.

The small number of instances constituting each biological dataset (three instances per dataset) prevents a clear statistical characterization of the performances of CRM2 in terms of accuracy. As general trend, we have observed that the accuracy approaches 100% in the majority of the biological datasets analyzed. Nevertheless, in a number of datasets this trend changes, leading the accuracy level to low values. Investigating the reason why this phenomenon arises and the possible corresponding remedies warrants additional analysis.

In the next Chapter we introduce another problem arising when the genotypes of a child and the ones of his two parents are known and a deletion polymorphism may be observed as an event on the child chromosome, which is one of the main factors giving rise to several human diseases.

A B R A N C H -A N D -P R I C E A L G O R I T H M F O R T H E PA R S I M O N I O U S L O S S O F H E T E R O Z Y G O S I T Y P R O B L E M
Single Nucleotide Polymorphisms (SNPs), together with the recombination process, constitute the predominant form of human variation and are as such source of human disease [START_REF] Li | Low nucleotide diversity in man[END_REF][START_REF] Wang | Largescale identification, mapping, and genotyping of single-nucleotide polymorphisms in the human genome[END_REF][START_REF] Cargill | Characterization of single-nucleotide polymorphisms in coding regions of human genes[END_REF][START_REF] Halushka | Patterns of single-nucleotide polymorphisms in candidate genes for bloodpressure homeostasis[END_REF][START_REF] Terwilliger | Linkage disequilibrium mapping of complex disease: fantasy or reality?[END_REF][START_REF] Hoehe | Sequence variability and candidate gene analysis in complex disease: association of µ opioid receptor gene variation with substance dependence[END_REF]Catanzaro et al. 2010a]. A large number of other types of variations exist in nature, including insertions, inversions, translocations. One type of variation being deletions, which occur when a subsequence of the human genome is present in a reference genome but is not in the genome of an individual being analyzed [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]. When the genotypes of a child and its two parents are known, a deletion polymorphism may be observed as a Loss of Heterozygosity (LOH) event on the child chromosome. In fact, the laws of Mendelian inheritance state that each individual inherits one copy of a chromosome from the father and one from the mother. Hence, for a given SNP, an individual can be either homozygous, i. e., the nucleotides of the parental DNA strands are equal, or heterozygous, i. e., the nucleotides of the parental DNA strands are different. For example, the first individual in Figure 6 is homozygous at the first SNP and heterozygous at the second SNP. When a deletion polymorphism occurs, an individual carries only a single copy of the chromosomal segment while the other is missing. As an example, the first individual in Figure 6 carries a deletion at the third SNP of the considered chromosome region (denoted by the symbol '-'). A deletion is de novo if it is only found in the child and not in one of its parents. Otherwise,

Population Molecular Sequence of a specific chromosome region

Father (1st DNA strand) These SNPs are located in the three variant sites. The symbol '-' represents a deletion, i. e., a lack of a nucleotide the deletion is said to be inherited i. e., passed from one of the two parents to the child. If the deletion event modifies the heterozygosity of an individual at a given site of a chromosomal region then we say that a LOH event occurred at that site [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]].

… A C A C G G T C … T T C G G G G T C … A G T C G A C C G … Father (2st DNA strand) … A C A C G G T C … T T C G A G G T C … A G T C -A C C G … Mother (1st DNA strand) … A C A T G G T C … T T C G G G G T C … A G T C A A C C G … Mother (2st DNA strand) … A C A C G G T C … T T C G G G G T C … A G T C A A C C G … Child (1st DNA strand) … A C A T G G T C … T T C G G G G T C … A G T C A A C C G … Child (2st DNA strand) … A C A C G G T C … T T C G A G G T C … A G T C G A C C G …

Deletions may have noxious effects on an individual as they may

give rise to several human diseases. For example, recent studies showed that schizophrenia [START_REF] Stefansson | Large recurrent microdeletions associated with schizophrenia[END_REF] Detecting deletions across the genome of individuals could be of fundamental assistance for the diagnosis and the treatment of certain human diseases, hence increasing research efforts have been dedicated

to this task in recent years [START_REF] Frazer | A second generation human haplotype map of over 3.1 million SNPs[END_REF].

A natural approach to perform the task consists of comparing the genomes of a given population of affected individuals with the genomes from a population of unaffected ones. However, the genomes of the Halldórsson et al. [2011a] observed that detecting deletions from genotype data may not be straightforward due to the limit of both genotyping technology and the presence of uncertainty in the genotyping process. In fact, SNP genotyping technology is not able to discern easily the difference between a homozygous site and a deletion, hence the output will always be a homozygous SNP even if the true genotype of the individual may carry only a single copy of the genotype. Moreover, even if a deletion polymorphism were observed in molecular data, such an event could be due either to the presence of real deletions or to genotyping errors, i. e., misreadings caused by the genotyping technology. [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] proposed a possible approach to address these major limitations. In particular, the authors introduced the parsimony criterion (see [START_REF] Albert | Parsimony, phylogeny, and genomics[END_REF]) to identify chromosomic regions from a population of individuals that underwent to a massive loss of heterozygosity events. Subsequently, they stated hour of computing time. Before describing the solution approach, in the next section we introduce some useful notation and definitions.

notation and problem formulation

Consider a trio t, i. e., a set of two parents and an offspring (see Figure 7), and denote s as a given SNP genotyped in t. Then, one of the following three situations may occur:

1. The SNP s can be Inconsistent whit a Loss of Heterozygosity (ILOH), a situation that occurs when the child is heterozygous. In this

Individuals SNPs

Father (1st DNA strand) case the two alleles must have been inherited from each parent.

… A A C … Father (2st DNA strand) … A A C … Mother (1st DNA strand) … G A A … Mother (2st DNA strand) … G C A … Child (1st DNA strand) … A A C … Child (2st DNA strand) … G A C …
For example, this is the case for the first highlighted column of the sequences of the trio shown in Figure 7.

2. The SNP s can be Consistent with a Loss of Heterozygosity (CLOH), a situation that occurs when a deletion may (but needs not) be introduced to explain the trio's inheritance pattern. Confer the second highlighted column of the sequences of the trio shown in Figure 7, the SNP of the child could be explained by means of a deletion of the paternal pattern.

3. The SNP s can show Evidence of a Loss of Heterozygosity (ELOH), a situation that occurs when a deletion or a genotyping error are the only possible explanation for the trio inheritance pattern.

Confer the third highlighted column of the sequences of the trio shown in Figure 7, the SNP of the child can be explained only by means of a deletion of the maternal pattern.

According to the above definitions, a trio genotyped at m SNPs can be encoded as a string of length m over an alphabet Σ = {1, 0, X},

where '1' codes for a SNP inconsistent with having a loss of heterozy-gosity; '0' codes for a SNP consistent with a loss of heterozygosity;

and 'X' codes for a SNP showing evidence of loss of heterozygosity [Halldórsson et al. 2011a]. For example, the string t = X100X0X010 represents a trio genotyped at 10 SNPs, thereof 5 consistent with having a loss of heterozygosity, 3 showing evidence of a loss of heterozygosity, and 2 inconsistent with having a loss of heterozygosity. As shown in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF], not all of the ELOH-SNPs in a trio are relevant for disease association studies. In particular, the ones that are relevant are only those that are in common with other trios. The identification of the relevant ELOH-SNPs in a set T of n trios genotyped at m SNPs can be performed by solving a NP-hard combinatorial optimization problem that can be stated as follows.

Consider a trio t ∈ T and let p be an ELOH-SNP in t. Since a trio t can also be seen as a string of length m over Σ, with a little abuse of notation we also denote p as the position of the p-th character in t. For example, given the trio t = X100X0X010 , p = 5 denotes both the second ELOH-SNP in t and the fifth character in the string t. Let I t p = {(l p , r p ), x p } be a point-interval for a ELOH-SNP p ∈ t, i. e., a pair constituted by an interval (l p , r p ) : 0 < l p < r p < m + 1 and a point

x p ∈ (l p , r p ). Let us associate a point-interval I t p to each ELOH-SNP in t, in such a way that x p denotes an ELOH-SNP located at the p-th character of t and r p (respectively l p ) represents the position of the closest ILOH-SNP on the right (respectively left) of x p . Whenever an ambiguity may arise, we shall write x t p to denote that the point x p is associated to the point-interval I t p . If such an ILOH-SNP on the right (respectively left) is missing, we set r p = m + 1 (respectively l p = 0).

For example, given the trio t = X100X0X010 we can associate to the first ELOH-SNP the point-interval I t 1 = {(0, 2), 1}; to the second ELOH-SNP the point-interval I t 5 = {(2, 9), 5}; and to the third ELOH-SNP the point-interval I t 7 = {(2, 9), 7}. The graph G is a Max-Point Tolerance Graph (MPTG) [START_REF] Catanzaro | Max-Point-Tolerance Graphs[END_REF], i. e., it is a generalization of Interval Graphs (IGs) (namely, intersection graphs of intervals on a line [START_REF] Booth | Testing for the consecutive ones property, interval graphs, and graph planarity using PQ-tree algorithms[END_REF]) in which each pair of intervals can "tolerate" a non-empty intersection (without forming an edge) as long as both of their distinguished points are not contained in this intersection. For example, Figure 9 shows the max-point-tolerance graph induced by the point-intervals shown in Figure 8. Combinatorial optimization problems that are NP-hard in general often turn out to be solvable in polynomial time on IGs (see e.g., [START_REF] Fishburn | Interval Orders and Interval Graphs:A Study of Partially Ordered Sets[END_REF] and [START_REF] Fulkerson | Incidence matrices and interval graphs[END_REF]). However, this fact does not hold in general on MPTGs (see [START_REF] Catanzaro | Max-Point-Tolerance Graphs[END_REF]).

In particular, the PLOHP has been shown to be NP-hard [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF].

In this specific application, the cost function f is such that f(S) = ψ(|S|), for any S ⊆ V and for some ψ : [0 . . . |V|] → R + 0 . In particular, both [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and Halldórsson et al. [2011a]

defined f as f(C) =          c 1 if |C| = 1 c 2 if |C| 2, ∀ C ∈ C(G) (11)
where c 1 and c 2 are two constants such that 0 < c 1 c 2 . Generally, there are different costs depending on the SNP and deletion should be considered. In fact, genotyping technologies are usually characterized by a high variability in the quality of the SNP genotypes produced Consortium 2003. A common method for dealing with these is to remove from analysis markers that show many ELOH events [START_REF] Conrad | A high-resolution survey of deletion polymorphism in the human genome[END_REF], this method however may remove most of the signal from the data in the preprocessing step. Similarly, different regions in the genome may have different propensity for carrying deletions [START_REF] Korbel | PEMer: a computational framework with simulation-based error models for inferring genomic structural variants from massive paired-end sequencing data[END_REF]. Based on Catanzaro et al. [2013] analysis about different values of c 1 and c 2 pursuit the assumption made by Halldórsson et al. [2011a] in which functions f(c) always assign the same cost to each ELOH and error, respectively. Remark that the approach developed here can be extended to any general linear cost set function.

The NP-hardness of the PLOHP justifies the development of solution approaches such as those described in the next section.

a branch-and-price algorithm for the plohp

In this section we describe an exact solution algorithm for the PLOHP based on Branch-and-Price (B&P). Throughout the thesis we assume, without loss of generality, that the cost function f in the PLOHP is a sizedefined submodular set function, i.e., it is such that f(S) = ψ(|S|), for any subset S ⊆ V and for some ψ : [0 . . . |V|] → R + 0 . This characterization of f is commonly assumed in loss of heterozigosity studies [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]] and allows to solve the PLOHP in polynomial time when the MPTG is reducible to an interval-graph [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF]]. As the number of maximal cliques in a MPTG is polynomially bounded and a polynomial time algorithm to determine all of them exists (see [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]), a polynomial sized formulation for the PLOHP can be obtained as follows. Consider a vertex v ∈ V and let

An integer linear programming formulation for the PLOHP

C v = {C ∈ C(G) : v ∈ C}. Let x v be a decision variable equal to 1 if vertex v ∈ V forms
z BIP = min C∈C(G) |C| 2 f(C)y C + v∈V f(v)x v (12a) s.t. C∈C v (G) |C| 2 y C + x v 1 ∀ v ∈ V (12b) x v ∈ {0, 1} ∀ v ∈ V (12c) y C ∈ {0, 1} ∀ C ∈ C(G), |C| 2. ( 12d
)
The

z RMP = min C∈C τ (G) |C| 2 f(C)y C + v∈V f(v)x v (13a) s.t. C∈C τ v (G) |C| 2 y C + x v 1 ∀ v ∈ V (13b) x v 0 ∀ v ∈ V (13c) y C 0 ∀ C ∈ C τ (G), |C| 2. ( 13d
)
Note that in the RMP constraints x v 1, for all v ∈ V, and y C 1, for all C ∈ C τ (G), |C| 2, have been omitted as they are redundant.

The RMP can be solved by column generation techniques. To this end, denote µ v as the dual variables associated with constraints (13b). The dual of the RMP is characterized by the following dual constraints:

µ v f(v) ∀ v ∈ V (14a) v∈V:v∈C µ v f(C) ∀C ∈ C τ (G) : |C| 2. ( 14b
)
A variable with negative reduced cost to be added to the RMP corresponds to a dual constraint violated by the current optimal solution.

Since all of the x v variables are present in the RMP, constraints (14a)

will never be violated. In contrast, constraints (14b) are violated if

∃ Ĉ ∈ C(G) : | Ĉ| 2 and v∈V:v∈ Ĉ µ v > f( Ĉ). ( 15
)
Checking whether (15) holds in the current optimal solution to the RMP involves solving a maximum node-weighted clique problem in G with weights {µ}. This task can be performed e.g., by means of Algorithm 2 that is described in Section 3.2.2. Algorithm 1 implements the column generation technique by using Algorithm 2 as pricing oracle. Solve the RMP 3:

Get the dual vector µ associated to (13b)

4:

Set weight µ v on each vertex v ∈ V

5:

Ĉ =PricingOracle(G,µ)

6: if Ĉ = ∅ then 7:
Add y Ĉ and its column to the RMP

if (p i p k p j ) ∧ (l k < p i ) ∧ (r k > p j ) then 6: C = C ∪ {p k } 7: for all v ∈ V : v ∈ C do 8: cost = (cost + µ v ) 9:
if cost > max then Each vertex v in the clique is connected to these two vertices and its corresponding point-interval is such that

p v ∈ [p l , p r ] ⊆ [l v , r v ].
Hence, a clique can be defined by the leftmost and rightmost vertices, respectively. Then, a possible approach to enumerate the polynomial bounded number of cliques in G consists of picking any pair of adjacent vertices in G (i.e., distinct point-intervals in I, say (I u , p u ) and

(I v , p v ), such that I u ∩ I v = (l u , r u ) ∩ (l v , r v ) ⊇ {p u , p v })
and searching for the largest set of points that fall inside

[p i = min{p u , p v }, p j = max{p u , p v }],
i.e., the largest set of point-intervals I k that satisfy con-

dition (p i p k p j ) ∧ (l k < p i ) ∧ (r k > p j ).
Algorithm 2 outlines the pseudo code necessary to perform the overall task. We note that the complexity of Algorithm 2 is O(|I| 3 ).

Initializing the columns of the RMP

In order to initialize the RMP it is necessary to find a subset C τ (G) ⊂ C(G) covering all of the vertices in G and such that any clique in C τ (G) has cardinality greater than or equal to 2. This is due to the presence of all the x v variable in the RMP. A possible approach to carry out this task is outlined in Algorithm 3. In particular, the algorithm takes it as input the set I in non-decreasing order. Subsequently if compatible(I i , I j , I k ) then 6:

Q = Q \ I k 7: C = C ∪ I k 8: C = C ∪ C 9: return C 10: function compatible(I u , I v , I z ) 11: return (p u p z p v ) ∧ (l z < p u ) ∧ (r z > p v )

Branching rules

The optimal solution to the RMP can be fractional. To ensure its integrality we embody the previously described column generation approach in a B&B algorithm, in which the branching is performed on both x and y variables. Specifically, concerning x variables, the algorithm branches according to the most fractional variable rule (see [START_REF] Nemhauser | Optimization[END_REF]). The same rule also applies to y variables with the exception that a y variable that has been set to 0 in a previous node of the search tree cannot be set to 1 in one of its descendants.

To ensure the application of this rule we perform the following steps:

(1) we store all the maximal clique generated in each node;

(2) during the column generation process, we filter, from among all of the maximal cliques, the ones associated to the y variables that have been set to 0 in the parent node;

(3) we propagate the filtered set to all the descendant nodes.

Preprocessing

In this section we develop a number of preprocessing techniques both to reduce the size of an instance of the PLOHP and to decompose an instance into a family of independent subproblems having a smaller size than the original one. We discuss the computational impact of these techniques in Section 3.3.

Removing Vertices

A basic approach to reduce the size of an instance of the PLOHP consists of removing variables from CLHF whose value can be known apriori and consequently remove the constraints associated with those variables. This is possible e.g., when dealing with isolated vertices in the graph. In fact, any feasible solution to CLHF is characterized by having x v = 1 for any isolated vertex v ∈ V. Isolated vertices can be efficiently detected by inspecting the set I for intervals having nocompatibility between each other. For example, the point-interval I v induces an isolated vertex in V if there does not exists in I a distinct point-interval I u such that (l u p v r u ) ∧ (l v p u r v ). It is easy to see that the removal of isolated vertices from an instance of the PLOHP can be performed in O(|I|(|I| -1)/2).

Interestingly, isolated vertices are not the only vertices that can be removed from an instance of the PLOHP. For example, given a point-

interval I v = [(l v , r v ), p v ]
∈ I, consider the situation in which there exists another interval

I u = [(l u , r u ), p u ] ∈ I such that l u = l v , p u = p v , r u = r v .
In such a case, we say that I u is a cloned interval of I v .

Let Γ I v be the set of cloned intervals of I v . Then, it is easy to see that we can remove from V the set of vertices induced by Γ I v and keep just I v as representative of the induced equivalence class Γ I v .

The identification of the cloned intervals can be performed in time O(|I|(|I| -1)/2) and can be carried out together with the removal of isolated vertices.

Decomposing a MPTG into connected components

The particular nature of point-intervals may give rise to MPTG with several connected components. For example, the point-intervals shown in Figure 8 leads to to a MPTG with by two connected components (see Figure 9). The potential existence of connected components in a MPTG can be exploited both to reduce the time necessary to identify a maximal clique in G and to decompose the problem into independent subproblems of smaller size. A possible way to quickly identify connected components in a MPTG consists of ordering the point-intervals according to a specific precedence rule. In particular, given two distinct point-intervals I u , I v ∈ I, we say that I v precedes I u , in symbols I v ≺ I u , if one of the following cases is verified: 1.

l v < l u ; 2. if l v = l u then p v < p u ; 3. if l v = l u and p v = p u then r v r u .
For example, Figure 10 shows the result of the application of the precedence rule on the point-intervals shown in Figure 8. It is easy to realize that, if the point-intervals have been sorted according to the above precedence rule, a connected component arises in G whenever there exist two I u = {(l u , r u ), p u } and

I v = {(l v , r v ), p v } such that p v l u .
Algorithm 4 provides a heuristic to identify connected components in a MPTG. In fact condition p i < l i+1 is sufficient to identify a set of nodes not connected with a given node, but not necessary. We opted for a heuristic approach rather then an exact one because in preliminary experiments we observed that the exact approach was slower than the heuristic one and did not prove particularly beneficial in terms of identification of connected components. The algorithm takes as input the set I sorted according to the above precedence rule. Subsequently, the algorithm visits each point-interval in I and constructs a connected component B of G by collecting all of the point-intervals I i such that I i does not satisfy the condition p i l i+1 . It is easy to see that the computational complexity of Algorithm 4 is O(|I|). Moreover, it is worth noting that the decomposition of G into connected components enables the parallelization of the RMP, an approach that may vastly shorten the solution time of an instance of the problem.

experimental results

In this section we analyze the performance of the B&P algorithm to solve instances of the parsimonious loss of heterozigosity problem.

Our experiments were motivated by a number of goals, namely: to compare the performance of the B&P algorithm with the ones obtained by CLHF [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF], which currently is the best exact algorithm for the PLOHP; to evaluate the benefits obtained by using the presolving and block decomposition strategies previously described; and finally, to allow the analysis of larger data sets with respect to the ones currently analyzed. As in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and Halldórsson et al. [2011b], we emphasize that our experiments simply aim to evaluate the computational performance of the exact algorithms. We neither attempt to study the efficiency of the B&P al-

Algorithm 4 | Connected Components Finder Input: sorted I Output: K: set of connected components of G 1: i = 0, K = ∅, K = ∅ 2: while i < |I| do 3: if p i > l i+1 then 4: K = K ∪ {I i } 5: else 6: K = K ∪ {K}, K = ∅ 7: return K
gorithm to predict LOH events across the genomes of a population of individuals nor to compare its accuracy versus LOH predictors that do not use the parsimony criterion. This analysis has been performed by Halldórsson et al. [2011b], [START_REF] Conrad | A high-resolution survey of deletion polymorphism in the human genome[END_REF], and [START_REF] Speed | Identification of Deletion Polymorphisms from Haplotypes[END_REF] and [START_REF] Mccarroll | Integrated detection and population-genetic analysis of SNPs and copy number variation[END_REF] and we refer the interested reader to these articles for further information.

Implementation

We implemented CLHF in ANSI C++ by using FICO Xpress 7.6, Optimizer libraries v26.01.04 (64-bit Hyper capacity). We implemented the B&P algorithm in ANSI C++ by using SCIP Optimization Suite

Calibrating the pricing oracle

Determining the number of variables with negative reduced cost to be added in the RMP is crucial for the performances of the proposed B&P. To perform this task, we generated a learning set of instances of the PLOHP constituted by 100 random instances. We generated such a set by using the same approach and parameters discussed in Halldórsson et al. [2011b]. Each random instance has between 35725 and 66670 point-intervals and between 16915 and 52522 cliques on the real segment [1, 1000]. Table 6 reports on the main characteristics of this set. We introduced a parametrized number of violated cliques to be added by the pricing oracle at runtime. Such a number is a fraction of the overall number of maximal cliques contained in a random instance and we set it to be equal to one of the following seven different percentages: 1%, 3%, 5%, 7%, 10%, 15% and 20%. rank test with Holmes correction for multiple tests [START_REF] Fay | Wilcoxon-Mann-Whitney or t-test? On assumptions for hypothesis tests and multiple interpretations of decision rules[END_REF]].

Table 7 reports on the p-values so obtained. Specifically, the presence of a value smaller then 0.05 in a cell means that the sample in the row is statistically significant smaller than the sample in the column. In contrast, the presence of a value bigger than 0.05 in a cell means that the sample in the row is statistically significant bigger than the sample in the column. In particular, Table 7 shows that adding the 1% of the violated cliques always provide a smaller solution time compared to the others approaches. For this reasons, we decided to use the setup in applying the B&P algorithm to the remaining experiments.

For the sake of completeness, we should say that we also performed tests using others criteria for selecting the variables to add at the RMP.

In particular, we tried to add at each iteration only the most violated clique but the performance results to be really poor. We also tried to randomly select one violated clique at each iteration avoiding the sort-ing of the variables by their negative reduced cost, but unfortunately also these test performed poorly.

Benchmark instances

In order to compare the performance of the B&P algorithm with the ones obtained by CLHF [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF], we considered the instances of the PLOHP described in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and Halldórsson et al. [2011b]. These instances, hereafter denoted as "gens-1", "gens-2", "gens-3", "gens-4" and "gens-5", are characterized by hav- and γ = [1, 3575], respectively. The authors obtained this way 15 instances of the PLOHP, that are hereafter grouped into three datasets denoted as "gens-x α ", "gens-x β " and "gens-x γ ", x ∈ {1, 5}, respec-1% 3% 5% 7% 10% 15% 20% 1% -1.97784 × 10 -18 1.97784 × 10 -18 1.97784 × 10 -18 1.97796 × 10 -18 1.97796 × 10 -18 1.97796 × 10 -18 3% 1 -2.95093 × 10 -17 1.97796 × 10 -18 1.97796 × 10 -18 1.97784 × 10 -18 1.97796 × 10 -18 tons ranging from 0 to 5 and over 70% of cloned point-intervals. The instances "gens-x β ", x ∈ {1, 5} are characterized by a number of connected components ranging from 33 to 40, a number of singletons ranging from 3 to 7 and over 70% of cloned point-intervals. Finally, the instances "gens-x γ ", x ∈ {1, 5} are characterized by a number of connected components ranging from 40 to 61, a number of singletons ranging from 2 to 13 and over 80% of cloned point-intervals.

5% 1 1 - 1.79421 × 10 -17 2.03865 × 10 -18 1.97796 × 10 -17 1.97784 × 10 -18 7% 1 1 1 - 6.37311 × 10 -18 1.97784 × 10 -18 1.97796 × 10 -18 10% 1 1 1 1 - 4.83732 × 10 -17 2.16549 × 10 -18 15% 1 1 1 1 1 - 2.4767 × 10 -18 20% 1 1 1 1 1 1 -
We also considered a set of 5 artificial datasets of the PLOHP, hereafter denoted as "Sim-1", "Sim-2", "Sim-3", "Sim-4" and "Sim-5", each containing 10 random instances of the problem and mainly differing from one another both by the length and the number of the point-intervals contained in each of the corresponding instances. We generated the instances in each dataset by first fixing the real segment [1, 6000] as a reference and by generating, by means of the Mersenne Twister library [START_REF] Matsumoto | Mersenne twister: A 623dimensionally equidistributed uniform pseudo-random number generator[END_REF]], a random number of triplets (l v , p v , r v ) such that 1 l v p v r v 6000. Table 9 summarizes the characteristics of the considered datasets. Finally, in order to evaluate the performance of the B&P algorithm also on biological instances of the PLOHP, we considered a biological instance provided by deCode Genetics and hereafter called "Bio22". This in-stance is constituted by 20022 point-intervals on the real segment [1, 6000], each having minimum length 2 and maximum length 6000.

The point-intervals in Bio22 refers to genetic fragments from chromosome 22 extracted from a population of 18360 individuals and constists of 6000 single nucleotide polymorphisms. All the instances used in our experiments can be downloaded at http://homepages. ulb.ac.be/~lporrett/PLOHP/Datasets.tar.bz2.

Computational performances

Table 10 shows the performances of CLHF and the B&P in solving the 5 real instances described in Halldórsson et al. [2011b] with and without the presolving strategies. In particular, the table is constituted by three main sections, each one corresponding to the datasets "gens-

x α ", "gens-x β " and "gens-x γ ", x ∈ {1, 5}, previously described. Moreover, the columns are divided in four sections, each of them showing the solution time (expressed in seconds) and the number of cliques needed to optimally solve a given instance of the PLOHP. Specifically, the first pair of columns is related to the execution of CLHF. The second pair of columns is related to the execution of the B&P. The third pair of columns is related to the execution of CLHF using the presolving strategies. Finally, the fourth pair of columns are related to the execution of the B&P using the presolving strategies. As a general trend, Table 10 shows that, independently of the presence (or absence) of the presolving strategies, the B&P is always faster than CLHF. This fact alone justifies and shows the importance of using inverse projection when dealing with very large instances of the PLOHP. Moreover, the table also shows the benefits of using the presolving strategies: the removal of over 70% of cloned intervals allows the removal of over 70% of x v variables; the search for the connected components enables the decomposition of the corresponding induced MPTG into at least 29 connected components (independent subproblems), which leads to a vast speed-up of the solution process. Finally, the table shows that, although the B&P is faster than CLHF, it proves unable to solve the instances "gens-x γ ", x ∈ {1, 5}, within the considered time limit.

This fact is mainly due to the pricing oracle, which proves particularly time-consuming in gens-x γ . percentage), i.e., the difference between the optimal value to a given instance of the PLOHP in a specific dataset and the objective function value of the linear programming relaxation at the root node of the respective search tree, divided by the optimal value; the nodes explored in the search tree respectively.

Computational performances on larger datasets of the PLOHP

In general, Table 11 shows that CLHF and the B&P, combined with presolving strategies, is able to solve each simulated instance in about half an hour, except for the parameter set "Sim-5" that required on average 969.5657 seconds for the B&P while CLHF exceeds the 64GB of memory available for the experiments.

Tables 12 report on the minimum, maximum and average number of connected components,intervals, singletons, edges and Cloned nodes, while Table 13 report on the minimum, maximum and average time (expressed in seconds) to filter the singletons and the cloned intervals, and to identify the connected components.

Table 14 shows the performance of the B&P in solving "Bio22" instance using the presolving strategies. The table reports the solution time (expressed in seconds) necessary to solve the biological instance of the PLOHP; the gap (expressed in percentage), i.e., the difference between the optimal value and the objective function value of the linear programming relaxation at the root node of the search tree, divided by the optimal value; the nodes explored in the search tree; and the overall number of connected components, intervals, cloned intervals and singletons, respectively. In general, Table 14 shows that the B&P, combined with presolving strategies, is able to solve the biological instance in 664.40 seconds.

Interestingly, Table 11 shows that the biological instance differs from the instances in Table 10 by having a fewer number of point-intervals and a relatively smaller number of cloned ones. In this case, filtering the cloned point-intervals reduces the number of x v variables only by the 20%, instead of the reduction of variables x v by 70% shown in Table 10. However, the relatively higher number of connected components in Bio 22, allows to decompose the problem into a bigger family of small independent subproblems, providing a remarkable speed up in the solution time.

In the next chapter we study a problem arising while analysing the whole genome re-sequencing. Instead of having only the information about the SNPs genotype, the whole genome data allows the detection of all kinds of genomic variation and not only the ones highlighted by the SNPs. In particular, we examine the polymorphism present in the human genome that is caused by the insertion/deletion of a specific family of copy-number variations.

A In this chapter we analyze a particular form of SINE variation called

Arthrobacter Luteus (ALU) insertions.

An ALU sequence is an approximately 300 basepair long sequence derived from 7SL RNA gene [START_REF] Ullu | Alu sequences are processed 7SL RNA genes[END_REF]. Although ALU elements do not encode genes, many studies suggest their functional importance. ALU elements are recognized to affect protein synthesis at the transcriptional and post-transcriptional level [START_REF] Sorek | Alu-Containing Exons are Alternatively Spliced[END_REF][START_REF] Kelley | Transposable elements modulate human RNA abundance and splicing via specific RNA-protein interactions[END_REF] as well as DNA methylation [START_REF] Andrade | Genetic and epigenetic variations contributed by Alu retrotransposition[END_REF]] and other cellular processes [START_REF] Prescott | Alu repeats and human disease[END_REF]. Furthermore, they are thought to be major drivers of genome evolution [START_REF] Hormozdiari | Rates and patterns of great ape retrotransposition[END_REF][START_REF] Salem | Alu elements and hominid phylogenetics[END_REF]] and assist in the creation of structural variation [START_REF] Wang | Whole genome computational comparative genomics: A fruitful approach for ascertaining Alu insertion polymorphisms[END_REF]]. The importance of ALU elements is further highlighted by the potential association with genetic instability, one of the principal cause factors in many disorders including cancer [START_REF] Prescott | Alu repeats and human disease[END_REF][START_REF] Zhang | Alu distribution and mutation types of cancer genes[END_REF][START_REF] Helman | Somatic retrotransposition in human cancer revealed by whole-genome and exome sequencing[END_REF].

ALU repeats are SINE that occur frequently in the human genome, as well as in other genomes. The ALU sequence family has been propagated to more than one million copies in primate genomes over the last 65 million years. ALU repeats are the largest family of mobile elements in the human genome and the ALU family comprises more than 10% of the human genome [START_REF] Cordaux | The impact of retrotransposons on human genome evolution[END_REF]]. Most ALU repeats were inserted early in primate evolution, where it is estimated that there was approximately one new ALU insertion in every primate birth [START_REF] Mark | Alu repeats and human genomic diversity[END_REF].

The current rate of ALU insertion is estimated to be of the order of one ALU insertion in every 200 births [START_REF] Prescott | Alu repeats and human disease[END_REF]. Some members of these young ALU subfamilies have been inserted into the human genome so recently that they are polymorphic with respect to the presence or absence of insertion in different human genomes.

Those relatively few elements that are present in the genomes of some individuals and absent from others are referred to as ALU-insertion polymorphisms.

In this chapter we consider the problem of detecting polymorphic ALU insertions from DNA sequence reads using high throughput pairedend sequencing data [START_REF] Fullwood | Next-generation DNA sequencing of paired-end tags (PET) for transcriptome and genome analyses[END_REF].

The problem was presented for the first time in [START_REF] Sveinbjörnsson | PAIR: polymorphic Alu insertion recognition[END_REF], where the authors provided a polynomial algorithm able to detect ALU insertions coming from a single individual. Simultaneously, they introduce the problem for detecting inserted ALUs over multiple individuals providing an heuristic approach but not an efficient and exact solution method. Specifically, we investigate in detail the problem arising in the multiple individuals scenario by analyzing the combinatorial structure of the problem and presenting a more

notation and problem statement

In this section, we introduce notation that will prove useful to state the problem of detecting ALU insertions. In particular, we first introduce a single individual scenario, i.e, insertions detected from DNA sequence reads coming from a single individual. Then, we extend the concepts to the general case in which multiple individuals share the same insertion.

To this purpose, consider an individual as a set P of read pairs.

A read pair consists of a read of a fixed length L, followed by a short spacing Y, followed by another read of the same length L. In rection ensures that: if the location of one of the reads is known then the location of the mate (the other read) is also known, up to Y (see Figure 12a). The spacing between the two reads can be assumed to be known a priori or to be easily estimated from the sequence reads [START_REF] Korbel | PEMer: a computational framework with simulation-based error models for inferring genomic structural variants from massive paired-end sequencing data[END_REF]. In particular, let define Y max as the maximal insert length between two reads.

Using the description provided above, we can define a read pair π as a pair of left read π l and right read π r , which are mates to each other, denoted as mate(π l ) = π r and mate(π r ) = π l . We use π n to denote either a left or a right read when the relative position in the pair is not relevant. Since the genome sequence of the individual(s) being sequenced is however not known a priori, but is highly similar to the reference genome, reads are mapped to the reference to infer the genome of the individual. In particular, if π n is mapped to a reference genome, we use begin(π n ) and end(π n ) to represent its start and end position in the mapped reference genome. We define the insert length of a read pair π, measured with respect to the reference genome, as ∆(π) = end(π r )begin(π l ). At some locations in the reference genome the genomes of the reference and the individual(s) being sequenced will diverge. Some of this divergence is due to the insertion of ALU polymorphisms. In particular, whenever Specifically, δ(π l ) = {p : begin(π l ) p end(π l ) + Y} and δ(π r ) = {p : begin(π r ) -Y p end(π r )} identify sets of locations in the reference genome candidate to contain an ALU insertion for a left and right read, respectively. By considering two reads π n i and π n j coming from two distinct improper read pairs i and j, we say that π n i and π n j are compatible with an ALU insertion whenever δ(π n i ) ∩ δ(π n j ) = ∅ (see Figure 12b). We can now define a Region Compatible with Insertions (RCI) as follows:

∆(π) > 2 × L + Y max ,
Definition 1. Given the set of location of a reference genome G and P ⊆ P a set of improper read pairs. Le be Π = δ(π l i ) ∪ δ(π r i ) such that π l i , π r i ∈ P, then a Region Compatible with Insertions is any subset X ⊆ G such that X = j∈Π δ(π n j ) and X = ∅.

Since reads mapping to ALU sequences will almost always have multiple places on the genome that have similar quality mapping, each read either gives evidence of an ALU insertion or to reads improperly mapped or read. In order to detect actual ALU insertions from false positives, denote R i and E i as the set of RCIs and the set of locations leading to errors induced by individual i respectively, and h(ρ) and g(η) represent the costs of detecting an RCI ρ ∈ R i and a sequencing error η ∈ E i , respectively. In particular, we can define the following optimization problem:

Problem 3 (Polymorphic ALU Insertions Recognition Problem (PAIRP)).

Given Pi a set of improper read pairs associated to individual i, minimize the overall cost

Z = ρ∈R i h(ρ) + η∈E i g(η)
such that each ρ ∈ P is either compatible with an RCI or is classified as an error.

A possible extension of the PAIRP to a multiple individuals scenario MIPAIRP, begins by considering that two read pairs π i ∈ P a and π j ∈ P b , coming from two different individuals a and b, provide locations compatible with insertion whenever they are mutually compatible.

In particular, let denote R = i=1,...,n R i and E = i=1,...,n E i as the set of RCIs and the set of locations leading to errors induced by n individuals respectively. We can define the following optimization problem:

Problem 4 (Multiple Individuals Polymorphic ALU Insertions Recognition Problem (MIPAIRP)). Given P a set of improper read pairs, minimize the overall cost

Z = ρ∈R h(ρ) + η∈E g(η)
such that each ρ ∈ P is either compatible with an RCI or is classified as an error.

The PAIRP and MIPAIRP are based on the parsimony principle [START_REF] Victor | Parsimony, Phylogeny, and Genomics[END_REF]]. This fact implies that the optimal solutions to the problem provide estimations of insertion events that, in the worst case, are lower bounds on the overall number of true insertion events occurred in the set of individuals being considered [START_REF] Catanzaro | The minimum evolution problem: Overview and classification[END_REF][START_REF] Aringhieri | Optimal solutions for the balanced minimum evolution problem[END_REF]. [START_REF] Sveinbjörnsson | PAIR: polymorphic Alu insertion recognition[END_REF] provided a polynomial algorithm for the PAIRP and proposed a heuristic approach for the MIPAIRP without investigating the complexity of the latter any further. In the next sections we shall address this major issue and provide an algorithm able to exactly solve instances of both general problems.

the alu graph

In this section we define a new class of graphs, which we term ALU Graph (ALUG) and will turn out to be useful in transforming the PAIRP, and by extension the MIPAIRP, into a particular version of the Minimum Clique Partition Problem (MCPP) [START_REF] Michael | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF].

In order to characterize such a class of graphs, consider a set of improper read pairs P and denote between two distinct intervals exist, i.e., an edge is related to the presence/absence of mutual compatibility between two distinct reads, then G is a variant of interval graph (see [START_REF] Fishburn | Interval Orders and Interval Graphs:A Study of Partially Ordered Sets[END_REF]), more simply, an ALU Graph (ALUG). The class of the ALUGs can be seen as a specialization of the class of the Interval Graphs (IGs). In fact, the following 13a). Hence, any ALUG is also an IG. To complete the proof, it is sufficient to show that the converse is not true. In fact, the class of ALUGs is characterized by always having an even number of intervals, which in turn implies that an ALUG does not contain graphs with an odd cardinality of vertices. In contrast, the class of the IGs may also include such graphs. 13 shows the mapping of improper reads to a segment of the reference genome, while Figure 13b shows the clique associated to the RCI represented in Figure 13a Given an instance of the PAIRP and its corresponding ALUG, G, we observe that by definition, the RCIs correspond to cliques in G, i.e., to complete subgraphs of G. This is because any RCIs represent a set of improper reads all mutually compatible with ALU insertions. We also recall that a maximal clique of a graph is a clique that is not a subset of a larger clique and a maximum clique is a clique of maximum size.

I l k = [begin(π l k ), end(π l k ) + Y max ] and I r k = [begin(π r k ) -Y max , end(π r k )]
[l k , r k ] of G into the interval [begin(π l k ), end(π l k ) + Y max ] or [begin(π r k ) -Y max , end(π r k )] of A (see e.g., Figure
Then, the following result holds for ALUGs: Proposition 5. Let G be an ALUG. Then G contains at most v maximal cliques.

Proof. The statement follows from observing that the class of the ALUG is contained in the class of interval graphs and that the number of maximal cliques in a interval graph is proven to be at most equal to the number of the vertices (see [START_REF] Fulkerson | Incidence matrices and interval graphs[END_REF])

The procedure is described in [START_REF] Fulkerson | Incidence matrices and interval graphs[END_REF] and consists in deleting a random vertex on the graph and simply list it together with its neighbours in the reduced graph. If the original graph has n vertices, this yields a set of n cliques. The maximal ones of these n cliques are the maximal ones in the graph and there can be at most n. This implies that the maximum clique problem can be solved in polynomial time in an ALUG.

A representation of MIPAIRP by ALUGs consist of considering multiple set of intervals I, each of them coming from a different individual.

Hence, let us define the set I as the union of all I i where i represent the individual and G(I) is an ALUG for the MIPAIRP. Moreover, it is trivial to see that RCI corresponds to cliques in G(I) and Proposition 5 is valid also for G(I).

In the next section we exploit Proposition 5 to develop an exact approach for the solution of the MIPAIRP based on integer programming.

the complexity of pairp and mipairp

Given an instance of the PAIRP/MIPAIRP and its corresponding ALUG, G, we note that in any optimal solution to the instance, a false positive read will always be a read that does not belong to any RCI selected.

Hence, in any optimal solution to the instance, a false positive read corresponds to a clique of the ALUG having cardinality '1'. This insight allows us to consider the PAIRP and the MIPAIRP as instance of the MCPP in a particular class of graphs [START_REF] Michael | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]]. The MCPP is known to be NP-hard in general [START_REF] Michael | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]; in this section we will show that the MCPP (i) remains hard even when restricted to the class of the ALUG, (ii) can be solved in polynomial time for the PAIRP version if the cost functions h and g satisfy some specific properties while (iii) the MIPAIRP can not be solved optimally in polynomial time unless P=NP.

Before proceeding, we introduce some notation that proves useful throughout the section. We say that a set function f is zero-cardinal if f(∅) = 0; non-negative if f assumes only non-negative values; and non-decreasing if f(T ) f(S) for any T ⊆ S ⊆ V. We say that f is sub-modular if it satisfies the following property [START_REF] Schrijver | Combinatorial optimization[END_REF]]:

f(S ∪ {u}) + f(T ) f(S) + f(T ∪ {u}) ∀ T ⊆ S ⊆ V, u ∈ V \ S
We say that f is a polymatroid rank function if it is zero-cardinal, nondecreasing, non-negative, and sub-modular. Moreover, similarly to [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF], we define a value-polymatroid set function f as a zero-cardinal, non-decreasing, non-negative set function that satisfies the following property:

f(S ∪ {u}) + f(T ) f(S) + f(T ∪ {u}) ∀ S, T ⊆ V : f(S) f(T ), u ∈ V \ (S ∪ T ).
Note that a value-polymatroid set function is also polymatroidal, but the converse is generally not true [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF]]. Finally, a set function f is size-defined sub-modular if there exists a function

Φ : [0 . . . |V|] → R + 0 such that f(S) = Φ(|S|)
, for any S ⊆ V. As shown in [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF]], a size-defined sub-modular set function f is both value-polymatroidal and polymatroidal. The fact that ALUGs are variants of the class of interval graphs and the previous definitions turn out to be useful to investigate the complexity of the PAIRP. Specifically, denote C(G) the set of cliques of G and set:

f(C) =              0 if C = ∅ g(C) if |C| = 1 h(C) if |C| 2 ∀C ∈ C(G) (16
l k ) = l k , end(π l k ) = r k , begin(π r k ) = l k + 3 * , end(π r k ) = r k + 3 * and Y max = λ.
The ALUG G induced by P has I(G ) = {I l ∪ I r } as the set of all the "left" and "right" intervals. By construction it is easy to se that G is the disjoint union of the graph induced by I l and I r and that G is isomorphic to both of them. Hence, solving the PAIRP on G has at least the same complexity of solving the MCPP on an interval graph.

Proposition 7. The decision version of the PAIRP is NP-complete even when the cost function f is restricted to polymatroidal set functions.

Proof. The statement follows Proposition 6 and that the MCPP on an interval graph is NP-complete when the cost function is polymatroidal [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF] In general, it is easy to realize that the decision version of the PAIRP is NP-complete for any cost function f(C) = ψ(|C|)σ(C) such that ψ : [0...|V|] → R + 0 and σ(C) is a generic function on C. In fact, such a case also includes the rooted-TSP cost function on a tree (see [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF]) which is trivially polymatroidal. Although Proposition 7 states that the decision version of the PAIRP is in general NP-complete, it is worth noting that in some special cases the problem can be still solved in polynomial time. For example, the following proposition holds:

Proposition 8. Let G = (V, E) be an ALUG and f : C(G) → R + 0 a valuepolymatroidal cost function. If G is also an interval graph then it is possible to compute a minimum cost partition into cliques of G in polynomial time.

Proof. The statement follows from the fact that the minimum clique partition problem on an interval graph can be solved in polynomial time when the cost function is a value-polymatroidal set function [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF].

Proposition 8 turns out useful to show that if G is an interval graph the PAIRP can be solved in polynomial-time when the following objective function is used:

f α (C) =              0 if C = ∅ c1 if |C| = 1 c2 if |C| 2 ∀C ∈ C(G) (17)
where c 1 and c 2 are two constants such that 0 that C ∪ {t} is a clique in G. Note that t belongs to some T ∈ P \ C and that since f α is non-decreasing, it holds that f α (C) f α (T ) f(T \ {t}). Observe also that since f α is α-value-polymatroidal, it holds that f α (C ∪ {t}) + f α (T \ {t}) f α (C) + f α (T ). Hence, it is possible to enlarge C until it becomes a maximal clique of G without getting worse than the cost of P .

< c 1 αc 1 < c 2 (α + 1)c 1 ,
({v}) = |C|c 1 < c 2 = f α (C) (α + 1)c 1 .
Based on Propositions 7,8 and 9 a recursive dynamic programming approach for the PAIRP is provided in [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF].

On the other hand, in the multiple individual scenario, the presence Based on Propositions 10 and cost function ( 18) the recursive dynamic programming approach provided in [START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF] can not be used for the MIPAIRP.

an integer programming model for the mipairp

The NP-hardness of the MIPAIRP justifies the development of exact and approximate solution approaches for the problem. In this section we present an integer programming model for the MIPAIRP. To this end, C∈ Ĉ(G)

y C + x v 1 ∀v ∈ V (21b) x v ∈ {0, 1} ∀v ∈ V (21c) y C ∈ {0, 1} ∀C ∈ Ĉ(G) (21d)
Formulation 5 has the benefit of being polynomial-sized, due to Proposition 5, hence in principle its relaxation does not require the use of column generation techniques to be solved.

Formulations 4 and 5 can be strengthened by adapting appropriately the inequalities described in [START_REF] Grötschel | Facets of the clique partitioning polytope[END_REF], [START_REF] Bandelt | Lifting theorems and facet characterization for a class of clique partitioning inequalities[END_REF] and [START_REF] Oosten | The clique partitioning problem: Facets and patching facets[END_REF]. Moreover, additional valid inequalities can be considered. Specifically, given a pair of distinct vertices v, w ∈ V, we say that v dominates w if Proposition 11. Let v be a dominating vertex in V. Then, the inequality

C v ∈ Ĉ(G)
y C 1 is valid for both Formulations 4 and 5.

experimental results

In this section we analyze the performance of our model in solving instances of the Multiple Individuals Polymorphic ALU Insertions Recognition Problem (MIPAIRP). Our experiments were motivated mainly to measure the runtime performances of our model in tackling real size instances of the MIPAIRP. We emphasize that our experiments simply aim to evaluate the computational performance of our model; we neither attempt to study the efficiency of it to predict ALU insertions across the genomes of a population of individuals nor to compare its accuracy versus ALU insertions predictors using an objective function that is different from the one used in [START_REF] Sveinbjörnsson | PAIR: polymorphic Alu insertion recognition[END_REF].

length 120 and insertion distance Y max fixed to 200. Hence, we also considered increasing size instances of the MIPAIRP obtained by augmenting the number of individuals from each "set-x", x ∈ {1, 3}, by α = 1 000, β = 2 000 and γ = 3 000, respectively. We generate for each "set-x", x ∈ {1, 3}, 15 random instances of the MIPAIRP for a total of 135, that are hereafter grouped into three datasets of 45 denoted as "Set-x α ", "Set-x β " and "Set-x γ ", x ∈ {1, 3}, respectively. Table 15 shows the main characteristics of the 135 instances. In particular, the table is composed by three main sections, each one corresponding to the datasets "Set-x α ", "Set-x β " and "Set-x γ ", x ∈ {1, 3}, previously described.

The first column of Table 15 reports the name of the datasets. The second column reports the number of reads generated for each individual. The third column reports the overall number of individuals analyzed in each given instance. Finally, the fourth column reports the average number of maximal cliques contained in all the 15 instances.

Computational performances

Table 16 shows the performances of Formultaion 5 in solving the 135 instances described is Section 4.5.2. In particular, the table is constituted by three main sections, each one corresponding to the datasets "Set-x α ", "Set-x β " and "Set-x γ ", x ∈ {1, 3}, previously described. The first column of Table 16 reports

perspectives

The primary goals of this thesis have been the study of problems arising in the field of biological systems under the OR point of view. Our Finally, even if all the problems faced in this thesis are NP-hard, we showed how to make some progress and be able to efficiently solve large scale instances by providing compact formulations and efficient decomposition methods. Moreover, we have made improvements in being able to efficiently solve large scale instances for well known problems in the domain. It is our belief that the proposed methods and solving techniques will be helpful in analysing real-size datasets to identify new variation sites associated with genetic disease, and the basis for the development of new mathematical models for the arising problems in the GWAS domain.

  problems in the domain of Genome-Wide Association Studies: (i) the Pure Parsimony Haplotyping under Uncertain Data Problem that consists in finding the minimum number of haplotypes necessary to explain a given set of genotypes containing possible reading errors; (ii) the Parsimonious Loss of Heterozygosity Problem that consists of partitioning suspected polymorphisms from a set of individuals into a minimum number of deletion areas; (iii) and the Multiple Individuals Polymorphic ALU Insertions Recognition Problem that consists of finding the set of locations in the genome where ALU sequences are inserted in some individual(s). All three problems are NP-hard combinatorial optimization problems. Therefore, we analyse their combinatorial structure and we propose an exact approach to solution for each of them. The proposed models are efficient, accurate, compact, polynomial-sized and usable in all those cases for which the parsimony criterion is well suited for estimation. vii R É S U M É Le domaine interdisciplinaire de la biologie des systèmes a évolué rapidement au cours des dernières années. Différentes disciplines ont contribué au développement de la branche expérimentale aussi bien que de la branche théorique Bien que la biologie computationnelle a été une activité en en croissance en informatique depuis plusde deux décennies, ce n'est que depuis quelques années que des modèles d'optimisation ont été de plus en plus développés et analysés par des chercheurs dont la spécialité de base est la recherche opérationnelle. Cette thèse vise à apporter une contribution au domaine de la biologie computationnelle en appliquant la programmation mathématique à certains problèmes de biologie moléculaire. Plus précisément, nous abordons trois problèmes dans le domaine de Genome-Wide Association Studies: (i) le problème appelé Pure Parsimony Haplotyping under Uncertain Data, qui consiste à trouver le nombre minimum d'haplotypes nécessaire pour expliquer un ensemble donné de génotypes contenant des erreurs de lecture potentielles; (ii) le problème appelé Parsimonious Loss of Heterozygosity Problem, qui consiste dans le partitionnement des polymorphismes soupçonnés à partir d'un ensemble d'individus en un nombre minimal de zones de suppression; (iii) et le troisième problème appelé Multiple Individuals Polymorphic ALU Insertions Recognition Problem, qui consiste à trouver l'ensemble des emplacements dans le génome où une séquence d'ALU est insérée dans certains individus. Les trois problèmes sont des problèmes d'optimisation combinatoire NP-difficile. Par conséquent, nous avons analysé leur structure combinatoire et proposé une approche exacte de résolution pour chacun d'entre eux. Les modèles proposés sont efficaces, précis, compacts, de taille polynomiale, et utilisables dans tous les cas pour lesquels le critère de parcimonie est bien adapté à l'estimation. x P U B L I C AT I O N S Some chapters of this thesis are partially based on articles that have been published or submitted for publication. Specifically: Catanzaro, D., M. Labbé, and L. Porretta. 2011. "A Class Representative Model for Pure Parsimony Haplotyping under Uncertain Data." PLOS ONE 6, no. 3 (). Porretta, L., D. Catanzaro, B. Fortz, and B. V. Haldòrsson. 2016. "A Branch&Price Algorithm for the Minimum Cost Clique Cover Problem in Max-Point Tolerance Graphs." Porretta, L., B. Fortz, and B. V. Haldòrsson. 2017. "An Integer Programming formulation of the Multiple Individual Polymorphic ALU Insertion Recognition Problem." "The isolated man does not develop any intellectual power. It is necessary for him to be immersed in an environment of other men, whose techniques he absorbs during the first twenty years of his life. He may then perhaps do a little research of his own and make a very few discoveries which are passed on to other men. From this point of view the search for new techniques must be regarded as carried out by the human community as a whole, rather than by individuals" -Alan Turing A C K N O W L E D G E M E N T S I want to convey my sincere gratitude to Professors Bernard Fortz, supervisor of this work, for offering me the possibility to write my thesis in G.O.M., for giving me the freedom to follow my research interests and finally for making sure that I could spend some of my time abroad.

  Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Genome-Wide Association Studies . . . . . . . . . . . . 1.3 Combinatorial Optimization . . . . . . . . . . . . . . . . 1.4 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . 2 a class representative model for pph-ud 2.1 Notation And Problem Statement . . . . . . . . . . . . . 2.2 Design and Implementation . . . . . . . . . . . . . . . . 2.2.1 Reducing model size . . . . . . . . . . . . . . . . 2.2.2 Accounting for perfect phylogeny constraints . 2.3 Experimental results . . . . . . . . . . . . . . . . . . . . 2.3.1 Implementation . . . . . . . . . . . . . . . . . . . 2.3.2 Benchmark instances . . . . . . . . . . . . . . . . 2.3.3 Computational performances . . . . . . . . . . . 3 a branch-and-price algorithm for the plohp 3.1 Notation and Problem Formulation . . . . . . . . . . . 3.2 A Branch-and-Price algorithm for the PLOHP . . . . . 3.2.1 An integer linear programming formulation for the Parsimonious Loss of Heterozygosity Problem (PLOHP) . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 Figure 2

 12 Figure 1The polyhedral representation of a Mixed Integer Linear Programming (MILP) problem . . . . . 7 Figure 2 Any two copies of the human genome differ from one another by approximately 0.1% of nucleotide sites. In this example, most of the DNA sequence is identical in these chromosomes, but there are three nucleotides where variation occurs. A pattern of DNA sequence variation defines a haplotype. . . . . . . . . . . . . . . . 17

Figure 3 Figure 4 Figure 5 Figure 6

 3456 Figure 3 Graphical representation of an instance of Pure Parsimony Haplotyping (PPH) and two alternative solutions . . . . . . . . . . . . . . . . . . . . 21 Figure 4 Bipartite graph representation of Solutions 1 and 2 of Figure 3. . . . . . . . . . . . . . . . . . 23 Figure 5 Examples of subsets of genotypes induced by the two alternative solutions of Figure 3. . . . . 24

Figure 8

 8 Figure 8An example of point-intervals . . . . . . . . . .

Figure 9

 9 Figure 9 The MPTG induced by the point-intervals shown in Figure 8. . . . . . . . . . . . . . . . . . . . . .

Figure 10

 10 Figure 10 An example of block decomposition of the pointintervals shown in Figure 8. . . . . . . . . . . .

Figure 11

 11 Figure 11 BoxPlot of the solution time taken by the Branchand-Price (B&P) algorithm to exactly solve the learning set of instances of the PLOHP when considering different percentages of violated cliques. . . . . . . . . . . . . . . . . . . . . . . .

Figure 12

 12 Figure 12 Example of improper read pairs aligned to a reference. Arrows show read directions. The blue part of the reads can be mapped to the reference and the red parts are clipped or mapped somewhere else in the reference. The grey part represents the read pair distance Y and the dark grey part represents the area where the read pair distance Y of the two reads overlaps. The interval of possible Arthrobacter Luteus (ALU) locations is highlighted in red over the reference. Note that read pair π 5 and π 6 contribute with two improper reads, while π 7 and π 8 just with one. . . . . . . . . . . . . . . . . . . . . . .

  FigureExample of ALUG for Multiple Individuals Polymorphic ALU Insertions Recognition Problem (MIPAIRP).

Figure

  Figure 15a shows the mapping of improper reads coming from multiple individuals, identified by different colours, to a segment of the reference genome, while Figure 14b shows the cliques associated to the RCIs represented in Figure 14a . . . . . . . . . . . . . . . . . . . . . 89 Figure Example of ALUG for MIPAIRP. Figure 15a shows the mapping of improper reads coming from multiple individuals to a segment of the reference genome, while Figure 15b shows the cliques associated to the RCIs represented in Figure 15a . . . . . . . . . . . . . . . . . . . . . 90 Figure Example of a graph where function(18) prove to be not value-polymatroid. This figure shows a graph composed of vertices belonging to two different individuals X and Y. In particular X = {a, b, c} and Y = {d, e, f, g}. . . . . . . . . . . . . 92

  are vectors c T = (c 1 , . . . , c n ) and h T = (h 1 , . . . , h p ), an m × n matrix A = (a ij ), an m × p matrix G = (g ij ) and a vector b = (b 1 , . . . , b m ) T . The entries of c, h, A, G, b are considered as rational. The vectors x = (x 1 , . . . , x n ) T and y = (y 1 , . . . , y p ) T contain the decision variables of the problem. The elements of x are non-negative integers while the elements of y are non-negative real-valued numbers.

Figure 1 :

 1 Figure 1: The polyhedral representation of a MILP problem

Figure 2 :

 2 Figure 2: Any two copies of the human genome differ from one another by approximately 0.1% of nucleotide sites. In this example, most of the DNA sequence is identical in these chromosomes, but there are three nucleotides where variation occurs. A pattern of DNA sequence variation defines a haplotype.

Figure 3 :

 3 Figure 3: Graphical representation of an instance of PPH and two alternative solutions

  . Haplotyping a set of genotypes G means finding the set of haplotypes resolving G.

Figure 3

 3 Figure3are depicted in Figure4aand 4b, respectively. The bipartite graph representation of a solution suggests that in a feasible solution to PPH the haplotypes induce a family of subsets of genotypes satisfying the following three properties: (i) each subset of genotypes shares one haplotype, (ii) each genotype belongs to exactly two subsets, and

Figure 4 :

 4 Figure 4: Bipartite graph representation of Solutions 1 and 2 of Figure 3.

Figure 5 :

 5 Figure 5: Examples of subsets of genotypes induced by the two alternative solutions of Figure 3.

  [2006]'s datasets for testing the performances of our model. Specifically, through Hudson's MS program[START_REF] Hudson | Generating samples under a Wright-Fisher neutral model of genetic variation[END_REF]], Brown and Harrower created two families of datasets (called the uniform and nonuniform datasets) by randomly pairing the resulting haplotypes. The distinction in the two simulated methods comes in how the random pairing is performed. In the uniform datasets the haplotypes are randomly paired by sampling uniformly from the set of distinct haplotypes. In the nonuniform datasets the haplotypes are sampled uniformly from the collection of haplotypes generated by the coalescent process. In this collection, haplotypes may not be unique, so some haplotypes are sampled with higher frequency than others. Both the uniform and non-uniform datasets consist of collections of 30 or 50 genotypes having 10, 30, 50, 75 or 100 SNPs each. Each dataset contains a number of instances variable between 15 and 50. The authors also considered biological data from chromosomes 10 and 21, over all four Hap-Map[START_REF] Consortium | The international hapmap project[END_REF]] populations. For each input the authors selected sequences having 30, 50, and 75 SNPs, respectively, giving a total of 8 datasets consisting of 3 instances each. Brown and Harrower's datasets are not subjected to uncertainty, for this reason we considered a three sets of random generated error mask matrices having an error ratio (i.e., the number of entries equal to 1 divided by mp) equal to 5%, 10%, and 15% respectively.Brown et al. [2006]'s datasets and the error mask matrices used in our experiments can be downloaded at the address: http://homepages.ulb.ac.be/~lporrett/PPHUD/PPHerr. zip.

Figure 6 :

 6 Figure 6: In this example, we compare the DNA sequences from a set of individuals correlated between each other (father, mother and child).

  such a criterion in terms of a combinatorial optimization problem, called the Parsimonious Loss of Heterozygosity Problem (PLOHP), consisting of finding a minimum clique covering on a particular class of interval graphs called Max-Point Tolerance Graph (MPTG) (see Section 3.1). The authors proved the NP-hardness of the PLOHP, investigated some fundamental properties of MPTGs, and proposed an Integer Linear Programming (ILP) formulation able to optimally solve instances of the PLOHP containing up to 2000 individuals and 2000 SNPs within 5 hours of computing time. Unfortunately, the proposed formulation proved unable to optimally solve instances of the PLOHP involving a larger number of individuals, a desirable requirement in disease association studies.Starting from[START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]' seminal work, in this chapter we present an ILP formulation for the PLOHP based on column generation. We introduce a number of preprocessing techniques to reduce the size of a given instance of the problem and we investigate possible decomposition strategies to divide a reduced instance into a family of independent subproblems having even smaller size. We develop a new efficient algorithm to find maximum node-weighted cliques in MPTGs and we embody this algorithm in a Branch-and-Price (B&P) algorithm for the PLOHP, that turned out to be 10-30x faster than the previous approach described in the literature. The new algorithm is freely available and enables the solution of very large practical instances of the PLOHP, containing over 6000 trios and SNPs, within 1

Figure 7 :

 7 Figure 7: An example of a trio, i. e., a set of two parents and their offspring, and their genotypes. SNPs are highlighted in red. Specifically, the first highlighted column represents a SNP inconsistent with a loss of heterozygosity; the second highlighted column represents a SNP consistent with a loss of heterozygosity; the third highlighted column represents a SNP showing an evidence of a loss of heterozygosity.

Figure 8 :Figure 9 :

 89 Figure 8: An example of point-intervals

  Consider a set I of point-intervals and let G = (V, E) be the MPTG induced by I. Then, the following proposition holds [Catanzaro et al. 2013]: Proposition 3. A MPTG contains at most |I|(|I| -1)/2 maximal cliques.

Formulation 2 .

 2 a singleton in a solution to the problem and 0 otherwise. Similarly, let y C be a decision variable equal to 1 if the maximal clique C ∈ C(G), with |C| 2, is selected in a solution to the problem and 0 otherwise. Then, a valid ILP formulation for the PLOHP is: Catanzaro-Labbé-Halldórsson ILP Formulation (CLHF)

  objective function (12a) minimizes the cost of finding a clique cover of G, and constraints (12b) impose that, for each vertex v ∈ V, the solution contains either a clique or a singleton in C(G) that covers v. Note that, as the overall number of maximal cliques in a MPTG is polynomially bounded, CLHF is polynomial-sized. However, when the instances of the PLOHP are very large (e.g., they include more than 250,000 point-intervals), the use of CLHF may become unpractical due to the large number of y variables. A possible approach to overcome such a limitation consists of reformulating CLHF by inverse projection[START_REF] Martin | Large Scale Linear and Integer Optimization: A Unified Approach[END_REF]]. To this end, denote C τ (G) as a subset of C(G) that includes τ < |C(G)| maximal cliques in G and set C τ v (G) = {C ∈ C τ (G) : v ∈ C}. Now, consider the following Reduced Master Problem (RMP) associated with the linear programming relaxation of CLHF: Formulation 3. CLHF -The Reduced Master Problem (RMP)

Specifically, Algorithm 1

 1 takes G as an input and iterates the following steps: line 2 solves the RMP; line 3 gets the dual values {µ} and

  8: until Ĉ = ∅ 9: return z * RMP associates a weight µ v to each vertex v ∈ V; line 5 calls the oracle and searches in C(G) for a maximum node-weighted clique Ĉ satisfying (15); if such a clique exists, line 7 adds variable y Ĉ (and its corresponding column) to the RMP and iterates lines 1-8; if, there are NO violated constraints (13b), the solution to the RMP is provably optimal and Algorithm 1 returns the optimal value to the RMP. Preliminary experiments showed that to decrease the running time of Algorithm 1 it is advisable to add not just the most violated node-weighted clique at each iteration of the pricing oracle, but possibly many cliques C having an overall weight larger than f(C). Section 3.3 explores, via computational experiments, the problem of deciding how many cliques should be added in a generic instance of the PLOHP.3.2.2 Finding a maximum node-weighted clique in a MPTGA critical step in Algorithm 1 consists of finding a maximum nodeweighted clique in G. Although this problem is NP-hard for general graphs, it can be solved in polynomial time when the graph is MPTG and the objective function is size-defined submodular. In fact, it is Algorithm 2 | Pricing Oracle Input: G, µ Output: C : a maximum node-weighted clique in G 1: set C = ∅, max = 0, get I from G 2: for all I i , I j ∈ I : i < j and adjacent(I i , I j ) do

  , as G is a MPTG, any subset of vertices that forms a maximal clique C in G is characterized by having two vertices whose corresponding points, say p l and p r , are the furthest to the left and to the right, respectively, in the set of point-intervals induced by C.

  iterates the following procedure: for each point-interval I k ∈ I (i) find a maximal clique C that cover (the vertex induced by) I k and (ii) remove from I the point-intervals (whose induced vertices are) in C.The procedure stops when I = ∅. The complexity of Algorithm 3 isO(|I| 3 ). Algorithm 3 | C τ (G)-Generator Input: Set of intervals I Output: C τ (G) 1: C = ∅, Q = I 2:for all I i , I j ∈ I : i < j and adjacent(I i , I j ) and Q = ∅ do

Figure 10 :

 10 Figure 10: An example of block decomposition of the point-intervals shown in Figure 8.

  3.1.0[START_REF] Achterberg | SCIP: Solving constraint integer programs[END_REF]] to handle the column generation and the B&P routines. Moreover, we used FICO Xpress Optimizer as linear programming solver. The experiments have been performed on an Intel Core i7-4930K CPU, 3.40GHz, equipped with 64 GByte RAM and operating system Ubuntu release 12.10 (kernel linux 3.5.0-41generic). During the runtime of CLHF we activated Xpress automatic cuts, Xpress presolving strategy, and Xpress primal heuristic to generate the first upper bound for the PLOHP. Finally, similar toHalldórsson et al. [2011b] and[START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF], we set the maximum runtime for the Formulations to 5 hours. The executable codes used in the experiments can be downloaded together with the instances at url: http://homepages.ulb.ac.be/~lporrett/PLOHP/Code.tar.bz2.

Figure 11 :

 11 Figure 11: BoxPlot of the solution time taken by the B&P algorithm to exactly solve the learning set of instances of the PLOHP when considering different percentages of violated cliques.

  ing 645572, 645552, 645510, 645543 and 659885 point-intervals on the real segment[1, 3575], respectively, each having minimum length 1 and maximum length 9. Catanzaro et al. observed that the ILP formulation presented in[START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] proved unable to solve any of these instances within 1 hour of computing time. Hence, the authors also considered smaller instances of the PLOHP obtained from the given ones by extracting from each "gens-x", x ∈ {1, 5}, the pointintervals contained in the real segments α = [1, 1500], β =[1, 2000]

  N I N T E G E R P R O G R A M M I N G F O R M U L AT I O N O F T H E M U LT I P L E I N D I V I D U A L P O LY M O R P H I C A L U I N S E R T I O N R E C O G N I T I O N P R O B L E MWhole genome re-sequencing using next generation sequencers rapidly become the sledgehammer of genome wide association studies and started to be preferable over SNP genotyping for association studies as it allows for the detection of all genomic variation and not only SNP variation. Even if SNPs remains the most abundant form of variation between two individuals, other forms of variation exist, such as inversions, copy-number variations, Long INterspersed Elements (LINE) and Short INterspersed Elements (SINE) elements.

  general predictive model able to solve it efficiently. In particular we show that the problem of detecting ALU insertion over multiple individuals can be formulated as a specific version of the clique partition problem in a particular class of graphs called interval graphs with particular properties. We prove that this problem is NP-hard in general and we provide a methodology to solve it based on Integer Programming (IP). Specifically, we present an IP formulation and strengthening valid inequalities to reduce the solution space. The work thus gives perspective on the mathematics of detecting ALU insertion in genotype data, provides for the first time methodology suitable for provably optimal solution of hard real instances, and suggests new directions on the development of future efficient exact solution approaches.

Figure 12 :

 12 Figure 12: Example of improper read pairs aligned to a reference. Arrows show read directions. The blue part of the reads can be mapped to the reference and the red parts are clipped or mapped somewhere else in the reference. The grey part represents the read pair distance Y and the dark grey part represents the area where the read pair distance Y of the two reads overlaps. The interval of possible ALU locations is highlighted in red over the reference. Note that read pair π 5 and π 6 contribute with two improper reads, while π 7 and π 8 just with one.

  a read pair is defined as improper, i.e., the ends of the two reads map to locations that are inconsistent with the read pair distance Y max . Improper read pair is one of the main signals indicating the presence of an ALU insertion and each left and right read gives partial information about the location of the ALU insertion[START_REF] Qian | PopAlu: populationscale detection of Alu polymorphisms[END_REF]].

  as the intervals induced by the left and right reads of the k-th read pair, respectively. Moreover, denote I = I l ∪ I r as the set of intervals induced by all the reads in P and set v = |I|. Consider a graph G having a vertex for each interval I l k and I r k , k = 1, . . . , | P|, and an edge between two vertices if an intersection

Proposition 4 .

 4 The class of the IGs strictly contains the class of the ALUGs. Proof. A generic interval [l k , r k ] of an interval graph is completely characterized by the left and right margins l k and r k , respectively. In a similar way, a generic interval of an ALUG A is completely characterized by the interval [L k , R k ] where L k = begin(π l k ) and R k = end(π l k ) + Y max for each left read and L k = begin(π r k ) -Y max and R k = end(π r k ) for each right read, respectively. We can transform any ALUG A into an IG G by mapping each interval

  Figure13: The ALUG is based on the reads showed in Figure12b. In particular, Figure13shows the mapping of improper reads to a segment

Proposition 6 .

 6 Solving the PAIRP on an ALUG has at least the same complexity of solving the MCPP on an interval graph. Proof. Given a generic interval graph G, let define I(G) the set of all the intervals of G. A generic interval [l k , r k ] of an interval graph is completely characterized by the left and right margins l k and r k , respectively. Let define Λ = k∈I(G) {l k , r k } as the union of the extreme points of each interval k. Therefore, let define λ = min x,y∈Λ |x -y| as the smallest space between intervals. Furthermore, let define = max k r k the rightmost extreme point of the rightmost interval. We can now generate the sets of improper read pairs P using the following transformation. For each interval k ∈ I(G) we generate an read pair π k having begin(π

Figure 14 :

 14 Figure 14: Example of ALUG for MIPAIRP. Figure 15a shows the mapping of improper reads coming from multiple individuals, identified by different colours, to a segment of the reference genome, while Figure 14b shows the cliques associated to the RCIs represented in Figure 14a

  of multiple individuals need to be taken in account in the definition of the objective function for the MIPAIRP. For example, by penalising all cases in which an RCI is identified only by improper reads belonging all to different individuals. Since an RCI in an ALUG is represented by a clique, these biased RCI may arise whenever the cardinality of the clique is equal to the number of individuals contributing in the same clique (see Figure14).

Figure 15 :

 15 Figure 15: Example of ALUG for MIPAIRP.Figure15ashows the mapping of improper reads coming from multiple individuals to a segment of the reference genome, while Figure15bshows the cliques associated to the RCIs represented in Figure15a

Figure 16 :

 16 Figure 16: Example of a graph where function(18) prove to be not valuepolymatroid. This figure shows a graph composed of vertices belonging to two different individuals X and Y. In particular X = {a, b, c} and Y = {d, e, f, g}.

  given a vertex v ∈ V, we denote C v = {C ∈ C(G) : v ∈ C}. Moreover, denote y C as a decision variable equal to 1 if the clique C ∈ C(G) is selected in the optimal solution to the problem and 0 otherwise. Then, Formulation 4 is a valid formulation for the MIPAIRP: ) impose that each vertex v ∈ V belongs to a selected clique C ∈ C(G) and constraints (19c) impose the integrality on variables y C . Formulation 4 is characterized by an exponential number of variables and constraints and its linear relaxation can be exactly solved by using column generation techniques. Specifically, observe that a variable with negative reduced cost in the linear relaxation of Formulation 4 corresponds to a dual constraint violated by the current dual solution. Denoted µ v as the dual variables associated with constraints (19b), the dual of the linear relaxation of Formulation 4, denoted by LP1, is characterized by the following constraints:v∈V:v∈C µ v f(C) ∀C ∈ C(G)(20)Constraints (20) are violated if there exists a clique C ∈ C(G) such that v∈V:v∈C µ v > f(C) ∀C ∈ C(G). The existence of such a clique can be checked in polynomial time by using Proposition 5 and this in turn implies that the linear relaxation of LP1 can be solved in polynomial time. Interestingly, if the cost function f is defined as in (18) then Formulation 4 can be rewritten as follows. Denote x v as a decision variable equal to 1 if vertex v ∈ V forms a clique of cardinality 1 in the optimal solution to the problem and 0 otherwise. Moreover, denote Ĉ(G) as the set of all maximal cliques in G having cardinality greater or equal to 2 and Ĉ(G) = {C ∈ C(G) : v ∈ C, |C| 2}.Then, Formulation 5 is a valid formulation for the MIPAIRP.

  (i)N(w) ⊂ N(v) and (ii) for any C v ∈ Ĉ(G) and C w ∈ Ĉ(G) it holds that |C v | > |C w |. In such a case, we also say that v is dominating and w is dominated. Dominated vertices are irrelevant to the clique partitioning of G, since in any optimal solution to the problem they will always be identified as cliques of cardinality one. Hence, both formulations can be strengthened by adding the following valid inequalities, whose proof trivially follows from Proposition 8 and the definition of domination:

  the name of the datasets. The second column reports the average number of constraints generated for each instance. The third column reports the average number of variables generated for each instance. The fourth column reports the average number of iteration of the dual simplex executed by the solver. Finally, the fifth column reports the average time spent by the solver finding genotypes are not accurate, i.e., when some single nucleotide polymorphisms are missing or affected by errors. We proposed an exact approach to solution of PPH-UD based on an extend version ofCatanzaro et al. [2010b] class representative model for PPH, possibly one of the best integer programming models described so far in the literature on PPH. We proposed a number of reduction rules to improve the performances of the model and showed how to embody the perfect phylogeny constrains in order to represent known familial relationships among members of the input population.Afterwards, in Chapter 3 we investigated the Parsimonious Loss of Heterozygosity Problem (PLOHP), a NP-hard combinatorial optimization problem consisting of finding the minimum cost clique cover problem on a MPTG. The optimal solution to the PLOHP has a remarkable importance in practice, as it enables the association of major human diseases with chromosomic regions from patients that exhibit loss of heterozigosity events. Practical instances of the PLOHP are usually characterized by a very large size and this fact usually prevents the use of exact solution approaches such as those described in Catanzaro et al.[2013] andHalldórsson et al. [2011b]. Moreover, we investigated ways to speed up the best known exact solution algorithm for the PLOHP as well as techniques to enlarge the size of the instances that can be optimally solved. In particular, we presented a Branchand-Price (B&P) algorithm for the PLOHP and we developed a number of preprocessing techniques and decomposition strategies to dramatically reduce the size of its instances. Computational experiments showed that the proposed algorithm is 10-30x faster than previous approaches described in the literature. The new algorithm enables the solution of very large practical instances of the PLOHP, containing over 6 000 trios and SNPs.Finally, in Chapter 4 we investigate the Multiple Individuals Polymorphic ALU Insertions Recognition Problem (MIPAIRP), i.e., the problem of partitioning suspected polymorphisms of a set of individuals into the minimum number of deletion areas. Specifically, we showed that the MIPAIRP can be formulate as particular instance of the clique partition problem in a undirected interval graph G, and we prove the general NP-hardness of the problem. Moreover, we extended the results described in[START_REF] Sveinbjörnsson | PAIR: polymorphic Alu insertion recognition[END_REF] by providing a state-of-the-art integer programming formulation and a possible strengthening valid inequalities able to exactly solve real instances of the MIPAIRP containing up to 3 000 individuals and 3 000 SNPs within 12 hour computing time. The overall performances could be further enhanced by implementing the preprocessing techniques presented in Chapter 3.

  contributions and the observations made in our work helped illustrate how profitably the sophisticated optimization techniques of OR can be applied to a non-mathematical domain such as GWAS. Moreover, our work also pose a number of interesting open questions for the specific research issues attacked in this thesis and for the research in the field of GWAS, in general.In particular, in Chapter 2 we have investigated the Pure Parsimony Haplotyping under Uncertain Data (PPH-UD) which take in account input genotypes that are not accurate. The lack of algorithm able to take in account bias and errors is one of the well known challenge in the analysis of biological data. In fact, several procedures have been developed for dealing with errors in linkage studies, forensic analy-ses and non-invasive genotyping before using predictive approaches[START_REF] Paolella | The Alu family repeat promoter has a tRNA-like bipartite structure[END_REF]. Our exact approach to solution of PPH-UD based on an extend version ofCatanzaro et al. [2010b] proposed in Chapter 2 should encourage the use of mathematical programming in solving these kind of problems. More so if such possible extensions have been carried out without compromising the performances in solving even larger instances. Furthermore, since the introduction in the last few years of error rate monitoring in the DNA sequencing phase [Manley et al. 2016], the use of a non-binary error mask matrix for the PLOHP deserve for sure future research efforts. Afterwards, in Chapter 3 we investigated the Parsimonious Loss of Heterozygosity Problem (PLOHP), a NP-hard combinatorial optimization problem consisting of finding the minimum cost clique cover problem on a MPTG. The decomposition methods presented allowed usto solve 10-30x faster than previous approaches described in the literature and very large practical instances of the PLOHP, containing over 6 000 trios and SNPs. Another avenue of improvement involves analysing instances that could be solved faster by using a polynomial algorithm. This can not be done in general since PLOHP is NP-hard, but if we are able in polynomial time to characterise a specific instance as an interval graph, instead of a MPTG, the we could use the polynomial algorithm proposed by[START_REF] Gijswijt | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF]. Hopefully, future research efforts will provide new insights on the combinatorics of MPTG and will enable the analysis of even larger genomic fragments.Moreover, in Chapter 4 we investigate the Multiple Individuals Polymorphic ALU Insertions Recognition Problem (MIPAIRP), a NP-hard combinatorial optimization problem consisting of partitioning suspected polymorphisms of a set of individuals into the minimum number of deletion areas. The IP formulation presented allowed to to exactly solve real instances of the MIPAIRP containing up to 3 000 individuals and 3 000 SNPs within 12 hour computing time. Our results give perspective on the development of future approaches to solution of the problem that may turn out to be useful in practical applications.

Table 4

 4 Performances of the CRM for PPH(RM version)

	Table 10	Computational performances comparison of Catanzaro-
		Labbé-Halldórsson Integer Linear Programming (ILP)
		Formulation (CLHF) and the B&P algorithm with
		and without presolving strategies . . . . . . . .
	Table 11	Computational performances comparison of CLHF
		and the B&P algorithm on artificial instances of

on

Brown et al. [2006]

's datasets. . . . . . . . . Table 5 Accuracy of the CRM for PPH-UD under different error ratios. . . . . . . . . . . . . . . . . . . Table 6 Summary of number of nodes, edges and cliques contained in the learning set of instances of the PLOHP used to calibrate the pricing oracle. . . . Table 7 Wilcoxon signed-rank test on solution time for each pair of percentages . . . . . . . . . . . . . Table 8 Statistics of gens-x α , gens-x β , gens-x γ . . . . Table 9 Parameter sets used to generate random instances of the PLOHP. . . . . . . . . . . . . . . . the PLOHP. . . . . . . . . . . . . . . . . . . . . . xx List of Tables xxi

Table

  Statistics of simulated instances Sim1, Sim2, Sim3, Sim4, Sim5. . . . . . . . . . . . . . . . . 74 Sim5. . . . . . . . . . . . . . . . . 75 Table Computational performances of the B&P on biological instance of the PLOHP. . . . . . . . . . 76 Table Statistics of Set-x α , Set-x β , Set-x γ . . . . . . . 95

			acronyms	xxiii
	LOH	Loss of Heterozygosity
	Table LP MCPP MILP I N T R O D U C T I O N The minimum, the maximum and the average time (expressed in seconds) both to filter the singletons and the cloned intervals, and to iden-tify the connected components in Sim1, Sim2, (alphabetical order) Linear Programming Minimum Clique Partition Problem Mixed Integer Linear Programming Sim3, Sim4, A C R O N Y M S ALU Arthrobacter Luteus MIPAIRP Multiple Individuals Polymorphic ALU Insertions Recognition 1.1 motivation
		Problem
		ALUG	ALU Graph
	MP Model-based optimization is a key methodology in engineering em-Master Problem
	B&B ployed in the design, analysis, construction, and operation of all kind Branch-and-Bound
	MPPH	Minimum Perfect Phylogeny Haplotyping problem
	SNP	Single Nucleotide Polymorphism
		xxii

Table Performances of Formulation 5 on Set-x α , Setx β , Set-x γ . . . . . . . . . . . . . . . . . . . . . 97 B&C Branch-and-Cut B&P Branch-and-Price CRM Class Representative Model CG Column Generation CLHF Catanzaro-Labbé-Halldórsson ILP Formulation CLOH Consistent with a Loss of Heterozygosity COP Combinatorial Optimization Problem DNA Deoxyribonucleic Acid ELOH Evidence of a Loss of Heterozygosity GWAS Genome-Wide Association Studies IG Interval Graph ILOH Inconsistent whit a Loss of Heterozygosity ILP Integer Linear Programming IP Integer Programming LINE Long INterspersed Elements MPTG Max-Point Tolerance Graph OR Operations Research PAIRP Polymorphic ALU Insertions Recognition Problem PPH Pure Parsimony Haplotyping PPH-UD Pure Parsimony Haplotyping under Uncertain Data PLOHP Parsimonious Loss of Heterozygosity Problem RCI Region Compatible with Insertions RMP Reduced Master Problem RNA Ribonucleic Acid SINE Short INterspersed Elements of devices. Since engineering approaches are playing a growing role in the rapid evolution of systems biology [Wolkenhauer et al. 2005; Sontag 2005; Doyle et al. 2006; Kremling et al. 2007], it is reasonable to expect that mathematical optimization methods will contribute in a significant way to advances in systems biology. In fact, optimization is already playing a key role as given the several examples of applications of optimization in systems biology, offered in Greenberg

  It unfolds through the identification of a number of individuals carrying a disease or a trait and proceeds by comparing these individuals to those that do not carry or are not known to carry that disease/trait. Both sets of individuals

are genotyped for a large number of genetic variants which are then tested for association with the disease/trait. Studies using tens of thousands of individuals are becoming common (or are widespreading) and are becoming increasingly the norm/standard in the study of association of genetic variants to disease

[START_REF] Styrkarsdottir | Multiple genetic loci for bone mineral density and fractures[END_REF][START_REF] Gudbjartsson | Many sequence variants affecting diversity of adult human height[END_REF][START_REF] Rivadeneira | Twenty loci associated with bone mineral density identified by large-scale metaanalysis of genome-wide association datasets[END_REF]

]. Since October 2010, 702 human GWASs have been published on 421 traits, the majority exhibiting medical relevance. "The National Human Genome Research Institute at the National Institutes of Health" weekly updates a catalogue of published GWAS results (http://www.genome.gov/gwastudies).

Table 1 :

 1 mances of the CRM for PPH-UD, we compared the numerical results of the model with the corresponding ones of the CRM for PPH(RM version, see Catanzaro et al. [2010b] running on the same datasets in absence of uncertainty. The performances of the CRM for PPH are shown in Table4. Moreover, in order to obtain a qualitative measure of the accuracy of the CRM for PPH-UD, we compared the optimal solutions provided by CRM for PPH in absence of uncertainty with the corresponding ones provided by CRM for PPH-UD. Specifically, fixed a generic instance of PPH-UD, we computed the accuracy of the CRM for PPH as the ratio between the number of matching haplotypes in the solutions provided by both the CRM for PPH-UD and the CRM for PPH divided the overall number of haplotypes in the solution provided by CRM for PPH. The accuracy of the CRM for PPH-UD under increasing error ratios it is shown in Table5. For sake of notation, in the following subsections we shall denote CRM1 and CRM2 as the CRM for PPH and the CRM for PPH-UD, respectively.The numerical experiments relatives to the uniform datasets showed that CRM2 takes sensibly more time than CRM1 to solveBrown et al. [2006]'s datasets, confirming the hardness of PPH-UD with respect to PPH. Specifically, while CRM1 took in average 8 seconds to solve the most difficult dataset having 10 SNPs, CRM2 took in average 10 seconds independently from the error ratio, and even longer on instances 03, 05, 06, 08 and 11 of dataset 50x10r4 where 19.657, 62.160, 34.429, 40.416, and 15.974 seconds were needed to find the optimum. This trend persists also in the instances having 30 SNPs, in which When considering instances having a larger number of SNPs, we experienced a generalized increment of the average solution time taken by CRM2, proportional to the increment of the error ratio. Performances of the CRM for PPH-UD when considering input data having an error ratio of 5%.

	2.3.3.1 Uniform Datasets

CRM1 took in average 11.772 seconds while CRM2 needed an average solution time of 43.273 seconds when considering an error ratio of 5%, with the exception of instances 02, 08, 09, 11, 13 and 14 which needed

62.182, 51.462, 60.079, 60.020, 122.443, and 58

.514 seconds, respectively. We experienced also a generalized decrease of the average solution time when considering an error ratio of 10% and, vice versa, an increase of the average solution time when considering an error ratio of 15%.

Table 2 :

 2 Performances of the CRM for PPH-UD when considering input data having an error ratio of 10%.

	Dataset		Time (sec.)			Gap (%)			Nodes
		Average	Max	Min	Average Max	Min Average Max Min
	Uniform
	50x10	10.092	57.006	1.560	0.000	0.000 0.000	1.000	1
	50x10r4	11.570	33.895	3.275	0.000	0.000 0.000	1.000	1
	50x10r16	9.209	18.627	2.790	1.068	8.333 0.000	1.467	5
	50x30	44.276	157.222	7.467	1.123	5.882 0.000	4.133	15
	30x50	39.772	222.039	6.187	1.006	8.333 0.000	42.700	1149
	30x75	49.558	160.645	8.048	1.409	6.250 0.000	19.500	67
	30x100	97.879	262.453 23.986	0.664	3.836 0.000	32.800	98
	Non-Uniform
	50x10	9.925	126.922	2.113	1.000	8.333 0.000	1.133	1
	50x30	38.689	84.860	4.000	0.400	5.882 0.000	11.267	265
	30x50	33.272	73.425	1.793	0.614	4.737 0.000	54.400	66
	30x75	88.030	325.864 32.260	0.800	4.348 0.000	95.400	1539
	30x100	631.495 7207.900 50.928	1.157	4.270 0.000 2371.200 32622	12
	Biological
	CHR10-CEU	27.340	68.299	1.669	0.000	0.000 0.000	53.667	155
	CHR21-CEU	13.455	36.601	0.522	0.000	0.000 0.000	6.667	12
	CHR10-HBC	87.770	248.825	2.228	4.347	7.143 0.015 384.333	1127
	CHR21-HBC	39.878	118.866	0.075	10.317	16.667 0.000	27.000	69
	CHR10-JPT	23.436	69.457	0.002	1.667	5.000 0.000	72.667	213
	CHR21-JPT	2403.091 7188.610 2.781	2.222	6.667 0.000 4406.333 13125
	CHR10-YRI	702.533 1777.450 62.414	0.000	0.000 0.000 583.000	1635	48
	CHR21-YRI	2545.198 7630.880 1.500	5.449	16.346 0.000 123.333	368

Table 3 :

 3 Performances of the CRM for PPH-UD when considering input data

	having an error ratio of 15%.

Table 4 :

 4 Performances of the CRM for PPH (RM version) on Brown et

	al. [2006]'s datasets.

Table 5 :

 5 Accuracy of the CRM for PPH-UD under different error ratios.

Table 2 ,

 2 CRM2 took in average 10 times more than the average solution time taken by CRM1

to solve instances having 10 SNPs, reaching a maximum solution time of 135.199 seconds when tackling instance 06 affected by an error ratio of 5%. Similarly, when tackling instances having 30 SNPs, CRM2 took in average 4 times more than the average solution time taken by CRM1, reaching a maximum solution time of 193.406 seconds when tackling instance 00 affected by an error ratio of 5%. When dealing with instances having more than 30 SNPs, CRM2 took sensibly more than CRM1 reaching a solution time of 7210.800 seconds when tackling the instance 100-30.03 affected by an error ratio of 5%.

Table 6 :

 6 Summary of number of nodes, edges and cliques contained in the learning set of instances of the PLOHP used to calibrate the pricing oracle.

		Nodes	Edges	Cliques
	Maximun	66 670	2 060 044 865	52 522
	Minimum	35 725	742 125 353	16 915
	Average	57 581.35 1 335 575 876.85 36 824.55
	Standard Deviation 8 056.94	313 44 967.49	8 011.84

Table 7 :

 7 Wilcoxon signed-rank test on solution time for each pair of percent-

	ages

Table 9 :

 9 Parameter sets used to generate random instances of the PLOHP.

Table 10 :

 10 Computational performances comparison of CLHF and the B&P algorithm with and without presolving strategies

Table 11 :

 11 Table11shows the performance of CLHF and the B&P in solving the Computational performances comparison of CLHF and the B&P algorithm on artificial instances of the PLOHP.

	"Sim" instances when using the presolving strategies. Both CLHF and
	the B&P algorithm without presolving strategies proved unable to
	solve any instance in this set. The table reports on the mean and
	standard deviation of: the solution time (expressed in seconds) neces-
	sary to solve a random instance of the PLOHP; the gap (expressed in

Table 12 :

 12 Statistics of simulated instances Sim1, Sim2, Sim3, Sim4, Sim5.

Table 11

 11 also shows that the solution time necessary to solve the artificial instances of the PLOHP is longer than the corresponding one necessary to solve any other instances. This fact is possibly due to

			Sim1			Sim2
		Min	Max	Average	Min	Max	Average
	Singleton	0.43	0.52	0.44	0.43	0.52	0.43
	Clone	0.42	0.43	0.42	0.42	0.50	0.42
	Conn. Comp. 6.2•10 -5 1.3•10 -5 6.2•10 -5 6.1•10 -5 1.7•10 -5 6.2•10 -5
			Sim3			Sim4
		Min	Max	Average	Min	Max	Average
	Singleton	0.43	0.52	0.43	0.41	0.42	0.41
	Clone	0.42	0.50	0.42	0.40	0.41	0.42
	Conn. Comp. 6.1•10 -5 6.2•10 -5 6.2•10 -5 6.0•10 -5 1.3•10 -5 6.1•10 -5
			Sim5
		Min	Max	Average
	Singleton	0.73	0.78	0.74
	Clone	0.65	0.66	0.65
	Conn. Comp. 7.7•10 -5 7.9•10 -5 7.8•10 -5

Table 13 :

 13 The minimum, the maximum and the average time (expressed in

	seconds) both to filter the singletons and the cloned intervals, and
	to identify the connected components in Sim1, Sim2, Sim3, Sim4,
	Sim5.

Table 14 :

 14 Computational performances of the B&P on biological instance of

	the PLOHP.
	the larger number of point-intervals present in the artificial instances
	which in turn increases the number of cliques in the corresponding
	MPTG.

  Consider a graph G = (V, E) and a cost function f α defined as in (17). Then, there exists a minimum cost partition into cliques of G, say P , such that none of the cliques in P have cardinality greater than 1 and smaller than α + 1. Moreover, if P contains cliques of cardinality greater than or equal to α + 1 then at least one of them is a maximal clique of G. Then, due to the nature of f α , it is possible to obtain a lower cost partition into cliques of G by just breaking C into |C| cliques of cardinality 1. In fact, in such a case we would have that v∈C f α

	Proposition 9. Proof. By contradiction, suppose there exists a clique C ∈ P such
	that 2	|C|	α.
			and α is a positive integer such that 2	α	|V| -1. In
			fact, in such a case it is easy to see that the set function f α (C) is
			size-defined submodular, hence value-polymatroidal; thus, if G(I) is
			an ALUG, by Proposition 8 the PAIRP can be solved in polynomial
			time. Moreover, it is worth noting that the optimal solution to the
			problem can be characterized when considering function (17). In fact,
			the following proposition holds:

tively.

Table 8 shows the main characteristics of the 15 instances. In particular, the table is composed by three main sections, each one corresponding to the datasets "gens-x α ", "gens-x β " and "gens-x γ ", x ∈ {1, 5}, previously described. The first three columns of Table 8 re Proof. The statement follows by observing that a RCI is present in

cardinality bigger or equal to 2. This is impossible if for all the cliques

where the cardinality of the clique C is equal to the number of individuals belonging to that clique.

Based on Proposition 10, the new following cost function seems much more well-suited:

where c1, c2 and c3 are three constants such that 0 < c 1 αc 1 < c 2 c 3 (α + 1)c 1 , α is a positive integer such that 2 α < |V| -1,

and ω C the number of individuals contained in the clique C ∈ C(G).

Unfortunately, in such a case it is easy to see in the following example that the new set function ( 18) is not size-defined submodular, hence value-polymatroidal, but a simple submodular rank function.

Given the graph G in Figure 16, S = {a, b, c}, T = {e, f, g}, u = {d} and c 1 = 1 and c 2 = c 3 = 2. The main property of the polymatroid set function for function ( 18) is not satisfied because:

gives as result that:

Implementation

We modelled the MIPAIRP with PYOMO [START_REF] Hart | Pyomooptimization modeling in python[END_REF]], a opensource software package that supports a diverse set of optimization capabilities for formulating, solving, and analyzing optimization models, using Gurobi Optimizer version 6.5.2 build v6.5.2rc1. The experiments were performed on an Intel Core i7-4930K CPU, 3.40GHz, equipped with 64 GByte RAM and operating system Ubuntu release 16.04.3 LTS (kernel GNU/Linux 4.4.0-51-generic).

Benchmark Instances

We tested our model on three datasets, hereafter denoted as "Set-1", "Set-2" and "Set-3", characterized by having 50 000, 100 000, 150 000,