File Coverage

blib/lib/Bio/Tools/Run/Alignment/MSAProbs.pm
Criterion Covered Total %
statement 54 190 28.4
branch 3 64 4.6
condition 0 20 0.0
subroutine 15 24 62.5
pod 8 8 100.0
total 80 306 26.1


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::Tools::Run::Alignment::MSAProbs
3             #
4             # Please direct questions and support issues to
5             #
6             # Cared for by Jessen Bredeson
7             #
8             # Copyright Jessen Bredeson
9             #
10             # You may distribute this module under the same terms as perl itself
11             #
12             # POD documentation - main docs before the code
13              
14             =head1 NAME
15              
16             Bio::Tools::Run::Alignment::MSAProbs - Object for the calculation of a
17             multiple sequence alignment (MSA) from a set of unaligned sequences using
18             the MSAProbs program
19              
20             =head1 SYNOPSIS
21              
22             # Build a msaprobs alignment factory
23             $factory = Bio::Tools::Run::Alignment::MSAProbs->new(@params);
24              
25             # Pass the factory a list of sequences to be aligned.
26             $inputfilename = 't/cysprot.fa';
27             # $aln is a SimpleAlign object.
28             $aln = $factory->align($inputfilename);
29              
30             # or where @seq_array is an array of Bio::Seq objects
31             $seq_array_ref = \@seq_array;
32             $aln = $factory->align($seq_array_ref);
33              
34             #There are various additional options and input formats available.
35             #See the DESCRIPTION section that follows for additional details.
36              
37             =head1 DESCRIPTION
38              
39             MSAProbs is Liu, Schmidt, and Maskell's (2010) alignment program using HMM
40             and partition function posterior probabilities. For more a more in-depth
41             description see the original publication:
42              
43             Liu, Y., Schmidt, B., and Maskell, D. L. (2010) MSAProbs: multiple
44             sequence alignment based on pair hidden Markov models and partition
45             function posterior probabilities. I 26(16): 1958-1964
46             doi:10.1093/bioinformatics/btq338
47            
48             -OR-
49              
50             http://bioinformatics.oxfordjournals.org/content/26/16/1958.abstract
51              
52             You can download the source code from
53             http://sourceforge.net/projects/msaprobs/
54              
55             It is recommended you use at least version 0.9; behaviour with earlier
56             versions is questionable.
57              
58             =head2 Helping the module find your executable
59              
60             You will need to help MSAProbs to find the 'msaprobs' executable. This can
61             be done in (at least) three ways:
62              
63             1. Make sure the msaprobs executable is in your path (i.e.
64             'which msaprobs' returns a valid program)
65             2. define an environmental variable MSAPROBSDIR which points to a
66             directory containing the 'msaprobs' app:
67             In bash
68             export MSAPROBSDIR=/home/progs/msaprobs or
69             In csh/tcsh
70             setenv MSAPROBSDIR /home/progs/msaprobs
71              
72             3. include a definition of an environmental variable MSAPROBSDIR
73             in every script that will
74             BEGIN {$ENV{MSAPROBSDIR} = '/home/progs/msaprobs'; }
75             use Bio::Tools::Run::Alignment::MSAProbs;
76              
77             =head1 FEEDBACK
78              
79             =head2 Mailing Lists
80              
81             User feedback is an integral part of the evolution of this and other
82             Bioperl modules. Send your comments and suggestions preferably to one
83             of the Bioperl mailing lists. Your participation is much appreciated.
84              
85             bioperl-l@bioperl.org - General discussion
86             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
87              
88             =head2 Support
89              
90             Please direct usage questions or support issues to the mailing list:
91              
92             I
93              
94             rather than to the module maintainer directly. Many experienced and
95             reponsive experts will be able look at the problem and quickly
96             address it. Please include a thorough description of the problem
97             with code and data examples if at all possible.
98              
99             =head2 Reporting Bugs
100              
101             Report bugs to the Bioperl bug tracking system to help us keep track
102             the bugs and their resolution. Bug reports can be submitted via the web:
103              
104             http://bugzilla.open-bio.org/
105              
106             =head1 AUTHOR - Jessen Bredeson
107              
108             Email jessenbredeson@berkeley.edu
109              
110             =head1 CONTRIBUTIONS
111              
112             This MSAProbs module was adapted from the Bio::Tools::Run::Alignment::Muscle
113             module, written by Jason Stajich and almost all of the credit should be given
114             to him.
115              
116             Email jason-at-bioperl-dot-org
117              
118             =head1 APPENDIX
119              
120             The rest of the documentation details each of the object
121             methods. Internal methods are usually preceded with a _
122              
123             =cut
124              
125             package Bio::Tools::Run::Alignment::MSAProbs;
126              
127 1         70 use vars qw($AUTOLOAD @ISA $PROGRAMNAME $PROGRAM %DEFAULTS %MSAPROBS_PARAMS
128             %MSAPROBS_SWITCHES %OK_FIELD
129 1     1   112708 );
  1         2  
130 1     1   4 use strict;
  1         1  
  1         16  
131 1     1   556 use Bio::Seq;
  1         48255  
  1         26  
132 1     1   508 use Bio::SeqIO;
  1         21890  
  1         31  
133 1     1   675 use Bio::SimpleAlign;
  1         28855  
  1         33  
134 1     1   488 use Bio::AlignIO;
  1         1100  
  1         24  
135 1     1   5 use Bio::Root::Root;
  1         1  
  1         14  
136 1     1   3 use Bio::Root::IO;
  1         1  
  1         15  
137 1     1   439 use Bio::Factory::ApplicationFactoryI;
  1         108  
  1         18  
138 1     1   9 use Bio::Tools::GuessSeqFormat;
  1         1  
  1         15  
139 1     1   415 use Bio::Tools::Run::WrapperBase;
  1         1  
  1         91  
140             @ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase
141             Bio::Factory::ApplicationFactoryI);
142              
143              
144             BEGIN {
145 1     1   6 %DEFAULTS = ( 'QUIET' => 1,
146             '_AFORMAT' => 'fasta',
147             '_CONSISTENCY' => 2,
148             '_ITERATIONS' => 10,
149             '_CLUSTALW' => 0,
150             '_ALIGNMENT_ORDER' => 0
151             );
152 1         2 %MSAPROBS_PARAMS = ( 'NUM_THREADS' => 'NUM_THREADS',
153             'CONSISTENCY' => 'C',
154             'ITERATIONS' => 'IR',
155             'ANNOT_FILE' => 'ANNOT'
156             );
157 1         3 %MSAPROBS_SWITCHES = ( 'CLUSTALW' => 'CLUSTALW',
158             'ALIGNMENT_ORDER' => 'A'
159             );
160            
161             # Authorize attribute fields
162 1         2 %OK_FIELD = map{ uc($_) => 1 } qw(INFILE OUTFILE VERBOSE QUIET VERSION),
  11         1419  
163             keys %MSAPROBS_PARAMS,
164             keys %MSAPROBS_SWITCHES;
165             }
166              
167             =head2 program_name
168              
169             Title : program_name
170             Usage : $factory->program_name()
171             Function: holds the program name
172             Returns: string
173             Args : None
174              
175             =cut
176              
177             sub program_name {
178 6     6 1 31 return 'msaprobs';
179             }
180              
181             =head2 program_dir
182              
183             Title : program_dir
184             Usage : $factory->program_dir(@params)
185             Function: returns the program directory, obtained from ENV variable.
186             Returns: string
187             Args :
188              
189             =cut
190              
191             sub program_dir {
192 3 50   3 1 15 return Bio::Root::IO->catfile($ENV{MSAPROBSDIR}) if $ENV{MSAPROBSDIR};
193             }
194              
195             =head2 version
196              
197             Title : version
198             Usage : exit if $prog->version() < 0.9.4
199             Function: Determine the version number of the program
200             Example :
201             Returns : float or undef
202             Args : none
203              
204            
205             =cut
206              
207             sub version {
208 0     0 1 0 my ($self) = @_;
209 0         0 my( $exe,$version );
210 0 0       0 return unless $exe = $self->executable;
211 0         0 my $string = `$exe -version 2>&1` ;
212              
213 0         0 $string =~ /MSAPROBS\s+VERSION\s+([\d\.]+)/i;
214 0 0       0 $version =~ s/\.(\d+)$/$1/ if ($version = $1);
215            
216 0   0     0 return $version || undef;
217             }
218              
219             =head2 new
220              
221             Title : new
222             Usage : my $msaprobs = Bio::Tools::Run::Alignment::MSAProbs->new();
223             Function: Constructor
224             Returns : Bio::Tools::Run::Alignment::MSAProbs
225             Args : -outfile => $outname
226              
227              
228             =cut
229              
230             sub new {
231 1     1 1 80 my ($class,@args) = @_;
232            
233 1         9 my $self = $class->SUPER::new(@args);
234            
235 1         34 my( @msap_args, @obj_args, $field );
236 1         3 while( my $arg = shift @args ) {
237 1         2 $field = uc $arg;
238 1         2 $field =~ s/^-//;
239 1 50       14 $arg = '-'.$arg if $arg !~ /^-/;
240             $self->throw("Invalid argument: $field")
241 1 50       4 unless $OK_FIELD{$field};
242 1         4 push @msap_args, lc($arg),shift @args;
243             }
244 1         3 map{ $self->{lc($_)} = $DEFAULTS{$_} } keys %DEFAULTS;
  6         10  
245            
246             $self->_set_from_args(\@msap_args,
247             -create => 1,
248             -case_sensitive => 1,
249 1         5 -methods => [map{lc($_);} keys %OK_FIELD]);
  11         18  
250            
251            
252 1         9 return $self;
253             }
254              
255             =head2 run
256              
257             Title : run
258             Usage : my $output = $application->run(\@seqs);
259             Function: Generic run of an application
260             Returns : Bio::SimpleAlign object
261             Args : Arrayref of Bio::PrimarySeqI objects or
262             a filename to run on
263              
264            
265             =cut
266              
267             sub run {
268 0     0 1   my( $self,$input ) = @_;
269 0   0       $input ||= $self->infile;
270 0           return $self->align($input);
271             }
272              
273             =head2 align
274              
275             Title : align
276             Usage :
277             $inputfilename = 't/data/cysprot.fa';
278             $aln = $factory->align($inputfilename);
279             or
280             $seq_array_ref = \@seq_array;
281             # @seq_array is array of Seq objs
282             $aln = $factory->align($seq_array_ref);
283             Function: Perform a multiple sequence alignment
284             Returns : Reference to a SimpleAlign object containing the
285             sequence alignment.
286             Args : Name of a file containing a set of unaligned fasta sequences
287             or else an array of references to Bio::Seq objects.
288              
289             Throws an exception if argument is not either a string (eg a
290             filename) or a reference to an array of Bio::Seq objects. If
291             argument is string, throws exception if file corresponding to string
292             name can not be found. If argument is Bio::Seq array, throws
293             exception if less than two sequence objects are in array.
294              
295            
296             =cut
297              
298             sub align {
299 0     0 1   my ($self,$input) = @_;
300             # Create input file pointer
301 0           $self->io->_io_cleanup();
302 0           my $infilename;
303 0 0         if( defined($input) ) {
    0          
304 0           $infilename = $self->_setinput($input);
305             } elsif( defined($self->infile) ) {
306 0           $infilename = $self->_setinput($self->infile);
307             } else {
308 0           $self->throw("No inputdata provided\n");
309             }
310 0 0         unless( $infilename ) {
311 0           $self->throw("Bad input data or less than 2 sequences in $infilename !");
312             }
313            
314 0           my $param_string = $self->_setparams();
315              
316             # run msaprobs
317 0           return &_run($self, $infilename, $param_string);
318             }
319              
320             =head2 error_string
321              
322             Title : error_string
323             Usage : $obj->error_string($newval)
324             Function: Where the output from the last analysus run is stored.
325             Returns : value of error_string
326             Args : newvalue (optional)
327              
328              
329             =cut
330              
331             sub error_string{
332 0     0 1   my ($self,$value) = @_;
333 0 0         if( defined $value) {
334 0           $self->{'error_string'} = $value;
335             }
336 0           return $self->{'error_string'};
337              
338             }
339              
340             =head2 infile
341              
342             Title : infile
343             Usage : $prog->infile($filename)
344             Function: get/set the fasta (and only a fasta) file to run on
345             or the array reference containing the Bio::SeqI objects
346             Returns : name of input sequence file or object array ref
347             Args : name of input sequence file or object array ref
348              
349            
350             =cut
351              
352             =head2 outfile
353              
354             Title : outfile
355             Usage : $prog->outfile($filename)
356             Function: get/set the file to save output to
357             Returns : outfile name if set
358             Args : newvalue (optional)
359              
360            
361             =cut
362              
363             =head2 annot_file
364              
365             Title : annot_file
366             Usage : $prog->annot_file($filename)
367             Function: get/set the file name to write the MSA annotation to
368             Returns : filename or undef
369             Args : filename (optional)
370              
371            
372             =cut
373              
374             =head2 num_threads
375              
376             Title : num_threads
377             Usage : $prog->num_threads($cores)
378             Function: get/set number of cores on your machine
379             Returns : integer
380             Args : integer (optional; executable auto-detects)
381              
382            
383             =cut
384              
385             =head2 consistency
386              
387             Title : consistency
388             Usage : $prog->consistency($passes)
389             Function: get/set the number of consistency transformation passes
390             Returns : integer
391             Args : integer 0..5, [default 2] (optional)
392              
393            
394             =cut
395              
396             =head2 iterations
397              
398             Title : iterations
399             Usage : $prog->iterations($passes)
400             Function: get/set the number of iterative-refinement passes
401             Returns : integer
402             Args : integer 0..1000, [default 10] (optional)
403              
404            
405             =cut
406              
407             =head2 alignment_order
408              
409             Title : alignment_order
410             Usage : $prog->alignment_order($bool)
411             Function: specify whether or not to output aligned sequences in
412             alignment order, not input order
413             Returns : boolean
414             Args : boolean [default: off] (optional)
415              
416            
417             =cut
418              
419             =head2 clustalw
420              
421             Title : clustalw
422             Usage : $prog->clustalw($bool)
423             Function: write output in clustalw format; makes no sense unless
424             outfile() is also specified
425             Returns : boolean
426             Args : boolean [default: off] (optional)
427              
428            
429             =cut
430              
431             =head1 Bio::Tools::Run::WrapperBase methods
432              
433             =cut
434              
435             =head2 no_param_checks
436              
437             Title : no_param_checks
438             Usage : $obj->no_param_checks($newval)
439             Function: Boolean flag as to whether or not we should
440             trust the sanity checks for parameter values
441             Returns : value of no_param_checks
442             Args : newvalue (optional)
443              
444              
445             =cut
446              
447             =head2 save_tempfiles
448              
449             Title : save_tempfiles
450             Usage : $obj->save_tempfiles($newval)
451             Function:
452             Returns : value of save_tempfiles
453             Args : newvalue (optional)
454              
455              
456             =cut
457              
458             =head2 outfile_name
459              
460             Title : outfile_name
461             Usage : my $outfile = $msaprobs->outfile_name();
462             Function: Get the name of the output file from a run
463             (if you wanted to do something special)
464             Returns : string
465             Args : none
466              
467              
468             =cut
469              
470              
471             =head2 tempdir
472              
473             Title : tempdir
474             Usage : my $tmpdir = $self->tempdir();
475             Function: Retrieve a temporary directory name (which is created)
476             Returns : string which is the name of the temporary directory
477             Args : none
478              
479              
480             =cut
481              
482             =head2 cleanup
483              
484             Title : cleanup
485             Usage : $msaprobs->cleanup();
486             Function: Will cleanup the tempdir directory
487             Returns : none
488             Args : none
489              
490              
491             =cut
492              
493             =head2 io
494              
495             Title : io
496             Usage : $obj->io($newval)
497             Function: Gets a L object
498             Returns : L
499             Args : none
500              
501              
502             =cut
503              
504             =head1 Private Methods
505              
506             =cut
507              
508             =head2 _run
509              
510             Title : _run
511             Usage : Internal function, not to be called directly
512             Function: makes actual system call to msaprobs program
513             Example :
514             Returns : nothing; msaprobs output is written to a
515             temporary file OR specified output file
516             Args : Name of a file containing a set of unaligned fasta sequences
517             and hash of parameters to be passed to msaprobs
518              
519              
520             =cut
521              
522             sub _run {
523 0     0     my ($self,$infilename,$params) = @_;
524 0           my $commandstring = $self->executable.' '.$infilename.$params;
525 0           $self->debug( "msaprobs command = $commandstring \n");
526              
527 0           my $status = system($commandstring);
528 0           my $outfile = $self->outfile_name;
529 0 0         if( !-s $outfile ) {
530 0           $self->warn( "MSAProbs call crashed: $? [command $commandstring]\n");
531 0           return undef;
532             }
533            
534 0 0         if( $self->clustalw ){
535 0           $outfile = $self->_clustalize($outfile);
536 0           $self->aformat('clustalw');
537             }
538 0           my $in = Bio::AlignIO->new(
539             '-file' => $outfile,
540             '-format' => $self->aformat,
541             '-displayname_flat' => 1
542             );
543 0           my $aln = $in->next_aln();
544 0           undef $in;
545            
546 0           return $aln;
547             }
548            
549             =head2 _setinput
550              
551             Title : _setinput
552             Usage : Internal function, not to be called directly
553             Function: Create input file for msaprobs program
554             Example :
555             Returns : name of file containing msaprobs data input AND
556             Args : Arrayref of Seqs or input file name
557              
558              
559             =cut
560              
561             sub _setinput {
562 0     0     my( $self,$input ) = @_;
563 0           my( $infilename,$outtemp,$tfh,@sequences );
564 0 0         if (! ref $input) {
    0          
565             # check that file exists or throw
566 0 0 0       return unless (-s $input && -r $input);
567             # let's peek and guess
568 0           $infilename = $input;
569 0 0         open(IN,$input) || $self->throw("Cannot open $input");
570 0           my $header;
571 0           while( defined ($header = ) ) {
572 0 0         last if $header !~ /^\s+$/; }
573 0           close(IN);
574 0 0         $header =~ /^>\s*\S+/ ||
575             $self->throw("Need to provide a FASTA-formatted file to msaprobs!");
576 0           my $inseqio = Bio::SeqIO->new(
577             -file => $input,
578             -format => 'fasta' );
579 0           while( my $seq = $inseqio->next_seq ){
580 0           push @sequences, $seq; }
581 0           undef $inseqio;
582             # have to check each seq for terminal '*', so
583             # continue below and write clean output to temp file
584             }elsif( ref($input) =~ /ARRAY/i ){ # $input may be an
585             # array of BioSeq objects...
586             # Open temporary file for both reading & writing of array
587 0 0         if( ! ref($input->[0]) ) {
    0          
588 0           $self->warn("passed an array ref which did not contain objects to _setinput");
589 0           return;
590             }elsif( $input->[0]->isa('Bio::PrimarySeqI') ){
591 0           @sequences = @$input;
592             }else{
593 0           $self->warn( "got an array ref with 1st entry ".
594             $input->[0].
595             " and don't know what to do with it\n");
596 0           return;
597             }
598             }else{
599 0           $self->warn("Got $input and don't know what to do with it\n");
600 0           return;
601             }
602            
603 0           ($tfh,$infilename) = $self->io->tempfile();
604 0           $outtemp = Bio::SeqIO->new('-fh' => $tfh,
605             '-format' => 'fasta');
606 0           my( @out,$string );
607 0           my $ct = 1;
608 0           while( my $seq = shift @sequences){
609 0 0 0       return unless ( ref($seq) &&
610             $seq->isa("Bio::PrimarySeqI") );
611 0 0 0       if( ! defined $seq->display_id ||
612             $seq->display_id =~ /^\s+$/){
613 0           $seq->display_id( "Seq".$ct++ ); }
614 0           $string = $seq->seq;
615 0           $string =~ s/\*$//;
616 0           $seq->seq($string);
617 0 0         if( $string =~ tr/~.-/~.-/ ){
618 0           $self->warn("These sequences may have already been aligned!");
619             }
620 0           push @out, $seq;
621             }
622 0           $outtemp->write_seq(@out);
623            
624 0           $outtemp->close();
625 0           undef $outtemp;
626 0           close($tfh);
627 0           $tfh = undef;
628 0           return $infilename;
629             }
630              
631              
632             =head2 _setparams
633              
634             Title : _setparams
635             Usage : Internal function, not to be called directly
636             Function: Create parameter inputs for msaprobs program
637             Example :
638             Returns : parameter string to be passed to msaprobs
639             during align
640             Args : name of calling object
641              
642            
643             =cut
644              
645             sub _setparams {
646 0     0     my ($self) = @_;
647 0           my ($attr,$method,$value,$param_string);
648 0           $param_string = '';
649            
650 0 0         unless( defined $self->outfile ){
651 0           $self->aformat($DEFAULTS{'AFORMAT'});
652 0           $self->clustalw(0);
653             }
654            
655             #put switches/params in format expected by MSAProbs
656 0           for $attr ( keys %MSAPROBS_PARAMS ){
657 0           $method = lc $attr;
658 0           $value = $self->$method();
659 0 0         next unless (defined $value);
660 0           my $attr_key = lc $MSAPROBS_PARAMS{$attr};
661 0           $attr_key = ' -'.$attr_key;
662 0           $param_string .= $attr_key.' '.$value;
663             }
664            
665 0           for $attr ( keys %MSAPROBS_SWITCHES ){
666 0           $method = lc $attr;
667 0           $value = $self->$method();
668 0 0         next unless $value;
669 0           my $attr_key = lc $MSAPROBS_SWITCHES{$attr};
670 0           $attr_key = ' -'.$attr_key;
671 0           $param_string .= $attr_key;
672             }
673              
674             # Set default output file if no explicit file specified
675             # or if a clustalw-formatted file is desired...
676 0 0 0       if( $self->clustalw || ! $self->outfile ) {
677 0           my ($tfh, $outfile) = $self->io->tempfile(-dir => $self->tempdir);
678 0           close($tfh);
679 0           undef $tfh;
680 0           $self->outfile_name($outfile);
681             }else{
682 0           $self->outfile_name($self->outfile);
683             }
684 0 0         my $null = ($^O =~ m/mswin/i) ? 'NUL' : '/dev/null';
685 0 0         $param_string .= ' -v' if $self->verbose > 0;
686 0           $param_string .= ' >'.$self->outfile_name;
687 0 0 0       $param_string .= " 2>$null" if $self->quiet &&
688             $self->verbose < 1;
689 0           $self->arguments($param_string);
690 0           return $param_string;
691             }
692              
693             =head2 aformat
694              
695             Title : aformat
696             Usage : my $alignmentformat = $self->aformat();
697             Function: Get/Set alignment format
698             Returns : string
699             Args : string
700              
701              
702             =cut
703              
704             sub aformat{
705 0     0 1   my $self = shift;
706 0 0         $self->{'_aformat'} = shift if @_;
707 0           return $self->{'_aformat'};
708             }
709              
710              
711             sub _clustalize {
712 0     0     my $self = shift;
713 0           my $infile = shift;
714 0           my $outfile = $self->outfile;
715            
716 0           local $/ = "\n";
717 0           my( $in,$out,$firstline,$line );
718 0           $in = Bio::Root::IO->new(-file => $infile);
719 0           $out = Bio::Root::IO->new(-file => '>'.$outfile);
720            
721 0           while( defined( $firstline = $in->_readline )) {
722 0 0         last if $firstline !~ /^\s*$/; }
723 0           $in->_pushback('CLUSTALW format, '.$firstline);
724            
725 0           while( defined( $line = $in->_readline )) {
726 0           $out->_print( $line ); }
727            
728 0           $out->close();
729 0           $in->close();
730 0           undef $out;
731 0           undef $in;
732 0           $self->debug($outfile);
733 0 0         return $outfile if -s $outfile;
734             }
735              
736              
737              
738             1; # Needed to keep compiler happy