PTree: patternbased, stochastic search for maximum parsimony phylogenies
 DOI
 10.7717/peerj.89
 Published
 Accepted
 Received
 Academic Editor
 Keith Crandall
 Subject Areas
 Bioinformatics, Computational Biology, Evolutionary Studies
 Keywords
 Phylogeny reconstruction, Maximum parsimony, Local search, Stochastic search
 Copyright
 © 2013 Gregor et al.
 Licence
 This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
 Cite this article
 2013. PTree: patternbased, stochastic search for maximum parsimony phylogenies. PeerJ 1:e89 https://doi.org/10.7717/peerj.89
Abstract
Phylogenetic reconstruction is vital to analyzing the evolutionary relationship of genes within and across populations of different species. Nowadays, with next generation sequencing technologies producing sets comprising thousands of sequences, robust identification of the tree topology, which is optimal according to standard criteria such as maximum parsimony, maximum likelihood or posterior probability, with phylogenetic inference methods is a computationally very demanding task. Here, we describe a stochastic search method for a maximum parsimony tree, implemented in a software package we named PTree. Our method is based on a new patternbased technique that enables us to infer intermediate sequences efficiently where the incorporation of these sequences in the current tree topology yields a phylogenetic tree with a lower cost. Evaluation across multiple datasets showed that our method is comparable to the algorithms implemented in PAUP* or TNT, which are widely used by the bioinformatics community, in terms of topological accuracy and runtime. We show that our method can process largescale datasets of 1,000–8,000 sequences. We believe that our novel patternbased method enriches the current set of tools and methods for phylogenetic tree inference. The software is available under: http://algbio.cs.uniduesseldorf.de/webapps/wadownload/.
Introduction
Phylogenetic analysis infers the evolutionary relationships among genes from within or across distinct populations of different species. As input, we are often given a set of aligned genetic sequences for which we wish to work out the pattern of ancestry. This analysis plays an important role in, for example, drug or vaccine development (McHardy & Adams, 2009; Russell et al., 2008). With the cost of sequencing decreasing rapidly due to next generation sequencing technologies (Metzker, 2010), more and more sequences are becoming available and deposited in sequence repositories such as GenBank (Benson et al., 2011). Therefore, new, fast and accurate methods capable of handling largescale datasets are required (Sanderson, 2007).
Commonly used methods for phylogenetic inference fall into four categories: distancebased methods, maximum parsimony, maximum likelihood and Bayesian methods. Distancebased methods are data clustering methods that consider only the pairwise measure of evolutionary distances among sequences. Maximum parsimony assumes that the correct phylogenetic tree is the one requiring the smallest number of evolutionary events to explain the input sequences. A maximum likelihood method requires a substitution model to assess the probability of particular phylogenetic trees, where it aims to find a tree with the highest likelihood with respect to the given substitution model. Bayesian methods rely on probabilistic models of sequence evolution, like maximum likelihood tree inference. Different from maximum likelihood estimates, however, trees are constructed for instance based on the consensus of a set of trees sampled from the highest probability regions of the posterior distribution over evolutionary model parameters and trees. There are several software packages that implement variants of these methods: PHYLIP (Felsenstein, 2005), PAUP* (Swofford, 2002), TNT (Goloboff, Farris & Nixon, 2008), PSODA (Carroll et al., 2009), POY (Varón, Vinh & Wheeler, 2010), MEGA5 (Tamura et al., 2011) and MRBAYES (Huelsenbeck & Ronquist, 2001).
Maximum likelihood and Bayesian methods are considered to be the most accurate, and are nowadays relatively fast due to the computational capacity of modern computers and methodological advances (Guindon & Gascuel, 2003; Huelsenbeck & Ronquist, 2001; Stamatakis, Ludwig & Meier, 2005; Zwickl, 2006). However, these methods depend on large numbers of parameters that have to be estimated, such as branch lengths, tree topology, parameters of the substitution model and sitespecific rate variations, which make them even more computationally demanding. Although maximum likelihood methods are widely used, Kück et al. (2012) showed that the high confidence in maximum likelihood trees is not always justified for certain tree shapes. Distancebased methods, such as the unweighted pair group method with arithmetic mean (UPGMA) and the neighborjoining (NJ) algorithm, are quite fast, with runtime depending only quadratically and cubically, respectively, on the number of input sequences (Saitou & Nei, 1987; Sokal & Michener, 1958). However, distancebased methods are less accurate than all other techniques. We here have focused on the maximum parsimony criterion for tree inference, since this approach has a reasonable tradeoff between speed and accuracy and is still considered as an important optimality criterion for the evaluation of phylogenetic trees (Steel & Penny, 2000), especially for datasets at lower evolutionary divergence (Guindon & Gascuel, 2003). Moreover, Goloboff, Catalano & Farris (2009) showed that maximum parsimony can be successfully employed in the analysis of a largescale dataset of 73,060 taxa. Additionally, they found that long branch attraction (Felsenstein, 1978), which can cause problems when using maximum parsimony, does not play an important role in determining the general structure of the tree when the taxonomic space is sufficiently covered by the sample. Furthermore, a tree computed by a maximum parsimony method can be used as a good starting solution for a subsequent run of a maximum likelihood method (Sundberg et al., 2008).
Finding the optimal solution to the maximum parsimony problem is NPhard (Graham, 1982). The research on exact maximum parsimony searches has resulted in two methods (Kasap & Benkrid, 2011; White & Holland, 2011). Although these methods make use of parallel hardware, they are applicable only on small datasets (<40 taxa). A conjecture as to whether it would be possible to build exact maximum parsimony trees from several trees on fewer taxa turned out to be invalid in general. Fischer (2012) showed that maximum parsimony trees are not hereditary, i.e., an instance of the maximum parsimony problem cannot, in general, be reduced to the same problem on fewer taxa. Therefore, approximation algorithms have been designed to identify good, near optimal solutions. In a commonly used approach, tree rearrangement operators such as nearest neighbor interchange (NNI), subtree pruning and regrafting (SPR), or tree bisection and reconnection (TBR) are applied to define a state space of tree topologies, which is then traversed in local searches with strategies such as best first, best first with backtracking, simulated annealing or genetic algorithms (Felsenstein, 2003). Operators introducing more extensive tree rearrangements such as TBR and SPR usually yield a better solution than less extensive tree rearrangement operators such as NNI; however, this comes with the cost of longer runtimes. Some work has also been devoted to exploring how to modify the size of a neighborhood of a tree that is being explored during one step of a local search algorithm to ensure that the algorithm converges to a near optimal solution and avoids local minima. Goëffon, Richer & Hao (2008) showed that the progressive neighborhood, where the size of the neighborhood is shrinking towards the end of a local search algorithm, results in faster convergence to a near optimal solution.
Genetic algorithms, methods that mimic the process of natural evolution, can be also employed for the maximum parsimony problem (Hill et al., 2005). Genetic algorithms can be also combined with local search methods such that each tree resulting from a recombination operation is subsequently subjected to a local search method that tries to improve its cost. Such algorithms are called memetic algorithms. One implementation of a memetic algorithm is described in Richer, Goëffon & Hao (2009).
Greedy algorithms, which build a solution step by step, can be also employed to reconstruct a phylogenetic tree; such algorithms add one sequence at a time. However, the performance of these algorithms strongly depends on the order in which individual sequences are added. They are thus often combined with other strategies (e.g., local searches). One such an implementation is described in Ribeiro & Vianna (2005).
Here, we describe PTree, a patternbased, stochastic search method for the maximum parsimony problem that can be applied for large datasets (∼8,000 taxa). The method is based on our patternbased tree reconstruction method which is a new technique that can be used to reconstruct phylogenetic trees from genetic sequences; it also suggests how to explore the neighborhood of a tree in the local search. To investigate the performance of the method, we have compared it with the local search techniques implemented in the two widely used software packages PAUP* and TNT.
Design and Methods
Method overview
In this section, we will describe our method in a bottomup fashion. First, we will describe the key idea of our method, which is the Repeated Substitution Pattern (RSP). Given a local tree topology, i.e., a particular internal node of a tree, its parent and children, the RSP defines how the sequences of potential new intermediate nodes, which improve the parsimony cost of the respective local tree topology, can be inferred. As a lot of intermediate nodes may be inferred for an internal node of a tree, we define how this number can be restricted in the subsection Intermediates Sampling. In particular, we aim to choose, primarily, intermediate nodes with the highest potential to lower the parsimony cost of a particular local tree topology. Then, we describe the PatternBased Tree Reconstruction method that reconstructs a phylogenetic tree topology from a sequence alignment, employing the RSP and Intermediates Sampling. In our software PTree the PatternBased Tree Reconstruction method is employed as an operation in a local search algorithm to identify the most parsimonious tree topology. The Randomization of several steps in PTree for further performance improvement is described in the last subsection. In the following, we will use terms node and sequence interchangeably, as in our method, each node of a tree represents one DNA sequence.
Repeated substitution pattern
The RSP defines how new intermediate sequences are inferred at an internal node of the tree, considering only a particular local tree topology. Let us consider the local tree topology depicted in Fig. 1, where A is an internal node, V_{0} its parent and V_{1}…V_{N} are its children, and let S_{i} be the substitution set of node V_{i} consisting of the differences between sequences of V_{i} and A (as seen from A to V_{i}). Note that the arrow that goes from node A to its parent V_{0} does not represent the direction in which the local tree topology is rooted but it serves only for the purpose of the definition of the RSP. A repeated substitution set is defined as the nonempty intersection (Y_{k} = S_{i}∩S_{j}) between two distinct substitution sets S_{i}, S_{j} for i, j∈[0…N] and i ≠ j. Based on this, we infer a new intermediate sequence I_{k} for all distinct nonempty repeated substitution sets Y_{k} at internal node A. A new intermediate node I_{k} originates from the sequence of A and subsequent application of the substitutions contained in the corresponding repeated substitution set Y_{k}. We say that a RSP was found at internal node A, if there is at least one nonempty repeated substitution set Y_{k} (i.e., the intersection between at least two distinct substitution sets is nonempty). If the RSP was found at a particular internal node, the output of the application of the RSP to the internal node of the tree is a nonempty list that contains all new intermediate nodes that can be inferred at this internal node. In the following, we will denote the application of the RSP to internal node A as a call of function rPattern(A) that returns a list of new intermediate nodes (I_{k}).
Figure 2A depicts an example of a local tree topology where node A is an internal node that represents DNA sequence (ACACA), V_{0} its parent that represents DNA sequence (ATACC). V_{1}…V_{3} are children of A that represent DNA sequences (TTACA), (ACCCA) and (ATAAA) respectively. Substitution sets are depicted at respective arrows as seen from internal node A. For instance, the differences in the sequences of A and V_{0} are described via substitution set S_{0}. Sequences of A and V_{0} differ in two positions, namely, C is substituted by T at the second position and A is substituted by C at the fifth position. According to the definition of the RSP, we can see that the intersection of at least two substitution sets (with distinct indices) is nonempty, namely S_{0}∩S_{1} = S_{0}∩S_{3} = S_{1}∩S_{3} = {C2T} = Y_{1}. Thus, we can create a sequence of intermediate node I_{1}, such that we apply the mutations in the intersection set Y_{1} = {C2T} to the DNA sequence of node A(ACACA), which results in intermediate sequence I_{1}(ATACA), i.e., rPattern(A) returns a list of length one that contains I_{1}. Figure 2B depicts an expected local tree topology after node I_{1} is added. After the new internal node I_{1} is added to the local tree topology, its cost decreases from 7 to 5. Thus, the application of the RSP enables the inference of new intermediate nodes that can refine the overall tree topology.
Intermediates sampling
As the application of the RSP, i.e., application of the function rPattern(A) to an internal node A of degree n may produce up to n(n−1)/2 intermediate nodes, we restrict the number of intermediate nodes that can be inferred for a particular internal node with a twostep procedure.
Step 1. We define the maximum number of intermediate nodes that can be inferred at some internal node as a function of the degree of an internal node: F(n) = max(round(c_{1} n), c_{2}) where c_{1} and c_{2} are fixed coefficients (e.g., c_{1} = 4.0 and c_{2} = 1). The purpose of this function is to allow the inference of more intermediate nodes at internal nodes of higher degree. Let I_{A} be the set of intermediate nodes that can be inferred at internal node A (i.e., I_{A}≔rPattern(A)) and M≔F (degree of A). To restrict I_{A} (i.e., the number of intermediate nodes that will be inferred at A), we choose at most M intermediate nodes out of I_{A} with the biggest cost decrease, i.e., we compute for each intermediate node from I_{A} the cost decrease if a node was included in the local tree topology regardless of the other candidate intermediate sequences. For instance, the cost decrease after node I_{1} is included in the local tree topology in Fig. 2 is two. Note that we perform the choice of intermediate nodes at random if we have more options for how to choose intermediate nodes with the biggest cost decrease. Given I_{A}, let us denote the set of at most M intermediate nodes with the biggest cost decrease selected in this step as I_{A_s1}.
Step 2. Given I_{A_s1}, the intermediate nodes selected in the first step, the corresponding internal node A, and its parent and children, we apply the following procedure:

Compute the distance matrix (DM) based on pairwise evolutionary distances among the sequences.

Compute the minimum spanning tree (MST) based on the DM.

Remove intermediate nodes of degree one or two.

Repeat steps (1) to (4) until no intermediate nodes can be removed in step (3).

Return the list of intermediate nodes that remain in the MST; let us denote this list as I_{A_s2}.
In the following, we will denote the application of intermediates sampling as a call of function iSample with a list of intermediate nodes and a respective internal node as its arguments. Thus, iSample(rPattern(A), A) is a list of restricted size I_{A_s2} that contains intermediate nodes with the highest potential to lower the parsimony cost of the local tree topology defined at internal node A.
Patternbased tree reconstruction
This method reconstructs a tree topology from a given sequence alignment that can represent both original input sequences that represent, for example, species, and the corresponding inferred intermediate ancestral sequences, where the latter are auxiliary and optional. Let seqSet be the set of the aligned input sequences. The method comprises seven steps:

Initialize the distance matrix (DM), i.e., compute a neighborjoining (NJ) tree based on pairwise evolutionary distances among the given input sequences seqSet. Set the DM to the path metric representing the distances in the inferred NJ tree, such that for each pair of the input sequences (i.e., leaf nodes of the NJ tree), the DM contains an entry that specifies their distance in the NJ tree. Here, we initialize the DM for the consequent minimum spanning tree (MST) computation in this way, since we observed that our method yields better results than if we computed the MST only based on the pairwise number of differences among the input sequences.

Compute an MST from the DM. The MST represents a tree with the lowest cost considering only the given input sequences seqSet. The MST is an approximation of a minimum Steiner tree (Hwang, Richards & Winter, 1992); therefore it is a good starting point for the intermediate inference, in which we further improve the tree topology (i.e., lower its overall parsimony cost) by adding new intermediate nodes.

Infer new candidate intermediate sequences based on the RSPs. Note that a candidate intermediate sequence denotes a sequence that can be added to the tree topology, while an intermediate sequence refers to a sequence that has been added to the current tree topology. All candidate intermediate sequences, inferred at this step, are collected by the depthfirst search (DFS) algorithm, which efficiently traverses the whole tree and visits all internal nodes to identify RSPs. Thus, we get the list of all candidate intermediate sequences cSeq by calling function iSample(rPattern(A), A) for each internal node A of the tree and by concatenating the resulting lists. More precisely, new candidate intermediate sequences are inferred under consideration of a particular local tree topology, consisting of an arbitrary internal node, its parent, and its children. New candidate intermediate sequences are inferred if a RSP can be found for this local tree topology. Intermediates sampling, the twostep approach, is then employed to get only a restricted number of candidate intermediate nodes at a particular internal node. After all new candidate intermediate sequences have been collected, we remove duplicate candidate intermediate sequences from cSeq (duplicates to each other, or to sequences represented by nodes in the current tree, i.e., in seqSet). Note that inferred intermediate sequences are subsequently added to the sequence set for the tree inference, and can be placed as either internal or leaf nodes in the subsequent tree inference steps. Therefore, the inference of new intermediate sequences plays an important role.

Merge the new candidate intermediate sequences cSeq, inferred in the previous step, with all sequences from the current tree seqSet. Let seqSet be the merged set. Iteratively recompute the distance matrix (DM), i.e., set the entries of the DM to specify the number of different characters between each pair of sequences from seqSet, recompute the MST according to the DM and remove intermediate nodes of degree one (i.e., leaf nodes) or two (superfluous internal nodes), as they are uninformative, until no (candidate) intermediate sequence can be removed. Note that the corresponding sequences are also removed from seqSet. The candidate intermediate sequences that remain in the current tree after this step become intermediate sequences and seqSet contains all sequences of the nodes of the current tree.

Repeat steps (2) to (5) until no new candidate intermediate sequences can be found in step (3).

For each internal node that represents an original input sequence (e.g., existing species), create a copy of the node and place it as a child node of the respective internal node that represents the same sequence, so that all original input sequences are represented by leaf nodes of the resulting tree.

Return the resulting phylogenetic tree topology.
Let us denote the application of this patternbased tree reconstruction method to the set of input sequences seqSet as a call of function pbTree(seqSet) that returns a tree topology.
PTree
Our software PTree employs the patternbased tree reconstruction method as an operation in a local search for the most parsimonious tree topology. PTree searches for an optimal solution to the maximum parsimony problem by traversing the tree space of feasible solutions and keeping only the last feasible solution with the lowest parsimony cost that has been found so far.
In addition to the set of aligned input genetic sequences seqSet, PTree requires specification of the parameters d, b and k: The parameter d denotes the percentages of internal nodes to be deleted from the tree. To allow refinement of the tree topology, the search for the most parsimonious tree topology cannot terminate before b iterations overall (the burnin phase) and before k iterations where no more parsimonious tree relative to the best tree found so far has been observed. For our experiments, we set d to 10%, b to 100 and k to 10.
PTree operates in three steps:

The algorithm initializes the search with a neighborjoining (NJ) tree inferred from the input sequence alignment seqSet. We set the currentTree to be this NJ tree. In PTree, we have incorporated a relaxed neighborjoining algorithm with computational complexity O(N^{2}logN) (Sheneman, Evans & Foster, 2006). Distance corrections based on the Jukes–Cantor (Jukes & Cantor, 1969) or the Kimura twoparameter model (Kimura, 1980) that are available in this implementation can optionally be performed, however, in our experiments we found no notable differences in results with different evolutionary models.

Updating the current tree topology, i.e., delete d% of randomly chosen internal nodes from the currentTree and return rSeq, the set of sequences associated with remaining nodes (internal and leaf nodes). The set rSeq is passed to our patternbased tree reconstruction method, i.e., the function pbTree(rSeq) returns a newTree. If the newTree has a lower parsimony cost than the currentTree, the newTree becomes the currentTree; else, the currentTree remains the same as it was.

Repeat steps (2) and (3) until the number of iterations is larger than b and no update of the currentTree has been performed for the last k iterations. Then, return the currentTree as the final phylogenetic tree.
Randomization
Randomization of some steps of our method substantially improves its accuracy. The MST algorithm, the NJ algorithm and the intermediates sampling are randomized. Since some edges of a graph can have the same weights and any of these edges can be added at a particular step of an MST algorithm, the MST is not unequivocally defined. Standard implementations of the MST algorithm usually return the same tree after two consecutive runs with the same input. Thus, we have modified the MST algorithm so that it is likely that two consecutive runs produce two different tree topologies with the same cost, using the same input sequences. At the point where the MST algorithm can add several distinct edges with the same cost (weight) to the current graph, it decides at random which edge will be added. There can be more options for choosing two clusters that will be merged at a particular step of an NJ algorithm, and thus the resulting NJ tree is also not unequivocally defined; therefore, we have employed a randomized version of the relaxed NJ algorithm (Sheneman, Evans & Foster, 2006). Different MST topologies facilitate inference of more various candidate intermediate sequences, which leads to better accuracy of the resulting phylogenetic trees, since we also explore the larger neighborhood of the current trees in the local search.
Results
To validate the performance of PTree, we compared it with widely used methods on reference datasets of PhyML (Guindon & Gascuel, 2003) and RAxML (Stamatakis, Ludwig & Meier, 2005), as well as on alignments of concatenated Human Immunodeficiency Virus (HIV) reverse transcriptase and polymerase sequences from EuResist (2010). We tested for topological accuracy, cost optimality and runtime. All computations were done on a machine with 2.8 GHz DualCore AMD Opteron Processor 8820. The Java Virtual Machine was given 4GB of the main memory for the datasets up to 4,000 taxa and 8GB for the datasets of 8,000 taxa, and all programs were run in a single thread.
Topological accuracy
Tests for topological accuracy were done in accordance to PhyML (Guindon & Gascuel, 2003) on 5,000 datasets with known topology where each dataset contained 40 taxa of 500bp length each. To compare the true and the reconstructed tree topologies, we used the RobinsonFoulds distance (Robinson & Foulds, 1981), although this measure is more conservative than other measures such as SPR or TBR, we chose it to allow a comparison with results shown for PhyML in Guindon & Gascuel (2003). Reference methods were PAUP* NJ, PAUP* maximum parsimony (with the TBR branch swapping option used in the heuristic search), PAUP* maximum likelihood (with the NNI branch swapping option used in the heuristic search), PhyML and TNT (with the SPR branch swapping option used in the heuristic search). The results are depicted in Fig. 3. As expected, the best results in terms of the Robinson–Foulds distance (Robinson & Foulds, 1981) between the inferred and true tree topology were achieved by the maximum likelihood methods, followed by the parsimony methods (including PTree). The worst results were yielded by NJ, the representative of a distancebased method. Furthermore, within the parsimony methods, PTree outperforms PAUP* used with the tree bisection and reconnection (TBR) heuristic, which is the most extensive and accurate search heuristic.
Cost optimality and runtime
To assess the performance of PTree in terms of parsimony costs and runtime, we compared it with PAUP* NNI, SPR, TBR and TNT SPR heuristics on the reference datasets of RAxML, as well as on a dataset of HIV reverse transcriptase and polymerase sequences. For each dataset, we created seven subsets of different size (125–8,000 sequences). Sequences of the HIV datasets were ∼1,600bp in length and the sequences of the RAxML datasets were ∼1,200bp. For the RAxML subsets, we took the arb_10000 dataset, selected the first 8,000 sequences and then subsequently selected the first half of the remaining sequences for a new dataset. The HIV subsets were created in a similar fashion: we downloaded ∼8,000 reverse transcriptase and polymerase sequences from EuResist (2010), aligned them and repeatedly selected the first half of the remaining sequences as a new dataset. As two consecutive runs of PAUP* (NNI, SPR, or TBR) or PTree can result in different results, we identified the mean (average over multiple runs), minimum, and maximum parsimony costs and runtimes. Thus, PAUP* was run six times (except for the tests with the SPR and TBR heuristics with the RAxML dataset with 4,000 sequences that were run just once; the test with the TBR heuristic with the RAxML dataset with 2,000 sequences that was run twice; and the test with the TBR heuristic with the HIV dataset with 4,000 sequences that was run twice as well); PTree (run with 100 iterations and without distance correction) was run 10 times for all datasets. Some tests were done with lower number of runs due to the limited computational resources. The comparison of the resulting average parsimony costs is shown in Tables 1 and 3, and the comparison of the corresponding average runtimes is shown in Tables 2 and 4. The comparisons of the corresponding minimum and maximum parsimony costs, as well as the comparison of the minimum and maximum runtimes are shown in Tables S1–S8. Furthermore, to enable the comparison of different methods proportionally to the resulting parsimony costs and runtimes of PTree, Tables S9–S20 show the results of Tables 1–4 and Tables S1–S8, respectively, in terms of percentages. In terms of the average values, PTree outperforms PAUP* with the NNI heuristic search in terms of the parsimony costs, except for the 250sequence HIV dataset, where PTree yields slightly higher costs. Moreover, PTree was also substantially faster than PAUP* NNI for large datasets with 1,000 or more sequences. Compared to PAUP* with the SPR or TBR heuristic and TNT with the SPR heuristic, PTree performs worse in all tests in terms of the parsimony costs. Nevertheless, both search heuristics implemented in PAUP* are usually unfeasible due to their high runtime on large datasets with more than 2,000 sequences. However, the SPR heuristic implemented in TNT is applicable for large datasets. While PTree was faster than TNT SPR with the HIV dataset with more than 4,000 sequences, PTree was slower than TNT SPR with the RAxML dataset. The best results in terms of runtime were gained by PAUP* NJ, which shows substantially decreased runtime compared to all tested methods. However, parsimony costs and topological accuracy suggest that NJ should be used only for some initial results.
Size of input dataset  

125  250  500  1,000  2,000  4,000  8,000  
Method  
NJ  6,298  11,789  21,911  42,383  79,567  152,546  289,472 
PAUP* (NNI)  6,132  11,557  21,559  41,636  78,421  150,316  287,530 
PTree  6,104  11,572  21,549  41,573  78,276  149,950  286,815 
TNT (SPR)  6,078  11,487  21,401  41,250  77,556  148,527  283,921 
PAUP* (SPR)  6,083  11,482  21,405  41,274  77,594  148,683  – 
PAUP* (TBR)  6,080  11,476  21,381  41,235  77,537  148,601  – 
Size of input dataset  

125  250  500  1,000  2,000  4,000  8,000  
Method  
NJ  0.2 s  0.2 s  0.5 s  2 s  9 s  35 s  9 m 13 s 
PAUP* (NNI)  2.4 s  11.8 s  1 m 54 s  23 m 44 s  1 h 46 m  11 h 26 m  92 h 52 m 
PTree  19.4 s  58.5 s  3 m 30 s  14 m 13 s  1 h 4 m  4 h 47 m  23 h 48 m 
TNT (SPR)  4 s  14 s  1 m 24 s  6 m 38 s  41 m 2 s  6 h 34 m  33 h 38 m 
PAUP* (SPR)  27.7 s  1 m 30 s  19 m 35 s  2 h 44 m  22 h  155 h 17 m  – 
PAUP* (TBR)  27.7 s  4 m 33 s  28 m 43 s  5 h 2 m  35 h 49 m  190 h 31 m  – 
Size of input dataset  

125  250  500  1,000  2,000  4,000  8,000  
Method  
NJ  957  2,675  6,045  12,704  21,698  71,883  151,678 
PAUP* (NNI)  946  2,650  5,970  12,489  21,188  70,391  148,718 
PTree  934  2,607  5,898  12,299  20,951  70,247  148,423 
TNT (SPR)  934  2,592  5,850  12,214  20,743  69,545  146,906 
PAUP* (SPR)  934  2,597  5,852  12,203  20,718  69,508  – 
PAUP* (TBR)  933  2,589  5,836  12,159  20,673  69,380  – 
Size of input dataset  

125  250  500  1,000  2,000  4,000  8,000  
Method  
NJ  0.1 s  0.1 s  0.1 s  1 s  8 s  1 m  8 m 28 s 
PAUP* (NNI)  1.4 s  11.4 s  1 m 30 s  17 m 7 s  2 h 51 m  28 h 33 m  139 h 49 m 
PTree  7.5 s  38.5 s  2 m 14 s  6 m 43 s  23 m 56 s  2 h 10 m  11 h 10 m 
TNT (SPR)  0.7 s  3 s  14 s  51 s  5 m 23 s  57 m 57 s  4 h 48 m 
PAUP* (SPR)  17.7 s  2 m 44 s  19 m 51 s  2 h 5 m  60 h 23 m  >1 month  – 
PAUP* (TBR)  41.2 s  4 m 56 s  39 m 22 s  4 h 38 m  252 h 42 m  >1 month  – 
Conclusion
PTree is a novel stochastic patternbased search method for maximum parsimony trees. We compared our method with maximum parsimony tree inference using different local search methods with PAUP* and TNT, which are widely used by the bioinformatics community. Our tests showed that PTree outperformed NNI searches in terms of the topological accuracy and cost optimality of identified trees. PTree is also substantially faster than PAUP* NNI searches and comparable to TNT SPR for large datasets (1,000–8,000 sequences). The implementation of our current method is in Java, which is a platform independent programming language (only the NJ algorithm (Sheneman, Evans & Foster, 2006) was added as a C/C++ library). To further improve runtimes, we may reimplement some parts of the algorithm in C and add it as a C/C++ library using the Java Native Interface in the future. We found that more elaborate search heuristics, PAUP* SPR, PAUP* TBR and TNT SPR, yield trees with better parsimony costs than PTree, however they are often more computationally demanding.
Our pattern based tree reconstruction method could be also combined with other algorithms that compute phylogenetic trees. For instance, it could be employed in a memetic algorithm in several ways: as an operation in a local search algorithm, as a mutation operator or as a recombination operator. To recombine two trees, the patternbased tree reconstruction method can be given all previously inferred internal nodes of both trees and all original input sequences as an input, which results in a new tree. Thus, our new patternbased method enriches the set of current phylogeny reconstruction tools and methods.