Genetics, Vol. 170, 419-431, May 2005, Copyright © 2005
doi:10.1534/genetics.103.025692

This Article Abstract Full Text (PDF) All Versions of this Article: genetics.103.025692v1 170/1/419    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Suchard, M. A. Search for Related Content PubMed PubMed Citation Articles by Suchard, M. A.

Stochastic Models for Horizontal Gene Transfer

Taking a Random Walk Through Tree Space

Department of Biomathematics, David Geffen School of Medicine, University of California, Los Angeles, California 90095-1766

¹ Address for correspondence: Department of Biomathematics, David Geffen School of Medicine, UCLA, 650 Charles Young Dr., Box 951766, Los Angeles, CA 90095-1766.
E-mail: msuchard{at}ucla.edu

	ABSTRACT

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
Horizontal gene transfer (HGT) plays a critical role in evolution across all domains of life with important biological and medical implications. I propose a simple class of stochastic models to examine HGT using multiple orthologous gene alignments. The models function in a hierarchical phylogenetic framework. The top level of the hierarchy is based on a random walk process in "tree space" that allows for the development of a joint probabilistic distribution over multiple gene trees and an unknown, but estimable species tree. I consider two general forms of random walks. The first form is derived from the subtree prune and regraft (SPR) operator that mirrors the observed effects that HGT has on inferred trees. The second form is based on walks over complete graphs and offers numerically tractable solutions for an increasing number of taxa. The bottom level of the hierarchy utilizes standard phylogenetic models to reconstruct gene trees given multiple gene alignments conditional on the random walk process. I develop a well-mixing Markov chain Monte Carlo algorithm to fit the models in a Bayesian framework. I demonstrate the flexibility of these stochastic models to test competing ideas about HGT by examining the complexity hypothesis. Using 144 orthologous gene alignments from six prokaryotes previously collected and analyzed, Bayesian model selection finds support for (1) the SPR model over the alternative form, (2) the 16S rRNA reconstruction as the most likely species tree, and (3) increased HGT of operational genes compared to informational genes.

The process of horizontal (or lateral) gene transfer (HGT) plays a critical role across all domains of life and in particular among prokaryotes (JAIN et al. 1999; KOONIN et al. 2001). For example, many prokaryotes are agile at quickly adapting to new environments. Often, this ability stems from the acquisition of new genes through HGT rather than through random mutation (LAWRENCE 1999). At least three mechanisms promote HGT in prokaryotes (JAIN et al. 2002). These include: (1) transformation in which free DNA sequences are absorbed from the environment, (2) conjugation between two different prokaryotic species, and (3) transduction of genetic material through viruses. Finally, HGT also has medical importance (BROWN 2003). In the field of infectious diseases, HGT among bacterial pathogens of antibiotic resistance genes has greatly contributed to the emergence of multidrug-resistant bacteria in clinical settings (LEVERSTEIN-VAN HALL et al. 2002). In the field of oncology, HGT may also affect tumor progression; BERGSMEDH et al. (2001) show that eukaryotic cells can transfer active oncogenes.

Three general methods have been employed to examine HGT. The first focuses on single genomes and identifies genes suspected to have been imported through HGT by examining variation in nucleotide base composition and codon usage patterns (LAWRENCE and OCHMAN, 1997). The latter two methods are comparative studies across species. One uses similarity approaches based on gene content to identify HGT (RAGAN 2001) and to propose average genome or species-level trees (SNEL et al. 1999), while the alternative method endorses phylogenetic reconstruction using orthologous genes (JAIN et al. 1999). Base composition and codon bias studies may perform poorly when compared to phylogenetic methods (KOSKI et al. 2001). Further, phylogenetic methods offer at least one advantage over similarity-based approaches. The reconstructed phylogenies have direct biological interpretability as descriptions of the underlying evolutionary histories of the different genes (DOOLITTLE 1999). If a reconstructed gene tree differs from the assumed phylogeny of the species being studied, then HGT is offered as a possible explanation (SYVANEN 1994). One intrinsic difficulty is that the true species tree is often itself unknown. Therefore, it is necessary to either fix the species tree to equal the inferred gene tree for a specially chosen gene, e.g., the 16S rRNA tree (WOESE 2000), or simultaneously estimate the species tree and gene trees given a biologically plausible model relating them. As a first step, several research groups have attacked the inverse problem of reconstructing a species tree given gene trees subject to HGT. Most notable are the parsimony-based reconciled tree work by Page and colleagues (e.g., PAGE 2000) and the algorithmic work of MIRKIN et al. (2003).

I propose a simple class of stochastic models for HGT that enable the simultaneous estimation of the underlying species tree relating a group of organisms and the gene trees subject to HGT for a set of orthologous gene alignments. These HGT models function in a hierarchical manner (SUCHARD et al. 2003a) in which standard Bayesian phylogenetic approaches (e.g., SINSHEIMER et al. 1996; YANG and RANNALA 1997; MAU et al. 1999; LI et al. 2000; HUELSENBECK et al. 2001) are used to reconstruct each gene tree from its corresponding gene alignment. Simultaneous to the reconstructions, the HGT models impose a second probabilistic distribution over the gene trees (MADDISON 1997). This hierarchical distribution describes the gene trees likelihoods given an unknown species tree and an unknown number of HGT events leading from that species tree to each gene tree. The model is fit in a Bayesian framework that naturally handles uncertainty in discrete parameters such as all the trees and the number of HGT events and compares various models using Bayes factors (SUCHARD et al. 2001). Stochastic models fit in statistical frameworks offer several advantages over parsimony approaches. First, parsimony may underestimate the number of HGT events linking the species tree to the gene trees. This consequence is similarly seen in parsimonious reconstructions of the tree themselves, in which the number of nucleotide substitutions is underestimated. Second, it is easier in a statistical framework to include measures of uncertainty and these levels may be high in the inferred gene trees given the sparse data from which they are reconstructed.

One additional advantage of building stochastic models for HGT is the ability to compare competing models and to incorporate possible differences in the stochastic processes across genes, while assessing the significance of these differences in a formal statistical framework. As one example of possible differences across genes, JAIN et al. (1999) propose the complexity hypothesis. Under this hypothesis, genes are divided into one of two classes, informational or operational genes. Between classes, the rates of HGT differ. It is suspected that rates are higher for operational genes than for informational genes. This hypothesis and others can be tested by integrating over all possible species trees and gene trees weighed by their posterior probabilities. This Bayesian model-averaging approach reduces the possible bias inherent in selecting a specific species tree, minimizes underestimation of the uncertainty associated with the hypotheses (TAYLOR et al. 1996), and eliminates the need for ad hoc analyses. Formal comparison of different models for HGT will help gather further insight into the underlying biological processes.

	MODEL

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
Within-gene reconstruction model:
I begin with a hierarchical framework for phylogenetic reconstruction using molecular sequence data Y (SUCHARD et al. 2003a). Data Y = (Y₁, ... , Y_K) consist of K naturally disjoint partitions. Partition data Y_k for k = 1, ... , K represent the aligned DNA sequences of length L_k from one specific gene per partition, sequenced from the same N taxa across all partitions. A hierarchical phylogenetic model enables the pooling of information across gene partitions to improve estimate precision in individual partitions, while permitting estimation and testing of tendencies in across-partition quantities. For HGT, such across-partition quantities include: (1) an overall species tree, (2) appropriate stochastic models from which to construct a probability distribution over individual gene trees given the species tree, and (3) the stochastic model parameters that may vary between different classes of genes.

To utilize standard Bayesian models for phylogenetic reconstruction (e.g., SINSHEIMER et al. 1996; YANG and RANNALA 1997; MAU et al. 1999; LI et al. 2000) within a gene partition, data Y_k further divide into ordered homologous sites Y_kl for l = 1, ... , L_k. Site data Y_kl = (Y_kl1, ... , Y_klN)^t contain one nucleotide from each taxon, such that Y_kln (A, G, C, T) or their ambiguous wildcards for n = 1, ... , N. I assume that sites within a partition are independent and identically distributed, and the likelihood of observing Y_kl is given by a multinomial distribution over the 4^N possible outcomes with ambiguous nucleotides being integrated over their possible realizations. The multinomial outcome probabilities become functions of an unknown tree _k that describes the relatedness of the N taxa, branch lengths t_k = (t_k1, ... , t_kB), and a model to describe nucleotide mutation along these branches, all within partition k.

I elect for a reversible, continuous-time Markov chain (CTMC) model for nucleotide substitution (FELSENSTEIN 1981) popularized by TAMURA and NEI (1993) (TN93). The TN93 model is further parameterized by two transition:transversion rate ratios, _k between purines A and G and _k between pyrimidines C and T, and the stationary distribution of the underlying Markov chain _k = (_kA, _kG, _kC, _kT). The final scale parameter in the TN93 model is fixed such that branch lengths measure the expected number of nucleotide substitutions between the nodes in _k that the branch connects. Because I assume a reversible model for nucleotide substitution and make no clock-like restrictions on branch lengths, the root of each tree is unidentifiable (FELSENSTEIN 1981). As a consequence, the descriptions of all trees to follow are unrooted with N – 2 internal nodes and B = 2N – 3 branches.

Across-gene hierarchical model:
Following the hierarchical framework of SUCHARD et al. (2003a), I take branch lengths t_k as exponentially distributed with unknown expected divergence µ_k within partition k and model

and

(1)

where V = (A, G, M)^t and = (_A, _G, _C, _T) are unknown across-partition-level expectations, variance-covariance matrix has diagonal form, and ^–2, ^–2, ^–2_µ, and N are unknown across-partition-level measures of precision. Leaving V, , , and N as unknowns specified only by hyperprior distributions and estimating these parameters simultaneously with the within-partition-level continuous parameters, _k, _k, and µ_k for all k, enables the borrowing of strength of information from one partition by another, producing more precise within-partition-level estimates. I assume conjugate (when possible) and flat or noninformative hyperpriors on these across-partition-level parameters, as discussed in SUCHARD et al. (2003a). While the development of hierarchical priors over the continuous within-partition-level parameters has been straightforward, constructing a hierarchical prior over gene trees _k that incorporates the stochastic nature of HGT is more involved. This is illustrated in the next section.

Horizontal gene transfer models:
To build a stochastic model for HGT, I first present a formal description of the set of all possible N-taxon trees, commonly referred to as "tree space" (BILLERA et al. 2001), as a mathematical graph and then discuss several possible random walks (D. ALDOUS and J. FILL, unpublished results) on this graph that mirror the observed effects of HGT.

There exist M = (2N – 5)!/2^N–3(N – 3)! possible trees relating N extant taxa (FELSENSTEIN 1981). On the basis of these M trees, I construct a graph = (, ) with vertex set and edge set . Each tree represents a different vertex, or node, in the graph, such that the size of the vertex set || = M. An edge uv of a graph describes a direct connection between two of the graph's vertices u, v . The number of edges emanating from a single vertex v defines its degree d(v). Two vertices that are joined together by a single edge are called adjacent. Restricting attention to simple graphs in which pairs of vertices may be connected to each other only by a single edge and no vertex is connected to itself by a looping edge, a single vertex v from graph may be adjacent from as few as zero to as many as M – 1 other vertices. The set of all vertices adjacent to v are its neighborhood (v) and the size of this neighborhood |(v)| = d(v). The specification of a neighborhood for each vertex completes the description of , and many choices are available.

Subtree-prune-regraft-based model:
One approach to defining neighborhoods for each possible tree stems from subtree transfer operations (ALLEN and STEEL 2001). Subtree transfer operators act on trees producing local rearrangements. Applying a subtree transfer operator to one tree results in the creation of one of several possible new topologies that differs from by an extent dependent on the operator. The collection of all trees one operation away from = v becomes its neighborhood (v) under that operator. Nearest-neighbor interchange (ROBINSON 1971), tree bisection and reconnection (SWOFFORD et al. 1996), and subtree prune and regraft (SPR) (HEIN 1990, 1993) are three examples. In light of the goals of this article, SPR offers an advantage over the former two operators because of its potential biological interpretation. Applying the SPR operator to = v with its resultant drawn from _SPR(v) mirrors the differences observed between a species tree and an individual gene tree affected by one HGT (or recombination) event (HEIN 1990, 1993; JAIN et al. 1999; ALLEN and STEEL 2001).

Figure 1 illustrates one realization of the SPR operator applied to a six-taxon tree. The operator works in two steps. The first step selects and cuts any branch in the initial tree, _initial. Cutting the branch prunes away a subtree, _subtree. This subtree then regrafts itself using the same cut branch to a new internal node obtained by subdividing a preexisting branch in _initial – _subtree.

View larger version (11K): [in this window] [in a new window] [Download PPT slide]   FIGURE 1.— Subtree prune and regraft operator applied to a six-taxon tree. (1) Operator selects and cuts any branch in the initial tree, pruning away a subtree. (2) Operator regrafts subtree by selecting and subdividing a preexisting branch in the remaining tree. (3) Resultant tree for this realization.

One straightforward stochastic process on any simple graph is an unweighted random walk. A random walk on proceeds from vertex to vertex along existing edges of the graph, generating a discrete-time Markov chain (DTMC), where the states of the chain are the visited vertices. As unweighted, the chain uniform randomly chooses its next vertex to visit from all neighbors of its current vertex. For this DTMC, the one-event transition probability matrix A has entries

(2)

defining the probability of u changing into v as a result of one random event. It should be noted that A is just the adjacency matrix of rescaled to be a stochastic matrix [i.e., _v(A)_uv = 1].

On the basis of K random walks on the graph _SPR induced by the SPR operator, I construct a hierarchical prior over the joint distribution of all gene trees _k. To accomplish this task, I assume:

• An unknown species tree exists.
  
• The vertex representing is the initial state of K Markov chains.
  
• The Markov chains are conditionally independent given and A.
  
• The vertex representing _k is the final state of the kth chain.
  
• And each chain is of unknown length 0 E_k < .
  

Figure 2 depicts one set of the possible paths of K = 4 Markov chains starting at species tree and ending at gene trees _k on a small portion of a representative graph. The lengths of paths E_k shown range from one to three. I illustrate no paths of length zero, but these realizations should be most likely. A parsimony-like analysis considering beginning and end points of the chains in Figure 2 would, for example, underestimate E₄ as zero instead of three.

View larger version (26K): [in this window] [in a new window] [Download PPT slide]   FIGURE 2.— One possible Markov chain realization on a simplified graph for the species tree and four gene trees 1, ... , 4. All chains begin at the same vertex representing the species tree (in white). The chain producing gene tree 1 has length E1 = 3 (blue), the chain for 2 has length E2 = 1 (red), the chain for 3 has length E3 = 2 (yellow), and the chain for 4 has length E4 = 3 (green). Note that this latter chain returns to its starting state; a parsimony-like analysis would estimate E4 = 0. Not depicted are chains with actual length zero; these are most probable a priori.

(3)

To complete the hierarchical specification, I assign a prior distribution over by letting

(4)

where z = (z₁, ... , z_M) are constants, the prior probabilities of the M possible N-taxon trees. When little or no information is available about , one reasonable choice is z₁ = ... = z_M = 1/M; alternately, one may choose z such that the prior odds of competing hypotheses regarding are one in a hypothesis-testing setting (SUCHARD et al. 2003a). A further choice is discussed later. I further assume a conditionally independent prior on all E_k,

(5)

where _k is the expected number of HGT events for gene k and is a deterministic function of across-gene-level parameters. This prior is conjugate to (3), allowing all E_k to be integrated out of the model, improving sampling efficiency (LIU 1994),

(6)

Letting q(_k = v| = u, _k) = (P)_uv, the multistep transition probability matrix,

(7)

where I is the M x M identity matrix and Q = P – I is the CTMC infinitesimal rate matrix representation of the HGT process. In this parameterization, _k are scaled as the expected number of HGT events per gene. Let = (₁, ... , _K). Then, recalling the conditional independence assumption between Markov chains, the joint distribution over all gene trees _k becomes

(8)

Calculating the probabilities in (8) requires numerical methods to determine the matrix exponential involving P_SPR. These methods involve calculating the complete set of eigenvalues and eigenvectors of P_SPR, requiring (M³) operations. Such procedures become quickly computationally prohibitive as N, and hence M, increases. As a consequence, numerical approximations may be necessary to develop weighted graph extensions to _SPR directly. The weights in these extended graphs would be functions of unknown parameters and sampling these parameters would necessitate repetitive diagonalization.

Random walks with analytic solutions:
An alternative to this computational barrier involves using random walks on graphs for which analytic solutions are known for any size M. To help find such solutions, Equation 7 demonstrates the close connection between a DTMC with a Poisson-distributed number of events and a CTMC. In fact, any such DTMC can be expressed as a unique CTMC, called the "continuized" version (D. ALDOUS and J. FILL, unpublished results). Analytic solutions for several weighted and unweighted CTMC processes on a complete graph are commonly used in phylogenetics. In a complete graph, all vertices are adjacent to all others. The most notable examples are the CTMC models for nucleotide substitution. The simplest model by JUKES and CANTOR (1969) is unweighted. In the APPENDIX, I present the multistep transition probability matrix P_GJC for a generalized Jukes-Cantor (GJC) model involving an arbitrary number of vertices M. Proposed by KIMURA (1980), the next most sophisticated model for a complete graph is weighted. This model presupposes that the vertices are divided into two disjoint sets, _{1 2} = , and that transitions within and between ₁ and ₂ occur at varying rates. In terms of HGT, such a weighted random walk may prove useful to model varying rates of HGT between different groups of taxa. Letting M₁ = |₁|, M₂ = |₂|, and R equal the ratio of within- to between-transition rates, I present the multistep transition probability matrix P_GK given M₁, M₂, and R for a generalized Kimura (GK) model in the APPENDIX.

Modeling differences across gene classes:
I incorporate potential differences across genes in the expected number of HGT events _k by employing a generalized linear model (GLM) approach (MCCULLAGH and NELDER 1983). GLMs link the mean response, in this case _k, to a set of linear predictors. First, I divide all K genes into one of C possible classes, where the definition of the classes depends on the specific research question at hand. To identify gene-class membership in the GLM, I construct a K x C design matrix D = (D_kc), where matrix elements D_k1 = 1 for all k, representing the baseline multiplier for the reference class, and

(9)

for c = 2, ... , C, representing the offset multipliers for the remaining classes. Such a design matrix is standard in regression problems involving categorical dependent variables. I model

(10)

where linear combinations of predictors = (₁, ... , _C) specify, on the log-scale, the expected number of HGT events for all classes. I complete the hierarchical prior specification by assuming

(11)

I set L = (–2, 0, ... , 0) and = diag(10, ... , 10). This provides a quite diffuse prior on , with the median expected number of HGT events per gene 0.14 (GARCIA-VALLVE et al. 2000) for all classes.

As an example of how this GLM construction functions, consider the C = 2 classes case. Then,

(12)

When ₂ = 0, no difference across classes exists. Likewise, when ₂ < 0, the expected number of HGT events per gene is smaller in class 2 than in class 1, and when ₂ > 0, the expected number is larger.

	STATISTICAL FRAMEWORK

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
Comparing the relative appropriateness of the various stochastic models for HGT proposed in preceding sections and testing for significant differences in the expected number of HGT events across genes can be accomplished using Bayesian model selection via Bayes factors. Bayes factors are the Bayesian analog of the likelihood-ratio test (LRT), but suffer from fewer difficulties than LRTs in discrete spaces, when comparing non-nested models and with sparse data (SUCHARD et al. 2001). A Bayes factor B₁₀ measures the relative change in the support of the data Y in favor of one statistical model M₁ over another model M₀ and equals the ratio of the marginal likelihood m(Y|M₁) of M₁ over the marginal likelihood m(Y|M₀) of M₀ (KASS and RAFTERY 1995). To calculate Bayes factors, frequently more efficient methods than estimating the multidimensional integrals hidden in the marginal likelihoods directly are available.

When models are nested, a relatively simple Bayes factor calculation is available via the Savage-Dickey ratio (VERDINELLI and WASSERMAN 1995) and involves generating a posterior sample from the larger model only (SUCHARD et al. 2003b). For example, to assess the significance of differences across gene classes in the expected number of HGT events, let M₁ represent the unrestricted model proposed above. Nested within M₁ exists M₀, the equal-rates model, where _c = 0 for c = 2, ... , C. Further, the GJC model is nested within the GK model, as both are equal when R = 1.

On the other hand, the GJC and SPR models are non-nested, but both possess zero free parameters in their respective P matrices. For two arbitrary models M₀ and M₁ in situations like this, it is possible to estimate the posterior probabilities p(M₀|Y) and p(M₁|Y) by constructing a mixture model over the joint space of M₀ and M₁. By applying the Bayes theorem,

(13)

where q(M₀) and q(M₁) are the prior probabilities of models M₀ and M₁ in the mixture. Generally, I assume equal prior probabilities, q(M₀) = q(M₁) = ¹/₂, when reporting posterior estimates. However, improved efficiency in estimating B₁₀ can be garnered by adjusting these prior probabilities such that p(M₀|Y) p(M₁|Y) (CARLIN and CHIB 1995; SUCHARD et al. 2002).

Models SPR and GK neither are nested nor contain the same number of free parameters. One might entertain constructing a reversible-jump Markov chain Monte Carlo (MCMC) sampler (GREEN 1995) over the joint space of these models to compute the Bayes factor in support of SPR over GK. However, a simpler algebraic solution exists given the two preceding Bayes factor calculations,

(14)

To estimate all model parameters and Bayes factors, I employ MCMC. I further develop this MCMC algorithm and discuss its performance in the APPENDIX.

	EXAMPLE

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
To illustrate these stochastic models for HGT and methods to test hypotheses about them, I examine a large set of orthologous, prokaryotic genes collected by JAIN et al. (1999). The data consist of K = 144 separate gene alignments. Each alignment contains orthologous copies of a single gene from six prokaryotes. These prokaryotes are: Aquifex aeolicus (Aa), an early branching thermophilic eubacterium; Escherichia coli (Ec), a proteobacterium; Synechocystis 6803 (S6), a cyanobacterium; Bacillus subtilis (Bs), a gram-positive bacterium; Methanococcus jannaschii (Mj), a methanogen; and Archaeoglobus fulgidus (Af), a thermophilic sulfate-reducing methanogen relative. The first four organisms are Eubacteria, while the last two are Archaea. JAIN et al. (1999) construct the gene alignments on the basis of amino acid translations, assuming a star tree to reduce alignment bias (LAKE 1991), and classify each gene into one of two distinct classes, informational and operational genes (RIVERA et al. 1998). Table 1 lists the functional characteristics of the genes that fall into each class. As a generalization, informational gene products interact in large complex systems; this is especially true of the translational and transcriptional apparatuses. On the other hand, most operational gene products function independently or in small protein assemblies. In total, JAIN et al. (1999) assign 56 genes as informational and 88 as operational, employ these genes to examine the complexity hypothesis, and find support for higher levels of HGT among the operational genes as compared to the informational genes.

Selection of stochastic model:
I begin by comparing the relative likelihoods of the three different stochastic models, SPR, GJC, and GK. For the GK model, I define my two disjoint sets of trees as (1) those that support a split between the four Eubacteria and the two Archaea, ₁, and (2) those that do not, ₂. These definitions offer a first approximation to modeling differing rates of HGT within life domains and across domains in this example. HGT events that start and end in set ₁ are within domain transfers, while events that start in ₁ and end in ₂, or vice versa, are across domain transfers.

The log₁₀ Bayes factor in favor of SPR over GJC and the log₁₀ Bayes factor in favor of GK over GJC are

(15)

Figure 4a illustrates the scaled regeneration quantile (SRQ) plot for estimating the relative posterior probabilities used to calculate log₁₀ B_SPR,GJC. No substantial deviation from the slope = 1 line implies the MCMC chain is mixing sufficiently to generate this estimate. Combining the results in (15), I calculate the log₁₀ Bayes factor in favor of SPR over GK as

(16)

View larger version (24K): [in this window] [in a new window] [Download PPT slide]   FIGURE 4.— Scaled regeneration quantile (SRQ) plots to assess MCMC sampler performance when estimating four relative posterior probabilities. Plot a was generated when comparing the SPR and GJC stochastic models. Plots b–d were generated when comparing the three most probable species trees. No substantial deviation in the slopes from 1 (dashed lines) implies that the chains are mixing well enough to produce stable estimates.

Estimating the species tree:
Figure 3 displays the currently accepted species tree relating the six prokaryotes studied here. The four Eubacteria and two Archaea form two distinct clades (FENG et al. 1997) and Aa is the earliest branching species of the Eubacteria studied (DECKERT et al. 1998). The branching order of the remaining three Eubacteria Ec, S6, and Bs is more ambiguous (GIOVANNONI et al. 1996). The three possible resolutions of this trifurcation are depicted on the right side of Figure 3. Much of the debate surrounding the trifurcation depends on data choice and reconstruction methodology. For example, the top resolution produces species tree _Ec-S6 that places Ec and S6 as nearest neighbors. Protein synthesis elongation factor (EF) Tu gene reconstructions support this tree (LAKE and RIVERA 1996) and JAIN et al. (1999) fix _Ec-S6 as their reference tree in their analysis. Reconstructions of 16S rRNA phylogeny support the middle resolution of species tree _Bs-S6 (COLE et al. 2003) with Bs and S6 as nearest neighbors. The final resolution of species tree _Ec-Bs gains support from reconstructions of phenylalanyl-tRNA synthetase (TEICHMANN and MITCHISON 1999). However, even these three critical genes are subject to HGT (WOLF et al. 1999; ZAP et al. 1999; KE et al. 2000) and their reconstructed phylogenies may inaccurately represent the true species tree.

View larger version (19K): [in this window] [in a new window] [Download PPT slide]   FIGURE 3.— Species tree relating six prokaryotes. Species are: Aquifex aeolicus (Aa), Escherichia coli (Ec), Synechocystis 6803 (S6), Bacillus subtilis (Bs), Methanococcus jannaschii (Mj), and Archaeoglobus fulgidus (Af). Branch order of three Eubacteria Ec, S6, and Bs is under debate, leading to three possible subtrees (shown on right).

(17)

Similar to the back calculation completed in previous section, I estimate

(18)

while direct calculation of log₁₀ B_Bs-S6,Ec-S6 using the sampler yields approximately the same result. Figure 4, b–d, depicts the SRQ plots relevant to these Bayes factor calculations. Again, the MCMC chain appears well mixing. Although the posterior support for _Ec-Bs and _Ec-S6 initially appears quite small, on a relative scale it is not; probabilities for the remaining 102 trees are >15 orders of magnitude smaller.

Data sets as large as the K = 144 gene alignments from JAIN et al. (1999) are currently rare. Consequentially, I examine via simulation the number of alignments necessary to identify the species tree under the SPR model. Under this simulation, I randomly sample without replacement a fixed number of gene alignments K and then estimate the posterior support for _Bs-S6, assuming this is the true species-tree. I repeat this simulation 20 times for each value of K. For K = 2, the expected posterior probability of _Bs-S6 = 0.14. This estimate is approaching its prior value, signifying appropriate MCMC sampling with limited data. Approximately K = 50 gene alignments are required to achieve an expected posterior probability 0.80 and K = 70 are required for 0.90.

Hierarchical estimates of evolutionary pressures:
Table 2 presents the posterior estimates of the across-gene-level parameters used to pool information about (_k, _k, µ_k, _k). The table also lists posterior estimates of

(19)

These transformed variables report the across-gene-level averages of the two transition:transversion ratios and expected divergence on their usual, instead of log, scale. As seen from Table 2, the average transition:transversion ratio for purines A' is significantly different from the ratio for pyrimidines G', as the ratios' 95% Bayesian credible intervals (BCIs) do not overlap, and both ratios are greater than one. This supports the use of the TN93 model for nucleotide substitution over a more restricted model. Estimates of A, G, and are consistent with a previous study using a subset of the data in a hierarchical framework (SUCHARD et al. 2003a). Also in comparison to this previous study, differences in estimates of M, ²_A, ²_G, ²_M, and N all trend in the correct directions given the increase in the number of taxa and genes fit here.

(20)

offers substantial support (KASS and RAFTERY 1995) in favor of differing rates.

View larger version (27K): [in this window] [in a new window] [Download PPT slide]   FIGURE 5.— Analysis of the complexity hypothesis. The top two plots depict the posterior distributions of linear predictors 1 and 2 using histograms and normal approximations (solid lines). Also shown are the predictors' prior densities (dashed lines). Greater prior than posterior density at 2 = 0 (dotted line) supports a difference in HGT rates between gene classes. The bottom plot depicts the posterior distributions of the expected number of HGT events per gene for informational genes (light shading) and operational genes (dark shading).

	REMARKS

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
In this article, I proposed a simple class of stochastic models for HGT. The models are based on a random walk process in tree space and allow for the development of a joint distribution over multiple gene trees given an unknown species tree. I consider two general forms of random walks. The first stems from subtree transfer operations, in particular the SPR operator that mirrors the observed effects that HGT has on an inferred tree. The second form is based on walks over complete graphs and offers numerically tractable solutions for increasing number of taxa. I fit these models using a Bayesian framework to data from six prokaryotes. I find strongest support for the species tree that places Bs and S6 as nearest neighbors. This tree is supported by 16S rRNA reconstructions, but differs from the EF-Tu tree assumed by JAIN et al. (1999). I demonstrate the flexibility of these stochastic models to test competing ideas about HGT by examining the complexity hypothesis and find support for increased HGT of operational genes compared to informational genes. This latter finding remains unchanged if I fix the species tree to equal the EF-Tu tree (data not shown).

The specific stochastic models for HGT developed in this article have important limitations. First and foremost, the random walks explore only the discrete, topological portion of tree space and do not consider changes in branch lengths between trees as part of the underlying HGT process. As a result, HGT between nearest neighbors in a tree remains unidentified as this process does not result in a change in the topological configuration of the tree. Model extensions that consider a continuous random drift process on the joint space of (, t) (BILLERA et al. 2001) may circumvent this shortfall. For a related problem involving coalescence, YANG (2002) shows that including branch lengths t into the probabilistic model across loci improves power. Additionally, I assume that the K DTMCs representing the random walks of the gene trees _k away from the species tree are conditionally independent given . This assumption implies that the evolutionary histories of all genes are unlinked, while evidence for the HGT of, at a minimum, complete operons abounds in prokaryotes (KOONIN et al. 2001). Possible modeling aspects include allowing for linked or partially linked genes.

HGT is not the only process that may cause incongruence between gene trees. Although the effects of lineage sorting should be minor given the extensive divergence between the species studied here, the inclusion of paralogous genes copies within the orthologous alignments may mislead inference. Also important, stochastic error due to sparse phylogenetic data, evolutionary model misspecification, and parallel/convergent evolution can falsely produce incongruence between trees (CAO et al. 1998). These effects should upwardly bias the inferred number of HGT events. However, I suspect this bias is less than one HGT event per gene as only a modest percentage of genes should be affected and the error should produce just minor changes in the inferred tree. There is no a priori reason to suspect that this bias differs between the informational and operational gene classes; so the bias does not affect the relative difference between classes in HGT rates and inference regarding the complexity hypothesis.

For the SPR model, numerical approximations to the matrix exponentials involving the multistep transition probability matrix P_SPR may offer promise in handling research problems with larger numbers of taxa N (MOLER and VAN LOAN 2003). As N increases, the square dimensions of P_SPR grow superexponential, while the size of the neighborhood of each vertex grows only as (N²). As a consequence, P_SPR becomes increasingly sparse. In this situation, the number of unique eigenvalues increases substantially slower than the matrix's dimension. Krylov subspace techniques (SIDJE and STEWART 1999) may stretch computational limits upward to N = 8 or more.

In spite of these limitations, these stochastic models for HGT offer several advantages over previous approaches to studying HGT using multiple orthologous gene alignments. Under these stochastic models, the species tree is an unknown parameter that may be either integrated out of the analysis as a nuisance parameter or estimated jointly with the multiple gene trees. Joint analysis decreases the possibility of bias introduced through fixing the species tree when knowledge about it is uncertain. A stochastic approach also overcomes the bias inherent in parsimony-like estimation. Further, the hierarchical framework in which the stochastic model sits enables the borrowing of strength in the estimation of all gene-partition-level estimates including the gene trees themselves. Finally, and most importantly, stochastic models lend themselves well to formal statistical testing, with no need for ad hoc procedures. The ability to compare differing models for HGT will continue to shed further insight into the underlying biological process.

	APPENDIX

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
Complete models:
To determine the multistep transition probability matrix P_GJC for the GJC model with M 2 states, I first recall that

(A1)

is generated from an unweighted complete graph. As a complete graph, it is trivially connected and, therefore, has a unique stationary distribution. This distribution is (1/M, ... , 1/M).

To determine the eigenvalues of Q_GJC, I write

(A2)

where Q_GJC is scaled such that _k is expressed in terms of the expected number of HGT events per gene, J is the M x M matrix of all ones, and I is the M x M identity matrix. Matrix J has a rank of one and, therefore, one nonzero eigenvalue that equals M/(M – 1). Given the eigenvalues of J and expression (A2), the M eigenvalues of Q_GJC become

(A3)

Like the standard Jukes-Cantor model, where M = 4, the GJC model for any M 2 continues to have only two distinct eigenvalues. Conceptually this results because the qualitative behavior of the underlying Markov chain does not change as the size of the state-space increases.

By letting _k , I see that the stationary distribution is the eigenvector corresponding to the 0 eigenvalue. By examining the other limiting case where _k = 0 and considering the initial conditions, algebraic rearrangement yields

(A4)

The state-space of the GK model is partitioned into two disjoint sets ₁ and ₂. Let M₁ = |₁| and M₂ = |₂|, where M₁ + M₂ = M, and let R be the ratio of rates for transitions within a structural set to transitions between sets. Then, following arguments similar to those above, one can find the multistep transition probability matrix P_GK for the GK model.

If u ₁, then

(A5)

where

(A6)

By symmetry, if u ₂, then

(A7)

where ₃ = M₂R + M₁. For R 1, note that there are four unique eigenvalues when M₁ M₂ and three unique eigenvalues otherwise. This is consistent with the standard Kimura model, in which M₁ = M₂ = 2 with three unique eigenvalues.

Sampling algorithm:
For each gene-partition k, let _k = (_k, t_k, _k, _k, µ_k, _k) and, then, assemble = (₁, ... , _K) to be the collection of all gene-level parameters. To specify the hierarchical prior parameters, let = (V, , , N, , ). Across-gene-level parameters also include R when considering the GK model and mixing parameter {0, 1} when comparing models SPR and GJC. I employ a MCMC approach to sample from each model's joint posterior distribution, p(, |Y). I generate samples from these posteriors using two nested Metropolis-within-Gibbs cycles, as laid out in SUCHARD et al. (2003a) for hierarchical phylogenetic models. The outer cycle first iterates over gene partitions k and then over the parameters in . Within each gene partition k, the inner cycle proceeds over the parameters in _k. With the exception of proposals for , , R, and , all parameter proposals follow those in SUCHARD et al. (2003a).

The multinomial prior placed on is conjugate to its sampling density. As a result, it is possible to Gibbs sample from its full conditional distribution for moderately small M. This full conditional distribution remains multinomial with M state probabilities given by

(A8)

where _k = v_k for all k and _– is the vector of all model parameters (, ) excluding . Similar to the reweighted prior approach to estimate , varying z can improve sampling efficiency when estimating the relative posterior probabilities of specific species trees .

I draw the transition ratio R and linear predictors via separate Metropolis-Hastings proposals. For R, I propose new parameter values by generating a normal random variate centered at the current value of R with a tunable variance s²_R. Given the high degree of correlation between column vectors in the design matrix D, I expect the posterior distribution of to also exhibit strong correlation. This expectation stems from a normal linear regression approximation to p(exp()|) that has a variance-covariance structure proportional to (D'D)^–1. As a consequence, component-by-component updating of _c in should lead to a slowly mixing MCMC chain (ROBERTS and SAHU 1997). To help ensure adequate mixing, I propose all _c simultaneously using a multivariate normal random variate centered at the current value of with a tunable variance-covariance matrix diag. I adjust the tunable variances such that proposals have acceptance rates of 30–40% (GELMAN et al. 1996) and fix the correlation matrix approximately equal to the posterior correlation of determined by a trial chain.

When comparing HGT models using a mixture approach, I sample the mixing parameter directly from its full conditional distribution in a Gibbs step,

(A9)

where

(A10)

for i = 0, 1, = u, and _k = v_k. Values a may be saved at each iteration and used to construct a Rao-Blackwellized estimator for p(M₁|Y) (SUCHARD et al. 2003a).

Finally, the inferred number of HGT events E_k for the SPR model can be recovered after posterior simulation. The full conditional distribution

(A11)

where = u and _k = v_k. Since _uvk 1,

(A12)

where

(A13)

As the full conditional distribution of E_k is bounded above, I can generate random draws from it using rejection sampling. Starting with a posterior sample , I draw one replicate E_k for each p = 1, ... , P. For each p, I first generate E* from a Poisson distribution and U from the uniform distribution. Then, if U (A^E*)_uv/(P)_uv, where ^(p) = u and , I set E_k^(p) = E*. Otherwise, I reject the current proposal and begin again by regenerating (E*, U).

MCMC performance:
I run my MCMC chains for 1.1 x 10⁵ outer Metropolis-within-Gibbs cycles, discard the first 10⁴ cycles as burn-in, and subsample every 10 cycles. This process retains P = 10⁴ posterior samples with decreased autocorrelation. The total chain length and burn-in time appear moderately longer than required by examining time-series plots of the model log-likelihood during simulation.

To assess the performance of the MCMC sampler, I employ scaled SRQ plots (MYKLAND et al. 1995; LI et al. 2000; SUCHARD et al. 2002). SRQ plots are useful to demonstrate adequate sampler mixing within discrete model parameters. For the primary measures in this study, two important discrete parameters are the species tree and the model mixture parameter . In particular, I use SRQ plots to assess mixing when comparing the relative probabilities of two possible species trees and of differing stochastic models for HGT. In these SRQ plots, the local slope around a given point depicts the ratio of the relative posterior probability estimate based on the entire MCMC chain to an estimate based on a short segment of the chain around that point. Substantial deviation of the slope from one implies that the sampler is slowly mixing and, as a result, the chain is not sufficiently long to generate stable estimates. For continuous model parameters and Bayes factors based on the Savage-Dickey ratio, I assess convergence by comparing posterior estimates obtained from simulations of at least five independent chains with starting values drawn directly from the model priors.

	ACKNOWLEDGEMENTS

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
I thank the Lake lab, in particular Jon Moore and Jim Lake, for stimulating my interest in HGT, for many provocative discussions, and for providing the SPR adjacency matrices and prokaryote sequences used in this study. I also thank John Huelsenbeck for his insights into HGT and Janet Sinsheimer and Vladimir Minin for commenting on this manuscript. Fred Fox and the National Science Foundation grant 9987641-sponsored University of California, Los Angeles, Training Program in Bioinformatics graciously made possible the computing facilities to fit all 144 alignments simultaneously. The complete data set is available to interested readers at http://www.biomath.medsch.ucla.edu/msuchard/datasets.html. I am supported in part by National Institutes of Health grants GM08042 and GM068955 and U.S. Public Health Service grant CA16042.

	LITERATURE CITED

TOP ABSTRACT MODEL STATISTICAL FRAMEWORK EXAMPLE REMARKS APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 

ALLEN, B., and M. STEEL, 2001 Subtree transfer operations and their induced matrices on evolutionary trees. Ann. Combinatorics 5: 1–15.

BERGSMEDH, A., A. SZELES, M. HENRIKSSON, A. BRATT, M. FOLKMAN et al., 2001 Horizontal transfer of oncogenes by uptake of apoptotic bodies. Proc. Natl. Acad. Sci. USA 98: 6407–6411.[Abstract/Free Full Text]

BILLERA, L., S. HOLMES and K. VOGTMANN, 2001 Geometry of the space of phylogenetic trees. Adv. Appl. Math. 27: 733–767.[CrossRef]

BROWN, J., 2003 Ancient horizontal gene transfer. Nat. Rev. Genet. 4: 121–132.[CrossRef][Medline]

CAO, Y., A. JANKE, P. WADDELL, M. WESTERMAN, O. TAKENAKA et al., 1998 Conflict among individual mitochondrial proteins in resolving the phylogeny of Eutherian ord