A mixed integer linear programming model to reconstruct phylogenies from single nucleotide polymorphism haplotypes under the maximum parsimony criterion

Problem 1

where T a phylogeny of$\mathcal{H},$$\mathcal{T}$the set of all possible phylogenies of$\mathcal{H},$$f:\mathcal{T}\to \mathbb{R}$ a function modeling the selected criterion of phylogeny estimation, and$g:\mathcal{H}\times \mathcal{T}\to \mathbb{R}$ is a characteristic function equal to one if a phylogeny T is compatible (according to the selected criterion) for the set$\mathcal{\text{H.}}$ A specific optimization problem is completely characterized by defining the functions f and g, and the phylogeny T^∗ that optimizes f and satisfies g is referred to as optimal.

One of the classic criteria for phylogeny estimation is the parsimony criterion, which assumes that one taxon evolves from another by means of small changes and that the most plausible phylogeny will be that requiring the smallest number of changes. That evolution proceeds by small rather than smallest changes is due to the fact that the neighborhood of possible alleles that are selected at each instant of the life of a taxon is finite and, perhaps more important, that the selective forces acting on the taxon may not be constant throughout its evolution[12, 13]. Over the long term (periods of environmental change, including the intracellular environment), the accumulation of small changes will not generally correspond to the smallest possible set of changes consistent with the observed final sequences. Nevertheless, it is plausible to believe, at least for well-conserved molecular regions where mutations are reasonably rare and unlikely to have occurred repeatedly at any given variant locus, that the process of approximating small changes with smallest changes could provide a good approximation to the true evolutionary process of the observed set of taxa[14]. Such an assumption is likely to be reasonable, for example, in intraspecies phylogenetics, where few generations have elapsed since the observed taxa shared a common ancestor and thus the expected number of mutations per locus is much less than one. When such assumtions hold, a phylogeny of$\mathcal{H}$is defined to be optimal under the parsimony criterion if it satisfies the following requirements: (i) it has the shortest length, i.e., the minimum sum of the edge weights, and (ii) it is such that, for each pair of distinct haplotypes$h_i,h_j\in \mathcal{H},$ the sum of the edge weights belonging to the path from h _i to h _j in T^∗ is not smaller than the observed number of changes between h _i and h _j [11]. The first condition imposes the assumption that the smallest number of mutations consistent with the observed sequences is a good approximation to the true accumulated set of mutations; the second condition correlates the edge weights to the observed data.

The parsimony assumption can be considered accurate in the limit of low mutation rates or short time scales, making it a reasonable model for situations such as analysis of intraspecies variation where little time is presumed to have elapsed since the existence of a common ancestor of all observed taxa. Maximum parsimony also remains valuable as a model for novel methodology development in phylogenetics because of its relative simplicity and amenity to theoretical analysis. Novel computational strategies, such as those developed in this paper, might therefore productively be developed and analyzed in the context of maximum parsimony before being extended to more complicated models of phylogenetics.

Finding the phylogeny that satisfies the parsimony criterion involves solving a specific version of the PEP, called the Most Parsimonious Phylogeny Estimation Problem (MPPEP). Some of the variants of the MPPEP, see e.g.,[15, 16], can be solved in polynomial time, however, in the most general case, the problem is$\mathcal{N}\mathcal{P}$-hard[11, 17] and this fact has justified the development of a number of exact and approximate solution approaches, such those described in[11, 17, 18]. Some recent versions of the MPPEP, such as the Most Parsimonious Phylogeny Estimation Problem from SNP haplotypes (MPPEP-SNP) investigated in this article, play a fundamental role in providing predictions of practical use in several medical bioinformatics applications, such as disease association studies[19] or reconstruction of tumor phylogenies[20, 21]. In these contexts, it would be highly desirable to have the most accurate inferences possible for specific applications, but this in turn would imply to have algorithms able to exactly solve instances of such versions. As regards the MPPEP-SNP, the literature describes some (rare) circumstances for which it is possible to solve the problem in polynomial time (see Section Methods). Unfortunately, in the general case the MPPEP-SNP is$\mathcal{N}\mathcal{P}$-hard and solving provably to optimality therefore generally requires the use of exact approaches based on implicit enumeration algorithms, similar to the mixed integer programming strategies described in[1, 2, 22].

In this article, we explore the prospects for improving on the implicit enumeration strategy of[1, 2] using a novel problem formulation and a series of additional constraints to more precisely bound the solution space and accelerate implicit enumeration of possible optimal phylogenies. We present a formulation for the problem based on an adaptation of[23]’s mixed integer formulation for the Steiner tree problem extended with a number of preprocessing techniques and reduction rules to further decrease its size. We then show that it is possible to exploit the high symmetry inherent in most instances of the problem, through a series of strengthening valid inequalities, to eliminate redundant solutions and reduce the practical search space. We demonstrate through a series of empirical tests on real and artificial data that these novel insights into the symmetry of the problem often leads to significant reductions in the gap between the optimal solution and its non-integral linear programming bound relative to the prior art as well as often substantially faster processing of moderately hard problem instances. More generally, the work provides an indication of the conditions under which such an optimal enumeration approach is likely to be feasible, suggesting that these strategies are usable for relatively large numbers of taxa, although with stricter limits on numbers of variable sites. The work thus provides methodology suitable for provably optimal solution of some harder instances that resist all prior approaches. In future work, it may provide useful guidance for strategies and prospects of similar optimization methods for other variants of phylogeny inference.

Notation and problem formulation

Before proceeding, we shall introduce some notation and definitions that will prove useful throughout the article. The human genome is divided in 23 pairs of chromosomes, i.e., organized structures of DNA that contain many genes, regulatory elements and other nucleotide sequences. When a nucleotide site of a specific chromosome region shows a variability within a population of individuals then it is called a Single Nucleotide Polymorphism (SNP). Specifically, a site is considered a SNP if for a minority of the population a certain nucleotide is observed (called the minor allele) while for the rest of the population a different nucleotide is observed (the major allele). The minor allele, or mutant type[24], is generally encoded as ‘1’, as opposed to the major allele, or wild type[24], generally encoded as ‘0’. A haplotype is a set of alleles, or more formally, a string of length m over an alphabet Σ = {0,1}[25].

Given a set$\mathcal{H}$of n haplotypes, denote$\mathcal{S}=\{1,2,\dots ,m\}$ as the set of alleles and h _i (s),$s\in \mathcal{S}$, as the s-th allele of haplotype$h_i\in \mathcal{H}$. Given two distinct haplotypes$h_i,h_j\in \mathcal{H},$ we denote${\mathcal{S}}_{h_i{h}_j}$ as the subset of different alleles between h _i and h _j ,$d_{h_i{h}_j}=\sum _{s\in {\mathcal{S}}_{h_i{h}_j}}\left|h_i\right(s)-h_j(s\left)\right|$ as the distance between h _i and h _j , and we say that h _i and h _j are adjacent if$d_{h_i{h}_j}=1$. From a biological point of view, the adjacency between a pair of distinct haplotypes means that one of the two haplotypes evolved from the other by mutation of a specific SNP over time.

Consider a graph$G=(\mathcal{H},E)$ having a vertex for each haplotype in$\mathcal{H}$and an edge for each pair of adjacent haplotypes$h_i,h_j\in \mathcal{\text{H.}}$ Then, a phylogeny T of$\mathcal{H}$is a spanning tree of G, i.e., an acyclic subgraph of G in which a pair of vertices$h_i,h_j\in \mathcal{H}$ is adjacent in T if$d_{h_i{h}_j}=1$. It is worth noting that, according to the definition of the edge set E, in general a phylogeny of$\mathcal{H}$may not exist as the graph$G=(\mathcal{H},E)$ might not be connected. To always ensure the existence of a phylogeny for$\mathcal{\text{H,}}$ we introduce the set${\mathcal{H}}^{\prime }$ which consists of the minimum number of haplotypes that should be added to$\mathcal{H}$in such a way that, defined$\overline{\mathcal{H}}=\mathcal{H}\cup {\mathcal{H}}^{\prime }$ and$\overline{E}=\left\{\right(h_i,h_j):h_i,h_j\in \overline{\mathcal{H}}\text{and}{d}_{h_i{h}_j}=1\}$, the graph$\overline{G}=(\overline{\mathcal{H}},\overline{E})$ is connected. From a biological point of view, the set${\mathcal{H}}^{\prime }$ represents the set of haplotypes that are ancestors of the observed ones but either had gone extinct or just were not observed in that sample (also called Steiner nodes).

Denote$\overline{T}$ as a phylogeny of$\overline{H}$,$\overline{E}\left(\overline{T}\right)$ as the edge set of$\overline{T}$, and$L\left(\overline{T}\right)$ as the length of the phylogeny$\overline{T}$, i.e., the sum of the distances$d_{h_i{h}_j}$, for all$(h_i,h_j)\in \overline{E}\left(\overline{T}\right)$. Then, the problem of finding a phylogeny of$\mathcal{H}$that satisfies the parsimony criterion consists of solving the following optimization problem:

A mixed integer programming model for the MPPEP-SNP

Denote u _i , i ∈ V, as a decision variable equal to 1 if the i-th vertex of V is considered in the optimal solution to the MPEPP-SNP and 0 otherwise;$x_i^s$ as a decision variable equal to 1 if in the optimal solution to the MPPEP-SNP the s-th coordinate of the vertex u _i , i ∈ V, is 1 and 0 otherwise;$z_{\mathit{\text{ij}}}^s$ as a decision variable equal to 1 if in the optimal solution to the MPPEP-SNP the pair of distinct vertices i,j ∈ V has a change at s-th coordinate, and 0 otherwise; and y _ij as a decision variable equal to 1 if the pair of distinct vertices i,j∈V is adjacent in the optimal solution to the MPPEP-SNP and 0 otherwise. Finally, let$V_{\mathcal{H}}=\{1,2,\dots ,n\}$,$V_{{\mathcal{H}}^{\prime }}=\{n+1,n+2,\dots ,n+\mathit{UB}\}$, and Q = {1,2,…,n + LB}. Then, a valid formulation for the MPPEP-SNP is the following:

Reducing model size

Formulation 1 is characterized by a polynomial number of variables and an exponential number of constraints. Its implementation can be performed by means of standard branch-and-cut approaches that use GSEC separation oracles such as those described in[32].

It is worth noting that variables$x_i^s$ and$z_{\mathit{\text{ij}}}^s$ can be relaxed in Formulation1c)-(1e) and the convexity constraint (1f) suffice to guarantee their integrality in any optimal solution to the problem. Moreover, Formulation 1 can be reduced in size by removing those variables that are redundant or whose value is known in the optimal solution to the problem. For example, variables y _ij can be removed from Formulation 1 as it is easy to realize that they are redundant. Similarly, all variables$z_{\mathit{\text{ij}}}^s$ such that$i,j\in V_{\mathcal{H}}$ and d _ij  > 1 do not need to be defined as their value will be always zero for any$s\in \mathcal{S}$ and in any feasible solution to the problem. Variables u _i , i ∈ Q, do not need to be declared as their value will be always 1 any feasible solution to the problem. Finally, variables$x_i^s$,$i\in V_{\mathcal{H}}$, can be removed as their value is univocally assigned by the input haplotype set$\mathcal{\text{H.}}$ The reduction process can be further combined with the preprocessing strategies described in[1] to obtain even smaller formulations. Such strategies allow one to remove alleles from the input haplotype set$\mathcal{H}$without altering the optimal solution to the problem. For example, suppose that the haplotype set$\mathcal{H}$is such that there exists an allele$ŝ\in \mathcal{S}$ such that$h_i\left(ŝ\right)=1$, for all$h_i\in \mathcal{H}$; then it is easy to realize that$ŝ$can be removed from$\mathcal{S}$as in any feasible solution to the problem the$\text{ŝ-th}$ coordinate of any vertex in the phylogeny would be characterized by having x i ŝ = 1. A similar situation occurs when there exists an allele$ŝ\in \mathcal{S}$ such that$h_i\left(ŝ\right)=0$, for all$h_i\in \mathcal{H}$. Analogously, suppose that the input haplotype set$\mathcal{H}$is characterized by equal alleles, i.e., by the existence of two alleles, say${ŝ}_1$ and${ŝ}_2$, such that$h_i\left({ŝ}_1\right)=h_i\left({ŝ}_2\right)$, for all$i\in \mathcal{S}$. Then it is easy to realize that if one removes one of the two alleles from$\mathcal{\text{S,}}$ say${ŝ}_2$, and multiplies the${ŝ}_1$-th coordinate by 2 does not alter neither the structure nor the value of the optimal solution to the problem. Describing all the preprocessing techniques for shrinking the input haplotype set$\mathcal{H}$is beyond the scope of the present article. The interested reader will find a systematic discussion of this topic in[1].

By applying the previously cited reduction strategies to Formulation 1 and denoting$Ŝ$as the set of alleles resulting from the application of the preprocessing strategies described in[1], w ^s as the number of alleles in$\mathcal{S}$equal to the s-th,$s\in Ŝ$, Z as the set for which variables$z_{\mathit{\text{ij}}}^s$ are defined, R = {n + LB + 1,n + LB + 2,…,n + UB}, and$C_{\mathcal{H}}=\{i\in C:i\in V_{\mathcal{H}}\}$, for any C ⊂ V, the following reduced formulation for the MPPEP-SNP can be obtained:

Strengthening valid inequalities

By exploiting the integrality of variables u _i ,$x_i^s$, and$z_{\mathit{\text{ij}}}^s$, a number of valid inequalities can be developed to strengthen Formulation 2.

Proposition 4

Constraints

$\begin{array}{ll}+& h_i\left(s_2\right)-x_j^{s_2}\le 2(1-z_{\mathit{\text{ij}}}^{s_1})-\sum _{\begin{array}{c}s\in Ŝ:\\ s\ne s_1\end{array}}{z}_{\mathit{\text{ij}}}^s\\ \forall s_1,s_2\in Ŝ:s_1\ne s_2,i\in V_{\mathcal{H}},j\in V_{{\mathcal{H}}^{\prime }}\end{array}$

(7)

$\begin{array}{ll}-& h_i\left(s_2\right)+x_j^{s_2}\le 2(1-z_{\mathit{\text{ij}}}^{s_1})-\sum _{\begin{array}{c}s\in Ŝ:\\ s\ne s_1\end{array}}{z}_{\mathit{\text{ij}}}^s\\ \forall s_1,s_2\in Ŝ:s_1\ne s_2,i\in V_{\mathcal{H}},j\in V_{{\mathcal{H}}^{\prime }}\end{array}$

(8)

are valid for Formulation 2.

Given an input haplotype set$\mathcal{H}$and a pair of non-adjacent haplotypes h _i and h _j , there exit multiple equivalent paths that may connect h _i and h _j in the unary hypercube. This characteristic constitutes the principal class of symmetries for the MPPEP-SNP and may lead to poor relaxations for the problem. For example, suppose that the input haplotype set$\mathcal{H}$is constituted by haplotypes h₁ = 〈00〉 and h₂ = 〈11〉. Then a possible solution to the instance may consist either of a star centered in haplotype h₃ = 〈10〉 or a star centered in haplotype h₃ = 〈01〉(see Figure1). Note that both solutions are feasible and optimal for the specific instance. A possible strategy to break this class of symmetries consists of imposing the following constraints:

Proof

In a feasible solution to the problem, the sum$\sum _{s\in Ŝ}{z}_{\mathit{\text{ij}}}^s$,$i,j\in V_{{\mathcal{H}}^{\prime }}$, i,j∈Z, can assume only value 0 or 1. If$\sum _{\begin{array}{c}j\in V:\\ j\in Z\end{array}}\sum _{s\in Ŝ}{z}_{(i+1)j}^s=0$, constraint (11) reduces to$\sum _{\begin{array}{c}j\in V:\\ j\in Z\end{array}}\sum _{s\in Ŝ}{z}_{\mathit{\text{ij}}}^s\ge 0$ which is trivially valid. Vice versa, If$\sum _{\begin{array}{c}j\in V:\\ j\in Z\end{array}}\sum _{s\in Ŝ}{z}_{(i+1)j}^s=1$, constraint (11) reduces to$\sum _{\begin{array}{c}j\in V:\\ j\in Z\end{array}}\sum _{s\in Ŝ}{z}_{\mathit{\text{ij}}}^s\ge 1$ which is again valid due to Propositions 1 and 2. □

Proposition 6 forces vertices in$V_{{\mathcal{H}}^{\prime }}$ to be sorted according to a decreasing degree order. In this way, it is possible to avoid the occurrence of symmetric solutions to the problem differing just for the degree of the Steiner vertices (see e.g., Figure2).

Let$Q_3=\{i,j\in V_{\mathcal{H}}:d_{\mathit{\text{ij}}}\ge 3\}$ and k ∈ V, k ∉ Q₃. Then the following proposition holds:

Implementation

We implemented Formulations 1 and 2 by means of Mosel 64 bit 3.2.0 of Xpress-MP, Optimizer version 22, 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 both formulations, we computed the overall solution time to solve a generic instance of the problem as the sum of the preprocessing time due to the implementation of[22]’s reduction rules plus the solution time taken by the Optimizer to find the optimal solution to the instance. In preliminary experiments, we observed that Formulation 2 has two main advantages with respect to Formulation 1, namely: it requires much less memory to load the model (at least 27% RAM less in the analyzed instances) and it is characterized by faster linear relaxations at each node of the search tree. As result, Formulation 2 allows potentially the analysis of larger instances than Formulation 1 and may be characterized by faster solution times. Hence, we preferred to use Formulation 2 in our experiments.

In preliminary experiments we observed that the Flow-RM outperforms the GESC-RM in terms of solution time. This fact is mainly due to the computational overhead introduced by the GSEC separation oracle which seems to be not compensated by the size of the analyzed instances. Hence, we did not consider the GESC-RM any further in our experiments.

During the runtime, we enabled the Xpress-MP automatic cuts and the Xpress-MP pre-solving strategy. Moreover, we also tested a number of generic primal heuristics for the Steiner tree problem to generate a first primal bound to the MPPEP-SNP (see, e.g.,[34]). Unfortunately, in preliminary experiments we observed that the use of such heuristics interferes negatively with the Xpress Optimizer, by delaying the solution time of the analyzed instances. Hence, we disabled the used of the generic primal heuristics and enabled the use of the Xpress-MP primal heuristic instead. The source code of the algorithm can be downloaded athttp://homepages.ulb.ac.be/~dacatanz/Site/Software_files/iMPPEP.zip.

Separation oracle for the forbidden path constraints

to the formulation; otherwise, we select a different pair of vertices in V∖Q _d and iterate this procedure until either 10 distinct paths have been generated or all possible pairs of vertices in V∖Q _d have been selected. If all vertices have been selected but less than 10 distinct paths have been generated, then we select a larger subset of V∖Q _d (e.g., a triplet of vertices) and we iterate again the previous steps. It is easy to realize that this procedure may potentially generate all the possible paths violating (13). However, we stop the procedure after generating 10 paths or after considering subset C containing more than 5 vertices as we observed in preliminary experiments that this strategy provides the best trade-off between speed and tightness of the cut.

Branching strategies

We observed that even better runtime performance can be obtained by sorting the coordinates of the input haplotypes in decreasing way according to the weights w ^s and by branching first on variables${\beta }_{\alpha }^s$, then on variables u _i , and subsequently on variables$x_i^s$ and finally on variables$z_{\mathit{ij}}^s$.

Performance analysis

In order to obtain a measure of the performance of the Flow-RM, we compared[1]’s polynomial-size formulation versus the Flow-RM on[1]’s real instances of the MPPEP-SNP, namely: Human chromosome Y constituted by 150 haplotypes having 49 SNPs each; bacterial DNA constituted by 17 haplotypes having 1510 SNPs each; Chimpanzee mitochondrial DNA constituted by 24 haplotypes having 1041 SNPs each; Chimpanzee chromosome Y constituted by 24 haplotypes having 1041 SNPs each; and a set of four Human mitochondrial DNA from HapMap[35] constituted by 40 haplotypes having 52 SNPs each, 395 haplotypes having 830 SNPs each, 13 haplotypes having 390 SNPs each, and 44 haplotypes having 405 SNPs each, respectively. Such instances consist only of non-recombining data (Y chromosome, mitochondrial, and bacterial DNA) and can be downloaded athttp://homepages.ulb.ac.be/~dacatanz/Site/Software_files/iMPPEP.zip.

Interestingly, sometimes in real applications the number of haplotypes can be much bigger than the number of SNPs. Hence, it is important to test the ability of an exact algorithm to tackle instances of the MPPEP-SNP containing e.g., hundreds haplotypes.[1] observed that their brute force enumeration algorithm is able to tackle instances of the problem containing up to 270 haplotypes having up to 9 SNPs each within 12 hours computing time. Unfortunately, the authors also observed that their algorithm is unable to solve larger instances of the MPPEP-SNP, no matter the maximum runtime considered. In this context, the Flow-RM makes the difference, being able to tackle instances of the MPPEP-SNP having up to 300 haplotypes and 10 SNPs within 3 hours computing time. To show this result, we considered a number of random instances of the problem containing 100, 150, 200, 250, and 300 haplotypes, respectively. Fixing the number of haplotypes n∈{100,150,200,250,300}, we created an instance of the problem by generating at random n strings of length 10 over the alphabet Σ={0,1}. During the generation process, we randomly selected the number of SNPs equal to 1 in a given haplotype, and subsequently we randomly chose the sites of the haplotype to be set to 1. We iterated the instance generation process 10 times for a fixed value of n, obtaining eventually an overall number of 50 random instances of the MPPEP-SNP downloadable athttp://homepages.ulb.ac.be/~dacatanz/Site/Software_files/iMPPEP.zip.

The results showed that the integrality gaps are usually very low, ranging from 0% to 4.63% and assuming in average a value about 1%, confirming once again the tightness of the Flow-RM and the efficiency of the strengthening valid inequalities.

Finally, we also tested the performance of the Flow-RM on a set of 5 HapMap Human mitochondrial DNA instances of the MPPEP-SNP that were not solvable by using[1]’s brute-force enumeration algorithm, namely: f1 constituted by 63 haplotypes having 16569 SNPs each, i2 constituted by 40 haplotypes having 977 SNPs each, k3 constituted by 100 haplotypes having 757 SNPs each, m4 constituted by 26 haplotypes having 48 SNPs each, and p5 constituted by 21 haplotypes having 16548 SNPs each. Such instances can be downloaded at the same address and consist only of non-recombining data (Y chromosome, mitochondrial, and bacterial DNA).

A part from m4, all the remaining instances gave rise to too large formulations (several hundreds Mbytes RAM) to be handled by the Xpress Optimizer. Hence, instead of analyzing entirely each instance we decomposed it into contiguous SNP blocks and analyzed the most difficult block separately. In more in detail, we define${\mathcal{H}}_r$ to be the haplotype matrix obtained by the application of[1]’s reduction rules, we sorted the columns of${\mathcal{H}}_r$ according to an increasing ordering of the weights w ^s ,$s\in Ŝ$; subsequently, we considered the submatrices obtained by taking k contiguous SNPs (or k-block) in$\text{Ŝ,}$k ∈ {10,13,15}. We did not consider greater values for k as we observed that k = 15 was already a threshold after which the haplotype submatrix gave rise to too large formulations. For each k-block$\mathcal{B}$in${\mathcal{H}}_r$ we considered the hamming distance$d_{h_i{h}_j}=\sum _{s\in \mathcal{B}}\left|h_i\right(s)-h_j(s\left)\right|$ between each pair of distinct haplotypes in${\mathcal{H}}_r$, and chose the k-block maximizing the sum$\sum _{h_i,h_j\in {\mathcal{H}}_r,h_i<h_j}{d}_{h_i{h}_j}$. Finally, we assumed three hours as maximum runtime per instance.