PTree: pattern-based, stochastic search for maximum parsimony phylogenies

Method overview

In this section, we will describe our method in a bottom-up 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 Pattern-Based Tree Reconstruction method that reconstructs a phylogenetic tree topology from a sequence alignment, employing the RSP and Intermediates Sampling. In our software PTree the Pattern-Based 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₀ its parent and V₁…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₀ 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 non-empty 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 non-empty 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 non-empty repeated substitution set Y_k (i.e., the intersection between at least two distinct substitution sets is non-empty). 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 non-empty 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₀ its parent that represents DNA sequence (ATACC). V₁…V₃ 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₀ are described via substitution set S₀. Sequences of A and V₀ 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 non-empty, namely S₀∩S₁ = S₀∩S₃ = S₁∩S₃ = {C2T} = Y₁. Thus, we can create a sequence of intermediate node I₁, such that we apply the mutations in the intersection set Y₁ = {C2T} to the DNA sequence of node A(ACACA), which results in intermediate sequence I₁(ATACA), i.e., rPattern(A) returns a list of length one that contains I₁. Figure 2B depicts an expected local tree topology after node I₁ is added. After the new internal node I₁ 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 two-step 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₁ n), c₂) where c₁ and c₂ are fixed coefficients (e.g., c₁ = 4.0 and c₂ = 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₁ 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:

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

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

3. Remove intermediate nodes of degree one or two.

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

5. 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.

Pattern-based 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:

1. Initialize the distance matrix (DM), i.e., compute a neighbor-joining (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.

2. 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.

3. 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 depth-first 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 two-step 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.

4. 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.

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

6. 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.

7. Return the resulting phylogenetic tree topology.

Let us denote the application of this pattern-based 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 pattern-based 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 burn-in 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:

1. The algorithm initializes the search with a neighbor-joining (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 neighbor-joining algorithm with computational complexity O(N²logN) (Sheneman, Evans & Foster, 2006). Distance corrections based on the Jukes–Cantor (Jukes & Cantor, 1969) or the Kimura two-parameter 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.

2. 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 pattern-based 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.

3. 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.

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 Robinson-Foulds 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 distance-based 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 250-sequence 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.