examples/TipDate.HIV2/README.txt

   1 README.txt\r
   2 18 September 2011, Ziheng Yang\r
   3 \r
   4 This example data file is for duplicating the analysis of Stadler and\r
   5 Yang (2012).  The method dates divergence times when sequences have\r
   6 sample dates, as in the case of viral sequences.  The data set was\r
   7 aligned and analyzed by Lemey et al. (2003).\r
   8 \r
   9 \r
  10 (A) Data format \r
  11 \r
  12 Have a look at the sequence alignment file HIV2ge.txt.  Note that the\r
  13 end of the sequence name has the sample date.  For example, P03h1995\r
  14 was isolated in 1995.\r
  15 \r
  16 \r
  17 (B) ML analysis (Stadler and Yang 2012 table #)\r
  18 \r
  19 To run the ML analysis, cd to the current folder and then run baseml.\r
  20 This uses the control file baseml.ctl.\r
  21 \r
  22 On Windows:\r
  23 \r
  24    cd paml4.6\r
  25    cd examples\TipDate\r
  26    ..\..\bin\baseml\r
  27 \r
  28 On Unix/Linux/Mac OSX:\r
  29 \r
  30    cd paml4.6\r
  31    cd examples/TipDate\r
  32    ../../bin/baseml\r
  33 \r
  34 The TipDate model is specified by the following line in baseml.ctl,\r
  35 where the second number is the time unit (100 years in this case):\r
  36 \r
  37        TipDate = 1 100  * TipDate (1) & time unit\r
  38 \r
  39 \r
  40 (C) Bayesian analysis (Stadler and Yang 2012 table #)\r
  41 \r
  42 The control file is mcmctree.ctl is for running mcmctree.  Right now\r
  43 the seqfile and treefile are specified as follows.  These are set up\r
  44 so that you can run two copies of mcmctree in two different folders\r
  45 (r1 and r2 inside the folder) at the same time.\r
  46 \r
  47         seqfile = ../HIV2ge.txt\r
  48        treefile = ../HIV2ge.tre\r
  49 \r
  50 Start two command terminals, and cd to r1 and r2 respectively.  Then in each window type\r
  51 \r
  52    ..\..\..\bin\mcmctree ..\mcmctree.ctl (on Windows)\r
  53 \r
  54    ..\..\..\bin\mcmctree ..\mcmctree.ctl (on Unix/Linux/MacOSX)\r
  55 \r
  56 You confirm that the two runs produce very similar results.  If not,\r
  57 you need increase nsample or burnin etc.\r
  58 \r
  59 Have a look at the control file mcmctree.ctl, the paml/mcmctree\r
  60 documentation in the doc/ folder.\r
  61 \r
  62 The following line in the control file specifies the TipDate model.\r
  63 The second number is the time unit.\r
  64 \r
  65        TipDate = 1 100  * TipDate (1) & time unit\r
  66 \r
  67 When you run the default analysis, you will see the following printout on the monitor.\r
  68 \r
  69 TipDate model\r
  70 Date range: (1995.00, 1982.00) => (0, 0.13). TimeUnit = 100.00.\r
  71 \r
  72 The program scans the sequence names, and take the last field in the\r
  73 sequence name as the sampling date.  The most recent sample date will\r
  74 be time zero, and the other times will be rescaled using the time\r
  75 unit.  For example the sequences in the example dataset are from 1982\r
  76 to 1995, with 1995 becoming time 0, and 1982 becoming 0.13 as one time\r
  77 unit is specified to be 100 years.\r
  78 \r
  79 You should specify a prior on the age of the root of the tree.\r
  80 \r
  81        RootAge = B(0.5, 2.0, 0.01, 0.02)\r
  82 \r
  83 The above means that the age of root is between 0.5 and 2.0 time units\r
  84 (that is, between 1945 and 1795), but those bounds are soft, so that\r
  85 the minimum- and maximum-age bounds are violated with probabilities 1%\r
  86 and 2%, respectively.\r
  87 \r
  88 Other things that are important include the prior on the substitution rate \r
  89 \r
  90    rgene_gamma = 2 10   * gamma prior G(alpha, beta) for rate for genes\r
  91 \r
  92 Note that the gamma distribution G(alpha, beta) has shape parameter\r
  93 alpha and scale parameter beta.  You specify alpha depending on how\r
  94 much confidence you have (with alpha = 1 or 2 to be diffuse priors and\r
  95 alpha = 5 or 10 to be informative priors, say), and then specify the\r
  96 scale parameter by having the mean alpha/beta in the right range for\r
  97 the data.  Here alpha = 2 means a fairly diffuse prior, while\r
  98 alpha/beta = 0.5 means 0.5 changes per site per time unit, or 0.005\r
  99 changes per site per year (since the time unit is 100 years).\r
 100 \r
 101 The clock model may also be important.\r
 102 \r
 103       clock = 1    * 1: global clock; 2: independent rates; 3: correlated rates\r
 104 \r
 105 References\r
 106 \r
 107 |P, Pybus OG, Wang B, Saksena NK, Salemi M, Vandamme AM. 2003. Tracing the origin and history of the HIV-2 epidemic. Proc Natl Acad Sci USA 100:6588-6592.\r
 108 \r
 109 Stadler, T. and Z. Yang.  2012  Dating phylogenies with sequentially sampled tips.  Syst Biol, submitted\r
 110 \r
 111 \r
 112 //end of file\r