view release on metacpan or  search on metacpan

lib/Bio/Tools/Run/Alignment/Clustalw.pm  view on Meta::CPAN

    # 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);