Bio/Search/Hit/HmmpfamHit.pm

   1 #
   2 # BioPerl module for Bio::Search::Hit::HmmpfamHit
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Sendu Bala <bix@sendu.me.uk>
   7 #
   8 # Copyright Sendu Bala
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Search::Hit::HmmpfamHit - A parser and hit object for hmmpfam hits
  17
  18 =head1 SYNOPSIS
  19
  20     # generally we use Bio::SearchIO to build these objects
  21     use Bio::SearchIO;
  22     my $in = Bio::SearchIO->new(-format => 'hmmer_pull',
  23                                                            -file   => 'result.hmmer');
  24
  25     while (my $result = $in->next_result) {
  26                 while (my $hit = $result->next_hit) {
  27                         print $hit->name, "\n";
  28                         print $hit->score, "\n";
  29                         print $hit->significance, "\n";
  30
  31                         while (my $hsp = $hit->next_hsp) {
  32                                 # process HSPI objects
  33                         }
  34                 }
  35     }
  36
  37 =head1 DESCRIPTION
  38
  39 This object implements a parser for hmmpfam hit output, a program in the HMMER
  40 package.
  41
  42 =head1 FEEDBACK
  43
  44 =head2 Mailing Lists
  45
  46 User feedback is an integral part of the evolution of this and other
  47 Bioperl modules. Send your comments and suggestions preferably to
  48 the Bioperl mailing list.  Your participation is much appreciated.
  49
  50   bioperl-l@bioperl.org                  - General discussion
  51   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  52
  53 =head2 Support
  54
  55 Please direct usage questions or support issues to the mailing list:
  56
  57 I<bioperl-l@bioperl.org>
  58
  59 rather than to the module maintainer directly. Many experienced and
  60 reponsive experts will be able look at the problem and quickly
  61 address it. Please include a thorough description of the problem
  62 with code and data examples if at all possible.
  63
  64 =head2 Reporting Bugs
  65
  66 Report bugs to the Bioperl bug tracking system to help us keep track
  67 of the bugs and their resolution. Bug reports can be submitted via the
  68 web:
  69
  70   https://github.com/bioperl/bioperl-live/issues
  71
  72 =head1 AUTHOR - Sendu Bala
  73
  74 Email bix@sendu.me.uk
  75
  76 =head1 APPENDIX
  77
  78 The rest of the documentation details each of the object methods.
  79 Internal methods are usually preceded with a _
  80
  81 =cut
  82
  83 # Let the code begin...
  84
  85 package Bio::Search::Hit::HmmpfamHit;
  86
  87 use strict;
  88
  89 use Bio::Search::HSP::HmmpfamHSP;
  90
  91 use base qw(Bio::Root::Root Bio::Search::Hit::PullHitI);
  92
  93 =head2 new
  94
  95  Title   : new
  96  Usage   : my $obj = Bio::Search::Hit::HmmpfamHit->new();
  97  Function: Builds a new Bio::Search::Hit::HmmpfamHit object.
  98  Returns : Bio::Search::Hit::HmmpfamHit
  99  Args    : -chunk    => [Bio::Root::IO, $start, $end] (required if no -parent)
 100            -parent   => Bio::PullParserI object (required if no -chunk)
 101            -hit_data => array ref with [name description score significance
 102                                                 num_hsps rank]
 103
 104            where the array ref provided to -chunk contains an IO object
 105            for a filehandle to something representing the raw data of the
 106            hit, and $start and $end define the tell() position within the
 107            filehandle that the hit data starts and ends (optional; defaults
 108            to start and end of the entire thing described by the filehandle)
 109
 110 =cut
 111
 112 sub new {
 113     my ($class, @args) = @_;
 114     my $self = $class->SUPER::new(@args);
 115
 116         $self->_setup(@args);
 117
 118         my $fields = $self->_fields;
 119         foreach my $field (qw( next_domain domains hsp_data )) {
 120                 $fields->{$field} = undef;
 121         }
 122
 123         my $hit_data = $self->_raw_hit_data;
 124         if ($hit_data && ref($hit_data) eq 'ARRAY') {
 125                 foreach my $field (qw(name description score significance num_hsps rank)) {
 126                         $fields->{$field} = shift(@{$hit_data});
 127                 }
 128         }
 129         $fields->{hit_start} = 1;
 130
 131         delete $self->_fields->{accession};
 132
 133         $self->_dependencies( { ( length => 'hsp_data' ) } );
 134
 135     return $self;
 136 }
 137
 138 #
 139 # PullParserI discovery methods so we can answer all HitI questions
 140 #
 141
 142 sub _discover_description {
 143         # this should be set when this object is created, but if it was undef as is
 144         # possible, this _discover method will be called: just return and keep the
 145         # return value undef
 146         return;
 147 }
 148
 149 sub _discover_hsp_data {
 150         my $self = shift;
 151         my $hsp_table = $self->get_field('hsp_table');
 152         my $hsp_data = $hsp_table->{$self->get_field('name')} || undef;
 153         if ($hsp_data) {
 154                 if (defined $hsp_data->{hit_length}) {
 155                         $self->_fields->{length} = $hsp_data->{hit_length};
 156                 }
 157
 158                 # rank query_start query_end hit_start hit_end score evalue
 159                 $self->_fields->{hsp_data} = $hsp_data->{hsp_data};
 160         }
 161 }
 162
 163 sub _discover_query_start {
 164         my $self = shift;
 165         my $hsp_data = $self->get_field('hsp_data') || return;
 166
 167         my ($this_hsp) = sort { $a->[1] <=> $b->[1] } @{$hsp_data};
 168         $self->_fields->{query_start} = $this_hsp->[1];
 169 }
 170
 171 sub _discover_query_end {
 172         my $self = shift;
 173         my $hsp_data = $self->get_field('hsp_data') || return;
 174
 175         my ($this_hsp) = sort { $b->[2] <=> $a->[2] } @{$hsp_data};
 176         $self->_fields->{query_end} = $this_hsp->[2];
 177 }
 178
 179 sub _discover_hit_start {
 180         my $self = shift;
 181         my $hsp_data = $self->get_field('hsp_data') || return;
 182
 183         my ($this_hsp) = sort { $a->[3] <=> $b->[3] } @{$hsp_data};
 184         $self->_fields->{hit_start} = $this_hsp->[3];
 185 }
 186
 187 sub _discover_hit_end {
 188         my $self = shift;
 189         my $hsp_data = $self->get_field('hsp_data') || return;
 190
 191         my ($this_hsp) = sort { $b->[4] <=> $a->[4] } @{$hsp_data};
 192         $self->_fields->{hit_end} = $this_hsp->[4];
 193 }
 194
 195 sub _discover_next_hsp {
 196         my $self = shift;
 197         my $hsp_data = $self->get_field('hsp_data') || return;
 198         unless (defined $self->{_next_hsp_index}) {
 199                 $self->{_next_hsp_index} = 0;
 200         }
 201         return if $self->{_next_hsp_index} == -1;
 202
 203         $self->_fields->{next_hsp} = Bio::Search::HSP::HmmpfamHSP->new(-parent => $self,
 204                                                                                                                                   -hsp_data => $hsp_data->[$self->{_next_hsp_index}++]);
 205
 206         if ($self->{_next_hsp_index} > $#{$hsp_data}) {
 207                 $self->{_next_hsp_index} = -1;
 208         }
 209 }
 210
 211 =head2 next_hsp
 212
 213  Title    : next_hsp
 214  Usage    : while( $hsp = $obj->next_hsp()) { ... }
 215  Function : Returns the next available High Scoring Pair
 216  Example  :
 217  Returns  : L<Bio::Search::HSP::HSPI> object or null if finished
 218  Args     : none
 219
 220 =cut
 221
 222 sub next_hsp {
 223     my $self = shift;
 224     my $hsp = $self->get_field('next_hsp');
 225         undef $self->_fields->{next_hsp};
 226         return $hsp;
 227 }
 228
 229 =head2 next_domain
 230
 231  Title   : next_domain
 232  Usage   : my $domain = $hit->next_domain();
 233  Function: An alias for L<next_hsp()>, this will return the next HSP
 234  Returns : L<Bio::Search::HSP::HSPI> object
 235  Args    : none
 236
 237 =cut
 238
 239 *next_domain = \&next_hsp;
 240
 241 =head2 hsps
 242
 243  Usage     : $hit_object->hsps();
 244  Purpose   : Get a list containing all HSP objects.
 245  Example   : @hsps = $hit_object->hsps();
 246  Returns   : list of L<Bio::Search::HSP::BlastHSP> objects.
 247  Argument  : none
 248
 249 =cut
 250
 251 sub hsps {
 252     my $self = shift;
 253         my $old = $self->{_next_hsp_index} || 0;
 254         $self->rewind;
 255         my @hsps;
 256         while (defined(my $hsp = $self->next_hsp)) {
 257                 push(@hsps, $hsp);
 258         }
 259         $self->{_next_hsp_index} =  @hsps > 0 ? $old : -1;
 260         return @hsps;
 261 }
 262
 263 =head2 domains
 264
 265  Title   : domains
 266  Usage   : my @domains = $hit->domains();
 267  Function: An alias for L<hsps()>, this will return the full list of hsps
 268  Returns : array of L<Bio::Search::HSP::HSPI> objects
 269  Args    : none
 270
 271 =cut
 272
 273 *domains = \&hsps;
 274
 275 =head2 hsp
 276
 277  Usage     : $hit_object->hsp( [string] );
 278  Purpose   : Get a single HSPI object for the present HitI object.
 279  Example   : $hspObj  = $hit_object->hsp;  # same as 'best'
 280            : $hspObj  = $hit_object->hsp('best');
 281            : $hspObj  = $hit_object->hsp('worst');
 282  Returns   : Object reference for a L<Bio::Search::HSP::HSPI> object.
 283  Argument  : String (or no argument).
 284            :   No argument (default) = highest scoring HSP (same as 'best').
 285            :   'best'  = highest scoring HSP.
 286            :   'worst' = lowest scoring HSP.
 287  Throws    : Exception if an unrecognized argument is used.
 288
 289 See Also   : L<hsps()|hsps>, L<num_hsps>()
 290
 291 =cut
 292
 293 sub hsp {
 294     my ($self, $type) = @_;
 295         $type ||= 'best';
 296         my $hsp_data = $self->get_field('hsp_data') || return;
 297
 298         my $sort;
 299         if ($type eq 'best') {
 300                 $sort = sub { $a->[6] <=> $b->[6] };
 301         }
 302         elsif ($type eq 'worst') {
 303                 $sort = sub { $b->[6] <=> $a->[6] };
 304         }
 305         else {
 306                 $self->throw("Unknown arg '$type' given to hsp()");
 307         }
 308
 309         my ($this_hsp) = sort $sort @{$hsp_data};
 310         return Bio::Search::HSP::HmmpfamHSP->new(-parent => $self, -hsp_data => $this_hsp);
 311 }
 312
 313 =head2 rewind
 314
 315  Title   : rewind
 316  Usage   : $result->rewind;
 317  Function: Allow one to reset the Hit iterator to the beginning, so that
 318            next_hit() will subsequently return the first hit and so on.
 319  Returns : n/a
 320  Args    : none
 321
 322 =cut
 323
 324 sub rewind {
 325         my $self = shift;
 326         my $hsp_data = $self->get_field('hsp_data') || return;
 327         $self->{_next_hsp_index} = @{$hsp_data} > 0 ? 0 : -1;
 328 }
 329
 330 # have p() a synonym of significance()
 331 sub p {
 332         return shift->significance;
 333 }
 334
 335 =head2 strand
 336
 337  Usage     : $sbjct->strand( [seq_type] );
 338  Purpose   : Gets the strand(s) for the query, sbjct, or both sequences.
 339            : For hmmpfam, the answers are always 1 (forward strand).
 340  Example   : $qstrand = $sbjct->strand('query');
 341            : $sstrand = $sbjct->strand('hit');
 342            : ($qstrand, $sstrand) = $sbjct->strand();
 343  Returns   : scalar context: integer '1'
 344            : array context without args: list of two strings (1, 1)
 345            : Array context can be "induced" by providing an argument of 'list'
 346                    : or 'array'.
 347  Argument  : In scalar context: seq_type = 'query' or 'hit' or 'sbjct' (default
 348            : = 'query') ('sbjct' is synonymous with 'hit')
 349
 350 =cut
 351
 352 sub strand {
 353     my ($self, $type) = @_;
 354         $type ||= (wantarray ? 'list' : 'query');
 355     $type = lc($type);
 356         if ($type eq 'list' || $type eq 'array') {
 357                 return (1, 1);
 358         }
 359         return 1;
 360 }
 361
 362 =head2 frac_aligned_query
 363
 364  Usage     : $hit_object->frac_aligned_query();
 365  Purpose   : Get the fraction of the query sequence which has been aligned
 366            : across all HSPs (not including intervals between non-overlapping
 367            : HSPs).
 368  Example   : $frac_alnq = $hit_object->frac_aligned_query();
 369  Returns   : undef (the length of query sequences is unknown in Hmmpfam reports)
 370  Argument  : none
 371
 372 =cut
 373
 374 # noop
 375 sub frac_aligned_query { }
 376
 377 1;