File Coverage

blib/lib/Bio/Tools/Run/Primate.pm
Criterion Covered Total %
statement 48 159 30.1
branch 6 54 11.1
condition 0 22 0.0
subroutine 13 22 59.0
pod 6 6 100.0
total 73 263 27.7


line stmt bran cond sub pod time code
1             # BioPerl module for Bio::Tools::Run::Primate
2             #
3             # Please direct questions and support issues to
4             #
5             # Cared for by
6             #
7             # Copyright Shawn Hoon
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             Wrapper for Primate, Guy Slater's near exact match finder for short sequence
15             tags.
16              
17             =head1 SYNOPSIS
18              
19             use Bio::Tools::Run::Primate;
20             use Bio::SeqIO;
21              
22             my $query = "primer.fa";
23             my $target = "contig.fa";
24              
25             my @params = ("query" => $query,"target" => $target,"m"=>0);
26             my $fact = Bio::Tools::Run::Primate->new(@params);
27              
28             my @feat = $fact->run;
29             foreach my $feat(@feat) {
30             print $feat->seqname."\t".$feat->primary_tag."\t".$feat->start.
31             "\t".$feat->end."\t".$feat->strand."\t".$feat->seq->seq."\n";
32             }
33              
34             =head1 DESCRIPTION
35              
36             Primate is available under to ensembl-nci package at
37             http://cvsweb.sanger.ac.uk/cgi-bin/cvsweb.cgi/ensembl-nci/?cvsroot=Ensembl
38              
39             =head1 FEEDBACK
40              
41             =head2 Mailing Lists
42              
43             User feedback is an integral part of the evolution of this and other
44             Bioperl modules. Send your comments and suggestions preferably to one
45             of the Bioperl mailing lists. Your participation is much appreciated.
46              
47             bioperl-l@bioperl.org - General discussion
48             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
49              
50             =head2 Support
51              
52             Please direct usage questions or support issues to the mailing list:
53              
54             I
55              
56             rather than to the module maintainer directly. Many experienced and
57             reponsive experts will be able look at the problem and quickly
58             address it. Please include a thorough description of the problem
59             with code and data examples if at all possible.
60              
61             =head2 Reporting Bugs
62              
63             Report bugs to the Bioperl bug tracking system to help us keep track
64             the bugs and their resolution. Bug reports can be submitted via the
65             web:
66              
67             http://redmine.open-bio.org/projects/bioperl/
68              
69             =head1 AUTHOR - Shawn Hoon
70              
71             Email shawnh@fugu-sg.org
72              
73             =head1 APPENDIX
74              
75             The rest of the documentation details each of the object
76             methods. Internal methods are usually preceded with a _
77              
78             =cut
79              
80              
81             package Bio::Tools::Run::Primate;
82 1         67 use vars qw($AUTOLOAD @ISA $PROGRAM $PROGRAMDIR @PRIMATE_PARAMS $PROGRAMNAME
83 1     1   103732 @OTHER_SWITCHES %OK_FIELD);
  1         1  
84 1     1   4 use strict;
  1         1  
  1         14  
85 1     1   470 use Bio::Root::Root;
  1         17545  
  1         23  
86 1     1   4 use Bio::Root::IO;
  1         1  
  1         13  
87 1     1   377 use Bio::Factory::ApplicationFactoryI;
  1         101  
  1         17  
88 1     1   433 use Bio::SeqIO;
  1         20793  
  1         24  
89 1     1   501 use Bio::SeqFeature::Generic;
  1         40477  
  1         27  
90 1     1   407 use Bio::Tools::Run::WrapperBase;
  1         2  
  1         60  
91              
92              
93             @ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase);
94              
95              
96             BEGIN {
97              
98 1     1   3 @PRIMATE_PARAMS = qw(V Q T M B QUERY TARGET OUTFILE PROGRAM EXECUTABLE);
99 1         2 @OTHER_SWITCHES = qw(QUIET VERBOSE);
100              
101             # Authorize attribute fields
102 1         2 foreach my $attr ( @PRIMATE_PARAMS,@OTHER_SWITCHES) { $OK_FIELD{$attr}++; }
  12         1117  
103             }
104              
105             =head2 program_name
106              
107             Title : program_name
108             Usage : $factory>program_name()
109             Function: holds the program name
110             Returns: string
111             Args : None
112              
113             =cut
114              
115             sub program_name {
116 6     6 1 25 return 'primate';
117             }
118              
119             =head2 program_dir
120              
121             Title : program_dir
122             Usage : $factory->program_dir(@params)
123             Function: returns the program directory, obtained from ENV variable.
124             Returns: string
125             Args :
126              
127             =cut
128              
129             sub program_dir {
130 3 50   3 1 55 return Bio::Root::IO->catfile($ENV{PRIMATEDIR}) if $ENV{PRIMATEDIR};
131             }
132              
133             =head2 new
134              
135             Title : new
136             Usage : my $obj = Bio::Tools::Run::Primate->new()
137             Function: Builds a new Bio::Tools::Run::Primate objet
138             Returns : Bio::Tools::Run::Primate
139             Args : query => the L object or a file path
140             target => the L object or a file path
141             m => the number of mismatches allowed, default 1(integer)
142             b => [TRUE|FALSE] find best match, default FALSE
143             executable=>where the program sits
144              
145             =cut
146              
147             sub new {
148 1     1 1 99 my ($class, @args) = @_;
149 1         9 my $self = $class->SUPER::new(@args);
150            
151 1         33 my ($attr, $value);
152              
153 1         3 while (@args) {
154 4         5 $attr = shift @args;
155 4         3 $value = shift @args;
156 4 50       8 next if( $attr =~ /^-/ ); # don't want named parameters
157 4 50       10 if($attr =~/^q$/i){
158 0         0 $self->query($value);
159             }
160 4 50       7 if($attr =~/^t$/i){
161 0         0 $self->target($value);
162             }
163 4         19 $self->$attr($value);
164             }
165 1         2 return $self;
166             }
167              
168             sub AUTOLOAD {
169 4     4   2 my $self = shift;
170 4         4 my $attr = $AUTOLOAD;
171 4         9 $attr =~ s/.*:://;
172 4         5 $attr = uc $attr;
173 4 50       7 $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
174 4 50       8 $self->{$attr} = shift if @_;
175 4         18 return $self->{$attr};
176             }
177              
178             =head2 version
179              
180             Title : version
181             Usage : $primate->version
182             Function: Determine the version number of the program
183             Returns : float or undef
184             Args : none
185              
186             =cut
187              
188             sub version {
189 0     0 1   my ($self) = @_;
190              
191 0           my $exe = $self->executable();
192 0 0         return undef unless defined $exe;
193 0           my $string = `$exe -v ` ;
194 0           $string =~ /\(([\d.]+)\)/;
195 0   0       return $1 || undef;
196             }
197              
198             =head2 search
199              
200             Title : search
201             Usage : DEPRECATED. Use $factory->run() instead
202             Function: Perform a primate search
203             Returns : Array of L
204             Args :
205              
206             =cut
207              
208             sub search {
209 0     0 1   return shift->run(@_);
210             }
211              
212              
213             =head2 run
214              
215             Title : run
216             Usage : @feat = $factory->run();
217             Function: Perform a primate search
218             Returns : Array of L
219             Args :
220              
221             =cut
222              
223             sub run{
224 0     0 1   my ($self,$target) = @_;
225 0   0       $target = $target ||$self->target;
226 0 0         $target || $self->throw("Need a target sequence");
227 0 0         $self->query || $self->throw("Need a query sequence");
228              
229             # Create input file pointer
230 0           my ($query_file,$target_file)= $self->_setinput($self->query,$target);
231 0 0 0       if (!($query_file && $target_file)) {$self->throw("Unable to create temp files for query and target !");}
  0            
232              
233             # Create parameter string to pass to primate program
234 0           my $param_string = $self->_setparams();
235              
236             # run primate
237 0           my @feats= $self->_run($query_file,$target_file,$param_string);
238 0           return @feats;
239             }
240              
241             #################################################
242             #INTERNAL METHODS
243              
244             =head2 _run
245              
246             Title : _run
247             Usage : Internal function, not to be called directly
248             Function: makes actual system call to dba program
249             Returns : array of L
250             Args : path to query and target file and parameter string
251              
252             =cut
253              
254             sub _run {
255 0     0     my ($self,$query_file,$target_file,$param_string) = @_;
256 0           my $instring;
257 0           $self->debug( "Program ".$self->executable."\n");
258 0           my ($tfh,$outfile) = $self->io->tempfile(-dir=>$self->tempdir);
259 0           close($tfh); # this is to make sure we don't have
260             # open filehandles
261 0           undef $tfh;
262 0           my $commandstring = $self->executable.
263             " $param_string -q $query_file -t $target_file > $outfile";
264 0           $self->debug( "primate command = $commandstring");
265 0           my $status = system($commandstring);
266 0 0         $self->throw( "primate call ($commandstring) crashed: $? \n") unless $status==0;
267              
268             #parse pff format and return a Bio::Search::HSP::GenericHSP array
269 0           my @feats = $self->_parse_results($outfile);
270              
271 0           return @feats;
272             }
273              
274             =head2 _parse_results
275              
276             Title : _parse_results
277             Usage : Internal function, not to be called directly
278             Function: Passes primate output
279             Returns : array of L
280             Args : the name of the output file
281              
282             =cut
283              
284             sub _parse_results {
285 0     0     my ($self,$outfile) = @_;
286 0 0         $outfile||$self->throw("No outfile specified");
287 0           my @feats;
288 0           my %query = $self->_query_seq();
289              
290 0           open(OUT,$outfile);
291 0           while(my $entry = ){
292 0           chomp($entry);
293 0 0         if($entry =~ /primate/ ) {
294 0           my ($dummy,$tagname, $seqname, $strand,$seq_end,$mismatch) = split(" " , $entry );
295             #map primate coordinates to Seq coordinates
296 0           my $seq_start = $seq_end- length($query{$tagname})+2;
297 0           $seq_end++;
298 0           my $feature = Bio::SeqFeature::Generic->new( -seq_id => $seqname,
299             -strand => $strand,
300             -score => $mismatch,
301             -start => $seq_start,
302             -end => $seq_end,
303             -frame => 1,
304             -source => 'primate',
305             -primary => $tagname);
306 0           $feature->attach_seq($self->_target_seq);
307 0           push @feats,$feature;
308             }
309             }
310              
311 0           return @feats;
312             }
313              
314              
315             =head2 _setinput()
316              
317             Title : _setinput
318             Usage : Internal function, not to be called directly
319             Function: Create input files for primate
320             Returns : name of file containing query and target
321             Args : query and target (either a filename or a L
322              
323             =cut
324              
325             sub _setinput {
326 0     0     my ($self, $query,$target) = @_;
327 0           my ($query_file,$target_file,$tfh1,$tfh2);
328              
329 0 0         my @query = ref ($query) eq "ARRAY" ? @{$query} : ($query);
  0            
330 0           foreach my $query(@query){
331              
332 0 0 0       if(ref($query)&& $query->isa("Bio::PrimarySeqI")){
    0          
333 0           ($tfh1,$query_file) = $self->io->tempfile(-dir=>$self->tempdir);
334 0           my $out1 = Bio::SeqIO->new(-fh=> $tfh1 , '-format' => 'fasta');
335 0           my %query;
336 0           $query{$query->primary_id} = $query->seq;
337 0           $self->_query_seq(\%query);
338 0 0         $out1->write_seq($query) || return 0;
339 0           close ($tfh1);
340 0           undef $tfh1;
341             }
342             elsif (-e $query){
343 0           my $in = Bio::SeqIO->new(-file => $query , '-format' => 'fasta');
344 0           ($tfh1,$query_file) = $self->io->tempfile(-dir=>$self->tempdir);
345 0           my $out1 = Bio::SeqIO->new(-fh=> $tfh1 , '-format' => 'fasta');
346 0           my %query;
347 0           while(my $seq1 = $in->next_seq()){
348 0 0         $out1->write_seq($seq1) || return 0;
349 0           $query{$seq1->primary_id} = $seq1->seq;
350             }
351 0           close($tfh1);
352 0           undef $tfh1;
353 0           $self->_query_seq(\%query);
354             }
355             else {
356 0           return 0;
357             }
358             }
359 0 0 0       if(ref($target) && $target->isa("Bio::PrimarySeqI")){
    0          
360 0           ($tfh2,$target_file) = $self->io->tempfile(-dir=>$self->tempdir);
361 0           my $out1 = Bio::SeqIO->new(-fh=> $tfh2 , '-format' => 'Fasta');
362 0 0         $out1->write_seq($target)|| return 0;
363 0           $self->_target_seq($target);
364 0           close($tfh2);
365 0           undef $tfh2;
366             }
367             elsif (-e $target){
368 0           my $in = Bio::SeqIO->new(-file => $target , '-format' => 'fasta');
369 0           ($tfh2,$target_file) = $self->io->tempfile(-dir=>$self->tempdir);
370 0           my $out = Bio::SeqIO->new(-fh=> $tfh2 , '-format' => 'fasta');
371 0   0       my $seq1 = $in->next_seq() || return 0;
372 0           $out->write_seq($seq1);
373 0           close($tfh2);
374 0           undef $tfh2;
375 0           $self->_target_seq($seq1);
376             }
377             else {
378 0           return 0;
379             }
380              
381 0           return $query_file,$target_file;
382             }
383              
384             =head2 _setparams()
385              
386             Title : _setparams
387             Usage : Internal function, not to be called directly
388             Function: Create parameter inputs for primate program
389             Returns : parameter string to be passed to primate
390             Args : the param array
391              
392             =cut
393              
394             sub _setparams {
395 0     0     my ($attr, $value, $self);
396              
397 0           $self = shift;
398              
399 0           my $param_string = "";
400 0           for $attr ( @PRIMATE_PARAMS ) {
401 0           $value = $self->$attr();
402 0 0         next unless (defined $value);
403              
404 0           my $attr_key = lc $attr; #put params in format expected by dba
405 0           $attr_key = ' -'.$attr_key;
406 0 0 0       if(($attr_key !~/QUERY/i) && ($attr_key !~/TARGET/i)){
407 0           $param_string .= $attr_key.' '.$value;
408             }
409             }
410              
411 0 0 0       if ($self->quiet() || $self->verbose() < 0) {
412 0 0         my $null = ($^O =~ m/mswin/i) ? 'NUL' : '/dev/null';
413 0           $param_string .= " >$null ";
414             }
415 0           return $param_string;
416             }
417              
418             =head2 _query_seq()
419              
420             Title : _query_seq
421             Usage : Internal function, not to be called directly
422             Function: get/set for the query sequence
423             Returns : a hash of seq with key the query tag
424             Args : optional
425              
426             =cut
427              
428             sub _query_seq {
429 0     0     my ($self,$seq) = @_;
430 0 0         if(defined $seq){
431 0           $self->{'_query_seq'} = $seq;
432             }
433 0           return %{$self->{'_query_seq'}};
  0            
434             }
435              
436             =head2 _target_seq()
437              
438             Title : _target_seq
439             Usage : Internal function, not to be called directly
440             Function: get/set for the target sequence
441             Returns : L
442             Args : optional
443              
444             =cut
445              
446             sub _target_seq {
447 0     0     my ($self,$seq) = @_;
448 0 0         if(defined $seq){
449 0           $self->{'_target_seq'} = $seq;
450             }
451 0           return $self->{'_target_seq'};
452             }
453              
454             1; # Needed to keep compiler happy