File Coverage

blib/lib/Bio/Tools/Run/Alignment/DBA.pm
Criterion Covered Total %
statement 54 236 22.8
branch 5 56 8.9
condition 0 16 0.0
subroutine 15 23 65.2
pod 5 5 100.0
total 79 336 23.5


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::Tools::Run::Alignment::DBA
3             #
4             # Copyright Shawn Hoon
5             #
6             # You may distribute this module under the same terms as perl itself
7             # POD documentation - main docs before the code
8              
9             =head1 NAME
10              
11             Bio::Tools::Run::Alignment::DBA - Object for the alignment of two
12             sequences using the DNA Block Aligner program.
13              
14             =head1 SYNOPSIS
15              
16             use Bio::Tools::Run::Alignment::DBA;
17              
18             # Build a dba alignment factory
19             my @params = ('matchA' => 0.75,
20             'matchB' => '0.55',
21             'dymem' =>'linear');
22             my $factory = Bio::Tools::Run::Alignment::DBA->new(@params);
23              
24             # Pass the factory a filename with 2 sequences to be aligned.
25             $inputfilename = 't/data/dbaseq.fa';
26             # @hsps is an array of GenericHSP objects
27             my @hsps = $factory->align($inputfilename);
28              
29             # or
30             my @files = ('t/data/dbaseq1.fa','t/data/dbaseq2.fa');
31             my @hsps = $factory->align(\@files);
32              
33             # or where @seq_array is an array of Bio::Seq objects
34             $seq_array_ref = \@seq_array;
35             my @hsps = $factory->align($seq_array_ref);
36              
37             =head1 DESCRIPTION
38              
39             DNA Block Aligner program (DBA) was developed by Ewan Birney. DBA
40             is part of the Wise package available at
41             L.
42              
43             You will need to enable dba to find the dba program. This can
44             be done in a few different ways:
45              
46             1. Define an environmental variable WISEDIR:
47              
48             export WISEDIR =/usr/local/share/wise2.2.0
49              
50             2. Include a definition of an environmental variable WISEDIR in
51             every script that will use DBA.pm:
52              
53             $ENV{WISEDIR} = '/usr/local/share/wise2.2.20';
54              
55             3. Make sure that the dba application is in your PATH.
56              
57             =head1 FEEDBACK
58              
59             =head2 Mailing Lists
60              
61             User feedback is an integral part of the evolution of this and other
62             Bioperl modules. Send your comments and suggestions preferably to one
63             of the Bioperl mailing lists. Your participation is much appreciated.
64              
65             bioperl-l@bioperl.org - General discussion
66             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
67              
68             =head2 Support
69              
70             Please direct usage questions or support issues to the mailing list:
71              
72             I
73              
74             rather than to the module maintainer directly. Many experienced and
75             reponsive experts will be able look at the problem and quickly
76             address it. Please include a thorough description of the problem
77             with code and data examples if at all possible.
78              
79             =head2 Reporting Bugs
80              
81             Report bugs to the Bioperl bug tracking system to help us keep track
82             the bugs and their resolution. Bug reports can be submitted via the
83             web:
84              
85             http://redmine.open-bio.org/projects/bioperl/
86              
87             =head1 AUTHOR - Shawn Hoon
88              
89             Email shawnh@fugu-sg.org
90              
91             =head1 APPENDIX
92              
93             The rest of the documentation details each of the object
94             methods. Internal methods are usually preceded with a _
95              
96             =cut
97              
98             package Bio::Tools::Run::Alignment::DBA;
99 1         69 use vars qw($AUTOLOAD @ISA $PROGRAM $PROGRAMDIR $PROGRAMNAME
100 1     1   111823 @DBA_SWITCHES @DBA_PARAMS @OTHER_SWITCHES %OK_FIELD);
  1         1  
101 1     1   4 use strict;
  1         1  
  1         14  
102 1     1   462 use Bio::SeqIO;
  1         38115  
  1         23  
103 1     1   628 use Bio::SimpleAlign;
  1         50612  
  1         34  
104 1     1   489 use Bio::AlignIO;
  1         959  
  1         23  
105 1     1   4 use Bio::Root::Root;
  1         1  
  1         12  
106 1     1   2 use Bio::Root::IO;
  1         2  
  1         14  
107 1     1   384 use Bio::Factory::ApplicationFactoryI;
  1         99  
  1         20  
108 1     1   581 use Bio::Search::HSP::GenericHSP;
  1         10211  
  1         34  
109 1     1   397 use Bio::Tools::Run::WrapperBase;
  1         2  
  1         70  
110              
111             @ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase);
112              
113             BEGIN {
114              
115 1     1   5 @DBA_PARAMS = qw(MATCHA MATCHB MATCHC MATCHD GAP BLOCKOPEN UMATCH SINGLE
116             NOMATCHN PARAMS KBYTE DYMEM DYDEBUG ERRORLOG);
117 1         1 @OTHER_SWITCHES = qw(OUTFILE);
118 1         1 @DBA_SWITCHES = qw(HELP SILENT QUIET ERROROFFSTD ALIGN LABEL);
119              
120             # Authorize attribute fields
121 1         2 foreach my $attr ( @DBA_PARAMS, @DBA_SWITCHES,
122 21         1585 @OTHER_SWITCHES) { $OK_FIELD{$attr}++; }
123             }
124              
125             =head2 program_name
126              
127             Title : program_name
128             Usage : $factory>program_name()
129             Function: holds the program name
130             Returns: string
131             Args : None
132              
133             =cut
134              
135             sub program_name {
136 6     6 1 23 return 'dba';
137             }
138              
139             =head2 program_dir
140              
141             Title : program_dir
142             Usage : $factory->program_dir(@params)
143             Function: returns the program directory, obtained from ENV variable.
144             Returns: string
145             Args :
146              
147             =cut
148              
149             sub program_dir {
150 3 50   3 1 10 return Bio::Root::IO->catfile($ENV{WISEDIR},"/src/bin") if $ENV{WISEDIR};
151             }
152              
153             sub new {
154 1     1 1 75 my ($class, @args) = @_;
155 1         9 my $self = $class->SUPER::new(@args);
156              
157 1         33 my ($attr, $value);
158              
159 1         3 while (@args) {
160 3         4 $attr = shift @args;
161 3         2 $value = shift @args;
162 3 50       6 next if( $attr =~ /^-/ ); # don't want named parameters
163 3 50       6 if ($attr =~/'PROGRAM'/i ) {
164 0         0 $self->executable($value);
165 0         0 next;
166             }
167 3         14 $self->$attr($value);
168             }
169 1         3 return $self;
170             }
171              
172             sub AUTOLOAD {
173 3     3   11 my $self = shift;
174 3         4 my $attr = $AUTOLOAD;
175 3         7 $attr =~ s/.*:://;
176 3         4 $attr = uc $attr;
177 3 50       6 $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
178 3 50       4 $self->{$attr} = shift if @_;
179 3         7 return $self->{$attr};
180             }
181              
182             =head2 version
183              
184             Title : version
185             Usage : exit if $prog->version() < 1.8
186             Function: Determine the version number of the program
187             Example :
188             Returns : float or undef
189             Args : none
190              
191             =cut
192              
193             sub version {
194 0     0 1   my ($self) = @_;
195              
196 0           my $exe = $self->executable();
197 0 0         return undef unless defined $exe;
198 0           my $string = `$exe -- ` ;
199 0           $string =~ /\(([\d.]+)\)/;
200 0   0       return $1 || undef;
201             }
202              
203             =head2 align
204              
205             Title : align
206             Usage :
207             $inputfilename = 't/data/seq.fa';
208             @hsps = $factory->align($inputfilename);
209             or
210             #@seq_array is array of Seq objs
211             $seq_array_ref = \@seq_array;
212             @hsps = $factory->align($seq_array_ref);
213             or
214             my @files = ('t/data/seq1.fa','t/data/seq2.fa');
215             @hsps = $factory->align(\@files);
216             Function: Perform a DBA alignment
217              
218              
219             Returns : An array of Bio::Search::HSP::GenericHSP objects
220             Args : Name of a file containing a set of 2 fasta sequences
221             or else a reference to an array to 2 Bio::Seq objects.
222             or else a reference to an array of 2 file
223             names containing 1 fasta sequence each
224              
225             Throws an exception if argument is not either a string (eg a
226             filename) or a reference to an array of 2 Bio::Seq objects. If
227             argument is string, throws exception if file corresponding to string
228             name can not be found. If argument is Bio::Seq array, throws
229             exception if less than two sequence objects are in array.
230              
231             =cut
232              
233             sub align {
234              
235 0     0 1   my ($self,$input) = @_;
236 0           my ($temp,$infile1, $infile2, $seq);
237 0           my ($attr, $value, $switch);
238              
239             # Create input file pointer
240 0           ($infile1,$infile2)= $self->_setinput($input);
241 0 0 0       if (!($infile1 && $infile2)) {$self->throw("Bad input data (sequences need an id ) or less than 2 sequences in $input !");}
  0            
242              
243             # Create parameter string to pass to dba program
244 0           my $param_string = $self->_setparams();
245              
246             # run dba
247 0           my @hsps = $self->_run($infile1,$infile2,$param_string);
248 0           return @hsps;
249             }
250              
251             #################################################
252              
253             =head2 _run
254              
255             Title : _run
256             Usage : Internal function, not to be called directly
257             Function: makes actual system call to dba program
258             Example :
259             Returns : nothing; dba output is written to a temp file
260             Args : Name of a file containing a set of unaligned fasta sequences
261             and hash of parameters to be passed to dba
262              
263             =cut
264              
265             sub _run {
266 0     0     my ($self,$infile1,$infile2,$param_string) = @_;
267 0           my $instring;
268 0           $self->debug( "Program ".$self->executable."\n");
269 0 0         unless( $self->outfile){
270 0           my ($tfh, $outfile) = $self->io->tempfile(-dir=>$self->tempdir);
271 0           close($tfh);
272 0           undef $tfh;
273 0           $self->outfile($outfile);
274             }
275 0           my $outfile = $self->outfile();
276 0           my $commandstring = $self->executable." $param_string -pff $infile1 $infile2 > $outfile";
277 0           $self->debug( "dba command = $commandstring");
278 0           my $status = system($commandstring);
279 0 0         $self->throw( "DBA call ($commandstring) crashed: $? \n") unless $status==0;
280             #parse pff format and return a Bio::Search::HSP::GenericHSP array
281 0           my $hsps = $self->_parse_results($outfile);
282              
283 0           return @{$hsps};
  0            
284             }
285              
286             =head2 _parse_results
287              
288             Title : __parse_results
289             Usage : Internal function, not to be called directly
290             Function: Parses dba output
291             Example :
292             Returns : an reference to an array of GenericHSPs
293             Args : the name of the output file
294              
295             =cut
296              
297             sub _parse_results {
298 0     0     my ($self,$outfile) = @_;
299 0 0         $outfile||$self->throw("No outfile specified");
300 0           my ($start,$end,$name,$seqname,$seq,$seqchar,$tempname,%align);
301 0           my $count = 0;
302 0           my @hsps;
303 0           open(OUT,$outfile);
304 0           my (%query,%subject);
305 0           while(my $entry = ){
306 0 0         if($entry =~ /^>(.+)/ ) {
307 0           $tempname = $1;
308 0 0         if( defined $name ) {
309 0 0         if($count == 0){
    0          
310 0           my @parse = split("\t",$name);
311 0           $query{seqname} = $parse[0];
312 0           $query{start} = $parse[3];
313 0           $query{end} = $parse[4];
314 0           $query{score} = $parse[5];
315 0 0         $query{strand} = ($parse[6] eq '+') ? 1 : -1;
316 0           my @tags = split(";",$parse[8]);
317 0           foreach my $tag(@tags){
318 0           $tag =~/(\S+)\s+(\S+)/;
319 0           $query{$1} = $2;
320             }
321 0           $query{seq} = $seqchar;
322 0           $count++;
323             }
324             elsif ($count == 1){
325 0           my @parse = split("\t",$name);
326 0           $subject{seqname} = $parse[0];
327 0           $subject{start} = $parse[3];
328 0           $subject{end} = $parse[4];
329 0           $subject{score} = $parse[5];
330 0 0         $subject{strand} = ($parse[6] eq '+') ? 1:-1;
331 0           my @tags = split(";",$parse[8]);
332 0           foreach my $tag(@tags){
333 0           $tag =~/(\S+)\s+(\S+)/;
334 0           $subject{$1} = $2;
335             }
336 0           $subject{seq} = $seqchar;
337             #create homology string
338 0           my $xor = $query{seq}^$subject{seq};
339 0           my $identical = $xor=~tr/\c@/*/;
340 0           $xor=~tr/*/ /c;
341             my $hsp= Bio::Search::HSP::GenericHSP->new(-algorithm =>'DBA',
342             -score =>$query{score},
343             -hsp_length =>length($query{seq}),
344             -query_gaps =>$query{gaps},
345             -hit_gaps =>$subject{gaps},
346             -query_name =>$query{seqname},
347             -query_start =>$query{start},
348             -query_end =>$query{end},
349             -hit_name =>$subject{seqname},
350             -hit_start =>$subject{start},
351             -hit_end =>$subject{end},
352             -hit_length =>length($self->_subject_seq->seq),
353             -query_length =>length($self->_query_seq->seq),
354             -query_seq =>$query{seq},
355             -hit_seq =>$subject{seq},
356 0           -conserved =>$identical,
357             -identical =>$identical,
358             -homology_seq =>$xor);
359 0           push @hsps, $hsp;
360 0           $count = 0;
361             }
362             }
363 0           $name = $tempname;
364 0           $seqchar = "";
365 0           next;
366             }
367 0           $entry =~ s/[^A-Za-z\.\-]//g;
368 0           $seqchar .= $entry;
369             }
370             #do for the last entry
371 0 0         if($count == 1){
372 0           my @parse = split("\t",$name);
373 0           $subject{seqname} = $parse[1];
374 0           $subject{start} = $parse[3];
375 0           $subject{end} = $parse[4];
376 0           $subject{score} = $parse[5];
377 0 0         $subject{strand} = ($parse[6] eq '+') ? 1:-1;
378 0           my @tags = split(";",$parse[8]);
379 0           foreach my $tag(@tags){
380 0           $tag =~/(\S+)\s+(\S+)/;
381 0           $subject{$1} = $2;
382             }
383 0           $subject{seq} = $seqchar;
384              
385             #create homology string
386            
387 0           my $xor = $query{seq}^$subject{seq};
388 0           my $identical = $xor=~tr/\c@/*/;
389 0           $xor=~tr/*/ /c;
390             my $hsp= Bio::Search::HSP::GenericHSP->new(-algorithm =>'DBA',
391             -score =>$query{score},
392             -hsp_length =>length($query{seq}),
393             -query_gaps =>$query{gaps},
394             -hit_gaps =>$subject{gaps},
395             -query_name =>$query{seqname},
396             -query_start =>$query{start},
397             -query_end =>$query{end},
398             -hit_name =>$subject{seqname},
399             -hit_start =>$subject{start},
400             -hit_end =>$subject{end},
401             -hit_length =>length($self->_subject_seq->seq),
402             -query_length =>length($self->_query_seq->seq),
403             -query_seq =>$query{seq},
404             -hit_seq =>$subject{seq},
405 0           -conserved =>$identical,
406             -identical =>$identical,
407             -homology_seq =>$xor);
408 0           push @hsps, $hsp;
409             }
410              
411              
412 0           return \@hsps;
413             }
414              
415             =head2 _setinput()
416              
417             Title : _setinput
418             Usage : Internal function, not to be called directly
419             Function: Create input file for dba program
420             Example :
421             Returns : name of file containing dba data input
422             Args : Seq or Align object reference or input file name
423              
424             =cut
425              
426             sub _setinput {
427 0     0     my ($self, $input, $suffix) = @_;
428 0           my ($infilename, $seq, $temp, $tfh1,$tfh2,$outfile1,$outfile2);
429              
430             #there is gotta be some repetition here...need to clean up
431              
432 0 0         if (ref($input) ne "ARRAY"){ #a single file containg 2 seqeunces
433 0           $infilename = $input;
434 0 0         unless(-e $input){return 0;}
  0            
435 0           my $in = Bio::SeqIO->new(-file => $infilename , '-format' => 'Fasta');
436 0           ($tfh1,$outfile1) = $self->io->tempfile(-dir=>$self->tempdir);
437 0           ($tfh2,$outfile2) = $self->io->tempfile(-dir=>$self->tempdir);
438            
439 0           my $out1 = Bio::SeqIO->new(-fh=> $tfh1 , '-format' => 'Fasta','-flush'=>1);
440 0           my $out2 = Bio::SeqIO->new(-fh=> $tfh2 , '-format' => 'Fasta','-flush'=>1);
441              
442 0   0       my $seq1 = $in->next_seq() || return 0;
443 0   0       my $seq2 = $in->next_seq() || return 0;
444 0           $out1->write_seq($seq1);
445 0           $out2->write_seq($seq2);
446 0           $self->_query_seq($seq1);
447 0           $self->_subject_seq($seq2);
448 0           $out1->close();
449 0           $out2->close();
450 0           close($tfh1);
451 0           close($tfh2);
452 0           undef $tfh1;
453 0           undef $tfh2;
454 0           return $outfile1,$outfile2;
455             }
456             else {
457 0 0         scalar(@{$input}) == 2 || $self->throw("dba alignment can only be run on 2 sequences not.");
  0            
458 0 0 0       if(ref($input->[0]) eq ""){#passing in two file names
    0          
459 0           my $in1 = Bio::SeqIO->new(-file => $input->[0], '-format' => 'fasta');
460 0           my $in2 = Bio::SeqIO->new(-file => $input->[1], '-format' => 'fasta');
461              
462 0           ($tfh1,$outfile1) = $self->io->tempfile(-dir=>$self->tempdir);
463 0           ($tfh2,$outfile2) = $self->io->tempfile(-dir=>$self->tempdir);
464              
465 0           my $out1 = Bio::SeqIO->new(-fh=> $tfh1 , '-format' => 'fasta');
466 0           my $out2 = Bio::SeqIO->new(-fh=> $tfh2 , '-format' => 'fasta');
467              
468 0   0       my $seq1 = $in1->next_seq() || return 0;
469 0   0       my $seq2 = $in2->next_seq() || return 0;
470 0           $out1->write_seq($seq1);
471 0           $out2->write_seq($seq2);
472 0           $self->_query_seq($seq1);
473 0           $self->_subject_seq($seq2);
474 0           close($tfh1);
475 0           close($tfh2);
476 0           undef $tfh1;
477 0           undef $tfh2;
478 0           return $outfile1,$outfile2;
479             }
480             elsif($input->[0]->isa("Bio::PrimarySeqI") && $input->[1]->isa("Bio::PrimarySeqI")) {
481 0           ($tfh1,$outfile1) = $self->io->tempfile(-dir=>$self->tempdir);
482 0           ($tfh2,$outfile2) = $self->io->tempfile(-dir=>$self->tempdir);
483              
484 0           my $out1 = Bio::SeqIO->new(-fh=> $tfh1 , '-format' => 'fasta');
485 0           my $out2 = Bio::SeqIO->new(-fh=> $tfh2 , '-format' => 'fasta');
486              
487 0           $out1->write_seq($input->[0]);
488 0           $out2->write_seq($input->[1]);
489 0           $self->_query_seq($input->[0]);
490 0           $self->_subject_seq($input->[1]);
491              
492 0           close($tfh1);
493 0           close($tfh2);
494 0           undef $tfh1;
495 0           undef $tfh2;
496 0           return $outfile1,$outfile2;
497             }
498             else {
499 0           return 0;
500             }
501             }
502 0           return 0;
503             }
504              
505             =head2 _setparams()
506              
507             Title : _setparams
508             Usage : Internal function, not to be called directly
509             Function: Create parameter inputs for dba program
510             Example :
511             Returns : parameter string to be passed to dba
512             during align or profile_align
513             Args : name of calling object
514              
515             =cut
516              
517             sub _setparams {
518 0     0     my ($attr, $value, $self);
519              
520 0           $self = shift;
521              
522 0           my $param_string = "";
523 0           for $attr ( @DBA_PARAMS ) {
524 0           $value = $self->$attr();
525 0 0         next unless (defined $value);
526             # next if $attr =~/outfile/i;
527              
528 0           my $attr_key = lc $attr; #put params in format expected by dba
529 0 0         if($attr_key =~ /match([ABCDabcd])/i){
530 0           $attr_key = "match".uc($1);
531             }
532 0           $attr_key = ' -'.$attr_key;
533 0           $param_string .= $attr_key.' '.$value;
534             }
535 0           for $attr ( @DBA_SWITCHES) {
536 0           $value = $self->$attr();
537 0 0         next unless ($value);
538 0           my $attr_key = lc $attr; #put switches in format expected by dba
539 0           $attr_key = ' -'.$attr_key;
540 0           $param_string .= $attr_key ;
541             }
542 0           return $param_string;
543             }
544              
545             =head2 _query_seq()
546              
547             Title : _query_seq
548             Usage : Internal function, not to be called directly
549             Function: get/set for the query sequence
550             Example :
551             Returns :
552             Args :
553              
554             =cut
555              
556             sub _query_seq {
557 0     0     my ($self,$seq) = @_;
558 0 0         if(defined $seq){
559 0           $self->{'_query_seq'} = $seq;
560             }
561 0           return $self->{'_query_seq'};
562             }
563              
564             =head2 _subject_seq()
565              
566             Title : _subject_seq
567             Usage : Internal function, not to be called directly
568             Function: get/set for the subject sequence
569             Example :
570             Returns :
571              
572             Args :
573              
574             =cut
575              
576             sub _subject_seq {
577 0     0     my ($self,$seq) = @_;
578 0 0         if(defined $seq){
579 0           $self->{'_subject_seq'} = $seq;
580             }
581 0           return $self->{'_subject_seq'};
582             }
583             1; # Needed to keep compiler happy
584