view release on metacpan or  search on metacpan

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

=head2  _run

 Title   :  _run
 Usage   :  Internal function, not to be called directly	
 Function:  makes actual system call to msaprobs program
 Example :
 Returns : nothing; msaprobs output is written to a
           temporary file OR specified output file
 Args    : Name of a file containing a set of unaligned fasta sequences
           and hash of parameters to be passed to msaprobs


=cut

sub _run {
    my ($self,$infilename,$params) = @_;
    my $commandstring = $self->executable.' '.$infilename.$params;
    $self->debug( "msaprobs command = $commandstring \n");

    my $status  = system($commandstring);
    my $outfile = $self->outfile_name; 
    if( !-s $outfile ) {
	$self->warn( "MSAProbs call crashed: $? [command $commandstring]\n");
	return undef;
    }
    
    if( $self->clustalw ){
    $outfile = $self->_clustalize($outfile);
    $self->aformat('clustalw');
    }
    my $in  = Bio::AlignIO->new(
        '-file'   => $outfile, 
        '-format' => $self->aformat,
        '-displayname_flat' => 1 
        );
    my $aln = $in->next_aln();
    undef $in;
    
    return $aln;
}
    
=head2  _setinput

 Title   :  _setinput
 Usage   :  Internal function, not to be called directly	
 Function:  Create input file for msaprobs program
 Example :
 Returns : name of file containing msaprobs data input AND
 Args    : Arrayref of Seqs or input file name


=cut

sub _setinput {
    my( $self,$input ) = @_;
    my( $infilename,$outtemp,$tfh,@sequences );
    if (! ref $input) {
        # check that file exists or throw
        return unless (-s $input && -r $input);
        # let's peek and guess
        $infilename = $input;
        open(IN,$input) || $self->throw("Cannot open $input");
        my $header;
        while( defined ($header = <IN>) ) {
            last if $header !~ /^\s+$/;  }
        close(IN);
        $header =~ /^>\s*\S+/ ||
        $self->throw("Need to provide a FASTA-formatted file to msaprobs!");
	    my $inseqio = Bio::SeqIO->new(
	        -file   => $input,
	        -format => 'fasta' );
        while( my $seq = $inseqio->next_seq ){
            push @sequences, $seq;  }
        undef $inseqio;
        # have to check each seq for terminal '*', so
        # continue below and write clean output to temp file
    }elsif( ref($input) =~ /ARRAY/i ){ #  $input may be an
        #  array of BioSeq objects...
        #  Open temporary file for both reading & writing of array
        if( ! ref($input->[0]) ) {
            $self->warn("passed an array ref which did not contain objects to _setinput");
            return;
        }elsif( $input->[0]->isa('Bio::PrimarySeqI') ){		
            @sequences = @$input;
        }else{ 
            $self->warn( "got an array ref with 1st entry ".
                $input->[0].
                " and don't know what to do with it\n");
            return;
        }
    }else{ 
        $self->warn("Got $input and don't know what to do with it\n");
        return;
    }
    
    ($tfh,$infilename) = $self->io->tempfile();
    $outtemp =  Bio::SeqIO->new('-fh'     => $tfh,
			    	            '-format' => 'fasta');
    my( @out,$string );
	my $ct = 1;
    while( my $seq = shift @sequences){
        return unless ( ref($seq) && 
            $seq->isa("Bio::PrimarySeqI") );
        if( ! defined $seq->display_id ||
            $seq->display_id =~ /^\s+$/){
            $seq->display_id( "Seq".$ct++ ); } 
        $string = $seq->seq;
        $string =~ s/\*$//;
        $seq->seq($string);
        if( $string =~ tr/~.-/~.-/ ){
        $self->warn("These sequences may have already been aligned!");
        }
        push @out, $seq;
    }
    $outtemp->write_seq(@out);
    
    $outtemp->close();
    undef $outtemp;
    close($tfh);
    $tfh = undef;