| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
© 2005 Society of Systematic Biologists
Weighted Least-Squares Likelihood Ratio Test for Branch Testing in Phylogenies Reconstructed from Distance Measures
Edited by Junhyong Kim: Associate Editor Chris Simon Editor
1 Institut Cavanilles de Biodiversitat i Biología Evolutiva, Universitat de València, Edifici d' Instituts de Paterna Apartat 2085 46071 València, Spain E-mail: rafael.sanjuan{at}uv.es borys.wrobel{at}uv.es (R.S.)
2 Department of Marine Genetics and Biotechnology, Institute of Oceanology, Polish Academy of Sciences Gdynia, Poland
 |
Abstract |
|---|
 Top Abstract METHODRESULTSDISCUSSIONREFERENCES |
|---|
A variety of analytical methods is available for branch testing in distance-based phylogenies. However, these methods are rarely used, possibly because the estimation of some of their statistics, especially the covariances, is not always feasible. We show that these difficulties can be overcome if some simplifying assumptions are made, namely distance independence. The weighted least-squares likelihood ratio test (WLS-LRT) we propose is easy to perform, using only the distances and some of their associated variances. If no variances are known, the use of the Felsenstein F-test, also based on weighted least squares, is discussed. Using simulated data and a data set of 43 mammalian mitochondrial sequences we demonstrate that the WLS-LRT performs as well as the generalized least-squares test, and indeed better for a large number of taxa data set. We thus show that the assumption of independence does not negatively affect the reliability or the accuracy of the least-squares approach. The results of the WLS-LRT are no worse than the results of the bootstrap methods, such as the Felsenstein bootstrap selection probability test and the Dopazo test. We also show that WLS-LRT can be applied in instances where other analytical methods are inappropriate. This point is illustrated by analyzing the relationships between human immunodeficiency virus type 1 (HIV-1) sequences isolated from various organs of different individuals.
Keywords: Analytical branch test; bootstrap; distance; generalized least-squares; maximum likelihood; phylogeny; weighted least-squares
Received August 12, 2003; Revised December 14, 2003; Accepted November 10, 2004
Bootstrapping (Efron, 1982; Felsenstein, 1985) has become the most frequently used approach to assess statistical significance in the field of phylogenetic reconstruction. Felsenstein (1985) bootstrap selection probability (BP) is commonly used to assign confidence levels for clades, and an internal branch test based on bootstrap has been proposed by Dopazo (1994). In the BP test, the frequency of a given clade among the trees constructed using the bootstrap pseudoreplicates is calculated. The Dopazo method estimates the distribution of the length of a given branch based on bootstrapping and then tests if the distribution includes the zero value. Resampling methods have the advantage of avoiding the assumption of underlying evolutionary models and are not limited to any number of taxa, as are some of the analytical tests described below. However, the statistical properties of the bootstrap are unclear and its ability to assess tree accuracy and repeatability has been challenged (Hillis and Bull, 1993; Felsenstein and Kishino, 1993; Berry and Gascuel, 1996; Efron et al., 1996; Li and Zharkikh, 1995).
Since the incorporation of bootstrapping to phylogenetics, alternative analytical methods have been proposed. For example, Nei et al. (1985) proposed a method that can be applied to UPGMA trees (Sneath and Sokal, 1973) obtained from sequence, electrophoretic, or restriction data, and Li (1989) and Tajima (1992) suggested two tests that are not specific for any particular tree-building method. In all three cases, the branch length variances are estimated from distance variances and covariances, which are themselves obtained from the distance matrix. However, the calculation of the covariances became too complicated for more than four or five taxa. Takahata and Tajima (1991) intended to solve this problem by calculating the maximum and minimum covariances and thus deducing an interval for branch length variances.
In the maximum likelihood framework, the likelihood ratio test can be used for comparing nested models, which makes it suitable for branch testing. In this framework, there are methods that can be used to assign clade support, such as the quartet puzzling algorithm (Strimmer and von Haesler, 1996) and the Bayesian (Markov Chain Monte Carlo) methods (Larget and Simon, 1999). For non-nested topologies there are also maximum likelihood methods that involve bootstrapping (parametric or nonparametric): the Kishino-Hasegawa test (Kishino and Hasegawa, 1989), modified later to take into account test multiplicity (Shimodaira and Hasegawa, 1999); the Swofford-Olsen-Waddel-Hillis test (Swofford et al., 1996; Goldman et al., 2000); the Approximately Unbiased Test (Shimodaira, 2002); and the Expected Likelihood Test (Strimmer and Rambaut, 2002).
In the least-squares (LS) methods (Cavalli-Sforza and Edwards, 1967; Fitch and Margoliash, 1967), the residual sum of squares between observed and expected distances is minimized. The generalized sum of squares is of the form 
i,j,k,lwij,kl(dij– eij)(dkl– ekl) where dijand eij stand, respectively, for the observed and estimated distances, and wij,kl for the weights assigned to each element of the sum. If distances are assumed to be independent, then this sum becomes i,jwij(dij – eij)2, which is a weighted sum of squares, and if all the weights are set equal, then it becomes the ordinary sum of squares. Given a topology, the generalized least-squares (GLS) statistic can be calculated analytically (Bulmer, 1991). LS have been incorporated into tree-building methods such as the Fitch-Margoliash (FM; Fitch and Margoliash, 1967) and minimum evolution (ME; Rzhetsky and Nei, 1993). Although not explicitly an LS method, the widely used neighbor-joining (NJ) algorithm (Saitou and Nei, 1987) is related to LS (Felsenstein, 1997) and usually finds a tree very close to the ME solution (Rzhetsky and Nei, 1992a). The recently proposed balanced minimum evolution (BME) method (Desper and Gascuel, 2002) is a special case of the LS approach. BME is both fast and accurate, and the fact that it does not output negative branch lengths is an additional advantage over the previous methods (Desper and Gascuel, 2004). It is worth noting that because a non-GLS criterion is used in FM, ME, and BME, all these reconstruction methods assume independence between distances (Swofford, 2001; Felsenstein, 1997). The GLS approach is rarely used because of the difficulties in estimating covariances, both theoretical and computational. Using the GLS for branch testing, two trees can be compared (Bulmer, 1991), one of which being a "contraction" of the other (that is, the length of the branch being tested is set to zero). Under the null hypothesis, both trees are equally likely, and thus the difference of the residual sum of squares for the two trees should follow the chi-square distribution with degrees of freedom equal to the difference in the number of branch lengths being estimated, that is, one. However, the Bulmer test is applicable only when Tajima and Nei (1984) nucleotide distances are used. More recently, Susko (2003) has extended the GLS applicability showing that, provided the distances satisfy the maximum likelihood criterion, their distribution is approximately multivariate normal, which allows their variances and covariances to be estimated. The GLS statistic follows then the chi-square distribution; given the data (from which the distance matrix is calculated), a probability can be assigned to each topology to allow the sorting of a set of competing topologies according to their probability and to establish confidence regions for topologies.
Distance-based analytical tests use only the distance matrix, whereas resampling methods such as the bootstrap use the primary data (the matrix of characters). As a consequence, analytical methods remain useful even if there is no access to the original data, when the original data are the distances themselves (e.g., DNA hybridization distances), when the distances are indexes that were not directly drawn from characters (e.g., Fst index), or when the distances are averages across different data sets. Resampling is difficult or impossible in such situations. Analytical methods require the estimation of the variances of the distances and the covariances between them and this has limited their application in two ways. First, the covariances are rarely part of the original data and are difficult to estimate. The solution to this problem, provided recently by Susko (2003), allows for the estimation of a large number of covariances, but estimating a large number of parameters can have an adverse effect on the accuracy. Second, analytical methods proposed thus far are based on models of molecular evolution, and their applicability is limited to sequence data.
The purpose of this work is to show that when some general simplifications are made, analytical tests such as the weighted least-squares (WLS) test we propose or the Felsenstein (1983)F-test do not lose their accuracy or reliability, and that they give results nearly as good as the results of the commonly used resampling methods. These simplified analytical tests are not restricted for sequence data and do not require estimating a large number of parameters.
|
METHOD |
|---|
|
TopAbstractMETHODRESULTSDISCUSSIONREFERENCES |
|---|
Weighted Least-Squares and Maximum Likelihood
Consider an unrooted tree B for which 2T – 3 branch lengths have been estimated by minimizing some function f (the desirable properties of f will be discussed further below). Under hypothesis B, each observed distance dij has a value drawn randomly from a distribution with expectation bij, the patristic distance. Assuming that this distribution is normal, N(bij,
ij), and considering ij constant (ij = , 
i, j), the likelihood of bij is P[bij()] = N(dij| bij, ), that is,
|
| (0001) |
If all dij were independent, after including all pair-wise distances, the log-likelihood of the tree B would be
|
| (0002A) |
. If the assumption of the constant variance does not hold, then Equation 2a becomes
|
| (0002B) |
and thus the maximum likelihood solution is equivalent to the WLS solution.
Choosing the Function f Related to the Maximum Likelihood
A test for topologies using the GLS statistic was first proposed by Bulmer (1991) and was then made more general by Susko (2003). The only underlying assumption necessary for the GLS calculation is distance normality. In the WLS method we propose, the distances are assumed to be normally distributed and independent (i.e., the covariances are ignored). The variances can be estimated analytically, for example, by the sample average method (Susko, 2003), or by nonparametric bootstrapping (Felsenstein, 1985; Susko, 2003). Alternatively, they can be part of the original data. In general, the relationship between distances (bij) and their associated variances (ij) can be expressed as 2ij = n = 0
2(n)bijn. If 2(n) = 0, n 
0, then 2ij = 2(0) and the function f, which is minimized to obtain the maximum likelihood solution, is the ordinary sum of squares. If 2(n) = 0, n p, then
|
| (0003) |
and f is determined by the weighted sum of squares, where p is called the power of this sum. As a reasonable approximation, terms bij in the denominator of the weighted sums of squares can be replaced by the observed distances (dij); this overcomes the difficulties in minimizing f when unknown parameters (i.e., the branch lengths) appear both in the numerator and the denominator (Felsenstein, 1983). If all or some of the variances associated with the distances are known, it is possible to find the optimal power value by fitting the equation
|
| (0004) |
to the data. This allows reducing the T(T–1)/2 variances to only two parameters and performing branch testing when only a subset of the variances is known, which might be the case when the variances are part of the original data (see the analysis of human immunodeficiency virus type 1 (HIV-1) organ compartmentalization below). Finally, the power value p estimated using Equation eq2 can be used to calculate the sum of squares as
|
| (0005) |
Weighted Least-Squares and the Interior Branch Testing
To test a given branch of the optimal tree B, it is possible to collapse this branch (make its length equal to zero) and to recalculate the length of the 2T – 4 remaining branches in this new, "contracted" or "collapsed" tree (C), as well as the associated sum of squares to obtain the tree likelihood (PC) as shown in Equation 2a, 2b. Then, taking the difference between the log-likelihood of trees B and C, we obtain:
|
| (0006) |
where SSB and SSC are the sums of squares for the optimal and collapsed trees, respectively, as defined in Equation eq3. Under the assumptions of distance normality and independency, 2ln [PB (,T)/PC (,T)] can be approximated by the chi-square distribution (Sokal and Rohlf, 2000). Hypothesis B employs 2T – 3 parameters, whereas C employs 2T – 4 parameters. Hence, the weighted least-squares likelihood ratio test (WLS-LRT) we propose has one degree of freedom and is a particular form of the GLS approach (Bulmer, 1991; Susko, 2003).
If we ignore all variances or if we do not trust the reliability of the estimations of p and 2(p), it is possible to use the Felsenstein (1983)F-test to compare two alternative trees. Under the null hypothesis that the topology B and C explain the data equally well, the fraction 
will distribute as an F(dfC – dfB, df}B) statistic, provided that distances are normally distributed and independent. The degrees of freedom associated with each topology are dfB = T(T – 5)/2 + 3 and dfC = T(T – 5)/2 + 4. In other words,
|
| (0007) |
The square root of this statistic will distribute as a one-tailed t-statistic with T(T–5)/2 + 3 degrees of freedom. As no estimation of 2(p) is needed, it is possible to test branches even when none of the variances associated with the distances are known. In such a case, however, it is not possible to select the optimal power value from the data.
Correcting for Multiple Tests
When testing a phylogeny implies testing a series of branches belonging to the same evolutionary hypothesis (in the extreme case, all the branches in the tree), a multiple test correction should be undertaken. As it is unclear how the significance threshold is affected in this situation, one can use Bonferroni or Dunn-Sidák (Sokal and Rholf, 2000) general corrections or the correction suggested by Felsenstein (1985) for phylogenies. Given that in phylogenetic trees, the multiple tests have an unknown degree of statistical dependence, Dunn-Sidák correction is the most appropriate because it simply assumes complete independency, leading in the worst situation to conservative results. Dunn-Sidák corrected significance threshold is 
' = 1 – (1 – )1/n, where n is the number of tests performed.
Implementation
The WLS-LRT, as well as the related Felsenstein (1983)F-test, have been implemented in a C program (WeightLESS) derived from the source code of the PHYLIP package (Felsenstein, 1993), available for download at www.iopan.gda.pl/~wrobel, where detailed instructions about its use can be found. The program provides an interface similar to that of PHYLIP package and takes as input the bifurcating tree in which branches are to be tested, a distance matrix, a power value, and a 2(p) value. WeightLESS calculates the significance values according to Dunn-
idák correction for the case in which all internal branches are going to be tested. If the user is interested in testing a lower number of branches, this number can also be introduced and, WeightLESS calculates the corresponding Dunn-idák significance threshold. In the output, P-values for every branch are listed, indicating which of them are significant. The program also outputs a tree in which interior branches are labelled with the P-values.
As described above, testing the branches involves collapsing them and recalculating the remaining branch lengths by the WLS method. In this second stage, we have used the algorithm implemented in FITCH (Felsenstein, 1993), which requires O(n2) time, where n is the number of taxa. Since this is repeated for every collapsed branch, our method requires O(n3) time.
WeightLESS allows also estimating parameters p and 2(p), either if a variance matrix or if several distance matrices are provided. In the latter case, the computation time at the parameter estimation stage depends linearly on the number of matrices and quadratically on the number of taxa. In practice, this may be quite lengthy if the number of matrices is high, as recommended if the matrices are obtained by bootstrapping. However, this step needs to be done only once for a given data set.
|
RESULTS |
|---|
|
TopAbstractMETHODRESULTSDISCUSSIONREFERENCES |
|---|
Simulations
To assess the accuracy of branch testing in the WLS framework, a four-taxa simulation scheme was employed. The topology was fixed and the length of the interior branch was either 0, 0.025, or 0.05 substitutions per site. For each interior branch length, EVOLVER (part of the PAML package; Yang, 1997) was used to simulate sequence alignments with 100, 300, 1000, 3000, and 15,000 characters for 100 four-taxa phylogenies in which the length of each of the four exterior branches was assigned a random value ranging from 0.03 to 0.07 substitutions per site. This was done because, in perfectly symmetric trees, all pair-wise distances would be identical and therefore no estimation of parameters p and
2(p) would be possible. Felsenstein (1984; Kishino and Hasegawa 1989) model was used (with the transition:transversion parameter set to 2.0). From the 100 simulated data sets, unrooted trees were reconstructed with the FM-modified FITCH method using programs DNADIST and FITCH (belonging to the PHYLIP package; Felsenstein, 1993); negative branch lengths were not allowed for reasons explained elsewhere (Kuhner and Felsenstein, 1994; Susko, 2003). We used the same model of sequence evolution to calculate distances and to simulate the data. To calculate the two parameters needed for the WLS-LRT, 1000 bootstrap pseudoreplicates for each simulated alignment were obtained, and Equation eq2 was fitted by ordinary linear regression (Sokal and Rohlf, 2000). Table 1 shows the results after testing the interior branch in the simulated phylogenies using the WLS-LRT, compared to the results obtained using other methods, such as the Felsenstein (1983)F-test, Susko's implementation of the GLS-LRT (Bulmer, 1991; Susko, 2003), Felsenstein's (1985) bootstrap selection probability test, and the Dopazo (1994) bootstrap internal branch test. A thousand bootstrap pseudoreplicates were performed in the resampling methods. It can be observed that both the WLS-LRT and GLS-LRT showed a very low proportion of false positives (i.e., assigning a significant value to a zero-length branch, that is, type I error). In contrast, type I error for the F-test ranged between 7% and 18%, depending on the length of the interior branch. High incidence of false positives was also observed when the BP threshold was set at 70%, a value commonly used in this test. In these simulations, the threshold needed to be set at 95% to keep the incidence of type I error low. As far as type II error is concerned, when the interior branch length was 0.025 substitutions per site, the sequence length needed to be 1000 nucleotides long to reduce the incidence of false negatives below 10% with all tests except for the F-test. When the interior branch was 0.05, 300 characters were enough for a roughly equal performance.
|
In all simulations the performance of the WLS-LRT and the GLS-LRT was very similar. We can observe in passing that for some simulated data sets with short sequence lengths the GLS statistic could not be calculated using the GLSDNA program (Susko, 2003). This was possibly occurring when the covariance matrix was singular, which has a higher chance of occurring for shorter sequences. We conclude that the assumption of distance independence, which differentiates WLS from GLS, is far from affecting either the reliability or the accuracy of the test. The power was only slightly better (lower type II error) for bootstrap methods than for LS-based tests. However, although the WRT-LRT performed well in the four-taxa problem with simulated data, phylogenetic analyses are usually more complex. The use of a data set of biological sequences taken from a higher number of taxa is therefore in order.
rRNA-tRNA Data
We have used a previously published data set (Stanhope et al., 1998) of 12S rRNA, tRNA-valine, and 16S rRNA sequences from 43 mammalian species, 2086 nucleotides long (TreeBase accession code M623). Information about each species can be retrieved at www.ncbi.nlm.nih.gov/entrez/query.fcgi?db = Taxonomy. Using ModelTest 3.1 (Posada and Crandall, 1998) and PAUP* 4.0b10 for Unix (Swofford, 2002), the substitution model that best explained the data under Akaike Information Criterion (Akaike, 1974) was the general time-reversible (GTR/REV) model with rate variation among sites following a gamma distribution plus an invariant sites category. One thousand bootstrap pseudomatrices were generated with SEQBOOT (part of the PHYLIP package, Felsenstein, 1993) and PAUP*. These pseudomatrices were used to estimate the variances associated with each interspecies distance. The linear regression of the log-variances on the log-distances is shown in Figure 1. Parameters needed for the WLS-LRT were estimated using Equation eq2: power value p = 1.823 and 2(p) = 0.0060. This solution fitted the data significantly better (r2 = 0.981) than the commonly used default values p = 2 (r2 = 0.972, partial F(1,42) = 18.947, P < 0.0001) or p = 1 (r2 = 0.927, partial F(1,42) = 119.368, P < 0.0001). Thus, p = 1.823 was used as the exponent value in subsequent calculations of the weighted sums of squares both at the tree inference and the branch testing levels. Three methods were employed for phylogenetic reconstruction: FITCH, minimum evolution (ME), and neighbor-joining (NJ), implemented in FITCH and NEIGHBOR programs (PHYLIP package). Again, negative branches were not allowed. FITCH, ME, and NJ trees are shown in Figures 2–4, respectively. The minimal SSB corresponded to the FITCH tree (SSB = 3.2925), and the value for the ME tree (SSB = 3.3639) was larger than the value obtained for the NJ tree (SSB = 3.3047). Once the three trees were constructed, the 40 interior branches were tested with WLS-LRT and Felsenstein (1983)F-test. The significance threshold was held at = 0.05 using Dunn-Sidák multiple test correction (0.05' = 1.282 x 10– 3). Bonferroni and Felsenstein (1985) corrections gave very similar 0.05' values (1.250 x 10– 3 and 1.292 x 10– 3, respectively) and did not alter the results presented hereafter. The results of the WLS-LRT were identical to those obtained with the F-test. Finally, setting the power to the suboptimal default value p = 2 did not affect the results of the WLS-LRT.
|
Figures 2–4 allow comparing the performances of the WLS-LRT and the Felsenstein (1985) bootstrap selection probability (BP) test in FITCH, ME, and NJ trees, respectively (1000 bootstrap pseudoreplicates were performed in the BP test). Overall, the significance threshold of the WLS-LRT corresponded to a BP of roughly 60% to 80%, with a slightly higher uncertainty for ME and NJ trees than for FITCH tree. Particularly, there was a branch having only 25.0% and 20.0% bootstrap support in ME and NJ trees, respectively, which was significant according to the WLS-LRT test. Conversely, there were few branches in all three trees with high (75.3% to 93.9%) bootstrap support, which were nonsignificant according to the WLS-LRT. Although the WLS-LRT probability value for these latter branches is generally close to the Dunn-Sidák corrected significance threshold, and therefore these results can be attributed to a lack of power of the WLS-LRT, an alternative explanation can be given for this discrepancy. As mentioned above, looking at Table 1 it is apparent that in a four-taxa topology, the appropriate bootstrap threshold is closer to 95% than to the commonly used 70%. It is thus possible that small subtrees in the 43-taxa phylogeny that were delimited by neighboring nodes with a very strong statistical support (i.e., BP close to 100%) can show this same effect. This can be easily understood, because in a small tree or highly delimited small subtree, the probability of a given branch to occur simply by chance in bootstrap pseudoreplicates is higher than in large trees, where more combinations are possible. Branches with high BP and nonsignificant WLS-LRT values seemed to follow this pattern (Figs. 2–4).
|
|
The other previously proposed bootstrap test for branch testing, the Dopazo test, gave results very similar to the WLS-LRT for the FITCH, ME, and NJ trees of rRNA-tRNA sequences (Table 2). Again, 1000 pseudoreplicates were used. Taking the results of the Dopazo test as a reference, for the FITCH tree the WLS-LRT produced three false negatives (the total number of branches tested was again 40), the P-value in these cases being relatively close to the significance limit (0.00223, 0.00501 and 0.09512), and no false positives. For the ME tree, there were only two putative false negatives (very close to the significance threshold) and no false positives. Analogous results were obtained for the NJ tree. Finally, for the above-mentioned branch, in which the WLS-LRT was significant and the BP was as low as 57.7, 25.0, and 20.0% in FITCH, ME, and NJ trees, Dopazo values were 100, 99.8, and 99.0%, respectively, reflecting discrepancies between the two bootstrap methods and the possibility that the true BP appropriate threshold might be much lower than 70% in this case.
|
The GLS statistic for the rRNA-tRNA distance matrix could not be calculated using the GLSDNA program (Susko, 2003). This is probably due to the fact that, with large number of taxa, some of them closely related, the covariance matrix can become almost singular. For such cases, the GLSDNA_EIG routine not described in the original publication is available from its author's homepage www.mathstat.dal.ca/~tsusko. Moreover, this latter method does not allow for highly complex substitution models such as the GTR/REV. Therefore, we used the suboptimal F84 (Felsenstein, 1993) model. As can be seen in Table 2, the GLS-LRT calculated using GLSDNA_EIG is extremely conservative for the rRNA-tRNA data set. This poor performance compared to the WLS-LRT is not attributable to the inadequacy of the substitution model employed in the GLSDNA_EIG test, since this same behavior was observed when a simple two-parameters model (Kimura, 1980) was used in all tests (data not shown). In this approximate method, the GLS is interpreted as a WLS for a set of linear transformations of the distances. To make the covariance matrix invertible, an arbitrary eigenvalue cutoff has to be selected, and linear transformations for which the variance is lower than the cutoff are ignored. These linear transformations are usually those with a better discriminatory power, and this is why the results are expected to be highly conservative. The greater the eigenvalue cutoff, the lower the statistical power will be. However, with values below 10– 9, no further power improvement was observed, and the results remained highly unsatisfactory (Table 2). We conclude that difficulties in obtaining the covariance matrix make the GLS-LRT hardly applicable for large data sets, and that the performance of the GLS-LRT by the above approximate methods is notably worse than the performance of the WLS method.
Performance of the Weighted Least-Squares Tests When the Number of Known Variances Is Small
So far, the whole set of T(T–1)/2 variances had been taken into account to estimate the power value and 2(p) (Fig. 1). In the rRNA-tRNA gene data set, with 43 species, the distance matrix contains 903 informative values. Obviously, the fewer variances that are known, the less reliable the estimates of the two parameters from Equation eq2 will be. We therefore sought to determine whether using values of the power and 2(p) away from the optimal would affect the reliability of the WLS-LRT. Let us consider the extreme case in which the two parameters were estimated from only two data points. Two distances were randomly picked from the complete rRNA-tRNA distance matrix, the power and 2(p) were estimated using Equation eq2, and the WLS-LRT was performed for the 40 branches. This procedure was repeated 50 times. Please notice that we used the tree built with the optimal power value, as the purpose here was to check the performance of the WLS-LRT and not the reliability of the tree building method itself. The performance of the WLS-LRT in this situation remained good: taking the results obtained with the full set of variances as a reference, type I error was 1.10 ± 0.23% and type II error was 1.55 ± 0.55% for the FITCH tree. Analogous results were obtained for ME and NJ trees, the percent of false negatives and false positives remaining always below 2.5%. In contrast, the estimation of the power value was highly inaccurate, ranging from 0.232 to 3.389, with a coefficient of variation as high as 29%, and the same is valid for 2(p), for which the coefficient of variation reached 113%. However, the power value and 2(p) were highly correlated (r2 = 0.989) in such a way that the lack of precision in their estimation had little impact on the WLS-LRT statistic. This is why the test performed well even when only two data points were used to estimate parameters in Equation eq2, and also why the WLS-LRT gave the same results when the power value p = 2 was used instead of the optimal p = 1.823.
In the still more extreme case in which no variances are known, no estimation of the power and 2(p) is possible. The Felsenstein (1983)F-test can be used to avoid the estimation of 2(p), but the power value is still needed to calculate the weighted sum of squares. It is interesting, therefore, to study the performance of the F-test when the suboptimal power values are used. Table 3 shows the results of applying this test with power values p = {0, 1, 2, 3, 4, 5}, plus the optimal p = 1.823. As mentioned before, the results of the F-test were the same as the results of the WLS-LRT when the optimal power value was used. However, Table 3 shows that using a lower power increased the incidence of false negatives. Therefore, this test might not be performed with a power value that might be far from the optimal. It has been established that this optimal value is close to p = 2 for pair-wise DNA distance data (Kuhner and Felsenstein, 1994), but this might not be the case for other kinds of data.
|
Analysis of the HIV-1 Organ Compartmentalization
Both in the simulations and in the analysis of the rRNA-tRNA data, the variances of the distances were estimated by bootstrapping. However, other methods to estimate the variances allow the use of the WLS-LRT. We would like to illustrate this point with an example concerning human immunodeficiency virus type 1 (HIV-1) organ compartmentalization, an important issue in virology (McGrath et al., 2001). This question was thus far addressed by making phylogenetic reconstructions using sequences data from different organs of the same individual, in order to assess whether viral isolates from the particular organs group together (Ball et al., 1994). Such an approach can be applied only for viral sequences isolated from the organs of a single individual, which makes it difficult to draw general conclusions. To solve this difficulty, the studies should go beyond the intra-patient level, but then the above approach becomes useless because sequences from different hosts group monophyletically. We sought to solve this problem by constructing a distance matrix in which intra-patient distances between pairs of organs were averaged over several patients.
We constructed a data set from previously published (Keys et al., 1993; Ball et al., 1994; Korber et al., 1994; van't Wout et al., 1998; Gatanaga et al., 1999; Morris et al., 1999; Gorry et al., 2001; Wang et al., 2001; Collins et al., 2002) sequences of the V3 hypervariable region of the HIV-1 isolated from various organs (bone marrow, brain, cerebrospinal fluid, kidney, liver, lung, lymphatic nodes, and spleen) of 22 individuals. Only for one patient were sequences of the virus isolated from all eight organs available (further information is given in the Appendix). Using ModelTest 3.1 (Posada and Crandall, 1998), Tamura-Nei model with rate variation among sites following a gamma distribution plus an invariant sites category was appropriate according to the Akaike Information Criterion (Akaike, 1974). Consequently, this model was used to estimate the distance between each pair of sequences isolated from different organs of the same individual, using MEGA 3 (Kumar et al. 2004); 71 distances were obtained in total. Then, the distances for each pair of organs were averaged over all patients, and thus an 8 x 7 matrix was constructed, which was used to obtain the FITCH tree shown in Figure 3. For 9 out of the 28 pairs of organs, the distances were averages of three or more (up to 16) values; for the remaining 17 pairs, only one distance was available.
|
|
Previous studies (van't Wout et al., 1998) pointed out that the nervous tissue harbors specific HIV-1 variants. This hypothesis can be tested by determining whether the branch grouping cerebrospinal and brain is statistically significant. One could manage to bootstrap the aligned sequence data set, but this would result only in the estimation of intra-patient variability, and inter-patient variability would be ignored. In contrast, the WLS-LRT does not have this limitation. The variances calculated using the 9 pairs of organs for which three or more distances were available were regressed against their corresponding average distances using Equation eq2, giving an estimation of the power p = 1.766 and
2(p) = 0.8589. The sum of squares for the FITCH tree was SSB = 7.6457, and the SSC = 16.1897. The WLS-LRT showed that this branch was significant (
12 = 9.948, P = 0.0016). To sum up, the WLS-LRT approach allowed for the simultaneous analysis of sequences isolated from 22 individuals, taking into account the inter-patient variability, thus providing a solid support to the hypothesis that nervous system harbors divergent HIV-1 variants.
|
DISCUSSION |
|---|
|
TopAbstractMETHODRESULTSDISCUSSIONREFERENCES |
|---|
The first step of any distance method is the reduction of the information contained in the primary data (characters) to a distance matrix. As a necessary consequence, using distance data for phylogenetic reconstruction unavoidably leads to a reduced statistical power compared to character-based methods such as the ML. Additionally, ML inference is generally more robust than classical inference based on mean distances (Huelsenbeck, 1995). However, distance-based methods remain commonly used for phylogenetic reconstruction because they are computationally efficient; they have been shown to be consistent, as opposed to maximum parsimony; and their degree of accuracy is usually close to ML (Kuhner and Felsenstein, 1994). Once the information of the primary data has been reduced to a distance matrix for phylogenetic reconstruction, it should be optimal to use analytical tests for assessing the statistical support of this distance-based reconstruction, for which it is not necessary to go back to the character matrix, as opposed to the bootstrap. Eliminating this step has obvious computational advantages. In this work, we have focused on the LS approach, a well-established criterion for tree inference (Fitch and Margoliash, 1967; Bulmer, 1991; Swofford, 2001; Felsenstein, 1997; Desper and Gascuel, 2004), as well as for assessing the statistical support of phylogenies (Bulmer, 1991; Susko, 2003). We have pointed out that the use of the GLS is restricted by the difficulty of estimating a large set of covariances. Previous simulation studies (Kuhner and Felsenstein, 1994; Desper and Gascuel, 2004) have shown that WLS provides a highly accurate and consistent approach for tree reconstruction. There is no reason why statistical tests based on the same criterion, such as the WLS-LRT presented in this work, should be inaccurate or biased. In the simulations done here, the WLS-LRT was highly reliable in terms of type I error and had a reasonably low type II error. In contrast, the GLS test showed a poor sensitivity in the rRNA-tRNA example. This loss of statistical power is attributable to the increased number of taxa, which makes it difficult to handle the covariance matrix (Susko, 2003), thus imposing the use of an approximate version of the GLS method, whose performance was notably worse than that of the WLS method.
Both GLS and WLS are based on the assumption that distances are normally distributed. In the case of nucleotide or amino acid sequences, it was traditionally believed that this assumption is not fulfilled (Takahata and Kimura, 1981). It has been recently shown, however, that, provided the distances are maximum likelihood estimates, their distribution is approximately multivariate normal (Susko, 2003). Moreover, the WLS-LRT is applicable to a wider range of distance measures, not necessarily derived from sequence data, and this is why it seems reasonable to assume normality in general. The assumption that makes the WLS different from the GLS is the independence between distances, a priori unfulfilled in phylogenetic reconstructions because of the common evolutionary history. However, our results suggest that this does not affect the performance of the test. The WLS-LRT is expected to perform well if the model ij2 = (p)2dijp provides a satisfactory description of the relationship between the distances and their variances. This seems to be the case at least for DNA sequence data (Fig. 1). A previous study (Kuhner and Felsenstein, 1994), using simulated nucleotide data, established p = 2 (the original FITCH) as the integer value for which the best results could be obtained regarding topological accuracy. Here, we provide a criterion for selecting the optimal power value, which does not need to be any of the commonly used integers 0, 1, and 2. For the rRNA-tRNA data, the model fitted the data accurately, and we were able to apply the WLS-LRT under optimal conditions. However, the conditions might not be optimal for two reasons. First, when the variances are part of the primary data rather than estimated, it is possible that only some would be known. In this scenario, estimates of the power value and 2(p) can be much less reliable, but this will have a negligible effect on the WLS statistic provided that Equation eq2 fits the data accurately. The second scenario, in which the model fit is poor, remains unexplored and would probably lead to weaker performance of the test. In some cases, if few or no variances are available, it may be necessary to use the Felsenstein F-test, for which the estimation of 2(p) is not necessary. We show that the results of this test need to be taken with caution. In our simulations, high incidence of type I error was observed. Moreover, the rRNA-tRNA data illustrate that a priori knowledge of the power value is necessary.
Contrary to the simulations, the rRNA-tRNA data set consisted of a number of sequences high enough to appreciate differences between FITCH, ME, and NJ trees: for the FITCH tree, the WLS-LRT and the Felsenstein F-test gave slightly closer results to the bootstrap methods than for the ME and the NJ trees. The FITCH method is directly based on the WLS criterion and thus gives the maximum likelihood solution under the assumptions of distance normality and independence. These are precisely the two assumptions that the WLS-LRT makes, and this is why FITCH trees are better suited to WLS-LRT. Some authors claimed before that the FM method is less efficient than the ME or NJ methods in finding the correct topology (Tateno et al., 1982; Saitou and Imanishi, 1989; Rzhetsky and Nei, 1992b); but more recently reported extensive simulation studies using FITCH modification indicate the contrary (Kuhner and Felsentein, 1994; Desper and Gascuel, 2002).
The non-parametric bootstrap does not require any underlying evolution model and can be used for many kinds of data, regardless of the number of taxa considered or the reconstruction method employed. It has been shown previously that the bootstrap consensus trees have a higher chance of reproducing the correct topology than unique optimal trees (Sitnikova et al., 1995; Berry and Gascuel, 1996). Despite these advantages, the statistical basis of the bootstrap is unclear (Hillis and Bull, 1993; Felsentein and Kishino, 1993; Li and Zharkikh, 1995; Berry and Gascuel, 1996; Efron et al., 1996). The results of Hillis and Bull (1993) are misconstrued by some authors to suggest that the appropriate threshold for clade selection corresponds to a BP of 70%. However, this value was established for maximum parsimony and depends on many factors, which include the tree symmetry, the divergence between taxa, and the constancy of the substitution rate among branches (Hillis and Bull, 1993). In our simulations, a BP of 70% did not lead to reliable results, and the threshold had to be pushed up to 95% in order keep the incidence of type I error low. In the analysis of the rRNA-tRNA data, with 43 rather than 4 sequences, a threshold roughly around 60–80% seemed more appropriate. However, there were some strong differences between the BP and the WLS-LRT, which did not seem to be randomly distributed along the topology: for those nodes that were delimited by neighboring highly significant branches, BP up to 95% were nonsignificant according to the WLS-LRT, and conversely, it seemed that in those nodes that were not delimited by neighboring solid branches, BP down to 20% were significant according to WLS-LRT. Although not fully proven here, the BP appropriate threshold might be down-biased in the former situation, whereas the BP could be up-biased in the latter. These results illustrate how complicated it is to establish a reliable BP threshold for clade selection, and argues in favor of using interior branch tests.
|
ACKNOWLEDGEMENTS |
|---|
We would like to thank Miguel Nebot and Olivier Gascuel for stimulating discussions and comments. RS receives a pre-doctoral fellowship from the Spanish MEC, while BW's stay at the University of Valencia was supported by a Marie-Curie fellowship. This work was financed by the project BMC2001-3096 of the MiCyT (Spain).
|
REFERENCES |
|---|
|
TopAbstractMETHODRESULTSDISCUSSIONREFERENCES |
|---|
-
Akaike H. A new look at the statistical model identification. IEEE Trans. Autom. Contr. (1974) 19:716–723.[CrossRef]
Ball J. K., Holmes E. C., Whitwell H., Desselberger U. Genomic variation of human immunodeficiency virus type 1 (HIV-1): Molecular analyses of HIV-1 in sequential blood samples and various organs obtained at autopsy. J. Gen. Virol. (1994) 75:67–79.[Medline]
Berry V., Gascuel O. On the interpretation of bootstrap trees: Appropriate threshold of clade selection and induced gain. Mol. Biol. Evol. (1996) 13:999–1011.
Bulmer M. Use of the method of generalized least-squares in reconstructing phylogenies from sequence data. Mol. Biol. Evol. (1991) 8:868–883.[Web of Science]
Cavalli-Sforza L. L., Edwards A. W. Phylogenetic analysis. Models and estimation procedures. Evolution (1967) 21:550–570.[CrossRef][Web of Science]
Collins K. R., Quiñones-Mateu M. E., Wu M., Luzze H., Johnson J. L., Hirsch C., Tossi Z., Arts E. J. Human immunodeficiency virus type 1 (HIV-1) quasispecies at the sites of Mycobacterium tuberculosis infection contribute to systemic HIV-1 heterogeneity. J. Virol. (2002) 76:1697–1706.
Desper R., Gascuel O. Fast and accurate phylogeny reconstruction algorithms based on the minimum-evolution principle. J. Comput. Biol. (2002) 19:687–705.
Desper R., Gascuel O. Theoretical foundation of the balanced minimum evolution method of phylogenetic inference and its relationship to weighted least-squares tree fitting. Mol. Biol. Evol. (2004) 21:587–598.
Dopazo J. Estimating errors and confidence intervals for branch lengths in phylogenetic trees by a bootstrap approach. J. Mol. Evol. (1994) 38:300–304.[Web of Science][Medline]
Efron B. The jackknife, the bootstrap, and other resampling plans. Conf. Board. Math. Sci. Soc. Ind. Appl. Math. (1982) 38:1–92.
Efron B., Halloran E., Holmes S. Bootstrap confidence levels for phylogenetic trees. Proc. Natl. Acad. Sci. USA. (1996) 93:7085–7090.
Felsenstein J. Distance methods for inferring phylogenies: A justification. Evolution. (1983) 38:783–791.
Felsenstein J. Distance methods for inferring phylogenies: A justification. Evolution (1984) 38:16–24.[CrossRef][Web of Science]
Felsenstein J. Confidence limits on phylogenies: An approach using the bootstrap. Evolution (1985) 39:783–791.[CrossRef][Web of Science]
Felsenstein J. PHYLIP (Phylogeny Inference Package) version 3.6c. Distributed by the author. (1993) Seattle: Department of Genetics, University of Washington.
Felsenstein J. An alternating least-squares approach to inferring phylogenies from pairwise distances. Syst. Biol. (1997) 46:101–111.
Felsenstein J., Kishino H. Is there something wrong with the bootstrap phylogenies? A reply to Hillis and Bull. Syst. Biol. (1993) 42:193–200.
Fitch W. M., Margoliash E. Construction of phylogenetic trees. Science (1967) 55:279–284.
Gatanaga H., Oka S., Ida S., Wakabayashi T., Shioda T., Iwamoto A. Active HIV-1 redistribution and replication in the brain with HIV encephalitis. Arch. Virol. (1999) 144:29–43.[CrossRef][Web of Science][Medline]
Goldman N., Anderson J. P., Rodrigo A. G. Likelihood-based tests of topologies in phylogenetics. Syst. Biol. (2000) 49:652–670.
Gorry P. R., Bristol G., Zack J. A., Ritola K., Swanstrom R., Birch C. J., Bell J. E., Bannert N., Crawford K., Wang H., Schols D., de Clercq E., Kunstman K., Wolinsky S. M., Gabuzda D. Macrophage tropism of human immunodeficiency virus type 1 isolates from brain and lymphoid tissues predicts neurotropism independent of coreceptor specificity. J. Virol. (2001) 75:10073–10089.
Hillis D. M., Bull J. J. An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis. Syst. Biol. (1993) 42:182–192.
Huelsenbeck J. P. The robustness of two phylogenetic methods: Four-taxon simulation reveal a slight superiority of maximum likelihood over neighbor-joining. Mol. Biol. Evol. (1995) 12:843–849.[Abstract]
Keys B., Karis J., Fadeel B., Valentin A., Norkrans G., Hagberg L., Chiodi F. V3 sequences of paired HIV-1 isolates from blood and cerebrospinal fluid cluster according to host and show variation related to the clinical stage of disease. Virology (1993) 196:475–483.[CrossRef][Web of Science][Medline]
Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. (1980) 16:111–120.[CrossRef][Web of Science][Medline]
Kishino H., Hasegawa M. Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea. J. Mol. Evol. (1989) 29:170–179.[CrossRef][Web of Science][Medline]
Korber B. T. M., Kunstman K. J., Patterson B. K., Furtado M., McEvilly M. M., Levy R., Wolinsky S. M. Genetic differences between blood- and brain-derived viral sequences from human immunodeficiency virus type 1-infected patients: Evidence of conserved elements in the V3 region of the envelope protein of brain-derived sequences. J. Virol. (1994) 68:7467–7481.
Kuhner M. K., Felsenstein J. A simulation comparison of phylogeny algorithms under unequal evolutionary rates. Mol. Biol. Evol. (1994) 11:459–468.[Abstract]
Kumar S., Tamura K., Nei M. MEGA3: Integrated software for molecular evolutionary genetics analysis and sequence alignment. Briefings in Bioinformatics (2004) 5:150–163.
Larget B., Simon D. L. Markov chain Monte Carlo algorithms for the Bayesian analysis of phylogenetic trees. Mol. Biol. Evol. (1999) 16:750–759.[Web of Science]
Li W.-H. A statistical test of phylogenies estimated from sequence data. Mol. Biol. Evol. (1989) 6:424–435.[Abstract]
Li W.-H., Zharkikh A. Statistical tests of DNA phylogenies. Syst. Biol. (1995) 44:49–63.
McGrath K. M., Hoffman N. G., Resch W., Nelson J. A. E., Swanstrom R. Using HIV-1 variability to explore virus biology. Virus Res. (2001) 76:137–160.[CrossRef][Web of Science][Medline]
Morris A., Marsden M., Halcrow K., Hughes E. S., Brettle R. P., Bell J. E., Simmonds P. Mosaic structure of the human immunodeficiency virus type 1 genome infecting lymphoid cells and the brain: Evidence for frequent in vivo recombination events in the evolution of regional populations. J. Virol. (1999) 73:8720–8731.
Nei M., Stephens J. C., Saitou N. Methods for computing the standard errors of branching points in an evolutionary tree and their application to molecular data from humans and apes. Mol. Biol. Evol. (1985) 2:66–85.[Abstract]
Posada D., Crandall K. A. Modeltest: Testing the model of DNA substitution. Bioinformatics (1998) 14:917–918.
Rzhetsky A., Nei M. A simple method for estimating and testing minimum-evolution trees. Mol. Biol. Evol. (1992a) 9:945–967.[Web of Science]
Rzhetsky A., Nei M. Statistical properties of the ordinary least-squares, generalized least-squares, and minimum-evolution methods of phylogenetic inference. J. Mol. Evol. (1992b) 35:367–375.[CrossRef][Web of Science][Medline]
Rzhetsky A, Nei M. Theoretical foundation of the minimum-evolution method of phylogenetic inference. Mol Biol Evol. (1993) 10:1073–1095.[Abstract]
Saitou N., Imanishi T. Relative efficiencies of the Fitch-Margoliash, maximum parsimony, maximum likelihood, minimum evolution, and the neighbor-joining methods of phylogenetic tree construction in obtaining the correct tree. Mol. Biol. Evol. (1989) 6:514–525.[Web of Science]
Saitou N., Nei M. The neighbor-joining method: A new method for reconstructing phylogenetic trees. Mol. Biol. Evol. (1987) 4:406–425.[Abstract]
Sanjuán R, Codoñer F. M., Moya A., Elena S. F. Natural selection and the organ-specific compartmentalization of HIV-1 V3 hypervariable region. Evolution (2004) 58:1185–1194.[CrossRef][Web of Science][Medline]
Shimodaira H. An approximately unbiased test of phylogenetic tree selection. Syst. Biol. (2002) 51:492–508.
Shimodaira H., Hasegawa M. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Mol. Biol. Evol. (1999) 16:1114–1116.[Web of Science]
Sitnikova T., Rzhetsky A., Nei M. Interior-branch and bootstrap tests of phylogenetic trees. Mol. Biol. Evol. (1995) 12:319–333.[Abstract]
Sneath P. H. A., Sokal R. R. Numerical taxonomy (1973) San Francisco: Freeman.
Sokal R. R., Rohlf F. J. Biometry: The principles and practice of statistics in biological research (2000) New York: Freeman. 686–697.
Stanhope M. J., Waddell V. G. W., Madsen O., de Jong W., Hedges S. B., Cleven G. C., Kao D., Springer M. S. Molecular evidence for multiple origins of Insectivora and for a new order of endemic African insectivore mammals. Proc. Natl. Acad. Sci. USA (1998) 95:9967–9972.
Strimmer K., Rambaut A. Inferring confidence sets of possibly misspecified gene trees. Proc. Royal. Soc. London Ser. B (2002) 269:137–142.
Strimmer K., van Haeseler A. Quartet puzzling: A quartet maximum-likelihood method for reconstructing tree topologies. Mol. Biol. Evol. (1996) 13:964–969.[Web of Science]
Susko E. Confidence regions and hypothesis testing for topologies using generalized least-squares. Mol. Biol. Evol. (2003) 20:862–868.
Swofford D. L. PAUP*.Phylogenetic analysis using parsimony (*and other methods). [4.0beta]. (2002) Sunderland, MA: Sinauer Associates.
Swofford D. L., Olsen G. J., Waddel P. J., Hillis D. M. Phylogenetic inference. In: Molecular systematics.—Hillis D. M., Mortiz C., Mable B. K., eds. (1996) 2nd edition. Sunderland, Massachusetts: Sinauer Associates. Pp. 407–514.
Tajima F. Statistical method for estimating the standard errors of branch lengths in a phylogenetic tree reconstructed without assuming equal rates of nucleotide substitution among different lineages. Mol. Biol. Evol. (1992) 9:168–181.[Abstract]
Tajima F, Nei M. Estimation of evolutionary distance between nucleotide sequences. Mol. Biol. Evol. (1984) 1:269–285.[Abstract]
Takahata N., Kimura M. A model of evolutionary base substitutions and its application with special reference to rapid change of pseudogenes. Genetics (1981) 98:641–657.
Takahata N., Tajima F. Sampling errors in phylogeny. Mol. Biol. Evol. (1991) 8:494–502.[Web of Science]
Tateno Y., Nei M., Tajima F. Accuracy of estimated phylogenetic trees from molecular data. I. Distantly related species. J. Mol. Evol. (1982) 18:387–404.[CrossRef][Web of Science][Medline]
van't Wout A. B., Ran I. J., Kuiken C. L., Kootstra N. A., Pals S. T., Schuitemaker H. Analysis of the temporal relationship between human immunodeficiency virus type 1 quasispecies in sequential blood samples and various organs obtained at autopsy. J. Virol. (1998) 72:488–496.
Wang T. H., Donaldson Y. K., Brettle R. P., Bell J. E., Simmonds P. Identification of shared populations of human immunodeficiency virus type 1 infecting microglia and tissue macrophages outside the central nervous system. J. Virol. (2001) 75:11686–11699.
Yang Z. PAML: A program package for phylogenetic analysis by maximum likelihood. Comput. Appl. Biosci. (1997) 13:555–556.
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||