CODEML parameter (p,w) estimate errors

[Yuri]

Hello,

I've been running some tests on synthetic data using CODEML and I've encountered some concerns with omega and category weight estimates when branch lengths are estimated. This problem is reproducible and consistent in PAML 4.10.10. From what I've tested, this does not seem to occur in PAML 4.10.9 with the same datasets. M2a and M8 seem to be most affected.

For example, running M2a (NSsites = 2) I can produce the following estimates from the same dataset, running the same .ctl file.

4.10.10

MLEs of dN/dS (w) for site classes (K=3)

p:   1.00000  0.00000  0.00000
w:   1.00000  1.00000 42.40770

4.10.9

MLEs of dN/dS (w) for site classes (K=3)

p:   0.40383  0.00000  0.59617
w:   0.82781  1.00000  1.99162

[Sandra AC]

Hi Yuri,

Thanks for your message! In order for us to be able to reproduce and troubleshoot this issue, could you please share the control files and input tree and sequence files you have used?

Thanks in advance!
Sandy

[Sandra AC]

Hi there,

Just in case other users face a similar situation than the one described in this thread, please find some details below as to how troubleshooting proceeded -- Yuri has confirmed that all issues have now been resolved!

I analysed Yuri's dataset with PAML v4.10.9 and v4.10.10 and using Yuri's control file and a simplified version (i.e., only variables used in positive selection analyses included such as this control file for site models). The MLEs when calculated with CODEML in both PAML versions and when using Yuri´s control file (purple) and the simplified control file (blue) were the following  (grep 'MLEs' -A4 */out*txt) :

```
paml_41010_myctl/out_4_10_10.txt:MLEs of dN/dS (w) for site classes (K=3)
paml_41010_myctl/out_4_10_10.txt-
paml_41010_myctl/out_4_10_10.txt-p:   0.15012  0.29047  0.55941
paml_41010_myctl/out_4_10_10.txt-w:   0.66042  1.00000  2.02596
paml_41010_myctl/out_4_10_10.txt-
--
paml_41010/out_4_10_10.txt:MLEs of dN/dS (w) for site classes (K=3)
paml_41010/out_4_10_10.txt-
paml_41010/out_4_10_10.txt-p:   0.40383  0.00000  0.59617
paml_41010/out_4_10_10.txt-w:   0.82781  1.00000  1.99162
paml_41010/out_4_10_10.txt-
--
paml_4109_myctl/out_4_10_9.txt:MLEs of dN/dS (w) for site classes (K=3)
paml_4109_myctl/out_4_10_9.txt-
paml_4109_myctl/out_4_10_9.txt-p:   0.15012  0.29047  0.55941
paml_4109_myctl/out_4_10_9.txt-w:   0.66043  1.00000  2.02596
paml_4109_myctl/out_4_10_9.txt-
--
paml_4109/out_4_10_9.txt:MLEs of dN/dS (w) for site classes (K=3)
paml_4109/out_4_10_9.txt-
paml_4109/out_4_10_9.txt-p:   0.40383  0.00000  0.59617
paml_4109/out_4_10_9.txt-w:   0.82781  1.00000  1.99162
paml_4109/out_4_10_9.txt-
```

The results were consistent across runs using the same control file and v4.10.9 or v4.10.10. It seems that the option set for variable "method" was the one driving the main differences:

• "method = 0" optimises all parameters simultaneously. The algorithm is slower per iteration but explores the parameter space more thoroughly and seems to be more reliable at finding a better optimum for site models like M2a. This is the option enabled by default if variable "method" is not included in the control file.
• "method = 1" optimises one branch at a time. This algorithm is much faster but can be more prone to stopping at a worse optimum (in this test, it converged to a lower log-likelihood [−185909.862  for "method = 1" vs  −185908.789 for "method = 0"] with the neutral site-class proportion collapsing to zero). This is the option that was set in Yuri's control file.

To test robustness in convergence, the best procedure would be to change the starting values for variable omega (e.g., set variable "omega" to different values) to check whether independent runs converge to the same point. For tests of positive selection, using "method = 0" (default, does not need to be specified in the control file) may help in reaching convergence.

Hope this helps!

All the best,
Sandy