 Software
 Open access
 Published:
PQ, a new program for phylogeny reconstruction
BMC Bioinformatics volume 19, Article number: 374 (2018)
Abstract
Background
Many algorithms and programs are available for phylogenetic reconstruction of families of proteins. Methods used widely at present use either a number of distancebased principles or characterbased principles of maximum parsimony or maximum likelihood.
Results
We developed a novel program, named PQ, for reconstructing protein and nucleic acid phylogenies following a new characterbased principle. Being tested on natural sequences PQ improves upon the results of maximum parsimony and maximum likelihood. Working with alignments of 10 and 15 sequences, it also outperforms the FastME program, which is based on one of the distancebased principles. Among all tested programs PQ is proved to be the least susceptible to long branch attraction. FastME outperforms PQ when processing alignments of 45 sequences, however. We confirm a recent result that on natural sequences FastME outperforms maximum parsimony and maximum likelihood. At the same time, both PQ and FastME are inferior to maximum parsimony and maximum likelihood on simulated sequences. PQ is open source and available to the public via an online interface.
Conclusions
The software we developed offers an opensource alternative for phylogenetic reconstruction for relatively small sets of proteins and nucleic acids, with up to a few tens of sequences.
Background
Phylogenetic reconstruction based on biological sequences is widely used in bioinformatics. Orthologous RNA and protein sequences are used to investigate the evolutionary relationships between taxonomic groups. Molecular biologists investigating protein families often reconstruct the phylogeny of these families to understand the evolutionary origins of important protein features, such as substrate specificity of enzymes.
Many software tools are available for phylogenetic reconstruction, and different tools often produce different results with the same input. At present, several types of phylogenetic algorithms are commonly used. The maximum parsimony (MP) criterion [1] informs the first type of algorithms; these algorithms rate trees using the number of mutations that are required to obtain a given set of sequences. The second class of algorithms are based on probabilistic models of sequence evolution and on the maximum likelihood (ML) criterion [2]. A specific variant of ML algorithms are quartet puzzle (QP) algorithms [3], where the criterion is not the likelihood itself, but the number of quartets of sequences such that the quartet topology induced by a given tree has the maximum likelihood among three possible topologies. The third class of algorithms uses evolutionary distance criteria. These distancebased algorithms vary widely, though the most popular are the neighborjoining algorithm [4] and algorithms based on several varieties of the minimum evolution (ME) criterion.
This paper presents a new characterbased algorithm based on a novel criterion PQ (for positionquartet) that resembles both MP and QP, but significantly differs from that. This new criterion is inspired by the fact that a correct tree often includes a number of branches that split sequences into groups with or without certain characters in certain alignment positions. It seems natural to count such branchcompatible positions and take their number as an optimality score for a tree.
However, the mentioned approach can hardly be applied as is, because branches close to the edges of the tree are more likely to produce a compatible position by chance, compared with branches more central to the tree. Thus, an optimization of this “positionbranch” score would give an advantage to certain tree topologies, namely those having less “deep” branches. Moreover, in alignments of a substantial number of sequences, completely compatible positions are rather rare and counting a small number of such positions is not informative.
With these considerations in mind, our method counts branchcompatible positions, not in the whole tree, but instead in its fourleaf subtrees, which have only one branch each. The topology of a tree is known to be unambiguously determined by the topologies of its fourleaf subtrees. At the same time, many branchcompatible positions should occur in a foursequence alignment. Hence, a correct tree should contain more alignment positions that support splits of the fourleaf subtrees, relative to an incorrect tree.
We propose the positionquartet (PQ) score, which counts the number of pairs of an alignment position and a quartet of sequences such that the position supports the subtree of the quartet. In the simplest variant (which is used for nucleic acid alignments) “support” means that one side of the quartet contains the same letter in the position and both letters on the other side are some other ones. The mentioned sides of any quartet are uniquely determined by the topology of the tree. If a position provides a “double” support (i.e., one letter in both sequences from one side and some other letter in both sequences from the other side of the quartet), then such positionquartet pair counts twice.
A refined version of the PQ score relies on the fact that in proteins, a specific feature of a clade may not be a single aminoacid residue at a certain position, but instead may represent a group of related residues at the position. This fact inspired us to use scoring matrices for amino acid residues. More precisely, a position supports a quartet, if the value of the scoring matrix on two letters on one side of the quartet is greater than on any two letters from different sides. The measure of support is the difference between the matrix value on the supported side and the maximum of matrix values across the split of the quartet. Again, if both sides of a quartet are supported by a position, the measure of support for such positionquartet pair is the sum of two measures. The overall score of a tree topology is the sum of these support measures over all positions of the alignment and all quartets of the tree.
In what follows, we report tests of our program with the BLOSUM62 matrix. We plan to compose a matrix designed especially for phylogenetic reconstruction with PQ, as BLOSUM62 was designed for protein alignment.
The PQ score resembles the parsimony score, as they are both summed over all positions of the alignment. They differ significantly, however, because the PQ score of a position is the sum of scores over all quartets of input sequences, while the parsimony score is the minimal number of mutations needed to produce the letters at a position via the given tree.
The criterion used in the quartetpuzzling (QP) method also resembles the PQ score. In the QP method, the main score is the number of quartets such that the treeinduced topology has the maximum likelihood among three possible quartet topologies. PQ and quartetpuzzling differ in two main respects: first, PQ uses the sum over all positions and all quartets instead of a simple count of quartets; second, PQ does not use the maximum likelihood criterion. In addition, the program TREEPUZZLE [5], which is the only available realization of the quartetpuzzling method, yields a tree as a majorityrule consensus of many trees obtained by stepwise addition in randomized orders of input sequences, while PQ produces the tree with the highest found score.
Our tests show that PQ, MP, and QP yield different results. TNT [6] (a realization of MP) and PQ both produce fully resolved trees, and in all our tests, species trees are more distant from MP trees than they are from PQ trees, on average. TREEPUZZLE (a realization of QP) usually produces unresolved trees, so it cannot be compared with PQ directly. Thus to compare PQ with QP we prepared a script that produces a resolved tree basing on draft trees generated by TREEPUZZLE.
To evaluate the quality of phylogenetic reconstructions performed with PQ, we used natural instead of simulated protein sequences. With the available models of protein evolution, simulated sequence alignments differ from natural alignments in many respects. In the RAxML manual [7], A. Stamatakis writes, “…the current methods available for generation of simulated alignments are not very realistic. …Typically, search algorithms execute significantly less (factor 5–10) topological moves on simulated data until convergence as opposed to real data, i.e. the number of successful Nearest Neighbor Interchanges (NNIs) or subtree rearrangements is lower” and later: “…a program that yields good topological RobinsonFoulds distances on simulated data can in fact perform much worse on real data than a program that does not perform well on simulated data” (p. 60). Our results support the last statement. For example, ME outperforms ML on natural data but is inferior to ML on simulated data.
We used sequence alignments of orthologous proteins for testing; one protein per organism. We compared the reconstructed trees with species trees. We recognize that the actual tree of a given set of orthologous proteins may differ from the species tree because of horizontal gene transfer (HGT) and/or the loss of paralogs, but these deviations should not lead to incorrect conclusions when comparing phylogeny reconstruction methods. If a method reconstructs the actual tree better than another method, then the result from the first method will be closer to the species tree, in most cases. Exceptions to this trend are possible because reconstruction errors can by chance partly compensate for the difference between the real and species trees. Such exceptions, however, will produce just random noise, which is equally likely to improve the results from both methods. If such exceptions are rare, the resulting noise will not influence the comparison significantly. If not, and thus the noise is sufficiently large, then the comparison will yield statistically insignificant results.
On all the sets of alignments we tested, PQ shows a statistically significant (p<0.001) advantage over ML and MP. This indicates that the deviations between actual for each protein and species trees do not significantly affect our conclusions about the results of program comparison.
Testing phylogenetic programs on natural nucleotide sequences is a much more complicated task. We performed just two small tests on extractions from alignments of ribosomal RNA. These tests show that PQ performs well on nucleotide sequences, too.
We also performed tests on simulated protein and nucleic acid alignments. On the simulations, PQ is inferior to ML and MP. Also on the simulations ME and QP have less accuracy than MP and ML, in contrast to our tests on natural sequences. In our opinion, this primarily demonstrates a low quality of simulations made with current mutation models.
Algorithm
Tree score
Consider a multiple alignment of protein sequences and an unrooted binary phylogenetic tree with leaves labeled with the sequences of the alignment. We assume that more than three sequences are present. Let us denote the letter (i.e., an amino acid residue or the gap symbol) in the cth column of the ith sequence of the alignment as a_{ic}. Each fourelement subset {i,j,k,l} of the sequences of the alignment can be divided into two twoelement subsets following by the tree topology. We assume that this division is {i,j}∪{k,l}, which means that the tree contains at least one branch that separates i and j from k and l. We also fix an amino acid substitution matrix S(a,b), such as BLOSUM62.
The tree score Q is calculated using the following formula:
where c accounts for all columns of the alignment, q accounts for all quartets {i,j,k,l} of sequences such that a_{ic}, a_{jc}, a_{kc}, a_{lc} are residues (not gaps), and Q_{cq} (called the positionquartet score or the PQ score) is given by the following formula:
where
For example, if the matrix S(a,b) is diagonal, with all diagonal elements equal to 1 and other elements equal to 0 (which is a natural choice for nucleic acid sequences), then the PQ score Q_{cq} is equal to:

0 if all four letters a_{ic},a_{jc},a_{kc},a_{lc} are different;

0 if the intersection of two sides of the split quartet, {a_{ic},a_{jc}} and {a_{kc},a_{lc}}, is not empty (particularly if all four letters are the same);

1 if a_{ic}=a_{jc} while a_{ic}≠a_{kc}, a_{kc}≠a_{lc}, and a_{ic}≠a_{lc};

1 if a_{kc}=a_{lc} while a_{ic}≠a_{kc}, a_{ic}≠a_{jc}, and a_{jc}≠a_{kc};

2 if a_{ic}=a_{jc} and a_{kc}=a_{lc}, but a_{ic}≠a_{kc}.
We also implemented a generalized variant of the PQ score. It is based on the idea that a quartet that has two pairs of similar letters of both its sides should “cost” more than just a sum of contributions of two sides. Thus it seems natural to multiple the score Q_{cq} of a positionquartet pair (c,q) by a certain number if both sides of the quartet contribute positively to the score.
More precisely, let α be any positive number. Replace the above formula (1) for Q_{cq} with the following:
This formula reduces to (1) if α=1.
Our implementation of PQ includes two ways of accounting gaps, in addition to the default variant in which gaps are ignored. The gap symbol is treated as an additional letter in both variants. One variant makes no difference between gaps and other letters, which denote amino acid residues or nucleotides, and the other accounts for Q_{cq} only if the quartet q in the position c includes one gap at most.
Normalized tree score
Together with the tree score described above, the normalized tree score is computed as follows. For each quartet of input sequences q and each position c the maximum positionquartet score \(Q^{m}_{cq}\) is calculated as the maximum value of the abovedescribed Q_{cq} scores among all three possible splits, regardless of the split realized in the tree. We define Q^{m} as the sum of all \(Q^{m}_{cq}\). Note that Q^{m} does not depend on tree topology, but depends only on the input alignment. Finally, we define the normalized tree score S as the ratio Q/Q^{m}. If the input alignment is fixed, then S is proportional to Q; S simultaneously gives a moreobjective indicator for the treereconstruction quality when considering various alignments. Indeed, Q depends on the total numbers of quartets and positions, while S is the fraction of positionquartet pairs that support the tree and thus does not directly depend on the size of an input alignment. Tests show that both values negatively correlate with distance from the inferred tree to the reference tree, but for all tested sets the correlation coefficient between S and the distance is higher in absolute value.
Search algorithms
For a given alignment, the tree with the highest score must be identified. An exact solution requires factorial time, so we used several standard heuristics to select a tree scored nearly the highest. It is possible that trees with several topologies have the same highest score, in this case, the program returns the one found first.
Stepwise addition
This heuristic fixes the order of the input sequences. For the first four sequences, it finds the tree with the best score, which only requires checking three trees. Then the fifth sequence is added, and the best tree is chosen from the trees with five leaves such that their subtrees with the first four leaves coincide with the tree found at the first step. Sequences are added in this manner until a tree corresponding to the entire set of sequences is obtained.
Multiple stepwise addition
The process of stepwise addition is repeated several times while changing the input order of sequences with random shuffling. The result is the bestscoring tree among all obtained trees.
NNI hill climbing
From an initial tree, such as the result of stepwise addition, this heuristic performs all possible nearestneighbor interchanges (NNI) [8], one by one. If the current NNI yields a tree with a higher score, then that tree is processed again. This heuristic repeats until all NNIs of the current tree yield trees with scores not greater than the score of the current tree.
NNI Monte Carlo optimization
An initial temperature T=T_{ini} is set, T_{ini}=1000 by default, and K=12000000. Only the ratio K/T is significant, so we set K to be large enough to allow T to be expressed as an integer. Then all possible NNIs are performed one by one in an initial tree. If the current NNI gives a tree with a score Q_{new} that is greater than the score Q_{old} of the current tree, then the procedure is repeated with the new tree. If Q_{new}<Q_{old}, then the new tree is next processed with the probability:
and with the probability 1−P the next NNI is performed on the old tree. T is reduced by T_{ini}/N after each step, where N is a parameter, N=1000 by default. The process stops when T reaches zero. The tree with the highest score among all tested is output.
SPR hill climbing
SPR hill climbing is analogous to NNI hill climbing, but uses subtree pruning and regrafting (SPR [9]) instead of NNI.
Materials and methods
Compared software
We compared results of our program with implementations of four wellknown algorithms for phylogenetic reconstruction. These algorithms are: maximum parsimony (MP) implemented in TNT 1.1 [6], maximum likelihood (ML) implemented in RAxML 8.2.8 [7, 10], balanced mimimum evolution (ME) implemented in FastME 2.1.5 [11] and quartet puzzle (QP) implemented in TREEPUZZLE 5.2 [5].
For MP the parameters are as follows:

Program: TNT

Result: RAxML_parsimonyTree

Search strategy: “mult”, which means several rounds of randomized stepwise addition of sequences followed by search using tree bisection and reconnection (TBR).
For our ML tests, we used the PROTGAMMAAUTO model of RAxML for amino acid sequences and GTRGAMMA model for nucleotide sequences. All other parameters remained set at default values. We took the socalled “bestTree” from the output of RAxML, as the result for comparison. The parameters for ML are as follows:

Program: RAxML 8.2.8

Result: RAxML_bestTree

Amino acid substitution model: PROTGAMMAAUTO. This involves automatic model choice and using the gamma distribution of rates; see [7] for details.

Nucleotide substitution model: GTRGAMMA.

Search strategy: starting with MP tree several SPR steps are performed with the radius (i.e. the number of nodes away from the original pruning position) determined automatically by RAxML.
For ME tests, we used FastME 2.1.5 with the default parameters:

Program: FastME v2.1.5.

Amino acid substitution model for distance calculation: LG, gamma rate variation parameter (alpha) equals 1, do not remove sites with gaps.

Initial tree: BIONJ (see [11] for details).

Search strategy: NNI and SPR postprocessing.
For QP tests, we used the program TREEPUZZLE 5.2. This program produces an unresolved tree in general case, which makes impossible a direct comparison with other programs producing resolved (binary) trees. Thus we implemented a script that takes socalled “puzzling step trees” generated by TREEPUZZLE and inputs it to the program consense of PHYLIP [12] package. The latter is able to produce a resolved consensus of a number of trees with socalled extended majority rule. The number of puzzling steps was set to 100, other parameters were by default:

Program: a pipeline from TREEPUZZLE 5.2 to consense.

Substitution model: auto; parameter estimates: approximate.

Rate of site heterogeneity: uniform.

Approximate quartet likelihood.

Number of puzzling steps: 100.

List puzzling step trees.

Consensus type: majority rule (extended)
Data sets of protein alignments
We used three sets of organisms: 25 Metazoa species, 45 Fungi species and 45 Proteobacteria species.
The fungal and proteobacterial species were selected trying to maximize the total number of common Pfam [13] families in their proteomes. Pfam families consist of evolutionary domains, which are segments of proteins whose evolution included only point mutations and small insertions or deletions, without large rearrangements. The evolution of these domains can be studied by analyzing their alignments.
The metazoan species were chosen with the NCBI taxonomy in mind: the goal was a set of popular organisms, with many sequenced proteins and a fully resolved taxonomic tree.
For each set we found as many orthologous groups of protein domains as was possible, using the procedure described in [14]. In brief, this procedure uses the following instructions.
From a set of species, take all Pfam families that are present in all species. For each family, take all sequences of protein domains of this family from all species. Then construct pairwise global alignments of the sequences from different species and compute the alignment scores. Finally, find the best bidirectional hits, which are pairs of domains from different species in which each member of the pair has the maximum alignment score with the other member when compared with all other domains of the same species. An orthologous group is defined as a set of domains, one from each species, such that each pair of the domains forms a best bidirectional hit.
The organisms are listed in Additional file 1, and the sequences of orthologous groups are in Additional file 2.
To examine the relative effectiveness of the programs when analyzing differently sized alignments, we used alignments of subsets of sequences from each orthologous group in addition to alignments of entire orthologous groups. We thus tested the programs on nine alignment datasets, as listed in Table 1.
Each metazoan orthologous group was randomly split into 10 and 15 sequences; each fungal or proteobacterial orthologous group was split into 15 and 30 sequences. All the sets of sequences so obtained were aligned using Muscle 3.8.31 [15].
An alignment was removed from the dataset if: (i) it contains two or more identical sequences, or (ii) the distance matrix (generated by the protdist program of the PHYLIP package) contains negative distances, meaning that some sequences are too distant so that the distance likelihood function has no maximum. This explains why, for example, the Metazoa25 dataset contains fewer alignments than the Metazoa15 dataset.
Comparison procedure for protein alignments
To compare two fullyresolved (binary) trees for the same set of species, we use the normalized Robinson–Foulds distance [16], which is the number of different splits in the two trees, divided by the total number of splits in the trees. This value remains between 0 and 1.
A species tree was created for each dataset. For the 25 metazoa, we designed this tree to be unique, as supported by the NCBI Taxonomy, since the 25 species were selected to ensure this. For the 45 fungi, the species tree is the consensus of all trees that were built by all four of the MP, ML, ME and QP methods using all alignments of 45 fungal domains. This consensus was created with the program consense of PHYLIP package using the “Majority rule extended” option, which yields a binary consensus tree. The same procedure was used for the 45 proteobacteria.
For the other datasets we studied, namely Metazoa10, Metazoa15, Fungi15, Fungi30, Proteobacteria15 and Proteobacteria30, species trees were obtained by restricting the corresponding complete tree to the appropriate subset of organisms.
All three of the complete species trees are included in Additional file 2 in Newick format and as PNG images.
For each alignment, we computed the normalized Robinson–Foulds distances between the corresponding species tree and the trees created by the five methods: PQ, MP, ML, ME and QP.
For each alignment dataset, we compared the results from PQ with results from MP, ML, ME, and QP. To compare PQ with, for instance, MP, we counted the number of alignments for which the distance from the PQ tree to the species tree is less than the distance from the corresponding MP tree to the species tree. We also counted the number of alignments for which the distance from the PQ tree is greater than the distance from the MP tree. These two numbers were then compared by the sign test. If the pvalue is less than 0.001, then one of the compared methods is judged to be more effective for the present dataset.
As a reference for fungal and proteobacterial alignments, we may use the consensus of trees created by any one program alone with almost the same results. All three consensus trees are close to each other. For Fungi, the maximum normalized Robinson–Foulds distance 2/42≈0.048 occurs between the MP and ML consensus trees, meaning that each tree contains two splits of 42 that are not presented in another tree. For Proteobacteria, the maximum distance 8/42≈0.19 occurs between the ME and QP consensus trees. The comparison results depend only slightly on the choice of the reference tree. For example, comparing PQ with ML on Proteobacteria30, the result is 430/186 using the overall consensus as a reference, i.e., in 430 cases the PQ reconstruction is closer to the reference and in 186 cases it is farther. Compare these values with 428/195 using the ML consensus, 429/181 using the ME consensus, 421/183 using the MP consensus, and 431/194 using the QP consensus; these are all quite close to each other.
Datasets of nucleic acid alignments and comparison procedure for them
To produce a good reference dataset of nucleic acid alignments is a much more complicated task comparing to the same one for protein alignments. We decided to perform a rather small test to check the ability of PQ to reconstruct phylogeny from a set of nucleic acid sequences.
For 45 fungi and 45 proteobacteria that are involved in the protein test, we downloaded their small ribosomal RNA from the database Silva [17]. We aligned these two sets of RNA sequences by Muscle, then excluded redundant sequences (there are two pairs of completely identical rRNA sequences in the fungal set), also, we removed all sites represented by only one sequence. The resulting alignments consist of 43 sequences and 1853 columns for Fungi and of 45 sequences and 1666 columns for Proteobacteria, these alignments are available in Additional file 3. Then 100 times for Fungi and 100 times for Proteobacteria we performed the following procedure: random selection of a number N from the range 300 to 800; random selection of 15 species and N columns from the alignment; composing an artificial alignment from these rows and columns. The resulting set of 200 artificial subalignments was used for testing programs. These subalignments and trees inferred from them are available in Additional file 3. We used restrictions of our species trees to corresponding species subsets as reference trees.
Simulated alignments
Amino acid simulated alignments were extracted from raw data to the paper [18] from Dryad Digital Repository [19]. From there we used 500 “reference” alignments from the folder “simulation/30taxa” in the archive rawData.zip. According to that paper, “30sequence multiple sequence alignments were simulated using Artificial Life Framework (ALF) [20]. The sequence length was drawn from a Gamma distribution with parameters k=2.78, θ=133.81. Sequences were evolved along 30taxa birth–death trees (with parameters λ=10μ) scaled such that the distance from root to deepest branch was 100 point accepted mutation (PAM) units. Characters were substituted according to WAG substitution matrices [21], and insertions and deletions were applied at a rate of 0.0001 event/PAM/site, with length following a Zipfian distribution with exponent 1.821 truncated to at most 50 characters (default ALF parameters).”
Five hundred nucleotide 15sequence alignments were simulated using phylosim R package [22]. The trees for simulations were created by rtree function from the phylosim package with parameters by default, which means branch lengths uniformly distributed in interval 0 to 100 PAM. The length of the initial sequence was chosen uniformly from 300 to 800. Characters were substituted according to GTR substitution model with a mutation rate heterogenity modeled according to a Gamma distribution with the shape parameter of 4.5 and the fraction of invariant sites of 0.5. Insertions and deletions were applied at a rate of 0.0045 event/PAM/site, with the maximum length of 4. The simulated nucleotide alignments, the trees used for simulations, and the trees inferred from the alignments are available in Additional file 4.
Implementation
We implemented PQ in a commandline application written in ANSI C. The source code, an executable file for Windows, and a brief user manual are available at http://mouse.belozersky.msu.ru/software/pq/.
The program takes an alignment in Fasta format as input and outputs an unrooted tree with no branch lengths in Newick format. Users may select a number of parameters, among them the file with the scoring matrix, the positive integer value of α, and the optimization strategy to be used. Further details are available in the online user manual.
A web interface is available at http://mouse.belozersky.msu.ru/tools/pq/. It allows the reconstruction of phylogeny from alignments of up to 100 sequences using any optimization strategy except for SPR. For user convenience, the web interface returns an unrooted tree without branch lengths along with a rooted phylogram that has the same topology. Branch lengths are computed by the program proml in the PHYLIP package. The resulting tree with branch lengths is rooted to its midpoint. The program drawgram in PHYLIP is used to generate an image of the tree.
Results and discussion
Time and memory complexity
The time complexity of PQ with parameters by default, i.e., using 10fold stepwise addition followed by gradient NNI search, is C_{1}N^{4}L+C_{2}N^{5}, where N is the number of sequences in the input alignment, L is the number of informative (not completely conserved) sites in the alignment, and C_{1} and C_{2} are coefficients that do not depend on N or L. During the stepwise addition, calculation of Q_{cq} for all alignment columns c and all quartets q requires O(N^{4}L) operations. After that \({N\choose 4}\) sums over columns can be stored in memory. Stepwise addition implies N−4 steps of O(N^{4}) operations each, because each step requires testing, in average, (2N−3)/2 branches and testing each branch requires calculations with O(N^{3}) quartets (not O(N^{4}) because the fourth member of each quartet is fixed, it is the added leaf). During the NNI search, each round implies testing N−3 branches, with calculations with O(N^{4}) stored quartets for each branch.
The memory complexity of the program is proportional to \({N\choose 4}\).
Testing on fungal alignments shows that the performance of PQ with default parameters takes for a 30sequence alignment in average 26 times more time and for a 45sequence alignment 223 times more time comparing with a 15sequence alignment. This approximately coincides with the N^{5} rule.
SPR requires more computation time than NNI and the difference grows dramatically with the number of sequences. For alignments of the Metazoa10 dataset, SPR takes on average of 1.3 times more time than NNI hill climbing and 2.5 times more time than single stepwise addition; for Proteobacteria45, the values are 30 times and 210 times, respectively. Theoretical considerations give the sixth power dependence of time with respect to the number of sequences for one round of an SPR search. However, the average number of the rounds also may grow with the sequence number and the rule of this growth is hard to predict theoretically.
Comparing with other programs, the fastest one is FastME. The work of FastME with one 45sequence alignment takes (at our computer) in average 0.13 s. For TNT this time is 0.23 s, for PQ (with parameters by default) is 12 s, for TREEPUZZLE is 100 s and for RAxML is 430 s. Among these programs, PQ has the worst time dependence on the number of sequences. A rough extrapolation shows that PQ would work faster than RAxML up to approximately 150 sequences in the input data.
Tree scores and distances to the species tree
Table 2 lists the mean normalized tree scores S, mean normalized Robinson–Foulds distances to the species trees D, and correlation coefficients: r_{SD} between the scores and the distances, r_{SL} between the scores and the lengths of alignments, and r_{DL} between the distances and the lengths. All data are for trees obtained through NNI hill climbing using the BLOSUM62 scoring matrix. The parameter α was equal to 1, and gaps were ignored. We also tested other values of α, namely 2, 3, 5 and 10, and we took gaps into account, but neither of those improved accuracy, so we omit those results from this paper.
Turning to an analysis of the distances between the reconstructed and species trees, first, notice the difference between fungi and proteobacteria datasets. Trees reconstructed from proteobacterial alignments are on average much more distant from the corresponding species tree than are trees reconstructed from fungal alignments. This divergence may be explained by HGT, which is rather frequent among bacteria. Due to HGT, the real phylogeny of a protein family may differ slightly from the phylogeny of the corresponding organisms, and this difference will increase the distances we consider. Other causes likely contribute to this divergence as well; the lower values of S for proteobacterial datasets hint that specific features of proteobacterial alignments make phylogeny reconstruction more difficult. The correlation r_{SD} between the normalized scores and distances to the species trees is rather stable for all fungal and metazoan datasets and is practically independent of the size of the alignments. For proteobacterial datasets, the values of r_{SD} are also stable with respect to alignment size, but they are significantly lower than those for eukaryotic datasets.
Optimization strategies
For all alignments, we reconstructed phylogenies with PQ using the following six optimization heuristics: single stepwise addition, stepwise addition with randomized order repeated tenfold, 100fold repeated stepwise addition, NNI hill climbing, NNI Monte Carlo search, and SPR hill climbing. Each NNI and SPR search started with the bestscoring result of the tenfold repeated stepwise addition. We measured the frequency at which each heuristic reaches the maximum tree score of the six trees, and how frequently the heuristic produces the minimum Robinson–Foulds distance to the species tree. The results are listed in Tables 3 and 4.
We expected and found that morecomplicated optimization algorithms are required to obtain a maximum possible tree score for alignments of more sequences. Less expected, we found that the difference between complicated and simple optimization algorithms is less for distance to the species tree than it is for tree scores. This likely indicates that the tree score well distinguishes a tree that is far enough from the real tree from a tree that is close to the real tree, but that the score often fails to choose among two nearly correct trees. This trend resembles results obtained by Takahashi and Nei [23] in tests with MP, ML, and ME scores using simulated data.
Analysis of the results presented in Table 3 suggests that proteobacterial alignments have some features that make phylogenetic reconstruction harder than it is with eukaryotic alignments. Note that the data in Table 3 is independent of the species tree and, therefore, does not depend directly on possible HGTs. Nevertheless, with prokaryotic alignments each search strategy reaches the highest tree score less frequently than with eukaryotic alignments of the same number of sequences. This result is in accordance with the lower normalized tree scores for proteobacterial alignments. HTG from taxons other than Proteobacteria may make tree topology more complicated, and this is one possible explanation of the phenomenon.
Another feature that complicates the reconstruction lies in the shorter average length of proteobacterial protein domains, as compared with eukaryotic protein domains. For example, the median alignment length in Fungi45 is 264, and in Proteobacteria45 is 160. The normalized tree score correlates well with the length of the alignment as it is shown in Table 2. But the domain length is not the only factor complicating reconstructions of proteobacterial phylogeny. To check this, we extracted alignments of medium length, namely all alignments of the length between 161 and 263, from Fungi45 and Proteobacteria45. These datasets include nearly equal numbers of such alignments: 222 from Fungi45 and 217 from Proteobacteria45. For these mediumlength alignments, the difference between Fungi and Proteobacteria is also impressive. For example, 100fold stepwise addition gives a maximum score among scores that can be reached with at least one of the heuristics for only 8, which is 3.7%, proteobacterial mediumlength alignments and for 96, which is 43.2%, fungal mediumlength alignments. It means that even working with alignments of the approximately same length, the simple search strategy produces the same result as more complicated strategies much less frequently in case of proteobacteria comparing with the case of fungi.
The behavior of the mean normalized score confirms this lengthindependent relative complexity of proteobacterial alignments. For fungal mediumlength alignments mean value of S is 0.9899, which is lower than that for the total set of fungal 45sequence alignments (0.9912) but higher than that for proteobacteial mediumlength alignments, 0.9795.
Comparison with other programs on protein alignments
We examined the results of NNI hill climbing to compare PQ with other software, and list the results in Tables 5, 6, 7, and 8.
Table 5 contains the average distances to species trees, for each dataset and each tested method.
Table 6 contains the numbers of alignments producing relatively good results. As thresholds for this “relative goodness” we chose the lower quartiles of RF distances among trees built by PQ for each particular dataset, thus these numbers for PQ are always close to 25% of the dataset volume. The percents are not equal to 25% exactly because RF distance takes a limited number of possible values. For example, for Metazoa10 the lower quartile of RF distances between PQ trees and reference trees is 1/7, i.e. the lowest possible nonzero value. Thus for this data set, the percent of good results is equal to the percent of perfect results, i.e. alignments for which the inferred phylogeny coincides with the real phylogeny. For 15species data sets, the percents of perfect results are much lower, 1.2 to 2.3% for Metazoa15, 1.3 to 4.1% for Fungi15 and 0 to 0.3% for Proteobacteria15. For other datasets, there are almost no perfect results of any program.
Table 7 contains the percents of alignments producing relatively bad results. Thresholds are the higher quartiles of RF distances among trees built by PQ for each dataset.
Table 8 contains the results of pairwise comparisons of PQ with ME, ML, MP, and QP, as detailed in Materials and Methods. We conclude from Table 8 that PQ reconstructs phylogeny more accurately than do ML and MP for all the datasets we tested. However, there is a significant point to note about relative accuracy of PQ and ML. The distances between ML trees and species trees correlate with lengths of alignments stronger, comparing with distances between PQ trees and species trees. For example, for Fungi30 the correlation coefficient is − 0.46 for ML trees and − 0.33 for PQ trees, for Proteobacteria30 − 0.22 and − 0.15, respectively. Regarding only alignments of Fungi45 with the length greater than 550, ML has a statistically significant advantage over PQ. Namely among 64 such alignments, for 47 the ML tree is closer to the species tree and only for 11 is more distant than the PQ tree. For all other sets the difference between ML and PQ for long (length > 550) alignments is not significant, but the ratio of two numbers, “ML better” to “PQ better” is always less for long alignments than for short ones. It is not completely clear if this effect is due to the alignment length itself or is related to some features of large proteins.
For sets with alignments of 10, 15, and 25 sequences, PQ is more accurate than ME. The same is correct for the Proteobacteria30 set. For two sets, Fungi30 and Proteobacteria45, the difference between PQ and ME is not statistically significant, and for Fungi45 ME outperforms PQ.
Note that the advantage of ME over both ML and MP accords with G. Gonnet’s results from only, as far as we know, testing phylogeny reconstruction methods on large natural datasets [24]. The commonly held opinion that ML is more accurate than distancebased methods is probably based on tests with simulated alignments, which may differ significantly from alignments of natural sequences.
PQ is more accurate than QP for all metazoan sets, and also for Proteobacteria30. For other sets, the difference between PQ and QP is not statistically significant, but PQ is always slightly better.
Comparison with other programs on nucleotide alignments
Tables 9 and 10 demonstrate results of the five programs on subalignments of rRNA sequences. All programs show medium results for subalignments of fungal 18S rRNA and poor results (average distance to reference is about 0.5) for proteobacterial subalignments. For both sets PQ shows slightly better results comparing with ME and QP and significantly better results comparing with ML and MP. For fungal subalignments ML shows a greater dependence on the subalignment length than other programs, which is in accordance with the same phenomenon for protein alignments.
Comparison with other programs on simulated alignments
Tables 11 and 12 demonstrate results of the five programs on simulated alignments. On amino acid simulations, the best results are demonstrated by ML, MP is much worse, PQ and QP are approximately equal and slightly worse than MP and the worst is ME.
On nucleic acid simulations, MP is the best, even better than ML. Here ME works slightly better than PQ, while QP becomes the worst method.
These results dramatically differ from the results on natural sequences. It means that the used simulation procedures produce alignments that are not realistic and cannot be used for comparison of phylogenetic programs. Probably the natural evolution of biological sequences possesses some properties that are not taken into account by standard algorithms for its computer simulation.
Longbranch attraction
Longbranch attraction (LBA) often occurs as an artifact in phylogenetic reconstruction [25, 26]. We attempted to investigate the frequency of LBA using our set of fungal alignments when reconstructing phylogeny with the four programs we tested. From 827 alignments of 45 fungal orthologous sequences, we selected alignments satisfying the following condition: each of five branches marked by letters in Fig. 1 was reconstructed by at least one of the four tested programs. Two hundred ninety five such alignments are available in our data.
Next, each of those 295 alignments was restricted to 18 sequences. First, we removed sequences from two poorly represented classes: Leotiomycetes and Dothidiomycetes. Then, in each alignment, we found a sequence among Saccharomycetes that has the maximal mean distance from sequences of Eurotiomycetes. The same was done for sequences from Sordariomycetes. The species set for each of the 295 selected orthologous groups consists of all five Basidiomycota, all 11 Eurotiomycetes, and the two mostrapidly evolving sequences, one from Saccharomycetes and one from Sordariomycetes.
To evaluate the degree of unevenness of the evolution rate in our data, for each of 295 selected protein families we computed the ratio of two values: the first is the average distance from the “fastest” sequence of Saccharomycetes to all sequences from Eurotiomycetes, the second is the average distance between sequences of these two classes. For different protein families, this ratio is proved to be between 1.03 and 3.89, in average 1.18.
We intended to observe the attraction of branches adjacent to the two lone and rapidly evolving species (Fig. 2). This effect cannot result from an erroneous selection of orthologs, because with the selected 45species alignments the Pezizomycotina branch was reconstructed.
We reconstructed phylogenies from the restricted alignments with PQ, ME, ML, MP, and QP. For each program, we counted the number of trees containing an erroneous split that separates two lone sequences from others (as in Fig. 2b). Such a split appears in 13 PQ trees, 17 ME trees, 20 ML trees, 20 MP trees, and 18 QP trees.
We repeated the same test, switching Sordariomycetes and Eurotiomycetes so that all Sordariomycetes sequences remained, and the mostrapidly evolving sequences were chosen from Saccharomycetes and Eurotiomycetes. The results of this test are close to those from the first test. An erroneous split appears in 7 PQ trees, 16 ME trees, 13 ML trees, 12 MP trees, and 14 QP trees.
The alignments and trees are available in Additional file 5.
To estimate the statistical significance of comparison of two programs, say A and B, we applied the following procedure. Let m be the number of orthologous groups (among the 295 selected) for which A in both described tests does not produce LBA while B at least once in the two tests produces a tree with LBA. Then let n be the number of orthologous groups with the opposite situation, namely A produces LBA at least once and B in both tests does not make this error. Now compare m and n with the sign test.
This procedure shows that PQ outperforms other tested methods with respect to susceptibility to LBA. Namely, with PQ as A and other programs as B, the results of the described procedure (m/n) are as follows:

B=ME: m/n=15/4,p=0.01

B=ML: m/n=17/2,p=3.6·10^{−4}

B=MP: m/n=20/5,p=2·10^{−3}

B=QP: m/n=13/0,p=1.2·10^{−4}
Conclusion
PQ effectively reconstructs phylogenetic trees following a new characterbased criterion. Our tests indicate that PQ, at least on alignments of 45 and less relatively short sequences, is more accurate than methods that use the maximum parsimony and maximum likelihood criteria. For sets of 10 or 15 sequences, PQ outperforms the FastME program, which is based on the minimumevolution criterion. A test on susceptibility to long branch attraction shows that PQ may be the algorithm least susceptible to this problem. PQ, therefore, provides an effective alternative for phylogenetic reconstruction in some situations.
Also, we confirmed the result of G. Gonnet that distancebased methods (in our case FastME) outperform maximum likelihood in accuracy on natural sequences. This result is not supported by simulation studies that suggest an unsatisfactory quality of the existing simulation algorithms.
Availability and requirements
The datasets supporting the conclusions of this article are included within the article and its additional files.
The described software is available online:

Project name: PQ

Project home page: http://mouse.belozersky.msu.ru/software/pq/

Operating systems: Platform independent

Programming language: C

License: GNU GPL
Abbreviations
 HGT:

Horizontal gene transfer
 LBA:

Longbranch attraction
 MP:

Maximum parsimony
 ML:

Maximum likelihood
 ME:

Minimum evolution
 NNI:

Nearest neighbor interchange
 QP:

Quartet puzzle
 SPR:

Subtree pruning and regrafting
 TBR:

Tree bisection and reconnection
References
Fitch WM. Toward defining the course of evolution: minimum change for a specific tree topology. Syst Zool. 1971; 20:406–16.
Felsenstein J. Evolutionary trees from DNA sequences: a maximum likelihood approach. J Mol Biol. 1981; 17:368–76.
Strimmer K, von Haeseler A. Quartet Puzzling: a Quartet maximumlikelihood method for reconstructing tree topologies. Mol Biol Evol. 1996; 13(7):964–69.
Saitou N, Nei M. The neighborjoining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987; 4:406–25.
Schmidt HA, Strimmer K, Vingron M, von Haeseler A. TREEPUZZLE: maximum likelihood phylogenetic analysis using quartets and parallel computing. Bioinformatics. 2002; 18(3):502–4.
Goloboff P, Farris J, Nixon C. TNT, a free program for phylogenetic analysis. Cladistics. 2008; 24:774–86.
Stamatakis A. The RAxML v8.2.X Manual. 2016. https://sco.hits.org/exelixis/resource/download/NewManual.pdf.
Camin JH, Sokal RR. A method for deducing branching sequences in phylogeny. Evolution. 1965; 19:311–26.
Swofford DL, Sullivan J. Phylogeny inference based on parsimony and other methods using PAUP*: Theory In: Salemi M, Vandamme AM, editors. The Phylogenetic Handbook. A Practical Approach to DNA and Protein Phylogeny. Cambridge Univ. Press: 2003. p. 160–206.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and postanalysis of large phylogenies. Bioinformatics. 2014; 30(9):1312–3.
Lefort V, Desper R, Gascuel O. FastME 2.0: A Comprehensive, Accurate, and Fast DistanceBased Phylogeny Inference. Mol Biol Evol. 2015; 32(10):2798–800.
Felsenstein J. PHYLIP (Phylogeny Inference Package) version 3.6. Distributed by the author. Seattle: Department of Genome Sciences, University of Washington; 2005. http://evolution.genetics.washington.edu/phylip.html.
Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, Potter SC, Punta M, Qureshi M, SangradorVegas A, Salazar GA, Tate J, Bateman A. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016; 44:D279–85.
Spirin S. A benchmark for evaluation of phylogeny reconstruction programs. PeerJ Preprints. 2016; 4:e2628v1.
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004; 32(5):1792–7.
Robinson DR, Foulds LR. Comparison of phylogenetic trees. Math Biosci. 1981; 53:131–47.
Yilmaz P, Parfrey LW, Yarza P, Gerken J, Pruesse E, Quast C, Schweer T, Peplies J, Ludwig W, Glöckner FO. The SILVA and “Allspecies Living Tree Project (LTP)” taxonomic frameworks. Nucleic Acids Res. 2014; 42:D643–8.
Tan G, Muffato M, Ledergerber C, Herrero J, Goldman N, Gil M, Dessimoz C. Current methods for automated filtering of multiple sequence alignments frequently worsen singlegene phylogenetic inference. Syst Biol. 2015; 64(5):778–91.
Tan G, Muffato M, Ledergerber C, Herrero J, Goldman N, Gil M, Dessimoz C. Data from: Current methods for automated filtering of multiple sequence alignments frequently worsen singlegene phylogenetic inference. Dryad Digit Repository. 2015. https://doi.org/10.5061/dryad.pc5j0.
Dalquen DA, Anisimova M, Gonnet GH, Dessimoz C. ALF — a simulation framework for genome evolution. Mol Biol Evol. 2012; 29:1115–23.
Whelan S, Goldman N. A general empirical model of protein evolution derived from multiple protein families using a maximumlikelihood approach. Mol Biol Evol. 2001; 18:691–9.
Sipos B, Massingham T, Jordan GE, Goldman N. PhyloSim — Monte Carlo simulation of sequence evolution in the R statistical computing environment. BMC Bioinformatics. 2011; 12:104.
Takahashi K, Nei M. Efficiencies of fast algorithms of phylogenetic inference under the criteria of maximum parsimony, minimum evolution, and maximum likelihood when a large number of sequences are used. Mol Biol Evol. 2000; 17(8):1251–8.
Gonnet GH. Surprising results on phylogenetic tree building methods based on molecular sequences. BMC Bioinformatics. 2012; 13:148.
Felsenstein J. Cases in which parsimony or compatibility methods will be positively misleading. Syst Zool. 1978; 27:401–10.
Bergsten J. A review of longbranch attraction. Cladistics. 2005; 21:163–93.
Acknowledgements
Authors are grateful to Andrei Alexeevski and Anna Karyagina for fruitful discussions and moral support. We acknowledge the Willi Hennig Society for sponsorhip of free availability of the program TNT.
Funding
SS and DP were supported by the Russian Science Foundation, grant no. 161410319. The funding body did not play any role in the design of the study and collection, analysis, or interpretation of data and in writing the manuscript.
Author information
Authors and Affiliations
Contributions
SS invented an idea of the algorithm. MK wrote the first version of the program and contributed to developing the algorithm and preparation of test sets. DP carried out the programming of the current version of the program. SS performed the comparison with other programs. SS wrote the first version of the manuscript. All authors contributed to improving the manuscript. All authors read and approved the final version of the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 2
Protein data. The archive Proteindata.tar.gz contains nine folders that each hold data of one data set used in this work. Each folder contains two subfolders called Alignments and Trees. Subfolder Alignments contains sequence alignments in fasta format. Names of the files are Pfam identifiers with additional figures, for example, the file PF00012_3.fasta contains an alignment of the sequences of protein domains from the third orthologous group of Pfam family PF00012. Names of the sequences in alignments are Uniprot organism mnemonics. Subfolder Trees contains five subfolders, PQ, MP, ML, ME, and QP with trees in Newick format reconstructed from the alignments with five methods. Names of the tree files correspond to names of alignment files. The subfolder Trees of folders Metazoa25, Fungi45 and Proteobacteria45 also contains the three species trees used as reference, in Newick format and as PNG images. On the metazoan tree image, all nontrivial branches are labeled with taxon names. On the fungal tree image, branches corresponding to phyla, subphyla, and classes of Pezizomycotina are labeled. On the proteobacterial tree image, branches corresponding to classes are labeled. (TAR 28,930 kb)
Additional file 3
Nucleic acid data. The archive Nucleicdata.tar.gz contains two folders called Fungi and Proteobacteria. Each folder contains two subfolders called Alignments and Trees. Subfolder Alignments contains the alignment of small ribosomal subunit RNA of corresponding organisms and 100 15sequence subalignments in fasta format. Names of the sequences in alignments are Uniprot organism mnemonics. Subfolder Trees contains five subfolders, PQ, MP, ML, ME, and QP with trees in Newick format reconstructed from the alignments with five methods. Names of the tree files correspond to names of alignment files. Also the subfolder Trees contains the species trees used as reference, in Newick format and as PNG images. On the fungal tree image, branches corresponding to phyla, subphyla, and classes of Pezizomycotina are labeled. On the proteobacterial tree image, branches corresponding to classes are labeled. (TAR 456 kb)
Additional file 4
Simulated data. The archive Simulateddata.tar.gz contains two folders called Alignments and Trees. Folder Alignments contains 500 simulated 15sequence nucleotide alignments in fasta format. Folder Trees contains six subfolders, PQ, MP, ML, ME, QP and Reference with trees in Newick format reconstructed from the alignments with five methods and with reference trees used for simulations. Names of the tree files correspond to names of alignment files. (TAR 1860 kb)
Additional file 5
Data for testing long branch attraction artifact. The archive LBA.tar.gz contains two folders called Test1 and Test2. Each folder contains two subfolders called Alignments and Trees. Their content are similar to the content of the corresponding files in Additional file 2. Alignments and trees in the folder Test1 each includes sequences from 18 fungal species: five Basidiomycota, 11 Eurotiomycetes, one from Sordariomucetes and one from Saccharomycetes. Alignments and trees in the folder Test2 each includes sequences from 19 fungal species: five Basidiomycota, 12 Sordariomycetes, one from Eurotiomycetes and one from Saccharomycetes. (TAR 1700 kb)
Additional file 6
Scripts. The archive Scripts.tar.gz contains Bash and Python 2.7 scripts used for obtaininh results of the paper and the file ReadMe.txt with their description. (TAR 19 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Penzar, D., Krivozubov, M. & Spirin, S. PQ, a new program for phylogeny reconstruction. BMC Bioinformatics 19, 374 (2018). https://doi.org/10.1186/s1285901823994
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1285901823994