Department of Biomathematics, UCLA School of Medicine
Department of Biostatistics, UCLA School of Public Health
Department of Human Genetics, UCLA School of Medicine
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Likelihood ratio tests for evolutionary models (for a review, see Huelsenbeck and Rannala 1997
) can be remiss in that the topology space of evolutionary relatedness is discrete, data are sparse, parameter estimates may lie on the boundaries, and standard likelihood asymptotics may not apply (Navidi, Churchill, and von Haeseler 1991
, 1993; Goldman 1993
; Sinsheimer, Lake, and Little 1996
; Lange 1997
; Whelan and Goldman 1999
). Using Markov chain Monte Carlo (MCMC) methods (Gilks, Richardson, and Spiegelhalter 1996
) to approximate posterior distributions allows us to broach evolutionary model selection in a computationally feasible manner, with topology determination as an application of reversible jump MCMC (Green 1995
). Although MCMC methods have previously been used in the reconstruction of evolutionary relatedness (Kuhner, Yamato, and Felsenstein 1995, 1998
; Rannala and Yang 1996
; Mau and Newton 1997
; Yang and Rannala 1997
; Larget and Simon 1999
; Mau, Newton, and Larget 1999
; Li, Pearl, and Doss 2000
), our methods differ from these in being fully Bayes with a proper, yet nontruncated and uninformative, prior in modeling assumptions, in likelihood computation, in proposal kernels, and in the range of hypotheses tested. The Bayesian hypothesis-testing approach we propose in this paper provides a framework to simultaneously infer evolutionary relationships and test a large set of modeling hypotheses, of which we illustrate only a few.
In the Materials and Methods section, we describe the data upon which evolutionary relatedness is determined and models for reconstructing evolutionary trees, we introduce a reversible jump MCMC approach to estimate these relationships, and we show that Bayes factor comparison of evolutionary models is possible using vague but proper priors and can be used without conditioning on a particular topology. To illustrate, in the Results section, we compare several hierarchical evolutionary models, examine the appropriateness of a molecular clock, and test the existence of multiple rate classes.
![]() |
Materials and Methods |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
We assumed that nucleotide sites were independent and identically
distributed (iid) within a set of sites evolving under the same
evolutionary constraints (rate class r). Consequently, the
likelihood of observing a given pattern
X1jrX2jr
XNjr within r was multinomially distributed, where the probability was determined by an unknown bifurcating topology
describing the evolutionary relatedness of the organisms, a set of branch lengths tb
T for b = 1, 2, ... , 2N - 3, and a Markovian model for evolutionary change along this topology (Sinsheimer, Lake, and Little 1996
). The set T did not necessarily maintain consistent definition between different topologies. As a result, N-taxon topologies were nonnested models, each supported on separate parameter spaces T(
).
Models of Evolution
A popular class of evolutionary models are continuous-time Markov chain models, parameterized in terms of a 4 x 4 infinitesimal rate matrix of nucleotide change Q and branch lengths that correspond to the expected number of changes between nodes per site. The matrix Q satisfies the condition Q1 = 0, leaving 12 nonnegative, off-diagonal parameters in the most general form of Q. The transition matrix
![]() |
We explore three nested evolutionary models which reduce the parameterization of Q. The most general, TN93 (Tamura and Nei 1993
), allows for differing evolutionary rates between purine-to-purine transitions (
), pyrimidine-to-pyrimidine transitions (
), and transversions (purine-to-pyrimidine or pyrimidine-to-purine; ß) and allows the general stationary distribution of the nucleotides (
) to vary subject to the constraints
m
0,
m = 1 for m
(A, G, C, T/U) and detailed balance, Q
=
. The resulting infinitesimal rate matrix is
|
Previous evolutionary reconstructions from nucleic acid data using some MCMC methods fix the stationary distribution at either an empirical estimate from the observed data (Li, Pearl, and Doss 2000) or at values determined by preliminary MCMC sampling (Mau, Newton, and Larget 1999
). Like Larget and Simon (1999)
, we have not adopted either of these approaches. Empirical estimates give equal weight to all taxa and may therefore be biased when taxon selection oversamples certain subgroups, while fixing parameters can lead to underestimation of the variance of other parameters. Instead, our MCMC approach samples all model parameters.
For all three models, the position of the root, the most recent common ancestor (MRCA) of all N taxa, is not estimable without further parameter restrictions (Felsenstein 1981
), such as a molecular clock. If we can identify an outgroup taxon on the same branch as the root, a molecular clock among the remaining N - 1 taxa is a nested submodel in our framework and thus can be tested. A molecular clock allows for a computationally advantageous parameterization (Mau, Newton, and Larget 1999
) and reduces by half the number of branch lengths to be estimated.
We extend the HKY85 parameterization to a mixture model containing R infinitesimal rate matrices Qr and R sets of branch lengths Tr, where r = 1, ... , R, R is the number of different site classes present in the data, and each site in the data is assigned a priori to belong to class r. We choose HKY85 as an example for comparison with previous work (Yang 1995
; Larget and Simon 1999
) and note that such mixture models are easily implemented for TN93 or K80 as well. This mixture model is applicable when the reading frame of the DNA sequence is known and rate assignments are based on codon position in the reading frame or when data from different genes known to evolve under different selective pressures are combined. The model is a generalization of the Bayesian computation of Larget and Simon (1999
), in which they estimate multiple Qr matrices but not multiple branch lengths, and of the work of Yang (1995
), where he assumes that the branch lengths between different classes are scalar multiples. Multiple Qr matrices allow for different transition/transversion ratios and stationary distributions across classes, and multiple Tr sets allow for varying rates of evolution both across classes and between species. Yang's (1995)
scalar multiple branch lengths assume that the relative rates of evolution between classes are constant across species.
Bayesian Computation
Priors
Priors must remain proper to estimate Bayes factors. We employ flat or vague but completely proper priors over the entire parameter space (,
(
)), where
(
) = (
,
,
, T, µ) and µ is a hyperparameter to help define the prior for tb in T. When we employ multiple site classes,
= (
1, ... ,
R),
= (
1, ... ,
R), T = (T1, ... , TR), and µ= (µ1, ... , µR). We assume a prior for all parameters that is independent of topology such that q(
(
) |
) = q(
) and that all components of
are, a priori, independent. For the TN93 model, we set
|
Computation
We employ a Metropolis-within-Gibbs (Tierney 1994
) sampler using reversible model jumping (Green 1995
) among topologies. The parameter space dimension remains constant across topology models, although interpretation of portions of T is model-dependent. Each new state of the Markov chain is proposed via a Gibbs cycle. In each step within the cycle, a single parameter block is updated conditional on the remaining blocks using a Metropolis-Hastings algorithm (Metropolis et al. 1953
; Hastings 1970
). We use the following update cycle:
|
Previous MCMC approaches have updated and T simultaneously (Larget and Simon 1999
; Mau, Newton, and Larget 1999
); however, these approaches consider, at most, proportional rates of evolution between rate classes. Here, we include an extra T-only block to improve mixing within different sets of branch lengths for each class. We give our transition kernels for each Metropolis-Hastings step in appendix A. Where possible, we employ transition kernels that are symmetric to decrease computational complexity and are supported on the same bounded or discrete space as the kernel's underlying parameters to increase acceptance probabilities.
Each run of our MCMC chain consists of 500,000 full update cycles, and we disregard the first 100,000 steps as burn-in. For the starting state, we draw , µ, T | µ,
, and
directly from the prior distributions and set
equal to observed nucleotide frequencies in each class. We estimate functions of the chain's posterior by subsampling every 40 steps after burn-in. Multiple chains are run to insure adequate convergence. We use D =
tb, representing the total divergence between all taxa, µ,
,
, and
, to assess convergence within and across topologies. These parameters retain their interpretation as the sampler moves between topologies and may be used effectively to monitor how well the MCMC sampler is performing (Brooks and Guidici 1999
).
We calculate the likelihood of the data given , T,
,
, and
by integrating out the unknown states of the internal nodes using the pruning algorithm of Felsenstein (1981)
.
Model Comparisons
We make comparisons among models using Bayes factors (Kass and Raftery 1995
). The Bayes factor in favor of model M1 against model M0, given the data Y, can be expressed as
|
If model M0 is nested within another model M1 such that the parameter space of M1 is 1 = (
,
) and the parameter space of M0 is
0 = (
0,
) where
0 is a known constant with q0(
)
q1(
=
0,
), then B10 may be estimated via posterior simulation of M1 using the Savage-Dickey ratio (Verdinelli and Wasserman 1995
),
|
Restrictions on Evolutionary Rates
The HKY85 model is a restriction of TN93, and the K80 model is a restriction of HKY85. To test the appropriateness of the restricted models, we generate a posterior sample of the joint (,
,
) using our MCMC sampler under our most general TN93 model. We then estimate the Bayes factors in favor of TN93 against HKY85 and in favor of TN93 against K80. We then generate a posterior sample of the joint (
,
) under HKY85 to estimate the Bayes factor in favor of HKY85 against K80.
We approximate the posterior densities of (,
,
) and (
,
) using a normal approximation with the estimated posterior mean and posterior covariance evaluated at the joint restriction
=
and at (
m = 1/4,
=
) (in the former case) and at
m = 1/4 (in the latter case). We directly calculate the appropriate prior densities at these restrictions. When testing
=
, we recall that this restriction is equivalent to
-
= 0 and that the difference of two Uniform[0, 1) random variables is triangularly distributed on [-1, 1] with a density of 1 at the restriction (Feller 1971
). We then form the respective Bayes factors using equation (6)
.
Multiple Classes
In data sets putatively containing multiple site classes, we estimate the Bayes factors in favor of multiple Qr matrices under HKY85 by first generating a posterior sample of (1, ... ,
R,
1, ... ,
R) using our MCMC sampler. We approximate the posterior density at the restriction (
1 = ... =
R,
1 = ... =
R) using a normal approximation based on the posterior mean and posterior covariance of the sample reparameterized as
|
We form the Bayes factor by dividing this posterior density estimate by the prior evaluated at the joint restriction. The induced prior equals q(1 -
R, ... ,
R-1 -
R) x q(
1 -
R, ... ,
R-1 -
R) by prior independence. We identify that
|
![]() |
Following a similar derivation,
![]() |
Molecular-Clock Restrictions
To test the appropriateness of a molecular clock, we condition on the posterior mode topology, identify a known outgroup, and reparameterize the branch lengths in terms of coalescent height differences, ij. These parameters measure the difference in the sums of the branch lengths between two contemporary taxa i and j and their MRCA. Under a molecular clock,
ij = 0.
Given the outgroup R in the topology illustrated in figure 1
, a molecular clock constrains
|
|
|
Induced Priors on Coalescent Height Differences
For two- or three-taxon rooted topologies and for any arbitrary pair of taxa considered marginally in an N-taxon topology given an outgroup, we derive exact expressions for the ordinate of the induced prior q(ij = 0) in appendix B. To consider constraints jointly for N > 3, we estimate
, the ordinate of the induced joint prior q(
ij = 0), via simulation. We draw n = 50,000 samples of µ and T | µ directly from our priors and form the appropriate
ij from T. We calculate
using the multidimensional density estimator
|
Conditioning on d and recalling that
is the sum of independent Bernoulli random variables, the finite sample variance is approximated by
|
As a diagnostic for these simulations and the estimator, we compare with the analytic results developed in appendix B for n = 2 and 3. For n = 2, we calculate
= 1.0 ± 0.1 (estimate ± SD), and the exact result equals 0.955. For n = 3,
= 0.87 ± 0.02, while the exact result equals 0.897. To evaluate the estimator in higher dimensions, we drew 50,000 samples from a 12-dimensional multivariate N(1, I), where 1 = (1, ... , 1)t and I is the identity matrix, and obtained
= 3.5 x 10-8 ± 0.5 x 10-8, while the theoretical density is 4.0 x 10-8. These results return the theoretical densities to within the same order of magnitude and show only small simulation error or bias.
![]() |
Results |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Simulated Data
To insure that our methods would support a molecular clock if one were present, we simulated sequences of length 1,500 under a molecular clock for four contemporary taxa (A, B, C, D) and an outgroup (R) using the topology in figure 1
. We imposed a molecular clock by assigning branch lengths so that the evolutionary distances from MRCAs and contemporary taxa were equal. The approximate posterior density of ij = 0 was 1,106.3. A log10 B10 value of 0 implies that both models are equally likely, while values greater than 2 represent very strong evidence in support of the general model and values less than -2 represent very strong evidence in support of the restricted model (Kass and Raftery 1995
). The induced prior for
ij in this topology was 0.52, yielding a log10 B10 value of -3.32, which favors a molecular clock.
Tree of Life
The TOL data set consisted of 15 16S ribosomal RNA sequences (rRNA) (Lake 1988
). There were 1,039 aligned nucleotides after removal of gaps, and
obs = (0.2408, 0.3157, 0.2464, 0.1971)t. The species were drawn from four major classes of living organisms: eukaryotes, eubacteria, halobacteria, and eocytes, and also included the chloroplastic sequence from a eukaryote, Zea mays (chl.). Figure 2
(left) shows the modal topology (86% ± 3%, posterior probability mean ± SD determined from 10 independent chains) and the conditional posterior mean branch lengths estimated under the TN93 model. The model correctly clustered the eukaryotes, eocytes, halobacteria, and eubacteria into their appropriate monophyletic groups (clades) based on organism morphology and clustered the chloroplastic sequence in the eubacterial clade. This result is consistent with the endosymbiotic hypothesis of the origins of eukaryotic cellular organelles (Margulis 1981
) and has been demonstrated previously using rRNA (Lake 1988
; Bhattacharya and Medlin 1995
). Table 1 lists the marginal posterior means and standard errors of
,
,
, µ, and D under TN93, HKY85, and K80.
|
|
|
HIV
The HIV data contained the complete HIV genomes of two subtype D isolates, eight subtype B isolates, and one ADI subtype recombinant, MAL (Korber et al. 1997
). The subtype B isolates JRCSF and JRFL were collected from the same patient. There were 7,969 sites after removing all gaps in the aligned genomes, and
obs = (0.3698, 0.2365, 0.1708, 0.2229)t. Figure 4
displays the two topologies that account for virtually 100% of the posterior. These topologies were drawn with their conditional posterior mean branch lengths under TN93. One internal branch within the subtype B clade was approximately zero. At zero, the two shown topologies become equivalent. The sampler placed JRFL and JRCSF as nearest neighbors and correctly clustered the D and B subtypes. Table 1
gives the estimated evolutionary parameters and divergence.
|
|
|
|
Also in table 4 , we calculate the joint posterior using a multivariate normal approximation and prior densities via simulation evaluated at 0, and by taking the ratio of these values, we report the joint log10 B10 for each data set. The TOL and HIV examples offered strong support against a molecular clock (log10 B10 = 31.8 and 12.3, respectively), while the primates offer weaker support against a molecular clock (log10 B10 = 1.3).
We examined three subsets of our examples identified as interesting by the marginal diagnostics: (1) the eight subtype B isolates, (2) the anthropoids, and (3) the eukaryotes. Table 4 displays the corresponding joint log10 posterior and prior densities and log10 B10 for these subsets. As an illustration, figure 2 (right) plots the marginal posterior and prior distributions of the three coalescent height differences among the eukaryotes. The eukaryotes continued to offer strong support against a molecular clock (log10 B10 = 14.0), while the much more closely related B subtype isolates and anthropoids offered strong support in favor of a local molecular clock (log10 B10 = -3.7 and -2.4, respectively).
![]() |
Discussion |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
In allowing the sampler to explore the posterior across topology models, we overcome a shortfall of traditional analysis used to compare different continuous-time Markov evolutionary models. One can use a likelihood ratio test by maximizing the likelihood of the general model and the likelihood of the restricted model conditional on the same topology; however, the topologies that maximize the likelihood may differ under the two models. Then, general and restricted evolutionary models are no longer nested, and formal inference under a likelihood ratio test is no longer possible. In effect, our reversible jump MCMC sampler integrates out the nonnested portions of the parameter space. The Bayesian approach also allows us to effectively incorporate the uncertainty in the topology into the variance of the parameter estimates. Frequentist inference is forced to condition on topology and therefore underestimates the uncertainty.
The TOL example offers strong evidence against the universal appropriateness of a molecular clock; however, the anthropoids and subtype B isolates demonstrate that a local molecular clock for closely related taxa is a reasonable model. This finding is quite insensitive to prior choice. A molecular clock was originally employed in MCMC methods for evolutionary reconstruction to reduce computation (e.g., Mau and Newton 1997
; Yang and Rannala 1997
; Mau, Newton, and Larget 1999
; Li, Pearl, and Doss 2000
), but numerous examples of restriction violations exist (Ayala, Barrio, and Kwiatowski 1996
; Leitner et al. 1996
; Hillis, Mable, and Moritz 1996
; Simon et al. 1996
; Holmes, Pybus, and Harvey 1999
; Richman and Kohn 1999
). Larget and Simon (1999)
show that eliminating the molecular clock by doubling the number of estimable branch lengths does not produce an intractable problem; we extend the computation to allow for multiple sets of branch lengths that are not constrained by a molecular clock. In doing so, we further the frameworks of Thorne, Kishino, and Painter (1998)
and Huelsenbeck, Larget, and Swofford (2000)
in several important ways. Thorne, Kishino, and Painter (1998)
introduce a Bayesian approach that does not impose a molecular clock by first assuming that the true relationship of the taxa under study is known with complete certainty and employing empirical Bayes priors that use the data twice. However, they do not provide a statistical test of the appropriateness of a molecular clock. Huelsenbeck, Larget, and Swofford (2000)
continue to assume that the true relationship is known and formulate a likelihood ratio test that may become difficult to interpret when the data are sparse. We overcome both shortfalls by providing a framework to statistically test the appropriateness of a molecular clock while not having to condition on a known a priori topology and simultaneously make inference about the appropriate parameterization of the infinitesimal rate matrix. Additionally, pairwise diagnostic Bayes factors we propose allow the researcher to conveniently identify portions of an evolutionary history that violate a molecular clock and portions that support a local molecular clock. On the other hand, Huelsenbeck, Larget, and Swofford (2000)
allow for the estimation of divergence times, while our approach does not eliminate the confounding of time and evolutionary rate.
To provide both large and small jumps among so many branch lengths, we choose a 50/50 mixture of two transition kernelsthe first updating all branch lengths simultaneously using a reflective normal driver with small variance, and the second randomly selecting and updating one branch length using a driver with large variance. This mixture removed initially poor convergence in the HIV data set that had small branch lengths as compared with the other two examples. We find quick convergence and sufficient mixing for up to at least 15-taxon topologies without a molecular clock.
|
![]() |
Acknowledgements |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
Footnotes |
---|
1 Keywords: phylogenetics
Markov chain Monte Carlo
nested hypothesis testing
Bayes factors
2 Address for correspondence and reprints: Janet S. Sinsheimer, Department of Human Genetics, UCLA School of Medicine, Los Angeles, California 90095-7088. janet{at}sunlab.ph.ucla.edu
![]() |
literature cited |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Ayala, F. J., E. Barrio, J. Kwiatowski. 1996. Molecular clock or erratic evolution?. A tale of two genes. Proc. Natl. Acad. Sci. USA. 93:1172911734
Bhattacharya, D., L. Medlin. 1995. The phylogeny of plastids: a review based on comparisons of small-subunit ribosomal RNA coding regions. J. Phycol. 31:489498[ISI]
Brooks, S. P., P. Guidici. 1999. Convergence assessment for reversible jump MCMC simulationsPp. 733742 in J. Bernardo, J. Berger, A. P. Dawid, and A. F. M. Smith, eds. Bayesian Statistics 6. Oxford University Press, Cambridge, Mass
Brown, W. M., E. M. Prager, A. Wang, A. C. Wilson. 1982. Mitochondrial DNA sequences of primates, tempo and mode of evolution. J. Mol. Evol. 18:225239[ISI][Medline]
Crandall, K. A. ed. 1999. The evolution of HIVJohns Hopkins University Press, Baltimore, Md
Durbin, R., S. Eddy, A. Krogh, G. Mitchinson. 1998. Biological sequence analysis: probabilistic models of proteins and nucleic acidsCambridge University Press, Cambridge, England
Feller, W.. 1971. An introduction to probability theory and its applicationsVol. 2, 2nd edition. John Wiley and Sons, New York
Felsenstein, J.. 1978. The number of evolutionary trees. Syst. Zool. 27:2733[ISI]
.1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368376
Gelman, A., G. O. Roberts, W. R. Gilks. 1996. Efficient Metropolis jumping rulesPp. 599608 in J. M. Bernardo, J. O. Berger, A. P. Dawid, and A. F. M. Smith, eds. Bayesian Statistics 5. Oxford University Press, Oxford, England
Gilks, W. R., S. Richardson, D. J. Spiegelhalter. 1996. Markov chain Monte CarloChapman and Hall, New York
Goldman, N.. 1993. Statistical tests of models of DNA substitution. J. Mol. Evol. 36:182198[ISI][Medline]
Green, P. J.. 1995. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika. 82:711732[ISI]
Hasegawa, M., H. Kishino, T. Yano. 1985. Dating the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22:160174[ISI][Medline]
Hastings, W. K.. 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. 57:97109[ISI]
Hayasaka, K., K. T. Gojobori, S. Horai. 1988. Molecular phylogeny and evolution of primate mitochondrial DNA. Mol. Biol. Evol. 5:626644[Abstract]
Hillis, D. M., B. K. Mable, C. Moritz. 1996. Applications of molecular systematics: the state of the field and a look to the futurePp. 515543 in D. M. Hillis, C. Moritz, and B. K. Mable, eds. Molecular systematics. 2nd edition. Sinauer, Sunderland, Mass
Holmes, E. C., O. G. Pybus, P. H. Harvey. 1999. The molecular population dynamics of HIV-1Pp. 177207 in K. A. Crandall, ed. The evolution of HIV. Johns Hopkins University Press, Baltimore, Md
Huelsenbeck, J. P., B. Rannala. 1997. Phylogenetic methods come of age: testing hypotheses in an evolutionary context. Science. 276:227232
Huelsenbeck, J. P., B. Larget, D. Swofford. 2000. A compound Poisson process for relaxing the molecular clock. Genetics. 154:18791892
Jeffreys, H.. 1998. Theory of probabilityOxford classic texts in the physical sciences. 3rd edition. Oxford University Press, New York
Jukes, T., C. Cantor. 1969. Evolution of protein moleculesPp. 21132 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press, New York
Kass, R. E., A. E. Raftery. 1995. Bayes factors and model uncertainty. J. Am. Stat. Assoc. 90:773795[ISI]
Kimura, M.. 1980. A simple model for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16:111120[ISI][Medline]
Korber, B., B. Hahn, B. Foley, J. W. Mellors, T. Leitner, G. Myers, F. McCutchan, C. L. Kuikeneds1997. Human retroviruses and AIDS 1997: a compilation and analysis of nucleic acid and amino acid sequencesTheoretical Biology and Biophysics Group, Los Alamos National Laboratory, Los Alamos, NM. (http://hiv-web.lanl.gov)
Kuhner, M., J. Yamato, J. Felsenstein. 1995. Estimating effective population size and mutation rate from sequence data using Metropolis-Hastings sampling. Genetics. 140:14211430
.1998. Maximum likelihood estimation of population growth rates based on the coalescent. Genetics. 149:429434
Lake, J. A.. 1988. Origin of the eukaryotic nucleus determined by rate-invariant analysis of rRNA sequences. Nature. 331:184186[ISI][Medline]
Lange, K.. 1997. Mathematical and statistical methods for genetic analysisSpringer, New York
Larget, B., D. L. Simon. 1999. Markov chain Monte Carlo algorithms for the Bayesian analysis of phylogenetic trees. Mol. Biol. Evol. 16:750759
Leitner, T., D. Escanilla, C. Franzn, M. Uhln, J. Albert. 1996. Accurate reconstruction of a known HIV-1 transmission history by phylogenetic tree analysis. Proc. Natl. Acad. Sci. USA. 93:1086410869
Li, S., D. K. Pearl, H. Doss. 2000. Phylogenetic tree construction using Markov chain Monte Carlo. J. Am. Stat. Assoc. 95:493508[ISI]
Loftsgaarden, D. O., C. P. Quesenberry. 1965. A nonparametric estimate of a multivariate density function. Ann. Math. Stat. 36:10491051[ISI]
McCabe, K. M., G. Khan, Y. H. Zhang, E. O. Mason, E. R. McCabe. 1995. Amplification of bacterial DNA using highly conserved sequences: automated analysis and potential for molecular triage of sepsis. Pediatrics. 95:165169[Abstract]
Margulis, L.. 1981. Symbiosis in cell evolution: life and its environment on the early earthW. H. Freeman, San Francisco
Mau, B., M. A. Newton. 1997. Phylogenetic inference for binary data on dendograms using Markov chain Monte Carlo. J. Comput. Graph. Stat. 6:122131[ISI]
Mau, B., M. A. Newton, B. Larget. 1999. Bayesian phylogenetic inference via Markov chain Monte Carlo methods. Biometrics. 55:112[ISI][Medline]
Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller. 1953. Equations of state calculations by fast computing machines. J. Chem. Phys. 21:10871092[ISI]
Navidi, W. C., G. A. Churchill, A. von Haeseler. 1991. Methods for inferring phylogenies from nucleic acid sequence data by using maximum likelihood and linear invariants. Mol. Biol. Evol. 8:128143[Abstract]
.1993. Phylogenetic inference: linear invariants and maximum likelihood. Biometrics. 49:54355
Nerurkar, V. R., H. T. Nguyen, W. M. Dashwood, P. R. Hoffmann, C. Yin, D. M. Morens, A. H. Kaplan, R. Detels, R. Yanagihara. 1996. HIV type 1 subtype E in commercial sex workers and injection drug users in southern Vietnam. AIDS Res. Hum. Retroviruses. 12:841843[ISI][Medline]
Rannala, B., Z. Yang. 1996. Probability distribution of molecular evolutionary trees: a new method of phylogenetic inference. J. Mol. Evol. 43:304311[ISI][Medline]
Relman, D. A., T. M. Schmidt, A. Gajadhar, M. Sogin, J. Cross, K. Yoder, O. Sethabutr, P. Echeverria. 1996. Molecular phylogenetic analysis of Cyclospora, the human intestinal pathogen, suggests that it is closely related to Eimeria species. J. Infect. Dis. 173:440445[ISI][Medline]
Richman, A. D., J. R. Kohn. 1999. Self-incompatibility alleles from Physalis: implications for historical inference from balanced genetic polymorphisms. Proc. Natl. Acad. Sci. USA. 96:168172
Rudolph, K. M., A. J. Parkinson, C. M. Black, L. W. Mayer. 1993. Evaluation of polymerase chain reaction for diagnosis of pneumococcal pneumonia. J. Clin. Microbiol. 31:26612666[Abstract]
Rzhetsky, A., T. Sitnikova. 1996. When is it safe to use an oversimplified substitution model in tree-making?. Mol. Biol. Evol. 13:12551265[Abstract]
Simon, C., L. Nigro, J. Sullivan, K. Holsinger, A. Martin, A. Grapputo, A. Franke, C. McIntosh. 1996. Large differences in substitutional pattern and evolutionary rate of 12S ribosomal RNA genes. Mol. Biol. Evol. 13:923932
Sinsheimer, J. S., J. A. Lake, R. J. Little. 1996. Bayesian hypothesis testing of four-taxon topologies using molecular sequence data. Biometrics. 52:193210[ISI][Medline]
Swofford, D. L., G. J. Olsen, P. J. Waddell, D. M. Hillis. 1996. Phylogenetic inferencesPp. 407514 in D. M. Hillis, C. Moritz, and B. K. Mable, eds. Molecular systematics. 2nd edition. Sinauer, Sunderland, Mass
Tamura, K., 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:512526[Abstract]
Thorne, J. L., H. Kishino, I. S. Painter. 1998. Estimating the rate of evolution of the rate of molecular evolution. Mol. Biol. Evol. 15:16471657
Tierney, L.. 1994. Markov chains for exploring posterior distributions (with discussion). Ann. Stat. 22:17011762[ISI]
Verdinelli, I., L. Wasserman. 1995. Computing Bayes factors using a generalization of the Savage-Dickey density ratio. J. Am. Stat. Assoc. 90:614618[ISI]
Whelan, S., N. Goldman. 1999. Distributions of statistics used for the comparison of models of sequence evolution in phylogenetics. Mol. Biol. Evol. 16:12921299
Yang, Z.. 1995. A space-time process model for the evolution of DNA sequences. Genetics. 139:9931005
Yang, Z., B. Rannala. 1997. Bayesian phylogenetic inference using DNA sequences: a Markov chain Monte Carlo method. Mol. Biol. Evol. 14:717724[Abstract]