mcmctree - issue with 95% HPD predicted

[Annie Lebreton]

Dear PAML community,

I've been calibrating a Raxml tree on a time scale with mcmctree, but strangely I don't have any range for the 95% HPD predicted (e.g.  [&95%={2.172, 2.172}] ). When run multiple times, there is quite a difference between the resulting FigTree; acceptable one in my opinion; but I expect having some range in the 95% HPD. I didn't spot anything special in the logs. I tried sampling my dataset (full dataset of 931 genes, tests with 100 genes and 5 genes), increasing the nsample (20000 and 100000), using only 1 or 2 calibration points out of my 3... I'm running out of ideas...  

Any hint on why I have this pattern ?

Here is a tree I obtained:

((Mutel1: 4.309702, Sclhys1: 4.309702) [&95%={4.310, 4.310}]: 0.229867, (Albpec1: 3.664656, (((Hetan2: 0.915896, Stehi1: 0.915896) [&95%={0.916, 0.916}]: 2.050906, (Hercor1: 2.171631, Laxbic1: 2.171631) [&95%={2.172, 2.172}]: 0.795172) [&95%={2.967, 2.967}]: 0.364899, (((Varmin1: 0.766213, Ricme1: 0.766213) [&95%={0.766, 0.766}]: 0.197918, (Amycha1: 0.771304, Echtin1: 0.771304) [&95%={0.771, 0.771}]: 0.192827) [&95%={0.964, 0.964}]: 0.098220, ((Clapy1: 0.643540, (Lenvul1: 0.275602, Aurvu1: 0.275602) [&95%={0.276, 0.276}]: 0.367938) [&95%={0.644, 0.644}]: 0.130248, (Glocon1: 0.716106, ((Lacqui1: 0.310237, (Lacpsa1: 0.257866, (Lacviv1: 0.195272, ((Lacaka1: 0.125330, (Lacpse1: 0.122350, Lachen1: 0.122350) [&95%={0.122, 0.122}]: 0.002981) [&95%={0.125, 0.125}]: 0.041457, ((Lachat1: 0.096229, Lacdel1: 0.096229) [&95%={0.096, 0.096}]: 0.024483, Lacsan1: 0.120712) [&95%={0.121, 0.121}]: 0.046076) [&95%={0.167, 0.167}]: 0.028484) [&95%={0.195, 0.195}]: 0.062594) [&95%={0.258, 0.258}]: 0.052371) [&95%={0.310, 0.310}]: 0.149408, (Muloch1: 0.393179, (((Rusbre1: 0.123679, Ruscom1: 0.123679) [&95%={0.124, 0.124}]: 0.046728, (((Rusvin1: 0.098713, Rusoch1: 0.098713) [&95%={0.099, 0.099}]: 0.000683, (Ruseme1: 0.091769, Rusrug1: 0.091769) [&95%={0.092, 0.092}]: 0.007627) [&95%={0.099, 0.099}]: 0.009729, Rusdis1: 0.109126) [&95%={0.109, 0.109}]: 0.061282) [&95%={0.170, 0.170}]: 0.112332, (Lacvol1: 0.250648, Lacsub1: 0.250648) [&95%={0.251, 0.251}]: 0.032091) [&95%={0.283, 0.283}]: 0.110440) [&95%={0.393, 0.393}]: 0.066465) [&95%={0.460, 0.460}]: 0.256462) [&95%={0.716, 0.716}]: 0.057682) [&95%={0.774, 0.774}]: 0.288563) [&95%={1.062, 1.062}]: 2.269350) [&95%={3.332, 3.332}]: 0.332955) [&95%={3.665, 3.665}]: 0.874913) [&95%={4.540, 4.540}]; 

Best,

Annie

[Ziheng]

yes, that is strange indeed.  could you copy the control file and the tree file here.  i wonder what kind of calibrations you are using.  perhaps the mcmc is stuck at one point, without moving.

best wishes, 

ziheng

[Annie Lebreton]

Hi Ziheng, 

Thanks for your help,

To get the tree, I first run mcmctree (usedata=3) then codeml, extract the rate estimation and run mcmctree again (usedata=2). Here is the control file and the tree I used for this second run on a dataset of 10 genes. 

<mcmctree_2ndRun.ctl.>

          seed = -1

       seqfile =xdc.phy

      treefile = ../tree_lactarius.tre

       outfile =xdc_mcmctree_2ndRun.out

         ndata =10

       seqtype = 2  * 0: nucleotides; 1:codons; 2:AAs

       usedata = 2    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV

         clock = 2    * 1: global clock; 2: independent rates

       RootAge = '<5.0'  * safe constraint on root age, used if no fossil for root.

     cleandata =1    * remove sites with ambiguity data (1:yes, 0:no)?

       BDparas = 1 1 0.1  * birth, death, sampling

   kappa_gamma = 6 2      * gamma prior for kappa

   alpha_gamma = 1 1      * gamma prior for alpha

   rgene_gamma = 1 13.7557

  sigma2_gamma = 1 4.5

      finetune = 1: .05 .08 .1 .01 .01 .1 * auto (0 or 1): times, musigma2, rates, mixing, paras, FossilErr

         print = 1

        burnin = 20000

      sampfreq = 30

       nsample = 2000

<tree_lactarius.tre> 

34 1

((Mutel1,Sclhys1),(Albpec1,(((Hetan2,Stehi1)'>1.06<0.8085',(Hercor1,Laxbic1)),(((Varmin1,Ricme1),(Amycha1,Echtin1)),((Clapy1,(Lenvul1,Aurvu1)'>0.534<0.37'),(Glocon1,((Lacqui1,(Lacpsa1,(Lacviv1,((Lacaka1,(Lacpse1,Lachen1)),((Lachat1,Lacdel1),Lacsan1))))),(Muloch1,(((Rusbre1,Ruscom1),(((Rusvin1,Rusoch1),(Ruseme1,Rusrug1)),Rusdis1)),(Lacvol1,Lacsub1)))))))'>1.51<1.2041')));

Best,

Annie

[Ziheng]

hi annie,

the > < signs for specifying fossil bounds in the tree file are facing the wrong directions, so you need to fix them.  for example, 

'>1.06<0.8085'

should be 

'<1.06>0.8085'

or 

'>0.8085<1.06'

right now the calibration says that the age is larger than 1.06 but smaller than 0.8085.

the tree in examples/DatingSoftBound/ has 

((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);

which means that the human-chimp divergence time is >0.06 and < 0.08.  

the version of the program i have prints out the following error message: 

Error: fossil bounds in tree incorrect.

but this >1.06<0.8085'did not happen with the version you are using.  perhaps the error message was added later.

sorry about the confusion.

ziheng

[Annie Lebreton]

Dear Ziheng, 

It fixed the problem,  thanks !

Annie