File Coverage

blib/lib/Bio/Tools/Run/Phylo/Raxml.pm

Criterion	Covered	Total	%
statement	42	116	36.2
branch	2	36	5.5
condition	1	17	5.8
subroutine	14	21	66.6
pod	7	7	100.0
total	66	197	33.5

#
# BioPerl module for Bio::Tools::Run::Phylo::Raxml
#
# Please direct questions and support issues to
#
# Copyright Brian Osborne
#
# You may distribute this module under the same terms as perl itself
#
# POD documentation - main docs before the code

=head1 NAME

Bio::Tools::Run::Phylo::Raxml

=head1 SYNOPSIS

  # Build a Raxml factory
  $factory = Bio::Tools::Run::Phylo::Raxml->new(-p  => 100);

  # Get an alignment
  my $alignio = Bio::AlignIO->new(
        -format => 'fasta',
        -file   => '219877.cdna.fasta');
  my $alnobj = $alignio->next_aln;

  # Analyze the aligment and get a Tree
  my $tree = $factory->run($alnobj);

=head1 DESCRIPTION

Get a Bio::Tree object using raxml given a protein or DNA alignment.

=head1 FEEDBACK

=head2 Mailing Lists

User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to one
of the Bioperl mailing lists.  Your participation is much appreciated.

  [e-mail]                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

=head2 Support

Please direct usage questions or support issues to the mailing list:

I

Do not contact the module maintainer directly. Many experienced experts
at bioperl-l will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution.  Bug reports can be submitted via the web:

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

=head1 AUTHOR -  Brian Osborne

Email [e-mail]

=head1 APPENDIX

The rest of the documentation details each of the object
methods. Internal methods are usually preceded with a _

=cut

package Bio::Tools::Run::Phylo::Raxml;

use strict;
use File::Basename;
use File::Spec;
use Bio::Seq;
use Bio::SeqIO;
use Bio::TreeIO;
use Bio::AlignIO;
use Bio::Root::IO;
use Cwd;

use base qw(Bio::Root::Root Bio::Tools::Run::WrapperBase);

our @Raxml_PARAMS =
  qw(s n m a A b B c e E f g G i I J o p P q r R S t T w W x z N);
our @Raxml_SWITCHES =
  qw(SSE3 PTHREADS PTHREADS-SSE3 HYBRID HYBRID-SSE3 F h k K M j U v X y C d D);
our $PROGRAM_NAME = 'raxml';

# Specify some model if none is specified
my $DEFAULTAAMODEL = 'PROTCATDAYHOFF';
my $DEFAULTNTMODEL = 'GTRCAT';

=head2 new

 Title   : new
 Usage   : my $treebuilder = Bio::Tools::Run::Phylo::Raxml->new();
 Function: Constructor
 Returns : Bio::Tools::Run::Phylo::Raxml
 Args    : Same as those used to run raxml. For example:

 $factory = Bio::Tools::Run::Phylo::Raxml->new(-p  => 100, -SSE3 => 1)

=cut

sub new {
    my ( $class, @args ) = @_;
    my $self = $class->SUPER::new(@args);

    $self->_set_from_args(
        \@args,
        -case_sensitive => 1,
        -methods => [ @Raxml_PARAMS, @Raxml_SWITCHES ],
        -create  => 1
    );

    my ($out,$quiet) = $self->SUPER::_rearrange( [qw(OUTFILE_NAME QUIET)], @args );

    $self->outfile_name( $out || '' );
    $self->quiet( $quiet ) if $quiet;

    $self;
}

=head2 program_name

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

=cut

sub program_name {
    $PROGRAM_NAME;
}

=head2 program_dir

 Title   : program_dir
 Usage   : $factory->program_dir(@params)
 Function: returns the program directory
 Returns:  string
 Args    :

=cut

sub program_dir {
    undef;
}

=head2 error_string

 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)

=cut

sub error_string {
    my ( $self, $value ) = @_;

    $self->{'error_string'} = $value if ( defined $value );
    $self->{'error_string'};
}

=head2  version

 Title   : version
 Usage   : exit if $prog->version() < 1.8
 Function: Determine the version number of the program
 Example :
 Returns : float or undef
 Args    : none

=cut

sub version {
    my ($self) = @_;
    my $exe;

    return undef unless $exe = $self->executable;
    my $string = `$exe -v 2>&1`;

    $string =~ /raxml\s+version\s+([\d\.]+)/i;
    return $1 || undef;
}

=head2  quiet

 Title   : quiet
 Usage   :
 Function: get or set value for 'quiet'
 Example :
 Returns :
 Args    : the value

=cut

sub quiet {
    my ( $self, $value ) = @_;

    $self->{'_quiet'} = $value if ( defined $value );
    $self->{'_quiet'};
}

=head2 run

 Title   : run
 Usage   : $factory->run($stockholm_file) OR
           $factory->run($align_object)
 Function: Runs Raxml to generate a tree
 Returns : Bio::Tree::Tree object
 Args    : File name for your input alignment in stockholm format, OR
           Bio::Align::AlignI compliant object (eg. Bio::SimpleAlign).

=cut

sub run {
    my ($self, $in) = @_;

    if (ref $in && $in->isa("Bio::Align::AlignI")) {
        $in = $self->_write_alignfile($in);
    }
    elsif (! -e $in) {
        $self->throw("When not supplying a Bio::Align::AlignI object, you must supply a readable filename");
    }

    $self->_run($in);
}

=head2 _run

 Title   : _run
 Usage   : Internal function, not to be called directly
 Function: Runs the application
 Returns : Tree object
 Args    : Alignment file name

=cut

sub _run {
    my ( $self, $file ) = @_;

    my $exe       = $self->executable || return;
    my $param_str = $self->arguments . " " . $self->_setparams($file);
    my $command   = "$exe $param_str";
    $self->debug("Raxml command = $command");

    my $status  = system($command);

    # raxml creates tree files with names like "RAxML_bestTree.ABDBxjjdfg3"
    # if rapid bootstrapping was enabled, also a tree with RAxML_bipartitions.ABDBxjjdfg3
    # with support values is created, which then should be returned
    my $outfile = $self->f() eq 'a' ? 'RAxML_bipartitions.' : 'RAxML_bestTree.';
    $outfile .= $self->outfile_name;

    $outfile = File::Spec->catfile( ($self->w), $outfile ) if $self->w;

    if ( !-e $outfile || -z $outfile ) {
        $self->warn("Raxml call had status of $status: $? [command $command] \n");
        return undef;
    }

    my $treeio = Bio::TreeIO->new( -file => $outfile );
    my $tree = $treeio->next_tree;

# if bootstraps were enabled, the bootstraps are the ids; convert to
# bootstrap and no id
# if ($self->boot) {
#     my @nodes = $tree->get_nodes;
#     my %non_internal = map { $_ => 1 } ($tree->get_leaf_nodes, $tree->get_root_node);
#     foreach my $node (@nodes) {
#         next if exists $non_internal{$node};
#         $node->bootstrap && next; # protect ourselves incase the parser improves
#         $node->bootstrap($node->id);
#         $node->id('');
#     }
# }

    $tree;
}

=head2 _write_alignfile

 Title   : _write_alignfile
 Usage   : Internal function, not to be called directly
 Function: Create an alignment file
 Returns : filename
 Args    : Bio::Align::AlignI

=cut

sub _write_alignfile {
    my ( $self, $align ) = @_;

    my ( $tfh, $tempfile ) = $self->io->tempfile( -dir => '.' );

    my $out = Bio::AlignIO->new(
        -file   => ">$tempfile",
        -format => 'phylip'
    );
    $out->write_aln($align);
    $out->close();
    undef($out);
    close($tfh);
    undef($tfh);

    die "Alignment file $tempfile was not created" if ( ! -e $tempfile );

    $tempfile;
}

=head2 _alphabet

 Title   : _alphabet
 Usage   : my $alphabet = $self->_alphabet;
 Function: Get the alphabet of the input alignment, defaults to 'dna'
 Returns : 'dna' or 'protein'
 Args    : Alignment file

=cut

sub _alphabet {
    my ( $self, $file ) = @_;

    if ($file) {
        if ( -e $file ) {
            my $in = Bio::AlignIO->new( -file => $file );
            my $aln = $in->next_aln;

            # arbitrary, the first one
            my $seq      = $aln->get_seq_by_pos(1);
            my $alphabet = $seq->alphabet;
            $self->{_alphabet} = $alphabet;
        }
        else {
            die "File $file can not be found";
        }
    }

    # default is 'dna'
    return $self->{'_alphabet'} || 'dna';
}

=head2  _setparams

 Title   :  _setparams
 Usage   :  Internal function, not to be called directly
 Function:  Create parameter inputs for Raxml program
 Example :
 Returns : parameter string to be passed to Raxml
 Args    : name of calling object

=cut

sub _setparams {
    my ( $self, $infile ) = @_;
    my $param_string = '';

    # If 'model' is not set with '-m' check the alphabet of the input,
    # then specify the default model
    if ( !$self->m ) {
        my $model =
          ( $self->_alphabet($infile) eq 'dna' )
          ? $DEFAULTNTMODEL
          : $DEFAULTAAMODEL;
        $self->m($model);
    }

    # Set default output file if no explicit output file has been given.
    # Raxml insists that the output file name not contain '/' and its
    # output directory is set using the '-w' argument.
    if ( !$self->outfile_name ) {
        my $dir = getcwd();
        $self->w($dir);

        my ( $tfh, $outfile ) = $self->io->tempfile( -dir => $dir );
        close($tfh);
        undef $tfh;
        $outfile = basename($outfile);
        $self->outfile_name($outfile);
    }

    for my $attr (@Raxml_PARAMS) {
        my $value = $self->$attr();
        next unless ( defined $value );
        $param_string .= ' -' . $attr . ' ' . $value . ' ';
    }

    for my $attr (@Raxml_SWITCHES) {
        my $value = $self->$attr();
        next unless ($value);
        $param_string .= ' -' . $attr . ' ';
    }

    $param_string .= "-s $infile -n " . $self->outfile_name;

    my $null = File::Spec->devnull();
    $param_string .= " > $null 2> $null"
      if ( $self->quiet() || $self->verbose < 0 );

    $param_string;
}

=head1 Bio::Tools::Run::BaseWrapper methods

=cut

=head2 no_param_checks

Title   : no_param_checks
Usage   : $obj->no_param_checks($newval)
Function: Boolean flag as to whether or not we should
          trust the sanity checks for parameter values
Returns : value of no_param_checks
Args    : newvalue (optional)

=cut

=head2 save_tempfiles

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

=cut

=head2 outfile_name

Title   : outfile_name
Usage   : my $outfile = $Raxml->outfile_name();
Function: Get/Set the name of the output file for this run
           (if you wanted to do something special)
Returns : string
Args    : [optional] string to set value to

=cut

=head2 tempdir

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

=cut

=head2 cleanup

Title   : cleanup
Usage   : $Raxml->cleanup();
Function: Will cleanup the tempdir directory
Returns : none
Args    : none

=cut

=head2 io

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

=cut

1;    # Needed to keep compiler happy

__END__