© Copyright 1986-2008 by the University of Washington. Written by Joseph Felsenstein. Permission is granted to copy this document provided that no fee is charged for it and that this copyright notice is not removed.
This program infers an unrooted phylogeny from protein sequences, using a new method intermediate between the approaches of Eck and Dayhoff (1966) and Fitch (1971). Eck and Dayhoff (1966) allowed any amino acid to change to any other, and counted the number of such changes needed to evolve the protein sequences on each given phylogeny. This has the problem that it allows replacements which are not consistent with the genetic code, counting them equally with replacements that are consistent. Fitch, on the other hand, counted the minimum number of nucleotide substitutions that would be needed to achieve the given protein sequences. This counts silent changes equally with those that change the amino acid.
The present method insists that any changes of amino acid be consistent with the genetic code so that, for example, lysine is allowed to change to methionine but not to proline. However, changes between two amino acids via a third are allowed and counted as two changes if each of the two replacements is individually allowed. This sometimes allows changes that at first sight you would think should be outlawed. Thus we can change from phenylalanine to glutamine via leucine in two steps total. Consulting the genetic code, you will find that there is a leucine codon one step away from a phenylalanine codon, and a leucine codon one step away from glutamine. But they are not the same leucine codon. It actually takes three base substitutions to get from either of the phenylalanine codons TTT and TTC to either of the glutamine codons CAA or CAG. Why then does this program count only two? The answer is that recent DNA sequence comparisons seem to show that synonymous changes are considerably faster and easier than ones that change the amino acid. We are assuming that, in effect, synonymous changes occur so much more readily that they need not be counted. Thus, in the chain of changes TTT (Phe) --> CTT (Leu) --> CTA (Leu) --> CAA (Glu), the middle one is not counted because it does not change the amino acid (leucine).
To maintain consistency with the genetic code, it is necessary for the program internally to treat serine as two separate states (ser1 and ser2) since the two groups of serine codons are not adjacent in the code. Changes to the ser1 some overlaying of changes.
One problem that can arise is that two or more of the species can be so dissimilar that the distance between them would have to be infinite, as the likelihood rises indefinitely as the estimated divergence time increases. For example, with the Jukes-Cantor model, if the two sequences differ in 75% or more of their positions then the estimate of divergence time would be infinite. Since there is no way to represent an infinite distance in the output file, the program regards this as an error, issues an error message indicating which pair of species are causing the problem, and stops. It might be that, had it continued running, it would have also run into the same problem with other pairs of species. If the Kimura distance is being used there may be no error message; the program may simply give a large distance value (it is iterating towards infinity and the value is just where the iteration stopped). Likewise some maximum likelihood estimates may also become large for the same reason (the sequences showing more divergence than is expected even with infinite branch length). I hope in the future to add more warning messages that would alert the user the this.
If the similarity table is selected, the table that is produced is not in a format that can be used as input to the distance matrix programs. It has a heading, and the species names are also put at the tops of the columns of the table (or rather, the first 8 characters of each species name is there, the other two characters omitted to save space). There is not an option to put the table into a format that can be read by the distance matrix programs, nor is there one to make it into a table of fractions of difference by subtracting the similarity values from 1. This is done deliberately to make it more difficult for to use these values to construct trees. The similarity values are not corrected for multiple changes, and their use to construct trees (even after converting them to fractions of difference) would be wrong, as it would lead to severe conflict between the distant pairs of sequences and the close pairs of sequences.
The constants that are available to be changed by the user at the beginning of the program include "maxcategories", the maximum number of site categories, "iterations", which controls the number of times the program iterates the EM algorithm that is used to do the maximum likelihood distance, "namelength", the length of species names in characters, and "epsilon", a parameter which controls the accuracy of the results of the iterations which estimate the distances. Making "epsilon" smaller will increase run times but result in more decimal places of accuracy. This should not be necessary.
The program spends most of its time doing real arithmetic. The algorithm, with separate and independent computations occurring for each pattern, lends itself readily to parallel processing.
5 13 Alpha AACGTGGCCACAT Beta AAGGTCGCCACAC Gamma CAGTTCGCCACAA Delta GAGATTTCCGCCT Epsilon GAGATCTCCGCCC |
(Note that when the options for displaying the input data are turned off, the output is in a form suitable for use as an input file in the distance matrix programs).
Nucleic acid sequence Distance Matrix program, version 3.69
5 species, 13 sites
F84 Distance
Transition/transversion ratio = 2.000000
Name Sequences
---- ---------
Alpha AACGTGGCCA CAT
Beta ..G..C.... ..C
Gamma C.GT.C.... ..A
Delta G.GA.TT..G .C.
Epsilon G.GA.CT..G .CC
Empirical Base Frequencies:
A 0.24615
C 0.36923
G 0.21538
T(U) 0.16923
5
Alpha 0.000000 0.303900 0.857544 1.158927 1.542899
Beta 0.303900 0.000000 0.339727 0.913522 0.619671
Gamma 0.857544 0.339727 0.000000 1.631729 1.293713
Delta 1.158927 0.913522 1.631729 0.000000 0.165882
Epsilon 1.542899 0.619671 1.293713 0.165882 0.000000
|