PAML branch-site model: Identical lnL value between null and alternative models?

Question

I have been running PAML, via codeml, to try to identify positive selection in certain branches of my tree for one gene (but will eventually be trying to automate the analysis over multiple genes). This was done by setting model = 2 and NSsites = 2 in the codeml.ctl file. I attempted to run both the null and alternative models by setting fix_omega = 1 in the null.ctl and fix_omega = 0 in the alternative.ctl.

However, after running both models, they have returned exactly identical log-likelihood (lnL) values. I wasn't expecting to find anything significant, but the fact that they were identical struck me as odd. Out of curiosity, I altered the branch of interest in the newick tree file and reattempted the analysis, but again, the lnL values remain the same between the null and alternative models (although slightly different to the value returned with the different branch selected).

I have since repeated this many times, choosing a different branch(es) each time to see if I can get a difference between the null and alternative, but no luck. Am I missing something obvious? Should I really be getting the same exact result between null and alternative models, regardless of where in the tree I place #1? To reiterate, I wasn't expecting to find significant results, but was at least expecting (non-significantly) different numbers for lnL...

Answer 1

The issue with PAML is that it doesn't throw many exceptions (programming jargon for trapping and reporting errors) if something is wrong, generally the "exceptions" are identified in the failure to report an expected value in the output.

Possible explanation fix_omega is not turning the omega parameter on or off it is simply defining what the value of omega will be. If the default value, which is set at .4, but can be changed, is very close to the estimated value the likelihoods will be similar - or even identical.

Goto debugging As computer/data scientists, the very first thing we should do is just check there is run time variation between the options you are assessing via,

time codeml ./file

The thing to remember is fix_omega=0 is where Omega is being estimated from the data, whilst =1 means it will be supplied with a value. There should be longer run time on the "0" option.

Nitty gritty debugging

1 Please check you have

seqtype = 3

If not, then that is your answer, omega was never being set.

2 For debugging only make sure you have,

model = 1

This is global dN/dS for all the data not for a specific clade within the tree, the reason for this is given in the final step below

3, Then try a different NSsites model,

NSsites = 8

This should produce real differences. The NSsites options are,

• 0:one w;1:neutral;2:selection; 3:discrete;4:freqs; * 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
• 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1; * 13:3normal>0

w here is for omega

4 The proof of your system is that when you

fix_omega = 1
omega = .4 * try lots of different values, e.g. 0.3, 0.4, 0,7, 1, 1.5, 2

This should produce likelihood variation

---

Finally

model = 2

Then go to your tree file and make sure that the nomenclature used follows the manual, I think from memory the # is deployed. The output file should change to specify that you have two values for dN/dS from your treefile, thus the output from model=2 should vary from the output of model=1 if not the delineation has been ignored. If you do see two rates, great and now you are into parameter exploration. Parameter exploration is basically the options 1 to 4 described above, plus the stuff you mentioned.

---

Oh you mentioned using Newick format. I haven't used PAML for a while (dNdS stuff is ancient) but it it only used phylip format when I used it. I suspect if it got the wrong tree format PAML will throw an exception, but it's worth checking.