Bio/Tools/Run/Hmmer.pm

   1 # You may distribute this module under the same terms as perl itself
   2 # POD documentation - main docs before the code
   3
   4 =head1 NAME
   5
   6 Bio::Tools::Run::Hmmer - Wrapper for local execution of hmmsearch
   7 ,hmmbuild, hmmcalibrate, hmmalign, hmmpfam
   8
   9 =head1 SYNOPSIS
  10
  11   #run hmmpfam|hmmalign|hmmsearch
  12   my $factory = Bio::Tools::Run::Hmmer->new('program'=>'hmmsearch','hmm'=>'model.hmm');
  13
  14   # Pass the factory a Bio::Seq object or a file name
  15
  16   # returns a Bio::SearchIO object
  17   my $search = $factory->run($seq);
  18
  19
  20   my @feat;
  21   while (my $result = $searchio->next_result){
  22    while(my $hit = $result->next_hit){
  23     while (my $hsp = $hit->next_hsp){
  24             print join("\t", ( $r->query_name,
  25                                $hsp->query->start,
  26                                $hsp->query->end,
  27                                $hit->name,
  28                                $hsp->hit->start,
  29                                $hsp->hit->end,
  30                                $hsp->score,
  31                                $hsp->evalue,
  32                                $hsp->seq_str,
  33                                )), "\n";
  34     }
  35    }
  36   }
  37
  38   #build a hmm using hmmbuild
  39   my $aio = Bio::AlignIO->new(-file=>"protein.msf",-format=>'msf');
  40   my $aln = $aio->next_aln;
  41   my $factory =  Bio::Tools::Run::Hmmer->new('program'=>'hmmbuild','hmm'=>'model.hmm');
  42   $factory->run($aln);
  43
  44   #calibrate the hmm
  45   my $factory =  Bio::Tools::Run::Hmmer->new('program'=>'hmmcalibrate','hmm'=>'model.hmm');
  46   $factory->run();
  47
  48 =head1 DESCRIPTION
  49
  50 Wrapper module for Sean Eddy's  HMMER suite of program to allow running of hmmsearch,hmmpfam,hmmalign, hmmbuild,hmmconvert. Binaries are available at http://hmmer.wustl.edu/
  51
  52 =head1 FEEDBACK
  53
  54 =head2 Mailing Lists
  55
  56 User feedback is an integral part of the evolution of this and other
  57 Bioperl modules. Send your comments and suggestions preferably to one
  58 of the Bioperl mailing lists.  Your participation is much appreciated.
  59
  60  bioperl-l@bioperl.org          - General discussion
  61  http://bio.perl.org/MailList.html             - About the mailing lists
  62
  63 =head2 Reporting Bugs
  64
  65 Report bugs to the Bioperl bug tracking system to help us keep track
  66 the bugs and their resolution.  Bug reports can be submitted via email
  67 or the web:
  68
  69  bioperl-bugs@bioperl.org
  70  http://bugzilla.bioperl.org/
  71
  72 =head1 AUTHOR - Shawn
  73
  74  Email: shawnh@gmx.net
  75
  76 =head1 CONTRIBUTORS
  77
  78  Shawn Hoon shawnh@gmx.net
  79
  80 =head1 APPENDIX
  81
  82 The rest of the documentation details each of the object
  83 methods. Internal methods are usually preceded with a _
  84
  85 =cut
  86
  87 package Bio::Tools::Run::Hmmer;
  88
  89 use vars qw($AUTOLOAD @ISA @HMMER_PARAMS @HMMER_SWITCHES %OK_FIELD);
  90 use strict;
  91 use Bio::SeqIO;
  92 use Bio::Root::Root;
  93 use Bio::SearchIO;
  94 use Bio::AlignIO;
  95 use Bio::Tools::Run::WrapperBase;
  96
  97 @ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase);
  98
  99 BEGIN {
 100        @HMMER_PARAMS=qw(HMM PROGRAM DB N A E T Z );
 101        @HMMER_SWITCHES=qw(N);
 102        foreach my $attr ( @HMMER_PARAMS,@HMMER_SWITCHES)
 103                         { $OK_FIELD{$attr}++; }
 104 }
 105
 106 =head2 program_name
 107
 108  Title   : program_name
 109  Usage   : $factory>program_name()
 110  Function: holds the program name
 111  Returns:  string
 112  Args    : None
 113
 114 =cut
 115
 116 sub program_name {
 117   my ($self) = shift;
 118   return $self->program(@_);
 119 }
 120
 121 =head2 program_dir
 122
 123  Title   : program_dir
 124  Usage   : $factory->program_dir(@params)
 125  Function: returns the program directory, obtiained from ENV variable.
 126  Returns:  string
 127  Args    :
 128
 129 =cut
 130
 131 sub program_dir {
 132   return Bio::Root::IO->catfile($ENV{HMMERDIR}) if $ENV{HMMERDIR};
 133 }
 134
 135 sub AUTOLOAD {
 136        my $self = shift;
 137        my $attr = $AUTOLOAD;
 138        $attr =~ s/.*:://;
 139        $attr = uc $attr;
 140        $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
 141        $self->{$attr} = shift if @_;
 142        return $self->{$attr};
 143 }
 144
 145 =head2 new
 146
 147  Title   : new
 148  Usage   : $HMMER->new(@params)
 149  Function: creates a new HMMER factory
 150  Returns:  Bio::Tools::Run::HMMER
 151  Args    :
 152
 153 =cut
 154
 155 sub new {
 156        my ($class,@args) = @_;
 157        my $self = $class->SUPER::new(@args);
 158        $self->io->_initialize_io();
 159        my ($attr, $value);
 160        while (@args)  {
 161            $attr =   shift @args;
 162            $value =  shift @args;
 163            next if( $attr =~ /^-/ ); # don't want named parameters
 164            $self->$attr($value);
 165        }
 166        return $self;
 167 }
 168
 169 =head2 run
 170
 171  Title   :   run
 172  Usage   :   $obj->run($seqFile)
 173  Function:   Runs HMMER and returns Bio::SearchIO
 174  Returns :   A Bio::SearchIO
 175  Args    :   A Bio::PrimarySeqI or file name
 176
 177 =cut
 178
 179 sub run{
 180     my ($self,@seq) = @_;
 181
 182     if  (ref $seq[0] && $seq[0]->isa("Bio::PrimarySeqI") ){# it is an object
 183
 184         my $infile1 = $self->_writeSeqFile(@seq);
 185         return  $self->_run($infile1);
 186     }
 187     elsif(ref $seq[0] && $seq[0]->isa("Bio::Align::AlignI")){
 188         my $infile1 = $self->_writeAlignFile(@seq);
 189         return  $self->_run($infile1);
 190     }
 191     else {
 192         return  $self->_run(@seq);
 193     }
 194
 195 }
 196
 197 =head2 _run
 198
 199  Title   :   _run
 200  Usage   :   $obj->_run()
 201  Function:   Internal(not to be used directly)
 202  Returns :   An array of Bio::SeqFeature::Generic objects
 203  Args    :
 204
 205 =cut
 206
 207 sub _run {
 208      my ($self,$file)= @_;
 209
 210      my $str = $self->executable;
 211      my $param_str = $self->_setparams;
 212      $str.=" $param_str ".$file;
 213
 214     $self->debug("HMMER command = $str");
 215
 216      if($self->program_name=~/hmmpfam|hmmsearch|hmmalign/){
 217        open(HMM,"$str |") || $self->throw("HMMER call ($str) crashed: $?\n");
 218        my $searchio= Bio::SearchIO->new(-fh=>\*HMM,-format=>"hmmer");
 219        return $searchio;
 220       }
 221       else { # for hmmbuild or hmmcalibrate
 222         my $status = system($str);
 223         $self->throw("HMMER call($str) crashed: $?\n") unless $status==0;
 224         return 1;
 225       }
 226
 227 }
 228
 229 =head2 _setparams
 230
 231  Title   :  _setparams
 232  Usage   :  Internal function, not to be called directly
 233  Function:  creates a string of params to be used in the command string
 234  Example :
 235  Returns :  string of params
 236  Args    :
 237
 238 =cut
 239
 240 sub _setparams {
 241     my ($self) = @_;
 242     my $param_string;
 243     foreach my $attr(@HMMER_PARAMS){
 244         next if $attr=~/HMM|PROGRAM|DB/i;
 245         my $value = $self->$attr();
 246         next unless (defined $value);
 247         my $attr_key = ' -'.(uc $attr);
 248         $param_string .= $attr_key.' '.$value;
 249     }
 250     foreach my $attr(@HMMER_SWITCHES){
 251         my $value = $self->$attr();
 252         next unless (defined $value);
 253         my $attr_key = ' -'.($attr);
 254         $param_string .=$attr_key;
 255     }
 256     my ($hmm) = $self->HMM || $self->DB || $self->throw("Need to specify either HMM file or Database");
 257     $param_string.=' '.$hmm;
 258
 259     return $param_string;
 260 }
 261
 262
 263 =head2 _writeSeqFile
 264
 265  Title   :   _writeSeqFile
 266  Usage   :   obj->_writeSeqFile($seq)
 267  Function:   Internal(not to be used directly)
 268  Returns :
 269  Args    :
 270
 271 =cut
 272
 273 sub _writeSeqFile{
 274     my ($self,@seq) = @_;
 275     my ($tfh,$inputfile) = $self->io->tempfile(-dir=>$self->tempdir);
 276     my $in  = Bio::SeqIO->new(-fh => $tfh , '-format' => 'Fasta');
 277     foreach my $s(@seq){
 278       $in->write_seq($s);
 279     }
 280     $in->close();
 281     $in = undef;
 282     close($tfh);
 283     undef $tfh;
 284     return $inputfile;
 285
 286 }
 287
 288 =head2 _writeAlignFile
 289
 290  Title   :   _writeAlignFile
 291  Usage   :   obj->_writeAlignFile($seq)
 292  Function:   Internal(not to be used directly)
 293  Returns :
 294  Args    :
 295
 296 =cut
 297
 298 sub _writeAlignFile{
 299     my ($self,@align) = @_;
 300     my ($tfh,$inputfile) = $self->io->tempfile(-dir=>$self->tempdir);
 301     my $in  = Bio::AlignIO->new(-fh => $tfh , '-format' => 'msf');
 302     foreach my $s(@align){
 303       $in->write_aln($s);
 304     }
 305     $in->close();
 306     $in = undef;
 307     close($tfh);
 308     undef $tfh;
 309     return $inputfile;
 310 }
 311 1;