vlanore / diffsel

A C++ program to detect convergent evolution using Bayesian phylogenetic codon models.
Other
6 stars 2 forks source link

MCMC convergence checking #3

Closed kanglizhu closed 5 years ago

kanglizhu commented 5 years ago

According to the updated diffsel README, I have used samhd1_short.ali and samhd1.tree labelled with three "conditions" to conduct MCMC convergence checking and post-analysis. But I met some problems with MCMC convergence checking.

According to the paper “Detecting adaptive convergent amino acid evolution”, I know we can use two different approaches to establish MCMC convergence: the original method compares two MCMC chains using the tracecomp program, and the new method uses the Raftery and Lewis’s Diagnostic implemented in the R package coda after 200 iterations to estimate the number of necessary iterations, then run as many iterations as 120% of the estimated number and finally perform the same diagnostic to check convergence (as written in ‘diffselMCMCConvergenceAnalysis.Rmd’ script). For the original “ tracecomp” method, the chain has converged well if all lines have a large effsize (>300) and a small rel_diff (<0.1). For the “Raftery and Lewis” method, no error message printed means the chain has converged (is it right? ).

Using the samhd1_short.ali and samhd1.tree with 3 conditions, I used both methods to establish MCMC convergence. When the iterations in “Raftery and Lewis” method works well, the effsize of all lines in “ tracecomp” method is only about 5. I’m confused about this result, and I don’t know which step is wrong. So I have attached my result, sequence file and tree file, can you please help me check my results? check mcmc convergece .pdf samhd1.zip

vlanore commented 5 years ago

Hello,

There is a bug in tracecomp where it outputs wrong results when no burn-in is specified via command line.

To work around this, just manually specify your burn-in with -x <burn-in>. For example in your case, this should do the trick:

tracecomp -x 240 samhdi_run_1200_1 samhdi_run_1200_2

Let me know if this fixes your problem.

I'm sorry about this bug. tracecomp is an old piece of code that we haven't really looked into for a long time.

kanglizhu commented 5 years ago

Hello,

There is a bug in tracecomp where it outputs wrong results when no burn-in is specified via command line.

To work around this, just manually specify your burn-in with -x <burn-in>. For example in your case, this should do the trick:

tracecomp -x 240 samhdi_run_1200_1 samhdi_run_1200_2

Let me know if this fixes your problem.

I'm sorry about this bug. tracecomp is an old piece of code that we haven't really looked into for a long time.

Hello,

When i specify burn-in via command line, it works well now. Thank you so much for your help.

kanglizhu commented 5 years ago

Hello,

When i specify burn-in via command line, it works well. Thank you so much for your help.

发件人: Vincent Lanore [mailto:notifications@github.com] 发送时间: 2019年8月20日 21:27 收件人: vlanore/diffsel diffsel@noreply.github.com 抄送: Kangli ZHU 朱康丽 zhukangli@westlake.edu.cn; Author author@noreply.github.com 主题: Re: [vlanore/diffsel] MCMC convergence checking (#3)

Hello,

There is a bug in tracecomp where it outputs wrong results when no burn-in is specified via command line.

To work around this, just manually specify your burn-in with -x . For example in your case, this should do the trick:

tracecomp -x 240 samhdi_run_1200_1 samhdi_run_1200_2

Let me know if this fixes your problem.

I'm sorry about this bug. tracecomp is an old piece of code that we haven't really looked into for a long time.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/vlanore/diffsel/issues/3?email_source=notifications&email_token=AM6NJRBJYBQZHGDMPKXMZVLQFPWILA5CNFSM4INWNC3KYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4WI6RA#issuecomment-523013956, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AM6NJRDR7N4GXBM4DVY6O33QFPWILANCNFSM4INWNC3A.