PnpProbs: a better multiple sequence alignment tool by better handling of guide trees

Background This paper describes a new MSA tool called PnpProbs, which constructs better multiple sequence alignments by better handling of guide trees. It classifies sequences into two types: normally related and distantly related. For normally related sequences, it uses an adaptive approach to construct the guide tree needed for progressive alignment; it first estimates the input’s discrepancy by computing the standard deviation of their percent identities, and based on this estimate, it chooses the better method to construct the guide tree. For distantly related sequences, PnpProbs abandons the guide tree and uses instead some non-progressive alignment method to generate the alignment. Results To evaluate PnpProbs, we have compared it with thirteen other popular MSA tools, and PnpProbs has the best alignment scores in all but one test. We have also used it for phylogenetic analysis, and found that the phylogenetic trees constructed from PnpProbs’ alignments are closest to the model trees. Conclusions By combining the strength of the progressive and non-progressive alignment methods, we have developed an MSA tool called PnpProbs. We have compared PnpProbs with thirteen other popular MSA tools and our results showed that our tool usually constructed the best alignments.

PnpProbs is based on an adaptive approach we proposed in [1], in which we observed that sequences having different similarities have different characteristics and structural properties, and by using some reliable measure to estimate the similarity of the input (we do not know the true similarity because we do not have the correct alignments), we may exploit the corresponding properties to help generate better alignments. To study the feasibility of this idea, we have modified the open source code of MSAProbs [2] and developed a new adaptive MSA tool called GLProbs. Roughly speaking, both tools construct the alignments in the following three stages: (1) Determine the substitution scores for pairwise sequences based on some pair-Hidden Markov model(s), and then refine the scores to make them consistent with all input sequences.
(2) Construct a guide tree and based on it align the input sequences progressively to generate the multiple sequence alignment. (3) Refine the alignment given by Stage (2) to a better alignment for the final output.
The major difference between MSAProbs and GLProbs is in the first stage: MSAProbs uses a single model to determine the substitution scores, while GLProbs determines the scores adaptively. GLProbs first estimates the similarity of the input sequences by computing its average PID (percent identity), which is defined as follows: the PID of two sequences is the percentage of identical columns in their optimal (pairwise) alignment, and the average PID of a sequence family is the average of the PIDs of every pair of sequences in the family. If the input's average PID is high, GLProbs uses the global pair-Hidden Markov model (pair-HMM) to determine the scores; otherwise, it uses some local pair-HMMs.
We have made thorough comparisons between GLProbs and a dozen other leading MSA tools, and GLProbs had the highest accuracy in many of the comparisons (see [1] for more details of our evaluation of GLProbs).
We have also evaluated the performance of PnpProbs on phylogenetic inferencing over two benchmarks, namely Yule-Harding tree simulated data [17] and SABmark empirical data [5]. In Section "Phylogenetic analysis", we compare PnpProbs with five other MSA tools, namely GLProbs, MSAProbs, PicXAA, MUSCLE and ClustalW, and our results showed that the phylogenetic trees generated from the outputs of PnpProbs are closer to the model phylogenetic trees than those constructed from the five other MSA tools.
For verification of our results, all test data can be accessed from [17,18], and PnpProbs can be downloaded via the link https://github.com/ytye/PnpProbs.

Ideas for improving GLProbs
We observe some new structural properties and believe that by exploiting them we can further improve GLProbs' accuracy in general, and improve its accuracy significantly for sequences with low similarity. We focus on improving the second stage of GLProbs. Based on the substitution scores given in Stage 1, this stage determines a guide tree, which is supposed to capture the phylogeny relationship of input sequences. Then, it generates an MSA by performing profile-to-profile alignment according to the order suggested by the guide tree. Unlike GLProbs, we will use an adaptive approach to construct the guide trees. We classify the input sequences into two types: (i) distantly related sequences, whose similarities (or more precisely, average PID) are smaller than some threshold (as suggested by our study in Section "Non-progressive alignment for distantly related sequences", we set it to be 18 %), and (ii) normally related sequences, whose similarities are no smaller than the threshold. PnpProbs handles these two types of sequences differently.
For normally related sequences, we exploit some structure property for better guide tree generation. To explain, we show in Fig. 1 two protein families, F and G, whose average PID's are both 0.33, but their structures are quite different. In particular, each sequence of G has two regions (which we have highlighted in brown colour) over which the sequences are identical, and the sequences are totally different elsewhere. We note that similar structures can be found in real protein families, especially from those with similarity around 20-30 %: ( †) Their sequences have a number of conserved regions over which the sequences are very similar, and the sequences are very different elsewhere.
We have two observations about this structure: (i) The average PID cannot help us discover ( †), but the standard deviation can. As shown in Fig. 1, while F and G have the same average PIDs, the standard deviation of their PID's are quite different: for F, the PIDs of its sequence pairs are 0.5 (1st and 2nd sequences), 0.5 (2nd and 3rd), and 0 (1st and 3rd), and their standard deviation is significantly greater than 0, and for G, the PIDs of its sequence pairs are all equal to 0.33 and their standard deviation is 0. This is not surprising because the sequences in G are identical over the two conserved regions, and are totally different elsewhere. In general, if a family has small standard deviation of PID, it may have structure ( †). (ii) When aligning a family G with structure ( †) to some other family F, we should aim at finding alignment that is good mainly over G 's conserved regions, because G 's sequences are quite different elsewhere and even biologists may not know how to align the * * * * M K G T * * * * G T W * * * * * * metasequence :

Fig. 1
Some structure property that we use for handling normally related sequences sequences correctly over there. Furthermore, since G 's sequences are very similar over the conserved regions, having a good alignment (over the conserved regions) for one sequence of G will essentially give us good alignments for all the others. This suggests that when aligning G to F, we may proceed as if we were aligning a single (meta)-sequence to F (or more concretely, assume that G has only one single sequence).
Observation (i) motivates the following strategy to determine whether a family has structure ( †): If the standard deviation of the PIDs of the family is sufficiently small, we bet that it has structure ( †). For ease of reference, we will say that such family has low PID discrepancy, or simply low discrepancy.
Observation (ii) motivates us to try a guide tree construction method different from GLprobs' when handling families with low discrepancy. Note that the method UPGMA [19] is used in GLProbs to construct guide trees. The method iteratively merges clusters of sequences into larger clusters, and the two closest pair of clusters are chosen and merged in each iteration. The distance d k between two clusters C and C k , where C k is obtained after merging the clusters C i and C j , is In this paper, we try another guide tree construction method for families with low discrepancy; we will use the WPGMA method [19], which is the same as the UPGMA method, but it uses the following definition of distance: Note that (2) is equal to (1) when |C i | = |C j | = 1, or when both C i and C j can be regarded as containing only one single sequence (meta-sequence), as suggested by Observation (ii) for families with low discrepancy.
For distantly related sequences, they are only similar at some local domains or motifs, and these homologous regions may be rather small and are hidden in some long divergent regions. This causes troubles for the progressive alignment method, which is based on global pairwise alignments to merge and align iteratively clusters of sequences together to construct the MSA, and the order of merging depends solely on the guide tree. By insisting global alignments for inputs that have only local similarity, the progressive method may introduce, even in the early stage of execution, many mis-aligned columns and other mistakes, and these early mistakes cannot be corrected and may be propagated [20] and create more mistakes. To improve the alignment quality for distantly related sequences, we forget about the progressive methods and instead, we try non-progressive ones.
There exist many non-progressive MSA methods. For example, the non-progressive sequence annealing technique described in [21,22] combines successively confident alignable regions to build up the multiple alignment; the most similar segments (even in small size) will be aligned first in order to preserve those conserved motifs or domains.
We use this sequence annealing technique to handle input of distantly related sequences. Recall that in Stage 1, we have used the adaptive method to determine substitution scores. During the process, we have also found, for every pair of sequences x and y in the family, and every 1 ≤ i ≤ |x| and 1 ≤ j ≤ |y|, the probability Pr(x i , y j ) of aligning the ith character of x and the jth character of y in the best alignment. To construct an MSA for distantly related sequences, we will first sort all the character pairs (x i , y j ) in descending order of Pr(x i , y j ). Then, starting from the first character pair in the sort list, which has the highest probability of being aligned at the same column, we follow the character pairs in the list and try to insert each pair to the alignment (or more precisely, make the two characters in the pair aligned at the same column) one by one. However, we will actually make the insertion only if the alignment is still consistent after the insertion.
For checking of consistency, we will maintain a collection of "same-column" sets, which contains all the characters that we have determined that they should be aligned at the same column. We will keep track of these sets using a graph, in which its nodes are the sets, and for any two same-column sets S and S , we have a direct edge (S, S ) in the graph if the column for S must precede that of S in the alignment (e.g., when S contains the 10th character of sequence x and S contains the 20th of x). When we insert a pair (x i , y j ) to the alignment, we will update the graph by either (i) introducing a new same-column set (when both x i and y j are not currently in any same-column set), or (ii) adding either x i or y j in some existing same-column set (e.g., if y j is already in some S, then we need to add x i to S after inserting (x i , y j )), or (iii) merge two same-column sets (e.g., if x i is already in S and y j in S , then after inserting (x i , y j ) we need to merge S and S together).
We also need to update the edge set of the graph to reflect the changes. Note that we will not actually make the insertion unless the updated graph is still acyclic, which means that the column constraints are still consistency. When we have finished processing all the character pairs in the sorted list, we topological-sort the graph to get a skeleton of the MSA. We obtain the final MSA by adding to it those characters not in the skeleton. See [16,[21][22][23] for more details.

Construct better guide trees for normally related sequences
PnpProbs uses an adaptive approach to generating guide trees for normally related sequences. As mentioned in Section "Ideas for improving GLProbs", we have two methods, the UPGMA and the WPGMA method to construct guide trees. To study which methods is better, we have modified GLProbs such that it uses the WPGMA method to construct guide trees. For ease of reference, we use GLProbs-UPGMA to refer the original GLProbs, and GLProb-WPGMA to refer the modified one. We used both tools to align the normally related families in SABmark, OXBench and BALiBASE, and compute the TC scores of the resulting alignments, which is one of the most commonly used performance measure for evaluating multiple sequence alignments; the higher the scores, the better. Figure 2 shows the accumulated differences of their TC scores. To explain the figure, let us denote by σ G (PID) the standard deviation of the PIDs over all pairs of sequences in family G (we will drop the subscript G when there is no confusion). The curve in the figure is constructed as follows.  • We first classify the input families according to their σ (PID)s, and for each group i, i.e., the group with σ (PID)= i, we compute the average TC scores TC WPGMA and TC UPGMA over the alignments returned by GLProbs-WPGMA and GLProbs-UPGMA for the families in this group, respectively. Then, we i.e., the accumulated differences up to the group with σ (PID)= h is k.
Note that if the curve is increasing at (h, k), we have h > 0 and GLProbs-WPGMA is doing better than GLProb-UPGMA. As shown in Fig. 2, the accumulated differences is mainly increasing until σ (PID) reaches around 11.5 %, and hence GLProbs-WPGMA is doing better up to this point. Afterwards, the curve is decreasing, which means GLProbs-UPGMA is doing better. Therefore, as default, PnpProbs decides that a family has low discrepancy if its σ (PID) is smaller than 11.5 %, and uses the WPGMA method to construct its guide tree.

Non-progressive alignment for distantly related sequences
Recall that PnpProbs uses a non-progressive method to generate the MSA for distantly related sequences. To get more insight into the relative strength of the progressive and non-progressive methods, we have compared the performance of GLProbs with that of another MSA tool, PicXAA, which uses the nonprogessive sequence annealing method. To make the comparison more meaningful, we have modified the first stage of PicXAA so that it uses the same adaptive approach as GLProbs for generating substitution scores. We call the modified tool PicXAA-AD. Figure 3 shows the accumulated TC score difference between PicXAA-AD and GLProbs for aligning families in the three benchmark databases, namely SABmark, OXBench and BAliBASE. Note that the accumulated differences is increasing until the point around 18 %, and then is decreasing afterwards. This means that the nonprogressive tool PicXAA-AD is doing better when the similarity of the input is less than 18 %, and the progressive tool GLProbs is doing better for the other inputs.

The algorithm of PnpProbs
Given an input family of sequences, PnpProbs constructs its MSA as follows.   If avg(PID) ≥ 18 % (this is the normally related sequences case.) (a) Compute the distance matrix for every pair sequences. (b) Determine the guide tree construction method based on some threshold τ on the standard deviation σ (PID) of the PIDs, whose default value is 11.5 % as suggested by our study in Section "Construct better guide trees for normally related sequences": If σ (PID) < τ, use the WPGMA method to construct the guide tree; otherwise, use the UPGMA method (c) Based on the constructed guide tree, perform the profile-to-profile alignments to construct the MSA.
5. Refine the MSA given in the previous step as follows: we iteratively divide the MSA into two groups by randomly assign each sequence one of them, and we re-align these two groups using standard profile-profile alignment method to see if any improvement can be made. We stop when either (i) we have made 2N iterations and still cannot make any improvement, or (ii) we have made 4N iterations. Here, N is the number of input sequences.

Results
To evaluate the performance of PnpProbs, we have compared it with thirteen other leading multiple sequence alignment tools on three popular benchmark databases. PnpProbs has the best performance in almost all cases, and it achieves significant improvements over the other tools on distantly related sequences. We have also studied its practicability by using it for phylogenetic analysis.

Benchmark comparison
We have compared PnpProbs with the following multiple sequence alignment tools, ten of them use the progressive method: To measure the accuracy of their alignments, we used the sum-of-pairs score (SP) and the total-column score (TC), which were commonly used in previous studies [2,10,11,13,15]. Table 1 compares the performance of the tools on OXBench. It is divided into four categories according to the similarities of the input families. For example, the category "ALL(0-100 %)" show the average SP and TC scores over all the input families used in the test, and the category "(0-20 %)" are for families with similarities between 0 and 20 %. Notice that PnpProbs achieved the overall highest SP and TC scores, and has big improvement for distantly The table shows the average SP and average TC score (multiplied by 100). The best and second best results in each column are marked in bold and with *, respectively. The last column shows the running time using a single CPU thread. Note that we use default parameters for all tools  Fig. 4 shows that for the category (0-20 %), PnpProbs achieved improvements over ClustalW by 36.5 %, over PicXAA by 12.9 % and over GLProbs by 8.4 %. Table 2 shows the average SP and TC scores for SABmark 1.65. The Twilight Zone contains sequences with less than or equal to 25 % similarity, and the Superfamily contains sequences with similarity mostly between 20-50 % similarity. Table 3 shows the average SP and TC scores for BAliBASE 3.0. RV11 contains distantly related sequences (with less than 20 % similarity) and RV12 contains medium to divergent sequences with similarities from 20 to 40 %. For these two benchmark databases, PnpProbs achieved the highest scores in most tests. Again, its improvement was more significant for distantly related sequences, i.e., Twilight Zone and RV11 subsets.
For the efficiency of PnpProbs, we note from the last column of Tables 1, 2 and 3 that even using a single CPU thread, the running time of PnpProbs is comparable to most other tools. Moreover, it is straightforward to "parallelize" Step 1, 2, 3, 4a, 4i, 4iii and 5 of the algorithm of PnpProbs, and thus we can speedup PnpProbs' execution easily by using multiple-cores CPUs. Figure 5 shows PnpProbs' speed when running on a platform of six i7-3930k dual-cores with 64G RAM for inputs with different number of sequences. We note that PnpProbs takes an average of half an hour to align 1000 sequences.

Phylogenetic analysis
To compare the practicability of PnpProbs with other existing tools, we have used it, as well as five other MSA tools, namely GLProbs, MSAProbs, PicXAA, MUSCLE and ClustalW, to construct phylogenetic trees. Given a set of sequences, we first used the six MSA tools to construct six MSAs, and used them as input to the Maximum Parsimony method [24] to infer six hypothesized phylogenetic trees. Then, for each of these hypothesized trees, we calculated the Robinson-Foulds(RF) distance [25] between the tree and the model phylogenetic tree; the smaller the distance, the closer the two trees, and hence the better the corresponding MSA. Our tests used input sequences  chosen from two benchmark databases, namely Yule-Harding tree simulated data [17] and SABmark empirical data [5]. Figure 6 shows the results for inputs chosen from the Yule-Harding tree simulated database, which contains, for every family of sequences, a reference tree and a reference alignment for the family. We used the provided reference trees as the model trees to calculate the RF distance. We also use the reference alignment given in the database to construct a phylogenetic tree, and we refer this tree as RefAln.

Simulated data
Note that the Yule-Harding database is divided into four categories according to the simulated branch length diameter, and the larger the branch length diameter, the more divergent the phylogeny. (i.e., 1 minus the RF distance) as the score for measuring the similarity of two trees. Note from the figure that in most cases, the hypothesized trees derived from PnpProbs' alignments achieve scores higher than that of the other tools, and we can argue that the alignments of PnpProbs are better. Furthermore, as shown in Fig. 7, the RF distance differences between RefAln and the other hypothesized trees become larger when the phylogenies are more divergent (i.e., with larger branch length). However, the differences for PnpProbs increase mildly and are smaller than those of the five other MSA tools. Figure 8 shows the result for inputs obtained from the SABmark empirical database, which contains, for every family of sequences, a reference alignment. However, the database has no reference trees; thus we used RefAln as the model tree to compute the RF distance. We note that the hypothesized trees derived from PnpProbs' alignments have the highest accuracy.

Discussion
Our MSA tool PnpProbs aims at combining the strength of progressive and non-progressive methods for multiple sequence alignment; it uses progressive method for normally related sequences, and uses non-progressive method for distantly related ones. In [1], we proposed to use the average percent identity to estimate the similarity of a family of sequences, and in this paper, we proposed to use the standard deviation of the percentage identity to estimate the discrepancy of a sequence family. For normally related sequences, PnpProbs uses different methods to construct guide trees depending on the discrepancy of the family. Our experimental results showed that PnpProbs has the best TC scores in all but one test. We have also evaluated PnpProbs' practicability, and our results suggested that PnpProbs will be a useful tool for downstream phylogenetic analysis. For possible future research direction, we note that most of the MSA tools try a certain range of sizes of components to assemble multiple sequence alignment. For example, the progressive alignment method uses big components of sequence profiles, and the non-progressive sequence annealing technique uses small components, e.g., alignable columns or residue pairs. A natural research direction is to consider multiple sizes of decomposed components in one algorithm to build up the MSA such that families of sequences with long conserved regions apply large components and those with small conserved patterns use small components.