File Coverage

blib/lib/Bio/Tools/Run/Alignment/Clustalw.pm

Criterion	Covered	Total	%
statement	38	225	16.8
branch	0	90	0.0
condition	0	35	0.0
subroutine	14	24	58.3
pod	10	10	100.0
total	62	384	16.1

package Bio::Tools::Run::Alignment::Clustalw;
$Bio::Tools::Run::Alignment::Clustalw::VERSION = '1.7.2';
use utf8;
use strict;
use warnings;

use Bio::Seq;
use Bio::SeqIO;
use Bio::SimpleAlign;
use Bio::AlignIO;
use Bio::TreeIO;
use Bio::Root::IO;

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

# ABSTRACT: Object for the calculation of a multiple sequence alignment from a set of unaligned sequences or alignments using the Clustalw program
# AUTHOR: Peter Schattner
# OWNER: Peter Schattner
# LICENSE: Perl_5

# AUTHOR: Jason Stajich
# AUTHOR: Sendu Bala


our @CLUSTALW_PARAMS = qw(output ktuple topdiags window pairgap fixedgap
                          floatgap matrix type transit dnamatrix outfile
                          gapopen gapext maxdiv gapdist hgapresidues pwmatrix
                          pwdnamatrix pwgapopen pwgapext score transweight
                          seed helixgap outorder strandgap loopgap terminalgap
                          helixendin helixendout strandendin strandendout program
                          reps outputtree seed bootlabels bootstrap);

our @CLUSTALW_SWITCHES = qw(help check options negative noweights endgaps
                            nopgap nohgap novgap kimura tossgaps
                            kimura tossgaps njtree);
our @OTHER_SWITCHES = qw(quiet);
our $PROGRAM_NAME = 'clustalw';
our $PROGRAM_DIR = $ENV{'CLUSTALDIR'} || $ENV{'CLUSTALWDIR'};


sub program_name {
  return $PROGRAM_NAME;
}


sub program_dir {
  return $PROGRAM_DIR;
}

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

    $self->_set_from_args(\@args, -methods => [@CLUSTALW_PARAMS,
                                               @CLUSTALW_SWITCHES,
                                               @OTHER_SWITCHES],
                                  -create => 1);

    return $self;
}


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

    return undef unless $self->executable;
    my $prog = $self->executable;
    my $string = `$prog -- 2>&1` ;
    $string =~ /\(?([\d.]+)\)?/xms;
    return $1 || undef;
}


sub run {
    my ($self,$input) = @_;
    my ($temp,$infilename, $seq);
    my ($attr, $value, $switch);

    $self->io->_io_cleanup();
    # Create input file pointer
    $infilename = $self->_setinput($input);
    $self->throw("Bad input data (sequences need an id) or less than 2 sequences in $input!") unless $infilename;

    # Create parameter string to pass to clustalw program
    my $param_string = $self->_setparams();

    # run clustalw
    return $self->_run('both', $infilename, $param_string);
}


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

    $self->io->_io_cleanup();

    # Create input file pointer
    my $infilename = $self->_setinput($input);
    $self->throw("Bad input data (sequences need an id ) or less than 2 sequences in $input !") unless $infilename;

    # Create parameter string to pass to clustalw program
    my $param_string = $self->_setparams();

    # run clustalw
    my $aln = $self->_run('align', $infilename, $param_string);
}



sub profile_align {
    my ($self,$input1,$input2) = @_;

    $self->io->_io_cleanup();

    # Create input file pointer
    my $infilename1 = $self->_setinput($input1, 1);
    my $infilename2 = $self->_setinput($input2, 2);
    if (!$infilename1 || !$infilename2) {$self->throw("Bad input data: $input1 or $input2 !");}
    unless ( -e $infilename1 and -e  $infilename2) {$self->throw("Bad input file: $input1 or $input2 !");}

    # Create parameter string to pass to clustalw program
    my $param_string = $self->_setparams();

    # run clustalw
    my $aln = $self->_run('profile-aln', $infilename1, $infilename2, $param_string);
}


sub add_sequences {

    my ($self,$input1,$input2) = @_;
    my ($temp,$infilename1,$infilename2,$input,$seq);

    $self->io->_io_cleanup();
    # Create input file pointer
    $infilename1 = $self->_setinput($input1,1);
    $infilename2 = $self->_setinput($input2,2);
    if (!$infilename1 || !$infilename2) {$self->throw("Bad input data: $input1 or $input2 !");}
    unless ( -e $infilename1 and -e  $infilename2) {$self->throw("Bad input file: $input1 or $input2 !");}


    # Create parameter string to pass to clustalw program
    my $param_string = $self->_setparams();
    # run clustalw
    my $aln = $self->_run('add_sequences', $infilename1,
    $infilename2, $param_string);

}


sub tree {
    my ($self,$input) = @_;

    $self->io->_io_cleanup();

    # Create input file pointer
    my $infilename = $self->_setinput($input);

    if (!$infilename) {$self->throw("Bad input data (sequences need an id ) or less than 2 sequences in $input !");}

    # Create parameter string to pass to clustalw program
    my $param_string = $self->_setparams();

    # run clustalw
    my $tree = $self->_run('tree', $infilename, $param_string);
}


sub footprint {
    my ($self, $in, $slice_size, $deviate) = @_;

    my ($simplealn, $tree) = $self->run($in);

    # total tree length?
    my $total_length = $tree->total_branch_length;

    # tree length along sliding window, picking regions that significantly
    # deviate from the average tree length
    $slice_size ||= 5;
    $deviate ||= 33;
    my $threshold = $total_length - (($total_length / 100) * $deviate);
    my $length = $simplealn->length;
    my $below = 0;
    my $found_minima = 0;
    my $minima = [$threshold, ''];
    my @results;
    for my $i (1..($length - $slice_size + 1)) {
        my $slice = $simplealn->slice($i, ($i + $slice_size - 1), 1);
        my $tree = $self->tree($slice);
        $self->throw("No tree returned") unless defined $tree;
        my $slice_length = $tree->total_branch_length;

        $slice_length <= $threshold ? ($below = 1) : ($below = 0);
        if ($below) {
            unless ($found_minima) {
                if ($slice_length < ${$minima}[0]) {
                    $minima = [$slice_length, $slice];
                }
                else {
                    push(@results, ${$minima}[1]);
                    $minima = [$threshold, ''];
                    $found_minima = 1;
                }
            }
        }
        else {
            $found_minima = 0;
        }
    }

    return @results;
}


sub _run {
    my ($self, $command, $infile1, $infile2, $param_string) = @_;

    my ($instring, $tree);
    my $quiet = $self->quiet() || $self->verbose() < 0;

    if ($command =~ /align|both/) {
        if ($^O eq 'dec_osf') {
            $instring = $infile1;
            $command = '';
        }
        else {
            $instring = " -infile=". '"' . $infile1 . '"';
        }
        $param_string .= " $infile2";
    }

    if ($command =~ /profile/) {
        $instring =  "-profile1=$infile1  -profile2=$infile2";
        chmod 0777, $infile1, $infile2;
        $command = '-profile';
    }

    if ($command =~ /add_sequences/) {
        $instring =  "-profile1=$infile1  -profile2=$infile2";
        chmod 0777, $infile1,$infile2;
        $command = '-sequences';
    }

    if ($command =~ /tree/) {
        if( $^O eq 'dec_osf' ) {
            $instring =  $infile1;
            $command = '';
        }
        else {
            $instring = " -infile=". '"' . $infile1 . '"';
        }
        $param_string .= " $infile2";

        $self->debug( "Program ".$self->executable."\n");
        my $commandstring = $self->executable."$instring"."$param_string";
        my $null = ($^O =~ m/mswin/i) ? 'NUL' : '/dev/null';
        $commandstring .= " 1>$null" if $quiet;
        $self->debug( "clustal command = $commandstring");

        my $status = system($commandstring);
        unless( $status == 0 ) {
            $self->warn( "Clustalw call ($commandstring) crashed: $? \n");
            return undef;
        }

        return $self->_get_tree($infile1, $param_string);
    }

    my $output = $self->output || 'gcg';
    $self->debug( "Program ".$self->executable."\n");
    my $commandstring = $self->executable." $command"." $instring"." -output=$output". " $param_string";
    $self->debug( "clustal command = $commandstring\n");

    open(my $pipe, "$commandstring |") || $self->throw("ClustalW call ($commandstring) failed to start: $? | $!");
    my $score;
    while (<$pipe>) {
        print unless $quiet;
        # Kevin Brown suggested the following regex, though it matches multiple
        # times: we pick up the last one
        $score = $1 if ($_ =~ /Score:(\d+)/);
        # This one is printed at the end and seems the most appropriate to pick
        # up; we include the above regex incase 'Alignment Score' isn't given
        $score = $1 if ($_ =~ /Alignment Score (-?\d+)/);
    }
    close($pipe) || ($self->throw("ClustalW call ($commandstring) crashed: $?"));

    my $outfile = $self->outfile();

    # retrieve alignment (Note: MSF format for AlignIO = GCG format of clustalw)
    my $format = $output =~ /gcg/i ? 'msf' : $output;
    if ($format =~ /clustal/i) {
        $format = 'clustalw'; # force clustalw incase 'clustal' is requested
    }
    my $in  = Bio::AlignIO->new(-file  => $outfile, -format=> $format);
    my $aln = $in->next_aln();
    $in->close;
    $aln->score($score);

    if ($command eq 'both') {
        $tree = $self->_get_tree($infile1, $param_string);
    }

    # Clean up the temporary files created along the way...
    # Replace file suffix with dnd to find name of dendrogram file(s) to delete
    unless ( $self->save_tempfiles ) {
        foreach my $f ($infile1, $infile2) {
            $f =~ s/\.[^\.]*$// ;
            unlink $f .'.dnd' if ($f ne '');
        }
    }

    if ($command eq 'both') {
        return ($aln, $tree);
    }
    return $aln;
}

sub _get_tree {
    my ($self, $treefile, $param_string) = @_;

    $treefile =~ s/\.[^\.]*$// ;

    if ($param_string =~ /-bootstrap/) {
        $treefile .= '.phb';
    }
    elsif ($param_string =~ /-tree/) {
        $treefile .= '.ph';
    }
    else {
        $treefile .= '.dnd';
    }

    my $in = Bio::TreeIO->new('-file'  => $treefile,
                             '-format'=> 'newick');

    my $tree = $in->next_tree;
    unless ( $self->save_tempfiles ) {
        foreach my $f ( $treefile ) {
            unlink $f if( $f ne '' );
        }
    }

    return $tree;
}


sub _setinput {
    my ($self, $input, $suffix) = @_;
    my ($infilename, $seq, $temp, $tfh);

    # suffix is used to distinguish alignment files If $input is not a
    # reference it better be the name of a file with the sequence/

    # alignment data...
    unless (ref $input) {
        # check that file exists or throw
        $infilename = $input;
        return unless -e $input;
        return $infilename;
    }

    # $input may be an array of BioSeq objects...
    if (ref($input) eq "ARRAY") {
        #  Open temporary file for both reading & writing of BioSeq array
        ($tfh,$infilename) = $self->io->tempfile(-dir=>$self->tempdir);
        $temp = Bio::SeqIO->new('-fh'=>$tfh, '-format' =>'Fasta');

        # Need at least 2 seqs for alignment
        return unless (scalar(@$input) > 1);

        foreach $seq (@$input) {
            return unless (defined $seq && $seq->isa("Bio::PrimarySeqI") and $seq->id());
            $temp->write_seq($seq);
        }
        $temp->close();
        close($tfh);
        undef $tfh;
        return $infilename;
    }

    # $input may be a SimpleAlign object.
    elsif (ref($input) eq "Bio::SimpleAlign") {
        # Open temporary file for both reading & writing of SimpleAlign object
        ($tfh,$infilename) = $self->io->tempfile(-dir=>$self->tempdir);
        $temp = Bio::AlignIO->new('-fh'=> $tfh, '-format' => 'fasta');
        $temp->write_aln($input);
        close($tfh);
        undef $tfh;
        return $infilename;
    }

    # or $input may be a single BioSeq object (to be added to a previous alignment)
    elsif (ref($input) && $input->isa("Bio::PrimarySeqI") && $suffix==2) {
        # Open temporary file for both reading & writing of BioSeq object
        ($tfh,$infilename) = $self->io->tempfile();
        $temp = Bio::SeqIO->new(-fh=> $tfh, '-format' =>'Fasta');
        $temp->write_seq($input);
        close($tfh);
        undef $tfh;
        return $infilename;
    }

    return;
}


sub _setparams {
    my $self = shift;

    my $param_string = $self->SUPER::_setparams(-params => \@CLUSTALW_PARAMS,
                                                -switches => \@CLUSTALW_SWITCHES,
                                                -dash => 1,
                                                -lc => 1,
                                                -join => '=');

    # Set default output file if no explicit output file selected
    unless ($param_string =~ /outfile/) {
        my ($tfh, $outfile) = $self->io->tempfile(-dir => $self->tempdir());
        close($tfh);
        undef $tfh;
        $self->outfile($outfile);
        $param_string .= " -outfile=\"$outfile\"" ;
    }

    $param_string .= ' 2>&1';

    return $param_string;
}

1;

__END__