Bioinformatics Advance Access originally published online on November 5, 2004
Bioinformatics 2005 21(7):975-980; doi:10.1093/bioinformatics/bti109
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
An alternative model of amino acid replacement
Department of Plant and Microbial Biology 111 Koshland Hall #3102 University of California Berkeley, CA 94720-3102, USA
*To whom correspondence should be addressed.
 |
Abstract |
|---|
 TOP Abstract INTRODUCTIONMATHEMATICAL DETAILSREFERENCES |
|---|
Motivation: The observed correlations between pairs of homologous protein sequences are typically explained in terms of a Markovian dynamic of amino acid substitution. This model assumes that every location on the protein sequence has the same background distribution of amino acids, an assumption that is incompatible with the observed heterogeneity of protein amino acid profiles and with the success of profile multiple sequence alignment.
Results: We propose an alternative model of amino acid replacement during protein evolution based upon the assumption that the variation of the amino acid background distribution from one residue to the next is sufficient to explain the observed sequence correlations of homologs. The resulting dynamical model of independent replacements drawn from heterogeneous backgrounds is simple and consistent, and provides a unified homology match score for sequence–sequence, sequence–profile and profile–profile alignment.
Contact: gec{at}compbio.berkeley.edu
|
INTRODUCTION |
|---|
|
TOPAbstractINTRODUCTIONMATHEMATICAL DETAILSREFERENCES |
|---|
During evolution, a protein's amino acid sequence is altered by the insertion and deletion of residues and by the replacement of one residue by another. In principle, the alignment of protein sequences and the subsequent detection of protein homologs and the inference of protein phylogenies requires a dynamical model of this sequence evolution. The most common and widely used residue replacement dynamics is the standard Dayhoff model, which assumes that the substitution probability during some time interval depends only on the identities of the initial and replacement residues and that the dynamics is otherwise homogeneous along the protein chain, between protein families and across evolutionary epochs. In other words, under this model the dynamics of amino acid substitution resembles a continuous time, first order Markov chain (Dayhoff et al., 1972, 1978; Gonnet et al., 1992; Jones et al., 1992; Müller and Vingron, 2000).
However, it has long been known that this widely used Markovian substitution model is fundamentally unsatisfactory. One major problem is that the short and long time substitution dynamics are incompatible (Gonnet et al., 1992; Benner et al., 1994; Müller and Vingron, 2000). Benner et al. (1994) suggest that this is because at short evolutionary times the patterns of substitution are influenced by single base mutations between neighboring codons, whereas for more diverged sequences the genetic code is irrelevant and the patterns of replacement are dominated by the selection of chemically and structurally compatible residues.
A more serious problem with the Dayhoff Markovian model is that it assumes that every residue in every protein has the same background distribution of amino acids and that protein sequences rapidly evolve to this uninteresting equilibrium. In actuality, the amino acid background distribution varies markedly from one residue position to the next, as can be seen, for example, in Figure 1. These large site-to-site variations are stable across relatively long evolutionary time-scales, and they account for the success of protein hidden Markov models and other profile based multiple sequence alignment methods. (See, for example, Sjölander et al., 1996; Durbin et al., 1998.) Profile methods can detect substantially more remote homologies than pairwise alignment (Park et al., 1998; Green and Brenner, unpublished data). In short, the dynamics of amino acid substitution are not Markovian, stationary, nor homogeneous, and the prediction of rapidly decaying sequence correlations is at odds with the success of profile based remote homology detection.
|
A natural solution to the limitations of the Markov model is to assume that residue replacement is governed by different Markov processes for each position, each process potentially possessing its own background distribution and substitution probabilities. The appropriate Markov matrix for a particular protein position is chosen based upon predictions of the protein structure, or directly from the sequence data (Goldman et al., 1996, 1998; Thorne et al., 1996; Topham et al., 1997; Koshi and Goldstein, 1998; Dimmic et al., 2000; Lartillot and Philippe, 2004). However, this approach is both computationally and conceptually complex.
Here, we propose that the observed sequence correlation between diverged homologs is principally due to the heterogeneous, stable, background distribution of each protein site and, therefore, that a Markovian amino acid replacement dynamics is overly complicated and possibly unnecessary for the accurate construction of protein sequence alignments and phylogenies. As an alternative, we construct a dynamical model of amino acid replacement that explicitly assumes that each protein site has a different equilibrium distribution of the 20 canonical amino acids (which we refer to as that site's amino acid background, 
) and that the residue distribution of each site rapidly (relative to evolutionary time-scales) relaxes to this local, site-specific equilibrium. These distributions themselves conform to a probability distribution of backgrounds, P(), which we may discover by studying many families of homologous proteins (Fig. 2). We do not need to model the short time dynamics with any great accuracy, since the alignment of highly conserved homologs is relatively straightforward. Therefore, we will assume that replacement residues are randomly sampled from the background distribution of that site. Consequentially (and in direct contrast to the Markov model) the replacement residue is conditionally independent of the initial amino acid at all times. Note that multiple substitutions are statistically equivalent to a single substitution (since a single mutation is sufficient to relax a site to local equilibrium) and that a residue can be replaced by the same amino acid type. The resulting dynamic is a site-specific, continuous time, zeroth order Markov chain, similar in spirit to Felsenstein's (1981) model of nucleic acid substitution. The crucial difference is that the initial and replacement residues are non-trivially correlated because both have been sampled from the same background distribution, whereas non-homologous residues are drawn from different backgrounds. Bruno (1996) has used essentially the same dynamical model discussed here, albeit without incorporating the background prior distribution, to find maximum likelihood estimates of site-specific amino acid frequencies.
|
Under our residue replacement model, the principal origin of sequence correlation between diverged homologs is the background distribution of each protein site. This is also the central idea underlying profile based multiple sequence alignment algorithms. Therefore, we are not proposing a radically different method for homolog detection or sequence alignment; rather we are proposing a concrete and consistent dynamics for the underlying evolutionary process. The implications of this dynamics can be readily extended to cover not only profile–sequence based alignment, but also profile–profile and pairwise sequence–sequence alignment. Moreover, when we consider pairwise, sequence alignment below, we find that our model is essentially equivalent to the standard pairwise alignment methods, as they are used in practice. This alternative dynamical model of amino acid replacement is biologically reasonable, conceptually straightforward and can adequately explain many of the observed patterns of homolog sequence correlation without invoking a Markovian dynamic.
The correlations between pairs of homologous residues can be summarized by an amino acid substitution matrix, S, whose entries represent the log probability of observing the homologous pair of amino acids qij in a properly aligned pair of homologous proteins, against the probability pi pj of independently observing the residues in unrelated sequences (Altschul, 1991):
 | (1) |
= 1/3, 
digits), although the scaling is arbitrary. Assuming conditionally independent replacements, we can directly construct the large time limit substitution matrix from the background probability distribution, P(). Fortunately, this large distribution has previously been investigated and parameterized by fitting many columns from multiple alignments of homologous protein sequences to a mixture of Dirichlet distributions (Karplus, 1995; Durbin et al., 1998). Figure 2 displays a projection of the dist.20comp parameterization (Karplus, 1995) and Figure 3 displays the corresponding substitution matrix. The mathematical details of matrix construction are given below.
|
bits, rounded to the nearest integer.At shorter evolutionary times there is a significant chance that no mutation has occurred at all, resulting in an enhanced probability of amino acid conservation. Let c be the probability of zero mutation events; then the substitution matrix, adjusted for the possibility of zero mutations, is
 | (2) |
ij is the Kronecker delta function. A reasonable default model for the conservation probability c would be to assume that substitutions are Poissonian. Then c = exp(–t/
), where is the mean time between replacements. Note that although replacement is Markovian (albeit zeroth order), the dynamic decay of residue correlations at a position is not, due to the heterogeneity of the amino acid background at that position (an unobserved hidden variable). Equation (2) gives the log odds of aligned residues, given the inter-sequence divergence. Conversely, given a prior on the parameter c and a fixed alignment we can invert Equation (2) and estimate the divergence between sequences. Conserved residues indicate small divergences and unconserved pairs argue for large divergences, although different pairs are weighted differently.
As evolution proceeds, the background distribution of a site may itself change, due, for example, to a change in structure of that part of the protein. This will result in a loss of homology signal under our model, and it may no longer be possible to align the diverged residues, nor to recognize them as homologs. This is schematically illustrated in Figure 4.
|
Various families of substitution matrices have been developed, including PAM, BLOSUM and VTML. Different members of the same family represent different degrees of sequence divergence. In principle, we should match the divergence inherent in the substitution matrix to the divergence of the pair of sequences we wish to align (Altschul, 1993). However, this is computationally expensive, and, in practice, a single matrix is chosen based on its ability to align remote homologs, on the grounds that matching close homologs is relatively easy (Brenner et al., 1998). Under the Markov model, the chosen matrix has no particular significance. On the other hand, under our model there is a natural, non-trivial, long time limit matrix [Fig. 2; Equation (2), c = 0]. This matrix represents the sequence correlations at any time after the first few mutations, and before the underlying amino acid background itself diverges (Fig. 4).
Figures 5 and 6 demonstrate that the dist.20comp matrix represents a similar level of evolutionary divergence, and similar patterns of substitution as BLOSUM62 and VTML160, two substitution matrices commonly used for pairwise sequence alignment and remote homology detection. These three matrices have been created using very different evolutionary models; the dist.20comp matrix is based upon our heterogeneous, background/independent substitution model, and the dist.20comp background distribution is, in turn, derived from the columns of many multiple alignments of homologous protein sequences; the popular BLOSUM matrices are empirically derived from the BLOCKS database of reliable protein sequence alignments (Henikoff and Henikoff, 1992; Henikoff et al., 2000); and the classic PAM (Dayhoff et al., 1978) and modern VTML (Müller et al., 2002) matrix families are explicitly based upon the Markovian model of amino acid replacement. In a recent evaluation of pairwise remote homology detection, the VTML160 matrix was found to be more effective than any other VTML, PAM or BLOSUM matrix (Green and Brenner, 2002). However, as can be seen in Figure 6, the difference in remote homology detection ability of the three matrices is relatively small.
|
|
In summary, the important sequence correlations can be adequately explained by assuming conditionally independent replacements drawn from background distributions that vary from site to site, but are stable over evolutionary time-scales. The standard, Markovian model of amino acid replacement is unnecessary, overly complicated and inconsistent with observed substitution patterns.
This alternative, heterogeneous background, independent substitution model may be particularly useful for simultaneous sequence alignment and phylogenetic tree reconstruction, since it is necessary to align pairs of close homologs at the leaves, and multiply align many remote homologs at the interior nodes of the tree. Therefore, a simple (yet realistic) evolutionary dynamic that is consistent across a wide range of divergence times, and that leads naturally to sequence–sequence, sequence–profile and profile–profile alignment algorithms, may be advantageous.
|
MATHEMATICAL DETAILS |
|---|
|
TOPAbstractINTRODUCTIONMATHEMATICAL DETAILSREFERENCES |
|---|
A collection of homologous residues can be represented by a 20-component canonical amino acid count vector, n = {n1,...,n20}. The total number of counts can be 1, if the observation is taken from a single sequence, or many if the collection represents an entire column of a multiple sequence alignment or some other related set of residues.
In general, we wish to estimate whether two collections of homologous residues are related, given that detectably homologous residues are drawn from the same background amino acid distribution. The appropriate test statistic is the log odds of sampling the two amino acid count vectors from the same, but unobserved, background distribution, against the probability of independently sampling the two count vectors from different distributions:
 | (3) |
= {1,...,20}, follows the multinomial distribution;
 | (4) |
The probability distribution of background distributions P() has been studied and measured by collating columns from many multiple protein sequence alignments. Since this is a very large, multi-modal probability it is necessary to parameterize the distribution into a convenient representation. Typically, a mixture of Dirichlet distributions is used (Karplus, 1995; Sjölander et al., 1996; Durbin et al., 1998):
 | (5) |
k sum to one. The kth Dirichlet distribution is itself parameterized by the 20-component (canonical amino acids) non-negative vector 
k,
 | (6) |
i i. Dirichlet mixtures are used to model the background probability partially because Dirichlet distributions are naturally conjugate to the multinomial distribution (and therefore mathematically convenient) and partially because a Dirichlet mixture can approximate the true distribution with a reasonably small number of parameters. The underlying assumption is that the probability distribution is smooth, but lumpy (Fig. 2).
If we do not know the particular background from which the observations have been drawn, then we must average over all backgrounds to find the probability of observing a particular count vector:
 | (7) |
k). Therefore, the integral over must be equal to the corresponding Dirichlet normalization constant, Z(n + k). The final result is a mixture of multivariate negative hypergeometric distributions (Johnson and Kotz, 1969). The negative hypergeometric is an under-appreciated distribution [e.g. Equation (11.23) of Durbin et al., 1998] which bears the same relation to the hypergeometric as the negative binomial does to the binomial distribution. The multivariant generalization appears in this case as the combination of a Dirichlet and a multinomial. Confusingly, the negative hypergeometric distribution is sometimes called the inverse hypergeometric, an entirely different distribution, and vice versa.
The probability of independently sampling two count vectors, n1 and n2, from the same undetermined background is
 | (8) |
 | (9) |
If both count vectors contain only a single observation, then this profile–profile score reduces to a pairwise substitution matrix. Note, given that 
and 
(where xj is a Kronecker delta function), then all but the jth element of the product Z(xj + k)/Z(k) cancels. Thus,
 | (10) |
Applying Equation (10) to the 20-component Dirichlet mixture, dist.20comp generates the pairwise substitution matrix illustrated in Figure 3.
An interesting feature of this model is that it provides a unified homology match score for sequence–sequence, sequence–profile and profile–profile alignment [Equation (9)]. As far as we are aware, this profile–profile score has not been evaluated in a profile–profile alignment algorithm, although it is a natural generalization of the established hidden Markov model profile–sequence score. However, in the large sample limit, Equation (9) reduces to the Jensen–Shannon divergence between the two empirical amino acid distributions, a measure that has shown some promise in profile–profile alignment (Yona and Levitt, 2002; Edgar and Sjölander, 2004; Marti-Renom et al., 2004).
|
Acknowledgments |
|---|
We would like to thank Richard E. Green, Marcus Zachariah, Robert C. Edgar and Emma Hill for helpful discussions and suggestions. This work was supported by the National Institutes of Health (1-K22-HG00056). GEC received funding from the Sloan/DOE postdoctoral fellowship in computational molecular biology. SEB is a Searle Scholar (1-L-110).
Received on October 2, 2004; revised on October 25, 2004; accepted on October 25, 2004
|
REFERENCES |
|---|
|
TOPAbstractINTRODUCTIONMATHEMATICAL DETAILSREFERENCES |
|---|
Altschul, S.F. (1991) Amino acid substitution matrices from an information theoretic perspective. J. Mol. Biol., 219, 555–565[CrossRef][Web of Science][Medline].
Altschul, S.F. (1993) A protein alignment scoring system sensitive at all evolutionary distances. J. Mol. Evol., 36, 290–300[CrossRef][Web of Science][Medline].
Benner, S.A., Cohen, M.A., Gonnet, G.H. (1994) Amino acid substitution during functionally constrained divergent evolution of protein sequences. Protein Eng., 7, 1323–1332
Brenner, S.E., Chothia, C., Hubbard, T.J.P. (1998) Assessing sequence comparison methods with reliable structurally identified distant evolutionary relationships. Proc. Natl Acad. Sci. USA, 95, 6073–6078
Brenner, S.E., Koehl, P., Levitt, M. (2000) The ASTRAL compendium for protein structure and sequence analysis. Nucleic Acids Res., 28, 254–256
Bruno, W.J. (1996) Modeling residue usage in aligned protein sequences via maximum likelihood. Mol. Biol. Evol., 13, 1368–1374[Abstract].
Crooks, G.E. and Brenner, S.E. (2004a) Protein secondary structure: entropy, correlations and prediction. Bioinformatics, 20, 1603–1611
Crooks, G.E. and Brenner, S.E. (2004b) Measurements of protein sequence-structure correlations. Proteins, 57, 804–810[CrossRef][Web of Science][Medline].
Crooks, G.E., Hon, G., Chandonia, J.M., Brenner, S.E. (2004) WebLogo: a sequence logo generator. Genome Res., 14, 1188–1190
Dayhoff, M.O., Eck, R.V., Park, C.M. (1972) A model of evolutionary change in proteins. Atlas Protein Sequences Structure, 5, 89–99.
Dayhoff, M.O., Schwartz, R.M., Orcutt, B.C. (1978) A model of evolutionary change in proteins. Atlas Protein Sequences Structure, 5, (Suppl 3), 345–352.
Dimmic, M.W., Mindell, D.P., Goldstein, R.A. (2000) Modeling evolution at the protein level using an adjustable amino acid fitness model. Pac. Symp. Biocomput., 18–29.
Durbin, R., Eddy, S., Krogh, A., Mitchison, G. (1998) Biological Sequence Analysis. , Cambridge Cambridge University Press.
Edgar, R.C. and Sjölander, K. (2004) A comparison of scoring functions for protein sequence profile alignment. Bioinformatics, 20, 1301–1308
Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol., 17, 368–376[CrossRef][Web of Science][Medline].
Goldman, N., Thorne, J.L., Jones, D.T. (1996) Using evolutionary trees in protein secondary structure prediction and other comparative sequence analyses. J. Mol. Biol., 263, 196–208[CrossRef][Web of Science][Medline].
Goldman, N., Thorne, J.L., Jones, D.T. (1998) Assessing the impact of secondary structure and solvent accessibility on protein evolution. Genetics,, 149, 445–458
Gonnet, G.H., Cohen, M.A., Benner, S.A. (1992) Exhaustive matching of the entire protein sequence database. Science, 256, 1443–1445
Green, R.E. and Brenner, S.E. (2002) Bootstrapping and normalization for enhanced evaluations of pairwise sequence comparison. Proc. IEEE, 90, 1834–1847[CrossRef].
Henikoff, J.G., Greene, E.A., Pietrokovski, S., Henikoff, S. (2000) Increased coverage of protein families with the blocks database servers. Nucleic Acids Res., 28, 228–230
Henikoff, S. and Henikoff, J.G. (1992) Amino-acid substitution matrices from protein blocks. Proc. Natl Acad. Sci. USA, 89, 10915–10919
Johnson, N.L. and Kotz, S. Discrete Distributions, (1969) , New York John Wiley.
Jones, D.T., Taylor, W.R., Thornton, J.M. (1992) The rapid generation of mutation data matrices from protein sequences. Comput. Appl. Biosci., 8, , pp. 275–282
Karplus, K. (1995) Regularizers for estimating distributions of amino acids from small samples. , Santa Cruz Technical report University of California.
Koshi, J.M. and Goldstein, R. (1998) Models of natural mutations including site heterogeneity. Proteins, 32, 289–295[CrossRef][Web of Science][Medline].
Lartillot, N. and Philippe, H. (2004) A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process. Mol. Biol. Evol., 21, 1095–1109
Marti-Renom, M.A., Madhusudhan, M.S., Sali, A. (2004) Alignment of protein sequences by their profiles. Protein Sci., 13, 1071–1087[CrossRef][Web of Science][Medline].
Müller, T., Spang, R., Vingron, M. (2002) Estimating amino acid substitution models: a comparison of Dayhoff's estimator, the resolvent approach and a maximum likelihood method. Mol. Biol. Evol., 19, 8–13
Müller, T. and Vingron, M. (2000) Modeling amino acid replacement. J. Comput. Biol., 7, 761–776[CrossRef][Web of Science][Medline].
Murzin, A.G., Brenner, S.E., Hubbard, T., Chothia, C. (1995) SCOP: a structural classification of proteins database for the investigation of sequences and structures. J. Mol. Biol., 247, 536–540[CrossRef][Web of Science][Medline].
Park, J., Karplus, K., Barrett, C., Hughey, R., Haussler, D., Hubbard, T., Chothia, C. (1998) Sequence comparisons using multiple sequences detect three times as many remote homologues as pairwise methods. J. Mol. Biol., 284, 1201–1210[CrossRef][Web of Science][Medline].
Schneider, T.D. and Stephens, R.M. (1990) Sequence logos: a new way to display consensus sequences. Nucleic Acids Res., 18, 6097–6100
Sjölander, K., Karplus, K., Brown, M., Hughey, R., Krogh, A., Mian, I.S., Haussler, D. (1996) Dirichlet mixtures: a method for improving detection of weak but significant protein sequence homology. Comput. Appl. Biosci., 12, 327–345
Smith, T.F. and Waterman, M.S. (1981) Identification of common molecular subsequences. J. Mol. Biol., 147, 195–197[CrossRef][Web of Science][Medline].
Thorne, J.L., Goldman, N., Jones, D.T. (1996) Combining protein evolution and secondary structure. Mol. Biol. Evol., 13, 666–673[Abstract].
Topham, C.M., Srinivasan, N., Blundell, T.L. (1997) Prediction of the stability of protein mutants based on structural environment-dependent amino acid substitution and propensity tables. Protein Eng., 10, 7–21
Yona, G. and Levitt, M. (2002) Within the twilight zone: a sensitive profile-profile comparison tool based on information theory. J. Mol. Biol., 315, 1257–1275[CrossRef][Web of Science][Medline].
Zachariah, M.A., Crooks, G.E., Holbrook, S.R., Brenner, S.E. (2005) A generalized affine gap model significantly improves protein sequence alignment accuracy. Proteins, 58, 329–338[CrossRef][Web of Science][Medline].
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  S. Q. Le, N. Lartillot, and O. Gascuel Phylogenetic mixture models for proteins Phil Trans R Soc B, December 27, 2008; 363(1512): 3965 - 3976. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 L. Si Quang, O. Gascuel, and N. Lartillot Empirical profile mixture models for phylogenetic reconstruction Bioinformatics, October 15, 2008; 24(20): 2317 - 2323. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 S. Blanquart and N. Lartillot A Site- and Time-Heterogeneous Model of Amino Acid Replacement Mol. Biol. Evol., May 1, 2008; 25(5): 842 - 858. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
G. E. Crooks, R. E. Green, and S. E. Brenner Pairwise alignment incorporating dipeptide covariation Bioinformatics, October 1, 2005; 21(19): 3704 - 3710. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||