MrBayes icon indicating copy to clipboard operation
MrBayes copied to clipboard

Can the MrBayes estimates of divergence time go out of prior bounds?

Open CalamusLiber opened this issue 8 years ago • 1 comments

Dear users,

I tried to use MrBayes 3.2.6 (MPI) to date a fixed phylogeny. I set the priors as uniform with lower and/or upper bounds. But some of the resulting ages (including 95% HPD) present lower than the lower bound of priors. I paste my commands and the resulting chronogram afterward. Hope someone can help me. Thanks!

===== MB commands (passing data block) ===== [data set]

[tree topology] begin trees; tree huichong = [&R] ((((((((((((((RAP_Hysterothylacium_fabri,RAP_Hysterothylacium_spn),RAP_Hysterothylacium_aduncum),RAP_Hysterothylacium_reliquens),((RAP_Raphidascaris_longispicula,RAP_Raphidascaris_lophii),RAP_Hysterothylacium_longilabrum)),((RAP_Hysterothylacium_liparis,RAP_Hysterothylacium_sinense),RAP_Hysterothylacium_thalassini)),RAP_Raphidascaroides_nipponensis),(((((RAP_Raphidascaroides_moraveci,RAP_Raphidascaroides_brasiliensis),RAP_Hysterothylacium_pelagicum),RAP_Maricostula_tetrapteri),(RAP_Hysterothylacium_deardorffoverstreetorum,RAP_Raphidascaris_acus)),(((RAP_Hysterothylacium_amoyense,RAP_Hysterothylacium_zhoushanense),RAP_Goezia_pelagia),RAP_Iheringascaris_inquies))),(((((ASC_Baylisascaris_schroederi,ASC_Baylisascaris_ailuri),ASC_Baylisascaris_transfuga),(ASC_Baylisascaris_columnaris,ASC_Baylisascaris_procyonis)),((ASC_Ascaris_lumbricoides,ASC_Ascaris_suum),ASC_Parascaris_equorum)),ASC_Toxascaris_leonina)),((((ACA_Mawsonascaris_spn,ACA_Mawsonascaris_myliobatum),ACA_Mawsonascaris_zhoui),(ACA_Acanthocheilus_rotundatus,ACA_Pseudanisakis_rajae)),(((ASC_Toxocara_vitulorum,ASC_Toxocara_cati),ASC_Toxocara_canis),(ASC_Porrocaecum_depressum,ASC_Porrocaecum_reticulatum)))),(((((((((ANI_Anisakis_berlandi,ANI_Anisakis_simplex),ANI_Anisakis_pegreffii),ANI_Anisakis_ziphidarum),ANI_Anisakis_typica),((ANI_Anisakis_brevispiculata,ANI_Anisakis_paggiae),ANI_Anisakis_physeteris)),ANI_Terranova_caballeroi),((ANI_Pseudoterranova_azarasi,ANI_Pseudoterranova_cattani),ANI_Pseudoterranova_decipiens)),ANI_Pulchrascaris_chiloscyllii),(((((ANI_Contracaecum_ogmorhini,ANI_Contracaecum_rudolphii),ANI_Contracaecum_bioccai),ANI_Contracaecum_multipapillatum),ANI_Contracaecum_microcephalum),((ANI_Contracaecum_osculatum,ANI_Phocascaris_cystophorae),ANI_Contracaecum_radiatum)))),(((HETC_Dujardinascaris_gigantea,HETC_Ortleppascaris_sinensis),HETC_Krefftascaris_sharpiloi),HETC_Heterocheilus_tunicatus)),OGCOS_Oxyascaris_sp),OGCOS_Cruzia_americana),OGHETK_Ascaridia_galli); end;

begin mrbayes; charset ITS = 1-1382; charset 28S = 1383-2152; charset 18S = 2153-3908; charset COICOII = 3909-4948 4949-5540; charset 12S = 5541-6069; partition matrices = 5: ITS, 28S, 18S, COICOII, 12S; set partition = matrices;

lset applyto=(1) nst=6 rates=invgamma;
lset applyto=(2) nst=2 rates=invgamma;
lset applyto=(3) nst=2 rates=invgamma;
prset applyto=(3) statefreqpr=fixed(equal);
lset applyto=(4) nst=2 rates=invgamma;
lset applyto=(5) nst=6 rates=invgamma;

prset applyto=(all) ratepr=variable;
unlink statefreq=(all) revmat=(all) shape=(all) pinvar=(all) tratio=(all);

[all the constrained clades are monophyletic nodes in the tree above]
constraint RAP = RAP_Goezia_pelagia RAP_Hysterothylacium_aduncum RAP_Hysterothylacium_amoyense RAP_Hysterothylacium_deardorffoverstreetorum RAP_Hysterothylacium_fabri RAP_Hysterothylacium_liparis RAP_Hysterothylacium_longilabrum RAP_Hysterothylacium_pelagicum RAP_Hysterothylacium_reliquens RAP_Hysterothylacium_sinense RAP_Hysterothylacium_spn RAP_Hysterothylacium_thalassini RAP_Hysterothylacium_zhoushanense RAP_Iheringascaris_inquies RAP_Maricostula_tetrapteri RAP_Raphidascaris_acus RAP_Raphidascaris_longispicula RAP_Raphidascaris_lophii RAP_Raphidascaroides_brasiliensis RAP_Raphidascaroides_moraveci RAP_Raphidascaroides_nipponensis;
constraint ACA_Mawsonascaris = ACA_Mawsonascaris_myliobatum ACA_Mawsonascaris_spn ACA_Mawsonascaris_zhoui;
constraint ANI_Anisakis = ANI_Anisakis_berlandi ANI_Anisakis_brevispiculata ANI_Anisakis_paggiae ANI_Anisakis_pegreffii ANI_Anisakis_physeteris ANI_Anisakis_simplex ANI_Anisakis_typica ANI_Anisakis_ziphidarum;
constraint ANI_Pseudoterranova = ANI_Pseudoterranova_azarasi ANI_Pseudoterranova_cattani ANI_Pseudoterranova_decipiens;

calibrate
    RAP = unif(151.2,250.0)
    ACA_Mawsonascaris = unif(170.3,1E4) [**NOTE: I used 1E4 instead of unlimited upper bound. The resulting age of this node (MRCA) is only ~30 Mya!**]
    ANI_Anisakis = unif(33.9,52.4)
    ANI_Pseudoterranova = unif(28.5,50)
    ;
   
prset brlenspr=clock:uniform;
prset treeagepr = unif(0,541);
prset nodeagepr=calibrated;
prset ratepr=variable topologypr=fixed(huichong);
prset clockratepr=lognorm(-7,0.6);
prset clockvarpr=igr;
prset igrvarpr=exp(10);
   
set usebeagle=yes beagledevice=cpu beagleprecision=double;
set beaglescaling=dynamic beaglesse=yes;

mcmc ngen= 10000000 printfreq=1000 samplefreq=1000 nruns=4 nchains=5 savebrlens=yes;

end;

======== resulting tree (red box show the clade with wrong estimates) image

CalamusLiber avatar Jan 16 '17 10:01 CalamusLiber

type help constraint in MrBayes, the message contains:

It is important to note that simply defining a constraint using this          
   command is not sufficient for the program to actually implement the           
   constraint in an analysis. You must also enforce the constraints using        
   'prset topologypr = constraints (<list of constraints>)'. For more infor-     
   mation on this, see the help on the 'prset' command.  

but topologypr=fixed(huichong) was used in the analysis instead, meaning that the constraints were not enforced. This explains why the node ages are out of bounds.

There seems no proper way to fix the topology and define constraints at the same time using the topologypr command. a workaround might be using the following commands:

prset topologypr=constraint(RAP,ACA_Mawsonascaris,ANI_Anisakis,ANI_Pseudoterranova);
mcmcp ngen= 10000000 printfreq=1000 samplefreq=1000 nruns=4 nchains=5 savebrlens=yes;
startvals Tau=huichong;  [starting tree]
propset ExtSprClock(Tau,V)$prob=0;  [disable topology moves]
propset NNIClock(Tau,V)$prob=0;
propset ParsSPRClock(Tau,V)$prob=0;
mcmc;

zhangchicool avatar Dec 04 '18 06:12 zhangchicool