[Bioperl-l] PAML

[Jason Stajich]

Please always mail the list, I can't always answer questions.  It is
also much more helpful if you were to send a result file.

I think you want:

if( my $result = $paml_parser->next_result() ) {
  for my $ns_result ( $result->get_NSSite_results ) {
    print "model num is ", $ns_result-> model_num, "\n";
    print "lnL: ", $ns_result->likelihood, "\n";
    print "kappa: ", $ns_result->kappa, "\n";
  }
}

 If this works can you add something to the wiki HOWTO for PAML of the
basic code you are using and how you are doing the loglikehood test
comparison - you can get the number of free params for the chisq from
another field in the files as well - I wrote a script that did all
this a while ago but seem to have lost it. Be good to get that sort of
working code back into the wiki and/or repo for people to use and play
with.
-jason

--
Jason Stajich
jason at bioperl.org
jason.stajich at gmail.com
http://bioperl.org/wiki/User:Jason
http://twitter.com/hyphaltip



On Wed, Apr 7, 2010 at 3:45 AM, Dave Burt <david.burt1955 at o2.co.uk> wrote:
> Hi Jason,
>
>
>
> I am using Bioperl and various PAML modules - all is working except I have
> one question.
>
>
>
> How do I retrieve the likelihood (lnL) and kappa values from an PAML
> analysis?
>
>
>
> For example, I would like to compare model 0 vs. model 1, then repeat this
> for 1000's of genes.
>
>
>
> Here's the ctl file from a model 1 analysis
>
>
>
> Thanks Dave
>
>
>
>
>
> seqfile = c:/cygwin/tmp/6Apr3EZ7Wa/9g9rUwnO2d
>
> outfile = mlc
>
> treefile = c:/cygwin/tmp/6Apr3EZ7Wa/WcFjpKcmzU
>
> aaDist = 0
>
> verbose = 1
>
> noisy = 3
>
> RateAncestor = 1
>
> kappa = 2
>
> model = 1
>
> ndata = 1
>
> aaRatefile = wag.dat
>
> Small_Diff = .5e-6
>
> CodonFreq = 2
>
> runmode = 0
>
> alpha = 0
>
> omega = 0.4
>
> fix_kappa = 0
>
> Mgene = 0
>
> method = 0
>
> fix_omega = 0
>
> getSE = 0
>
> NSsites = 0
>
> seqtype = 1
>
> cleandata = 1
>
> icode = 0
>
> fix_alpha = 1
>
> clock = 0
>
> ncatG = 8
>
> Malpha = 0
>
> fix_blength = 0
>
>
>
>