← Back to branch summary

~ubuntu-branches/ubuntu/lucid/bioperl/lucid

« back to all changes in this revision

Viewing changes to Bio/Search/Hit/GenericHit.pm

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Search/Hit/GenericHit.pm
1
 
# $Id: GenericHit.pm,v 1.5 2002/02/07 22:12:41 jason Exp $
 
1
# $Id: GenericHit.pm,v 1.29 2003/11/25 18:46:14 jason Exp $
2
2
#
3
3
# BioPerl module for Bio::Search::Hit::GenericHit
4
4
#
16
16
 
17
17
=head1 SYNOPSIS
18
18
 
19
 
  {
20
19
    use Bio::Search::Hit::GenericHit;
21
20
    my $hit = new Bio::Search::Hit::GenericHit(-algorithm => 'blastp');
22
21
 
23
 
  }
 
22
    # typically one gets HitI objects from a SearchIO stream via a ResultI
 
23
    use Bio::SearchIO;
 
24
    my $parser = new Bio::SearchIO(-format => 'blast', -file => 'result.bls');
 
25
 
 
26
    my $result = $parser->next_result;
 
27
    my $hit    = $result->next_hit;
 
28
 
 
29
# TODO: Describe how to configure a SearchIO stream so that it generates
 
30
#       GenericHit objects.
24
31
 
25
32
=head1 DESCRIPTION
26
33
 
27
34
This object handles the hit data from a Database Sequence Search such
28
35
as FASTA or BLAST.
29
36
 
 
37
Unless you're writing a parser, you won't ever need to create a
 
38
GenericHit or any other HitI-implementing object. If you use
 
39
the SearchIO system, HitI objects are created automatically from
 
40
a SearchIO stream which returns Bio::Search::Hit::HitI objects.
 
41
 
 
42
For documentation on what you can do with GenericHit (and other HitI
 
43
objects), please see the API documentation in
 
44
L<Bio::Search::Hit::HitI|Bio::Search::Hit::HitI>.
 
45
 
30
46
=head1 FEEDBACK
31
47
 
32
48
=head2 Mailing Lists
45
61
email or the web:
46
62
 
47
63
  bioperl-bugs@bioperl.org
48
 
  http://bioperl.org/bioperl-bugs/
 
64
  http://bugzilla.bioperl.org/
49
65
 
50
66
=head1 AUTHOR - Jason Stajich and Steve Chervitz
51
67
 
73
89
 
74
90
use Bio::Root::Root;
75
91
use Bio::Search::Hit::HitI;
 
92
require Bio::Search::SearchUtils;
76
93
 
77
94
@ISA = qw(Bio::Root::Root Bio::Search::Hit::HitI );
78
95
 
87
104
           -accession    => Accession number (optional)
88
105
           -length       => Length of the Hit (optional)
89
106
           -score        => Raw Score for the Hit (optional)
 
107
           -bits         => Bit Score for the Hit (optional)
90
108
           -significance => Significance value for the Hit (optional)
91
109
           -algorithm    => Algorithm used (BLASTP, FASTX, etc...)
92
110
           -hsps         => Array ref of HSPs for this Hit. 
 
111
           -found_again  => boolean, true if hit appears in a 
 
112
                            "previously found" section of a PSI-Blast report.
93
113
 
94
114
=cut
95
115
 
97
117
  my($class,@args) = @_;
98
118
 
99
119
  my $self = $class->SUPER::new(@args);
100
 
  my ($hsps, $name,$desc, $acc, $length,
101
 
      $score,$algo,$signif) = $self->_rearrange([qw(HSPS NAME DESCRIPTION
102
 
                                                   ACCESSION
103
 
                                                   LENGTH SCORE ALGORITHM 
104
 
                                                   SIGNIFICANCE)], @args);
 
120
  my ($hsps, $name,$query_len,$desc, $acc, $locus, $length,
 
121
      $score,$algo,$signif,$bits,
 
122
      $rank) = $self->_rearrange([qw(HSPS 
 
123
                                     NAME 
 
124
                                     QUERY_LEN
 
125
                                     DESCRIPTION
 
126
                                     ACCESSION
 
127
                                     LOCUS
 
128
                                     LENGTH SCORE ALGORITHM 
 
129
                                     SIGNIFICANCE BITS
 
130
                                     RANK )], @args);
 
131
  
 
132
  defined $query_len && $self->query_length($query_len);
105
133
 
106
134
  if( ! defined $name ) { 
107
135
      $self->throw("Must have defined a valid name for Hit");
108
136
  } else { 
109
137
      $self->name($name);
110
138
  }  
 
139
 
111
140
  defined $acc    && $self->accession($acc);
 
141
  defined $locus  && $self->locus($locus);
112
142
  defined $desc   && $self->description($desc);
113
143
  defined $length && $self->length($length);
114
144
  defined $algo   && $self->algorithm($algo);
115
145
  defined $signif && $self->significance($signif);
116
146
  defined $score  && $self->raw_score($score);
117
 
  
 
147
  defined $bits   && $self->bits($bits);
 
148
  defined $rank   && $self->rank($rank);
 
149
 
118
150
  $self->{'_iterator'} = 0;
119
151
  $self->{'_hsps'} = [];
120
152
  if( defined $hsps  ) {
121
153
      if( ref($hsps) !~ /array/i ) {
122
 
          $self->warn("Did not specify a valid array ref for the param HSPS ($hsps)");
 
154
          $self->warn("Did not specify a valid array ref for the param HSPS ($hsps)");
123
155
      } else {
124
 
          while( @$hsps ) { 
125
 
              $self->add_hsp(shift @$hsps );
126
 
          }
 
156
          while( @$hsps ) { 
 
157
              $self->add_hsp(shift @$hsps );
 
158
          }
127
159
      }
128
160
  } 
129
161
  return $self;
143
175
sub add_hsp {
144
176
   my ($self,$hsp) = @_;
145
177
   if( !defined $hsp || ! $hsp->isa('Bio::Search::HSP::HSPI') ) { 
146
 
       $self->warn("Must provide a valid Bio::Search::HSP::HSPI object to object: $self method: add_hsp");
 
178
       $self->throw("Must provide a valid Bio::Search::HSP::HSPI object to object: $self method: add_hsp value: $hsp");
147
179
       return undef;
148
180
   }
 
181
#   print STDERR "GenericHit::add_hsp()\n";
149
182
   push @{$self->{'_hsps'}}, $hsp;
150
183
   return scalar @{$self->{'_hsps'}};
151
184
}
170
203
    my ($self,$value) = @_;
171
204
    my $previous = $self->{'_name'};
172
205
    if( defined $value || ! defined $previous ) {
173
 
        $value = $previous = '' unless defined $value;
174
 
        $self->{'_name'} = $value;
 
206
        $value = $previous = '' unless defined $value;
 
207
        $self->{'_name'} = $value;
175
208
    } 
176
209
    return $previous;
177
210
}
190
223
    my ($self,$value) = @_;
191
224
    my $previous = $self->{'_accession'};
192
225
    if( defined $value || ! defined $previous ) { 
193
 
        $value = $previous = '' unless defined $value;
194
 
        $self->{'_accession'} = $value;
 
226
        $value = $previous = '' unless defined $value;
 
227
        $self->{'_accession'} = $value;
195
228
    } 
196
 
    return $previous;
 
229
        return $previous;
197
230
}
198
231
 
199
232
=head2 description
210
243
    my ($self,$value) = @_;
211
244
    my $previous = $self->{'_description'};
212
245
    if( defined $value || ! defined $previous ) { 
213
 
        $value = $previous = '' unless defined $value;
214
 
        $self->{'_description'} = $value;
 
246
        $value = $previous = '' unless defined $value;
 
247
        $self->{'_description'} = $value;
215
248
    } 
216
249
    return $previous;
217
250
}
230
263
    my ($self,$value) = @_;
231
264
    my $previous = $self->{'_length'};
232
265
    if( defined $value || ! defined $previous ) { 
233
 
        $value = $previous = 0 unless defined $value;
234
 
        $self->{'_length'} = $value;
 
266
        $value = $previous = 0 unless defined $value;
 
267
        $self->{'_length'} = $value;
235
268
    } 
236
269
    return $previous;
237
270
}
255
288
    my ($self,$value) = @_;
256
289
    my $previous = $self->{'_algorithm'};
257
290
    if( defined $value || ! defined $previous ) { 
258
 
        $value = $previous = '' unless defined $value;
259
 
        $self->{'_algorithm'} = $value;
 
291
        $value = $previous = '' unless defined $value;
 
292
        $self->{'_algorithm'} = $value;
260
293
    } 
261
294
    return $previous;
262
295
}
277
310
    my ($self,$value) = @_;
278
311
    my $previous = $self->{'_score'};
279
312
    if( defined $value || ! defined $previous ) { 
280
 
        $value = $previous = '' unless defined $value;
281
 
        $self->{'_score'} = $value;
 
313
        $value = $previous = '' unless defined $value;
 
314
        $self->{'_score'} = $value;
282
315
    } 
283
316
    return $previous;
284
317
}
285
318
 
 
319
=head2 score
 
320
 
 
321
Equivalent to L<raw_score()|raw_score>
 
322
 
 
323
=cut
 
324
 
 
325
sub score { shift->raw_score(@_); }
 
326
 
286
327
=head2 significance
287
328
 
288
329
 Title   : significance
299
340
sub significance {
300
341
    my ($self,$value) = @_;
301
342
    my $previous = $self->{'_significance'};
302
 
    if( defined $value || ! defined $previous ) { 
303
 
        $value = $previous = '' unless defined $value;
304
 
        $self->{'_significance'} = $value;
305
 
    } 
 
343
    if( defined $value ) { 
 
344
        $self->{'_significance'} = $value;
 
345
    } elsif ( ! defined $previous ) {
 
346
        unless( defined $self->{'_hsps'}->[0] ) {
 
347
            $self->warn("No HSPs for this Hit (".$self->name.")");
 
348
            return undef;
 
349
        }
 
350
        # Set the significance of the Hit to that of the top HSP.
 
351
        $previous = $self->{'_significance'} = ($self->hsps)[0]->significance;
 
352
    }
 
353
 
 
354
    return $previous;
 
355
}
 
356
 
 
357
=head2 bits
 
358
 
 
359
 Usage     : $hit_object->bits();
 
360
 Purpose   : Gets the bit score of the best HSP for the current hit.
 
361
 Example   : $bits = $hit_object->bits();
 
362
 Returns   : Integer or undef if bit score is not set
 
363
 Argument  : n/a
 
364
 Comments  : For BLAST1, the non-bit score is listed in the summary line.
 
365
 
 
366
See Also   : L<score()|score>
 
367
 
 
368
=cut
 
369
 
 
370
#---------
 
371
sub bits { 
 
372
#---------
 
373
    my ($self,$value) = @_; 
 
374
    my $previous = $self->{'_bits'};
 
375
    if( defined $value ) { 
 
376
        $self->{'_bits'} = $value;
 
377
    } elsif ( ! defined $previous ) {
 
378
        # Set the significance of the Hit to that of the top HSP.
 
379
        unless( defined $self->{'_hsps'}->[0] ) {
 
380
            $self->warn("No HSPs for this Hit (".$self->name.")");
 
381
            return undef;
 
382
        }
 
383
        $previous = $self->{'_bits'} = ($self->hsps)[0]->bits;
 
384
    }    
306
385
    return $previous;
307
386
}
308
387
 
324
403
    return $self->{'_hsps'}->[$self->{'_iterator'}++];    
325
404
}
326
405
 
 
406
 
 
407
=head2 hsps
 
408
 
 
409
 Usage     : $hit_object->hsps();
 
410
 Purpose   : Get a list containing all HSP objects.
 
411
           : Get the numbers of HSPs for the current hit.
 
412
 Example   : @hsps = $hit_object->hsps();
 
413
           : $num  = $hit_object->hsps();  # alternatively, use num_hsps()
 
414
 Returns   : Array context : list of Bio::Search::HSP::BlastHSP.pm objects.
 
415
           : Scalar context: integer (number of HSPs).
 
416
           :                 (Equivalent to num_hsps()).
 
417
 Argument  : n/a. Relies on wantarray
 
418
 Throws    : Exception if the HSPs have not been collected.
 
419
 
 
420
See Also   : L<hsp()|hsp>, L<num_hsps()|num_hsps>
 
421
 
 
422
=cut
 
423
 
 
424
#---------
 
425
sub hsps {
 
426
#---------
 
427
   my $self = shift;
 
428
   
 
429
   if (not ref $self->{'_hsps'}) {
 
430
       $self->throw("Can't get HSPs: data not collected.");
 
431
   }
 
432
   
 
433
   return wantarray 
 
434
       #  returning list containing all HSPs.
 
435
       ? @{$self->{'_hsps'}}
 
436
   #  returning number of HSPs.
 
437
   : scalar(@{$self->{'_hsps'}});
 
438
}
 
439
 
 
440
=head2 num_hsps
 
441
 
 
442
 Usage     : $hit_object->num_hsps();
 
443
 Purpose   : Get the number of HSPs for the present Blast hit.
 
444
 Example   : $nhsps = $hit_object->num_hsps();
 
445
 Returns   : Integer
 
446
 Argument  : n/a
 
447
 Throws    : Exception if the HSPs have not been collected.
 
448
 
 
449
See Also   : L<hsps()|hsps>
 
450
 
 
451
=cut
 
452
 
 
453
#-------------
 
454
sub num_hsps {
 
455
    my $self = shift;
 
456
    
 
457
    if (not defined $self->{'_hsps'}) {
 
458
        $self->throw("Can't get HSPs: data not collected.");
 
459
    }
 
460
 
 
461
    return scalar(@{$self->{'_hsps'}});
 
462
}
 
463
 
327
464
=head2 rewind
328
465
 
329
466
 Title   : rewind
330
467
 Usage   : $hit->rewind;
331
 
 Function: Allow one to reset the HSP iteration to the beginning
 
468
 Function: Allow one to reset the HSP iterator to the beginning
332
469
           Since this is an in-memory implementation
333
470
 Returns : none
334
471
 Args    : none
340
477
   $self->{'_iterator'} = 0;
341
478
}
342
479
 
 
480
=head2 ambiguous_aln
 
481
 
 
482
 Usage     : $ambig_code = $hit_object->ambiguous_aln();
 
483
 Purpose   : Sets/Gets ambiguity code data member.
 
484
 Example   : (see usage)
 
485
 Returns   : String = 'q', 's', 'qs', '-'
 
486
           :   'q'  = query sequence contains overlapping sub-sequences 
 
487
           :          while sbjct does not.
 
488
           :   's'  = sbjct sequence contains overlapping sub-sequences 
 
489
           :          while query does not.
 
490
           :   'qs' = query and sbjct sequence contains overlapping sub-sequences
 
491
           :          relative to each other.
 
492
           :   '-'  = query and sbjct sequence do not contains multiple domains 
 
493
           :          relative to each other OR both contain the same distribution
 
494
           :          of similar domains.
 
495
 Argument  : n/a
 
496
 Throws    : n/a
 
497
 Comment   : Note: "sbjct" is synonymous with "hit"
 
498
 
 
499
=cut
 
500
 
 
501
#--------------------
 
502
sub ambiguous_aln { 
 
503
#--------------------
 
504
    my $self = shift;
 
505
    if(@_) { $self->{'_ambiguous_aln'} = shift; }
 
506
    $self->{'_ambiguous_aln'} || '-';
 
507
}
 
508
 
 
509
=head2 overlap
 
510
 
 
511
See documentation in L<Bio::Search::Hit::HitI::overlap()|Bio::Search::Hit::HitI>
 
512
 
 
513
=cut
 
514
 
 
515
#-------------
 
516
sub overlap { 
 
517
#-------------
 
518
    my $self = shift; 
 
519
    if(@_) { $self->{'_overlap'} = shift; }
 
520
    defined $self->{'_overlap'} ? $self->{'_overlap'} : 0;
 
521
}
 
522
 
 
523
 
 
524
=head2 n
 
525
 
 
526
 Usage     : $hit_object->n();
 
527
 Purpose   : Gets the N number for the current hit.
 
528
           : This is the number of HSPs in the set which was ascribed
 
529
           : the lowest P-value (listed on the description line).
 
530
           : This number is not the same as the total number of HSPs.
 
531
           : To get the total number of HSPs, use num_hsps().
 
532
 Example   : $n = $hit_object->n();
 
533
 Returns   : Integer
 
534
 Argument  : n/a
 
535
 Throws    : Exception if HSPs have not been set (BLAST2 reports).
 
536
 Comments  : Note that the N parameter is not reported in gapped BLAST2.
 
537
           : Calling n() on such reports will result in a call to num_hsps().
 
538
           : The num_hsps() method will count the actual number of
 
539
           : HSPs in the alignment listing, which may exceed N in
 
540
           : some cases.
 
541
 
 
542
See Also   : L<num_hsps()|num_hsps>
 
543
 
 
544
=cut
 
545
 
 
546
#-----
 
547
sub n { 
 
548
#-----
 
549
    my $self = shift; 
 
550
 
 
551
    # The check for $self->{'_n'} is a remnant from the 'query' mode days
 
552
    # in which the sbjct object would collect data from the description 
 
553
    # line only.
 
554
 
 
555
    my ($n);
 
556
    if(not defined($self->{'_n'})) {
 
557
        if( $self->hsp ) {
 
558
            $n = $self->hsp->n;
 
559
        }
 
560
    } else {
 
561
        $n = $self->{'_n'}; 
 
562
    } 
 
563
    $n ||= $self->num_hsps;
 
564
 
 
565
    return $n;
 
566
}
 
567
 
 
568
=head2 p
 
569
 
 
570
 Usage     : $hit_object->p( [format] );
 
571
 Purpose   : Get the P-value for the best HSP of the given BLAST hit.
 
572
           : (Note that P-values are not provided with NCBI Blast2 reports).
 
573
 Example   : $p =  $sbjct->p;
 
574
           : $p =  $sbjct->p('exp');  # get exponent only.
 
575
           : ($num, $exp) =  $sbjct->p('parts');  # split sci notation into parts
 
576
 Returns   : Float or scientific notation number (the raw P-value, DEFAULT).
 
577
           : Integer if format == 'exp' (the magnitude of the base 10 exponent).
 
578
           : 2-element list (float, int) if format == 'parts' and P-value
 
579
           :                is in scientific notation (See Comments).
 
580
 Argument  : format: string of 'raw' | 'exp' | 'parts'
 
581
           :    'raw' returns value given in report. Default. (1.2e-34)
 
582
           :    'exp' returns exponent value only (34)
 
583
           :    'parts' returns the decimal and exponent as a 
 
584
           :            2-element list (1.2, -34) (See Comments).
 
585
 Throws    : Warns if no P-value is defined. Uses expect instead.
 
586
 Comments  : Using the 'parts' argument is not recommended since it will not
 
587
           : work as expected if the P-value is not in scientific notation.
 
588
           : That is, floats are not converted into sci notation before
 
589
           : splitting into parts.
 
590
 
 
591
See Also   : L<expect()|expect>, L<signif()|signif>, L<Bio::Search::SearchUtils::get_exponent()|Bio::Search::SearchUtils>
 
592
 
 
593
=cut
 
594
 
 
595
#--------
 
596
sub p { 
 
597
#--------
 
598
# Some duplication of logic for p(), expect() and signif() for the sake of performance.
 
599
    my ($self, $fmt) = @_;
 
600
 
 
601
    my $val = $self->{'_p'};
 
602
 
 
603
    # $val can be zero.
 
604
    if(not defined $val) {
 
605
        # P-value not defined, must be a NCBI Blast2 report.
 
606
        # Use expect instead.
 
607
        $self->warn( "P-value not defined. Using expect() instead.");
 
608
        $val = $self->{'_expect'};
 
609
    }
 
610
 
 
611
    return $val if not $fmt or $fmt =~ /^raw/i;
 
612
    ## Special formats: exponent-only or as list.
 
613
    return &Bio::Search::SearchUtils::get_exponent($val) if $fmt =~ /^exp/i;
 
614
    return (split (/eE/, $val)) if $fmt =~ /^parts/i;
 
615
 
 
616
    ## Default: return the raw P-value.
 
617
    return $val;
 
618
}
 
619
 
 
620
=head2 hsp
 
621
 
 
622
 Usage     : $hit_object->hsp( [string] );
 
623
 Purpose   : Get a single HSPI object for the present HitI object.
 
624
 Example   : $hspObj  = $hit_object->hsp;  # same as 'best'
 
625
           : $hspObj  = $hit_object->hsp('best');
 
626
           : $hspObj  = $hit_object->hsp('worst');
 
627
 Returns   : Object reference for a Bio::Search::HSP::BlastHSP.pm object.
 
628
 Argument  : String (or no argument).
 
629
           :   No argument (default) = highest scoring HSP (same as 'best').
 
630
           :   'best' or 'first' = highest scoring HSP.
 
631
           :   'worst' or 'last' = lowest scoring HSP.
 
632
 Throws    : Exception if the HSPs have not been collected.
 
633
           : Exception if an unrecognized argument is used.
 
634
 
 
635
See Also   : L<hsps()|hsps>, L<num_hsps>()
 
636
 
 
637
=cut
 
638
 
 
639
#----------
 
640
sub hsp {
 
641
#----------
 
642
    my( $self, $option ) = @_;
 
643
    $option ||= 'best';
 
644
    
 
645
    if (not ref $self->{'_hsps'}) {
 
646
        $self->throw("Can't get HSPs: data not collected.");
 
647
    }
 
648
 
 
649
    my @hsps = @{$self->{'_hsps'}};
 
650
    
 
651
    return $hsps[0]      if $option =~ /best|first|1/i;
 
652
    return $hsps[$#hsps] if $option =~ /worst|last/i;
 
653
 
 
654
    $self->throw("Can't get HSP for: $option\n" .
 
655
                 "Valid arguments: 'best', 'worst'");
 
656
}
 
657
 
 
658
=head2 logical_length
 
659
 
 
660
 Usage     : $hit_object->logical_length( [seq_type] );
 
661
           : (mostly intended for internal use).
 
662
 Purpose   : Get the logical length of the hit sequence.
 
663
           : If the Blast is a TBLASTN or TBLASTX, the returned length 
 
664
           : is the length of the would-be amino acid sequence (length/3).
 
665
           : For all other BLAST flavors, this function is the same as length().
 
666
 Example   : $len    = $hit_object->logical_length();
 
667
 Returns   : Integer 
 
668
 Argument  : seq_type = 'query' or 'hit' or 'sbjct' (default = 'query')
 
669
             ('sbjct' is synonymous with 'hit')
 
670
 Throws    : n/a
 
671
 Comments  : This is important for functions like frac_aligned_query()
 
672
           : which need to operate in amino acid coordinate space when dealing
 
673
           : with [T]BLAST[NX] type reports.
 
674
 
 
675
See Also   : L<length()|length>, L<frac_aligned_query()|frac_aligned_query>, L<frac_aligned_hit()|frac_aligned_hit>
 
676
 
 
677
=cut
 
678
 
 
679
#--------------------
 
680
sub logical_length {
 
681
#--------------------
 
682
    my $self = shift;
 
683
    my $seqType = shift || 'query';
 
684
    $seqType = 'sbjct' if $seqType eq 'hit';
 
685
 
 
686
    my $length;
 
687
 
 
688
    # For the sbjct, return logical sbjct length
 
689
    if( $seqType eq 'sbjct' ) {
 
690
        $length = $self->length;
 
691
        # Adjust length based on BLAST flavor.
 
692
        if($self->algorithm =~ /TBLAST[NX]/ ) {
 
693
            $length /= 3;
 
694
        }
 
695
    } else {
 
696
        # Otherwise, return logical query length
 
697
        $length = $self->query_length();
 
698
        $self->throw("Must have defined query_len") unless ( $length );
 
699
 
 
700
        # Adjust length based on BLAST flavor.
 
701
        if($self->algorithm =~ /T?BLASTX/ ) {
 
702
            $length /= 3;
 
703
        }
 
704
    }
 
705
    return int($length);
 
706
}
 
707
 
 
708
=head2 length_aln
 
709
 
 
710
 Usage     : $hit_object->length_aln( [seq_type] );
 
711
 Purpose   : Get the total length of the aligned region for query or sbjct seq.
 
712
           : This number will include all HSPs
 
713
 Example   : $len    = $hit_object->length_aln(); # default = query
 
714
           : $lenAln = $hit_object->length_aln('query');
 
715
 Returns   : Integer 
 
716
 Argument  : seq_Type = 'query' or 'hit' or 'sbjct' (Default = 'query')
 
717
             ('sbjct' is synonymous with 'hit')
 
718
 Throws    : Exception if the argument is not recognized.
 
719
 Comments  : This method will report the logical length of the alignment,
 
720
           : meaning that for TBLAST[NX] reports, the length is reported
 
721
           : using amino acid coordinate space (i.e., nucleotides / 3).
 
722
           : 
 
723
           : This method requires that all HSPs be tiled. If they have not
 
724
           : already been tiled, they will be tiled first automatically..
 
725
           : If you don't want the tiled data, iterate through each HSP
 
726
           : calling length() on each (use hsps() to get all HSPs).
 
727
 
 
728
See Also   : L<length()|length>, L<frac_aligned_query()|frac_aligned_query>, L<frac_aligned_hit()|frac_aligned_hit>, L<gaps()|gaps>, L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>, L<Bio::Search::HSP::BlastHSP::length()|Bio::Search::HSP::BlastHSP>
 
729
 
 
730
=cut
 
731
 
 
732
#---------------'
 
733
sub length_aln {
 
734
#---------------
 
735
    my( $self, $seqType, $num ) = @_;
 
736
    
 
737
    $seqType ||= 'query';
 
738
    $seqType = 'sbjct' if $seqType eq 'hit';
 
739
 
 
740
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
741
 
 
742
    if( defined $num) {
 
743
        return $self->{'_length_aln_'.$seqType} = $num;
 
744
    }
 
745
 
 
746
    my $data = $self->{'_length_aln_'.$seqType};
 
747
    
 
748
    ## If we don't have data, figure out what went wrong.
 
749
    if(!$data) {
 
750
        $self->throw("Can't get length aln for sequence type \"$seqType\". " . 
 
751
                     "Valid types are 'query', 'hit', 'sbjct' ('sbjct' = 'hit')");
 
752
    }                
 
753
    return $data;
 
754
}    
 
755
 
 
756
=head2 gaps
 
757
 
 
758
 Usage     : $hit_object->gaps( [seq_type] );
 
759
 Purpose   : Get the number of gaps in the aligned query, hit, or both sequences.
 
760
           : Data is summed across all HSPs.
 
761
 Example   : $qgaps = $hit_object->gaps('query');
 
762
           : $hgaps = $hit_object->gaps('hit');
 
763
           : $tgaps = $hit_object->gaps();    # default = total (query + hit)
 
764
 Returns   : scalar context: integer
 
765
           : array context without args: two-element list of integers  
 
766
           :    (queryGaps, hitGaps)
 
767
           : Array context can be forced by providing an argument of 'list' or 'array'.
 
768
           :
 
769
           : CAUTION: Calling this method within printf or sprintf is arrray context.
 
770
           : So this function may not give you what you expect. For example:
 
771
           :          printf "Total gaps: %d", $hit->gaps();
 
772
           : Actually returns a two-element array, so what gets printed 
 
773
           : is the number of gaps in the query, not the total
 
774
           :
 
775
 Argument  : seq_type: 'query' | 'hit' or 'sbjct' | 'total' | 'list'  (default = 'total')
 
776
             ('sbjct' is synonymous with 'hit')
 
777
 Throws    : n/a
 
778
 Comments  : If you need data for each HSP, use hsps() and then interate
 
779
           : through each HSP object.
 
780
           : This method requires that all HSPs be tiled. If they have not
 
781
           : already been tiled, they will be tiled first automatically..
 
782
           : Not relying on wantarray since that will fail in situations 
 
783
           : such as printf "%d", $hit->gaps() in which you might expect to 
 
784
           : be printing the total gaps, but evaluates to array context.
 
785
 
 
786
See Also   : L<length_aln()|length_aln>
 
787
 
 
788
=cut
 
789
 
 
790
#----------
 
791
sub gaps {
 
792
#----------
 
793
    my( $self, $seqType, $num ) = @_;
 
794
 
 
795
    $seqType ||= (wantarray ? 'list' : 'total');
 
796
    $seqType = 'sbjct' if $seqType eq 'hit';
 
797
 
 
798
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
799
 
 
800
    $seqType = lc($seqType);
 
801
 
 
802
    if( defined $num ) {
 
803
        $self->throw("Can't set gaps for seqType '$seqType'. Must be 'query' or 'hit'\n") unless ($seqType eq 'sbjct' or $seqType eq 'query');
 
804
 
 
805
        return $self->{'_gaps_'.$seqType} = $num;
 
806
    }
 
807
    elsif($seqType =~ /list|array/i) {
 
808
        return ($self->{'_gaps_query'}, $self->{'_gaps_sbjct'});
 
809
    }
 
810
    elsif($seqType eq 'total') {
 
811
        return ($self->{'_gaps_query'} + $self->{'_gaps_sbjct'}) || 0;
 
812
    } else {
 
813
        return $self->{'_gaps_'.$seqType} || 0;
 
814
    }
 
815
}    
 
816
 
 
817
 
 
818
=head2 matches
 
819
 
 
820
See documentation in L<Bio::Search::Hit::HitI::matches()|Bio::Search::Hit::HitI>
 
821
 
 
822
=cut
 
823
 
 
824
#---------------
 
825
sub matches {
 
826
#---------------
 
827
    my( $self, $arg1, $arg2) = @_;
 
828
    my(@data,$data);
 
829
 
 
830
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
831
 
 
832
    unless( $arg1 ) {
 
833
        @data = ($self->{'_totalIdentical'}, $self->{'_totalConserved'});
 
834
 
 
835
        return @data if @data;
 
836
    } else {
 
837
 
 
838
        if( defined $arg2 ) {
 
839
            $self->{'_totalIdentical'} = $arg1;
 
840
            $self->{'_totalConserved'} = $arg2;
 
841
            return ( $arg1, $arg2 );
 
842
        }
 
843
        elsif($arg1 =~ /^id/i) { 
 
844
            $data = $self->{'_totalIdentical'};
 
845
        } else {
 
846
            $data = $self->{'_totalConserved'};
 
847
        }
 
848
        return $data if $data;
 
849
    }
 
850
    
 
851
    ## Something went wrong if we make it to here.
 
852
    $self->throw("Can't get identical or conserved data: no data.");
 
853
}
 
854
 
 
855
 
 
856
=head2 start
 
857
 
 
858
 Usage     : $sbjct->start( [seq_type] );
 
859
 Purpose   : Gets the start coordinate for the query, sbjct, or both sequences
 
860
           : in the BlastHit object. If there is more than one HSP, the lowest start
 
861
           : value of all HSPs is returned.
 
862
 Example   : $qbeg = $sbjct->start('query');
 
863
           : $sbeg = $sbjct->start('hit');
 
864
           : ($qbeg, $sbeg) = $sbjct->start();
 
865
 Returns   : scalar context: integer 
 
866
           : array context without args: list of two integers (queryStart, sbjctStart)
 
867
           : Array context can be "induced" by providing an argument of 'list' or 'array'.
 
868
 Argument  : In scalar context: seq_type = 'query' or 'hit' or 'sbjct' (default = 'query')
 
869
             ('sbjct' is synonymous with 'hit')
 
870
 Throws    : n/a
 
871
 Comments  : This method requires that all HSPs be tiled. If there is more than one
 
872
           : HSP and they have not already been tiled, they will be tiled first automatically..
 
873
           : Remember that the start and end coordinates of all HSPs are 
 
874
           : normalized so that start < end. Strand information can be
 
875
           : obtained by calling $hit->strand().
 
876
 
 
877
See Also   : L<end()|end>, L<range()|range>, L<strand()|strand>, 
 
878
             L<Bio::Search::HSP::BlastHSP::start|Bio::Search::HSP::BlastHSP>
 
879
 
 
880
=cut
 
881
 
 
882
#----------
 
883
sub start {
 
884
#----------
 
885
    my ($self, $seqType, $num) = @_;
 
886
 
 
887
    $seqType ||= (wantarray ? 'list' : 'query');
 
888
    $seqType = 'sbjct' if $seqType eq 'hit';
 
889
 
 
890
    if( defined $num ) {
 
891
        $seqType = "_\L$seqType\E";
 
892
        return $self->{$seqType.'Start'} = $num;
 
893
    }
 
894
 
 
895
    # If there is only one HSP, defer this call to the solitary HSP.
 
896
    if($self->num_hsps == 1) {
 
897
        return $self->hsp->start($seqType);
 
898
    } else {
 
899
        &Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
900
        if($seqType =~ /list|array/i) {
 
901
            return ($self->{'_queryStart'}, $self->{'_sbjctStart'});
 
902
        } else {
 
903
            ## Sensitive to member name changes.
 
904
            $seqType = "_\L$seqType\E";
 
905
            return $self->{$seqType.'Start'};
 
906
        }
 
907
    }
 
908
}
 
909
 
 
910
 
 
911
=head2 end
 
912
 
 
913
 Usage     : $sbjct->end( [seq_type] );
 
914
 Purpose   : Gets the end coordinate for the query, sbjct, or both sequences
 
915
           : in the BlastHit object. If there is more than one HSP, 
 
916
             the largest end
 
917
           : value of all HSPs is returned.
 
918
 Example   : $qend = $sbjct->end('query');
 
919
           : $send = $sbjct->end('hit');
 
920
           : ($qend, $send) = $sbjct->end();
 
921
 Returns   : scalar context: integer
 
922
           : array context without args: list of two integers 
 
923
           : (queryEnd, sbjctEnd)
 
924
           : Array context can be "induced" by providing an argument 
 
925
           : of 'list' or 'array'.
 
926
 Argument  : In scalar context: seq_type = 'query' or 'sbjct'
 
927
           :  (case insensitive). If not supplied, 'query' is used.
 
928
 Throws    : n/a
 
929
 Comments  : This method requires that all HSPs be tiled. If there is 
 
930
           : more than one HSP and they have not already been tiled, 
 
931
           : they will be tiled first automatically..
 
932
           : Remember that the start and end coordinates of all HSPs are 
 
933
           : normalized so that start < end. Strand information can be
 
934
           : obtained by calling $hit->strand().
 
935
 
 
936
See Also   : L<start()|start>, L<range()|range>, L<strand()|strand>
 
937
 
 
938
=cut
 
939
 
 
940
#----------
 
941
sub end {
 
942
#----------
 
943
    my ($self, $seqType, $num) = @_;
 
944
 
 
945
    $seqType ||= (wantarray ? 'list' : 'query');
 
946
    $seqType = 'sbjct' if $seqType eq 'hit';
 
947
 
 
948
    if( defined $num ) {
 
949
        $seqType = "_\L$seqType\E";
 
950
        return $self->{$seqType.'Stop'} = $num;
 
951
    }
 
952
 
 
953
    # If there is only one HSP, defer this call to the solitary HSP.
 
954
    if($self->num_hsps == 1) {
 
955
        return $self->hsp->end($seqType);
 
956
    } else {
 
957
        Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
958
        if($seqType =~ /list|array/i) {
 
959
            return ($self->{'_queryStop'}, $self->{'_sbjctStop'});
 
960
        } else {
 
961
            ## Sensitive to member name changes.
 
962
            $seqType = "_\L$seqType\E";
 
963
            return $self->{$seqType.'Stop'};
 
964
        }
 
965
    }
 
966
}
 
967
 
 
968
=head2 range
 
969
 
 
970
 Usage     : $sbjct->range( [seq_type] );
 
971
 Purpose   : Gets the (start, end) coordinates for the query or sbjct sequence
 
972
           : in the HSP alignment.
 
973
 Example   : ($qbeg, $qend) = $sbjct->range('query');
 
974
           : ($sbeg, $send) = $sbjct->range('hit');
 
975
 Returns   : Two-element array of integers 
 
976
 Argument  : seq_type = string, 'query' or 'hit' or 'sbjct'  (default = 'query')
 
977
             ('sbjct' is synonymous with 'hit')
 
978
 Throws    : n/a
 
979
 
 
980
See Also   : L<start()|start>, L<end()|end>
 
981
 
 
982
=cut
 
983
 
 
984
#----------
 
985
sub range {
 
986
#----------
 
987
    my ($self, $seqType) = @_;
 
988
    $seqType ||= 'query';
 
989
    $seqType = 'sbjct' if $seqType eq 'hit';
 
990
    return ($self->start($seqType), $self->end($seqType));
 
991
}
 
992
 
 
993
 
 
994
=head2 frac_identical
 
995
 
 
996
 Usage     : $hit_object->frac_identical( [seq_type] );
 
997
 Purpose   : Get the overall fraction of identical positions across all HSPs.
 
998
           : The number refers to only the aligned regions and does not
 
999
           : account for unaligned regions in between the HSPs, if any.
 
1000
 Example   : $frac_iden = $hit_object->frac_identical('query');
 
1001
 Returns   : Float (2-decimal precision, e.g., 0.75).
 
1002
 Argument  : seq_type: 'query' | 'hit' or 'sbjct' | 'total'
 
1003
           : default = 'query' (but see comments below).
 
1004
           : ('sbjct' is synonymous with 'hit')
 
1005
 Throws    : n/a
 
1006
 Comments  : Different versions of Blast report different values for the total
 
1007
           : length of the alignment. This is the number reported in the
 
1008
           : denominators in the stats section:
 
1009
           : "Identical = 34/120 Positives = 67/120".
 
1010
           : NCBI BLAST uses the total length of the alignment (with gaps)
 
1011
           : WU-BLAST uses the length of the query sequence (without gaps).
 
1012
           :
 
1013
           : Therefore, when called with an argument of 'total',
 
1014
           : this method will report different values depending on the
 
1015
           : version of BLAST used. Total does NOT take into account HSP
 
1016
           : tiling, so it should not be used.
 
1017
           :
 
1018
           : To get the fraction identical among only the aligned residues,
 
1019
           : ignoring the gaps, call this method without an argument or 
 
1020
           : with an argument of 'query' or 'hit'.
 
1021
           :
 
1022
           : If you need data for each HSP, use hsps() and then iterate
 
1023
           : through the HSP objects.
 
1024
           : This method requires that all HSPs be tiled. If they have not
 
1025
           : already been tiled, they will be tiled first automatically.
 
1026
 
 
1027
See Also   : L<frac_conserved()|frac_conserved>, L<frac_aligned_query()|frac_aligned_query>, L<matches()|matches>, L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>
 
1028
 
 
1029
=cut
 
1030
 
 
1031
#------------------
 
1032
sub frac_identical {
 
1033
#------------------
 
1034
    my ($self, $seqType) = @_;
 
1035
    $seqType ||= 'query';
 
1036
    $seqType = 'sbjct' if $seqType eq 'hit';
 
1037
 
 
1038
    ## Sensitive to member name format.
 
1039
    $seqType = lc($seqType);
 
1040
 
 
1041
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1042
 
 
1043
    my $ident = $self->matches('id');
 
1044
    my $total = $self->length_aln($seqType);
 
1045
    my $ratio = $ident / $total;
 
1046
    my $ratio_rounded = sprintf( "%.3f", $ratio);
 
1047
 
 
1048
    # Round down iff normal rounding yields 1 (just like blast)
 
1049
    $ratio_rounded = 0.999 if (($ratio_rounded == 1) && ($ratio < 1));
 
1050
    return $ratio_rounded;
 
1051
}
 
1052
 
 
1053
 
 
1054
=head2 frac_conserved
 
1055
 
 
1056
 Usage     : $hit_object->frac_conserved( [seq_type] );
 
1057
 Purpose   : Get the overall fraction of conserved positions across all HSPs.
 
1058
           : The number refers to only the aligned regions and does not
 
1059
           : account for unaligned regions in between the HSPs, if any.
 
1060
 Example   : $frac_cons = $hit_object->frac_conserved('hit');
 
1061
 Returns   : Float (2-decimal precision, e.g., 0.75).
 
1062
 Argument  : seq_type: 'query' | 'hit' or 'sbjct' | 'total'
 
1063
           : default = 'query' (but see comments below).
 
1064
           : ('sbjct' is synonymous with 'hit')
 
1065
 Throws    : n/a
 
1066
 Comments  : Different versions of Blast report different values for the total
 
1067
           : length of the alignment. This is the number reported in the
 
1068
           : denominators in the stats section:
 
1069
           : "Positives = 34/120 Positives = 67/120".
 
1070
           : NCBI BLAST uses the total length of the alignment (with gaps)
 
1071
           : WU-BLAST uses the length of the query sequence (without gaps).
 
1072
           :
 
1073
           : Therefore, when called with an argument of 'total',
 
1074
           : this method will report different values depending on the
 
1075
           : version of BLAST used. Total does NOT take into account HSP
 
1076
           : tiling, so it should not be used.
 
1077
           :
 
1078
           : To get the fraction conserved among only the aligned residues,
 
1079
           : ignoring the gaps, call this method without an argument or 
 
1080
           : with an argument of 'query' or 'hit'.
 
1081
           :
 
1082
           : If you need data for each HSP, use hsps() and then interate
 
1083
           : through the HSP objects.
 
1084
           : This method requires that all HSPs be tiled. If they have not
 
1085
           : already been tiled, they will be tiled first automatically.
 
1086
 
 
1087
See Also   : L<frac_identical()|frac_identical>, L<matches()|matches>, L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>
 
1088
 
 
1089
=cut
 
1090
 
 
1091
#--------------------
 
1092
sub frac_conserved {
 
1093
#--------------------
 
1094
    my ($self, $seqType) = @_;
 
1095
    $seqType ||= 'query';
 
1096
    $seqType = 'sbjct' if $seqType eq 'hit';
 
1097
 
 
1098
    ## Sensitive to member name format.
 
1099
    $seqType = lc($seqType);
 
1100
 
 
1101
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1102
 
 
1103
    my $consv = $self->matches('cons');
 
1104
    my $total = $self->length_aln($seqType);
 
1105
    my $ratio = $consv / $total;
 
1106
    my $ratio_rounded = sprintf( "%.3f", $ratio);
 
1107
 
 
1108
    # Round down iff normal rounding yields 1 (just like blast)
 
1109
    $ratio_rounded = 0.999 if (($ratio_rounded == 1) && ($ratio < 1));
 
1110
    return $ratio_rounded;
 
1111
}
 
1112
 
 
1113
 
 
1114
 
 
1115
 
 
1116
=head2 frac_aligned_query
 
1117
 
 
1118
 Usage     : $hit_object->frac_aligned_query();
 
1119
 Purpose   : Get the fraction of the query sequence which has been aligned
 
1120
           : across all HSPs (not including intervals between non-overlapping
 
1121
           : HSPs).
 
1122
 Example   : $frac_alnq = $hit_object->frac_aligned_query();
 
1123
 Returns   : Float (2-decimal precision, e.g., 0.75).
 
1124
 Argument  : n/a
 
1125
 Throws    : n/a
 
1126
 Comments  : If you need data for each HSP, use hsps() and then interate
 
1127
           : through the HSP objects.
 
1128
           : To compute the fraction aligned, the logical length of the query
 
1129
           : sequence is used, meaning that for [T]BLASTX reports, the 
 
1130
           : full length of the query sequence is converted into amino acids
 
1131
           : by dividing by 3. This is necessary because of the way 
 
1132
           : the lengths of aligned sequences are computed.
 
1133
           : This method requires that all HSPs be tiled. If they have not
 
1134
           : already been tiled, they will be tiled first automatically.
 
1135
 
 
1136
See Also   : L<frac_aligned_hit()|frac_aligned_hit>, L<logical_length()|logical_length>, L<length_aln()|length_aln>,  L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>
 
1137
 
 
1138
=cut
 
1139
 
 
1140
#----------------------
 
1141
sub frac_aligned_query {
 
1142
#----------------------
 
1143
    my $self = shift;
 
1144
 
 
1145
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1146
    sprintf( "%.2f", $self->length_aln('query') / 
 
1147
             $self->logical_length('query'));
 
1148
}
 
1149
 
 
1150
 
 
1151
 
 
1152
=head2 frac_aligned_hit
 
1153
 
 
1154
 Usage     : $hit_object->frac_aligned_hit();
 
1155
 Purpose   : Get the fraction of the hit (sbjct) sequence which has been aligned
 
1156
           : across all HSPs (not including intervals between non-overlapping
 
1157
           : HSPs).
 
1158
 Example   : $frac_alnq = $hit_object->frac_aligned_hit();
 
1159
 Returns   : Float (2-decimal precision, e.g., 0.75).
 
1160
 Argument  : n/a
 
1161
 Throws    : n/a
 
1162
 Comments  : If you need data for each HSP, use hsps() and then interate
 
1163
           : through the HSP objects.
 
1164
           : To compute the fraction aligned, the logical length of the sbjct
 
1165
           : sequence is used, meaning that for TBLAST[NX] reports, the 
 
1166
           : full length of the sbjct sequence is converted into amino acids
 
1167
           : by dividing by 3. This is necessary because of the way 
 
1168
           : the lengths of aligned sequences are computed.
 
1169
           : This method requires that all HSPs be tiled. If they have not
 
1170
           : already been tiled, they will be tiled first automatically.
 
1171
 
 
1172
See Also   : L<frac_aligned_query()|frac_aligned_query>, L<matches()|matches>, , L<logical_length()|logical_length>, L<length_aln()|length_aln>,  L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>
 
1173
 
 
1174
=cut
 
1175
 
 
1176
#--------------------
 
1177
sub frac_aligned_hit {
 
1178
#--------------------
 
1179
    my $self = shift;
 
1180
 
 
1181
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1182
    sprintf( "%.2f", $self->length_aln('sbjct') / $self->logical_length('sbjct'));
 
1183
}
 
1184
 
 
1185
 
 
1186
## These methods are being maintained for backward compatibility. 
 
1187
 
 
1188
=head2 frac_aligned_sbjct
 
1189
 
 
1190
Same as L<frac_aligned_hit()|frac_aligned_hit>
 
1191
 
 
1192
=cut
 
1193
 
 
1194
*frac_aligned_sbjct = \&fract_aligned_hit;
 
1195
 
 
1196
=head2 num_unaligned_sbjct
 
1197
 
 
1198
Same as L<num_unaligned_hit()|num_unaligned_hit>
 
1199
 
 
1200
=cut
 
1201
 
 
1202
*num_unaligned_sbjct = \&num_unaligned_hit;
 
1203
 
 
1204
 
 
1205
=head2 num_unaligned_hit
 
1206
 
 
1207
 Usage     : $hit_object->num_unaligned_hit();
 
1208
 Purpose   : Get the number of the unaligned residues in the hit sequence.
 
1209
           : Sums across all all HSPs.
 
1210
 Example   : $num_unaln = $hit_object->num_unaligned_hit();
 
1211
 Returns   : Integer
 
1212
 Argument  : n/a
 
1213
 Throws    : n/a
 
1214
 Comments  : See notes regarding logical lengths in the comments for frac_aligned_hit().
 
1215
           : They apply here as well.
 
1216
           : If you need data for each HSP, use hsps() and then interate
 
1217
           : through the HSP objects.
 
1218
           : This method requires that all HSPs be tiled. If they have not
 
1219
           : already been tiled, they will be tiled first automatically..
 
1220
 
 
1221
See Also   : L<num_unaligned_query()|num_unaligned_query>,  L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>, L<frac_aligned_hit()|frac_aligned_hit>
 
1222
 
 
1223
=cut
 
1224
 
 
1225
#---------------------
 
1226
sub num_unaligned_hit {
 
1227
#---------------------
 
1228
    my $self = shift;
 
1229
 
 
1230
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1231
 
 
1232
    my $num = $self->logical_length('sbjct') - $self->length_aln('sbjct');
 
1233
    ($num < 0 ? 0 : $num );
 
1234
}
 
1235
 
 
1236
 
 
1237
=head2 num_unaligned_query
 
1238
 
 
1239
 Usage     : $hit_object->num_unaligned_query();
 
1240
 Purpose   : Get the number of the unaligned residues in the query sequence.
 
1241
           : Sums across all all HSPs.
 
1242
 Example   : $num_unaln = $hit_object->num_unaligned_query();
 
1243
 Returns   : Integer
 
1244
 Argument  : n/a
 
1245
 Throws    : n/a
 
1246
 Comments  : See notes regarding logical lengths in the comments for frac_aligned_query().
 
1247
           : They apply here as well.
 
1248
           : If you need data for each HSP, use hsps() and then interate
 
1249
           : through the HSP objects.
 
1250
           : This method requires that all HSPs be tiled. If they have not
 
1251
           : already been tiled, they will be tiled first automatically..
 
1252
 
 
1253
See Also   : L<num_unaligned_hit()|num_unaligned_hit>, L<frac_aligned_query()|frac_aligned_query>,  L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>
 
1254
 
 
1255
=cut
 
1256
 
 
1257
#-----------------------
 
1258
sub num_unaligned_query {
 
1259
#-----------------------
 
1260
    my $self = shift;
 
1261
 
 
1262
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1263
 
 
1264
    my $num = $self->logical_length('query') - $self->length_aln('query');
 
1265
    ($num < 0 ? 0 : $num );
 
1266
}
 
1267
 
 
1268
 
 
1269
 
 
1270
=head2 seq_inds
 
1271
 
 
1272
 Usage     : $hit->seq_inds( seq_type, class, collapse );
 
1273
 Purpose   : Get a list of residue positions (indices) across all HSPs
 
1274
           : for identical or conserved residues in the query or sbjct sequence.
 
1275
 Example   : @s_ind = $hit->seq_inds('query', 'identical');
 
1276
           : @h_ind = $hit->seq_inds('hit', 'conserved');
 
1277
           : @h_ind = $hit->seq_inds('hit', 'conserved', 1);
 
1278
 Returns   : Array of integers 
 
1279
           : May include ranges if collapse is non-zero.
 
1280
 Argument  : [0] seq_type  = 'query' or 'hit' or 'sbjct'  (default = 'query')
 
1281
           :                 ('sbjct' is synonymous with 'hit')
 
1282
           : [1] class     = 'identical' or 'conserved' (default = 'identical')
 
1283
           :              (can be shortened to 'id' or 'cons')
 
1284
           :              (actually, anything not 'id' will evaluate to 'conserved').
 
1285
           : [2] collapse  = boolean, if non-zero, consecutive positions are merged
 
1286
           :             using a range notation, e.g., "1 2 3 4 5 7 9 10 11" 
 
1287
           :             collapses to "1-5 7 9-11". This is useful for 
 
1288
           :             consolidating long lists. Default = no collapse.
 
1289
 Throws    : n/a.
 
1290
 
 
1291
See Also   : L<Bio::Search::HSP::BlastHSP::seq_inds()|Bio::Search::HSP::BlastHSP>
 
1292
 
 
1293
=cut
 
1294
 
 
1295
#-------------
 
1296
sub seq_inds {
 
1297
#-------------
 
1298
    my ($self, $seqType, $class, $collapse) = @_;
 
1299
 
 
1300
    $seqType  ||= 'query';
 
1301
    $class ||= 'identical';
 
1302
    $collapse ||= 0;
 
1303
 
 
1304
    $seqType = 'sbjct' if $seqType eq 'hit';
 
1305
 
 
1306
    my (@inds, $hsp);
 
1307
    foreach $hsp ($self->hsps) {
 
1308
        # This will merge data for all HSPs together.
 
1309
        push @inds, $hsp->seq_inds($seqType, $class);
 
1310
    }
 
1311
    
 
1312
    # Need to remove duplicates and sort the merged positions.
 
1313
    if(@inds) {
 
1314
        my %tmp = map { $_, 1 } @inds;
 
1315
        @inds = sort {$a <=> $b} keys %tmp;
 
1316
    }
 
1317
 
 
1318
    $collapse ?  &Bio::Search::SearchUtils::collapse_nums(@inds) : @inds; 
 
1319
}
 
1320
 
 
1321
 
 
1322
=head2 strand
 
1323
 
 
1324
See documentation in L<Bio::Search::Hit::HitI::strand()|Bio::Search::Hit::HitI>
 
1325
 
 
1326
=cut
 
1327
 
 
1328
#----------'
 
1329
sub strand {
 
1330
#----------
 
1331
    my ($self, $seqType, $strnd) = @_;
 
1332
 
 
1333
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1334
 
 
1335
    $seqType ||= (wantarray ? 'list' : 'query');
 
1336
    $seqType = 'sbjct' if $seqType eq 'hit';
 
1337
 
 
1338
    $seqType = lc($seqType);
 
1339
 
 
1340
    if( defined $strnd ) {
 
1341
        $self->throw("Can't set strand for seqType '$seqType'. Must be 'query' or 'hit'\n") unless ($seqType eq 'sbjct' or $seqType eq 'query');
 
1342
 
 
1343
        return $self->{'_strand_'.$seqType} = $strnd;
 
1344
    }
 
1345
 
 
1346
    my ($qstr, $hstr);
 
1347
    # If there is only one HSP, defer this call to the solitary HSP.
 
1348
    if($self->num_hsps == 1) {
 
1349
        return $self->hsp->strand($seqType);
 
1350
    } 
 
1351
    elsif( defined $self->{'_strand_query'}) {
 
1352
        # Get the data computed during hsp tiling.
 
1353
        $qstr = $self->{'_strand_query'};
 
1354
        $hstr = $self->{'_strand_sbjct'}
 
1355
    }
 
1356
    else {
 
1357
        # otherwise, iterate through all HSPs collecting strand info.
 
1358
        # This will return the string "-1/1" if there are HSPs on different strands.
 
1359
        # NOTE: This was the pre-10/21/02 procedure which will no longer be used,
 
1360
        # (unless the above elsif{} is commented out).
 
1361
        my (%qstr, %hstr);
 
1362
        foreach my $hsp( $self->hsps ) {
 
1363
            my ( $q, $h ) = $hsp->strand();
 
1364
            $qstr{ $q }++;
 
1365
            $hstr{ $h }++;
 
1366
        }
 
1367
        $qstr = join( '/', sort keys %qstr);
 
1368
        $hstr = join( '/', sort keys %hstr);
 
1369
    }
 
1370
 
 
1371
    if($seqType =~ /list|array/i) {
 
1372
        return ($qstr, $hstr);
 
1373
    } elsif( $seqType eq 'query' ) {
 
1374
        return $qstr;
 
1375
    } else {
 
1376
        return $hstr;
 
1377
    }
 
1378
}
 
1379
 
 
1380
=head2 frame
 
1381
 
 
1382
See documentation in L<Bio::Search::Hit::HitI::frame()|Bio::Search::Hit::HitI>
 
1383
 
 
1384
=cut
 
1385
 
 
1386
#----------'
 
1387
sub frame { 
 
1388
#----------
 
1389
    my( $self, $frm ) = @_;
 
1390
 
 
1391
    Bio::Search::SearchUtils::tile_hsps($self) unless $self->tiled_hsps;
 
1392
 
 
1393
    if( defined $frm ) {
 
1394
        return $self->{'_frame'} = $frm;
 
1395
    }
 
1396
 
 
1397
    # The check for $self->{'_frame'} is a remnant from the 'query' mode days
 
1398
    # in which the sbjct object would collect data from the description line only.
 
1399
 
 
1400
    my ($frame);
 
1401
    if(not defined($self->{'_frame'})) {
 
1402
        $frame = $self->hsp->frame;
 
1403
    } else {
 
1404
        $frame = $self->{'_frame'}; 
 
1405
    } 
 
1406
    return $frame;
 
1407
}
 
1408
 
 
1409
=head2 rank
 
1410
 
 
1411
 Title   : rank
 
1412
 Usage   : $obj->rank($newval)
 
1413
 Function: Get/Set the rank of this Hit in the Query search list
 
1414
           i.e. this is the Nth hit for a specific query
 
1415
 Returns : value of rank
 
1416
 Args    : newvalue (optional)
 
1417
 
 
1418
 
 
1419
=cut
 
1420
 
 
1421
sub rank{
 
1422
    my $self = shift;
 
1423
    return $self->{'_rank'} = shift if @_;
 
1424
    return $self->{'_rank'} || 1;
 
1425
}
 
1426
 
 
1427
=head2 locus
 
1428
 
 
1429
 Title   : locus
 
1430
 Usage   : $locus = $hit->locus();
 
1431
 Function: Retrieve the locus (if available) for the hit
 
1432
 Returns : a scalar string (empty string if not set)
 
1433
 Args    : none
 
1434
 
 
1435
=cut
 
1436
 
 
1437
sub locus {
 
1438
    my ($self,$value) = @_;
 
1439
    my $previous = $self->{'_locus'};
 
1440
    if( defined $value || ! defined $previous ) { 
 
1441
      unless (defined $value) {
 
1442
        if ($self->{'_name'} =~/(gb|emb|dbj|ref)\|(.*)\|(.*)/) {
 
1443
                  $value = $previous = $3;
 
1444
                } else {
 
1445
          $value = $previous = '';
 
1446
        }
 
1447
      }
 
1448
          $self->{'_locus'} = $value;
 
1449
    } 
 
1450
        return $previous;
 
1451
}
 
1452
 
 
1453
=head2 each_accession_number
 
1454
 
 
1455
 Title   : each_accession_number
 
1456
 Usage   : @each_accession_number = $hit->each_accession_number();
 
1457
 Function: Get each accession number listed in the description of the hit.
 
1458
           If there are no alternatives, then only the primary accession will 
 
1459
           be given
 
1460
 Returns : list of all accession numbers in the description
 
1461
 Args    : none
 
1462
 
 
1463
=cut
 
1464
 
 
1465
sub each_accession_number {
 
1466
    my ($self,$value) = @_;
 
1467
    my $desc = $self->{'_description'};
 
1468
    #put primary accnum on the list
 
1469
    my @accnums;
 
1470
    push (@accnums,$self->{'_accession'});
 
1471
    if( defined $desc )  { 
 
1472
      while ($desc =~ /(\b\S+\|\S*\|\S*\s?)/g) {
 
1473
        my $id = $1;
 
1474
        my ($acc, $version);
 
1475
        if ($id =~ /(gb|emb|dbj|sp|pdb|bbs|ref|tp[gde])\|(.*)\|(.*)/) {
 
1476
            ($acc, $version) = split /\./, $2; 
 
1477
        } elsif ($id =~ /(pir|prf|pat|gnl)\|(.*)\|(.*)/) {
 
1478
            ($acc, $version) = split /\./, $3;  
 
1479
        } elsif( $id =~ /(gim|gi|bbm|bbs|lcl)\|(\d*)/) {
 
1480
            $acc = $id;
 
1481
        } elsif( $id =~ /(oth)\|(.*)\|(.*)\|(.*)/ ) { # discontinued...
 
1482
            ($acc,$version) = ($2);
 
1483
        } else {
 
1484
                     #punt, not matching the db's at ftp://ftp.ncbi.nih.gov/blast/db/README
 
1485
                     #Database Name                     Identifier Syntax
 
1486
          #============================      ========================
 
1487
          #GenBank                           gb|accession|locus
 
1488
          #EMBL Data Library                 emb|accession|locus
 
1489
          #DDBJ, DNA Database of Japan       dbj|accession|locus
 
1490
          #NBRF PIR                          pir||entry
 
1491
          #Protein Research Foundation       prf||name
 
1492
          #SWISS-PROT                        sp|accession|entry name
 
1493
          #Brookhaven Protein Data Bank      pdb|entry|chain
 
1494
          #Patents                           pat|country|number 
 
1495
          #GenInfo Backbone Id               bbs|number 
 
1496
          #General database identifier           gnl|database|identifier
 
1497
          #NCBI Reference Sequence           ref|accession|locus
 
1498
          #Local Sequence identifier         lcl|identifier
 
1499
              $acc=$id;
 
1500
            }
 
1501
            push(@accnums, $acc);
 
1502
          }
 
1503
    }  
 
1504
    return @accnums;
 
1505
}
 
1506
 
 
1507
=head2 tiled_hsps
 
1508
 
 
1509
See documentation in L<Bio::Search::SearchUtils::tile_hsps()|Bio::Search::SearchUtils>
 
1510
 
 
1511
=cut
 
1512
 
 
1513
sub tiled_hsps { 
 
1514
    my $self = shift;
 
1515
    return $self->{'_tiled_hsps'} = shift if @_;
 
1516
    return $self->{'_tiled_hsps'};
 
1517
}
 
1518
 
 
1519
=head2 query_length
 
1520
 
 
1521
 Title   : query_length
 
1522
 Usage   : $obj->query_length($newval)
 
1523
 Function: Get/Set the query_length
 
1524
 Returns : value of query_length (a scalar)
 
1525
 Args    : on set, new value (a scalar or undef, optional)
 
1526
 
 
1527
 
 
1528
=cut
 
1529
 
 
1530
sub query_length{
 
1531
    my $self = shift;
 
1532
 
 
1533
    return $self->{'_query_length'} = shift if @_;
 
1534
    return $self->{'_query_length'};
 
1535
}
 
1536
 
 
1537
 
343
1538
1;