import paml4.8

[paml.git] / doc / pamlHistory.txt

History and bug fixes of PAML (Phylogenetic Analysis by Maximum Likelihood)\r
\r
Ziheng Yang\r
\r
\r
If you find problems or have questions, please visit the PAML\r
discussion site\r
https://groups.google.com/forum/#!forum/pamlsoftware.\r
\r
\r
\r
\r
Version 4.8a, August 2014\r
\r
(*) baseml (nhomo = 3 or 4).  I have added the nonhomogeneous GTR\r
(REV) model with a set of abcde parameters for every branch.  The\r
specifications are as follows.\r
\r
model=7 nhomo=4 fix_kappa=0: one set of rate parameters abcde in GTR for the whole tree.\r
model=7 nhomo=4 fix_kappa=1: one set of rate parameters abcde in GTR for every branch.\r
\r
In version 4.8, only the second option (one set of abcde parameters\r
for the whole tree) is implemented, no matter whether you choose 0 or\r
1 for fix_kappa.  After this change, the use of fix_kappa in the nhomo\r
models is consistent between HKY and GTR.  For HKY, the specifications\r
have always been as follows.\r
\r
model=4 nhomo=4 fix_kappa=0: one kappa in HKY for the whole tree.\r
model=4 nhomo=4 fix_kappa=1: one kappa in HKY for every branch.\r
 \r
Note that nhomo=4 means one set of frequency parameters for every\r
branch on the tree.\r
\r
I have also added a method for calculating the expected number of\r
changes from any nucleotide i to nucleotide j along every branch under\r
the nonhomogeneous models (nhomo=3 or 4).  This accounts for the fact\r
that the base frequencies are changing over time, and so does the rate\r
of substitution.\r
\r
I tested only the nhomo=4 option carefully and the options nhomo=3 and\r
5 are not well tested.\r
\r
\r
\r
Version 4.8, March 2014\r
\r
(*) mcmctree.  The prior for rates for loci (with ndata = 2 or more)\r
has changed in this version, and the gamma-Dirichlet prior is used.\r
The old i.i.d. prior is replaced and unavailable anymore.  The new\r
gamma-Dirichlet prior is decribed in a paper \r
dos Reis M, Zhu T, Yang Z. 2014. The impact of the rate prior on\r
Bayesian estimation of divergence times with multiple loci. Syst\r
Biol. in press.\r
\r
(*) codeml.  Added a Bayesian method for estimating distance t and dN/dS\r
ratio w in pairwise comparisons.  The option is \r
      runmode = -3  * -3: pairwise Bayesian\r
      runmode = -3 1.1 1.1 1.1 2.2   * -3: pairwise Bayesian\r
The four numbers are the parameters in the gamma priors for t and w (a_t b_t a_w b_w).\r
The method is described in a manuscript\r
Angelis K, dos Reis M, Yang Z. 2014. Bayesian estimation of nonsynonymous/synonymous rate \r
ratios for pairwise sequence comparisons.  Mol Biol Evol.\r
\r
(*) baseml, codeml, mcmctree).  In version 4.7a, I disabled the option\r
of reading species indexes in the tree file (for example, the number\r
19, instead of the sequence name, may be used in the tree file to mean\r
the 19th sequence in the sequence alignment).  I have added this back.\r
Version 4.7 does not have this problem.\r
\r
(*) codeml in paml 4.6 and 4.7 has a bug if your sequences have stop\r
codons.  Stop codons have always been disallowed in the codon models\r
implemented in codeml.  In versions 4.6-4.7, if there is a stop codon\r
in any sequence, the codons in all sequences at that site (column in\r
alignment) are changed into ??? and then the sequences are processed\r
in the usual way.  The likelihood should then be the same as the\r
likelihood if that column is removed from the alignment.  However, the\r
implementation of the idea was incorrect and gives wrong results if\r
the sequences do not contain any ambiguity characters and if you\r
choose cleandata = 0.  I have now fixed this.  The idea is not really\r
useful, so the good advice is that you delete the column of stop\r
codons before the data are analyzed using codeml.\r
 \r
\r
Version 4.7a, October 2013\r
\r
(*) codeml.  The mutation-selection model of Yang and Nielsen (2008)\r
is noted to work for the site and branch-site models but not for the\r
branch models.  This is now fixed.\r
\r
(*) mcmctree.  I fixed and modified the infinitesites program, which\r
generates the limiting distribution when the amount of sequence data\r
approaches infinity.  This works for both the clock (clock=1) and\r
relaxed clock models (clock=2 or 3).  The mcmctree tutorial explains\r
how to run this program.\r
\r
\r
Version 4.7, January 2013\r
\r
(*) mcmctree.  This version includes a method for dating divergences\r
using sampled tips, as in viral datasets.  The prior for node ages in\r
the given rooted tree is generated using a\r
birth-death-sequential-sampling model, described in Stadler and Yang\r
(2012 Syst Biol, submitted).  The SIV/HIV-2 example dataset analyzed\r
in the paper is in the examples/TipDate/ folder.  Look at the readme\r
file there.\r
\r
(*) codeml.  Version 4.6 crashes if you have a lot of ambiguity codons\r
in the alignment, or precisely if the total number of fully determined\r
codons (61 in the standard code) and the ambiguous codons exceed 127.\r
I was using "unsigned char" to represent the codons, which go from 0\r
to 255 and can represent 256 distinct codons.  Somehow I changed this\r
to "char", which on some systems goes from -128 to 127.  Can't\r
remember why I made the change.  Anyway this is now fixed.  Versions\r
4.5 or earlier do not have this problem.  Thanks to David Yu for\r
pointing out the problem.\r
\r
\r
\r
Version 4.6, September 2012\r
\r
(*) baseml. The model of autocorrelated-rates among sites (Yang 1995\r
Genetics 139:993-1005) is broken in paml versions 4.3 to 4.5\r
inclusive.  Earlier versions were apparently fine.  When you select\r
the auto-discrete-gamma model (alpha > 0, rho > 0) or the\r
nonparametric autocorrelated-rates model (fix_alpha = 0, ncatG = 3,\r
nparK = 4), the program aborts with the following error message:\r
\r
"Error: use 10, 20, 32, 64, 128, 512, 1024 for npoints for legendre.."\r
\r
This is now fixed.  The option of nparK = 3 seems to have problems in\r
all versions, and the program prints out the following warning\r
message.  \r
\r
"nparK==3, double stochastic, may not work.  Use nparK=4?"\r
\r
This option is supposed to fit a Markov-dependent model with equal\r
proportions in the rate categies.  I have not got time to investigate\r
this option, but a glance at table 3 in that paper indicates that this\r
model was not used in the table and probably never properly\r
implemented, so please follow the advice and use nparK = 4 instead.\r
\r
(*) codeml.  I have now added the site model M2a_rel of Weadick & Chang\r
(2012 Mol. Biol. Evol.), in which w2 > 0.  This is specified using\r
NSsites = 22 while model M2a is specified as NSsites = 2 (with w2 >\r
1).  M2a_rel is the null model for the likelihood ratio test based on\r
clade model C, suggested by Weadick & Chang (2012 Mol. Biol. Evol.).\r
Note that the very old site model M2 from Nielsen & Yang (1998\r
Genetics) assumes w0 = 0.  Yang, Wong, Nielsen (2005 MBE) changed the\r
model so that 0 < w0 < 1 and renamed the model M2a.  My tests using\r
the versions around that time reveals the following:\r
\r
3.13d (M2):  w0 = 0, w2 > 0\r
3.14a (M2a): w0 > 0, w2 >= 1\r
3.14b (M2a): w0 > 0, w2 >= 1\r
\r
If seems that when we changed M2 into M2a, we also changed the\r
constraint on w2.  I don't have a copy of version 3.14 to check.\r
\r
(*) basseml & codeml.  To following option\r
    bootstrap = 1000\r
\r
should produce 1000 bootstrap samples in the file boot.txt.  A bug is\r
introduced in versions 4.4 and 4.5, so that garbiage is created\r
instead of bootstrap datasets.  Versions 4.3 and earlier were correct.\r
This bug is now fixed.\r
\r
(*) evolver.  The simulation of codon sequences (evolver 6\r
MCcodon.dat) is incorrect when a genetic code different from the\r
universal code is used.  This is a conceptual error I made and it has\r
been in all previous versions when simulation of codon sequences is\r
allowed.  I have assumed incorrectly that setting the frequencies of\r
stop codons to 0 should be sufficient to get correct results.  The\r
zero frequencies allow the program to identify the stop codons\r
correctly, but one still needs the code table to know which changes\r
are synonymous and which nonsynonymous.  evolver has been using the\r
universal code to simulate codon sequences.  Simulation using the\r
universal code has been correct, and simulation using the\r
mitochondrial code (for example) is affcted and the results should be\r
incorrect although the impact may be slight.  This bug is now fixed.\r
A variable is added in the control file MCcodon.dat right below the\r
codon frequencies, which tells the program which genetic code to use.\r
\r
1   * genetic code (0:universal; 1:mammalian mt; 2-10:see below)\r
\r
The program also checks that the stop codons have frequency 0 and that\r
the frequencies of the sense codons sum to 1.  Thanks to Mario dos\r
Reis for pointing this out.\r
\r
(*) mcmctree.  A bug is found in the calculation of the prior on rates\r
under the correlated-rates model (clock=3).  This bug has been in\r
version 4 (July 2007) to version 4.5 (December 2011), inclusive.  The\r
prior is calculated by multiplying the conditional densities of\r
equation (7) in Rannala & Yand (2007) over all interior nodes on the\r
master tree (see also Figure 1 in that paper).  However, the rA used\r
in the program was not for the current branch: it was instead for the\r
ancestral branch.  The correct formula used was correct for the root\r
of the master tree, but incorrect for other interior nodes.  The\r
impact of the bug on posterior time estimates may be expected to be\r
small, but if the clock is seriously violated and rates change over\r
neighboring branches, this may be a concern.\r
\r
\r
\r
Version 4.5, December 2011\r
\r
baseml.  I have added the joint ancestral reconstruction under the\r
nonhomogeneous models (nhomo = 3 or 4).  This works for the model of\r
one rate for all sites, and does not work for the model of gamma rates\r
for sites or the partition models.  Also I implemented the joint\r
reconstruction only and not the marginal reconstruction.\r
\r
mcmctree.  The program now prints out a file called FigTree.tre in the\r
current working directory that is readable by FigTree.  The branch\r
lengths show the posterior means of times while the bars generated by\r
FigTree will indicate the 95% CIs.  I have also added an option to\r
deal with TipDate data, such as viral sequences with dates.  See\r
examples/TipDate/README.txt for more notes.\r
\r
\r
\r
Version 4.4e, May 2011\r
\r
mcmctree.  I have added an option of automatically adjusting the step\r
lengths for proposals in the MCMC algorithm.  The program will use the\r
burn in to automatically adjust the step lengths, using those values\r
specified in the control file as initial values.\r
\r
codeml: The iteration method that updates one branch length at a time\r
(method = 1) is broken for the free-ratios model (model = 1) from\r
versions 4.3 to 4.4d inclusive.  When codeml runs, you will see\r
messages like the following on the screen: \r
Warning rounding error? b=3 cycle=0 lnL=6036.9413757 != 6070.1069284 \r
Version 4.2 and earlier are correct.  This bug is now fixed.  \r
Thanks to Tae-Kun Seo for reporting the bug.\r
\r
\r
\r
Version 4.4d, March 2011\r
\r
yn00.  The option of using weighting (weighting = 1) in the\r
calculation of dN and dS by the YN00 method (Yang and Nielsen 2000)\r
has been broken since Version 4.4 (February 2010).  A bug was\r
introduced when I changed the character coding in Version 4.4, so that\r
this bug affects Versions 4.4, 4.4a, 4.4b, and 4.4c.  The bug\r
typically causes a crash.  The result for weighting = 0 is not\r
affected.  Versions before 4.4 are correct.  This bug is now fixed.\r
Thanks for Charlie\r
(http://www.ucl.ac.uk/discussions/viewtopic.php?t=8726&highlight=yn00).\r
\r
\r
\r
Version 4.4c, August 2010\r
\r
mcmctree.  A bug causes the exotic calibrations (gamma, skew normal\r
and skew t) to be sometimes ignored, if you have multiple loci\r
(partitions) and if some species/sequences are missing at some loci.\r
If all calibrations are bounds (minimum, maximum, and joint), or if\r
the partitions have the same number of species/sequences, the\r
calculation is still correct.  This bug is in versions 4.1-4.4b.  This\r
is now fixed.  Thanks to Mike Steiper.\r
\r
\r
\r
Version 4.4b, July 2010\r
\r
mcmctree.  A serious bug has been discovered which affects the option\r
of approximate likelihood calculation (usedata = 2).  If the first\r
child of the root in the tree for any locus is a tip, the MLEs of\r
branch lengths in the in.BV file are assigned to wrong branches, so\r
that the results are entirely incorrect, and also very different from\r
those obtained using the exact likelihood calculation (usedata = 1).\r
Suppose the tree for a locus is (A, B), with a bifurcation at the\r
root, and with A and B representing species or clades.  If A is a tip,\r
the program is incorrect, while if A is an internal node, the program\r
is correct.  This bug is in versions 4.2b - 4.4a.  Hopefully this is\r
now corrected.  Thanks to Mario dos Reis.\r
\r
mcmctree.  I also changed the specification on L, U, and B bounds to\r
allow the user to specify the tail probabilities (which used to be\r
fixed at 2.5%) for each fossil calibration.  See table 8 in\r
pamlDOC.pdf for details.  You will still have to run the program\r
without data to get the effective prior of times used by the program.\r
\r
\r
\r
Version 4.4a, June 2010\r
\r
codeml (RateAncestor=1).  For ancestral reconstruction under codon\r
models, the program prints out the translated amino acids when it\r
prints the codons.  A bug caused the translation to be incorrect,\r
so instead of the amino acid, the program prints something like h or\r
H.  This bug was introduced in Version 4.4, February 2010 when I\r
changed the way that the characters are coded.  This is fixed now.\r
\r
baseml (clock = 1, 2, or 3 Mgene option).  If you assign different\r
rates for different genes (or site parations) and the model assumes\r
global clock or local clock, the program fails to print out the\r
relative rates for genes.  This bug seems to be introduced between\r
3.13 and 3.14 and has since.  The log likelihood value and estimates\r
of parameters are all correct, and the problem is only that MLEs of\r
rates for genes are not printed in the output file properly.\r
\r
\r
\r
Version 4.4a, April 2010\r
\r
In versions 4.2, 4.3, 4.4, Clade model D does not print the w\r
estimates for the site classes correctly.  For the example run in\r
examples/datasets2/, the incorrect output looks like the following\r
\r
site class             0         1         2\r
proportion        0.48590   0.10718   0.40692\r
branch type 0:    0.10557   4.17186   4.17186\r
branch type 1:    0.10557   4.17186   3.05340\r
\r
while the correct results (from version 4.1) are \r
\r
site class         0         1         2\r
proportion       0.48590  0.10718  0.40692\r
background w     0.10557  4.17186  3.05340\r
foreground w     0.10557  4.17186  0.95081\r
\r
The ratios w2 and w3 for site class 2 are incorrectly printed, with w2\r
equal to w1.  The branch lengths and lnL values etc. are all correct.\r
This bug affects versions 4.2, 4.3, and 4.4, and version 4.1 and\r
earlier versions appear to be correct.  Also the bug affects clade\r
model D, and clade model C appears to be fine.  This is fixed now.\r
Thanks to cajawe for reporting the bug (see\r
http://www.ucl.ac.uk/discussions/viewtopic.php?p=28240#28240).\r
\r
\r
\r
Version 4.4, February 2010\r
\r
I changed the way that the characters are coded in the programs, which\r
means that essentially every program in the package is modified.  It\r
is likely that some bugs are introduced during the process.  \r
\r
codeml.  The iteration algorithm specified by method = 1 when applied\r
to amino acid model of gamma rates among sites (seqtype = 2, alpha >\r
0, method = 1) is broken in versions 4.2-4.3b.  The bug causes the\r
iteration to fail to converge, with messages like the following\r
printed on the monitor: \r
"Warning rounding error? b=5 cycle=0 lnL=1718.5667685 != 1749.4724545".\r
\r
mcmctree.  If the gamma prior for rates has alpha <1/3, the initial\r
values are generated incorrectly.  This is now fixed.\r
\r
\r
\r
Version 4.3b, November 2009\r
\r
mcmctree.  A bug was introduced in version 4.2, which causes the\r
program to ignore the minimum age constraint on the root age and uses\r
the maximum bound specified by RootAge in the control file, if there\r
is a calibration on the root and it is a minimum-age constraint.  If\r
you have no calibration on the root, or if the calibration is an\r
maximum bound, both minimum and maximum bound or gamma distribution,\r
there is no problem.  The error affects versions 4.2 and 4.3, and\r
version 4.1 is correct.  This is now fixed.  Thanks to Mathieu\r
Groussin.\r
\r
evolver.  I have now added back the option of printing out the site\r
pattern counts instead of the sequences (specified at the beginning of\r
the control file such as MCbase.dat).  Use of pattern counts takes\r
less disk space.  The option was removed in version 4.2b after the\r
algorithm for collaspsing sites into patterns was rewritten.\r
\r
\r
\r
Version 4.3a, October 2009\r
\r
mcmctree.  A bug was introduced in version 4.3, which causes the root\r
age to be proposed incorrectly.  The effect of the bug appears to be\r
subtle and is hard to assess.  It means that the results from all\r
three clock models (clock = 0, 1, 2) are incorrect.  Here is a bit more\r
detail about the bug.  The internal maximum node age was set too high,\r
and caused loss of significant digits in a smart reflection algorithm\r
which reflects proposed values into the feasible interval.  This\r
update fixes the bug.  Please use this version to reanalyze your data\r
if you have used version 4.3.  Version 4.2 does not have this bug.\r
\r
mcmctree.  I now supplies compiling options so that you can use hard\r
minimum hard maximum, or hard minimum soft maximum.  The executables\r
are called mcmctreeHH and mcmctreeHS.  See the commands for compiling\r
at the beginning of the file mcmctree.c.\r
\r
\r
\r
Version 4.3, August 2009\r
\r
mcmctree.  The approximate likelihood calculation was unreliable,\r
because the Hessian matrix was not calculated reliably, especially if\r
the maximum likelihood estimates of some branch lengths are 0.  I have\r
now implemented the approach of Seo, Kishino and Thorne (2004), which\r
appears to be more stable.\r
\r
baseml & codeml: I changed these two programs to read fasta alignments\r
directly, working out the number of sequences and the sequence length\r
automatically.  However, this does not work when you have multiple\r
alignments in one file, that is, if ndata = 2 or higher.\r
\r
\r
\r
Version 4.2b, April 2009\r
\r
mcmctree. There is a bug that leads to crash when usedata = 3 and when\r
the data at some loci have ambiguity characters while other loci do\r
not.  usedata = 3 means that the program invokes baseml or codeml to\r
calculate the maximum likelihood estimates (MLEs) of branch lengths\r
and the variances of those MLEs.  The temporary sequence alignment\r
files for the loci (tmp?.txt) generated by mcmctree have garbiage and\r
are unreadable by baseml (or codeml).  This bug is now corrected.\r
(You can replace the wrong alignment file for the locus by your\r
original alignment and run baseml outside mcmctree.)\r
\r
\r
All programs.  I rewrote the code for collapsing sites into patterns.\r
This should have no impact on user files, except that evolver does not\r
produce data files in the "pattern-count" format anymore.  Now two\r
formats (paml/phylip format and PAUP nexus format) are used, specified\r
by the first line of the control data file (e.g., MCbase.dat,\r
MCcodon.dat, MCaa.dat).  This pattern-count format is still accepted\r
as input by baseml, codeml, etc. (See the section on "Site pattern\r
counts" in the documentation).\r
\r
baseml & codeml.  A bug was introduced to the local clock model\r
(clock=6), so that the program generates the following error message:\r
"need fossils for this locus".  For example, running codeml on the\r
control file /examples/MouseLemurs/codonml2.ctl will do so.  This bug\r
is now fixed.\r
\r
codeml: The bug introduced in version 4.2 concerning the free-ratios\r
model (Yang 1998 MBE) was not fixed properly in version 4.2a (see the\r
notes for 4.2a below).  The results are incorrect when there are\r
exactly 7 branch types (for example when you run the free-ratios model\r
on a 5-species tree) but are correct for other numbers of branch\r
types.  I hope this is now finally fixed.\r
\r
\r
\r
Version 4.2a, February 2009\r
\r
\r
codeml: Two bugs were known to have been introduced in version 4.2\r
when I rewrote the code for the branch-site and clade models.  The\r
first is the null hypothesis in the branch-site test (model = 2\r
NSsites = 2, fix_omega = 1 omega = 1).  The reference for this model\r
is Yang, Wong & Nielsen (2005), and Zhang et al. (2005).  The\r
alternative model (fix_omega = 0) is fine.\r
\r
The second bug is in the free-ratios model (Yang 1998 MBE).\r
\r
The results from those two models are incorrect and should be ignored\r
if you used version 4.2.  Both mistakes are fixed in this update.\r
\r
Thanks to W.P. Zhou and BYLee.\r
\r
\r
\r
Version 4.2, December 2008\r
\r
MCMCtree: The implementation of the lower bound and of the upper bound\r
is changed.  The improper density of equation (15) of Yang and Rannala\r
(2006) for the lower bound is noted to push up divergence times, a\r
feature we consider undesirable.  Now a truncated Cauchy distribution\r
is used instead.  This will be described in more detail somewhere.\r
\r
The default file name for the MCMC samples is mcmc.out.  To change\r
this, add the following line in the control file mcmctree.ctl:\r
      mcmcfile = mcmc.out \r
\r
codeml: Clade models C and D are extended to accommodate more than two\r
branch types.  The old models accept two branch types only, referred\r
to as foreground and background branches.  The structure of the new\r
models are below.  You use branch/node numbers to specify the branch\r
types (with #0 to be the default, followed by #1, #2, #3, ...).  As\r
before, the BEB calculation is implemented for clade model C only, and\r
is not available for clade model D.  The BEB calculation is expensive,\r
so I expect clade model C will work if you have only a few branch\r
types.  For each additional branch type, the BEB computation is 10\r
times more expensive.\r
\r
Site class   Proportion    w ratio\r
\r
  0             p0         0 < w0 < 1\r
  1             p1         w1 (w1 = 1 is fixed in model C)\r
  2             p2         w2, w3, w4, ...\r
\r
I also added an option of fixing the last w to 1 in clade model C,\r
specified in the control file as follows.  For example, if there are\r
three branch types in the tree (see the table above), the last w fixed\r
will be w4.\r
\r
       fix_omega = 1\r
       omega = 1\r
\r
The two versions of the clade model can be compared to test whether\r
the last w is > 1.\r
\r
\r
codeml: branch model with amino acid chemical properties, specified\r
using model = 2 and aaDist = 7.  Under this model, the types of\r
branches are specified just like the branch model (by labelling\r
branches).  For each branch type, the program estimates a few w's,\r
depending on specifications in the file OmegaAA.dat.  This model seems\r
to be broken for a long time since its implementation.  Look at the\r
readme.txt file in examples/mtCDNAape/.\r
\r
\r
\r
Version 4.1, August 2008\r
\r
MCMCtree: A few changes have been introduced to this program.  The\r
format for specifying fossil calibrations in the tree file has been\r
changed, so that the gamma distribution is specified as 'G(alpha,\r
beta)' and the old format '>0.6=0.7<0.8' does not work anymore.  The\r
bounds can be specified using either the old form '<0.8' or the new\r
form 'U(0.8)'.  Two exotic distributions, skew normal and skew t, are\r
added, using the format 'SN(location, scale, shape)' and 'ST(location,\r
scale, shape, df)'.  I edited the documentation pamlDOC.pdf and\r
rewrote the part for mcmctree.  The format of the MCMC output file\r
mcmc.out has been modified so that it is readable by Andrew Rambaut's\r
Tracer program.  Note that MCMCtree starts taking samples and writing\r
into the file only after burnin.\r
\r
codeml. The mutation-selection models of Yang and Nielsen (2008\r
Mol. Biol. Evol. 25, 568-579) are specified using the options\r
\r
    CodonFreq = 6 or 7  * 6:FMutSel0, 7:FMutSel\r
      estFreq = 1\r
\r
Use examples/mtCDNAape/codeml.HC.ctl to duplicate results in table 1\r
in the paper.  It seems that some results are in the mlc file while\r
some other results are in the rst and rst1 files.\r
\r
\r
codeml (branch model): When codon sequences are analyzed under the\r
branch model (model = 1 or 2, NSsites = 0), the main output file\r
includes a tree under the line "w ratios as labels for TreeView:".\r
This line was introduced in version 3.15.  However, the branch lengths\r
in the tree had the dN values, while I intended them to be the\r
original branch lengths (which are t, the number of nucleotide\r
substitutions per codon).  This bug may also affect the results of\r
ancestral reconstruction under the same branch model.  This bug is\r
fixed now.\r
\r
\r
baseml and codeml (discrete-gamma model): My intention has been to use\r
the mean and not the median to represent all rates in each category\r
when the discrete-gamma model is implemented.  See Yang (1994\r
J. Mol. Evol. 39:306-314) for the distinction.  However, a bug caused\r
a few recent versions to use the median instead.  This affects baseml\r
versions 3.14c and 3.15, while version 3.14b or earlier and version 4\r
or later are correct.  codeml in version 4a & 4b were incorrect while\r
other versions are correct.  I hope both baseml and codeml are correct\r
from version 4c.\r
\r
\r
baseml and codeml (branch and clade labels): The local clock models\r
(clock = 2) in the two programs and also the branch or branch-site\r
models in codeml use labels to identify branches.  When nested clade\r
label ($) are used, the program may get confused so that the branches\r
are not labeled correctly.  I think this is not fixed in version 4c.\r
\r
\r
Version 4b, January 2008\r
\r
mcmctree.  Changed a variable name in the control file from MaxRootAge\r
to RootAge.  This should be specified in any of the following formats: \r
\r
       RootAge = <1.0  * constraint on root age, used if no fossil for root.\r
       RootAge = '<1.0'  * constraint on root age, used if no fossil for root.\r
       RootAge = <1.0>0.8  * constraint on root age, used if no fossil for root.\r
       RootAge = '<1.0>0.8'  * constraint on root age, used if no fossil for root.\r
\r
Added some text around the output, to make the results more comprehensible.\r
\r
evolver (option 8).  option 8 for calculating bipartition distances\r
between trees is broken.  I changed the program to prepare for a\r
figure and forgot to remove the changes afterwards.  This is now\r
fixed.\r
\r
\r
Version 4, July 2007\r
\r
mcmctree.  The MCMC algorithm of Rannala and Yang (2007) dealing with\r
rate drift is added.  The section of manual on mcmctree is now\r
rewritten.  An example data set is included in the\r
examples/DatingSoftBound/ folder.\r
\r
codeml. This now implements the mutation-selection models of codon\r
substitution of Yang and Nielsen (2008).  I should include more\r
details after the models are better tested, for example, after the\r
paper is written.\r
\r
baseml.  The log likelihood under the gamma-rates model (that is, when\r
alpha > 0 in the control file) is not calculated correctly, or not as\r
intended.  Version 3.14a is correct, but versions 3.14c and 3.15 are\r
incorrect.  I am not sure about version 3.14b, as I have not kept a\r
copy.  The bug causes the median to be used to represent all rates in\r
each category of the discre-gamma model whereas the mean was intended.\r
See the discussion around equation (10) and the paragraph below it in\r
Yang (1994. J. Mol. Evol. 39:306-314).  The bug is not expected to\r
have much effect on tree topologies or branch lengths.  Thanks to\r
Simon Whelan for reporting the bug.\r
\r
\r
\r
Version 3.15, March 2006\r
\r
evolverNSbranches (simulation program to generate data under the\r
branch models of codon substitution).  A bug was introduced after\r
version 3.14a so that versions 3.14b (?), 3.14c, and 3.15 have a bug\r
that makes the program abort prematurely.  This is now corrected.\r
Thanks to N. Clark for reporting the bug.\r
\r
\r
Version 3.15, January 2006\r
\r
evolver.  The simulation option of the evolver program in version\r
3.14b (?), 3.14c, and 3.15 has a serious bug that makes the program\r
generate nonsensible sequences when the user tree has 0 internal\r
branch lengths.  The bug was introduced between versions 3.14a (which\r
is correct) and 3.14c (which is wrong).  I don't seem to have kept\r
3.14b to check whether it is in error.  In this case, the smart idea\r
that led to the bug was the thought that no evolution occurs if the\r
branch length is 0 (or less than 1e-20, to be more precise).  If the\r
branch length is 0, there is no need to evolve the sequence along the\r
branch, so I decided to copy the sequence at the start of the branch\r
to the node at the end of the branch.  This would be fine, except that\r
I introduced a bug at the same time that caused the program to skip\r
generating sequences at all nodes that are descendents of the\r
zero-length branch.  As a result, sequences at tips that are\r
descendents of the zero-length branch have no data.  On some systems\r
the uninitialized space is initialized as 0, and then the sequences\r
will be printed out as consisting of all Ts (T is the first nucleotide\r
in my program).\r
\r
\r
At the same time, I also turned on the option of simulating data on\r
random trees in evolver.  You change fixtree=1 into 0, and recompile.\r
Then use MCbaseRandomTree.dat to simulate data instead of MCbase.dat:\r
\r
    evolver 5 MCbaseRandomTree.dat\r
\r
\r
\r
Version 3.15,  November 2005\r
\r
mcmctree: This implements the MCMC algorithm of Yang and Rannala (2005\r
MBE) for estimating species divergence times using soft bounds.  The\r
example file is in examples/SoftBound/.  Look at the readme file and\r
see whether you can duplicate our results in the paper.  Note that\r
this program works with DNA sequence data only, and assumes the\r
molecular clock (clock = 1).  Note that assumption of the clock can\r
lead to seriously biased results if the clock is violated.  Work is\r
under way to relax the clock assumption.  \r
\r
yn00: Added three methods for estimating dS and dN: LWL85 (Li, Wu &\r
Luo 1985), LPB93 (Li 1993; Pamilo and Bianchi 1993), and LWL85m (Yang\r
2006 Computaional Molecular Evolution, Oxford University Press.  Eqs.\r
2.12 & 2.13).\r
\r
codonml.  Added example files for clade model C, in folder\r
examples/CladeModelC.  Thanks for Joe Bielawski for preparing these.\r
\r
Added a sequence data format, in which the numbers or frequencies of\r
site patterns are inputed instead of the sequences.  The format is\r
indicated by the option variable P on the first line of the data file,\r
used in the same way as option variables I (for interleaved format)\r
and S (for sequential format, which is the default).  See the section\r
on sequence data format in the paml documentation for details.  The\r
old control variable readpattf for baseml and codeml is now deleted.\r
You can also use evolver to simulate data using this format.  Look at\r
the control files MCbase.dat, MCcodon.dat, and MCaa.dat.  This format\r
is useful for very large data sets with >1M sites, when collapsing\r
sites into patterns can take a long time.  \r
\r
\r
\r
Version 3.14b, May 2005\r
\r
pamp crashes when there are more than 127 or 255 (depending on the\r
system) species in the data.  The bug is due to my declaration of the\r
variable visit[] in the routine PathwayMP1 in pamp.c as char rather\r
than int.  Another problem is that the formula used to calculate\r
distances under REV (GTR) is sometimes inapplicable and the program\r
prints out messages like "err: DistanceREV".  I introduced some ad hoc\r
way of dealing with the problem.\r
\r
\r
\r
Version 3.14b, April 2005\r
\r
aaml (codeml with seqtype = 2): \r
\r
Under the REV model (model = 9), codeml failed to print out the\r
estimated amino acid substitution rate matrix.  This is fixed.\r
\r
\r
aaml (codeml with seqtype = 2 or 3):\r
\r
If you have multiple data sets (ndata > 1) and analysing protein data\r
sets under the empirical models such as wag, jtt, dayhoff.  The\r
results for the first data set are correct, but all later data sets\r
are analyzed incorrectly under the corresponding +F models, that is,\r
wag+F, jtt+F, dayhoff+F.  A bug in the program means that for the\r
second and later data sets, the equilibrium amino acid frequencies are\r
from the real data and not correctly set to those specified by the\r
empirical models.  Thanks to Ben Evans.\r
\r
\r
Version 3.14b, February 2005\r
\r
codonml and aaml:\r
\r
There is a bug that affects the joint reconstruction of ancestral\r
sequences when some branch lengths are equal to each other and when\r
those branches are visited one after the other in the program.  For\r
example, if some branch lengths are estimated to be 0, such a problem\r
might occur.  I think this bug is in the program for a long time,\r
although I note that joint ancestral reconstruction is turned off in\r
some versions for other reasons.  The marginal reconstruction is not\r
affected by this bug.  The update fixes this bug.\r
\r
\r
baseml and codeml:\r
\r
The program crashes on some systems for the local clock models (clock\r
= 5 and 6) after the calculation is finished.  If you manage to run\r
the model, the results should be fine.  There are bugs in the program\r
when freeing memory at the end of the calculation.  They are fixed in\r
this update.  See the discussion site at\r
\r
http://www.rannala.org/phpBB2/viewtopic.php?t=228\r
\r
Thanks to fjvera.\r
\r
\r
\r
Version 3.14a, December 2004\r
\r
codonml (codeml for codons):\r
\r
The mechanisctic amino acid substitution models (table 3 in Yang et\r
al. 1998 MBBE 15: 1600-1611) appear to be broken in paml/codeml\r
versions 3.13 and 3.14.  They were correct in version 3.12.  This is\r
now fixed.  If you use those models, please run the examples in the\r
folder examples/mtCDNA/ and compare results with those published in\r
the paper to confirm the program.  The results in the paper are\r
correct, but the program implementation may be broken due to lack of\r
maintenance.\r
\r
codonml (codeml for codons)\r
\r
When you fit the branch models with three or more branch types (omega\r
ratios), the program aborts with an error message saying that only two\r
omega ratios are allowed.  This is due to a bug in the program and is\r
now fixed.\r
\r
\r
Version 3.14, September 2004\r
\r
\r
1. codonml (codeml for codons): \r
\r
(1a) NSsites models M1 (neutral), M2 (selection), and M8 (beta& ) have\r
been changed.  The references are two manuscripts that are not yet\r
published (Wong, Yang, Goldman & Nielsen, 2004 Genetics 168:\r
1041-1051; Yang, Wong, and Nielsen 2005 MBE).  The details are as in the\r
following table.\r
\r
[table missing or never written.]\r
\r
The old models M1 and M2 are noted to be rather unrealistic as they do\r
not account for sites with 0 < w0 < 1.  Thus in the new models 0 < w0 < 1\r
is estimated from the data.  Also insisting on a class of sites with 0\r
= 1 appears to help avoid false positives as sites under weak\r
selection will be absorbed in this neutral class rather than being\r
claimed to be under positive selection.  Under M8, the constraint s >\r
1 is now automatically enforced to help avoid multiple local peaks and\r
also to make the model a positive-selection model.  Note that the\r
newer models replace the old ones, which are available in older\r
versions of paml/codeml only.\r
\r
After the changes, the models are nested as follows: M0 (one ratio) <\r
M1a (NearlyNeutral) < M2a (PositiveSelection) < M3 (discrete), and M7\r
(beta) < M8 (beta& ), where "<" means the model on the left is a\r
special case of the model on the right:.  Two LRTs can be constructed\r
using those models.  The first one compares M1a against M2a, using a\r
chi square with df 2, and the second one compares M7 against M8, using\r
a chi square with d.f. = 2.\r
\r
(1b) Similarly branch-site model A (Yang and Nielsen 2002) is changed.\r
The old model assumes 0 = 0 and is unrealistic.  This is replaced by 0\r
< 0 < 1, estimated from the data.  The new model is still called\r
branch-site model A.  It can be compared with the new M1a\r
(NearlyNeutral) to form a likelihood ratio test, with d.f.  2.  This\r
is called test 1.  Another test, called test 2, uses a null model A\r
but with 2 = 1 fixed (use fix_omega = 1 and omega = 1).  Test 2 is\r
very conservative, but test 1 might mistake relaxed selective\r
constraint on the foreground branch as positive selection.\r
\r
\r
(1c) Similarly the clade model C ( see also Forsberg and Christiansen\r
2003; Bielawski and Yang 2004) is changed as follows.  The new model C\r
replaces 0 = 0 by 0 < 0 < 1 and has five parameters in the\r
distribution: p0, p1, 0, 2, and 3.  The new model C can be compared\r
with the new M1a (NearlyNeutral), with d.f.  3.\r
\r
(1d) The empirical Bayes procedure for calculating posterior\r
probabilities of site classes under the site models, branch-site\r
models, and clade models are called naïve empirical Bayes (NEB).  They\r
are naïve because they use the MLEs of parameters without accounting\r
for sampling errors in them.  This problem is now fixed using a\r
procedure called Bayes empirical Bayes (BEB) (Yang, Wong, and Nielsen\r
2005).  The BEB procedure is implemented for the site models M2a and\r
M8, the (modified) branch-site model A, and the (modified) clade model\r
C.  The old NEB calculation is kept in the output and the new BEB\r
results follow the NEB result.  Perhaps I will remove the NEB\r
calculation in later versions.\r
\r
\r
baseml/codeml: rewrote likelihood clock and local clock models.\r
Implemented models for combined analysis of multiple genes\r
incorporating multiple calibration points (Yang and Yoder 2003).  The\r
ad hoc rate smoothing procedure of Yang (2004) is implemented for\r
nucleotide, amino acid, and codon substitution models, and the\r
implementation also deals with missing species at some loci.  The\r
option variable clock in the control files is now used differently\r
from before.  See later in this documentation and also the readme and\r
readme2 files in the folder examples/MouseLemurs/.\r
\r
baseml/codeml local clock models.  clock = 5 and 6 are added in\r
3.14beta5.\r
\r
codeml: added branch-site models C and D (Bielawski and Yang in press).\r
\r
mcmctree: this program is disabled in this release.  The old program\r
died and the new program is still under construction.\r
\r
The main body of this documentation may not be up to date. \r
\r
\r
\r
Bug fixes and minor corrections\r
\r
\r
baseml/codeml.  The output in the file rates under the gamma model\r
lists inferred rate for each site.  The output is incorrect if the\r
tree is large and scaling is used to avoid underflow.  The use of\r
scaling is indicated by two lines of output on the monitor like the\r
following:\r
\r
"2 node(s) used for scaling (Yang 2000 J Mol Evol 51:423-432):\r
 155 350".  \r
\r
The results are clearly wrong as the probabilities are much greater\r
than 1, and the rates are many orders of magnitude too large, etc.\r
Also under the same problematic condition the expected numbers of\r
sites with certain site patterns in the file lnf are incorrect.  When\r
scaling is not used, the results should be o.k. and look reasonable.\r
This affects versions 3.13 and 3.14beta1-3.  This is fixed in version\r
3.14.  Thanks for Nick Goldman for pointing out the error.\r
\r
yn00.  The program crashes for large datasets with many codons due to\r
a memory allocation error.  This affects versions 3.13 and some\r
versions of 3.14beta1.  The problem is fixed in version 3.14.\r
\r
aaml (codeml for amino acids).  Model REVaa_0 is broken due to lack of\r
maintenance in versions 3.1, 3.11, 3,12, 3,13, and 3,14beta1-4.  It\r
seems to be working in version 3.  This is fixed in 3.14beta5.  2\r
April 2004.\r
\r
codonml (codeml for codons) in branch-site models (model = 2 NSsites =\r
2 or 3) prior to 3.14beta3 have a scaling problem, which makes the\r
length of the foreground branch to be incorrectly estimated.  Other\r
branch lengths and substitution parameters, such as the parameters in\r
the distribution, are corrected calculated, as well as the log\r
likelihood values and the posterior probabilities.  27 March 2004\r
\r
pamp: The gamma parameter for variable rates among sites using the\r
method of Yang & Kumar (1996) is not estimated correctly in 3.14beta1,\r
beta2, beta3.  Versions 3.13d and earlier seem fine.\r
\r
codonml (codeml for codons) in version 3.14beta and 3.14beta1 does not\r
work when fix_omega = 1 is used for branch models.  The results are\r
incorrect.  Please use newer versions 3.14beta3 or later.  Also\r
version 3.13 is fine.\r
\r
baseml: Models TN93 and REV are wrong when used with Mgene = 3 or 4.\r
This seems correct in version 3.12 but went wrong in 3.13 and 3.14\r
since I inserted TN92 between HKY85 and TN93.  Thanks for Lee Bofkin\r
for reporting the error.\r
\r
codonml (codeml for codons): When codon models are used to reconstruct\r
ancestral sequences (with RateAncestor = 1), the program lists\r
synonymous and nonsynonymous changes at each site under the heading\r
"Changes at sites (syn nonsyn)."  This listing is incorrect due to a\r
bug in the program.  This bug affects versions 3.12, 3.13 and 3.14,\r
and version 3.11 seems correct.  Thanks to Joe Bielawski.\r
\r
baseml/codeml: The SE's for divergence times under the clock models\r
are calculated incorrectly.  This happens when you use clock = 1 or 2,\r
supply fossil date to calculate absolute times, and request standard\r
errors for times.  The estimates of times themselves are correct, but\r
standard errors for times are wrong.  The SEs for times and for the\r
rate under the TipDate model (clock = 3) are wrong as well.  The\r
programs print out the variances after the +-, instead of their square\r
roots.  This error was introduced in version 3.13.  Versions prior to\r
3.13 are correct.  I posted an update 3.13d to fix this bug.\r
\r
evolver: the simulation program can now accept species names in the\r
tree.  Note that the data file formats for evolver (MCbase.dat,\r
MCcodon.dat, MCaa.dat) have changed and you might have to change your\r
own data files.  Look at the files included in the package.\r
\r
The documentation in paml v3.13 from August 2002 - 12 December 2002\r
had a mistake about the critical values for the newer test using a\r
modified M8.  The critical values for the test are 2.71 at the 5%\r
significance level and 5.41 at the 1% level, rather than 1.95 and 3.32\r
as in the documentation.\r
\r
\r
\r
\r
\r
Version 3.14beta2, November 2003\r
\r
(1) baseml: Models TN93 and REV are wrong when used with Mgene = 3 or\r
    4.  This seems correct in version 3.12 but went wrong in 3.13 and\r
    3.14 since I inserted TN92 between HKY85 and TN93.  Thanks for Lee\r
    Bofkin for reporting the error.\r
\r
(2) codonml (codeml for codons): When codon models are used to\r
    reconstruct ancestral sequences (with RateAncestor = 1), the\r
    program lists synonymous and nonsynonymous changes at each site\r
    under the heading "Changes at sites (syn nonsyn)."  This listing\r
    is incorrect due to a bug in the program.  This bug affects\r
    versions 3.12, 3.13 and 3.14, and version 3.11 seems correct.\r
    Thanks to Joe Bielawski.\r
\r
\r
\r
Version 3.14beta, 20 July 2003\r
\r
(1) baseml/codeml: rewrote likelihood clock and local clock models.\r
   Implemented models for combined analysis of multiple genes using\r
   multiple calibration points (Yang and Yoder 2003).  The option\r
   variable clock in the control files is now used differently from\r
   before.  See documentation below and also the readme and example\r
   files in the folder paml/examples/MouseLemurs/.\r
\r
(2) codeml: added branch-site models C and D (Bielawski and Yang\r
   submitted).\r
\r
(3) baseml/codeml: The SE's for divergence times under the clock\r
   models are calculated incorrectly.  This happens when you use clock\r
   = 1 or 2, supply fossil date to calculate absolute times, and\r
   request standard errors for times.  The estimates of times\r
   themselves are correct, but standard errors for times are wrong.\r
   The SEs for times and for the rate under the TipDate model (clock =\r
   3) are wrong as well.  The programs print out the variances after\r
   the +-, instead of their square roots.  This error was introduced\r
   in version 3.13.  Versions prior to 3.13 are correct.  I posted an\r
   update 3.13d to fix this bug.\r
\r
(4) mcmctree: this program is decommissioned in this release.  The old\r
   program died and the new program is still under construction.\r
\r
(5) evolver: the simulation program can now accept species names in\r
   the tree.  Note that the data file formats for evolver (MCbase.dat,\r
   MCcodon.dat, MCaa.dat) have changed and you might have to change\r
   your own data files.  Look at the files included in the package.\r
\r
(6) Improved the iteration algorithm a little.  No change to the\r
   interface.\r
\r
(7) The documentation in paml v3.13 from August 2002 - 12 December\r
   2002 had a mistake about the critical values for the newer test\r
   using a modified M8.  The critical values for the test are 2.71 at\r
   the 5% significance level and 5.41 at the 1% level, rather than\r
   1.95 and 3.32 as in the documentation.\r
\r
(8) Edited the manual pamlDOC.pdf.\r
\r
\r
\r
\r
15 May 2003\r
paml3.13d\r
\r
baseml/codeml: The SE's for times under the clock models are\r
calculated incorrectly.  This happens when you use clock = 1 or 2,\r
supply fossil date to calculate absolute times, and request standard\r
errors for times.  The estimates of times themselves are correct, but\r
standard errors for times are wrong.  The SEs for times and for the\r
rate under the TipDate model (clock = 3) are wrong as well.  The\r
programs print out the variances after the +-, instead of their square\r
roots.  \r
\r
This error was introduced in version 3.13.  Versions prior to 3.13 are\r
correct.\r
\r
\r
12 December 2002 correction of documentation\r
\r
The paml documentation in paml v3.13 from August 2002 - 12 December 2002\r
has a mistake about the critical values for the newer test using a\r
modified M8.\r
\r
The critical values for the test are 2.71 at the 5% significance level\r
and 5.41 at the 1% level, rather than 1.95 and 3.32 as in the\r
documentation.  I made a mistake and quickly discovered it.  I thought\r
I corrected this a few months ago, but apparently the document\r
included in the paml release had this mistake.  I am sorry.  I'll\r
correct it now.\r
\r
\r
\r
Version 3.13, August 2002\r
\r
(1) baseml: nhomo = 1 and model = 6 (for both the one-rate and gamma-rates\r
    models) in version 3.12 gives incorrect results.  The error was introduced \r
    in version 3.12 and is now fixed.  Sorry for the trouble and thanks to\r
    Joshua Scott Rest.\r
\r
(2) baseml: Tamura's (1992) model is inserted between HKY85 and TN93.  \r
    Unfortunately this means that the model code for TN93 and REV are shifted, \r
    and you have to change your control file for those two models.  I also \r
    added user-specified REV or UNREST models.\r
    0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85, 5:T92, 6:TN93,\r
    7:REV, 8:UNREST, 9:REVu; 10:UNRESTu\r
\r
(4) I implemented the nonhomogeneous models of Galtier and Gouy\r
    (1998).  The options are model = 5 (T92) and nhomo = 3 (N1) or 4\r
    (N2) or 5 (user).  nhomo = 5 is added for the model of Yang and\r
    Roberts (1995) as well.  It means that the user specifies, in the\r
    tree, how many sets of base frequencies (or how many GC content\r
    parameters in the model of Galtier and Gouy 1998) should be used\r
    and which branches has which sets.\r
\r
(5) codeml: For pairwise amino acid sequence comparison (seqtype = 2\r
    and runmode = -2), the program crashes.  This is fixed now.\r
    Thanks to Dennis Paul Wall.  \r
\r
(6) codeml: The mechanistic models of amino acid substitution (model = 6 or 7)\r
    are broken for quite some time due to lack of maintenance.  I have\r
    now fixed them.  See the file examples/mtCDNA/readme.txt.\r
\r
(7) codemlsites is now removed, which is for batch run of multiple NSsites\r
    models.  This is now done by specifying several NSsites models in\r
    the control file codeml.ctl for codeml, for example: NSsites = 0 1\r
    7 8.\r
\r
(8) baseml & codeml used to ask the user what to do if the tree has branch \r
    lengths.  I have now introduced a variable fix_blength in the control \r
    file that does the same thing. \r
    fix_blength = -1  * 0: ignore, -1: random, 1: initial, 2: fixed\r
\r
\r
\r
Version 3.12, March 2002\r
\r
(1) baseml and codeml, Mgene=2 or 4 (different pi's) and RateAncestor\r
    = 1 Ancestral reconstruction does not work properly although the\r
    ML iteration is fine.  The output likelihood values (in the form\r
    "Node 20: lnL = -1047.965624") are different from the optimum\r
    likelihood already calculated.  The option was fine in version\r
    3.0e and 3.1, and was broken in version 3.11.\r
\r
(2) codeml stops with an error message "pi=0" when some amino acids are \r
    missing in the sequences.  Earlier versions of the program are correct.\r
\r
(3) yn00 now accepts multiple data sets (specified by ndata).\r
\r
\r
Version 3.11, September 2001\r
\r
(1) baseml and codeml (clock = 1 or clock = 2) sometimes print out\r
    ridiculous (e.g., large negative) tree lengths and branch lengths.\r
    I thought this was fixed in version 3.1, but it was not.\r
    Hopefully it is now fixed.\r
\r
(2) baseml and codeml: Malpha (a separate alpha for each site\r
    partition) does not work with Mgene = 0.  The output will give you\r
    the same alpha value for each gene as you specify in the control\r
    file.  This is now fixed.\r
\r
(3) baseml: When the REV or UNREST models (model = 6 or 7) are used\r
    together with Mgene = 2, 3, 4, the program outputs one Q matrix,\r
    when each site partition should have a different Q matrix.  The ML\r
    parameter estimates are correct, but I did not bother to format\r
    the output correctly.  The single Q matrix in the output is the\r
    one for the first partition.  I have added the output for the\r
    other matrices now.\r
\r
(4) baseml and codeml.  Ancestral reconstruction was disabled for\r
    large trees with more than 100 or 200 sequences in previous\r
    versions.  It is now allowed.  On large trees, the multiplication\r
    of transition probabilities across branches in the tree cause the\r
    product to become too small to be represented in the computer.\r
    This is known as underflows.  I implemented a scaling algorithm to\r
    avoid underflows some time ago but did not bother to do it with\r
    the ancestral reconstruction algorithm at the time.  The algorithm\r
    should now work, up to the upper limit set in the program, that\r
    is, 1000 or 2000 sequences.\r
\r
\r
\r
Version 3.1, June 2001\r
\r
Changed the way of counting base or amino acid frequencies so that\r
this version of baseml or codeml might generate different results from\r
previous versions under the following model specifications.  \r
\r
   . baseml & aaml: if you have multiple genes and the model uses the same\r
   frequencies in the model for the different site partitions.  \r
\r
   . codonml:\r
   if you have ambiguity characters at all.  Ambiguity characters were\r
   ignored when counting frequencies in previous versions but are used now.\r
\r
There is no difference if your data do not contain any alignment gaps\r
or ambiguity characters or if you choose to remove them before any\r
analysis (cleandata = 1).  The early versions used two ad hoc methods\r
in dealing with such ambiguity characters in different places of the\r
programs.  The new version uses another ad hoc method and in all\r
places.  I think it will be fine which ever version you use, but do\r
not use different versions in the same analysis.\r
\r
codonml (codeml with seqtype = 1) with Mgene = 1 does not work\r
correctly if the sequences have alignment gaps or ambiguity\r
characters.  Mgene = 1 means separate analysis.  The bug causes data\r
for the genes except the first to be read out of frames, so that most\r
likely you will see a lot of messages "stop codon XXX" before the\r
program aborts.  This is now fixed.\r
\r
\r
codonml (codeml with seqtype = 1) with fix_alpha=0: In the final\r
output for the codon + gamma model, parameter alpha is not correctly\r
identified.  If you have fix_kappa = 0, fix_omega =0, and fix_alpha=0,\r
the last three numbers should be kappa, omega (dn/ds), and alpha.  The\r
detailed output did not identify alpha correctly, although the\r
likelihood value and the outputted rates and frequencies (r and f) are\r
correct.  I believe this is now corrected.\r
\r
\r
\r
Version 3.0e, 14 May 2001\r
\r
codonml (codeml with seqtype = 1) with model = 1: The program crashes\r
when you have a small tree with <= 5 branches.  This is a new error\r
introduced in version 3.0c.  It is now fixed.\r
\r
baseml and codeml (clock = 1 or clock = 2) sometimes print out ridiculous\r
tree lengths and branch lengths.  \r
\r
yn00 in version 3.0c ignores icode and uses the universal code only.\r
It was correct in vesion 3.0b.  This is fixed in version 3.0d.\r
\r
\r
Versions 3.0c-3.0d, 20 February 2001\r
\r
We bought a Pentium III cluster running redhat linux 6.2.  The gcc\r
compiler 2.61 seems to have bugs and codeml does not run for NSsites =\r
7 or 8.  The program might output "points out of order" or "Det goes\r
to 0" error messages.  3.0d seems to be safer for these models.\r
\r
codonml (codeml with seqtype = 1) with NSsites = 3 might generate a\r
message "error: sortwM3." and then stop with no output after the\r
iteration.  This is now fixed.  During the iteration, there was no\r
restriction on the w values, so that w0 can be greater than w1 and so\r
on.  At the end of the iteration, I sort the three w values so that w0\r
< w1 < w2 before outputting.  3.0c has a bug at sorting.  If you don't\r
see the error message, everything is fine.\r
\r
\r
10 November 2000\r
\r
The probabilities for Joint ancestral reconstruction are not correct.\r
Many values are larger than 1.  Fixed in version 3.0c.\r
\r
\r
\r
Version 3.0c, 30 October 2000\r
\r
(1) Changed parameterization of clock models (clock = 1, 2, or 3), so that\r
the output now lists node ages, and the SEs are for node ages as well.\r
Node ages are measured by the expected number of nucleotide or amino\r
acid substitutions per site (nucleotide, amino acid, or codon) from\r
the node to the present time, and are proportional to the divergence\r
times.  In earlier versions, the output list proportions.\r
\r
(2) codonml: problem with NSsites model with multiple gene data\r
(Mgene=1) is now fixed(?).  In the earlier version, the likelihood\r
calculation was o.k., but identification of sites under positive\r
selection was incorrect as the program attempted to list sites in\r
another gene when it was analyzing data for one gene.  (1) Mgene\r
options for codons and amino acids:\r
\r
\r
\r
Version 3.0b, 6 October 2000\r
\r
The local clock models (clock = 2) are not implemented correctly in paml3.0 or 3.0a. \r
\r
Yoder, A. D., and Z. Yang.  2000.  Estimation of primate speciation\r
dates using local molecular clocks.  Molecular Biology and Evolution\r
17 (7): 1081-1090.\r
\r
I introduced a serious error before I made the models available to the\r
public, in ver3.0 (May 2000).  Two symptoms are noted right now: (1)\r
crashes (which are good as they make you realise that something is\r
wrong immediately) and (2) the local clock model (clock=2) having a\r
worse likelihood than the global clock model (clock = 1).  The results\r
in the paper are correct except for a typo ("37.30 - 48.48" in table 3\r
for the cetacean-whale divergence using the C3 calibration should be\r
"57.30 - 48.48").  I now include the example data sets and results\r
from that paper to help you prepare your own data files.\r
\r
Sorry for any damages done.\r
\r
Many thanks to Toby Johnson and Philip Awadalla for spotting the error.\r
\r
\r
\r
Version 3, 4 February 2000\r
\r
The program evolver may not generate sequences correctly if branch\r
lengths are very large, say >10.  \r
\r
Scaling by nodes to avoid underflow.  When there are many sequences in\r
the data (say > 200 or 300), the probability of observing data at a\r
site can become too small to represent in the machine.  This is\r
because the probabilities are multiplied along branches in the tree\r
and a large tree has many branches.  Since version 2.0, this underflow\r
is dealt with by scaling, and it works with both the one-rate and\r
gamma-rates models.  It did not work with the NSsites models in\r
codonml, but probably nobody has analysed large data sets to notice\r
this.  Anyway, this version hopefully fixes that problem, when I was\r
analysing the data set of Bush et al. (1999 MBE 16:1457-1465), which\r
has about 350 sequences.  If you use the program on very large data\r
sets, watch out for anything strange and let me know.  I suppose with\r
the option of using supplied branch lengths, data sets of such sizes\r
will become manageable, and this scaling will become important.\r
\r
\r
\r
Version 2.0h, Oct 14, 1999\r
\r
Fixed the following errors, all simple errors.\r
\r
The RemoveIndel routine in 2.0e and 20.f has a bug that removes most\r
sites in the sequence except those with all T's and C's for baseml.\r
The same routine in 2.0g has a different bug that removes most sites\r
except those will all T's and C's in the sequence for aaml.  In either\r
case, the results will be grossly wrong.  The bugs affect the analysis\r
only if you have gaps or ambiguity characters in the data and use\r
runmode = 1 for aaml in 2.0e and 2.0f, and runmode = 1 for codonml for\r
2.0g.\r
\r
The newly introduced option for amino acid reconstruction by baseml\r
(icode) does not work.  The program aborts with an error message.\r
This is now fixed.\r
\r
\r
\r
\r
Version 2.0g, Sept 21, 1999\r
\r
Fixed a bug in 2.0f that causes crash in codeml when getSE = 1.  \r
\r
If the tree has branch lengths, codeml and baseml in 2.0f will simply\r
use the branch lengths rather than estimating them by iteration.  I\r
have added a question for the user to choose to ignore the branch\r
lengths, use them as fixed or use them as initial values for the\r
iteration.\r
\r
\r
\r
Version 2.0d, June 30, 1999\r
Fixed some convergence problems in ML pairwise estimation of dS and \r
dN, in particular, in presence of data partitions (the G option).\r
\r
Fixed a bug in evolver.  The bug means that no sensible data are\r
generated if the tree has branch lengths smaller than 1e-7.  At some\r
point, I did something "clever" (taking advantage of the fact that at\r
such small branch lengths, no evolution would occur), which destroyed\r
the calculation.  \r
\r
\r
\r
Version 2.0c, May 12, 1999\r
Corrected a few problems with the evolver program, and with running \r
multiple data sets using baseml.  Added REV in MCbase.dat.\r
Thanks to John Mercer.\r
\r
Note--Uncomment ndata in baseml or codeml to analyse multiple\r
data sets generated from evolver.\r
\r
\r
\r
April 16, 1999\r
\r
Although posted only a few days ago, paml2.0 has an error in the\r
program codeml, which invalidates most analyses based on codon\r
substitution models.  The symptom is relatively easy to spot and is\r
indicated by one of the parameters (say omega) not being estimated at\r
all and in the output you will have exactly the same number as you\r
specified in the control file.\r
\r
Amino acid-based analysis or pairwise comparison of codon sequences\r
are not affected.  Use paml2.0a to redo the analysis.  Apologies for\r
the trouble.\r
\r
\r
\r
version 2.0, 31 March 1999\r
Added back the faster calculation when there is no missing data.\r
Enabled the programs/analysis that broke in the temporary versions, \r
such as pamp, joint reconstruction.\r
\r
Renamed listtree as evolver, and added simulation options.\r
\r
\r
\r
PAML 1.4 Temporary versions\r
Dec98 & Jan99 \r
\r
  pamlTMP.Jan99.notes\r
  ===================\r
  Ziheng Yang\r
  Fri Jan 15 10:01:42 GMT 1999\r
\r
  .Type make to compile the programs.  If you got error messages, see\r
  whether you know how to change the file Makefile.  If it does not work\r
  either, look at the files paml.cc, paml.gcc, etc.\r
\r
  .At the suggestion of David Posada, I added ML estimation of pairwise \r
  distances between amino acid sequences in the program codeml.  \r
\r
  The relavant variables in the control file codeml.ctl are \r
\r
   runmode = -2 for pairwise distance calculation\r
     aaRatefile = jones.dat * only used for aa seqs with model=empirical(_F)\r
              * dayhoff.dat, jones.dat, mtmam.dat, or your own\r
     model = 3 \r
     alpha = 0\r
\r
  model specifies the amino acid substitution model.  You should choose\r
  a value from 0, 1, 2, 3.  If you want to estimate the rate matrix from\r
  your data, you should do that using many sequences and then move the\r
  estimated substitution rate matrix (in the output file mlc) into a\r
  file and change the variable aaRatefile.  \r
\r
  Choosing alpha>0 means using gamma distance with that specified alpha\r
  parameter.  If you see anything strange, please let me know.\r
\r
  .I have included only three programs baseml, codeml, and listtree.\r
  I have not tested the others and so they are not included.\r
\r
\r
\r
  pamlTMP.Dec98.notes\r
  ===================\r
  Ziheng Yang\r
  Wed Dec  9 10:44:37 GMT 1998\r
\r
\r
  This temporary version of paml has involved a number of changes, and\r
  is quite likely to contain errors.  Please let me know if you notice\r
  anything strange.  Some of the options do not work for now.  Details\r
  follow.  I will compile the Win32 and PowerMac versions after\r
  everything is fixed and tested.\r
\r
  . Missing data: baseml and codeml now handle missing data.  Choose\r
  DelMissing = 1 in the control file to remove sites with ambiguity\r
  characters or alignment gaps, as did previous versions.\r
\r
  . Clock: The clock is now fixed.  It may still crash, but the option\r
  is not worse than the no-clock options.  The definitions of the node\r
  times or branch lengths are now different from the documentation\r
  (pamlDOC.html), and so check the branch lengths in the output tree\r
  topology.  Thanks to Jeff Thorne.\r
\r
  . Nexus file formats: Paml programs now read sequence files and tree\r
  structure files in the Nexus format used by paup and MacClade.  Only\r
  the data or tree are read, and everything else in the file is ignored.\r
\r
  .  Ancestral reconstruction: Marginal ancestral reconstructions are\r
  now done under the gamma-rates model as well as the constant-rate\r
  model.  For protein-coding genes, reconstructions can now be done at\r
  the nucleotide level with baseml (in particular, by using models that\r
  account for differences at the three codon positions), at the amino\r
  acid level with aaml, and at the codon level with aaml.  According to\r
  Belinda Chang, those different reconstructions tend to be similar.\r
  Results for ancestral reconstruction are now listed by site, and not\r
  by site pattern as in earlier versions.  Missing data are now handled\r
  in those analysis.  See instructions in this doc about baseml.ctl for\r
  more options.  Thanks to Belinda Chang.  I would like to take this\r
  opportunity to point out that the reconstructed ancestral sequences\r
  are not the real observed sequences and may contain errors.\r
\r
  . I have disabled the option for variable rates among sites in\r
  combination with variable dN/dS rate ratios among lineages in codonml.\r
  This option has never been implemented in the program but was not\r
  disabled.\r
\r
  . Distance matrices from codonml.  The program codonml outputs\r
  matrices of synonymous and nonsynonymous rates in all pairwise\r
  comparisons using the method of Nei & Gojobori (1986) into two files\r
  named DistanceNG.dN & DistanceNG.dN.  The matrices are lower-diagonal\r
  and can be used with programs, such as neighbor and fitch, in phylip\r
  (and possibly paup* too) for tree reconstruction or branch length\r
  estimation by distance methods.  If you choose runmode = -2 for ML\r
  pairwise comparison, the program will also output two matrices of ML\r
  estimates into two files named DistanceML.dS & DistanceML.dN, for\r
  synonymous and nonsynonymous rates.  Since many users seem interested\r
  in looking at dN/dS rate ratios among lineages, examination of the\r
  tree shapes indicated by branch lengths calculated from the two rates\r
  may be interesting although the analysis is intuitive.  Obviously, you\r
  should NOT name your own data files with those names; otherwise they\r
  will get overwritten.  If your species name has more than 10\r
  characters, you will have to edit the files and cut the names short\r
  before you can use Phylip programs.  I have decided to leave this work\r
  to the user.  I plan to create some other distance matrix files for\r
  nucleotides and amino acids as well.\r
\r
  . An minor prompt error when running codonml with model = 2 (variable\r
  dN/dS ratios among lineages) is fixed.  The example data set\r
  lysozymeSmall.nuc is included for the user to duplicate my results in\r
  Yang (1998 Molecular Biology and Evolution 15: 568-573).  The control\r
  file codemlLysozyme.ctl explains in more detail how to specify the\r
  options.  I would like to remind you that it is not correct to use the\r
  option model = 1 to estimate dN/dS ratios for branches to find out\r
  which ratios are greater than one, and then to use model = 2 to test\r
  whether that ratio is significantly greater than one.  The problem\r
  with such an analysis is that the hypothesis being tested is derived\r
  from the data which you use to test the hypothesis, and as a result,\r
  you tend to get significantly results too often.  Strictly speaking,\r
  the hypothesis should be formulated before the data are analysed.\r
\r
  Things that do not work in this version include\r
\r
  . All parsimony-based analyses are now broken.  The program pamp is\r
  not included.  Programs baseml and codeml used to calculate the MP\r
  score for each tree before doing an ML analysis on that tree.  Now a\r
  -1 is used to indicate this is not possible.\r
\r
\r
\r
Version 1.4 (UCL)\r
July 1998 \r
\r
     1. Most of the changes are in the program codeml (aaml for amino\r
  acid sequences and codonml for codon sequences). A few new models of\r
  codon substitution are implemented (see Yang and Nielsen 1998; Nielsen\r
  and Yang 1998; Yang 1998).\r
     2. Some problems relating to ancestral sequence reconstruction (in\r
  baseml and codeml) are fixed. The earlier algorithm does not work\r
  properly if the data contain more than several sequences. The\r
  algorithm makes a guess at the likely character states at the interior\r
  nodes of the tree, and then uses those to generate reconstructions\r
  (joint reconstructions, see eq. 2 in Yang et al. 1995) to be\r
  evaluated. Sometimes this strategy misses important reconstructions,\r
  and as a result, the probabilities for reconstructed characters at\r
  specific interior nodes (marginal reconstructions; see eq. 4 in Yang\r
  et al. 1995) are substantially underestimated if the number of\r
  sequences is not very small. I have now written separate codes to do\r
  the marginal reconstruction, evaluating all possible character states\r
  for each interior node. This works also under the gamma model of\r
  substitution rates among sites, a feature that was not implemented\r
  before. The joint reconstruction works with one rate for all sites\r
  only and has the old problem of possibly missing important\r
  reconstructions. However, the posterior probabilities for those joint\r
  reconstructions that do get evaluated are accurate. I have added an\r
  option (choose RateAncestor = 2 rather than 1) for the user to specify\r
  the reconstruction to be evaluated. The two approaches are expected to\r
  produce very similar results, and since the marginal reconstruction is\r
  always reliable, perhaps it can always be used for data analysis.\r
  (Thanks to Belinda Chang.)\r
     3. The output file for estimated rates for sites under the\r
  variable-rates models is now "Rates.out" instead of "rst".\r
     4. I have implemented ancestral sequence reconstructions based on\r
  codon-substitution models (codonml). This has not been tested\r
  carefully, and I would appreciate your comments if you use it.\r
     5. A number of minor changes have been made. I have fixed several\r
  floating exception errors that seem to occur on DEC Alpha machines\r
  only.\r
     6. The documentation is now changed into the html format. \r
\r
\r
\r
Version 1.3a,b,c (UCL)\r
June 1998 \r
\r
  1. Basemlg was put back at the request of a user.\r
\r
  2. Some changes are made to codeml concerning codon-substitution\r
  models (Yang and Nielsen 1998 JME; Yang 1998 MBE; Nielsen and Yang\r
  1998 Genetics)\r
\r
\r
\r
PAML Version 1.3: (UC Berkeley)\r
July 1997 \r
\r
  1. Starting with this version, the program basemlg will not be\r
  included in the package.  This program implements the (continuous)\r
  gamma model of substitution rate variation among nucleotide sites\r
  (Yang 1993).  It involves intensive computation and has been\r
  superseded by the discrete-gamma model implemented in the program\r
  baseml.  If you want to use this program, you should use previous\r
  versions of paml, which can be obtained from me if nowhere else.\r
\r
  2. The simulation program mcml.c is removed from the package.  This is\r
  now superced by the program listtree.c, which does different things\r
  besides simulating data sets.\r
\r
  3. I have implemented the stepwise addition algorithm in the programs\r
  baseml and codeml.  This algorithm is faster than the\r
  star-decomposition algorithm implemented before, but my implementation\r
  is crude and inefficient.  The NNI perturbation is implemented too.\r
  These work without the molecular clock.\r
\r
  4. Many changes are made with the program codonml (codeml with\r
  seqtype=1).  A few new codon-based models are implemented.  One\r
  assumes different nonsynonymous/synonymous (dN/dS) rate ratios among\r
  branches in the phylogeny, and may be useful for detecting episodic\r
  positive selection.  Another allows the dN/dS ratio to vary among\r
  sites and may be useful for testing neutrality and detecting\r
  positively-selected sites in the sequence.  These are based on my\r
  collaborative work with Rasmus Nielsen.  I made some changes to the\r
  output of estimates of dS and dN in pairwise comparisons (codonml with\r
  runmode=-2), so that the results are easy to understand.  Another\r
  change with this program is in the codon-substitution model of Goldman\r
  and Yang (1994).  We used the amino acid distance matrix of Grantham\r
  (1974) to modify the rate of nonsynonymous substitutions, but\r
  unfortunately the model does not fit data well.  It does not always\r
  fit the data better than if we simply assume equal distance between\r
  any two amino acids.  So I have added the option of ignoring the amino\r
  acid differences.  Goldman and Yang's original formula is still\r
  available but not recommended for use.  As a result of these changes,\r
  the formulation of the Goldman and Yang model is now somewhat\r
  different from the paper and previous versions of the program.  See\r
  Section 6.2 for details.\r
\r
  5. I have again made some changes with the user interface, based on my\r
  guess of what the "user" would like to see.\r
\r
\r
\r
Version 1.1a-d: (University of California, Berkeley)\r
1995-1996 \r
\r
1.1a No records of changes remain.\r
\r
1.1b November, 1995:\r
\r
The following line at the beginning of the file tools.h \r
      #define __cplusplus\r
is removed.  I added this line to avoid some (unjustified) warning\r
messages from a compiler I used to use and forgot to remove it when I\r
was preparing paml1.1.  This line made some compilers unable to\r
compile.\r
\r
A bug that caused baseml and basemlg to give wrong estimates of kappa\r
under the F84 model is now fixed.  Estimates of kappa under the F84\r
models listed in the paml 1.0 and 1.1 documents are all incorrect and\r
are a bit too small.  For example, the estimate of kappa under the F84\r
model for the mtDNA data of Brown et al. (1982) should be 4.344 while\r
the wrong estimate given in the large table of the document is 3.420.\r
Likelihood value, branch lengths and other parameters seem to have\r
been correctly calculated by paml 1.0 and 1.1.  Versions prior to 1.0\r
do not have this bug, and results in Yang (J. Mol. Evol. 1994,\r
39:306-314) and Yang, Lauder, and Lin (J. Mol. Evol. 1995,\r
41:587-596), where the F84 model was used, are correct.\r
\r
1.1c March 1996 \r
\r
Representation of tree topologies using branches (see the document)\r
are now allowed (restored) in the tree structure file.  This is for\r
coping with cases where some taxa are ancestors of others.  If this\r
representation is used, the line starts with a dollar sign.  So the\r
following line in the file trees.5s (the tree file for 5 species)\r
\r
$ 5      6 3  6 4  6 5  3 1  3 2\r
\r
represents the tree ((12)34), with 5 to be the ancestor of 1 and 2.\r
\r
Alignment gaps are now allowed when option G is used (for multiple\r
genes).  The gaps are removed before analysis.\r
\r
A small program listtree lists all the bifurcating trees for a given\r
(small) number of species.\r
\r
codonml (codeml.c with seqtype=1) now has an option for calculating\r
the number of synonymous substitutions per synonymous site (Ks) and\r
the number of nonsynonymous substitutions per nonsynonymous site (Ka)\r
and their ratio for each pairwise comparison, using the method of\r
Goldman and Yang (1994).  This is implemented with\r
runmode=-2. Estimates of Ks and Ka are collected in the file rst and\r
estimates of model parameters in the file mlc.  This option produces\r
table 2 in the Goldman & Yang paper.\r
\r
The interior nodes were incorrectly identified when the ancestral\r
sequences were listed in the file rst.  They were all one less than\r
their correct numbers.  Nodes 7, 8, 9, 10 for the stewart.aa data were\r
incorrectly identified as nodes 6, 7, 8, 9.  This is corrected now.\r
\r
1.1d. July 1996\r
\r
An unintended version of codeml.c was included in version 1.1c, which\r
causes the program to produce different results under the codon-based\r
model of Goldman and Yang (1994).  The small number in baseml.c was\r
set too small, making the program do unnecessary iterations; this is\r
fixed now.  I also changed the program pamp so that it will estimate\r
the substitution rate matrix without using a tree.\r
\r
I have included a few other genetic code tables in this version of\r
codeml; see the file codeml.ctl.\r
\r
\r
PAML Version 1.2: November 1996 (UC Berkeley)\r
\r
  1. I have done some changes with the user interface.  Name of files\r
(sequence data file, control file, tree structure file) are now\r
specified in the control files (The default are baseml.ctl,\r
codeml.ctl, pamp.ctl, mcmctree.ctl).  The command line for executing a\r
program is now ProgramName ControlFileName with ControlFileName\r
optional.  I have added some text in the output file of baseml so that\r
the results are eassier to understand.  Choose verbose=0 in the\r
control file if you do not want the detailed output.\r
\r
  2. A new program mcmctree is included, which performs Bayesian\r
estimation of phylogenies using Markov chain Monte Carlo methods\r
(Rannala and Yang 1996; Yang and Rannala 1997).\r
\r
  3. A number of minor changes are made, such as using species names\r
in the parenthesis tree representation.\r
\r
\r
\r
July 1995 \r
Version 1.1: (IMEG PennState)\r
\r
  Models for analyzing data of multiple genes (Yang 1996 JME) are\r
  implemented in baseml and codeml.  Two more programs are included:\r
  pamp for parsimony-based analysis implementing the methods of Yang and\r
  Kumar (1996 MBE) and mcml for Monte Carlo simulation (prepared for\r
  Arndt von Haeseler but apparently not used).\r
\r
\r
February 1995 \r
Version 1.0 (IMEG, PennState)\r
\r
  This is the first official version of paml.  Three main programs are\r
  included: baseml, basemlg and codeml.  Control files are used for all\r
  of them.  The program codeml is formed by merging two programs codonml\r
  and aaml, for codon sequences and amino acid sequences, respectively.\r
  One rate, gamma rates, and auto-discrete-gamma rates for sites are\r
  implemented for nucleotide, amino acid, and codon-based models.  The\r
  document is provided as postscript files (paml1_3.ps).\r
\r
\r
October 1994 \r
Version 0.12, no name for package (BAU and IMEG PennState) \r
\r
  The bracket notation of tree topologies is introduced.  Tree\r
  search by star decomposition is implemented in baseml and codonml.\r
  The options for baseml are moved into an option file baseml.ctl.  The\r
  molecular clock seems to work with the automatic tree search model\r
  (runmode=1 or 2).  The old readme file is renamed manual.\r
\r
\r
March 1994 \r
Version 0.10, no name for package: (NHM, England)\r
\r
  Programs tripml and comptree are renamed codonml and rell.  The readme\r
  file is expanded and the numerical optimization routine, ming1, used\r
  by dnaml and codonml, is improved.  The discrete-gamma model of Yang\r
  (1994 JME) is implemented in both baseml and codonml.\r
\r
  Preliminary version 0.11, no name for package: April 1994 (NHM,\r
  England) Programs dnaml and dnamls are renamed baseml and basemlg to\r
  avoid confusion with Joe Felsenstein's programs.  The\r
  auto-discrete-gamma model and the nonparametric models of rate\r
  variation and correlation mong sites of Yang (1995 Genetics) are added\r
  in baseml.  Nonhomogeneous-process models of Yang and Roberts (1995\r
  MBE) are implemented in baseml.\r
\r
\r
April 1993 \r
Version 0.01, no name for package (Cambridge).  Three programs are included in \r
the package:\r
\r
  1. dnamls implements the model of gamma rates among sites of Yang\r
  (1993) and works with arbitrary topologies and substitution models\r
  (JC69, K80, F81, HKY85).\r
\r
  2. dnaml implements the model of a single rate for all sites. (The\r
  same name as Joe Felsenstein's program was used.)\r
\r
  3. tripml implements the codon-based model of Goldman and Yang (1994).\r
  Trees are represented by the "branch notation".\r
\r
// end of file\r