Skip to main content
3 of 3
added 255 characters in body
M__
  • 13k
  • 5
  • 29
  • 46

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.

M__
  • 13k
  • 5
  • 29
  • 46