File Coverage

blib/lib/Bio/Tools/Run/Pseudowise.pm
Criterion Covered Total %
statement 40 136 29.4
branch 6 40 15.0
condition 0 17 0.0
subroutine 10 20 50.0
pod 6 6 100.0
total 62 219 28.3


line stmt bran cond sub pod time code
1             # BioPerl module for Bio::Tools::Run::Pseudowise
2             #
3             # Please direct questions and support issues to
4             #
5             # Cared for by
6             #
7             # Copyright Kiran
8             #
9             # You may distribute this module under the same terms as perl itself
10             # POD documentation - main docs before the code
11              
12             =head1 NAME
13              
14             Bio::Tools::Run::Pseudowise - Object for prediting pseudogenes in a
15             given sequence given a protein and a cdna sequence
16              
17             =head1 SYNOPSIS
18              
19             # Build a pseudowise alignment factory
20             my $factory = Bio::Tools::Run::Pseudowise->new();
21              
22             # Pass the factory 3 Bio:SeqI objects (in the order of query
23             # peptide and cdna and target_genomic)
24             # @genes is an array of GenericSeqFeature objects
25             my @genes = $factory->run($seq1, $seq2, $seq3);
26              
27             =head1 DESCRIPTION
28              
29             Pseudowise is a pseudogene predition program developed by Ewan Birney
30             http://www.sanger.ac.uk/software/wise2.
31              
32             =head1 FEEDBACK
33              
34             =head2 Mailing Lists
35              
36             User feedback is an integral part of the evolution of this and other
37             Bioperl modules. Send your comments and suggestions preferably to one
38             of the Bioperl mailing lists. Your participation is much appreciated.
39              
40             bioperl-l@bioperl.org - General discussion
41             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
42              
43             =head2 Support
44              
45             Please direct usage questions or support issues to the mailing list:
46              
47             I
48              
49             rather than to the module maintainer directly. Many experienced and
50             reponsive experts will be able look at the problem and quickly
51             address it. Please include a thorough description of the problem
52             with code and data examples if at all possible.
53              
54             =head2 Reporting Bugs
55              
56             Report bugs to the Bioperl bug tracking system to help us keep track
57             the bugs and their resolution. Bug reports can be submitted via the
58             web:
59              
60             http://redmine.open-bio.org/projects/bioperl/
61              
62             =head1 AUTHOR - Kiran
63              
64             Email kiran@fugu-sg.org
65              
66             =head1 APPENDIX
67              
68             The rest of the documentation details each of the object
69             methods. Internal methods are usually preceded with a _
70              
71             =cut
72              
73              
74             package Bio::Tools::Run::Pseudowise;
75 1         70 use vars qw($AUTOLOAD @ISA $PROGRAM_NAME $PROGRAM_DIR
76             @PSEUDOWISE_SWITCHES @PSEUDOWISE_PARAMS
77 1     1   95610 @OTHER_SWITCHES %OK_FIELD);
  1         2  
78 1     1   5 use strict;
  1         1  
  1         16  
79 1     1   321 use Bio::SeqIO;
  1         35713  
  1         27  
80 1     1   286 use Bio::Tools::Run::WrapperBase;
  1         3  
  1         23  
81 1     1   290 use Bio::Tools::Pseudowise;
  1         46591  
  1         110  
82              
83             @ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase);
84              
85             # You will need to enable pseudowise to find the pseudowise program. This
86             # can be done in (at least) two ways:
87             #
88             # 1. define an environmental variable WISEDIR
89             # export WISEDIR =/usr/local/share/wise2.2.0
90             # where the wise2.2.20 package is installed
91             #
92             # 2. include a definition of an environmental variable WISEDIR in
93             # every script that will use DBA.pm
94             # $ENV{WISEDIR} = '/usr/local/share/wise2.2.20';
95              
96             BEGIN {
97 1     1   4 $PROGRAM_NAME = 'pseudowise';
98 1 50       4 $PROGRAM_DIR = Bio::Root::IO->catfile($ENV{WISEDIR},"src","bin") if $ENV{WISEDIR};
99 1         3 @PSEUDOWISE_PARAMS = qw(SPLICE_MAX_COLLAR SPLICE_MIN_COLLAR
100             SPLICE_SCORE_OFFSET
101             GENESTATS NOMATCHN PARAMS KBYTE
102             DYMEM DYDEBUG PALDEBUG
103             ERRORLOG);
104              
105 1         2 @PSEUDOWISE_SWITCHES = qw(HELP SILENT QUIET ERROROFFSTD);
106            
107             # Authorize attribute fields
108 1         2 foreach my $attr ( @PSEUDOWISE_PARAMS, @PSEUDOWISE_SWITCHES,
109 15         982 @OTHER_SWITCHES) { $OK_FIELD{$attr}++; }
110              
111             }
112              
113             =head2 program_name
114              
115             Title : program_name
116             Usage : $factory>program_name()
117             Function: holds the program name
118             Returns: string
119             Args : None
120              
121             =cut
122              
123             sub program_name {
124 6     6 1 25 return $PROGRAM_NAME;
125             }
126              
127             =head2 program_dir
128              
129             Title : program_dir
130             Usage : $factory->program_dir(@params)
131             Function: returns the program directory, obtained from ENV variable.
132             Returns: string
133             Args :
134              
135             =cut
136              
137             sub program_dir {
138 3     3 1 8 return $PROGRAM_DIR;
139             }
140              
141             sub new {
142 1     1 1 87 my ($class, @args) = @_;
143 1         10 my $self = $class->SUPER::new(@args);
144            
145 1         40 my ($attr, $value);
146 1         3 while (@args) {
147 3         6 $attr = shift @args;
148 3         4 $value = shift @args;
149 3 100       8 next if( $attr =~ /^-/ ); # don't want named parameters
150 2 50       6 if ($attr =~/'PROGRAM'/i) {
151 0         0 $self->executable($value);
152 0         0 next;
153             }
154 2         11 $self->$attr($value);
155             }
156 1         3 return $self;
157             }
158              
159              
160             sub AUTOLOAD {
161 1     1   2 my $self = shift;
162 1         2 my $attr = $AUTOLOAD;
163 1         5 $attr =~ s/.*:://;
164 1         2 $attr = uc $attr;
165 1 50       5 $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
166 1 50       3 $self->{$attr} = shift if @_;
167 1         3 return $self->{$attr};
168             }
169              
170             =head2 version
171              
172             Title : version
173             Usage : exit if $prog->version() < 1.8
174             Function: Determine the version number of the program
175             Example :
176             Returns : float or undef
177             Args : none
178              
179             =cut
180              
181             sub version {
182 0     0 1   my ($self) = @_;
183              
184 0 0         return undef unless $self->executable;
185 0           my $string = `pseudowise -- ` ;
186 0           $string =~ /\(([\d.]+)\)/;
187 0   0       return $1 || undef;
188              
189             }
190              
191             =head2 predict_genes
192              
193             Title : predict_genes
194             Usage : DEPRECATED. Use $factory->run instead
195             Function: Predict pseudogenes
196             Returns : An array of Bio::Seqfeature::Generic objects
197             Args : Name of a file containing a set of 3 fasta sequences in the order of
198             peptide, cdna and genomic sequences
199             or else 3 Bio::Seq objects.
200              
201             Throws an exception if argument is not either a string (eg a filename)
202             or 3 Bio::Seq objects. If arguments are strings, throws exception if
203             file corresponding to string name can not be found.
204              
205             =cut
206              
207             sub predict_genes {
208 0     0 1   return shift->run(@_);
209             }
210              
211             =head2 run
212              
213             Title : run
214             Usage : my @feats = $factory->run($seq1, $seq2, $seq3);
215             Function: Executes pseudogene binary
216             Returns : An array of Bio::Seqfeature::Generic objects
217             Args : Name of a file containing a set of 3 fasta sequences in the order of
218             peptide, cdna and genomic sequences
219             or else 3 Bio::Seq objects.
220              
221             Throws an exception if argument is not either a string (eg a filename)
222             or 3 Bio::Seq objects. If arguments are strings, throws exception if
223             file corresponding to string name can not be found.
224              
225             =cut
226              
227             sub run {
228 0     0 1   my ($self,@args)=@_;
229 0           my ($attr, $value, $switch);
230            
231             # Create input file pointer
232 0           my @files = $self->_setinput(@args);
233 0 0 0       if( @files !=3 || grep { !defined } @files ) {
  0            
234 0           $self->throw("Bad input data (sequences need an id ) ");
235             }
236            
237 0           my $prot_name = $args[0]->display_id;
238 0           return $self->_run($prot_name, @files);
239             }
240              
241              
242             =head2 _run
243              
244             Title : _run
245             Usage : Internal function, not to be called directly
246             Function: makes actual system call to a pseudowise program
247             Example :
248             Returns : nothing; pseudowise output is written to a
249             temporary file $TMPOUTFILE
250             Args : Name of a files containing 3 sequences in the order of peptide, cdna and genomic
251              
252             =cut
253              
254             sub _run {
255 0     0     my ($self,$prot_name, @files) = @_;
256 0           my $instring;
257 0           $self->debug( "Program ".$self->executable."\n");
258 0           my ($tfh1,$outfile) = $self->io->tempfile(-dir=>$self->tempdir);
259 0           my $paramstring = $self->_setparams;
260 0           my $commandstring = sprintf("%s %s %s > %s",
261             $self->executable,
262             $paramstring,
263             join(" ", @files),
264             $outfile);
265            
266 0 0 0       if($self->silent || $self->quiet || ($self->verbose < 1)){
      0        
267 0 0         my $null = ($^O =~ m/mswin/i) ? 'NUL' : '/dev/null';
268 0           $commandstring .= " 2> $null";
269             }
270 0           $self->debug( "pseudowise command = $commandstring\n");
271             # my $status = system($commandstring);
272 0           `$commandstring`;
273             # $self->throw( "Pseudowise call ($commandstring) crashed: $? \n")
274             # unless $status == 0;
275             #parse the outpur and return a Bio::Seqfeature array
276 0           my $genes = $self->_parse_results($prot_name,$outfile);
277 0           close($tfh1);
278 0           undef $tfh1;
279 0 0         if( $self->verbose > 0 ) {
280 0 0         open($tfh1,$outfile) || die $!;
281 0           while(<$tfh1>) {
282 0           $self->debug ($_);
283             }
284             }
285 0           return @{$genes};
  0            
286             }
287              
288             =head2 _parse_results
289              
290             Title : __parse_results
291             Usage : Internal function, not to be called directly
292             Function: Parses pseudowise output
293             Example :
294             Returns : an reference to an array of Seqfeatures
295             Args : the name of the output file
296              
297             =cut
298              
299             sub _parse_results {
300 0     0     my ($self,$prot_name,$outfile) = @_;
301 0 0         $outfile||$self->throw("No outfile specified");
302 0           my $filehandle;
303 0 0         if (ref ($outfile) !~ /GLOB/i ) {
304 0 0         open ($filehandle, "<".$outfile)
305             or $self->throw ("Couldn't open file ".$outfile.": $!\n");
306             } else {
307 0           $filehandle = $outfile;
308             }
309            
310 0           my @genes;
311             #The big parsing loop - parses exons and predicted peptides
312 0           my $parser = Bio::Tools::Pseudowise->new(-verbose => $self->verbose,
313             -fh => $filehandle);
314 0           while( my $f = $parser->next_feature ) {
315 0           push @genes, $f;
316             }
317 0           return \@genes;
318             }
319              
320             =head2 _setinput()
321              
322             Title : _setinput
323             Usage : Internal function, not to be called directly
324             Function: Create input files for pseudowise program
325             Example :
326             Returns : name of file containing dba data input
327             Args : Seq objects in the order of query protein and cdna and target genomic sequence
328              
329             =cut
330              
331             sub _setinput {
332 0     0     my ($self, $seq1, $seq2, $seq3) = @_;
333 0           my ($tfh1,$tfh2,$tfh3,$outfile1,$outfile2,$outfile3);
334              
335 0 0 0       if(!($seq1->isa("Bio::PrimarySeqI") && $seq2->isa("Bio::PrimarySeqI") &&
      0        
336             $seq2->isa("Bio::PrimarySeqI")))
337 0           { $self->throw("One or more of the sequences are nor Bio::PrimarySeqI objects\n"); }
338 0           my $tempdir = $self->tempdir();
339 0           ($tfh1,$outfile1) = $self->io->tempfile(-dir=>$tempdir);
340 0           ($tfh2,$outfile2) = $self->io->tempfile(-dir=>$tempdir);
341 0           ($tfh3,$outfile3) = $self->io->tempfile(-dir=>$tempdir);
342              
343 0           my $out1 = Bio::SeqIO->new(-fh => $tfh1 ,'-format' => 'Fasta');
344 0           my $out2 = Bio::SeqIO->new(-fh => $tfh2, '-format' => 'Fasta');
345 0           my $out3 = Bio::SeqIO->new(-fh => $tfh3, '-format' => 'Fasta');
346            
347 0           $out1->write_seq($seq1);
348 0           $out2->write_seq($seq2);
349 0           $out3->write_seq($seq3);
350 0           $self->_query_pep_seq($seq1);
351 0           $self->_query_cdna_seq($seq2);
352 0           $self->_subject_dna_seq($seq3);
353              
354 0           close($tfh1);
355 0           close($tfh2);
356 0           close($tfh3);
357 0           undef ($tfh1);
358 0           undef ($tfh2);
359 0           undef ($tfh3);
360 0           return ($outfile1,$outfile2,$outfile3);
361             }
362              
363             sub _setparams {
364 0     0     my ($self) = @_;
365 0           my $param_string;
366 0           foreach my $attr(@PSEUDOWISE_PARAMS){
367 0           my $value = $self->$attr();
368 0 0         next unless (defined $value);
369 0           my $attr_key = ' -'.(lc $attr);
370 0           $param_string .=$attr_key.' '.$value;
371             }
372              
373 0           foreach my $attr(@PSEUDOWISE_SWITCHES){
374 0           my $value = $self->$attr();
375 0 0         next unless (defined $value);
376 0           my $attr_key = ' -'.(lc $attr);
377 0           $param_string .=$attr_key;
378             }
379              
380 0           return $param_string;
381             }
382            
383              
384              
385             =head2 _query_pep_seq()
386              
387             Title : _query_pep_seq
388             Usage : Internal function, not to be called directly
389             Function: get/set for the query sequence
390             Example :
391             Returns :
392             Args :
393              
394             =cut
395              
396             sub _query_pep_seq {
397 0     0     my ($self,$seq) = @_;
398 0 0         if(defined $seq){
399 0           $self->{'_query_pep_seq'} = $seq;
400             }
401 0           return $self->{'_query_pep_seq'};
402             }
403              
404             =head2 _query_cdna_seq()
405              
406             Title : _query_cdna_seq
407             Usage : Internal function, not to be called directly
408             Function: get/set for the query sequence
409             Example :
410             Returns :
411             Args :
412              
413             =cut
414              
415             sub _query_cdna_seq {
416 0     0     my ($self,$seq) = @_;
417 0 0         if(defined $seq){
418 0           $self->{'_query_cdna_seq'} = $seq;
419             }
420 0           return $self->{'_query_cdna_seq'};
421             }
422              
423             =head2 _subject_dna_seq()
424              
425             Title : _subject_dna_seq
426             Usage : Internal function, not to be called directly
427             Function: get/set for the subject sequence
428             Example :
429             Returns :
430              
431             Args :
432              
433             =cut
434              
435             sub _subject_dna_seq {
436 0     0     my ($self,$seq) = @_;
437 0 0         if(defined $seq){
438 0           $self->{'_subject_dna_seq'} = $seq;
439             }
440 0           return $self->{'_subject_dna_seq'};
441             }
442              
443             1; # Needed to keep compiler happy