File Coverage

blib/lib/Bio/Tools/Run/Phylo/PAML/Yn00.pm
Criterion Covered Total %
statement 28 125 22.4
branch 2 48 4.1
condition 0 18 0.0
subroutine 9 18 50.0
pod 9 9 100.0
total 48 218 22.0


line stmt bran cond sub pod time code
1             # $Id$
2             #
3             # BioPerl module for Bio::Tools::Run::Phylo::PAML::Yn00
4             #
5             # Please direct questions and support issues to
6             #
7             # Cared for by Jason Stajich
8             #
9             # Copyright Jason Stajich
10             #
11             # You may distribute this module under the same terms as perl itself
12              
13             # POD documentation - main docs before the code
14              
15             =head1 NAME
16              
17             Bio::Tools::Run::Phylo::PAML::Yn00 - Wrapper aroud the PAML program yn00
18              
19             =head1 SYNOPSIS
20              
21             use Bio::Tools::Run::Phylo::PAML::Yn00;
22             use Bio::AlignIO;
23             my $alignio = Bio::AlignIO->new(-format => 'phylip',
24             -file => 't/data/gf-s85.phylip');
25             my $aln = $alignio->next_aln;
26              
27             my $yn = Bio::Tools::Run::Phylo::PAML::Yn00->new();
28             $yn->alignment($aln);
29             my ($rc,$parser) = $yn->run();
30             while( my $result = $parser->next_result ) {
31             my @otus = $result->get_seqs();
32             my $MLmatrix = $result->get_MLmatrix();
33             # 0 and 1 correspond to the 1st and 2nd entry in the @otus array
34             my $dN = $MLmatrix->[0]->[1]->{dN};
35             my $dS = $MLmatrix->[0]->[1]->{dS};
36             my $kaks =$MLmatrix->[0]->[1]->{omega};
37             print "Ka = $dN Ks = $dS Ka/Ks = $kaks\n";
38             }
39              
40             =head1 DESCRIPTION
41              
42             This is a wrapper around the yn00 (method of Yang and Nielsen, 2000)
43             program of PAML (Phylogenetic Analysis by Maximum Likelihood) package
44             of Ziheng Yang. See http://abacus.gene.ucl.ac.uk/software/paml.html
45             for more information.
46              
47             This module will generate a proper yn00.ctl file and will run the
48             program in a separate temporary directory to avoid creating temp files
49             all over the place and will cleanup after itself.
50              
51             =head1 FEEDBACK
52              
53             =head2 Mailing Lists
54              
55             User feedback is an integral part of the evolution of this and other
56             Bioperl modules. Send your comments and suggestions preferably to
57             the Bioperl mailing list. Your participation is much appreciated.
58              
59             bioperl-l@bioperl.org - General discussion
60             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
61              
62             =head2 Support
63              
64             Please direct usage questions or support issues to the mailing list:
65              
66             I
67              
68             rather than to the module maintainer directly. Many experienced and
69             reponsive experts will be able look at the problem and quickly
70             address it. Please include a thorough description of the problem
71             with code and data examples if at all possible.
72              
73             =head2 Reporting Bugs
74              
75             Report bugs to the Bioperl bug tracking system to help us keep track
76             of the bugs and their resolution. Bug reports can be submitted via the
77             web:
78              
79              
80             =head1 AUTHOR - Jason Stajich
81              
82             Email jason-at-bioperl.org
83              
84             =head1 CONTRIBUTORS
85              
86             Additional contributors names and emails here
87              
88             =head1 APPENDIX
89              
90             The rest of the documentation details each of the object methods.
91             Internal methods are usually preceded with a _
92              
93             =cut
94              
95              
96             # Let the code begin...
97              
98              
99             package Bio::Tools::Run::Phylo::PAML::Yn00;
100 1     1   795 use vars qw(@ISA %VALIDVALUES $MINNAMELEN $PROGRAMNAME $PROGRAM);
  1         1  
  1         58  
101 1     1   3 use strict;
  1         2  
  1         15  
102 1     1   2 use Cwd;
  1         1  
  1         40  
103 1     1   3 use Bio::Root::Root;
  1         1  
  1         12  
104 1     1   3 use Bio::AlignIO;
  1         1  
  1         11  
105 1     1   2 use Bio::TreeIO;
  1         1  
  1         11  
106 1     1   3 use Bio::Tools::Run::WrapperBase;
  1         1  
  1         16  
107 1     1   2 use Bio::Tools::Phylo::PAML;
  1         2  
  1         113  
108              
109             @ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase);
110              
111              
112             =head2 Default Values
113              
114             See the L module for
115             documentation of the default values.
116              
117             =cut
118              
119             BEGIN {
120              
121 1     1   2 $MINNAMELEN = 25;
122 1 50       5 $PROGRAMNAME = 'yn00' . ($^O =~ /mswin/i ?'.exe':'');
123 1 50       3 if( defined $ENV{'PAMLDIR'} ) {
124 0         0 $PROGRAM = Bio::Root::IO->catfile($ENV{'PAMLDIR'},$PROGRAMNAME);
125             }
126             # valid values for parameters, the default one is always
127             # the first one in the array
128             # much of the documentation here is lifted directly from the codeml.ctl
129             # example file provided with the package
130             %VALIDVALUES = (
131 1         813 'noisy' => [ 0..3,9],
132             'verbose' => [ 0,1,2], # 0:concise, 1:detailed, 2:too much
133              
134             'weighting' => [0,1], # weighting pathways between codons
135             'commonf3x4' => [0,1], # use same f3x4 for all sites
136              
137             # (icode) genetic code
138             # 0:universal code
139             # 1:mamalian mt
140             # 2:yeast mt
141             # 3:mold mt,
142             # 4:invertebrate mt
143             # 5:ciliate nuclear
144             # 6:echinoderm mt
145             # 7:euplotid mt
146             # 8:alternative yeast nu.
147             # 9:ascidian mt
148             #10:blepharisma nu
149             # these correspond to 1-11 in the genbank transl table
150            
151             'icode' => [ 0..10],
152             'ndata' => [1..10],
153             );
154             }
155              
156              
157             =head2 program_name
158              
159             Title : program_name
160             Usage : $yn00->program_name()
161             Function: holds the program name
162             Returns: string
163             Args : None
164              
165             =cut
166              
167             sub program_name {
168 0     0 1   return $PROGRAMNAME;
169             }
170              
171             =head2 program_dir
172              
173             Title : program_dir
174             Usage : $yn00->program_dir()
175             Function: returns the program directory, obtained from ENV variable.
176             Returns: string
177             Args :
178              
179             =cut
180              
181             sub program_dir {
182 0 0   0 1   return Bio::Root::IO->catfile($ENV{PAMLDIR}) if $ENV{PAMLDIR};
183             }
184              
185             =head2 new
186              
187             Title : new
188             Usage : my $obj = Bio::Tools::Run::Phylo::PAML::Yn00->new();
189             Function: Builds a new Bio::Tools::Run::Phylo::PAML::Yn00 object
190             Returns : Bio::Tools::Run::Phylo::PAML::Yn00
191             Args : -alignment => the L object
192             -save_tempfiles => boolean to save the generated tempfiles and
193             NOT cleanup after onesself (default FALSE)
194              
195             =cut
196              
197             sub new {
198 0     0 1   my($class,@args) = @_;
199              
200 0           my $self = $class->SUPER::new(@args);
201 0           my ($aln,$st) = $self->_rearrange([qw(ALIGNMENT SAVE_TEMPFILES)],
202             @args);
203 0 0         defined $aln && $self->alignment($aln);
204 0 0         defined $st && $self->save_tempfiles($st);
205            
206 0           $self->set_default_parameters();
207 0           return $self;
208             }
209              
210             =head2 run
211              
212             Title : run
213             Usage : $yn->run();
214             Function: run the yn00 analysis using the default or updated parameters
215             the alignment parameter must have been set
216             Returns : 3 values,
217             $rc = 1 for success, 0 for errors
218             hash reference of the Yang calculated Ka/Ks values
219             this is a set of pairwise observations keyed as
220             sequencenameA->sequencenameB->datatype
221             hash reference same as the previous one except it for the
222             Nei and Gojobori calculated Ka,Ks,omega values
223             Args : none
224              
225              
226             =cut
227              
228             sub run{
229 0     0 1   my ($self,$aln) = @_;
230 0   0       ($aln) ||= $self->alignment();
231 0 0         if( ! $aln ) {
232 0           $self->warn("must have supplied a valid alignment file in order to run yn00");
233 0           return 0;
234             }
235 0           my ($tmpdir) = $self->tempdir();
236 0           my ($tempseqFH,$tempseqfile);
237 0 0 0       if( ! ref($aln) && -e $aln ) {
238 0           $tempseqfile = $aln;
239             } else {
240 0 0         ($tempseqFH,$tempseqfile) = $self->io->tempfile
241             ('-dir' => $tmpdir,
242             UNLINK => ($self->save_tempfiles ? 0 : 1));
243 0 0         my $alnout = Bio::AlignIO->new('-format' => 'phylip',
244             '-fh' => $tempseqFH,
245             '-interleaved' => 0,
246             '-idlength' => $MINNAMELEN > $aln->maxdisplayname_length() ? $MINNAMELEN : $aln->maxdisplayname_length() +1);
247            
248 0           $alnout->write_aln($aln);
249 0           $alnout->close();
250 0           undef $alnout;
251 0           close($tempseqFH);
252 0           undef $tempseqFH;
253             }
254             # now let's print the yn.ctl file.
255             # many of the these programs are finicky about what the filename is
256             # and won't even run without the properly named file. Ack
257            
258 0           my $yn_ctl = "$tmpdir/yn00.ctl";
259 0 0         open(YN, ">$yn_ctl") or $self->throw("cannot open $yn_ctl for writing");
260 0           print YN "seqfile = $tempseqfile\n";
261              
262 0           my $outfile = $self->outfile_name;
263              
264 0           print YN "outfile = $outfile\n";
265 0           my %params = $self->get_parameters;
266 0           while( my ($param,$val) = each %params ) {
267 0           print YN "$param = $val\n";
268             }
269 0           close(YN);
270 0           my ($rc,$parser) = (1);
271             {
272 0           my $cwd = cwd();
  0            
273 0           my $exit_status;
274 0           chdir($tmpdir);
275 0           my $ynexe = $self->executable();
276 0 0         $self->throw("unable to find executable for 'yn'") unless $ynexe;
277 0           open(RUN, "$ynexe |");
278 0           my @output = ;
279 0           $exit_status = close(RUN);
280 0           $self->error_string(join('',@output));
281 0 0 0       if( (grep { /\berr(or)?: /io } @output) || !$exit_status ) {
  0            
282 0           $self->warn("There was an error - see error_string for the program output");
283 0           $rc = 0;
284             }
285 0           eval {
286 0           $parser = Bio::Tools::Phylo::PAML->new(-file => "$tmpdir/mlc",
287             -dir => "$tmpdir");
288              
289             };
290 0 0         if( $@ ) {
291 0           $self->warn($self->error_string);
292             }
293 0           chdir($cwd);
294             }
295 0 0         if( $self->verbose > 0 ) {
296 0           open(IN, "$tmpdir/mlc");
297 0           while() {
298 0           $self->debug($_);
299             }
300             }
301            
302 0 0         unless ( $self->save_tempfiles ) {
303 0           unlink("$yn_ctl");
304 0           $self->cleanup();
305             }
306 0           return ($rc,$parser);
307             }
308              
309             =head2 error_string
310              
311             Title : error_string
312             Usage : $obj->error_string($newval)
313             Function: Where the output from the last analysus run is stored.
314             Returns : value of error_string
315             Args : newvalue (optional)
316              
317              
318             =cut
319              
320             sub error_string{
321 0     0 1   my ($self,$value) = @_;
322 0 0         if( defined $value) {
323 0           $self->{'error_string'} = $value;
324             }
325 0           return $self->{'error_string'};
326              
327             }
328              
329             =head2 alignment
330              
331             Title : alignment
332             Usage : $codeml->align($aln);
333             Function: Get/Set the L object
334             Returns : L object
335             Args : [optional] L
336             Comment : We could potentially add support for running directly on a file
337             but we shall keep it simple
338             See also: L
339              
340             =cut
341              
342             sub alignment{
343 0     0 1   my ($self,$aln) = @_;
344 0 0         if( defined $aln ) {
345 0 0 0       if( !ref($aln) || ! $aln->isa('Bio::Align::AlignI') ) {
346 0           $self->warn("Must specify a valid Bio::Align::AlignI object to the alignment function");
347 0           return undef;
348             }
349 0           $self->{'_alignment'} = $aln;
350             }
351 0           return $self->{'_alignment'};
352             }
353              
354             =head2 get_parameters
355              
356             Title : get_parameters
357             Usage : my %params = $self->get_parameters();
358             Function: returns the list of parameters as a hash
359             Returns : associative array keyed on parameter names
360             Args : none
361              
362              
363             =cut
364              
365             sub get_parameters{
366 0     0 1   my ($self) = @_;
367             # we're returning a copy of this
368 0           return %{ $self->{'_codemlparams'} };
  0            
369             }
370              
371              
372             =head2 set_parameter
373              
374             Title : set_parameter
375             Usage : $codeml->set_parameter($param,$val);
376             Function: Sets a codeml parameter, will be validated against
377             the valid values as set in the %VALIDVALUES class variable.
378             The checks can be ignored if on turns of param checks like this:
379             $codeml->no_param_checks(1)
380             Returns : boolean if set was success, if verbose is set to -1
381             then no warning will be reported
382             Args : $paramname => name of the parameter
383             $value => value to set the parameter to
384             See also: L
385              
386             =cut
387              
388             sub set_parameter{
389 0     0 1   my ($self,$param,$value) = @_;
390 0 0         if( ! defined $VALIDVALUES{$param} ) {
391 0           $self->warn("unknown parameter $param will not set unless you force by setting no_param_checks to true");
392 0           return 0;
393             }
394 0 0 0       if( ref( $VALIDVALUES{$param}) =~ /ARRAY/i &&
395 0           scalar @{$VALIDVALUES{$param}} > 0 ) {
396            
397 0 0         unless ( grep {$value} @{ $VALIDVALUES{$param} } ) {
  0            
  0            
398 0           $self->warn("parameter $param specified value $value is not recognized, please see the documentation and the code for this module or set the no_param_checks to a true value");
399 0           return 0;
400             }
401             }
402 0           $self->{'_codemlparams'}->{$param} = $value;
403 0           return 1;
404             }
405              
406             =head2 set_default_parameters
407              
408             Title : set_default_parameters
409             Usage : $codeml->set_default_parameters(0);
410             Function: (Re)set the default parameters from the defaults
411             (the first value in each array in the
412             %VALIDVALUES class variable)
413             Returns : none
414             Args : boolean: keep existing parameter values
415              
416              
417             =cut
418              
419             sub set_default_parameters{
420 0     0 1   my ($self,$keepold) = @_;
421 0 0         $keepold = 0 unless defined $keepold;
422            
423 0           while( my ($param,$val) = each %VALIDVALUES ) {
424             # skip if we want to keep old values and it is already set
425 0 0 0       next if( defined $self->{'_codemlparams'}->{$param} && $keepold);
426 0 0         if(ref($val)=~/ARRAY/i ) {
427 0           $self->{'_codemlparams'}->{$param} = $val->[0];
428             } else {
429 0           $self->{'_codemlparams'}->{$param} = $val;
430             }
431             }
432             }
433              
434              
435             =head1 Bio::Tools::Run::Wrapper methods
436              
437             =cut
438              
439             =head2 no_param_checks
440              
441             Title : no_param_checks
442             Usage : $obj->no_param_checks($newval)
443             Function: Boolean flag as to whether or not we should
444             trust the sanity checks for parameter values
445             Returns : value of no_param_checks
446             Args : newvalue (optional)
447              
448              
449             =cut
450              
451             =head2 save_tempfiles
452              
453             Title : save_tempfiles
454             Usage : $obj->save_tempfiles($newval)
455             Function:
456             Returns : value of save_tempfiles
457             Args : newvalue (optional)
458              
459              
460             =cut
461              
462             =head2 outfile_name
463              
464             Title : outfile_name
465             Usage : my $outfile = $codeml->outfile_name();
466             Function: Get/Set the name of the output file for this run
467             (if you wanted to do something special)
468             Returns : string
469             Args : [optional] string to set value to
470              
471              
472             =cut
473              
474              
475             =head2 tempdir
476              
477             Title : tempdir
478             Usage : my $tmpdir = $self->tempdir();
479             Function: Retrieve a temporary directory name (which is created)
480             Returns : string which is the name of the temporary directory
481             Args : none
482              
483              
484             =cut
485              
486             =head2 cleanup
487              
488             Title : cleanup
489             Usage : $codeml->cleanup();
490             Function: Will cleanup the tempdir directory after a PAML run
491             Returns : none
492             Args : none
493              
494              
495             =cut
496              
497             =head2 io
498              
499             Title : io
500             Usage : $obj->io($newval)
501             Function: Gets a L object
502             Returns : L
503             Args : none
504              
505              
506             =cut
507              
508             1;