Bio/Search/HSP/PullHSPI.pm

   1 #
   2 # BioPerl module for Bio::Search::HSP::PullHSPI
   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 # You may distribute this module under the same terms as perl itself
   9
  10 # POD documentation - main docs before the code
  11
  12 =head1 NAME
  13
  14 Bio::Search::HSP::PullHSPI - Bio::Search::HSP::HSPI interface for pull parsers.
  15
  16 =head1 SYNOPSIS
  17
  18   # This is an interface and cannot be instantiated
  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           while (my $hsp = $hit->next_hsp) {
  28               $r_type = $hsp->algorithm;
  29               $pvalue = $hsp->p();
  30               $evalue = $hsp->evalue();
  31               $frac_id = $hsp->frac_identical( ['query'|'hit'|'total'] );
  32               $frac_cons = $hsp->frac_conserved( ['query'|'hit'|'total'] );
  33               $gaps = $hsp->gaps( ['query'|'hit'|'total'] );
  34               $qseq = $hsp->query_string;
  35               $hseq = $hsp->hit_string;
  36               $homo_string = $hsp->homology_string;
  37               $len = $hsp->length( ['query'|'hit'|'total'] );
  38               $len = $hsp->length( ['query'|'hit'|'total'] );
  39               $rank = $hsp->rank;
  40           }
  41       }
  42   }
  43
  44
  45 =head1 DESCRIPTION
  46
  47 PullHSP is for fast implementations that only do parsing work on the hsp
  48 data when you actually request information by calling one of the HSPI
  49 methods.
  50
  51 Many methods of HSPI are implemented in a way suitable for inheriting classes
  52 that use Bio::PullParserI. It only really makes sense for PullHSP modules to be
  53 created by (and have as a -parent) PullHit modules.
  54
  55 In addition to the usual -chunk and -parent, -hsp_data is all you should supply
  56 when making a PullHSP object. This will store that data and make it accessible
  57 via _raw_hsp_data, which you can access in your subclass. It would be best to
  58 simply provide the data as the input -chunk instead, if the raw data is large
  59 enough.
  60
  61 =head1 SEE ALSO
  62
  63 This module inherits methods from these other modules:
  64
  65 L<Bio::SeqFeatureI>,
  66 L<Bio::SeqFeature::FeaturePair>
  67 L<Bio::SeqFeature::SimilarityPair>
  68
  69 Please refer to these modules for documentation of the
  70 many additional inherited methods.
  71
  72 =head1 FEEDBACK
  73
  74 =head2 Mailing Lists
  75
  76 User feedback is an integral part of the evolution of this and other
  77 Bioperl modules. Send your comments and suggestions preferably to
  78 the Bioperl mailing list.  Your participation is much appreciated.
  79
  80   bioperl-l@bioperl.org                  - General discussion
  81   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  82
  83 =head2 Support
  84
  85 Please direct usage questions or support issues to the mailing list:
  86
  87 I<bioperl-l@bioperl.org>
  88
  89 rather than to the module maintainer directly. Many experienced and
  90 reponsive experts will be able look at the problem and quickly
  91 address it. Please include a thorough description of the problem
  92 with code and data examples if at all possible.
  93
  94 =head2 Reporting Bugs
  95
  96 Report bugs to the Bioperl bug tracking system to help us keep track
  97 of the bugs and their resolution. Bug reports can be submitted via the
  98 web:
  99
 100   https://github.com/bioperl/bioperl-live/issues
 101
 102 =head1 AUTHOR - Sendu Bala
 103
 104 Email bix@sendu.me.uk
 105
 106 =head1 COPYRIGHT
 107
 108 Copyright (c) 2006 Sendu Bala. All Rights Reserved.
 109
 110 =head1 DISCLAIMER
 111
 112 This software is provided "as is" without warranty of any kind.
 113
 114 =head1 APPENDIX
 115
 116 The rest of the documentation details each of the object methods.
 117 Internal methods are usually preceded with a _
 118
 119 =cut
 120
 121 # Let the code begin...
 122
 123 package Bio::Search::HSP::PullHSPI;
 124
 125
 126 use strict;
 127
 128 use base qw(Bio::Search::HSP::HSPI Bio::PullParserI);
 129
 130 =head2 _setup
 131
 132  Title   : _setup
 133  Usage   : $self->_setup(@args)
 134  Function: Implementers should call this to setup common fields and deal with
 135            common arguments to new().
 136  Returns : n/a
 137  Args    : @args received in new().
 138
 139 =cut
 140
 141 sub _setup {
 142     my ($self, @args) = @_;
 143
 144         # fields most subclasses probably will want
 145         $self->_fields( { ( hsp_length => undef,
 146                         identical => undef,
 147                         percent_identity => undef,
 148                         conserved => undef,
 149                         hsp_gaps => undef,
 150                         query_gaps => undef,
 151                         hit_gaps => undef,
 152                                                 evalue => undef,
 153                                                 pvalue => undef,
 154                                                 score => undef,
 155                                                 query_start => undef,
 156                                                 query_end => undef,
 157                                                 query_string => undef,
 158                                                 hit_start => undef,
 159                                                 hit_end => undef,
 160                                                 hit_string => undef,
 161                                                 homology_string => undef,
 162                                                 rank => undef,
 163                         seq_inds => undef,
 164                         hit_identical_inds => undef,
 165                         hit_conserved_inds => undef,
 166                         hit_nomatch_inds => undef,
 167                         hit_gap_inds => undef,
 168                         query_identical_inds => undef,
 169                         query_conserved_inds => undef,
 170                         query_nomatch_inds => undef,
 171                         query_gap_inds => undef ) } );
 172
 173         my ($parent, $chunk, $hsp_data) = $self->_rearrange([qw(PARENT
 174                                                                                                                     CHUNK
 175                                                                                                                         HSP_DATA)], @args);
 176
 177     $self->throw("Need -parent or -chunk to be defined") unless defined $parent || $chunk;
 178
 179         $self->parent($parent) if $parent;
 180
 181     if ($chunk) {
 182         my ($io, $start, $end) = (undef, 0, undef);
 183         if (ref($chunk) eq 'ARRAY') {
 184             ($io, $start, $end) = @{$chunk};
 185         }
 186         else {
 187             $io = $chunk;
 188         }
 189         $self->chunk($io, -start => $start, -end => $end);
 190     }
 191
 192         $self->_raw_hsp_data($hsp_data) if $hsp_data;
 193
 194     return $self;
 195 }
 196
 197 sub _raw_hsp_data {
 198         my $self = shift;
 199         if (@_) {
 200                 $self->{_raw_hsp_data} = shift;
 201         }
 202         return $self->{_raw_hsp_data};
 203 }
 204
 205 #
 206 # Some of these methods are written explitely to avoid HSPI throwing not
 207 # implemented or the wrong ancestor class being used to answer the method;
 208 # if it didn't do that then PullParserI AUTOLOAD would have cought them.
 209 #
 210
 211 =head2 algorithm
 212
 213  Title   : algorithm
 214  Usage   : my $r_type = $hsp->algorithm
 215  Function: Obtain the name of the algorithm used to obtain the HSP
 216  Returns : string (e.g., BLASTP)
 217  Args    : none
 218
 219 =cut
 220
 221 sub algorithm {
 222         return shift->get_field('algorithm');
 223 }
 224
 225 =head2 pvalue
 226
 227  Title   : pvalue
 228  Usage   : my $pvalue = $hsp->pvalue();
 229  Function: Returns the P-value for this HSP or undef
 230  Returns : float or exponential (2e-10)
 231  Args    : none
 232
 233 =cut
 234
 235 sub pvalue {
 236         return shift->get_field('pvalue');
 237 }
 238
 239 =head2 evalue
 240
 241  Title   : evalue
 242  Usage   : my $evalue = $hsp->evalue();
 243  Function: Returns the e-value for this HSP
 244  Returns : float or exponential (2e-10)
 245  Args    : none
 246
 247 =cut
 248
 249 sub evalue {
 250         return shift->get_field('evalue');
 251 }
 252
 253 *expect = \&evalue;
 254
 255 =head2 frac_identical
 256
 257  Title   : frac_identical
 258  Usage   : my $frac_id = $hsp->frac_identical( ['query'|'hit'|'total'] );
 259  Function: Returns the fraction of identitical positions for this HSP
 260  Returns : Float in range 0.0 -> 1.0
 261  Args    : 'query' = num identical / length of query seq (without gaps)
 262            'hit'   = num identical / length of hit seq (without gaps)
 263            'total' = num identical / length of alignment (with gaps)
 264            default = 'total'
 265
 266 =cut
 267
 268 sub frac_identical {
 269         my ($self, $type) = @_;
 270
 271         $type = lc $type if defined $type;
 272         $type = 'hit' if (defined $type && $type =~ /subject|sbjct/);
 273         $type = 'total' if (! defined $type || $type eq 'hsp' || $type !~ /query|hit|subject|sbjct|total/);
 274
 275         my $ratio = $self->num_identical($type) / $self->length($type);
 276     return sprintf( "%.4f", $ratio);
 277 }
 278
 279 =head2 frac_conserved
 280
 281  Title    : frac_conserved
 282  Usage    : my $frac_cons = $hsp->frac_conserved( ['query'|'hit'|'total'] );
 283  Function : Returns the fraction of conserved positions for this HSP.
 284             This is the fraction of symbols in the alignment with a
 285             positive score.
 286  Returns : Float in range 0.0 -> 1.0
 287  Args    : 'query' = num conserved / length of query seq (without gaps)
 288            'hit'   = num conserved / length of hit seq (without gaps)
 289            'total' = num conserved / length of alignment (with gaps)
 290            default = 'total'
 291
 292 =cut
 293
 294 sub frac_conserved {
 295     my ($self, $type) = @_;
 296
 297         $type = lc $type if defined $type;
 298         $type = 'hit' if (defined $type && $type =~ /subject|sbjct/);
 299         $type = 'total' if (! defined $type || $type eq 'hsp' || $type !~ /query|hit|subject|sbjct|total/);
 300
 301         my $ratio = $self->num_conserved($type) / $self->length($type);
 302     return sprintf( "%.4f", $ratio);
 303 }
 304
 305 =head2 num_identical
 306
 307  Title   : num_identical
 308  Usage   : $obj->num_identical($newval)
 309  Function: returns the number of identical residues in the alignment
 310  Returns : integer
 311  Args    : integer (optional)
 312
 313 =cut
 314
 315 sub num_identical {
 316     my $self = shift;
 317         return scalar($self->seq_inds('hit', 'identical'));
 318 }
 319
 320 =head2 num_conserved
 321
 322  Title   : num_conserved
 323  Usage   : $obj->num_conserved($newval)
 324  Function: returns the number of conserved residues in the alignment
 325  Returns : inetger
 326  Args    : integer (optional)
 327
 328 =cut
 329
 330 sub num_conserved {
 331     my $self = shift;
 332         return scalar($self->seq_inds('hit', 'conserved-not-identical'));
 333 }
 334
 335 =head2 gaps
 336
 337  Title    : gaps
 338  Usage    : my $gaps = $hsp->gaps( ['query'|'hit'|'total'] );
 339  Function : Get the number of gap characters in the query, hit, or total alignment.
 340  Returns  : Integer, number of gap characters or 0 if none
 341  Args     : 'query', 'hit' or 'total'; default = 'total'
 342
 343 =cut
 344
 345 sub gaps {
 346     my ($self, $type) = @_;
 347     $type = lc $type if defined $type;
 348     $type = 'total' if (! defined $type || $type eq 'hsp' || $type !~ /query|hit|subject|sbjct|total/);
 349     $type = 'hit' if $type =~ /sbjct|subject/;
 350
 351         if ($type eq 'total') {
 352                 return scalar($self->seq_inds('hit', 'gap')) + scalar($self->seq_inds('query', 'gap'));
 353         }
 354         return scalar($self->seq_inds($type, 'gap'));
 355 }
 356
 357 =head2 query_string
 358
 359  Title   : query_string
 360  Usage   : my $qseq = $hsp->query_string;
 361  Function: Retrieves the query sequence of this HSP as a string
 362  Returns : string
 363  Args    : none
 364
 365 =cut
 366
 367 sub query_string {
 368         return shift->get_field('query_string');
 369 }
 370
 371 =head2 hit_string
 372
 373  Title   : hit_string
 374  Usage   : my $hseq = $hsp->hit_string;
 375  Function: Retrieves the hit sequence of this HSP as a string
 376  Returns : string
 377  Args    : none
 378
 379 =cut
 380
 381 sub hit_string {
 382         return shift->get_field('hit_string');
 383 }
 384
 385 =head2 homology_string
 386
 387  Title   : homology_string
 388  Usage   : my $homo_string = $hsp->homology_string;
 389  Function: Retrieves the homology sequence for this HSP as a string.
 390          : The homology sequence is the string of symbols in between the
 391          : query and hit sequences in the alignment indicating the degree
 392          : of conservation (e.g., identical, similar, not similar).
 393  Returns : string
 394  Args    : none
 395
 396 =cut
 397
 398 sub homology_string {
 399         return shift->get_field('homology_string');
 400 }
 401
 402 =head2 length
 403
 404  Title    : length
 405  Usage    : my $len = $hsp->length( ['query'|'hit'|'total'] );
 406  Function : Returns the length of the query or hit in the alignment (without gaps)
 407             or the aggregate length of the HSP (including gaps;
 408             this may be greater than either hit or query )
 409  Returns  : integer
 410  Args     : 'query' = length of query seq (without gaps)
 411             'hit'   = length of hit seq (without gaps)
 412             'total' = length of alignment (with gaps)
 413             default = 'total'
 414  Args    : none
 415
 416 =cut
 417
 418 sub length {
 419     my ($self, $type) = @_;
 420     $type = 'total' unless defined $type;
 421     $type = lc $type;
 422
 423     if ($type =~ /^q/i) {
 424         return $self->query->length;
 425     }
 426         elsif ($type =~ /^(hit|subject|sbjct)/) {
 427         return $self->hit->length;
 428     }
 429         else {
 430         return $self->hit->length + $self->gaps('hit');
 431         }
 432 }
 433
 434 =head2 hsp_length
 435
 436  Title   : hsp_length
 437  Usage   : my $len = $hsp->hsp_length()
 438  Function: shortcut  length('hsp')
 439  Returns : floating point between 0 and 100
 440  Args    : none
 441
 442 =cut
 443
 444 sub hsp_length {
 445         return shift->length('total');
 446 }
 447
 448 =head2 percent_identity
 449
 450  Title   : percent_identity
 451  Usage   : my $percentid = $hsp->percent_identity()
 452  Function: Returns the calculated percent identity for an HSP
 453  Returns : floating point between 0 and 100
 454  Args    : none
 455
 456 =cut
 457
 458 sub percent_identity{
 459         my ($self) = @_;
 460         return $self->frac_identical('hsp') * 100;
 461 }
 462
 463 =head2 get_aln
 464
 465  Title   : get_aln
 466  Usage   : my $aln = $hsp->get_aln
 467  Function: Returns a Bio::SimpleAlign representing the HSP alignment
 468  Returns : Bio::SimpleAlign
 469  Args    : none
 470
 471 =cut
 472
 473 sub get_aln {
 474         my $self = shift;
 475
 476     require Bio::LocatableSeq;
 477     require Bio::SimpleAlign;
 478     my $aln = Bio::SimpleAlign->new();
 479     my $hs = $self->seq('hit');
 480     my $qs = $self->seq('query');
 481         if ($hs && $qs) {
 482                 $aln->add_seq($hs);
 483                 $aln->add_seq($qs);
 484                 return $aln;
 485         }
 486         return;
 487 }
 488
 489 =head2 seq_inds
 490
 491  Title   : seq_inds
 492  Purpose   : Get a list of residue positions (indices) for all identical
 493            : or conserved residues in the query or sbjct sequence.
 494  Example   : @s_ind = $hsp->seq_inds('query', 'identical');
 495            : @h_ind = $hsp->seq_inds('hit', 'conserved');
 496            : @h_ind = $hsp->seq_inds('hit', 'conserved', 1);
 497  Returns   : List of integers
 498            : May include ranges if collapse is true.
 499  Argument  : seq_type  = 'query' or 'hit' or 'sbjct'  (default = query)
 500               ('sbjct' is synonymous with 'hit')
 501              class     = 'identical' or 'conserved' or 'nomatch' or 'gap'
 502                           (default = identical)
 503                           (can be shortened to 'id' or 'cons')
 504                                   Note that 'conserved' includes identical unless you
 505                                   use 'conserved-not-identical'
 506
 507              collapse  = boolean, if true, consecutive positions are merged
 508                          using a range notation, e.g., "1 2 3 4 5 7 9 10 11"
 509                          collapses to "1-5 7 9-11". This is useful for
 510                          consolidating long lists. Default = no collapse.
 511  Throws    : n/a.
 512  Comments  :
 513
 514 See Also   : L<Bio::Search::BlastUtils::collapse_nums()|Bio::Search::BlastUtils>, L<Bio::Search::Hit::HitI::seq_inds()|Bio::Search::Hit::HitI>
 515
 516 =cut
 517
 518 sub seq_inds {
 519     my ($self, $seqType, $class, $collapse) = @_;
 520
 521     $seqType ||= 'query';
 522     $class ||= 'identical';
 523     $collapse ||= 0;
 524     $seqType = lc($seqType);
 525     $class = lc($class);
 526     $seqType = 'hit' if $seqType eq 'sbjct';
 527
 528     my $t = substr($seqType,0,1);
 529     if ($t eq 'q') {
 530         $seqType = 'query';
 531     }
 532     elsif ($t eq 's' || $t eq 'h') {
 533         $seqType = 'hit';
 534     }
 535     else {
 536         $self->warn("unknown seqtype $seqType using 'query'");
 537         $seqType = 'query';
 538     }
 539
 540     $t = substr($class,0,1);
 541     if ($t eq 'c') {
 542         if ($class eq 'conserved-not-identical') {
 543             $class = 'conserved';
 544         }
 545         else {
 546             $class = 'conservedall';
 547         }
 548     }
 549     elsif ($t eq 'i') {
 550         $class = 'identical';
 551     }
 552     elsif ($t eq 'n') {
 553         $class = 'nomatch';
 554     }
 555     elsif ($t eq 'g') {
 556         $class = 'gap';
 557     }
 558     else {
 559         $self->warn("unknown sequence class $class using 'identical'");
 560         $class = 'identical';
 561     }
 562
 563     $seqType .= '_';
 564     $class .= '_inds';
 565
 566     my @ary;
 567     if ($class eq 'conservedall_inds') {
 568                 my %tmp = map { $_, 1 } @{$self->get_field($seqType.'conserved_inds')},
 569                                                                 @{$self->get_field($seqType.'identical_inds')};
 570                 @ary = sort {$a <=> $b} keys %tmp;
 571     }
 572     else {
 573         @ary = @{$self->get_field($seqType.$class)};
 574     }
 575
 576     return $collapse ? &Bio::Search::SearchUtils::collapse_nums(@ary) : @ary;
 577 }
 578
 579 =head2 Inherited from L<Bio::SeqFeature::SimilarityPair>
 580
 581 These methods come from L<Bio::SeqFeature::SimilarityPair>
 582
 583 =head2 query
 584
 585  Title   : query
 586  Usage   : my $query = $hsp->query
 587  Function: Returns a SeqFeature representing the query in the HSP
 588  Returns : Bio::SeqFeature::Similarity
 589  Args    : [optional] new value to set
 590
 591
 592 =head2 hit
 593
 594  Title   : hit
 595  Usage   : my $hit = $hsp->hit
 596  Function: Returns a SeqFeature representing the hit in the HSP
 597  Returns : Bio::SeqFeature::Similarity
 598  Args    : [optional] new value to set
 599
 600
 601 =head2 significance
 602
 603  Title   : significance
 604  Usage   : $evalue = $obj->significance();
 605            $obj->significance($evalue);
 606  Function: Get/Set the significance value (see Bio::SeqFeature::SimilarityPair)
 607  Returns : significance value (scientific notation string)
 608  Args    : significance value (sci notation string)
 609
 610 =cut
 611
 612 sub significance {
 613         return shift->get_field('evalue');
 614 }
 615
 616 =head2 score
 617
 618  Title   : score
 619  Usage   : my $score = $hsp->score();
 620  Function: Returns the score for this HSP or undef
 621  Returns : numeric
 622  Args    : [optional] numeric to set value
 623
 624 =cut
 625
 626 sub score {
 627         return shift->get_field('score');
 628 }
 629
 630 =head2 bits
 631
 632  Title   : bits
 633  Usage   : my $bits = $hsp->bits();
 634  Function: Returns the bit value for this HSP or undef
 635  Returns : numeric
 636  Args    : none
 637
 638 =cut
 639
 640 sub bits {
 641         return shift->get_field('bits');
 642 }
 643
 644 # override
 645
 646 =head2 strand
 647
 648  Title   : strand
 649  Usage   : $hsp->strand('query')
 650  Function: Retrieves the strand for the HSP component requested
 651  Returns : +1 or -1 (0 if unknown)
 652  Args    : 'hit' or 'subject' or 'sbjct' to retrieve the strand of the subject
 653            'query' to retrieve the query strand (default)
 654            'list' or 'array' to retrieve both query and hit together
 655
 656 =cut
 657
 658 sub strand {
 659     my $self = shift;
 660     my $val = shift;
 661     $val = 'query' unless defined $val;
 662     $val =~ s/^\s+//;
 663
 664     if ($val =~ /^q/i) {
 665         return $self->query->strand(@_);
 666     }
 667     elsif ($val =~ /^hi|^s/i) {
 668         return $self->hit->strand(@_);
 669     }
 670     elsif ($val =~ /^list|array/i) {
 671         return ($self->query->strand(@_), $self->hit->strand(@_) );
 672     }
 673     else {
 674         $self->warn("unrecognized component '$val' requested\n");
 675     }
 676     return 0;
 677 }
 678
 679 =head2 start
 680
 681  Title   : start
 682  Usage   : $hsp->start('query')
 683  Function: Retrieves the start for the HSP component requested
 684  Returns : integer, or list of two integers (query start and subject start) in
 685            list context
 686  Args    : 'hit' or 'subject' or 'sbjct' to retrieve the start of the subject
 687            'query' to retrieve the query start (default)
 688
 689 =cut
 690
 691 sub start {
 692     my $self = shift;
 693     my $val = shift;
 694     $val = (wantarray ? 'list' : 'query') unless defined $val;
 695     $val =~ s/^\s+//;
 696
 697     if ($val =~ /^q/i) {
 698         return $self->query->start(@_);
 699     }
 700     elsif ($val =~ /^(hi|s)/i) {
 701         return $self->hit->start(@_);
 702     }
 703     elsif ($val =~ /^list|array/i) {
 704         return ($self->query->start(@_), $self->hit->start(@_) );
 705     }
 706     else {
 707         $self->warn("unrecognized component '$val' requested\n");
 708     }
 709     return 0;
 710 }
 711
 712 =head2 end
 713
 714  Title   : end
 715  Usage   : $hsp->end('query')
 716  Function: Retrieves the end for the HSP component requested
 717  Returns : integer, or list of two integers (query end and subject end) in
 718            list context
 719  Args    : 'hit' or 'subject' or 'sbjct' to retrieve the end of the subject
 720            'query' to retrieve the query end (default)
 721
 722 =cut
 723
 724 sub end {
 725     my $self = shift;
 726     my $val = shift;
 727     $val = (wantarray ? 'list' : 'query') unless defined $val;
 728     $val =~ s/^\s+//;
 729
 730     if ($val =~ /^q/i) {
 731         return $self->query->end(@_);
 732     }
 733     elsif ($val =~ /^(hi|s)/i) {
 734         return $self->hit->end(@_);
 735     }
 736     elsif ($val =~ /^list|array/i) {
 737         return ($self->query->end(@_), $self->hit->end(@_) );
 738     }
 739     else {
 740         $self->warn("unrecognized end component '$val' requested\n");
 741     }
 742     return 0;
 743 }
 744
 745 =head2 seq
 746
 747  Usage     : $hsp->seq( [seq_type] );
 748  Purpose   : Get the query or sbjct sequence as a Bio::Seq.pm object.
 749  Example   : $seqObj = $hsp->seq('query');
 750  Returns   : Object reference for a Bio::LocatableSeq object.
 751  Argument  : seq_type = 'query' or 'hit' or 'sbjct' (default = 'query').
 752            : ('sbjct' is synonymous with 'hit')
 753            : default is 'query'
 754
 755 =cut
 756
 757 sub seq {
 758     my ($self, $seqType) = @_;
 759     $seqType ||= 'query';
 760     $seqType = 'hit' if $seqType eq 'sbjct';
 761     if ($seqType =~ /^(m|ho)/i ) {
 762         $self->throw("cannot call seq on the homology match string, it isn't really a sequence, use get_aln to convert the HSP to a Bio::AlignIO and generate a consensus from that.");
 763     }
 764
 765     my $str = $self->seq_str($seqType) || return;
 766     require Bio::LocatableSeq;
 767     my $id = ($seqType =~ /^q/i) ? $self->query->seq_id : $self->hit->seq_id;
 768     return Bio::LocatableSeq->new(  -ID        => $id,
 769                                     -SEQ       => $str,
 770                                     -START     => $self->start($seqType),
 771                                     -END       => $self->end($seqType),
 772                                     -STRAND    => $self->strand($seqType),
 773                                     -FORCE_NSE => $id ? 0 : 1,
 774                                     -DESC      => "$seqType sequence " );
 775
 776 }
 777
 778 =head2 seq_str
 779
 780  Usage     : $hsp->seq_str( seq_type );
 781  Purpose   : Get the full query, sbjct, or 'match' sequence as a string.
 782            : The 'match' sequence is the string of symbols in between the
 783            : query and sbjct sequences.
 784  Example   : $str = $hsp->seq_str('query');
 785  Returns   : String
 786  Argument  : seq_Type = 'query' or 'hit' or 'sbjct' or 'match'
 787            :  ('sbjct' is synonymous with 'hit')
 788            : default is 'query'
 789  Throws    : Exception if the argument does not match an accepted seq_type.
 790  Comments  :
 791
 792 See Also   : L<seq()|seq>, L<seq_inds()|seq_inds>, B<_set_match_seq()>
 793
 794 =cut
 795
 796 sub seq_str {
 797     my $self = shift;
 798     my $type = shift || 'query';
 799
 800     if ($type =~ /^q/i) {
 801         return $self->query_string(@_);
 802     }
 803     elsif ($type =~ /^(s)|(hi)/i) {
 804         return $self->hit_string(@_);
 805     }
 806     elsif ($type =~ /^(ho)|(ma)/i) {
 807         return $self->homology_string(@_);
 808     }
 809     else {
 810         $self->warn("unknown sequence type $type");
 811     }
 812     return '';
 813 }
 814
 815 =head2 rank
 816
 817  Usage     : $hsp->rank( [string] );
 818  Purpose   : Get the rank of the HSP within a given Blast hit.
 819  Example   : $rank = $hsp->rank;
 820  Returns   : Integer (1..n) corresponding to the order in which the HSP
 821              appears in the BLAST report.
 822
 823 =cut
 824
 825 sub rank {
 826     return shift->get_field('rank');
 827 }
 828
 829 =head2 matches
 830
 831  Usage     : $hsp->matches(-seq   => 'hit'|'query',
 832                            -start => $start,
 833                            -stop  => $stop);
 834  Purpose   : Get the total number of identical and conservative matches
 835            : in the query or sbjct sequence for the given HSP. Optionally can
 836            : report data within a defined interval along the seq.
 837  Example   : ($id,$cons) = $hsp_object->matches(-seq   => 'hit');
 838            : ($id,$cons) = $hsp_object->matches(-seq   => 'query',
 839                                                 -start => 300,
 840                                                 -stop  => 400);
 841  Returns   : 2-element array of integers
 842  Argument  : (1) seq_type = 'query' or 'hit' or 'sbjct' (default = query)
 843            :  ('sbjct' is synonymous with 'hit')
 844            : (2) start = Starting coordinate (optional)
 845            : (3) stop  = Ending coordinate (optional)
 846
 847 =cut
 848
 849 sub matches {
 850     my ($self, @args) = @_;
 851     my($seqType, $beg, $end) = $self->_rearrange([qw(SEQ START STOP)], @args);
 852     $seqType ||= 'query';
 853     $seqType = 'hit' if $seqType eq 'sbjct';
 854
 855     my @data;
 856     if ((!defined $beg && !defined $end) || ! $self->seq_str('match')) {
 857         push @data, ($self->num_identical, $self->num_conserved);
 858     }
 859         else {
 860         $beg ||= 0;
 861         $end ||= 0;
 862         my ($start, $stop) = $self->range($seqType);
 863
 864         if ($beg == 0) {
 865                         $beg = $start;
 866                         $end = $beg+$end;
 867                 }
 868         elsif ($end == 0) {
 869                         $end = $stop;
 870                         $beg = $end-$beg;
 871                 }
 872
 873         if ($end >= $stop) {
 874                         $end = $stop;
 875                 }
 876         else {
 877                         $end += 1;
 878                 }
 879         if ($beg < $start) {
 880                         $beg = $start;
 881                 }
 882
 883         my $seq = substr($self->seq_str('homology'), $beg-$start, ($end-$beg));
 884
 885         if (!CORE::length $seq) {
 886             $self->throw("Undefined sub-sequence ($beg,$end). Valid range = $start - $stop");
 887         }
 888         ## Get data for a substring.
 889         $seq =~ s/ //g;  # remove space (no info).
 890         my $len_cons = CORE::length $seq;
 891         $seq =~ s/\+//g;  # remove '+' characters (conservative substitutions)
 892         my $len_id = CORE::length $seq;
 893         push @data, ($len_id, $len_cons);
 894     }
 895
 896     return @data;
 897 }
 898
 899 =head2 n
 900
 901  Usage     : $hsp_obj->n()
 902  Purpose   : Get the N value (num HSPs on which P/Expect is based).
 903  Returns   : Integer or null string if not defined.
 904  Argument  : n/a
 905  Throws    : n/a
 906  Comments  : The 'N' value is listed in parenthesis with P/Expect value:
 907            : e.g., P(3) = 1.2e-30  ---> (N = 3).
 908            : Not defined in NCBI Blast2 with gaps.
 909            : This typically is equal to the number of HSPs but not always.
 910
 911 =cut
 912
 913 sub n {
 914     return shift->get_field('num_hsps');
 915 }
 916
 917 =head2 range
 918
 919  Usage     : $hsp->range( [seq_type] );
 920  Purpose   : Gets the (start, end) coordinates for the query or sbjct sequence
 921            : in the HSP alignment.
 922  Example   : ($query_beg, $query_end) = $hsp->range('query');
 923            : ($hit_beg, $hit_end) = $hsp->range('hit');
 924  Returns   : Two-element array of integers
 925  Argument  : seq_type = string, 'query' or 'hit' or 'sbjct'  (default = 'query')
 926            :  ('sbjct' is synonymous with 'hit')
 927  Throws    : n/a
 928  Comments  : This is a convenience method for constructions such as
 929              ($hsp->query->start, $hsp->query->end)
 930
 931 =cut
 932
 933 sub range {
 934     my ($self, $seqType) = @_;
 935     $seqType ||= 'query';
 936
 937     my ($start, $end);
 938     if ($seqType eq 'query') {
 939         $start = $self->query->start;
 940         $end = $self->query->end;
 941     }
 942     else {
 943         $start = $self->hit->start;
 944         $end = $self->hit->end;
 945     }
 946     return ($start, $end);
 947 }
 948
 949 #*** would want cigar stuff from GenericHSP - move to HSPI?
 950
 951 1;
 952