File Coverage

Bio/Search/HSP/PullHSPI.pm
Criterion Covered Total %
statement 137 200 68.5
branch 59 118 50.0
condition 17 56 30.3
subroutine 28 32 87.5
pod 28 28 100.0
total 269 434 61.9


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::Search::HSP::PullHSPI
3             #
4             # Please direct questions and support issues to
5             #
6             # Cared for by Sendu Bala
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,
66             L
67             L
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
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 2     2   8 use strict;
  2         3  
  2         48  
127              
128 2     2   7 use base qw(Bio::Search::HSP::HSPI Bio::PullParserI);
  2         1  
  2         927  
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 179     179   257 my ($self, @args) = @_;
143            
144             # fields most subclasses probably will want
145 179         2142 $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 179         497 my ($parent, $chunk, $hsp_data) = $self->_rearrange([qw(PARENT
174             CHUNK
175             HSP_DATA)], @args);
176            
177 179 50 33     427 $self->throw("Need -parent or -chunk to be defined") unless defined $parent || $chunk;
178            
179 179 50       511 $self->parent($parent) if $parent;
180            
181 179 100       274 if ($chunk) {
182 145         142 my ($io, $start, $end) = (undef, 0, undef);
183 145 50       233 if (ref($chunk) eq 'ARRAY') {
184 145         137 ($io, $start, $end) = @{$chunk};
  145         268  
185             }
186             else {
187 0         0 $io = $chunk;
188             }
189 145         397 $self->chunk($io, -start => $start, -end => $end);
190             }
191            
192 179 100       312 $self->_raw_hsp_data($hsp_data) if $hsp_data;
193            
194 179         270 return $self;
195             }
196              
197             sub _raw_hsp_data {
198 68     68   45 my $self = shift;
199 68 100       91 if (@_) {
200 34         52 $self->{_raw_hsp_data} = shift;
201             }
202 68         65 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 4     4 1 12 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 0     0 1 0 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 16     16 1 49 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 18     18 1 2338 my ($self, $type) = @_;
270            
271 18 50       47 $type = lc $type if defined $type;
272 18 50 33     84 $type = 'hit' if (defined $type && $type =~ /subject|sbjct/);
273 18 100 66     154 $type = 'total' if (! defined $type || $type eq 'hsp' || $type !~ /query|hit|subject|sbjct|total/);
      66        
274            
275 18         33 my $ratio = $self->num_identical($type) / $self->length($type);
276 18         219 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 0     0 1 0 my ($self, $type) = @_;
296            
297 0 0       0 $type = lc $type if defined $type;
298 0 0 0     0 $type = 'hit' if (defined $type && $type =~ /subject|sbjct/);
299 0 0 0     0 $type = 'total' if (! defined $type || $type eq 'hsp' || $type !~ /query|hit|subject|sbjct|total/);
      0        
300            
301 0         0 my $ratio = $self->num_conserved($type) / $self->length($type);
302 0         0 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 20     20 1 25 my $self = shift;
317 20         38 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 2     2 1 814 my $self = shift;
332 2         6 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 0     0 1 0 my ($self, $type) = @_;
347 0 0       0 $type = lc $type if defined $type;
348 0 0 0     0 $type = 'total' if (! defined $type || $type eq 'hsp' || $type !~ /query|hit|subject|sbjct|total/);
      0        
349 0 0       0 $type = 'hit' if $type =~ /sbjct|subject/;
350            
351 0 0       0 if ($type eq 'total') {
352 0         0 return scalar($self->seq_inds('hit', 'gap')) + scalar($self->seq_inds('query', 'gap'));
353             }
354 0         0 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 14     14 1 960 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 17     17 1 785 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 11     11 1 790 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 36     36 1 49 my ($self, $type) = @_;
420 36 50       62 $type = 'total' unless defined $type;
421 36         45 $type = lc $type;
422              
423 36 100       137 if ($type =~ /^q/i) {
    100          
424 10         30 return $self->query->length;
425             }
426             elsif ($type =~ /^(hit|subject|sbjct)/) {
427 10         28 return $self->hit->length;
428             }
429             else {
430 16         43 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 2     2 1 731 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 4     4 1 8 my ($self) = @_;
460 4         13 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 4     4 1 792 my $self = shift;
475            
476 4         23 require Bio::LocatableSeq;
477 4         1826 require Bio::SimpleAlign;
478 4         22 my $aln = Bio::SimpleAlign->new();
479 4         17 my $hs = $self->seq('hit');
480 4         13 my $qs = $self->seq('query');
481 4 50 33     20 if ($hs && $qs) {
482 4         12 $aln->add_seq($hs);
483 4         7 $aln->add_seq($qs);
484 4         10 return $aln;
485             }
486 0         0 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, L
515              
516             =cut
517              
518             sub seq_inds {
519 164     164 1 349 my ($self, $seqType, $class, $collapse) = @_;
520            
521 164   50     272 $seqType ||= 'query';
522 164   50     243 $class ||= 'identical';
523 164   50     513 $collapse ||= 0;
524 164         207 $seqType = lc($seqType);
525 164         171 $class = lc($class);
526 164 50       294 $seqType = 'hit' if $seqType eq 'sbjct';
527            
528 164         312 my $t = substr($seqType,0,1);
529 164 100 33     422 if ($t eq 'q') {
    50          
530 102         150 $seqType = 'query';
531             }
532             elsif ($t eq 's' || $t eq 'h') {
533 62         73 $seqType = 'hit';
534             }
535             else {
536 0         0 $self->warn("unknown seqtype $seqType using 'query'");
537 0         0 $seqType = 'query';
538             }
539            
540 164         190 $t = substr($class,0,1);
541 164 100       454 if ($t eq 'c') {
    100          
    100          
    50          
542 71 100       116 if ($class eq 'conserved-not-identical') {
543 36         59 $class = 'conserved';
544             }
545             else {
546 35         48 $class = 'conservedall';
547             }
548             }
549             elsif ($t eq 'i') {
550 54         56 $class = 'identical';
551             }
552             elsif ($t eq 'n') {
553 35         42 $class = 'nomatch';
554             }
555             elsif ($t eq 'g') {
556 4         6 $class = 'gap';
557             }
558             else {
559 0         0 $self->warn("unknown sequence class $class using 'identical'");
560 0         0 $class = 'identical';
561             }
562            
563 164         243 $seqType .= '_';
564 164         183 $class .= '_inds';
565            
566 164         137 my @ary;
567 164 100       244 if ($class eq 'conservedall_inds') {
568 5504         4706 my %tmp = map { $_, 1 } @{$self->get_field($seqType.'conserved_inds')},
  35         93  
569 35         31 @{$self->get_field($seqType.'identical_inds')};
  35         102  
570 35         698 @ary = sort {$a <=> $b} keys %tmp;
  36700         21345  
571             }
572             else {
573 129         125 @ary = @{$self->get_field($seqType.$class)};
  129         412  
574             }
575            
576 164 50       3298 return $collapse ? &Bio::Search::SearchUtils::collapse_nums(@ary) : @ary;
577             }
578              
579             =head2 Inherited from L
580              
581             These methods come from L
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 4     4 1 10 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 8     8 1 29 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 8     8 1 1421 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 retreive both query and hit together
655              
656             =cut
657              
658             sub strand {
659 14     14 1 16 my $self = shift;
660 14         13 my $val = shift;
661 14 50       26 $val = 'query' unless defined $val;
662 14         29 $val =~ s/^\s+//;
663              
664 14 100       38 if ($val =~ /^q/i) {
    50          
    0          
665 6         11 return $self->query->strand(@_);
666             }
667             elsif ($val =~ /^hi|^s/i) {
668 8         19 return $self->hit->strand(@_);
669             }
670             elsif ($val =~ /^list|array/i) {
671 0         0 return ($self->query->strand(@_), $self->hit->strand(@_) );
672             }
673             else {
674 0         0 $self->warn("unrecognized component '$val' requested\n");
675             }
676 0         0 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 14     14 1 1424 my $self = shift;
693 14         12 my $val = shift;
694 14 0       22 $val = (wantarray ? 'list' : 'query') unless defined $val;
    50          
695 14         32 $val =~ s/^\s+//;
696              
697 14 100       42 if ($val =~ /^q/i) {
    50          
    0          
698 6         13 return $self->query->start(@_);
699             }
700             elsif ($val =~ /^(hi|s)/i) {
701 8         17 return $self->hit->start(@_);
702             }
703             elsif ($val =~ /^list|array/i) {
704 0         0 return ($self->query->start(@_), $self->hit->start(@_) );
705             }
706             else {
707 0         0 $self->warn("unrecognized component '$val' requested\n");
708             }
709 0         0 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 14     14 1 17 my $self = shift;
726 14         12 my $val = shift;
727 14 0       26 $val = (wantarray ? 'list' : 'query') unless defined $val;
    50          
728 14         29 $val =~ s/^\s+//;
729              
730 14 100       45 if ($val =~ /^q/i) {
    50          
    0          
731 6         12 return $self->query->end(@_);
732             }
733             elsif ($val =~ /^(hi|s)/i) {
734 8         15 return $self->hit->end(@_);
735             }
736             elsif ($val =~ /^list|array/i) {
737 0         0 return ($self->query->end(@_), $self->hit->end(@_) );
738             }
739             else {
740 0         0 $self->warn("unrecognized end component '$val' requested\n");
741             }
742 0         0 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 12     12 1 15 my ($self, $seqType) = @_;
759 12   50     21 $seqType ||= 'query';
760 12 50       23 $seqType = 'hit' if $seqType eq 'sbjct';
761 12 50       48 if ($seqType =~ /^(m|ho)/i ) {
762 0         0 $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 12   100     28 my $str = $self->seq_str($seqType) || return;
766 10         48 require Bio::LocatableSeq;
767 10 100       39 my $id = ($seqType =~ /^q/i) ? $self->query->seq_id : $self->hit->seq_id;
768 10 100       33 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, L, B<_set_match_seq()>
793              
794             =cut
795              
796             sub seq_str {
797 19     19 1 766 my $self = shift;
798 19   100     42 my $type = shift || 'query';
799              
800 19 100       77 if ($type =~ /^q/i) {
    100          
    50          
801 7         23 return $self->query_string(@_);
802             }
803             elsif ($type =~ /^(s)|(hi)/i) {
804 10         26 return $self->hit_string(@_);
805             }
806             elsif ($type =~ /^(ho)|(ma)/i) {
807 2         5 return $self->homology_string(@_);
808             }
809             else {
810 0         0 $self->warn("unknown sequence type $type");
811             }
812 0         0 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 4     4 1 12 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 0     0 1 0 my ($self, @args) = @_;
851 0         0 my($seqType, $beg, $end) = $self->_rearrange([qw(SEQ START STOP)], @args);
852 0   0     0 $seqType ||= 'query';
853 0 0       0 $seqType = 'hit' if $seqType eq 'sbjct';
854            
855 0         0 my @data;
856 0 0 0     0 if ((!defined $beg && !defined $end) || ! $self->seq_str('match')) {
857 0         0 push @data, ($self->num_identical, $self->num_conserved);
858             }
859             else {
860 0   0     0 $beg ||= 0;
861 0   0     0 $end ||= 0;
862 0         0 my ($start, $stop) = $self->range($seqType);
863            
864 0 0       0 if ($beg == 0) {
    0          
865 0         0 $beg = $start;
866 0         0 $end = $beg+$end;
867             }
868             elsif ($end == 0) {
869 0         0 $end = $stop;
870 0         0 $beg = $end-$beg;
871             }
872            
873 0 0       0 if ($end >= $stop) {
874 0         0 $end = $stop;
875             }
876             else {
877 0         0 $end += 1;
878             }
879 0 0       0 if ($beg < $start) {
880 0         0 $beg = $start;
881             }
882            
883 0         0 my $seq = substr($self->seq_str('homology'), $beg-$start, ($end-$beg));
884            
885 0 0       0 if (!CORE::length $seq) {
886 0         0 $self->throw("Undefined sub-sequence ($beg,$end). Valid range = $start - $stop");
887             }
888             ## Get data for a substring.
889 0         0 $seq =~ s/ //g; # remove space (no info).
890 0         0 my $len_cons = CORE::length $seq;
891 0         0 $seq =~ s/\+//g; # remove '+' characters (conservative substitutions)
892 0         0 my $len_id = CORE::length $seq;
893 0         0 push @data, ($len_id, $len_cons);
894             }
895            
896 0         0 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 4     4 1 1389 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 4     4 1 1412 my ($self, $seqType) = @_;
935 4   50     18 $seqType ||= 'query';
936            
937 4         5 my ($start, $end);
938 4 50       9 if ($seqType eq 'query') {
939 4         11 $start = $self->query->start;
940 4         9 $end = $self->query->end;
941             }
942             else {
943 0         0 $start = $self->hit->start;
944 0         0 $end = $self->hit->end;
945             }
946 4         12 return ($start, $end);
947             }
948              
949             #*** would want cigar stuff from GenericHSP - move to HSPI?
950              
951             1;
952