README.txt
Ziheng Yang, January 2004

(A)  Overview

This folder includes the files used in Yang (2004) and this file
describes the use of the AHRS algorithm described in that paper.  Two
options (clock = 5 and 6) are implemented in baseml and codeml.  clock
= 5 is for the global clock, and clock = 6 implements the heuristic
rate smoothing (AHRS) algorithm combined with maximum likelihood
estimation of divergence dates.  The mouse lemur files are used as an
example.  The control files for the nucleotide, amino acid, and codon
based analysis published in Yang (2004) are baseml2.ctl, aaml2.ctl,
and codonml2.ctl.  You can open those files to check them.

The models in Yang (2004) have two features compared with the previous
models (clock = 1, 2, 3), which were described in Yang and Yoder
(2003).  The first is that the new models deal with multiple-loci data
sets in which some species are missing at some loci.  The second is
that the AHRS algorithm attempts to automatically assign branches to
rate groups, which are then used for maximum likelihood estimation of
divergence dates using the method of Yang and Yoder (2003).  So the
new global-clock option (clock = 5) is identical to the old global
clock option (clock = 3) if every locus has all the species.  Note
however that there are some differences in the data formats between
the new models and the old ones, due to the need to deal with missing
species at some loci.  See section (B) for details.


(B) Data Format

Have a look at the sequence data file MouseLemurs123.nuc and the tree
file MouseLemurs.trees in this folder.

The new models, clock = 5 (global clock or simply clock) and 6 (local
clock) are for combined analysis of data from multiple loci (or site
partitions).  An important difference from the old options, clock = 0,
1, 2, 3, is that some species can be missing at some loci.  Specify
the number of loci using the variable ndata in the control file
baseml.ctl or codeml.ctl.  Prepare a tree file, including a tree for
all species that occur in the sequence data.  This will be called the
species tree or master tree.  The beginning of the tree file should
have the number of species and the number of trees (which is usually
1) on the first line.  Also prepare a sequence data file, with the
sequence alignments listed one after another for the loci.  The
example file MouseLemurs123.nuc lists data from the 3 codon
positions as if they were three different genes.  When the program
runs, it reads the master tree first, and then the sequence alignments
locus by locus.  For each locus, it constructs (extracts) the sub-tree
from the master tree.

For example, suppose the master tree is for six species with two
fossil calibration points:

      (((S1, S2), S3) @0.10, ((S4, S5) @0.20, S6));

Suppose we have sequences for 4 species at the first locus and 5
species at the second locus (with ndata = 2).  The data file looks as
follows.

     4 100
     S1       TCTATGTTATATGTATGAGTATGA ...
     S3       TCGATATTATGTGTATGAGTATGA ...
     S4       TCTATATTACATGTATGAGTATGA ...
     S6       TCCATATTAAATGTATGAGTATGA ...

     5 200
     S1       GTTATATGTATGAGTATGATCTAT ...
     S2       GTTATATGCGTGAGTATGAACTAT ...
     S4       ATTATGTGTATGAGTATGATCGAT ...
     S5       GCTACATGTATGAGTATGATCTAT ...
     S6       ATTAAATGTATGAGTATGATCCAT ...


The program will then construct the "gene tree" for locus 1 to be 
      ((S1, S2), S3) @0.10;

You should have at least one calibration node for each locus.  If two
"sister" species never occur simultaneously at one locus so that there
is no way of identifying the date of their divergence, you might want
to consider them as one "species".  For example, if donkey is
sequenced at some loci and horse at others and at no locus both are
sequenced, you can rename the two species equus in the tree and
sequence data files.


(C) Models clock = 5 (clock) and clock = 6 (local clock)

When you run clock = 5, a global clock is assumed.  If all species are
present in every locus, this analysis can be achieved by using clock =
3 (combined analysis) as in Yang & Yoder (2003).  This is the case for
the included example data file MouseLemurs123.nuc.  Note again however
that the data formats are different.  clock = 3 uses MouseLemurs.nuc
while clock = 5 uses MouseLemurs124.nuc.  If some species are missing
at some loci in your data, you can't use clock = 3 anymore and clock =
5 is the only choice.

The AHRS algorithm (clock = 6) works in 3 steps.  Here is a brief
description.  See Yang (2004) for more details.

   Step 1: Estimate branch lengths on each gene tree without the
           clock.

   Step 2: Heuristic rate smoothing, to estimate one set of divergence
           times for the nodes in the species tree and as many sets of
           branch rates as the number of loci.  Each branch in the
           gene tree is assigned and estimated one rate.  For the
           above toy example, the program estimates 3 divergence times
           (five ancestral nodes minus two calibration nodes), plus 6
           rates (six branches in a rooted tree of four species) for
           the first locus and 8 rates for the second locus (eight
           branches in a rooted tree of five species).  A smoothing
           parameter called \nu is also estimatd for each locus.  Rate
           estimation is achieved by minimizing the discrepancy
           between the predicted branch lengths and the estimates
           obtained in Step 1, and by minimizing rate changes over
           time (across lineages).  The estimated rates for branches
           in the gene trees are then collapsed into a few branch rate
           categories.  The number of rate categories is set by a
           variable called nbrate=4 in the routine DatingHeteroData()
           in the file treesub.c.  You can change this and recompile.
           The program also allows you to manually classify the rates
           into groups.  When it asks for the number of rate groups,
           you can reply with 0: meaning no change, 1: meaning one
           rate (same as using clock = 5), or say 5 (in which case the
           program will further ask to input 4 cutting points to
           separate the rates into 5 groups).

   Step 3: Maximum likelihood estimation of divergence times using the
           assigned branch rate groups.  Divergence times are estimated
           simultaneously with the rates for branch rate groups.  For
           the toy example mentioned above, this estimates 3 = 5 - 2
           node ages in the species tree, plus 3 rates at each locus,
           with 9 time and rate parameters in all.


(D) Output

The output is fairly self-explanatory.  For clock = 6, look at the
output on the screen at each of the three steps:

   Step 1: Estimating branch lengths under no clock.
   Step 2: Ad hoc rate smoothing to estimate branch rates.
   Step 3: ML estimation of times and rates.

In particular, read in the output at the end of step 2 and check the
branch labels (say, using Rod Page's TreeView) and see whether the
assignment is reasonable and whether you want to use your own
classification scheme.  The programs print out a file named
RateDist.txt, which contains a distance matrix for each locus, with
the distance calculated as the difference between two rates.  You can
use an algorithm such as UPGMA to cluster the rates into groups
(clades) to help with the classification.  I used the neighbor program
in J. Felsenstein's phylip package to do this (figure 5b and the
4-rate manual (4RM) model in Yang 2004).


(E)  Notes

 1.  The example control files are baseml2.ctl, aaml2.ctl, codonml2.ctl.

 2.  In early versions of baseml and codeml, the order of the control
     variables in the control files does not matter to the
     specification.  Unfortunately that is not the case anymore for
     clock = 5 or 6.  The variables clock and ndata should be
     specified in the control file before some other variables such as
     fix_kappa, kappa, fix_alpha, alpha etc.  See the example files
     baseml2.ctl and codonml2.ctl etc. included in the same folder.
     The reason for this is that you can fix the kappa values for
     different genes at certain values during Step 1, and for the
     program to know how many kappa values to read, it needs to know
     the number of genes (ndata); see note 5 below.

 3.  The options are implemented in baseml for nucleotide sequences
     and in codeml for amino acid and codon sequences.  

 4.  Only rooted binary trees are accepted in the tree file.

 5.  Estimating and fixing substitution parameters in step 1.  If the
     model involves substitution parameters such as kappa, alpha, you
     can run the program multiple times to guarantee getting correct
     estimates for them Step 1 of the algorithm and then use them as
     fixed constants in both Step 1 (when you rerun the program) and
     Step 3.  (Step 2 does not use a substitution model and so those
     parameters are irrelevant.)  Step 1 in the clock = 6 model always
     uses the algorithm of iterating one branch at a time, and may be
     problematic if the gamma shape parameter (for rate variation
     among sites) is estimated at the same time.  Note that the
     calculation in step 1 is just the old clock = 0 analysis and is
     in older versions of the programs, if the coccrect tree files are
     supplied for each locus.  In that case you can use clock = 0
     method = 0 to try and estimate the substitution parameters and
     then fix them when you run clock = 6.  If you have different
     numbers of species at different loci, the tree for each locus is
     different and it would be awkward to use clock = 0 to complete
     step 1.  In that case you can use clock = 6 so that the correct
     trees are extracted and used for each locus, but you run the
     algorithms multiple times to make sure the correct MLEs are
     obtained.  Either way, when step 1 is successful, you can fix the
     parameters in the control file, by say,

                  data = 3
             fix_alpha = 1
                 alpha = 0.29169  0.16392  1.24726

     which means that the alpha parameter is fixed at the three values
     above for the three genes (codon positions in this case).  Other
     parameters that can be fixed this way and that can vary for site
     partitions or genes or proteins include kappa and alpha for
     baseml, alpha and aaRatefile for amino acid sequences, and icode
     (genetic code table) kappa, and omega for codon-based analsis.
 
 6.  I notice that minor differences in the rate estimates might lead
     to assignment of one or two branches in a different rate group,
     leading to differences in final time estimates.  As a result,
     multiple runs using the same setting may lead to slightly
     different time estimates.  You should watch out for failure of
     the iteration algorithm in any of the three steps, which can
     cause differences between runs.

 7.  For protein sequences, you can use different substitution rate
     matrices.  In the following, three nuclear proteins (using wag)
     are followed by one mitochondrial protein (using mtmam).

        aaRatefile = wag.dat wag.dat wag.dat mtmam.dat

     For codon sequences, you can use different genetic codes.  In the
     following, three nuclear genes are followed by one mitochondrial
     gene.

             icode = 0 0 0 1

 8.  Separating codon positions into different files.  You can use
     Mgene = 1 option in baseml to separate the three codon positions
     into three different files (named Gene1.seq, Gene2.seq, etc.).
     You need to use verbose = 1 to do this, I think.  Then you can
     copy and paste those files into one sequence data file, which
     will be in the format required by the options here.

Good luck.

References

Yang, Z. 2004. A heuristic rate smoothing procedure for maximum likelihood estimation of species divergence times. Acta Zoologica Sinica 50:645-656.

//end of file
