use Bio::Tools::Run::Phylo::Hyphy::SLAC;
  use Bio::AlignIO;
  use Bio::TreeIO;

  my $alignio = Bio::AlignIO->new(-format => 'fasta',
  			         -file   => 't/data/hyphy1.fasta');

  my $aln = $alignio->next_aln;

  my $treeio = Bio::TreeIO->new(
      -format => 'newick', -file => 't/data/hyphy1.tree');

  my $slac = Bio::Tools::Run::Phylo::Hyphy::SLAC->new();
  $slac->alignment($aln);
  $slac->tree($tree);
  my ($rc,$results) = $slac->run();

     @VALIDVALUES = 
        (
         {'geneticCode' => [ "Universal","VertebratemtDNA","YeastmtDNA","Mold/ProtozoanmtDNA",
                             "InvertebratemtDNA","CiliateNuclear","EchinodermmtDNA","EuplotidNuclear",
                             "Alt.YeastNuclear","AscidianmtDNA","FlatwormmtDNA","BlepharismaNuclear"]},
         {'New/Restore' => [ "New Analysis", "Restore"]},
         {'tempalnfile' => undef }, # aln file goes here
         {'Model Options' => [ { "Custom" => '010010' },
                               { "Default" => undef } ] 
         },
         {'temptreefile' => undef }, # tree file goes here
         {'Model Fit Results' => [ '/dev/null'] }, # this will not work under Windows
         {'dN/dS bias parameter' => [ { "Estimate dN/dS only" => undef },
                                      { "Neutral" => undef },
                                      { "Estimate" => undef },
                                      { "Estimate + CI" => undef },
                                      { "User" => '3' } ] },
         {'Ancestor Counting' => [ 'Single Ancestor Counting','Weighted Ancestor Counting',
                                  'Sample Ancestal States','Process Sampled Ancestal States',
                                  'One rate FEL','Two rate FEL','Rate Distribution',
                                  'Full site-by-site LRT','Multirate FEL'] },
         {'SLAC Options' => ['Full tree','Tips vs Internals'] },
         {'Treatment of Ambiguities' => ['Resolved','Averaged'] },
         {'Test Statistic' => ['Approximate','Simulated Null'] },
         {'Significance level' => '0.05' },
         {'Output options' => 'Export to File' }, #we force a tsv file here
         {'outfile' => undef }, # outfile goes here
         {'Rate class estimator' => [ 'Skip','Count'] },
        );

  my($class,@args) = @_;

  my $self = $class->SUPER::new(@args);
  my ($aln, $tree, $st, $params, $exe, 
      $ubl) = $self->_rearrange([qw(ALIGNMENT TREE SAVE_TEMPFILES 
				    PARAMS EXECUTABLE)],
				    @args);
  defined $aln && $self->alignment($aln);
  defined $tree && $self->tree($tree);
  defined $st  && $self->save_tempfiles($st);
  defined $exe && $self->executable($exe);
  $self->set_default_parameters();
  if( defined $params ) {
      if( ref($params) !~ /HASH/i ) { 
	  $self->warn("Must provide a valid hash ref for parameter -FLAGS");
      } else {
	  map { $self->set_parameter($_, $$params{$_}) } keys %$params;
      }
  }

  return $self;

   my ($self,$aln,$tree) = @_;

   $self->prepare($aln,$tree) unless (defined($self->{'_prepared'}));
   my ($rc,$results) = (1);
   {
       my $commandstring;
       my $exit_status;
       my $tempdir = $self->tempdir;
       my $slacexe = $self->executable();
       $self->throw("unable to find or run executable for 'HYPHY'") unless $slacexe && -e $slacexe && -x _;
       $commandstring = $slacexe . " BASEPATH=" . $self->program_dir . " " . $self->{'_wrapper'};
       open(RUN, "$commandstring |") or $self->throw("Cannot open exe $slacexe");
       my @output = <RUN>;
       $exit_status = close(RUN);
       $self->error_string(join('',@output));
       if( (grep { /\berr(or)?: /io } @output)  || !$exit_status) {
	   $self->warn("There was an error - see error_string for the program output");
	   $rc = 0;
       }
       my $outfile = $self->outfile_name;
       eval {
	   open(OUTFILE, "$outfile") or $self->throw("cannot open $outfile for reading");
           my $readed_header = 0;
           my @elems;
           while (<OUTFILE>) {
               if ($readed_header) {
                   # FEL results are tsv
                   my @values = split("\t",$_);
                   for my $i (0 .. (scalar(@values)-1)) {
                       $elems[$i] =~ s/\n//g;
                       push @{$results->{$elems[$i]}}, $values[$i];
                   }
               } else {
                   @elems = split("\t",$_);
                   $readed_header = 1;
               }
           }
       };
       if( $@ ) {
	   $self->warn($self->error_string);
       }
   }
   unless ( $self->save_tempfiles ) {
       unlink($self->{'_wrapper'});
      $self->cleanup();
   }
   return ($rc,$results);

   my $self = shift;

   my $batchfile = 'QuickSelectionDetection.bf';
   $self->SUPER::create_wrapper($batchfile);

   my ($self,$keepold) = @_;
   $keepold = 0 unless defined $keepold;
   foreach my $elem (@VALIDVALUES) {
       my ($param,$val) = each %$elem;
       # skip if we want to keep old values and it is already set
       if (ref($val)=~/ARRAY/i ) {
           unless (ref($val->[0])=~/HASH/i) {
               push @{ $self->{'_orderedparams'} }, {$param, $val->[0]};
           } else {
               $val = $val->[0];
           }
       } 
       if ( ref($val) =~ /HASH/i ) { 
           my $prevparam;
           while (defined($val)) {
               last unless (ref($val) =~ /HASH/i);
               last unless (defined($param));
               $prevparam = $param;
               ($param,$val) = each %{$val};
               push @{ $self->{'_orderedparams'} }, {$prevparam, $param};
               push @{ $self->{'_orderedparams'} }, {$param, $val} if (defined($val));
           }
       } elsif (ref($val) !~ /HASH/i && ref($val) !~ /ARRAY/i) { 
           push @{ $self->{'_orderedparams'} }, {$param, $val};
       }
   }

    my $self= shift;
    unless ( $self->save_tempfiles ) {
	$self->cleanup();
    }
    $self->SUPER::DESTROY();
}

1;

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  http://redmine.open-bio.org/projects/bioperl/

 Title   : program_name
 Usage   : $factory->program_name()
 Function: holds the program name
 Returns:  string
 Args    : None

 Title   : program_dir
 Usage   : ->program_dir()
 Function: returns the program directory, obtained from ENV variable.
 Returns:  string
 Args    :

 Title   : prepare
 Usage   : my $rundir = $slac->prepare($aln);
 Function: prepare the slac analysis using the default or updated parameters
           the alignment parameter must have been set
 Returns : value of rundir
 Args    : Bio::Align::AlignI object,
	   Bio::Tree::TreeI object [optional]

 Title   : error_string
 Usage   : $obj->error_string($newval)
 Function: Where the output from the last analysus run is stored.
 Returns : value of error_string
 Args    : newvalue (optional)

 Title   : alignment
 Usage   : $slac->align($aln);
 Function: Get/Set the Bio::Align::AlignI object
 Returns : Bio::Align::AlignI object
 Args    : [optional] Bio::Align::AlignI
 Comment : We could potentially add support for running directly on a file
           but we shall keep it simple
 See also: Bio::SimpleAlign

 Title   : tree
 Usage   : $slac->tree($tree, %params);
 Function: Get/Set the Bio::Tree::TreeI object
 Returns : Bio::Tree::TreeI 
 Args    : [optional] $tree => Bio::Tree::TreeI,
           [optional] %parameters => hash of tree-specific parameters:

 Comment : We could potentially add support for running directly on a file
           but we shall keep it simple
 See also: Bio::Tree::Tree

 Title   : get_parameters
 Usage   : my %params = $self->get_parameters();
 Function: returns the list of parameters as a hash
 Returns : associative array keyed on parameter names
 Args    : none

 Title   : set_parameter
 Usage   : $slac->set_parameter($param,$val);
 Function: Sets a slac parameter, will be validated against
           the valid values as set in the %VALIDVALUES class variable.  
           The checks can be ignored if one turns off param checks like this:
             $slac->no_param_checks(1)
 Returns : boolean if set was success, if verbose is set to -1
           then no warning will be reported
 Args    : $param => name of the parameter
           $value => value to set the parameter to
 See also: no_param_checks()

 Title   : update_ordered_parameters
 Usage   : $slac->update_ordered_parameters(0);
 Function: (Re)set the default parameters from the defaults
           (the first value in each array in the 
	    %VALIDVALUES class variable)
 Returns : none
 Args    : boolean: keep existing parameter values

 Title   : save_tempfiles
 Usage   : $obj->save_tempfiles($newval)
 Function: 
 Returns : value of save_tempfiles
 Args    : newvalue (optional)

 Title   : tempdir
 Usage   : my $tmpdir = $self->tempdir();
 Function: Retrieve a temporary directory name (which is created)
 Returns : string which is the name of the temporary directory
 Args    : none

 Title   : cleanup
 Usage   : $slac->cleanup();
 Function: Will cleanup the tempdir directory after a run
 Returns : none
 Args    : none

 Title   : io
 Usage   : $obj->io($newval)
 Function:  Gets a Bio::Root::IO object
 Returns : Bio::Root::IO
 Args    : none