File Coverage

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

Criterion	Covered	Total	%
statement	71	258	27.5
branch	12	102	11.7
condition	0	41	0.0
subroutine	19	27	70.3
pod	10	10	100.0
total	112	438	25.5

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

use vars qw($AUTOLOAD @ISA $PROGRAM_NAME $PROGRAM_DIR %DEFAULTS
            @TCOFFEE_PARAMS @TCOFFEE_SWITCHES @OTHER_SWITCHES %OK_FIELD
            );

use Cwd;
use Bio::Seq;
use Bio::SeqIO;
use Bio::SimpleAlign;
use Bio::AlignIO;
use Bio::Root::IO;
use Bio::Factory::ApplicationFactoryI;
use Bio::Tools::Run::WrapperBase;

@ISA = qw(Bio::Tools::Run::WrapperBase
          Bio::Factory::ApplicationFactoryI);

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


# You will need to enable TCoffee to find the tcoffee program. This can be done
# in (at least) twp ways:
#  1. define an environmental variable TCOFFEE:
#        export TCOFFEEDIR=/home/progs/tcoffee   or
#  2. include a definition of an environmental variable TCOFFEEDIR
#      in every script that will
#     use Bio::Tools::Run::Alignment::TCoffee.pm.
#        BEGIN {$ENV{TCOFFEEDIR} = '/home/progs/tcoffee'; }

BEGIN {
    $PROGRAM_NAME = 't_coffee';
    $PROGRAM_DIR = $ENV{'TCOFFEEDIR'};
    %DEFAULTS = ( 'MATRIX' => 'blosum',
                  'OUTPUT' => 'clustalw',
                  'AFORMAT'=> 'msf',
                  'METHODS' => [qw(lalign_id_pair clustalw_pair)]
                );


    @TCOFFEE_PARAMS = qw(IN TYPE PARAMETERS DO_NORMALISE EXTEND
                         DP_MODE KTUPLE NDIAGS DIAG_MODE SIM_MATRIX
                         MATRIX GAPOPEN GAPEXT COSMETIC_PENALTY TG_MODE
                         WEIGHT SEQ_TO_ALIGN NEWTREE USETREE TREE_MODE
                         OUTFILE OUTPUT CASE CPU OUT_LIB OUTORDER SEQNOS
                         RUN_NAME CONVERT
                         );

    @TCOFFEE_SWITCHES = qw(QUICKTREE);

    @OTHER_SWITCHES = qw(QUIET ALIGN KEEPDND);

# Authorize attribute fields
    foreach my $attr ( @TCOFFEE_PARAMS, @TCOFFEE_SWITCHES, @OTHER_SWITCHES ) {
        $OK_FIELD{$attr}++; }
}


sub program_name {
    return $PROGRAM_NAME;
}


sub program_dir {
    return $PROGRAM_DIR;
}

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

    while (@args) {
        $attr =   shift @args;
        $value =  shift @args;
        next if( $attr =~ /^-/); # don't want named parameters
        $self->$attr($value);
    }
    $self->matrix($DEFAULTS{'MATRIX'}) unless( $self->matrix );
    $self->output($DEFAULTS{'OUTPUT'}) unless( $self->output );
    $self->methods($DEFAULTS{'METHODS'}) unless $self->methods;
    # $self->aformat($DEFAULTS{'AFORMAT'}) unless $self->aformat;

    return $self;
}

sub AUTOLOAD {
    my $self = shift;
    my $attr = $AUTOLOAD;
    $attr =~ s/.*:://;
    $attr = uc $attr;
    # aliasing
    $attr = 'OUTFILE' if $attr eq 'OUTFILE_NAME';
    $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};

    $self->{$attr} = shift if @_;
    return $self->{$attr};
}


sub error_string{
   my ($self,$value) = @_;
   if( defined $value) {
      $self->{'error_string'} = $value;
    }
    return $self->{'error_string'};

}


sub version {
    my ($self) = @_;
    my $exe;
    return undef unless $exe = $self->executable;
    my $string = `$exe -quiet=stdout 2>&1` ;
    $string =~ /Version_([\d.]+)/;
    return $1 || undef;

}


sub run{
   my ($self,@args) = @_;
   my ($type,$seq,$profile) = $self->_rearrange([qw(TYPE
                                                    SEQ
                                                    PROFILE)],
                                                @args);
   if( $type =~ /align/i ) {
       return $self->align($seq);
   } elsif( $type =~ /profile/i ) {
       return $self->profile_align($profile,$seq);
   } else {
       $self->warn("unrecognized alignment type $type\n");
   }
   return undef;
}


sub align {
    my ($self,$input) = @_;
    # Create input file pointer
    $self->io->_io_cleanup();
    my ($infilename,$type) = $self->_setinput($input);
    if (!$infilename) {
        $self->throw("Bad input data or less than 2 sequences in $input !");
    }

    my $param_string = $self->_setparams();

    # run tcoffee
    return $self->_run('align', [$infilename,$type], $param_string);
}

#################################################


sub profile_align {
    my $self = shift;
    my $input1 = shift;
    my $input2 = shift;

    my ($temp,$infilename1,$infilename2,$type1,$type2,$input,$seq);

    $self->io->_io_cleanup();
    # Create input file pointers
    ($infilename1,$type1) = $self->_setinput($input1);
    ($infilename2,$type2) = $self->_setinput($input2);
    unless ($type1) {
        $self->throw("Unknown type for first argument");
    }
    unless ($type2) {
        $self->throw("Unknown type for second argument")
    }

    if (!$infilename1 || !$infilename2) {
     $self->throw("Bad input data: $input1 or $input2 !");
    }

    my $param_string = $self->_setparams();
    # run tcoffee
    my $aln = $self->_run('profile-aln',
                          [$infilename1,$type1],
                          [$infilename2,$type2],
                          $param_string)
            ;
}
#################################################


sub _run {
    my ($infilename, $infile1,$infile2) = ('','','');
    my $self = shift;
    my $command = shift;
    my $instring;
    if ($command =~ /align/) {
        my $infile = shift ;
        my $type;
        ($infilename,$type) = @$infile;
        $instring =  '-in='.join(',',($infilename, 'X'.$self->matrix,
                          $self->methods));
    }
    if ($command =~ /profile/) {
        my $in1 = shift ;
        my $in2 = shift ;
        my ($type1,$type2);
        ($infile1,$type1) = @$in1;
        ($infile2,$type2) = @$in2;
        # in later versions (tested on 5.72 and 7.54) the API for profile
        # alignment changed. This attempts to do the right thing for older
        # versions but corrects for newer ones
        if ($self->version && $self->version < 5) {
            # this breaks severely on newer TCoffee (>= v5)
            unless (($self->matrix =~ /none/i) || ($self->matrix =~ /null/i) ) {
                $instring = '-in='.join(',',
                                ($type2.$infile2),
                                'X'.$self->matrix,
                                (map {'M'.$_} $self->methods)
                                );
                $instring .= ' -profile='.$infile1;
            } else {
                $instring = '-in='.join(',',(
                            $type1.$infile1,
                            $type2.$infile2,
                            (map {'M'.$_} $self->methods)
                                            )
                                        );
            }
        } else {
            if ($type2 eq 'S') {
                # second infile is a sequence, not an alignment
                $instring .= ' -profile='.join(',',$infile1);
                $instring .= ' -seq='.join(',',$infile2);
            } else {
                $instring .= ' -profile='.join(',',$infile1,$infile2);
            }
            $instring .= ' -matrix='.$self->matrix unless (($self->matrix =~ /none/i) || ($self->matrix =~ /null/i)) ;
            $instring .= ' -method='.join(',',$self->methods) if ($self->methods) ;
        }
    }
    my $param_string = shift;
#    my ($paramfh,$parameterFile) = $self->io->tempfile;
#    print $paramfh ( "$instring\n-output=gcg$param_string") ;
#    close($paramfh);
#    my $commandstring = "t_coffee -output=gcg -parameters $parameterFile" ;    ##MJL

    my $commandstring = $self->executable." $instring $param_string";

    #$self->debug( "tcoffee command = $commandstring \n");

    my $status = system($commandstring);
    my $outfile = $self->outfile();

    if( !-e $outfile || -z $outfile ) {
        $self->warn( "TCoffee call crashed: $? [command $commandstring]\n");
        return undef;
    }

    # retrieve alignment (Note: MSF format for AlignIO = GCG format of
    # tcoffee)

    my $in  = Bio::AlignIO->new('-file' => $outfile, '-format' =>
                                $self->output);
    my $aln = $in->next_aln();

    # Replace file suffix with dnd to find name of dendrogram file(s) to delete
    if( ! $self->keepdnd ) {
        foreach my $f ( $infilename, $infile1, $infile2 ) {
            next if( !defined $f || $f eq '');
            $f =~ s/\.[^\.]*$// ;
            # because TCoffee writes these files to the CWD
            if( $Bio::Root::IO::PATHSEP ) {
                my @line = split(/$Bio::Root::IO::PATHSEP/, $f);
                $f = pop @line;
            } else {
                (undef, undef, $f) = File::Spec->splitpath($f);
            }
            unlink $f .'.dnd' if( $f ne '' );
        }
    }
    return $aln;
}



sub _setinput {
    my ($self,$input) = @_;
    my ($infilename, $seq, $temp, $tfh);
    #  If $input is not a reference it better be the name of a
    # file with the sequence/ alignment data...
    my $type = '';
    if (! ref $input) {
        # check that file exists or throw
        $infilename = $input;
        unless (-e $input) {return 0;}
        # let's peek and guess
        open(my $IN,$infilename) || $self->throw("Cannot open $infilename");
        my $header = <$IN>;
        if( $header =~ /^\s+\d+\s+\d+/ ||
            $header =~ /Pileup/i ||
            $header =~ /clustal/i ) { # phylip
            $type = 'A';
        }

        # On some systems, having filenames with / in them (ie. a file in a
        # directory) causes t-coffee to completely fail. It warns on all systems.
        # The -no_warning option solves this, but there is still some strange
        # bug when doing certain profile-related things. This is magically solved
        # by copying the profile file to a temp file in the current directory, so
        # it its filename supplied to t-coffee contains no /
        # (It's messy here - I just do this to /all/ input files to most easily
        #  catch all variants of providing a profile - it may only be the last
        #  form (isa("Bio::PrimarySeqI")) that causes a problem?)

        my (undef, undef, $adjustedfilename) = File::Spec->splitpath($infilename);
        if ($adjustedfilename ne $infilename) {
            my ($fh, $tempfile) = $self->io->tempfile(-dir => cwd());
            seek($IN, 0, 0);
            while (<$IN>) {
                print $fh $_;
            }
            close($fh);
            (undef, undef, $tempfile) = File::Spec->splitpath($tempfile);
            $infilename = $tempfile;
            $type = 'S';
        }

        close($IN);
        return ($infilename,$type);
    } elsif (ref($input) =~ /ARRAY/i ) {
        #  $input may be an array of BioSeq objects...
        #  Open temporary file for both reading & writing of array
        ($tfh,$infilename) = $self->io->tempfile(-dir => cwd());
        (undef, undef, $infilename) = File::Spec->splitpath($infilename);
        if( ! ref($input->[0]) ) {
            $self->warn("passed an array ref which did not contain objects to _setinput");
            return undef;
        } elsif( $input->[0]->isa('Bio::PrimarySeqI') ) {
            $temp =  Bio::SeqIO->new('-fh' => $tfh,
                         '-format' => 'fasta');
            my $ct = 1;
            foreach $seq (@$input) {
                return 0 unless ( ref($seq) &&
                          $seq->isa("Bio::PrimarySeqI") );
                if( ! defined $seq->display_id ||
                    $seq->display_id =~ /^\s+$/) {
                    $seq->display_id( "Seq".$ct++);
                }
                $temp->write_seq($seq);
            }
            $temp->close();
            undef $temp;
            close($tfh);
            $tfh = undef;
            $type = 'S';
        } elsif( $input->[0]->isa('Bio::Align::AlignI' ) ) {
            $temp =  Bio::AlignIO->new('-fh' => $tfh,
                           '-format' => $self->aformat);
            foreach my $aln (@$input) {
                next unless ( ref($aln) &&
                          $aln->isa("Bio::Align::AlignI") );
                $temp->write_aln($aln);
            }
            $temp->close();
            undef $temp;
            close($tfh);
            $tfh = undef;
            $type = 'A';
        }  else {
            $self->warn( "got an array ref with 1st entry ".
                 $input->[0].
                 " and don't know what to do with it\n");
        }
        return ($infilename,$type);
        #  $input may be a SimpleAlign object.
    } elsif ( $input->isa("Bio::Align::AlignI") ) {
        #  Open temporary file for both reading & writing of SimpleAlign object
        ($tfh, $infilename) = $self->io->tempfile(-dir => cwd());
        (undef, undef, $infilename) = File::Spec->splitpath($infilename);
        $temp =  Bio::AlignIO->new(-fh=>$tfh,
                       '-format' => 'clustalw');
        $temp->write_aln($input);
        close($tfh);
        undef $tfh;
        return ($infilename,'A');
    }

    #  or $input may be a single BioSeq object (to be added to
    # a previous alignment)
    elsif ( $input->isa("Bio::PrimarySeqI")) {
        #  Open temporary file for both reading & writing of BioSeq object
        ($tfh,$infilename) = $self->io->tempfile(-dir => cwd());
        (undef, undef, $infilename) = File::Spec->splitpath($infilename);
        $temp =  Bio::SeqIO->new(-fh=> $tfh, '-format' =>'Fasta');
        $temp->write_seq($input);
        $temp->close();
        close($tfh);
        undef $tfh;
        return ($infilename,'S');
    } else {
        $self->warn("Got $input and don't know what to do with it\n");
    }
    return 0;
}



sub _setparams {
    my ($self) = @_;
    my ($attr, $value,$param_string);
    $param_string = '';
    my $laststr;
    for  $attr ( @TCOFFEE_PARAMS ) {
        $value = $self->$attr();
        next unless (defined $value);
        my $attr_key = lc $attr;
        if( $attr_key =~ /matrix/ ) {
            $self->{'_in'} = [ "X".lc($value) ];
        } else {
            $attr_key = ' -'.$attr_key;
            $param_string .= $attr_key .'='.$value;
        }
   }
    for  $attr ( @TCOFFEE_SWITCHES) {
        $value = $self->$attr();
        next unless ($value);
        my $attr_key = lc $attr; #put switches in format expected by tcoffee
        $attr_key = ' -'.$attr_key;
        $param_string .= $attr_key ;
    }

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

    if ($self->quiet() || $self->verbose < 0) { $param_string .= ' -quiet';}

    # -no_warning is required on some systems with certain versions or failure
    # is guaranteed
    if ($self->version >= 4 && $self->version < 4.7) {
        $param_string .= ' -no_warning';
    }

    return $param_string;
}


sub aformat{
    my $self = shift;
    $self->{'_aformat'} = shift if @_;
    return $self->{'_aformat'};
}



sub methods{
   my ($self) = shift;

   @{$self->{'_methods'}} = @{shift || []} if @_;
   return @{$self->{'_methods'} || []};
}










1;

__END__