* Max-Planck-Institut für evolutionäre Anthropologie, Leipzig, Germany
Heinrich-Heine-Universität Düsseldorf, Germany
Forschungszentrum Jülich, Germany
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Key Words: sequence evolution site-specific substitution rates maximum likelihood tree reconstruction human mitochondrial DNA
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
To assess site-specific rate heterogeneity, several approaches were developed. They can be separated into three categories: Parsimonious attempts (Hasegawa et al. 1993; Wakeley 1993), maximum likelihoodbased methods (Felsentein 1981; Olsen, Pracht, and Overbeek, personal communication; Yang 1993; Nielsen 1997; Excoffier and Yang 1999; Meyer, Weiss, and von Haeseler 1999), and empirical "pairwise" methods (Van de Peer et al. 1993, 2000; Pesole and Saccone 2001).
In the parsimonious approach the number of substitutions each site requires in a most parsimonious tree is counted, and those sites that exceed a cut-off valuei.e., if a site falls in the upper 10% quantile of the distribution of substitutions (Wakeley 1993)are coined "fast sites." However, as the assumptions of parsimony are not always valid, this method tends to underestimate the amount of rate variation. Likewise, skewed base composition and biased substitutions among nucleotides (or amino acids), which both are known to mimic effects of rate variation among sites, cannot be taken fully into account. Thus, only approximate estimates of site-specific variability are inferred.
Maximum likelihood methods in contrast have a statistically well-defined basis and can cope with recurrent mutations, skewed base composition, and biased variation in rate by specifying a model of sequence evolution (Yang 1993). However, the drawback of such methods is that they are computationally intensive, such that the estimation of site-specific rates is not possible without assuming either a tree topology including branch lengths (Olsen, Pracht, and Overbeek, personal communication; Kelly and Rice 1996; Nielsen 1997) or a specific distribution of rates a priori (Felsenstein 1981; Yang 1995; Kelly and Rice 1996; Excoffier and Yang 1999; Meyer, Weiss, and von Haeseler 1999). Typically one assumes that rates are distributed according to the Gamma distribution. Yang (1994) suggested the discrete Gamma distribution, which has been used to obtain site-specific rate estimates (Excoffier and Yang 1999; Meyer, Weiss, and von Haeseler 1999). From a biological point of view, there is no reason to assume that sites evolve with "discrete" rates. Moreover, the discrete approach has the tendency to underestimate the rate of fast-evolving sites. Thus, the use of a continuous distribution would, in principle, be favorable (Meyer, Weiss, and von Haeseler 1999). Unfortunately, this is at present computationally infeasible.
Empirical "pairwise" methods circumvent the construction of a tree by focusing on pairs of sequences (Van de Peer et al. 1993, 2000; Pesole and Saccone 2001). In each sequence pair, a position contains either a pair of different nucleotides or the identical nucleotides, where the probability to observe the events depends on the site-specific substitution rate and the evolutionary distance separating the sequences. Van de Peer et al. (1993) used this relation to infer site-specific rates of alignment positions. Although their approach yields sensible results, their statistical basis is not well understood. Thus, we introduce a maximum likelihood framework for estimating site-specific substitution rates from pairs of sequences based on the ideas of Van de Peer et al. (1993). Moreover, we suggest an iterative maximum likelihood procedure to compute site-specific rates and the phylogenetic tree simultaneously. The performance of the method is evaluated by simulations. As an illustrative example, we applied the estimation procedure to the dataset of 53 complete mitochondrial DNA (mtDNA) sequences of humans (Ingman et al. 2000) and obtained the full rate spectrum.
![]() |
Modeling Substitutions |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
|
|
We will use d = r · t for the expected number of substitutions between two sequences.
Rate Heterogeneity
If we assume that the entries Qij of Q are the same for all sites in a sequence of length , then all sites are evolving according to the same model of sequence evolution and with the same rate of evolution. The latter assumption can be relaxed by introducing a site-specific scaling factor fl, l = 1, ...,
that defines the rate per site as
|
|
Inferring the Evolutionary Distance
The rate matrix Q specifies the transition probabilities from one nucleotide to another if d substitutions per site are expected by the relation
|
|
|
Inferring Site-Specific Rates
If the model of sequence evolution Q is known and is the same for each position in the sequence, then we estimate the site-specific scaling factor fl (l = 1, ..., ) as follows. For a given position in a sequence, we analyze (independent) pairs of sequences at that position. Instead of studying two sequences that accumulated d substitutions, we now consider k independent pairs of sequences at position l in a multiple sequence alignment. Let (
,
), ..., (
,
) denote the k pairs, that are d1, d2, ..., dk substitutions apart. From this data set, we can estimate the rate-specific factor fl by maximizing
|
Analyzing Real Data
To estimate f = ( f1, ..., f), equation 9 requires the knowledge of the d1, d2, ..., dk. Several approaches are possible: One can compute the pairwise distances for each pair separately and then estimate the site-specific rates. Another option is to use the known phylogeny of the 2k sequences to identify independent sequence pairs in the tree and to deduce their distances from the branch lengths in the tree. Here, independence means that the path i, (i = 1, ..., k) that connects the sequence pair (Xi, Yi) in the tree is disjoint from the remaining k - 1 paths in the tree. This approach will be abbreviated IP-method (independent pairs). A second approach is simply taking all possible distance pairs from a known phylogeny, ignoring the fact that the pairs are not independent. We will call this the AP-method (all pairs). Both methods require the knowledge of the phylogeny and the evolutionary model Q. Note also that the AP-method is a generalization of the method suggested by Van de Peer et al. (1993).
If, however, no information on the phylogenetic relationship of the sequences is available, we suggest an iterative procedure, which includes the estimation of the model of sequence evolution, the phylogeny, and heterogeneous rates among sites. To compute the maximum likelihood estimates of the tree, either PUZZLE (Strimmer and von Haeseler 1996) or DNAML from PHYLIP is used (Felsenstein 1981).
Algorithm 1: Iterative Procedure
Step 3 of the iterative procedure puts constraints on the rates and ensures that we are analyzing relative rates with respect to the average evolutionary rate.
![]() |
Efficiency |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
In real applications of the method, sequences are related by a tree, and estimates are therefore dependent. To test the performance of equation 9 to estimate f for a collection of aligned sequences, we employed computer simulations, in which we compared the estimated site-specific rates with the rates modeled ("true" rates).
To this end, random tree topologies based on the coalescence (Hudson 1991) were generated. Sequences were evolved according to the Jukes-Cantor model (Jukes and Cantor 1969) with Gamma-distributed rates using Seq-Gen (Rambaut and Grassly 1997), where the shape parameter is an indicator of the amount of rate heterogeneity. A small value of
implies a pronounced rate heterogeneity.
To test the influence of the number of sequence pairs on the reliability of the estimates, data sets with n = 25, 50, 100, 250 sequences were generated. Because of the complexity of the problem, we confine ourselves to three types of simulations. The simulation conditions are summarized in table 1.
|
Simulation 1IP versus AP
Here we wanted to elucidate the influence of the usage of all possible n(n - 1)/2 pairs of sequences and their respective distances on the estimation of f as compared to the usage of only n/2 pairs of sequences, if independent pairs were used.
Figure 1 shows the scatter plots for a sample of n = 25 and n = 100 sequences. For both approaches the accuracy of the estimates increases if the number of sequences is increased, as can be seen from the reduced variation of the dot plot. For n = 25 it sometimes happened that the maximization of equation 9 did not converge. In those instances an arbitrary large rate of 75 was assigned. For n = 100 the maximization always worked. Moreover, for large data sets the effect of observing horizontal lines of rate estimates (see fig. 1a) disappears. Thus a much finer resolution of rate estimates seems possible.
|
|
Figure 2 summarizes the results. The overall performance of the estimates increases as the number of sequences is increased, but the AP-method performs best if strong rate heterogeneity is present. If = 5.0 then an average correlation coefficient of 0.5212 is observed for n = 250 sequences, which increases to 0.9531 if
= 0.1 (see fig. 2). In the latter case the correlation coefficients range between 0.85 and 0.95, indicating a high degree of correlation between estimated rates and true rates. Moreover, the range of the distribution of the correlation coefficients is narrower if rate heterogeneity is strong. Also, the variance of the empirical distribution of correlation coefficients is reduced as n is increased.
|
Because the results are averaged over 100 random trees, our conclusions are independent of the underlying tree. Thus, the method seems to work with good accuracy as long as the data set is sufficiently large and rate heterogeneity is present.
Simulation 3The Iterative Procedure
Finally, we investigated the performance of the iterative procedure (algorithm 1). As maximum likelihood tree reconstruction is the time-limiting factor (step 3 of the algorithm), we employed the subsampling strategy for the simulated data with n = 100 and 250 sequences (Meyer, Weiss, and von Haeseler 1999). Here, 10 random subsamples of 50 sequences were drawn. For each subsample, site-specific rates were estimated using the iterative procedure, and the estimates at each site were averaged.
The iterative procedure stopped in all instances (see table 1 Sim3) after seven iterations. The analysis that follows is for the resulting tree and the corresponding rate vector. The crosses in figure 2 display the results of the iterative method. If we can compute the overall maximum likelihood tree (n 50), then the correlation coefficients of the iterative method fall within the range of the coefficients obtained when the tree is given (Sim 2). However, if n = 100 or 250, then we observe a reduced correlation compared to Sim 2. The correlation coefficients overlap with the coefficients one would obtain if 50 sequences were analyzed. Thus the reduction in correlation between true and estimated rates may be attributed to the fact that we have employed the subsampling strategy. If we were able to compute maximum likelihood trees for large data sets in reasonable time, this reduction in correlation would disappear.
Simulation 2 and Simulation 3 both show the phenomenon of a high correlation coefficient when strong rate heterogeneity is present and a reduced correlation coefficient when the sequence positions evolve with weak heterogeneity. Figure 3 shows the scatter plots of one simulation of the rate estimates after convergence of algorithm 1 for 50 sequences, assuming = 5.0 (fig. 3a),
= 0.8 (fig. 3b), and
= 0.1 (fig. 3c). Whereas the slope equals 1.0147 (
= 5.0), 1.0868 (
= 0.8), and 1.0195 (
= 0.1), the correlation coefficients vary from 0.4844, 0.7746, to 0.9164. Thus, while the slope is for all examples close to 1, we observe a reduced correlation for the weak rate heterogeneity case. An explanation may be the lack of fast-evolving positions, which are present for
= 0.1, say (see fig. 3c); these fast-evolving positions (true rates larger than 30), can be reliably estimated and cause the high correlation. In summary, our simulations show that we are able to estimate site-specific rates in reasonable computation time. Moreover, they show that knowledge of the phylogeny is not necessary to infer the rates.
|
![]() |
Site-Specific Rates of Human Mitochondrial DNA |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
![]() |
Discussion |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
With the proposed method, continuous rate estimates are obtained. This is a major advantage when pronounced rate heterogeneity is present, and thus the suggested method seems superior to other approaches, where evolutionary rates are more or less arbitrarily pooled in discrete categories (Excoffier and Yang 1999; Meyer, Weiss, and von Haeseler 1999).
If rate heterogeneity is low ( = 5.0), the correlation of true and estimated rates is small. This finding is explained by the fact that the absolute differences between rates are small. In fact, if
= 5.0, approximately 90% of the sites evolve with relative rates between 0.25 and 1.75. In contrast, for
= 0.1, only 14% of the sites are in that range, and the absolute differences in rates between sites are larger. Thus, we find sites evolving with relative rates between 0 and 50 (see fig. 3c). Our results show that the method is able to distinguish between fast- and slow-evolving sites with high reliability. This is true even when we have to estimate the tree with the iterative procedure. Here our method faces a bottleneck, which can be overcome only by faster maximum likelihood tree reconstruction programs (e.g., Schmidt et al. 2002). To estimate site-specific rates reliably, we need large data sets; however, at present it is impossible to estimate maximum likelihood trees for more than 50 sequences. More work is necessary to develop further approximation procedures.
Further extensions of our approach are straightforward: we need to incorporate more complex models of sequence evolution and to introduce a statistical framework to estimate the significance of estimates obtained. Moreover, the method makes the assumption that the distribution of variable sites is the same in all sequences. We do not expect that this is a problem for the population sequence data we have studied. However, if very distantly related sequences are studied, it is possible that the distribution of sites free to vary does change across the phylogeny (e.g., Fitch 1971; Miyamoto and Fitch 1995; Lockhart et al. 2000; Galtier 2001; Lopez, Casane, and Philippe 2002; Huelsenbeck 2002). At present it is not clear how to include rate changes at a single position in our analyses.
![]() |
Appendix |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
Acknowledgements |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
Footnotes |
---|
William Martin, Associate Editor
![]() |
Literature Cited |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Anderson, S., A. T. Bankier, B. G. Barrell, et al. (14 co-authors). 1981. Sequence and organization of the human mitochondrial genome. Nature 290:457-465.[ISI][Medline]
Aris-Brosou, S., and L. Excoffier. 1996. The impact of population expansion and mutation rate heterogeneity on DNA sequence polymorphism. Mol. Biol. Evol. 13:494-504.[Abstract]
Clayton, D. A. 2000. Transcription and replication of mitochondrial DNA. Hum. Reprod. 15:11-17.
De Rijk, P., Y. Van de Peer, I. Van den Broeck, and R. De Wachter. 1995. Evolution according to large ribosomal subunit RNA. J. Mol. Evol. 41:366-375.[ISI][Medline]
Excoffier, L., and Z. Yang. 1999. Substitution rate variation among sites in mitochondrial hypervariable region i of humans and chimpanzees. Mol. Biol. Evol. 16:1357-1368.[Abstract]
Felsenstein, J. 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368-376.[ISI][Medline]
Fitch, W. M. 1971. The nonidentity of invariable positions in the cytochomes c of different species. Biochem. Genet. 5:231-241.[CrossRef][ISI][Medline]
Galtier, N. 2001. Maximum-likelihood phylogenetic analysis under a covarion-like model. Mol. Biol. Evol. 18:866-873.
Gu, X., Y.-X. Fu, and W.-H. Li. 1995. Maximum likelihood estimation of the heterogeneity of substitution rate among nucleotide sites. Mol. Biol. Evol. 12:546-557.[Abstract]
Hasegawa, M., A. D. Rienzo, T. D. Kocher, and A. C. Wilson. 1993. Toward a more accurate time scale for the human mitochondrial DNA tree. J. Mol. Evol. 37:347-354.[ISI][Medline]
Hudson, R. R. 1991. Gene genealogies and the coalescent process. Oxf. Surv. Evol. Biol. 7:1-44.
Huelsenbeck, J. P. 2002. Testing a covariotide model of DNA substitution. Mol. Biol. Evol. 19:698-707.
Ingman, M., H. Kaessmann, S. Pääbo, and U. Gyllensten. 2000. Mitochondrial genome variation and the origin of modern humans. Nature 408:708-712.[CrossRef][ISI][Medline]
Jukes, T. H., and C. R. Cantor. 1969. Evolution of protein molecules. Pp. 21132 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press, New York.
Kelly, C., and J. Rice. 1996. Modelling nucleotide evolution: a heterogeneous rate analysis. Math. Biosci. 133:85-109.[CrossRef][ISI][Medline]
Kogelnik, A. M., M. T. Lott, M. D. Brown, S. B. Navathe, and D. C. Wallace. 1998. Mitomap: a human mitochondrial genome database1998 update. Nucleic Acids Res. 26:112-115.
Lockhart, P. J., D. Huson, U. Maier, M. J. Fraunholz, Y. V. D. Peer, A. C. Barbrook, C. J. Howe, and M. A. Steel. 2000. How molecules evolve in eubacteria. Mol. Biol. Evol. 17:835-838.
Lopez, P., D. Casane, and H. Philippe. 2002. Heterotachy, an important process of protein evolution. Mol. Biol. Evol. 19:1-7.
Meyer, S., G. Weiss, and A. von Haeseler. 1999. Pattern of nucleotide substitution and rate heterogeneity in the hypervariable regions I and II of human mtDNA. Genetics 152:1103-1110.
Miyamoto, M. M., and W. M. Fitch. 1995. Testing the covarion hypothesis of molecular evolution. Mol. Biol. Evol. 12:503-513.[Abstract]
Nielsen, R. 1997. Site-by-site estimation of the rate of substitution and the correlation of rates in mitochondrial DNA. Syst. Biol. 46:346-353.[ISI][Medline]
Pesole, G., and C. Saccone. 2001. A novel method for estimating substitution rate variation among sites in a large dataset of homologous DNA sequences. Genetics 157:859-865.
Rambaut, A., and N. C. Grassly. 1997. Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees. Comput. Appl. Biosci. 13:235-238.[Abstract]
Schmidt, H. A., K. Strimmer, M. Vingron, and A. von Haeseler. 2002. TREE-PUZZLE: maximum likelihood phylogenetic analysis using quartets and parallel computing. Bioinformatics 18:502-504.
Strimmer, K., and A. von Haeseler. 1996. Quartet puzzling: a quartet maximum likelihood method for reconstructing tree topologies. Mol. Biol. Evol. 13:964-969.
Swofford, D. L., G. J. Olsen, P. J. Waddell, and D. M. Hillis. 1996. Phylogenetic inference. Pp. 407514 in D. M. Hillis, C. Moritz, and B. K. Mable, eds. Molecular systematics. Sinauer Associates, Sunderland, Mass.
Tamura, K., and M. Nei. 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol. Biol. Evol. 10:512-526.[Abstract]
Tavaré, S. 1986. Some probabilistic and statistical problems in the analysis of DNA sequences. Pp. 5786 in M. S. Waterman, ed. Some mathematical questions in biology: DNA sequence analysis. The American Mathematical Society, Providence, R.I.
Van de Peer, Y., S. L. Baldauf, W. F. Doolittle, and A. Meyer. 2000. An updated and comprehensive rRNA phylogeny of (crown) eukaryotes based on rate-calibrated evolutionary distances. J. Mol. Evol. 51:565-576.[ISI][Medline]
Van de Peer, Y., J. M. Neefs, P. D. Rijk, and R. D. Wachter. 1993. Reconstructing evolution from eukaryotic small-ribosomal-subunit RNA sequences: calibration of the molecular clock. J. Mol. Evol. 37:221-232.[ISI][Medline]
Wakeley, J. 1993. Substitution rate variation among sites in hypervariable region 1 of human mitochondrial DNA. J. Mol. Evol. 37:613-623.[ISI][Medline]
Yang, Z. 1993. Maximum likelihood estimation of phylogeny from DNA sequences when substitution rates differ over sites. Mol. Biol. Evol. 10:1396-1401.[Abstract]
Yang, Z. 1994. Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J. Mol. Evol. 39:306-314.[ISI][Medline]
Yang, Z. 1995. A space-time process model for the evolution of DNA sequences. Genetics 139:993-1005.
Yang, Z. 1996. Among-site rate variation and its impact on phylogenetic analyses. Trends Ecol. Evol. 11:367-372.[CrossRef][ISI]
Yang, Z., and T. Wang. 1995. Mixed model analysis of DNA sequence evolution. Biometrics 51:552-561.[ISI][Medline]