File Coverage

blib/lib/Bio/Tools/Run/Infernal.pm
Criterion Covered Total %
statement 65 324 20.0
branch 17 178 9.5
condition 4 78 5.1
subroutine 13 32 40.6
pod 21 21 100.0
total 120 633 18.9


line stmt bran cond sub pod time code
1             #
2             # You may distribute this module under the same terms as perl itself
3             #
4             # POD documentation - main docs before the code
5             #
6             # _history
7             #
8             # March 2007 - first full implementation; needs some file IO tweaking between
9             # runs but works for now
10             # April 2008 - add 0.81 parameters (may be removed in the 1.0 release)
11             #
12             # July 2009 - updated for v1.0. No longer supporting pre-1.0 Infernal
13              
14             =head1 NAME
15              
16             Bio::Tools::Run::Infernal - Wrapper for local execution of cmalign, cmbuild,
17             cmsearch, cmscore
18              
19             =head1 SYNOPSIS
20              
21             # parameters which are switches are set with any value that evals TRUE,
22             # others are set to a specific value
23              
24             my $factory = Bio::Tools::Run::Infernal->new(@params);
25              
26             # run cmalign|cmbuild|cmsearch|cmscore|cmemit directly as a wrapper method
27             # this resets the program flag if previously set
28              
29             $factory->cmsearch(@seqs); # searches Bio::PrimarySeqI's based on set cov. model
30             # saves output to optional outfile_name, returns
31             # Bio::SearchIO
32              
33             # only values which are allowed for a program are set, so one can use the same
34             # wrapper for the following...
35              
36             $factory->cmalign(@seqs); # aligns Bio::PrimarySeqI's to a set cov. model,
37             # --merge option allows two alignments generated
38             # from the same CM to be merged.
39             # output to outfile_name, returns Bio::AlignIO
40             $factory->cmscore(); # scores set cov. model against Bio::PrimarySeqI,
41             # output to outfile_name/STDOUT.
42             $factory->cmbuild($aln); # builds covariance model based on alignment
43             # CM to outfile_name or model_file (one is required
44             # here), output to STDOUT.
45             $factory->cmemit(); # emits sequence from specified cov. model;
46             # set one if no file specified. output to
47             # outfile_name, returns Bio::SeqIO or (if -a is set)
48             # Bio::AlignIO
49             $factory->cmcalibrate($file); # calibrates specified cov. model; output to
50             # STDOUT
51             $factory->cmstat($file); # summary stats for cov. model; set one if no file
52             # specified; output to STDOUT
53              
54             # run based on the setting of the program parameter
55              
56             my $factory = Bio::Tools::Run::Infernal->new(-program => 'cmsearch',
57             @params);
58             my $search = $factory->run($seq);
59              
60             # using cmsearch returns a Bio::SearchIO object
61              
62             while (my $result = $searchio->next_result){
63             while(my $hit = $result->next_hit){
64             while (my $hsp = $hit->next_hsp){
65             print join("\t", ( $r->query_name,
66             $hit->name,
67             $hsp->hit->start,
68             $hsp->hit->end,
69             $hsp->meta,
70             $hsp->score,
71             )), "\n";
72             }
73             }
74             }
75              
76             =head1 DESCRIPTION
77              
78             Wrapper module for Sean Eddy's Infernal suite of programs. The current
79             implementation runs cmsearch, cmcalibrate, cmalign, cmemit, cmbuild, cmscore,
80             and cmstat. cmsearch will return a Bio::SearchIO, cmemit a Bio::SeqIO/AlignIO,
81             and cmalign a Bio::AlignIO. All others send output to STDOUT. Optionally,
82             any program's output can be redirected to outfile_name.
83              
84             We HIGHLY suggest upgrading to Infernal 1.0. In that spirit, this wrapper now
85             supports parameters for Infernal 1.0 only; for wrapping older versions of
86             Infernal we suggest using the version of Bio::Tools::Run::Infernal that came
87             with previous versions of BioPerl-run.
88              
89             NOTE: Due to conflicts in the way Infernal parameters are now formatted vs.
90             subroutine naming in Perl (specifically the inclusion of hyphens) and due to the
91             very large number of parameters available, setting and resetting parameters via
92             set_parameters() and reset_parameters() is required. All valid parameters can
93             be set, but only ones valid for the executable set via program()/program_name()
94             are used for calling the executables, the others are silently ignored.
95              
96             =head1 FEEDBACK
97              
98             =head2 Mailing Lists
99              
100             User feedback is an integral part of the evolution of this and other
101             Bioperl modules. Send your comments and suggestions preferably to one
102             of the Bioperl mailing lists. Your participation is much appreciated.
103              
104             bioperl-l@bioperl.org - General discussion
105             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
106              
107             =head2 Support
108              
109             Please direct usage questions or support issues to the mailing list:
110              
111             I
112              
113             rather than to the module maintainer directly. Many experienced and
114             reponsive experts will be able look at the problem and quickly
115             address it. Please include a thorough description of the problem
116             with code and data examples if at all possible.
117              
118             =head2 Reporting Bugs
119              
120             Report bugs to the Bioperl bug tracking system to help us keep track
121             the bugs and their resolution. Bug reports can be submitted via the
122             web:
123              
124             http://redmine.open-bio.org/projects/bioperl/
125              
126             =head1 AUTHOR - Chris Fields
127              
128             Email: cjfields-at-uiuc-dot-edu
129              
130             =head1 CONTRIBUTORS
131              
132             cjfields-at-uiuc-dot-edu
133              
134             =head1 APPENDIX
135              
136             The rest of the documentation details each of the object
137             methods. Internal methods are usually preceded with a _
138              
139             =cut
140              
141             package Bio::Tools::Run::Infernal;
142              
143 1     1   105253 use strict;
  1         2  
  1         26  
144 1     1   3 use warnings;
  1         1  
  1         30  
145 1     1   4 use base qw(Bio::Root::Root Bio::Tools::Run::WrapperBase Bio::ParameterBaseI);
  1         1  
  1         450  
146              
147 1     1   664 use Bio::SeqIO;
  1         21618  
  1         31  
148 1     1   491 use Bio::SearchIO;
  1         8029  
  1         40  
149 1     1   423 use Bio::AlignIO;
  1         55126  
  1         42  
150 1     1   7 use Data::Dumper;
  1         2  
  1         4733  
151              
152             # yes, these are the current parameters
153             our %INFERNAL_PARAMS = (
154             'A' => ['switch', '-', qw(cmbuild)],
155             'E' => ['param', '-', qw(cmsearch cmstat)],
156             'F' => ['switch', '-', qw(cmbuild)],
157             'Lmax' => ['param', '--', qw(cmscore)],
158             'Lmin' => ['param', '--', qw(cmscore)],
159             'T' => ['param', '-', qw(cmsearch cmstat)],
160             'Wbeta' => ['param', '--', qw(cmbuild)],
161             'Z' => ['param', '-', qw(cmsearch cmstat)],
162             'a' => ['switch', '-', qw(cmbuild cmemit cmscore)],
163             'afile' => ['param', '--', qw(cmstat)],
164             'ahmm' => ['param', '--', qw(cmemit)],
165             'all' => ['switch', '--', qw(cmstat)],
166             'aln-hbanded' => ['switch', '--', qw(cmsearch)],
167             'aln-optacc' => ['switch', '--', qw(cmsearch)],
168             'aln2bands' => ['switch', '--', qw(cmscore cmsearch)],
169             'banddump' => ['param', '--', qw(cmalign)],
170             'begin' => ['param', '--', qw(cmemit)],
171             'beta' => ['param', '--', qw(cmalign cmscore cmsearch cmstat)],
172             'betae' => ['param', '--', qw(cmscore)],
173             'betas' => ['param', '--', qw(cmscore)],
174             'bfile' => ['param', '--', qw(cmstat)],
175             'binary' => ['switch', '--', qw(cmbuild)],
176             'bits' => ['switch', '--', qw(cmstat)],
177             'bottomonly' => ['switch', '--', qw(cmsearch)],
178             'c' => ['switch', '-', qw(cmemit)],
179             'call' => ['switch', '--', qw(cmbuild)],
180             'cdump' => ['param', '--', qw(cmbuild)],
181             'cfile' => ['param', '--', qw(cmbuild)],
182             'checkfb' => ['switch', '--', qw(cmalign)],
183             'checkpost' => ['switch', '--', qw(cmalign)],
184             'cmL' => ['param', '--', qw(cmstat)],
185             'cmaxid' => ['param', '--', qw(cmbuild)],
186             'cmtbl' => ['param', '--', qw(cmbuild)],
187             'corig' => ['switch', '--', qw(cmbuild)],
188             'ctarget' => ['param', '--', qw(cmbuild)],
189             'cyk' => ['switch', '--', qw(cmalign cmbuild cmsearch)],
190             'devhelp' => ['switch', '--', qw(cmalign cmbuild cmcalibrate cmemit cmscore cmsearch)],
191             'dlev' => ['param', '--', qw(cmalign)],
192             'dna' => ['switch', '--', qw(cmalign cmemit cmsearch)],
193             'eX' => ['param', '--', qw(cmbuild)],
194             'eent' => ['switch', '--', qw(cmbuild)],
195             'efile' => ['param', '--', qw(cmstat)],
196             'ehmmre' => ['param', '--', qw(cmbuild)],
197             'elself' => ['param', '--', qw(cmbuild)],
198             'emap' => ['param', '--', qw(cmbuild)],
199             'emit' => ['switch', '--', qw(cmscore)],
200             'end' => ['param', '--', qw(cmemit)],
201             'enone' => ['switch', '--', qw(cmbuild)],
202             'ere' => ['param', '--', qw(cmbuild)],
203             'exp' => ['param', '--', qw(cmemit)],
204             'exp-T' => ['param', '--', qw(cmcalibrate)],
205             'exp-beta' => ['param', '--', qw(cmcalibrate)],
206             'exp-cmL-glc' => ['param', '--', qw(cmcalibrate)],
207             'exp-cmL-loc' => ['param', '--', qw(cmcalibrate)],
208             'exp-ffile' => ['param', '--', qw(cmcalibrate)],
209             'exp-fract' => ['param', '--', qw(cmcalibrate)],
210             'exp-gc' => ['param', '--', qw(cmcalibrate)],
211             'exp-hfile' => ['param', '--', qw(cmcalibrate)],
212             'exp-hmmLn-glc' => ['param', '--', qw(cmcalibrate)],
213             'exp-hmmLn-loc' => ['param', '--', qw(cmcalibrate)],
214             'exp-hmmLx' => ['param', '--', qw(cmcalibrate)],
215             'exp-no-qdb' => ['switch', '--', qw(cmcalibrate)],
216             'exp-pfile' => ['param', '--', qw(cmcalibrate)],
217             'exp-qqfile' => ['param', '--', qw(cmcalibrate)],
218             'exp-random' => ['switch', '--', qw(cmcalibrate)],
219             'exp-sfile' => ['param', '--', qw(cmcalibrate)],
220             'exp-tailn-cglc' => ['param', '--', qw(cmcalibrate)],
221             'exp-tailn-cloc' => ['param', '--', qw(cmcalibrate)],
222             'exp-tailn-hglc' => ['param', '--', qw(cmcalibrate)],
223             'exp-tailn-hloc' => ['param', '--', qw(cmcalibrate)],
224             'exp-tailp' => ['param', '--', qw(cmcalibrate)],
225             'exp-tailxn' => ['param', '--', qw(cmcalibrate)],
226             'fil-E-hmm' => ['param', '--', qw(cmsearch)],
227             'fil-E-qdb' => ['param', '--', qw(cmsearch)],
228             'fil-F' => ['param', '--', qw(cmcalibrate)],
229             'fil-N' => ['param', '--', qw(cmcalibrate)],
230             'fil-Smax-hmm' => ['param', '--', qw(cmsearch)],
231             'fil-T-hmm' => ['param', '--', qw(cmsearch)],
232             'fil-T-qdb' => ['param', '--', qw(cmsearch)],
233             'fil-aln2bands' => ['switch', '--', qw(cmcalibrate)],
234             'fil-beta' => ['param', '--', qw(cmsearch)],
235             'fil-dfile' => ['param', '--', qw(cmcalibrate)],
236             'fil-gemit' => ['switch', '--', qw(cmcalibrate)],
237             'fil-no-hmm' => ['switch', '--', qw(cmsearch)],
238             'fil-no-qdb' => ['switch', '--', qw(cmsearch)],
239             'fil-nonbanded' => ['switch', '--', qw(cmcalibrate)],
240             'fil-tau' => ['param', '--', qw(cmcalibrate)],
241             'fil-xhmm' => ['param', '--', qw(cmcalibrate)],
242             'fins' => ['switch', '--', qw(cmalign cmbuild)],
243             'forecast' => ['param', '--', qw(cmcalibrate cmsearch)],
244             'forward' => ['switch', '--', qw(cmscore cmsearch)],
245             'g' => ['switch', '-', qw(cmsearch cmstat)],
246             'ga' => ['switch', '--', qw(cmsearch cmstat)],
247             'gapthresh' => ['param', '--', qw(cmalign cmbuild)],
248             'gcfile' => ['param', '--', qw(cmsearch)],
249             'ge' => ['switch', '--', qw(cmstat)],
250             'gfc' => ['switch', '--', qw(cmstat)],
251             'gfi' => ['switch', '--', qw(cmstat)],
252             'gibbs' => ['switch', '--', qw(cmbuild)],
253             'gtbl' => ['param', '--', qw(cmbuild)],
254             'gtree' => ['param', '--', qw(cmbuild)],
255             'h' => ['switch', '-', qw(cmalign cmbuild cmcalibrate cmemit cmscore cmsearch cmstat)],
256             'hbanded' => ['switch', '--', qw(cmalign cmscore cmsearch)],
257             'hmm-W' => ['param', '--', qw(cmsearch)],
258             'hmm-cW' => ['param', '--', qw(cmsearch)],
259             'hmmL' => ['param', '--', qw(cmstat)],
260             'hsafe' => ['switch', '--', qw(cmalign cmscore)],
261             'ignorant' => ['switch', '--', qw(cmbuild)],
262             'iins' => ['switch', '--', qw(cmbuild)],
263             'infile' => ['param', '--', qw(cmscore)],
264             'informat' => ['param', '--', qw(cmalign cmbuild cmsearch)],
265             'inside' => ['switch', '--', qw(cmalign cmscore cmsearch)],
266             'l' => ['switch', '-', qw(cmalign cmbuild cmemit cmscore)],
267             'lambda' => ['param', '--', qw(cmsearch)],
268             'le' => ['switch', '--', qw(cmstat)],
269             'lfc' => ['switch', '--', qw(cmstat)],
270             'lfi' => ['switch', '--', qw(cmstat)],
271             'm' => ['switch', '-', qw(cmstat)],
272             'matchonly' => ['switch', '--', qw(cmalign)],
273             'merge' => ['switch', '--', qw(cmalign)],
274             'mpi' => ['switch', '--', qw(cmalign cmcalibrate cmscore cmsearch)],
275             'mxsize' => ['param', '--', qw(cmalign cmbuild cmcalibrate cmscore cmsearch)],
276             'n' => ['param', '-', qw(cmbuild cmemit cmscore)],
277             'nc' => ['switch', '--', qw(cmsearch cmstat)],
278             'no-null3' => ['switch', '--', qw(cmalign cmcalibrate cmscore cmsearch)],
279             'no-qdb' => ['switch', '--', qw(cmsearch)],
280             'noalign' => ['switch', '--', qw(cmsearch)],
281             'nobalance' => ['switch', '--', qw(cmbuild)],
282             'nodetach' => ['switch', '--', qw(cmbuild)],
283             'nonbanded' => ['switch', '--', qw(cmalign cmbuild cmscore)],
284             'null' => ['param', '--', qw(cmbuild)],
285             'null2' => ['switch', '--', qw(cmsearch)],
286             'o' => ['param', '-', qw(cmalign cmsearch)],
287             'old' => ['switch', '--', qw(cmscore)],
288             'onepost' => ['switch', '--', qw(cmalign)],
289             'optacc' => ['switch', '--', qw(cmalign)],
290             'outfile' => ['param', '--', qw(cmscore)],
291             'p' => ['switch', '-', qw(cmalign cmsearch)],
292             'pad' => ['switch', '--', qw(cmscore)],
293             'pbegin' => ['param', '--', qw(cmalign cmcalibrate cmemit cmscore cmsearch)],
294             'pbswitch' => ['param', '--', qw(cmbuild)],
295             'pebegin' => ['switch', '--', qw(cmalign cmcalibrate cmemit cmscore cmsearch)],
296             'pend' => ['param', '--', qw(cmalign cmcalibrate cmemit cmscore cmsearch)],
297             'pfend' => ['param', '--', qw(cmalign cmcalibrate cmemit cmscore cmsearch)],
298             'prior' => ['param', '--', qw(cmbuild)],
299             'q' => ['switch', '-', qw(cmalign)],
300             'qdb' => ['switch', '--', qw(cmalign cmscore)],
301             'qdbboth' => ['switch', '--', qw(cmscore)],
302             'qdbfile' => ['param', '--', qw(cmstat)],
303             'qdbsmall' => ['switch', '--', qw(cmscore)],
304             'random' => ['switch', '--', qw(cmscore)],
305             'rdump' => ['param', '--', qw(cmbuild)],
306             'refine' => ['param', '--', qw(cmbuild)],
307             'regress' => ['param', '--', qw(cmalign cmbuild cmscore)],
308             'resonly' => ['switch', '--', qw(cmalign)],
309             'rf' => ['switch', '--', qw(cmalign cmbuild)],
310             'rna' => ['switch', '--', qw(cmalign cmemit cmsearch)],
311             'rsearch' => ['param', '--', qw(cmbuild)],
312             'rtrans' => ['switch', '--', qw(cmsearch)],
313             's' => ['param', '-', qw(cmalign cmbuild cmcalibrate cmemit cmscore)],
314             'sample' => ['switch', '--', qw(cmalign)],
315             'scoreonly' => ['switch', '--', qw(cmscore)],
316             'search' => ['switch', '--', qw(cmscore cmstat)],
317             'seqfile' => ['param', '--', qw(cmstat)],
318             'sfile' => ['param', '--', qw(cmstat)],
319             'shmm' => ['param', '--', qw(cmemit)],
320             'small' => ['switch', '--', qw(cmalign)],
321             'stall' => ['switch', '--', qw(cmalign cmscore cmsearch)],
322             'sub' => ['switch', '--', qw(cmalign cmbuild cmscore)],
323             'sums' => ['switch', '--', qw(cmalign cmsearch)],
324             'tabfile' => ['param', '--', qw(cmsearch)],
325             'tau' => ['param', '--', qw(cmalign cmbuild cmscore cmsearch)],
326             'taue' => ['param', '--', qw(cmscore)],
327             'taus' => ['param', '--', qw(cmscore)],
328             'tc' => ['switch', '--', qw(cmsearch cmstat)],
329             'tfile' => ['param', '--', qw(cmalign cmbuild cmemit cmscore)],
330             'toponly' => ['switch', '--', qw(cmsearch cmstat)],
331             'u' => ['switch', '-', qw(cmemit)],
332             'v' => ['switch', '-', qw(cmbuild cmcalibrate)],
333             'viterbi' => ['switch', '--', qw(cmalign cmscore cmsearch)],
334             'wblosum' => ['switch', '--', qw(cmbuild)],
335             'wgiven' => ['switch', '--', qw(cmbuild)],
336             'wgsc' => ['switch', '--', qw(cmbuild)],
337             'wid' => ['param', '--', qw(cmbuild)],
338             'withali' => ['param', '--', qw(cmalign)],
339             'withpknots' => ['switch', '--', qw(cmalign)],
340             'wnone' => ['switch', '--', qw(cmbuild)],
341             'wpb' => ['switch', '--', qw(cmbuild)],
342             'x' => ['switch', '-', qw(cmsearch)],
343             'xfile' => ['param', '--', qw(cmstat)],
344             );
345              
346             our %INFERNAL_PROGRAM = (
347             'cmalign' => "cmalign [-options] \n".
348             'cmalign [-options] --merge ',
349             'cmbuild' => 'cmbuild [-options] ',
350             'cmcalibrate' => 'cmcalibrate [-options] ',
351             'cmemit' => 'cmemit [-options] ',
352             'cmscore' => 'cmscore [-options] ',
353             'cmsearch' => 'cmsearch [-options] ',
354             'cmstat' => 'cmstat [-options] ',
355             );
356              
357             # this is a simple lookup for easy validation for passed methods
358             our %LOCAL_PARAMS = map {$_ => 1} qw(program outfile tempfile model);
359              
360             =head2 new
361              
362             Title : new
363             Usage : my $wrapper = Bio::Tools::Run::Infernal->new(@params)
364             Function: creates a new Infernal factory
365             Returns: Bio::Tools::Run::Infernal wrapper
366             Args : list of parameters
367              
368             =cut
369              
370             sub new {
371 1     1 1 77 my ($class,@args) = @_;
372 1         9 my $self = $class->SUPER::new(@args);
373            
374             # these are specific parameters we do not want passed on to set_parameters
375 1         39 my ($program, $model, $validate, $q, $o1, $o2) =
376             $self->_rearrange([qw(PROGRAM
377             MODEL_FILE
378             VALIDATE_PARAMETERS
379             QUIET
380             OUTFILE_NAME
381             O)], @args);
382            
383 1 50 33     28 if ($o1 && $o2) {
384 0         0 $self->warn("Only assign to either -outfile_name or -o, not both;");
385             }
386 1   33     7 my $out = $o1 || $o2;
387 1         3 $self->validate_parameters($validate);
388 1 50       2 $q && $self->quiet($q);
389 1 50       3 $program && $self->program($program);
390 1 50       3 $model && $self->model_file($model);
391 1   50     4 $out ||= '';
392 1         7 $self->outfile_name($out);
393 1         6 $self->io->_initialize_io();
394 1         21 $self->set_parameters(@args);
395 1         2 return $self;
396             }
397              
398             =head2 program
399              
400             Title : program
401             Usage : $obj->program()
402             Function: Set the program called when run() is used. Synonym of
403             program_name()
404             Returns : String (program name)
405             Args : String (program name)
406             Status : Unstable (may delegate to program_name, which is the interface method)
407              
408             =cut
409              
410             sub program {
411 2     2 1 3 my $self = shift;
412 2         3 return $self->program_name(@_);
413             }
414              
415             =head2 program_name
416              
417             Title : program_name
418             Usage : $factory>program_name()
419             Function: holds the program name
420             Returns: string
421             Args : None
422              
423             =cut
424              
425             sub program_name {
426 8     8 1 6 my ($self) = shift;
427 8 100       13 if (@_) {
428 1         1 my $p = shift;
429             $self->throw("Program '$p' not supported")
430 1 50       15 if !exists $INFERNAL_PROGRAM{lc $p};
431 1         2 $self->{'_program'} = lc $p;
432             # set up cache of valid parameters
433 1         6 while (my ($p, $data) = each %INFERNAL_PARAMS) {
434 190         116 my %in_exe = map {$_ => 1} @$data[2..$#{$data}];
  290         376  
  190         186  
435 190 100       469 $self->{valid_params}->{$p} = 1 if exists $in_exe{$self->{'_program'}};
436             }
437             }
438 8         27 return $self->{'_program'};
439             }
440              
441             =head2 model_file
442              
443             Title : model_file
444             Usage : $obj->model_file()
445             Function: Set the model file used when run() is called.
446             Returns : String (file location of covariance model)
447             Args : String (file location of covariance model)
448              
449             =cut
450              
451             sub model_file {
452 0     0 1 0 my $self = shift;
453 0 0       0 return $self->{'_model_file'} = shift if @_;
454 0         0 return $self->{'_model_file'};
455             }
456              
457             =head2 program_dir
458              
459             Title : program_dir
460             Usage : $factory->program_dir(@params)
461             Function: returns the program directory, obtained from ENV variable.
462             Returns: string
463             Args :
464              
465             =cut
466              
467             sub program_dir {
468 3     3 1 3 my ($self, $dir) = @_;
469 3 50       7 if ($dir) {
470 0         0 $self->{_program_dir} = $dir;
471             }
472 3   50     9 return Bio::Root::IO->catfile($ENV{INFERNALDIR}) || '';
473             }
474              
475             =head2 version
476              
477             Title : version
478             Usage : $v = $prog->version();
479             Function: Determine the version number of the program (uses cmsearch)
480             Example :
481             Returns : float or undef
482             Args : none
483              
484             =cut
485              
486             sub version {
487 0     0 1 0 my ($self) = @_;
488 0 0       0 return unless $self->executable;
489 0         0 my $exe = $self->executable;
490 0         0 my $string = `$exe -h 2>&1`;
491 0         0 my $v;
492 0 0       0 if ($string =~ m{Infernal\s([\d.]+)}) {
493 0         0 $v = $1;
494 0 0       0 $self->deprecated(-message => "Only Infernal 1.0 and above is supported.",
495             -version => 1.006001) if $v < 1;
496             }
497 0   0     0 return $self->{'_progversion'} = $v || undef;
498             }
499              
500             =head2 run
501              
502             Title : run
503             Usage : $obj->run($seqFile)
504             Function: Runs Infernal and returns Bio::SearchIO
505             Returns : A Bio::SearchIO
506             Args : A Bio::PrimarySeqI or file name
507              
508             =cut
509              
510             # TODO: update to accept multiple seqs, alignments
511             sub run {
512 0     0 1 0 my ($self,@seq) = @_;
513 0 0 0     0 if (ref $seq[0] && $seq[0]->isa("Bio::PrimarySeqI") ){# it is an object
    0 0        
514 0         0 my $infile1 = $self->_writeSeqFile(@seq);
515 0         0 return $self->_run($infile1);
516             } elsif (ref $seq[0] && $seq[0]->isa("Bio::Align::AlignI") ){ # it is an object
517 0         0 my $infile1 = $self->_writeAlignFile(@seq);
518 0         0 return $self->_run($infile1);
519             } else {
520 0         0 return $self->_run(@seq);
521             }
522             }
523              
524             =head1 Specific program interface methods
525              
526             =head2 cmsearch
527              
528             Title : cmsearch
529             Usage : $obj->cmsearch($seqFile)
530             Function: Runs Infernal cmsearch and returns Bio::SearchIO
531             Returns : A Bio::SearchIO
532             Args : Bio::PrimarySeqI or file name
533              
534             =cut
535              
536             sub cmsearch {
537 0     0 1 0 my ($self,@seq) = @_;
538 0         0 $self->program('cmsearch');
539 0 0 0     0 if (ref $seq[0] && $seq[0]->isa("Bio::PrimarySeqI") ){# it is an object
540 0         0 my $infile1 = $self->_writeSeqFile(@seq);
541 0         0 return $self->_run(-seq_files => [$infile1]);
542             } else {
543 0         0 return $self->_run(-seq_files => \@seq);
544             }
545             }
546              
547             =head2 cmalign
548              
549             Title : cmalign
550             Usage : $obj->cmalign($seqFile)
551             Function: Runs Infernal cmalign and returns Bio::AlignIO
552             Returns : A Bio::AlignIO
553             Args : Bio::PrimarySeqI or file name
554              
555             =cut
556              
557             sub cmalign {
558 0     0 1 0 my ($self,@seq) = @_;
559 0         0 $self->program('cmalign');
560 0 0       0 if (ref $seq[0]) { # it is an object
561 0 0       0 if ($seq[0]->isa("Bio::PrimarySeqI") ){
    0          
562 0         0 my $infile1 = $self->_writeSeqFile(@seq);
563 0         0 return $self->_run(-seq_files => [$infile1]);
564             } elsif ( $seq[0]->isa("Bio::Align::AlignI") ) {
565 0 0       0 if (scalar(@seq) != 2) {
566 0         0 $self->throw("")
567             }
568 0         0 my $infile1 = $self->_writeAlignFile($seq[0]);
569 0         0 my $infile2 = $self->_writeAlignFile($seq[1]);
570 0         0 return $self->_run(-align_files => [$infile1, $infile2]);
571             }
572             } else {
573             # we can maybe add a check for the file extension and try to DTRT
574 0         0 my %params = $self->get_parameters('valid');
575 0 0       0 $params{merge} ? return $self->_run(-align_files => \@seq):
576             return $self->_run(-seq_files => \@seq);
577 0         0 return $self->_run(-seq_files => \@seq);
578             }
579             }
580              
581             =head2 cmemit
582              
583             Title : cmemit
584             Usage : $obj->cmemit($modelfile)
585             Function: Runs Infernal cmemit and returns Bio::AlignIO
586             Returns : A Bio::AlignIO
587             Args : None; set model_file() to use a specific model
588              
589             =cut
590              
591             sub cmemit {
592 0     0 1 0 my ($self) = shift;
593 0         0 $self->program('cmemit');
594 0         0 return $self->_run(@_);
595             }
596              
597             =head2 cmbuild
598              
599             Title : cmbuild
600             Usage : $obj->cmbuild($alignment)
601             Function: Runs Infernal cmbuild and saves covariance model
602             Returns : 1 on success (no object for covariance models)
603             Args : Bio::AlignIO with structural information (such as from Stockholm
604             format source) or alignment file name
605              
606             =cut
607              
608             sub cmbuild {
609 0     0 1 0 my ($self,@seq) = @_;
610 0         0 $self->program('cmbuild');
611 0 0 0     0 if (ref $seq[0] && $seq[0]->isa("Bio::Align::AlignI") ){# it is an object
612 0         0 my $infile1 = $self->_writeAlignFile(@seq);
613 0         0 return $self->_run(-align_files => [$infile1]);
614             } else {
615 0         0 return $self->_run(-align_files => \@seq);
616             }
617             }
618              
619             =head2 cmscore
620              
621             Title : cmscore
622             Usage : $obj->cmscore($seq)
623             Function: Runs Infernal cmscore and saves output
624             Returns : None
625             Args : None; set model_file() to use a specific model
626              
627             =cut
628              
629             sub cmscore {
630 0     0 1 0 my ($self,@seq) = @_;
631 0         0 $self->program('cmscore');
632 0         0 return $self->_run();
633             }
634              
635             =head2 cmcalibrate
636              
637             Title : cmcalibrate
638             Usage : $obj->cmcalibrate('file')
639             Function: Runs Infernal calibrate on specified CM
640             Returns : None
641             Args : None; set model_file() to use a specific model
642              
643             =cut
644              
645             sub cmcalibrate {
646 0     0 1 0 my ($self,@seq) = @_;
647 0         0 $self->program('cmcalibrate');
648 0         0 return $self->_run();
649             }
650              
651             =head2 cmstat
652              
653             Title : cmstat
654             Usage : $obj->cmstat($seq)
655             Function: Runs Infernal cmstat and saves output
656             Returns : None
657             Args : None; set model_file() to use a specific model
658              
659             =cut
660              
661             sub cmstat {
662 0     0 1 0 my ($self,@seq) = @_;
663 0         0 $self->program('cmstat');
664 0         0 return $self->_run();
665             }
666              
667             =head1 Bio::ParameterBaseI-specific methods
668              
669             These methods are part of the Bio::ParameterBaseI interface
670              
671             =cut
672              
673             =head2 set_parameters
674              
675             Title : set_parameters
676             Usage : $pobj->set_parameters(%params);
677             Function: sets the parameters listed in the hash or array
678             Returns : None
679             Args : [optional] hash or array of parameter/values. These can optionally
680             be hash or array references
681             Note : This only sets parameters; to set methods use the method name
682             =cut
683              
684             sub set_parameters {
685 1     1 1 1 my $self = shift;
686             # circumvent any issues arising from passing in refs
687 0         0 my %args = (ref($_[0]) eq 'HASH') ? %{$_[0]} :
688 1 50       6 (ref($_[0]) eq 'ARRAY') ? @{$_[0]} :
  0 50       0  
689             @_;
690             # set the parameters passed in, but only ones supported for the program
691 1         2 my ($prog, $validate) = ($self->program, $self->validate_parameters);
692              
693             # parameter cleanup
694 1         4 %args = map { my $a = $_;
  1         1  
695 1         3 $a =~ s{^-}{};
696 1         3 lc $a => $args{$_}
697             } sort keys %args;
698            
699 1         4 while (my ($key, $val) = each %args) {
700 1 50       3 if (exists $INFERNAL_PARAMS{$key}) {
701 0         0 my ($type, $prefix) = @{$INFERNAL_PARAMS{$key}}[0..1];
  0         0  
702 0         0 @{$self->{parameters}->{$key}} = ($type, $prefix);
  0         0  
703 0 0 0     0 unshift @{$self->{parameters}->{$key}},
  0 0       0  
704             $type eq 'param' ? $val :
705             $type eq 'switch' && $val ? 1 : 0;
706 0 0       0 if ($validate) {
707 0         0 my %in_exe = map {$_ => 1} @{$INFERNAL_PARAMS{$key}}[2..$#{$INFERNAL_PARAMS{$key}}];
  0         0  
  0         0  
  0         0  
708 0 0       0 $self->warn("Parameter $key not used for $prog") if !exists $in_exe{$key};
709             }
710             } else {
711 1 50       4 $self->warn("Parameter $key does not exist") if ($validate);
712             }
713             }
714             }
715              
716             =head2 reset_parameters
717              
718             Title : reset_parameters
719             Usage : resets values
720             Function: resets parameters to either undef or value in passed hash
721             Returns : none
722             Args : [optional] hash of parameter-value pairs
723              
724             =cut
725              
726             sub reset_parameters {
727 0     0 1 0 my $self = shift;
728 0         0 delete $self->{parameters};
729 0 0       0 if (@_) {
730 0         0 $self->set_parameters(@_);
731             }
732             }
733              
734             =head2 validate_parameters
735              
736             Title : validate_parameters
737             Usage : $pobj->validate_parameters(1);
738             Function: sets a flag indicating whether to validate parameters via
739             set_parameters() or reset_parameters()
740             Returns : Bool
741             Args : [optional] value evaluating to True/False
742             Note : Optionally implemented method; up to the implementation on whether
743             to automatically validate parameters or optionally do so
744              
745             =cut
746              
747             sub validate_parameters {
748 2     2 1 3 my ($self) = shift;
749 2 100       3 if (@_) {
750 1 50       4 $self->{validate_params} = defined $_[0] ? 1 : 0;
751             }
752 2         2 return $self->{validate_params};
753             }
754              
755             =head2 parameters_changed
756              
757             Title : parameters_changed
758             Usage : if ($pobj->parameters_changed) {...}
759             Function: Returns boolean true (1) if parameters have changed
760             Returns : Boolean (0 or 1)
761             Args : None
762             Note : This module does not run state checks, so this always returns True
763              
764             =cut
765              
766 0     0 1   sub parameters_changed { 1 }
767              
768             =head2 available_parameters
769              
770             Title : available_parameters
771             Usage : @params = $pobj->available_parameters()
772             Function: Returns a list of the available parameters
773             Returns : Array of parameters
774             Args : [optional] name of executable being used; defaults to returning all
775             available parameters
776              
777             =cut
778              
779             sub available_parameters {
780 0     0 1   my ($self, $exec) = @_;
781 0           my @params;
782 0 0         if ($exec) {
783 0 0         $self->throw("$exec is not part of the Infernal package") if !exists($INFERNAL_PROGRAM{$exec});
784 0           for my $p (sort keys %INFERNAL_PARAMS) {
785 0 0         if (grep { $exec eq $_ }
  0            
786 0           @{$INFERNAL_PARAMS{$p}}[2..$#{$INFERNAL_PARAMS{$p}}])
  0            
787             {
788 0           push @params, $p;
789             }
790             }
791             } else {
792 0           @params = (sort keys %INFERNAL_PARAMS, sort keys %LOCAL_PARAMS);
793             }
794 0           return @params;
795             }
796              
797             =head2 get_parameters
798              
799             Title : get_parameters
800             Usage : %params = $pobj->get_parameters;
801             Function: Returns list of set key-value pairs, parameter => value
802             Returns : List of key-value pairs
803             Args : [optional]
804             'full' - this option returns everything associated with the parameter
805             as an array ref value; that is, not just the value but also
806             the value, type, and prefix. Default is value only.
807             'valid'- same a 'full', but only returns the grouping valid for the
808             currently set executable
809              
810             =cut
811              
812             sub get_parameters {
813 0     0 1   my ($self, $option) = @_;
814 0   0       $option ||= ''; # no option
815 0           my %params;
816 0 0         if (exists $self->{parameters}) {
817             %params = (ref $option eq 'ARRAY') ? (
818 0           map {$_ => $self->{parameters}{$_}[0]}
819 0           grep { exists $self->{parameters}{$_} } @$option) :
820             (lc $option eq 'full') ?
821 0           (%{$self->{parameters}}) :
822             (lc $option eq 'valid') ?
823 0           (map {$_ => $self->{parameters}{$_}}
824 0           grep { exists $self->{valid_params}->{$_} } keys %{$self->{parameters}}) :
  0            
825 0 0         (map {$_ => $self->{parameters}{$_}[0]} keys %{$self->{parameters}});
  0 0          
  0 0          
826             } else {
827 0           %params = ();
828             }
829 0           return %params;
830             }
831              
832             =head1 to_* methods
833              
834             All to_* methods are implementation-specific
835              
836             =cut
837              
838             =head2 to_exe_string
839              
840             Title : to_exe_string
841             Usage : $string = $pobj->to_exe_string;
842             Function: Returns string (command line string in this case)
843             Returns : String
844             Args :
845              
846             =cut
847              
848             sub to_exe_string {
849 0     0 1   my ($self, @passed) = @_;
850 0           my ($seqs, $aligns) = $self->_rearrange([qw(SEQ_FILES ALIGN_FILES)], @passed);
851 0 0 0       if ($seqs || $aligns) {
852 0 0 0       $self->throw("Seqs or alignments must be an array reference") unless
      0        
      0        
853             ($seqs && ref($seqs) eq 'ARRAY') || ($aligns && ref($aligns) eq 'ARRAY' );
854             }
855            
856 0           my %args = map {$_ => []} qw(switch param input redirect);
  0            
857 0           my %params = $self->get_parameters('valid');
858            
859 0           my ($exe, $prog, $model, $outfile) = ($self->executable,
860             $self->program_name,
861             $self->model_file,
862             $self->outfile_name);
863            
864 0 0         $self->throw("Executable not found") unless defined($exe);
865            
866 0 0         delete $params{o} if exists $params{o};
867              
868 0 0 0       if (!defined($model) && $prog ne 'cmbuild') {
869 0           $self->throw("model_file() not defined")
870             }
871            
872 0   0       $outfile ||= '';
873            
874 0           for my $p (sort keys %params) {
875 0 0         if ($params{$p}[0]) {
876 0 0         my $val = $params{$p}[1] eq 'param' ? ' '.$params{$p}[0] : '';
877 0           push @{$args{$params{$p}[1]}}, $params{$p}[2].$p.$val;
  0            
878             }
879             }
880            
881             # TODO: not sure what happens when we pass in multiple seq or alignment
882             # filenames, may need checking
883 0 0 0       if ($prog eq 'cmscore' || $prog eq 'cmstat' || $prog eq 'cmcalibrate') {
    0 0        
    0          
    0          
    0          
884 0 0         push @{$args{'redirect'}}, "> $outfile" if $outfile;
  0            
885 0           push @{$args{'input'}}, $model;
  0            
886             } elsif ($prog eq 'cmsearch') {
887 0 0         if (!defined $seqs) {
888 0           $self->throw('cmsearch requires a sequence file name');
889             }
890 0 0         push @{$args{'param'}}, "-o $outfile" if $outfile;
  0            
891 0           push @{$args{'input'}}, ($model, @$seqs);
  0            
892             } elsif ($prog eq 'cmalign') {
893 0 0         if ($params{'merge'}) {
894 0 0 0       $self->throw('cmalign with --merge option requires two alignment files') if
895             !defined($aligns) || @$aligns < 2;
896 0           push @{$args{'input'}}, ($model, @$aligns);
  0            
897             } else {
898 0 0         $self->throw('cmalign requires a sequence file') if
899             !defined $seqs;
900 0           push @{$args{'input'}}, ($model, @$seqs);
  0            
901             }
902 0 0         push @{$args{'param'}}, "-o $outfile" if $outfile;
  0            
903             } elsif ($prog eq 'cmbuild') {
904 0 0         $self->throw('cmbuild requires one alignment file') if
905             !defined($aligns);
906 0 0         if ($model) {
907 0           push @{$args{'input'}}, ($model, @$aligns);
  0            
908 0 0         push @{$args{'redirect'}}, "> $outfile" if $outfile;
  0            
909             } else {
910 0           push @{$args{'input'}}, ($outfile, @$aligns);
  0            
911             }
912             } elsif ($prog eq 'cmemit') {
913 0 0         if (!$outfile) {
914 0           $self->throw('cmemit requires an outfile_name; tempfile support not implemented yet');
915             } else {
916 0           push @{$args{'input'}}, ($model, ,$outfile);
  0            
917             }
918             }
919            
920             # quiet!
921 0 0 0       if ($self->quiet && $prog ne 'cmsearch') {
922 0 0         if ($prog eq 'cmalign') {
923 0 0         push @{$args{switch}}, '-q' if !exists $params{q};
  0            
924             } else {
925 0 0         my $null = ($^O =~ m/mswin/i) ? 'NUL' : '/dev/null';
926 0           push @{$args{redirect}}, "> $null";
  0            
927             }
928             }
929            
930 0           my $string = "$exe ".join(' ',(@{$args{switch}},
931 0           @{$args{param}},
932 0           @{$args{input}},
933 0           @{$args{redirect}}));
  0            
934              
935 0           $string;
936             }
937              
938             ############### PRIVATE ###############
939              
940             #=head2 _run
941             #
942             # Title : _run
943             # Usage : $obj->_run()
944             # Function: Internal(not to be used directly)
945             # Returns :
946             # Args :
947             #
948             #=cut
949              
950             {
951             my %ALLOWED = map {$_ => 1} qw(run cmsearch cmalign cmemit cmbuild
952             cmcalibrate cmstat cmscore);
953              
954             sub _run {
955 0     0     my ($self)= shift;
956            
957 0           my ($prog, $model, $out, $version) = ($self->program,
958             $self->model_file,
959             $self->outfile_name,
960             $self->version);
961            
962 0 0         if (my $caller = (caller(1))[3]) {
963 0           $caller =~ s{.*::(\w+)$}{$1};
964 0 0         $self->throw("Calling _run() from disallowed method") unless exists $ALLOWED{$caller};
965             } else {
966 0           $self->throw("Can't call _run directly");
967             }
968            
969             # a model and a file must be defined for all but cmemit; cmemit must have a
970             # file or model defined (using $file if both are defined)
971            
972             # relevant files are passed on to the string builder
973 0           my $str = $self->to_exe_string(@_);
974 0           $self->debug("Infernal command: $str\n");
975            
976 0           my %has = $self->get_parameters('valid');
977            
978             my $obj =
979             ($prog eq 'cmsearch') ? Bio::SearchIO->new(-format => 'infernal',
980             -version => $version,
981             -model => $model) :
982             ($prog eq 'cmalign' ) ? Bio::AlignIO->new(-format => 'stockholm') :
983 0 0 0       ($prog eq 'cmemit' && $has{a}) ? Bio::AlignIO->new(-format => 'stockholm') :
    0          
    0          
    0          
984             ($prog eq 'cmemit') ? Bio::SeqIO->new(-format => 'fasta') :
985             undef;
986 0           my @args;
987             # file output
988 0 0         if ($out) {
989 0           my $status = system($str);
990 0 0 0       if($status || !-e $out || -z $out ) {
      0        
991 0 0         my $error = ($!) ? "$! Status: $status" : "Status: $status";
992 0           $self->throw( "Infernal call crashed: $error \n[command $str]\n");
993 0           return undef;
994             }
995 0 0 0       if ($obj && ref($obj)) {
996 0           $obj->file($out);
997 0           @args = (-file => $out);
998             }
999             # fh-based (no outfile)
1000             } else {
1001 0 0         open(my $fh,"$str |") || $self->throw("Infernal call ($str) crashed: $?\n");
1002 0 0 0       if ($obj && ref($obj)) {
1003 0           $obj->fh($fh);
1004 0           @args = (-fh => $fh);
1005             } else {
1006             # dump to debugging
1007 0           my $io;
1008 0           while(<$fh>) {$io .= $_;}
  0            
1009 0           close($fh);
1010 0 0         $self->debug($io) if $io;
1011 0           return 1;
1012             }
1013             }
1014 0 0 0       $obj->_initialize_io(@args) if $obj && ref($obj);
1015 0   0       return $obj || 1;
1016             }
1017              
1018             }
1019              
1020             =head2 _writeSeqFile
1021              
1022             Title : _writeSeqFile
1023             Usage : obj->_writeSeqFile($seq)
1024             Function: Internal(not to be used directly)
1025             Returns :
1026             Args :
1027              
1028             =cut
1029              
1030             sub _writeSeqFile {
1031 0     0     my ($self,@seq) = @_;
1032 0           my ($tfh,$inputfile) = $self->io->tempfile(-dir=>$self->tempdir);
1033 0           my $in = Bio::SeqIO->new(-fh => $tfh , '-format' => 'Fasta');
1034 0           foreach my $s(@seq){
1035 0           $in->write_seq($s);
1036             }
1037 0           $in->close();
1038 0           $in = undef;
1039 0           close($tfh);
1040 0           undef $tfh;
1041 0           return $inputfile;
1042             }
1043              
1044             =head2 _writeAlignFile
1045              
1046             Title : _writeAlignFile
1047             Usage : obj->_writeAlignFile($seq)
1048             Function: Internal(not to be used directly)
1049             Returns :
1050             Args :
1051              
1052             =cut
1053              
1054             sub _writeAlignFile{
1055 0     0     my ($self,@align) = @_;
1056 0           my ($tfh,$inputfile) = $self->io->tempfile(-dir=>$self->tempdir);
1057 0           my $in = Bio::AlignIO->new('-fh' => $tfh ,
1058             '-format' => 'stockholm');
1059 0           foreach my $s(@align){
1060 0           $in->write_aln($s);
1061             }
1062 0           $in->close();
1063 0           $in = undef;
1064 0           close($tfh);
1065 0           undef $tfh;
1066 0           return $inputfile;
1067             }
1068              
1069             # this is a private sub used to regenerate the class data structures,
1070             # dumped to STDOUT
1071              
1072             # could probably add in a description field if needed...
1073              
1074             sub _dump_params {
1075 0     0     my %params;
1076             my %usage;
1077 0           for my $exec (qw(cmalign cmbuild cmcalibrate cmemit cmscore cmsearch cmstat)) {
1078 0           my $output = `$exec --devhelp`;
1079 0 0         if ($?) {
1080 0           $output = `$exec -h`;
1081             }
1082 0           my @lines = split("\n",$output);
1083            
1084 0           for my $line (@lines) {
1085 0 0         next if $line =~ /^#/;
1086 0 0         if ($line =~ /^\s*(-{1,2})(\S+)\s+(<\S+>)?/) {
    0          
1087 0           my %data;
1088 0 0         ($data{prefix}, my $p, $data{arg}) = ($1, $2, $3 ? 'param' : 'switch');
1089 0 0         if (exists $params{$p}) {
1090 0 0         if ($data{prefix} ne $params{$p}{prefix}) {
1091 0           warn("$data{prefix} for $p in $exec doesn't match prefix for same parameter in ".$params{$p}{exec}[-1].":".$params{$p}{prefix});
1092             }
1093 0 0         if ($data{arg} ne $params{$p}{arg}) {
1094 0           warn("$data{arg} for $p in $exec doesn't match arg for same parameter in ".$params{$p}{exec}[-1].":".$params{$p}{arg});
1095             }
1096             }
1097            
1098 0           while (my ($key, $val) = each %data) {
1099 0           $params{$p}->{$key} = $val;
1100             }
1101 0           push @{$params{$p}->{exec}}, $exec;
  0            
1102             } elsif ($line =~ /Usage:\s*(.+)$/) {
1103 0           push @{$usage{$exec}}, $1;
  0            
1104             } else {
1105             #print "$line\n";
1106             }
1107             }
1108             }
1109              
1110             # generate data structure
1111 0           print "our %INFERNAL_PARAMS = (\n";
1112 0           for my $k (sort keys %params) {
1113 0           printf(" %-17s => [","'$k'");
1114 0           for my $sub (qw(arg prefix exec)) {
1115             my $str = (ref($params{$k}{$sub}) eq 'ARRAY') ?
1116 0           "qw(".join(' ', @{$params{$k}{$sub}}).")" :
1117 0 0         "'".$params{$k}{$sub}."',";
1118 0           printf("%-10s", $str);
1119             }
1120 0           print "],\n";
1121             }
1122 0           print ");\n\n";
1123            
1124             # generate usage data structure
1125 0           print "our %INFERNAL_PROGRAM = (\n";
1126 0           for my $k (sort keys %usage) {
1127 0           printf(" %-17s => [\n","'$k'");
1128 0           print ' '.join(",\n ", map {"'$_'"} @{$usage{$k}})."\n";
  0            
  0            
1129 0           print " ],\n";
1130             }
1131 0           print ");\n";
1132              
1133             }
1134              
1135             1;