Bio.Phylo.Paml.Codeml'S Results Parser Quietly Fails To Read All The Data

Question

[a1ultima]

Biopython comes with methods to interface with the PAML package for phylogenetic analysis.

In particular I am using Bio.Phylo.PAML to run analyses using PAML's codeml.exe program which in my case does Ka/Ks (dN/dS) ratio analysis on pairs of orthologous gene sequences.

After running the analysis using results = cml.run() I can see it has successfully generated result.out files that look about right. Most importantly is the final line of the file that I need to parse into Python:

t= 0.2173  S=   703.9  N=  1489.1  dN/dS=  0.2247  dN = 0.0344  dS = 0.1529

What I am most interested in is dN/dS = 0.2247

According to Biopython's PAML wiki this value can be obtained from Python by doing results = cml.run() which generates a dictionary with a set of values I am interested in after running the analysis. The wiki claims I can find the values I need in a key called 'parameters'. But this only returns one of the values I need t= 0.2173, look:

>>> results.values()
['Fcodon', 'One dN/dS ratio for branches, ', '4.7b', {0: {'description': 'one-ratio', 'parameters': {'t': 0.1982}}}, {'htlv': {}, 'stlv': {}}]

Notice, that my parameters key only contains the t= 0.2173 and has omitted S= 703.9 N= 1489.1 dN/dS= 0.2247 dN = 0.0344 dS = 0.1529

Could anybody with codeml experience explain to me why the parser fails to yield most of the parameters (values) I am interested in?

Thank you

[a1ultima]

Been hacking away at this problem for a while now, what is quite mysterious is that the PAML manual is given as version 4.7a (2013 May) whilst the codeml program I am using states it is PAML 4.7b (2013 Sep)... hmm perhaps this could be a very new bug caused by Biopython not yet being updated to deal with the new codeml version?

[Peter]

Could be a bug - why don't you report this to Biopython with a sample file and snippet of code to show the problem? https://github.com/biopython/biopython/issues

[a1ultima]

Don, see issues being solved here: https://github.com/biopython/biopython/issues/483#issuecomment-78052449

[Peter]

Did you report a bug? Could you post a link to it?

[a1ultima]

Yep here (sorry for the late reply):

https://github.com/biopython/biopython/issues/483

Answer 1

[AW]

Hi,

You can obtain omega values for each branch using the following commands:

from Bio.Phylo.PAML import codeml

results = codeml.read(paml_outputfile)

print results["NSsites"][0]["parameters"]["omega"]

This gives you a list of omega (dn/ds) for each branch

I am using version 4.7b and biopython-1.6.3

Hope that helps!

A