class Bio::PAML::Codeml::Report

Description¶ ↑

Run PAML codeml and get the results from the output file. The Codeml::Report object is returned by Bio::PAML::Codeml.query. For example

codeml = Bio::PAML::Codeml.new('codeml', :runmode => 0, 
    :RateAncestor => 1, :alpha => 0.5, :fix_alpha => 0)
result = codeml.query(alignment, tree)

where alignment and tree are Bioruby objects. This class assumes we have a buffer containing the output of codeml.

References¶ ↑

Phylogenetic Analysis by Maximum Likelihood (PAML) is a package of programs for phylogenetic analyses of DNA or protein sequences using maximum likelihood. It is maintained and distributed for academic use free of charge by Ziheng Yang. Suggestion citation

Yang, Z. 1997
PAML: a program package for phylogenetic analysis by maximum likelihood 
CABIOS 13:555-556

abacus.gene.ucl.ac.uk/software/paml.html

Examples¶ ↑

Invoke Bioruby’s PAML codeml parser, after having read the contents of the codeml result file into buf (for example using File.read)

>> c = Bio::PAML::Codeml::Report.new(buf)

Do we have two models?

>> c.models.size
=> 2
>> c.models[0].name
=> "M0"
>> c.models[1].name
=> "M3"

Check the general information

>> c.num_sequences
=> 6
>> c.num_codons
=> 134
>> c.descr
=> "M0-3"

Test whether the second model M3 is significant over M0

>> c.significant
=> true

Now fetch the results of the first model M0, and check its values

>> m0 = c.models[0]
>> m0.tree_length
=> 1.90227
>> m0.lnL
=> -1125.800375
>> m0.omega
=> 0.58589
>> m0.dN_dS
=> 0.58589
>> m0.kappa
=> 2.14311
>> m0.alpha
=> nil

We also have a tree (as a string)

>> m0.tree
=> "((((PITG_23265T0: 0.000004, PITG_23253T0: 0.400074): 0.000004, PITG_23257T0: 0.952614): 0.000004, PITG_23264T0: 0.445507): 0.000004, PITG_23267T0: 0.011814, PITG_23293T0: 0.092242);"

Check the M3 and its specific values

>> m3 = c.models[1]
>> m3.lnL
=> -1070.964046
>> m3.classes.size
=> 3
>> m3.classes[0]
=> {:w=>0.00928, :p=>0.56413}

And the tree

>> m3.tree
=> "((((PITG_23265T0: 0.000004, PITG_23253T0: 0.762597): 0.000004, PITG_23257T0: 2.721710): 0.000004, PITG_23264T0: 0.924326): 0.014562, PITG_23267T0: 0.000004, PITG_23293T0: 0.237433);"

Next take the overall posterior analysis

>> c.nb_sites.size
=> 44
>> c.nb_sites[0].to_a
=> [17, "I", 0.988, 3.293]

or by field

>> codon = c.nb_sites[0]
>> codon.position
=> 17
>> codon.probability
=> 0.988
>> codon.dN_dS
=> 3.293

with aliases

>> codon.p
=> 0.988
>> codon.w
=> 3.293

Now we generate special string ‘graph’ for positive selection. The following returns a string the length of the input alignment and shows the locations of positive selection:

>> c.nb_sites.graph[0..32]
=> "                **    *       * *"

And with dN/dS (high values are still an asterisk *)

>> c.nb_sites.graph_omega[0..32]
=> "                3*    6       6 2"

We also provide the raw buffers to adhere to the principle of unexpected use. Test the raw buffers for content:

>> c.header.to_s =~ /seed/
=> 1
>> m0.to_s =~ /one-ratio/
=> 3
>> m3.to_s =~ /discrete/
=> 3
>> c.footer.to_s =~ /Bayes/
=> 16

Finally we do a test on an M7+M8 run. Again, after loading the results file into buf

Invoke Bioruby’s PAML codeml parser

>> c = Bio::PAML::Codeml::Report.new(buf78)

Do we have two models?

>> c.models.size
=> 2
>> c.models[0].name
=> "M7"
>> c.models[1].name
=> "M8"

Assert the results are significant

>> c.significant
=> true

Compared to M0/M3 there are some differences. The important ones are the parameters and the full Bayesian result available for M7/M8. This is the naive Bayesian result:

>> c.nb_sites.size
=> 10

And this is the full Bayesian result:

>> c.sites.size
=> 30
>> c.sites[0].to_a
=> [17, "I", 0.672, 2.847]
>> c.sites.graph[0..32]
=> "                **    *       * *"

Note the differences of omega with earlier M0-M3 naive Bayesian analysis:

>> c.sites.graph_omega[0..32]
=> "                24    3       3 2"

The locations are the same, but the omega differs.