File Coverage

blib/lib/Bio/Tools/Run/Promoterwise.pm
Criterion Covered Total %
statement 40 117 34.1
branch 5 38 13.1
condition 0 14 0.0
subroutine 11 18 61.1
pod 5 5 100.0
total 61 192 31.7


line stmt bran cond sub pod time code
1             #
2             # Please direct questions and support issues to
3             #
4             # Cared for by Shawn Hoon
5             #
6             # Copyright Shawn Hoon
7             #
8             # You may distribute this module under the same terms as perl itself
9             # POD documentation - main docs before the code
10              
11             =head1 NAME
12              
13             Bio::Tools::Run::Promoterwise - Wrapper for aligning two sequences using
14             promoterwise
15              
16             =head1 SYNOPSIS
17              
18             # Build a Promoterwise alignment factory
19             my @params = ('-s'=>1,'-query_start'=>10,'-dymem'=>1);
20             my $factory = Bio::Tools::Run::Promoterwise->new(@params);
21              
22             my (@fp)= $factory->run($seq1,$seq2);
23              
24             # each feature pair is a group of hsps
25             foreach my $fp(@fp){
26             print "Hit Length: ".$fp->feature1->length."\n";
27             print "Hit Start: ".$fp->feature1->start."\n";
28             print "Hit End: ".$fp->feature1->end."\n";
29             print "Hsps: \n";
30             my @first_hsp = $fp->feature1->sub_SeqFeature;
31             my @second_hsp = $fp->feature2->sub_SeqFeature;
32             for ($i..$#first_hsp){
33             print $first_hsp[$i]->seq." ".$second_hsp[$i]->seq."\n";
34             }
35             }
36             print "end: ". $fp->feature2->start."\t".$fp->feature2->end."\n";
37              
38             #Available parameters include:
39             #( S T U V QUERY_START QUERY_END TARGET_START
40             #TARGET_END LHWINDOW LHSEED LHALN LHSCORE LHREJECT
41             #LHREJECT LHMAX DYMEM KBYTE DYCACHE)
42             #For an explanation of these parameters, please see documentation
43             #from the Wise package.
44              
45             =head1 DESCRIPTION
46              
47             Promoterwise is an alignment algorithm that relaxes the constraint
48             that local alignments have to be co-linear. Otherwise it provides a
49             similar model to DBA, which is designed for promoter sequence
50             alignments by Ewan Birney. It is part of the wise2 package available
51             at: http://www.sanger.ac.uk/software/wise2.
52              
53             =head1 FEEDBACK
54              
55             =head2 Mailing Lists
56              
57             User feedback is an integral part of the evolution of this and other
58             Bioperl modules. Send your comments and suggestions preferably to one
59             of the Bioperl mailing lists. Your participation is much appreciated.
60              
61             bioperl-l@bioperl.org - General discussion
62             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
63              
64             =head2 Support
65              
66             Please direct usage questions or support issues to the mailing list:
67              
68             I
69              
70             rather than to the module maintainer directly. Many experienced and
71             reponsive experts will be able look at the problem and quickly
72             address it. Please include a thorough description of the problem
73             with code and data examples if at all possible.
74              
75             =head2 Reporting Bugs
76              
77             Report bugs to the Bioperl bug tracking system to help us keep track
78             the bugs and their resolution. Bug reports can be submitted via
79             the web:
80              
81             http://redmine.open-bio.org/projects/bioperl/
82              
83             =head1 AUTHOR - Shawn Hoon
84              
85             Email: shawnh@fugu-sg.org
86              
87              
88             =head1 APPENDIX
89              
90             The rest of the documentation details each of the object
91             methods. Internal methods are usually preceded with a _
92              
93             =cut
94              
95             package Bio::Tools::Run::Promoterwise;
96 1         70 use vars qw($AUTOLOAD @ISA $PROGRAM $PROGRAMDIR $PROGRAMNAME
97             @PROMOTERWISE_SWITCHES @PROMOTERWISE_PARAMS
98 1     1   98390 @OTHER_SWITCHES %OK_FIELD);
  1         2  
99 1     1   448 use Bio::SeqIO;
  1         38028  
  1         24  
100 1     1   5 use Bio::Root::Root;
  1         1  
  1         13  
101 1     1   431 use Bio::Tools::Run::WrapperBase;
  1         1  
  1         22  
102 1     1   415 use Bio::Tools::Promoterwise;
  1         42691  
  1         35  
103 1     1   8 use strict;
  1         1  
  1         71  
104              
105             @ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase );
106              
107             # Two ways to run the program .....
108             # 1. define an environmental variable WISEDIR
109             # export WISEDIR =/usr/local/share/wise2.2.0
110             # where the wise2.2.20 package is installed
111             #
112             # 2. include a definition of an environmental variable WISEDIR in
113             # every script that will use DBA.pm
114             # $ENV{WISEDIR} = '/usr/local/share/wise2.2.20';
115              
116             BEGIN {
117 1     1   4 @PROMOTERWISE_PARAMS = qw( S T U V QUERY_START QUERY_END TARGET_START
118             TARGET_END LHWINDOW LHSEED LHALN LHSCORE LHREJECT
119             LHREJECT LHMAX DYMEM KBYTE DYCACHE);
120            
121              
122 1         2 @OTHER_SWITCHES = qw(SILENT QUIET ERROROFFSTD);
123              
124             # Authorize attribute fields
125 1         2 foreach my $attr ( @PROMOTERWISE_PARAMS, @PROMOTERWISE_SWITCHES,
126 21         746 @OTHER_SWITCHES) { $OK_FIELD{$attr}++; }
127             }
128              
129             =head2 program_name
130              
131             Title : program_name
132             Usage : $factory>program_name()
133             Function: holds the program name
134             Returns: string
135             Args : None
136              
137             =cut
138              
139             sub program_name {
140 6     6 1 25 return 'promoterwise';
141             }
142              
143             =head2 program_dir
144              
145             Title : program_dir
146             Usage : $factory->program_dir(@params)
147             Function: returns the program directory, obtained from ENV variable.
148             Returns: string
149             Args :
150              
151             =cut
152              
153             sub program_dir {
154 3 50   3 1 12 return Bio::Root::IO->catfile($ENV{WISEDIR},"/src/bin/") if $ENV{WISEDIR};
155             }
156              
157             sub new {
158 1     1 1 73 my ($class, @args) = @_;
159 1         12 my $self = $class->SUPER::new(@args);
160              
161 1         33 my ($attr, $value);
162 1         3 while (@args) {
163 3         4 $attr = shift @args;
164 3         4 $value = shift @args;
165 3 100       7 next if( $attr =~ /^-/ ); # don't want named parameters
166 2         12 $self->$attr($value);
167             }
168 1         2 return $self;
169             }
170              
171              
172             sub AUTOLOAD {
173 1     1   2 my $self = shift;
174 1         2 my $attr = $AUTOLOAD;
175 1         4 $attr =~ s/.*:://;
176 1         1 $attr = uc $attr;
177 1 50       5 $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
178 1 50       3 $self->{$attr} = shift if @_;
179 1         2 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 0         return undef unless $self->executable;
197 0           my $prog = $self->executable;
198 0           my $string = `$prog -version`;
199 0           $string =~ /(Version *)/i;
200 0   0       return $1 || undef;
201              
202             }
203              
204              
205             =head2 run
206              
207             Title : run
208             Usage : 2 sequence objects
209             @fp = $factory->run($seq1, $seq2);
210             Function: run
211             Returns : An array of
212             Args : Name of a file containing a set of 2 fasta sequences
213             or else 2 Bio::Seq objects.
214              
215             Throws an exception if argument is not either a string (eg a filename)
216             or 2 Bio::Seq objects. If arguments are strings, throws exception if
217             file corresponding to string name can not be found.
218              
219             =cut
220              
221             sub run{
222 0     0 1   my ($self, $seq1, $seq2)=@_;
223 0           my ($attr, $value, $switch);
224 0           $self->io->_io_cleanup();
225             # Create input file pointer
226 0           my ($infile1,$infile2)= $self->_setinput($seq1, $seq2);
227 0 0 0       if (!($infile1 && $infile2))
228 0           {$self->throw("Bad input data (sequences need an id ) ");}
229              
230              
231             # run promoterwise
232 0           my @fp = $self->_run($infile1,$infile2);
233 0           return @fp;
234             }
235              
236              
237             =head2 _run
238              
239             Title : _run
240             Usage : Internal function, not to be called directly
241             Function: makes actual system call to a promoterwise program
242             Example :
243             Returns : L
244             Args : Name of a files containing 2 sequences in the order of peptide and genomic
245              
246             =cut
247              
248             sub _run {
249 0     0     my ($self,$infile1,$infile2) = @_;
250 0           my $instring;
251 0           $self->debug( "Program ".$self->executable."\n");
252 0 0         unless ( $self->executable ) {
253 0           $self->throw("Cannot run Promoterwise unless the executable is found.".
254             " Check your environment variables or make sure ".
255             "promoterwise is in your path.");
256             }
257 0           my $paramstring = $self->_setparams;
258            
259 0           my $commandstring = $self->executable." $infile1 $infile2 $paramstring";
260             # this is to capture STDERR messages which leak out when you run programs
261             # with open(FH, "... |");
262 0 0 0       if( ( $self->silent && $self->quiet) &&
      0        
263             ($^O !~ /os2|dos|amigaos/) ) {
264 0 0         my $null = ($^O =~ m/mswin/i) ? 'NUL' : '/dev/null';
265 0           $commandstring .= " -quiet -silent -erroroffstd 2> $null";
266             }
267 0           $self->debug( "promoterwise command = $commandstring");
268 0 0         open(PW, "$commandstring |") ||
269             $self->throw( "Promoterwise call ($commandstring) crashed: $? \n");
270 0           my $pw_parser = Bio::Tools::Promoterwise->new(-fh=>\*PW,
271             -query1_seq=>$self->_query1_seq,
272             -query2_seq=>$self->_query2_seq);
273 0           my @fp;
274 0           while (my $fp = $pw_parser->next_result){
275 0           push @fp,$fp;
276             }
277              
278 0           return @fp;
279             }
280              
281             sub _setinput {
282 0     0     my ($self, $seq1, $seq2) = @_;
283 0           my ($tfh1,$tfh2,$outfile1,$outfile2);
284              
285 0 0 0       $self->throw("calling with not enough arguments") unless $seq1 && $seq2;
286              
287             # Not going to set _query_pep/_subject_dna_seq
288             # if you pass in a filename
289              
290 0 0         unless( ref($seq1) ) {
291 0 0         unless( -e $seq1 ) {
292 0           $self->throw("Sequence1 is not a Bio::PrimarySeqI object nor file\n");
293             }
294 0           $outfile1 = $seq1;
295             }
296             else {
297 0           ($tfh1,$outfile1) = $self->io->tempfile(-dir=>$self->tempdir);
298 0           my $out1 = Bio::SeqIO->new('-fh' => $tfh1,
299             '-format' => 'fasta');
300 0           $out1->write_seq($seq1);
301 0           $self->_query1_seq($seq1);
302             # Make sure you close things - this is what creates
303             # Out of filehandle errors
304 0           close($tfh1);
305 0           undef $tfh1;
306             }
307 0 0         unless( ref($seq2) ) {
308 0 0         unless( -e $seq2 ) {
309 0           $self->throw("Sequence2 is not a Bio::PrimarySeqI object nor file\n");
310             }
311 0           $outfile2 = $seq2;
312             }
313             else {
314 0           ($tfh2,$outfile2) = $self->io->tempfile(-dir=>$self->tempdir);
315 0           my $out2 = Bio::SeqIO->new('-fh' => $tfh2,
316             '-format' => 'fasta');
317 0           $out2->write_seq($seq2);
318              
319 0           $self->_query2_seq($seq2);
320             # Make sure you close things - this is what creates
321             # Out of filehandle errors
322 0           close($tfh2);
323 0           undef $tfh2;
324             }
325 0           return ($outfile1,$outfile2);
326             }
327              
328             =head2 _setparams
329              
330             Title : _setparams
331             Usage : Internal function, not to be called directly
332             Function: creates a string of params to be used in the command string
333             Example :
334             Returns : string of params
335             Args :
336              
337             =cut
338              
339             sub _setparams {
340 0     0     my ($self) = @_;
341 0           my $param_string;
342 0           foreach my $attr(@PROMOTERWISE_PARAMS){
343 0           my $value = $self->$attr();
344 0 0         next unless (defined $value);
345 0           my $attr_key = ' -'.(lc $attr);
346 0           $param_string .= $attr_key.' '.$value;
347             }
348 0           foreach my $attr(@PROMOTERWISE_SWITCHES){
349 0           my $value = $self->$attr();
350 0 0         next unless (defined $value);
351 0           my $attr_key = ' -'.(lc $attr);
352 0           $param_string .=$attr_key;
353             }
354              
355 0           $param_string = $param_string." -hitoutput tab"; #specify the output option
356 0           return $param_string;
357             }
358              
359             =head2 _query_pep_seq
360              
361             Title : _query_pep_seq
362             Usage : Internal function, not to be called directly
363             Function: get/set for the query sequence
364             Example :
365             Returns :
366             Args :
367              
368             =cut
369              
370             sub _query1_seq {
371 0     0     my ($self,$seq) = @_;
372 0 0         if(defined $seq){
373 0           $self->{'_query1_seq'} = $seq;
374             }
375 0           return $self->{'_query1_seq'};
376             }
377              
378             =head2 _subject_dna_seq
379              
380             Title : _subject_dna_seq
381             Usage : Internal function, not to be called directly
382             Function: get/set for the subject sequence
383             Example :
384             Returns :
385              
386             Args :
387              
388             =cut
389              
390             sub _query2_seq{
391 0     0     my ($self,$seq) = @_;
392 0 0         if(defined $seq){
393 0           $self->{'_query2_seq'} = $seq;
394             }
395 0           return $self->{'_query2_seq'};
396             }
397             1;