Mercurial > repos > willmclaren > ensembl_vep

comparison variant_effect_predictor/Bio/Tools/Blast/HSP.pm @ 0:21066c0abaf5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author willmclaren
date Fri, 03 Aug 2012 10:04:48 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:21066c0abaf5
1 #----------------------------------------------------------------------------
2 # PACKAGE : Bio::Tools::Blast::HSP
3 # AUTHOR : Steve Chervitz (sac@bioperl.org)
4 # CREATED : March 1996
5 # STATUS : Alpha
6 # REVISION: $Id: HSP.pm,v 1.18 2002/10/22 07:38:48 lapp Exp $
7 #
8 # For the latest version and documentation, visit the distribution site:
9 # http://genome-www.stanford.edu/perlOOP/bioperl/blast/
10 #
11 # To generate documentation, run this module through pod2html
12 # (preferably from Perl v5.004 or better).
13 #
14 # Copyright (c) 1996-2000 Steve Chervitz. All Rights Reserved.
15 # This module is free software; you can redistribute it and/or
16 # modify it under the same terms as Perl itself.
17 #----------------------------------------------------------------------------
18
19 package Bio::Tools::Blast::HSP;
20
21 use Bio::Root::Global qw(:devel);
22 use Bio::Root::Object ();
23 #use Bio::Root::Err qw(:std);
24
25 @ISA = qw( Bio::Root::Object);
26
27 use strict;
28 use vars qw($ID $GAP_SYMBOL @SCORE_CUTOFFS $Revision %STRAND_SYMBOL);
29 $ID = 'Bio::Tools::Blast::HSP';
30 $Revision = '$Id: HSP.pm,v 1.18 2002/10/22 07:38:48 lapp Exp $'; #'
31
32 $GAP_SYMBOL = '-'; # Need a more general way to handle gap symbols.
33 @SCORE_CUTOFFS = ( 100, 30 ); # Bit score cutoffs (see homol_score()).
34 %STRAND_SYMBOL = ('Plus' => 1, 'Minus' => -1);
35
36 ## POD Documentation:
37
38 =head1 NAME
39
40 Bio::Tools::Blast::HSP - Bioperl BLAST High-Scoring Segment Pair object
41
42 =head1 SYNOPSIS
43
44 =head2 Object Creation
45
46 The construction of HSP objects is handled by Bio::Tools::Blast:: Sbjct.pm.
47 You should not need to use this package directly. See L<_initialize()|_initialize>
48 for a description of constructor parameters.
49
50 require Bio::Tools::Blast::HSP;
51
52 $hspObj = eval{ new Bio::Tools::Blast::HSP(-DATA =>\@hspData,
53 -PARENT =>$sbjct_object,
54 -NAME =>$hspCount,
55 -PROGRAM =>'TBLASTN',
56 );
57 };
58
59 @hspData includes the raw BLAST report data for a specific HSP,
60 and is prepared by Bio::Tools::Blast::Sbjct.pm.
61
62 =head1 INSTALLATION
63
64 This module is included with the central Bioperl distribution:
65
66 http://bio.perl.org/Core/Latest
67 ftp://bio.perl.org/pub/DIST
68
69 Follow the installation instructions included in the README file.
70
71 =head1 DESCRIPTION
72
73 The Bio::Tools::Blast::HSP.pm module encapsulates data and methods for
74 manipulating, parsing, and analyzing HSPs ("High-scoring Segment Pairs")
75 derived from BLAST sequence analysis.
76
77 This module is a utility module used by the B<Bio::Tools::Blast::Sbjct.pm>
78 and is not intended for separate use. Please see documentation for
79 B<Bio::Tools::Blast.pm> for some basic information about using
80 HSP objects (L<Links:>).
81
82 =over 0
83
84 =item * Supports BLAST versions 1.x and 2.x, gapped and ungapped.
85
86 =back
87
88 Bio::Tools::Blast::HSP.pm has the ability to extract a list of all
89 residue indices for identical and conservative matches along both
90 query and sbjct sequences. Since this degree of detail is not always
91 needed, this behavior does not occur during construction of the HSP
92 object. These data will automatically be collected as necessary as
93 the HSP.pm object is used.
94
95 =head1 DEPENDENCIES
96
97 Bio::Tools::Blast::HSP.pm is a concrete class that inherits from
98 B<Bio::Root::Object.pm> and relies on B<Bio::Tools::Sbjct.pm> as a
99 container for HSP.pm objects. B<Bio::Seq.pm> and B<Bio::UnivAln.pm>
100 are employed for creating sequence and alignment objects,
101 respectively.
102
103
104 =head2 Relationship to UnivAln.pm & Seq.pm
105
106 HSP.pm can provide the query or sbjct sequence as a B<Bio::Seq.pm>
107 object via the L<seq()|seq> method. The HSP.pm object can also create a
108 two-sequence B<Bio::UnivAln.pm> alignment object using the the query
109 and sbjct sequences via the L<get_aln()|get_aln> method. Creation of alignment
110 objects is not automatic when constructing the HSP.pm object since
111 this level of functionality is not always required and would generate
112 a lot of extra overhead when crunching many reports.
113
114
115 =head1 FEEDBACK
116
117 =head2 Mailing Lists
118
119 User feedback is an integral part of the evolution of this and other
120 Bioperl modules. Send your comments and suggestions preferably to one
121 of the Bioperl mailing lists. Your participation is much appreciated.
122
123 bioperl-l@bioperl.org - General discussion
124 http://bio.perl.org/MailList.html - About the mailing lists
125
126 =head2 Reporting Bugs
127
128 Report bugs to the Bioperl bug tracking system to help us keep track
129 the bugs and their resolution. Bug reports can be submitted via email
130 or the web:
131
132 bioperl-bugs@bio.perl.org
133 http://bugzilla.bioperl.org/
134
135 =head1 AUTHOR
136
137 Steve Chervitz, E<lt>sac@bioperl.orgE<gt>
138
139 =head1 SEE ALSO
140
141 Bio::Tools::Blast::Sbjct.pm - Blast hit object.
142 Bio::Tools::Blast.pm - Blast object.
143 Bio::Seq.pm - Biosequence object
144 Bio::UnivAln.pm - Biosequence alignment object.
145 Bio::Root::Object.pm - Proposed base class for all Bioperl objects.
146
147 =head2 Links:
148
149 http://bio.perl.org/Core/POD/Tools/Blast/Sbjct.pm.html
150
151 http://bio.perl.org/Projects/modules.html - Online module documentation
152 http://bio.perl.org/Projects/Blast/ - Bioperl Blast Project
153 http://bio.perl.org/ - Bioperl Project Homepage
154
155 =head1 COPYRIGHT
156
157 Copyright (c) 1996-98 Steve Chervitz. All Rights Reserved.
158 This module is free software; you can redistribute it and/or
159 modify it under the same terms as Perl itself.
160
161 =cut
162
163
164
165 #
166 ##
167 ###
168 #### END of main POD documentation.
169 ###
170 ##
171 #
172
173 =head1 APPENDIX
174
175 Methods beginning with a leading underscore are considered private
176 and are intended for internal use by this module. They are
177 B<not> considered part of the public interface and are described here
178 for documentation purposes only.
179
180 =cut
181
182 #####################################################################################
183 ## CONSTRUCTOR ##
184 #####################################################################################
185
186 =head2 _initialize
187
188 Usage : n/a; automatically called by Bio::Root::Object::new()
189 : Bio::Tools::Blast::HSP.pm objects are constructed
190 : automatically by Bio::Tools::Sbjct.pm, so there is no need
191 : for direct consumption.
192 Purpose : Initializes HSP data and calls private methods to extract
193 : the data for a given HSP.
194 : Calls superclass constructor first (Bio::Root::Object.pm).
195 Returns : n/a
196 Argument : Named parameters passed from new():
197 : All tags must be uppercase (does not call _rearrange()).
198 : -DATA => array ref containing raw data for one HSP.
199 : -PARENT => Sbjct.pm object ref.
200 : -NAME => integer (1..n).
201 : -PROGRAM => string ('TBLASTN', 'BLASTP', etc.).
202
203 See Also : L<_set_data()|_set_data>, B<Bio::Root::Object::new()>, B<Bio::Tools::Blast::Sbjct::_set_hsps()>
204
205 =cut
206
207 #----------------
208 sub _initialize {
209 #----------------
210 my( $self, %param ) = @_;
211
212 $self->SUPER::_initialize( %param );
213
214 # The gapped and program booleans may be needed after the HSP object
215 # is built.
216 # $self->{'_gapped'} = $param{-GAPPED} || 0;
217 $self->{'_prog'} = $param{-PROGRAM} || 0;
218 $self->_set_data( @{$param{-DATA}} );
219 }
220
221 #####################################################################################
222 ## ACCESSORS ##
223 #####################################################################################
224
225
226 =head2 _set_data
227
228 Usage : n/a; called automatically during object construction.
229 Purpose : Sets the query sequence, sbjct sequence, and the "match" data
230 : which consists of the symbols between the query and sbjct lines
231 : in the alignment.
232 Argument : Array (all lines from a single, complete HSP, one line per element)
233 Throws : Propagates any exceptions from the methods called ("See Also")
234
235 See Also : L<_set_seq()|_set_seq>, L<_set_residues()|_set_residues>, L<_set_score_stats()|_set_score_stats>, L<_set_match_stats()|_set_match_stats>, L<_initialize()|_initialize>
236
237 =cut
238
239 #--------------
240 sub _set_data {
241 #--------------
242 my $self = shift;
243 my @data = @_;
244 my @queryList = (); # 'Query' = SEQUENCE USED TO QUERY THE DATABASE.
245 my @sbjctList = (); # 'Sbjct' = HOMOLOGOUS SEQUENCE FOUND IN THE DATABASE.
246 my @matchList = ();
247 my $matchLine = 0; # Alternating boolean: when true, load 'match' data.
248 my @linedat = ();
249
250 $DEBUG and print STDERR "$ID: set_data()\n";
251
252 my($line, $aln_row_len, $length_diff);
253 $length_diff = 0;
254
255 # Collecting data for all lines in the alignment
256 # and then storing the collections for possible processing later.
257 #
258 # Note that "match" lines may not be properly padded with spaces.
259 # This loop now properly handles such cases:
260 # Query: 1141 PSLVELTIRDCPRLEVGPMIRSLPKFPMLKKLDLAVANIIEEDLDVIGSLEELVIXXXXX 1200
261 # PSLVELTIRDCPRLEVGPMIRSLPKFPMLKKLDLAVANIIEEDLDVIGSLEELVI
262 # Sbjct: 1141 PSLVELTIRDCPRLEVGPMIRSLPKFPMLKKLDLAVANIIEEDLDVIGSLEELVILSLKL 1200
263
264 foreach $line( @data ) {
265 next if $line =~ /^\s*$/;
266
267 if( $line =~ /^ ?Score/ ) {
268 $self->_set_score_stats( $line );
269 } elsif( $line =~ /^ ?(Identities|Positives|Strand)/ ) {
270 $self->_set_match_stats( $line );
271 } elsif( $line =~ /^ ?Frame = ([\d+-]+)/ ) {
272 # Version 2.0.8 has Frame information on a separate line.
273 $self->{'_frame'} = $1;
274 } elsif( $line =~ /^(Query:?[\s\d]+)([^\s\d]+)/ ) {
275 push @queryList, $line;
276 $self->{'_match_indent'} = CORE::length $1;
277 $aln_row_len = (CORE::length $1) + (CORE::length $2);
278 $matchLine = 1;
279 } elsif( $matchLine ) {
280 # Pad the match line with spaces if necessary.
281 $length_diff = $aln_row_len - CORE::length $line;
282 $length_diff and $line .= ' 'x $length_diff;
283 push @matchList, $line;
284 $matchLine = 0;
285 } elsif( $line =~ /^Sbjct/ ) {
286 push @sbjctList, $line;
287 }
288 }
289
290 # Storing the query and sbjct lists in case they are needed later.
291 # We could make this conditional to save memory.
292 $self->{'_queryList'} = \@queryList;
293 $self->{'_sbjctList'} = \@sbjctList;
294
295 # Storing the match list in case it is needed later.
296 $self->{'_matchList'} = \@matchList;
297
298 if(not defined ($self->{'_numIdentical'})) {
299 $self->throw("Can't parse match statistics.",
300 "Possibly a new or unrecognized Blast format.");
301 }
302
303 if(!scalar @queryList or !scalar @sbjctList) {
304 $self->throw("Can't find query or sbjct alignment lines.",
305 "Possibly unrecognized Blast format.");
306 }
307 }
308
309
310
311 =head2 _set_score_stats
312
313 Usage : n/a; called automatically by _set_data()
314 Purpose : Sets various score statistics obtained from the HSP listing.
315 Argument : String with any of the following formats:
316 : blast2: Score = 30.1 bits (66), Expect = 9.2
317 : blast2: Score = 158.2 bits (544), Expect(2) = e-110
318 : blast1: Score = 410 (144.3 bits), Expect = 1.7e-40, P = 1.7e-40
319 : blast1: Score = 55 (19.4 bits), Expect = 5.3, Sum P(3) = 0.99
320 Throws : Exception if the stats cannot be parsed, probably due to a change
321 : in the Blast report format.
322
323 See Also : L<_set_data()|_set_data>
324
325 =cut
326
327 #--------------------
328 sub _set_score_stats {
329 #--------------------
330 my ($self, $data) = @_;
331
332 my ($expect, $p);
333
334 if($data =~ /Score = +([\d.e+-]+) bits \(([\d.e+-]+)\), +Expect = +([\d.e+-]+)/) {
335 # blast2 format n = 1
336 $self->{'_bits'} = $1;
337 $self->{'_score'} = $2;
338 $expect = $3;
339 } elsif($data =~ /Score = +([\d.e+-]+) bits \(([\d.e+-]+)\), +Expect\((\d+)\) = +([\d.e+-]+)/) {
340 # blast2 format n > 1
341 $self->{'_bits'} = $1;
342 $self->{'_score'} = $2;
343 $self->{'_n'} = $3;
344 $expect = $4;
345
346 } elsif($data =~ /Score = +([\d.e+-]+) \(([\d.e+-]+) bits\), +Expect = +([\d.e+-]+), P = +([\d.e-]+)/) {
347 # blast1 format, n = 1
348 $self->{'_score'} = $1;
349 $self->{'_bits'} = $2;
350 $expect = $3;
351 $p = $4;
352
353 } elsif($data =~ /Score = +([\d.e+-]+) \(([\d.e+-]+) bits\), +Expect = +([\d.e+-]+), +Sum P\((\d+)\) = +([\d.e-]+)/) {
354 # blast1 format, n > 1
355 $self->{'_score'} = $1;
356 $self->{'_bits'} = $2;
357 $expect = $3;
358 $self->{'_n'} = $4;
359 $p = $5;
360
361 } else {
362 $self->throw("Can't parse score statistics: unrecognized format.", "$data");
363 }
364
365 $expect = "1$expect" if $expect =~ /^e/i;
366 $p = "1$p" if defined $p and $p=~ /^e/i;
367
368 $self->{'_expect'} = $expect;
369 $self->{'_p'} = $p || undef;
370
371 }
372
373
374
375 =head2 _set_match_stats
376
377 Usage : n/a; called automatically by _set_data()
378 Purpose : Sets various matching statistics obtained from the HSP listing.
379 Argument : blast2: Identities = 23/74 (31%), Positives = 29/74 (39%), Gaps = 17/74 (22%)
380 : blast2: Identities = 57/98 (58%), Positives = 74/98 (75%)
381 : blast1: Identities = 87/204 (42%), Positives = 126/204 (61%)
382 : blast1: Identities = 87/204 (42%), Positives = 126/204 (61%), Frame = -3
383 : WU-blast: Identities = 310/553 (56%), Positives = 310/553 (56%), Strand = Minus / Plus
384 Throws : Exception if the stats cannot be parsed, probably due to a change
385 : in the Blast report format.
386 Comments : The "Gaps = " data in the HSP header has a different meaning depending
387 : on the type of Blast: for BLASTP, this number is the total number of
388 : gaps in query+sbjct; for TBLASTN, it is the number of gaps in the
389 : query sequence only. Thus, it is safer to collect the data
390 : separately by examining the actual sequence strings as is done
391 : in _set_seq().
392
393 See Also : L<_set_data()|_set_data>, L<_set_seq()|_set_seq>
394
395 =cut
396
397 #--------------------
398 sub _set_match_stats {
399 #--------------------
400 my ($self, $data) = @_;
401
402 if($data =~ m!Identities = (\d+)/(\d+)!) {
403 # blast1 or 2 format
404 $self->{'_numIdentical'} = $1;
405 $self->{'_totalLength'} = $2;
406 }
407
408 if($data =~ m!Positives = (\d+)/(\d+)!) {
409 # blast1 or 2 format
410 $self->{'_numConserved'} = $1;
411 $self->{'_totalLength'} = $2;
412 }
413
414 if($data =~ m!Frame = ([\d+-]+)!) {
415 $self->{'_frame'} = $1;
416 }
417
418 # Strand data is not always present in this line.
419 # _set_seq() will also set strand information.
420 if($data =~ m!Strand = (\w+) / (\w+)!) {
421 $self->{'_queryStrand'} = $1;
422 $self->{'_sbjctStrand'} = $2;
423 }
424
425 # if($data =~ m!Gaps = (\d+)/(\d+)!) {
426 # $self->{'_totalGaps'} = $1;
427 # } else {
428 # $self->{'_totalGaps'} = 0;
429 # }
430 }
431
432
433
434 =head2 _set_seq_data
435
436 Usage : n/a; called automatically when sequence data is requested.
437 Purpose : Sets the HSP sequence data for both query and sbjct sequences.
438 : Includes: start, stop, length, gaps, and raw sequence.
439 Argument : n/a
440 Throws : Propagates any exception thrown by _set_match_seq()
441 Comments : Uses raw data stored by _set_data() during object construction.
442 : These data are not always needed, so it is conditionally
443 : executed only upon demand by methods such as gaps(), _set_residues(),
444 : etc. _set_seq() does the dirty work.
445
446 See Also : L<_set_seq()|_set_seq>
447
448 =cut
449
450 sub _set_seq_data {
451 my $self = shift;
452
453 $self->_set_seq('query', @{$self->{'_queryList'}});
454 $self->_set_seq('sbjct', @{$self->{'_sbjctList'}});
455
456 # Liberate some memory.
457 @{$self->{'_queryList'}} = @{$self->{'_sbjctList'}} = ();
458 undef $self->{'_queryList'};
459 undef $self->{'_sbjctList'};
460
461 $self->{'_set_seq_data'} = 1;
462 }
463
464
465
466 =head2 _set_seq
467
468 Usage : n/a; called automatically by _set_seq_data()
469 : $hsp_obj->($seq_type, @data);
470 Purpose : Sets sequence information for both the query and sbjct sequences.
471 : Directly counts the number of gaps in each sequence (if gapped Blast).
472 Argument : $seq_type = 'query' or 'sbjct'
473 : @data = all seq lines with the form:
474 : Query: 61 SPHNVKDRKEQNGSINNAISPTATANTSGSQQINIDSALRDRSSNVAAQPSLSDASSGSN 120
475 Throws : Exception if data strings cannot be parsed, probably due to a change
476 : in the Blast report format.
477 Comments : Uses first argument to determine which data members to set
478 : making this method sensitive data member name changes.
479 : Behavior is dependent on the type of BLAST analysis (TBLASTN, BLASTP, etc).
480 Warning : Sequence endpoints are normalized so that start < end. This affects HSPs
481 : for TBLASTN/X hits on the minus strand. Normalization facilitates use
482 : of range information by methods such as match().
483
484 See Also : L<_set_seq_data()|_set_seq_data>, L<matches()|matches>, L<range()|range>, L<start()|start>, L<end()|end>
485
486 =cut
487
488 #-------------
489 sub _set_seq {
490 #-------------
491 my $self = shift;
492 my $seqType = shift;
493 my @data = @_;
494 my @ranges = ();
495 my @sequence = ();
496 my $numGaps = 0;
497
498 foreach( @data ) {
499 if( m/(\d+) *(\D+) *(\d+)/) {
500 push @ranges, ( $1, $3 ) ;
501 push @sequence, $2;
502 } else {
503 $self->warn("Bad sequence data: $_");
504 }
505 }
506
507 (scalar(@sequence) and scalar(@ranges)) || $self->throw("Can't set sequence: missing data",
508 "Possibly unrecognized Blast format.");
509
510 # Sensitive to member name changes.
511 $seqType = "_\L$seqType\E";
512 $self->{$seqType.'Start'} = $ranges[0];
513 $self->{$seqType.'Stop'} = $ranges[ $#ranges ];
514 $self->{$seqType.'Seq'} = \@sequence;
515
516 $self->{$seqType.'Length'} = abs($ranges[ $#ranges ] - $ranges[0]) + 1;
517
518 # Adjust lengths for BLASTX, TBLASTN, TBLASTX sequences
519 # Converting nucl coords to amino acid coords.
520
521 my $prog = $self->{'_prog'};
522 if($prog eq 'TBLASTN' and $seqType eq '_sbjct') {
523 $self->{$seqType.'Length'} /= 3;
524 } elsif($prog eq 'BLASTX' and $seqType eq '_query') {
525 $self->{$seqType.'Length'} /= 3;
526 } elsif($prog eq 'TBLASTX') {
527 $self->{$seqType.'Length'} /= 3;
528 }
529
530 $self->{$seqType.'Strand'} = 'Plus' if $prog =~ /BLAST[NX]/;
531
532 # Normalize sequence endpoints so that start < end.
533 # Reverse complement or 'minus strand' HSPs get flipped here.
534 if($self->{$seqType.'Start'} > $self->{$seqType.'Stop'}) {
535 ($self->{$seqType.'Start'}, $self->{$seqType.'Stop'}) =
536 ($self->{$seqType.'Stop'}, $self->{$seqType.'Start'});
537 $self->{$seqType.'Strand'} = 'Minus';
538 }
539
540 ## Count number of gaps in each seq. Only need to do this for gapped Blasts.
541 # if($self->{'_gapped'}) {
542 my $seqstr = join('', @sequence);
543 $seqstr =~ s/\s//g;
544 my $num_gaps = CORE::length($seqstr) - $self->{$seqType.'Length'};
545 $self->{$seqType.'Gaps'} = $num_gaps if $num_gaps > 0;
546 # }
547 }
548
549
550 =head2 _set_residues
551
552 Usage : n/a; called automatically when residue data is requested.
553 Purpose : Sets the residue numbers representing the identical and
554 : conserved positions. These data are obtained by analyzing the
555 : symbols between query and sbjct lines of the alignments.
556 Argument : n/a
557 Throws : Propagates any exception thrown by _set_seq_data() and _set_match_seq().
558 Comments : These data are not always needed, so it is conditionally
559 : executed only upon demand by methods such as seq_inds().
560 : Behavior is dependent on the type of BLAST analysis (TBLASTN, BLASTP, etc).
561
562 See Also : L<_set_seq_data()|_set_seq_data>, L<_set_match_seq()|_set_match_seq>, L<seq_inds()|seq_inds>
563
564 =cut
565
566 #------------------
567 sub _set_residues {
568 #------------------
569 my $self = shift;
570 my @sequence = ();
571
572 $self->_set_seq_data() unless $self->{'_set_seq_data'};
573
574 # Using hashes to avoid saving duplicate residue numbers.
575 my %identicalList_query = ();
576 my %identicalList_sbjct = ();
577 my %conservedList_query = ();
578 my %conservedList_sbjct = ();
579
580 my $aref = $self->_set_match_seq() if not ref $self->{'_matchSeq'};
581 $aref ||= $self->{'_matchSeq'};
582 my $seqString = join('', @$aref );
583
584 my $qseq = join('',@{$self->{'_querySeq'}});
585 my $sseq = join('',@{$self->{'_sbjctSeq'}});
586 my $resCount_query = $self->{'_queryStop'} || 0;
587 my $resCount_sbjct = $self->{'_sbjctStop'} || 0;
588
589 my $prog = $self->{'_prog'};
590 if($prog !~ /^BLASTP|^BLASTN/) {
591 if($prog eq 'TBLASTN') {
592 $resCount_sbjct /= 3;
593 } elsif($prog eq 'BLASTX') {
594 $resCount_query /= 3;
595 } elsif($prog eq 'TBLASTX') {
596 $resCount_query /= 3;
597 $resCount_sbjct /= 3;
598 }
599 }
600
601 my ($mchar, $schar, $qchar);
602 while( $mchar = chop($seqString) ) {
603 ($qchar, $schar) = (chop($qseq), chop($sseq));
604 if( $mchar eq '+' ) {
605 $conservedList_query{ $resCount_query } = 1;
606 $conservedList_sbjct{ $resCount_sbjct } = 1;
607 } elsif( $mchar ne ' ' ) {
608 $identicalList_query{ $resCount_query } = 1;
609 $identicalList_sbjct{ $resCount_sbjct } = 1;
610 }
611 $resCount_query-- if $qchar ne $GAP_SYMBOL;
612 $resCount_sbjct-- if $schar ne $GAP_SYMBOL;
613 }
614 $self->{'_identicalRes_query'} = \%identicalList_query;
615 $self->{'_conservedRes_query'} = \%conservedList_query;
616 $self->{'_identicalRes_sbjct'} = \%identicalList_sbjct;
617 $self->{'_conservedRes_sbjct'} = \%conservedList_sbjct;
618
619 }
620
621
622
623
624 =head2 _set_match_seq
625
626 Usage : n/a. Internal method.
627 : $hsp_obj->_set_match_seq()
628 Purpose : Set the 'match' sequence for the current HSP (symbols in between
629 : the query and sbjct lines.)
630 Returns : Array reference holding the match sequences lines.
631 Argument : n/a
632 Throws : Exception if the _matchList field is not set.
633 Comments : The match information is not always necessary. This method
634 : allows it to be conditionally prepared.
635 : Called by _set_residues>() and seq_str().
636
637 See Also : L<_set_residues()|_set_residues>, L<seq_str()|seq_str>
638
639 =cut
640
641 #-------------------
642 sub _set_match_seq {
643 #-------------------
644 my $self = shift;
645
646 ## DEBUGGING CODE:
647 # if($self->parent->name eq '1AK5_' and $self->parent->parent->name eq 'YAR073W') {
648 # print "\n_set_match_seq() called for HSP ", $self->name, " of hit ${\$self->parent->name} in query ${\$self->parent->parent->name}"; <STDIN>;
649 # }
650
651 ref($self->{'_matchList'}) || $self->throw("Can't set HSP match sequence: No data");
652
653 my @data = @{$self->{'_matchList'}};
654
655 my(@sequence);
656 foreach( @data ) {
657 chomp($_);
658 ## Remove leading spaces; (note: aln may begin with a space
659 ## which is why we can't use s/^ +//).
660 s/^ {$self->{'_match_indent'}}//;
661 push @sequence, $_;
662 }
663 # Liberate some memory.
664 @{$self->{'_matchList'}} = undef;
665 $self->{'_matchList'} = undef;
666
667 $self->{'_matchSeq'} = \@sequence;
668
669 ## DEBUGGING CODE:
670 # if($self->parent->name eq '1AK5_' and $self->parent->parent->name eq 'YAR073W') {
671 # print "RETURNING: $self->{'_matchSeq'}:\n @{$self->{'_matchSeq'}}";<STDIN>;
672 # }
673
674 $self->{'_matchSeq'};
675 }
676
677
678
679 =head2 score
680
681 Usage : $hsp_obj->score()
682 Purpose : Get the Blast score for the HSP.
683 Returns : Integer
684 Argument : n/a
685 Throws : n/a
686
687 See Also : L<bits()|bits>
688
689 =cut
690
691 #---------
692 sub score { my $self = shift; $self->{'_score'}; }
693 #---------
694
695
696
697 =head2 bits
698
699 Usage : $hsp_obj->bits()
700 Purpose : Get the Blast score in bits for the HSP.
701 Returns : Float
702 Argument : n/a
703 Throws : n/a
704
705
706 See Also : L<score()|score>
707
708 =cut
709
710 #--------
711 sub bits { my $self = shift; $self->{'_bits'}; }
712 #--------
713
714
715
716 =head2 n
717
718 Usage : $hsp_obj->n()
719 Purpose : Get the N value (num HSPs on which P/Expect is based).
720 : This value is not defined with NCBI Blast2 with gapping.
721 Returns : Integer or null string if not defined.
722 Argument : n/a
723 Throws : n/a
724 Comments : The 'N' value is listed in parenthesis with P/Expect value:
725 : e.g., P(3) = 1.2e-30 ---> (N = 3).
726 : Not defined in NCBI Blast2 with gaps.
727 : This typically is equal to the number of HSPs but not always.
728 : To obtain the number of HSPs, use Bio::Tools::Blast::Sbjct::num_hsps().
729
730 See Also : L<score()|score>
731
732 =cut
733
734 #-----
735 sub n { my $self = shift; $self->{'_n'} || ''; }
736 #-----
737
738
739
740 =head2 frame
741
742 Usage : $hsp_obj->frame()
743 Purpose : Get the reading frame number (-/+ 1, 2, 3) (TBLASTN/X only).
744 Returns : Integer or null string if not defined.
745 Argument : n/a
746 Throws : n/a
747
748 =cut
749
750 #---------
751 sub frame { my $self = shift; $self->{'_frame'} || ''; }
752 #---------
753
754
755
756 =head2 signif()
757
758 Usage : $hsp_obj->signif()
759 Purpose : Get the P-value or Expect value for the HSP.
760 Returns : Float (0.001 or 1.3e-43)
761 : Returns P-value if it is defined, otherwise, Expect value.
762 Argument : n/a
763 Throws : n/a
764 Comments : Provided for consistency with Sbjct::signif()
765 : Support for returning the significance data in different
766 : formats (e.g., exponent only), is not provided for HSP objects.
767 : This is only available for the Sbjct or Blast object.
768
769 See Also : L<p()|p>, L<expect()|expect>, B<Bio::Tools::Blast::Sbjct::signif()>
770
771 =cut
772
773 #-----------
774 sub signif {
775 #-----------
776 my $self = shift;
777 my $val ||= defined($self->{'_p'}) ? $self->{'_p'} : $self->{'_expect'};
778 $val;
779 }
780
781
782
783 =head2 expect
784
785 Usage : $hsp_obj->expect()
786 Purpose : Get the Expect value for the HSP.
787 Returns : Float (0.001 or 1.3e-43)
788 Argument : n/a
789 Throws : n/a
790 Comments : Support for returning the expectation data in different
791 : formats (e.g., exponent only), is not provided for HSP objects.
792 : This is only available for the Sbjct or Blast object.
793
794 See Also : L<p()|p>
795
796 =cut
797
798 #----------
799 sub expect { my $self = shift; $self->{'_expect'}; }
800 #----------
801
802
803
804 =head2 p
805
806 Usage : $hsp_obj->p()
807 Purpose : Get the P-value for the HSP.
808 Returns : Float (0.001 or 1.3e-43) or undef if not defined.
809 Argument : n/a
810 Throws : n/a
811 Comments : P-value is not defined with NCBI Blast2 reports.
812 : Support for returning the expectation data in different
813 : formats (e.g., exponent only) is not provided for HSP objects.
814 : This is only available for the Sbjct or Blast object.
815
816 See Also : L<expect()|expect>
817
818 =cut
819
820 #-----
821 sub p { my $self = shift; $self->{'_p'}; }
822 #-----
823
824
825 =head2 length
826
827 Usage : $hsp->length( [seq_type] )
828 Purpose : Get the length of the aligned portion of the query or sbjct.
829 Example : $hsp->length('query')
830 Returns : integer
831 Argument : seq_type: 'query' | 'sbjct' | 'total' (default = 'total')
832 Throws : n/a
833 Comments : 'total' length is the full length of the alignment
834 : as reported in the denominators in the alignment section:
835 : "Identical = 34/120 Positives = 67/120".
836 : Developer note: when using the built-in length function within
837 : this module, call it as CORE::length().
838
839 See Also : L<gaps()|gaps>
840
841 =cut
842
843 #-----------
844 sub length {
845 #-----------
846 my( $self, $type ) = @_;
847 $type ||= 'total';
848
849 $type ne 'total' and $self->_set_seq_data() unless $self->{'_set_seq_data'};
850
851 ## Sensitive to member name format.
852 $type = "_\L$type\E";
853 $self->{$type.'Length'};
854 }
855
856
857
858 =head2 gaps
859
860 Usage : $hsp->gaps( [seq_type] )
861 Purpose : Get the number of gaps in the query, sbjct, or total alignment.
862 : Also can return query gaps and sbjct gaps as a two-element list
863 : when in array context.
864 Example : $total_gaps = $hsp->gaps();
865 : ($qgaps, $sgaps) = $hsp->gaps();
866 : $qgaps = $hsp->gaps('query');
867 Returns : scalar context: integer
868 : array context without args: (int, int) = ('queryGaps', 'sbjctGaps')
869 Argument : seq_type: 'query' | 'sbjct' | 'total'
870 : (default = 'total', scalar context)
871 : Array context can be "induced" by providing an argument of 'list' or 'array'.
872 Throws : n/a
873
874 See Also : L<length()|length>, L<matches()|matches>
875
876 =cut
877
878 #---------
879 sub gaps {
880 #---------
881 my( $self, $seqType ) = @_;
882
883 $self->_set_seq_data() unless $self->{'_set_seq_data'};
884
885 $seqType ||= (wantarray ? 'list' : 'total');
886
887 if($seqType =~ /list|array/i) {
888 return (($self->{'_queryGaps'} || 0), ($self->{'_sbjctGaps'} || 0));
889 }
890
891 if($seqType eq 'total') {
892 return ($self->{'_queryGaps'} + $self->{'_sbjctGaps'}) || 0;
893 } else {
894 ## Sensitive to member name format.
895 $seqType = "_\L$seqType\E";
896 return $self->{$seqType.'Gaps'} || 0;
897 }
898 }
899
900
901
902 =head2 matches
903
904 Usage : $hsp->matches([seq_type], [start], [stop]);
905 Purpose : Get the total number of identical and conservative matches
906 : in the query or sbjct sequence for the given HSP. Optionally can
907 : report data within a defined interval along the seq.
908 : (Note: 'conservative' matches are called 'positives' in the
909 : Blast report.)
910 Example : ($id,$cons) = $hsp_object->matches('sbjct');
911 : ($id,$cons) = $hsp_object->matches('query',300,400);
912 Returns : 2-element array of integers
913 Argument : (1) seq_type = 'query' | 'sbjct' (default = query)
914 : (2) start = Starting coordinate (optional)
915 : (3) stop = Ending coordinate (optional)
916 Throws : Exception if the supplied coordinates are out of range.
917 Comments : Relies on seq_str('match') to get the string of alignment symbols
918 : between the query and sbjct lines which are used for determining
919 : the number of identical and conservative matches.
920
921 See Also : L<length()|length>, L<gaps()|gaps>, L<seq_str()|seq_str>, B<Bio::Tools::Blast::Sbjct::_adjust_contigs()>
922
923 =cut
924
925 #-----------
926 sub matches {
927 #-----------
928 my( $self, %param ) = @_;
929 my(@data);
930 my($seqType, $beg, $end) = ($param{-SEQ}, $param{-START}, $param{-STOP});
931 $seqType ||= 'query';
932
933 if(!defined $beg && !defined $end) {
934 ## Get data for the whole alignment.
935 push @data, ($self->{'_numIdentical'}, $self->{'_numConserved'});
936 } else {
937 ## Get the substring representing the desired sub-section of aln.
938 $beg ||= 0;
939 $end ||= 0;
940 my($start,$stop) = $self->range($seqType);
941 if($beg == 0) { $beg = $start; $end = $beg+$end; }
942 elsif($end == 0) { $end = $stop; $beg = $end-$beg; }
943
944 if($end >= $stop) { $end = $stop; } ##ML changed from if (end >stop)
945 else { $end += 1;} ##ML moved from commented position below, makes
946 ##more sense here
947 # if($end > $stop) { $end = $stop; }
948 if($beg < $start) { $beg = $start; }
949 # else { $end += 1;}
950
951 # my $seq = substr($self->seq_str('match'), $beg-$start, ($end-$beg));
952
953 ## ML: START fix for substr out of range error ------------------
954 my $seq = "";
955 if (($self->{'_prog'} eq 'TBLASTN') and ($seqType eq 'sbjct'))
956 {
957 $seq = substr($self->seq_str('match'),
958 int(($beg-$start)/3), int(($end-$beg+1)/3));
959
960 } elsif (($self->{'_prog'} eq 'BLASTX') and ($seqType eq 'query'))
961 {
962 $seq = substr($self->seq_str('match'),
963 int(($beg-$start)/3), int(($end-$beg+1)/3));
964 } else {
965 $seq = substr($self->seq_str('match'),
966 $beg-$start, ($end-$beg));
967 }
968 ## ML: End of fix for substr out of range error -----------------
969
970
971 ## ML: debugging code
972 ## This is where we get our exception. Try printing out the values going
973 ## into this:
974 ##
975 # print STDERR
976 # qq(*------------MY EXCEPTION --------------------\nSeq: ") ,
977 # $self->seq_str("$seqType"), qq("\n),$self->name,",( index:";
978 # print STDERR $beg-$start, ", len: ", $end-$beg," ), (HSPRealLen:",
979 # CORE::length $self->seq_str("$seqType");
980 # print STDERR ", HSPCalcLen: ", $stop - $start +1 ," ),
981 # ( beg: $beg, end: $end ), ( start: $start, stop: stop )\n";
982 ## ML: END DEBUGGING CODE----------
983
984 if(!CORE::length $seq) {
985 $self->throw("Undefined sub-sequence ($beg,$end). Valid range = $start - $stop");
986 }
987 ## Get data for a substring.
988 # printf "Collecting HSP subsection data: beg,end = %d,%d; start,stop = %d,%d\n%s<---\n", $beg, $end, $start, $stop, $seq;
989 # printf "Original match seq:\n%s\n",$self->seq_str('match');
990 $seq =~ s/ //g; # remove space (no info).
991 my $len_cons = CORE::length $seq;
992 $seq =~ s/\+//g; # remove '+' characters (conservative substitutions)
993 my $len_id = CORE::length $seq;
994 push @data, ($len_id, $len_cons);
995 # printf " HSP = %s\n id = %d; cons = %d\n", $self->name, $len_id, $len_cons; <STDIN>;
996 }
997 @data;
998 }
999
1000
1001
1002 =head2 frac_identical
1003
1004 Usage : $hsp_object->frac_identical( [seq_type] );
1005 Purpose : Get the fraction of identical positions within the given HSP.
1006 Example : $frac_iden = $hsp_object->frac_identical('query');
1007 Returns : Float (2-decimal precision, e.g., 0.75).
1008 Argument : seq_type: 'query' | 'sbjct' | 'total'
1009 : default = 'total' (but see comments below).
1010 Throws : n/a
1011 Comments : Different versions of Blast report different values for the total
1012 : length of the alignment. This is the number reported in the
1013 : denominators in the stats section:
1014 : "Identical = 34/120 Positives = 67/120".
1015 : BLAST-GP uses the total length of the alignment (with gaps)
1016 : WU-BLAST uses the length of the query sequence (without gaps).
1017 : Therefore, when called without an argument or an argument of 'total',
1018 : this method will report different values depending on the
1019 : version of BLAST used.
1020 :
1021 : To get the fraction identical among only the aligned residues,
1022 : ignoring the gaps, call this method with an argument of 'query'
1023 : or 'sbjct'.
1024
1025 See Also : L<frac_conserved()|frac_conserved>, L<num_identical()|num_identical>, L<matches()|matches>
1026
1027 =cut
1028
1029 #-------------------
1030 sub frac_identical {
1031 #-------------------
1032 # The value is calculated as opposed to storing it from the parsed results.
1033 # This saves storage and also permits flexibility in determining for which
1034 # sequence (query or sbjct) the figure is to be calculated.
1035
1036 my( $self, $seqType ) = @_;
1037 $seqType ||= 'total';
1038
1039 if($seqType ne 'total') {
1040 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1041 }
1042 ## Sensitive to member name format.
1043 $seqType = "_\L$seqType\E";
1044
1045 sprintf( "%.2f", $self->{'_numIdentical'}/$self->{$seqType.'Length'});
1046 }
1047
1048
1049 =head2 frac_conserved
1050
1051 Usage : $hsp_object->frac_conserved( [seq_type] );
1052 Purpose : Get the fraction of conserved positions within the given HSP.
1053 : (Note: 'conservative' positions are called 'positives' in the
1054 : Blast report.)
1055 Example : $frac_cons = $hsp_object->frac_conserved('query');
1056 Returns : Float (2-decimal precision, e.g., 0.75).
1057 Argument : seq_type: 'query' | 'sbjct'
1058 : default = 'total' (but see comments below).
1059 Throws : n/a
1060 Comments : Different versions of Blast report different values for the total
1061 : length of the alignment. This is the number reported in the
1062 : denominators in the stats section:
1063 : "Identical = 34/120 Positives = 67/120".
1064 : BLAST-GP uses the total length of the alignment (with gaps)
1065 : WU-BLAST uses the length of the query sequence (without gaps).
1066 : Therefore, when called without an argument or an argument of 'total',
1067 : this method will report different values depending on the
1068 : version of BLAST used.
1069 :
1070 : To get the fraction conserved among only the aligned residues,
1071 : ignoring the gaps, call this method with an argument of 'query'
1072 : or 'sbjct'.
1073
1074 See Also : L<frac_conserved()|frac_conserved>, L<num_conserved()|num_conserved>, L<matches()|matches>
1075
1076 =cut
1077
1078 #--------------------
1079 sub frac_conserved {
1080 #--------------------
1081 # The value is calculated as opposed to storing it from the parsed results.
1082 # This saves storage and also permits flexibility in determining for which
1083 # sequence (query or sbjct) the figure is to be calculated.
1084
1085 my( $self, $seqType ) = @_;
1086 $seqType ||= 'total';
1087
1088 if($seqType ne 'total') {
1089 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1090 }
1091
1092 ## Sensitive to member name format.
1093 $seqType = "_\L$seqType\E";
1094
1095 sprintf( "%.2f", $self->{'_numConserved'}/$self->{$seqType.'Length'});
1096 }
1097
1098
1099 =head2 num_identical
1100
1101 Usage : $hsp_object->num_identical();
1102 Purpose : Get the number of identical positions within the given HSP.
1103 Example : $num_iden = $hsp_object->num_identical();
1104 Returns : integer
1105 Argument : n/a
1106 Throws : n/a
1107
1108 See Also : L<num_conserved()|num_conserved>, L<frac_identical()|frac_identical>
1109
1110 =cut
1111
1112 #-------------------
1113 sub num_identical {
1114 #-------------------
1115 my( $self) = shift;
1116
1117 $self->{'_numIdentical'};
1118 }
1119
1120
1121 =head2 num_conserved
1122
1123 Usage : $hsp_object->num_conserved();
1124 Purpose : Get the number of conserved positions within the given HSP.
1125 Example : $num_iden = $hsp_object->num_conserved();
1126 Returns : integer
1127 Argument : n/a
1128 Throws : n/a
1129
1130 See Also : L<num_identical()|num_identical>, L<frac_conserved()|frac_conserved>
1131
1132 =cut
1133
1134 #-------------------
1135 sub num_conserved {
1136 #-------------------
1137 my( $self) = shift;
1138
1139 $self->{'_numConserved'};
1140 }
1141
1142
1143
1144 =head2 range
1145
1146 Usage : $hsp->range( [seq_type] );
1147 Purpose : Gets the (start, end) coordinates for the query or sbjct sequence
1148 : in the HSP alignment.
1149 Example : ($qbeg, $qend) = $hsp->range('query');
1150 : ($sbeg, $send) = $hsp->range('sbjct');
1151 Returns : Two-element array of integers
1152 Argument : seq_type = string, 'query' or 'sbjct' (default = 'query')
1153 : (case insensitive).
1154 Throws : n/a
1155
1156 See Also : L<start()|start>, L<end()|end>
1157
1158 =cut
1159
1160 #----------
1161 sub range {
1162 #----------
1163 my ($self, $seqType) = @_;
1164
1165 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1166
1167 $seqType ||= 'query';
1168 ## Sensitive to member name changes.
1169 $seqType = "_\L$seqType\E";
1170
1171 return ($self->{$seqType.'Start'},$self->{$seqType.'Stop'});
1172 }
1173
1174 =head2 start
1175
1176 Usage : $hsp->start( [seq_type] );
1177 Purpose : Gets the start coordinate for the query, sbjct, or both sequences
1178 : in the HSP alignment.
1179 Example : $qbeg = $hsp->start('query');
1180 : $sbeg = $hsp->start('sbjct');
1181 : ($qbeg, $sbeg) = $hsp->start();
1182 Returns : scalar context: integer
1183 : array context without args: list of two integers
1184 Argument : In scalar context: seq_type = 'query' or 'sbjct'
1185 : (case insensitive). If not supplied, 'query' is used.
1186 : Array context can be "induced" by providing an argument of 'list' or 'array'.
1187 Throws : n/a
1188
1189 See Also : L<end()|end>, L<range()|range>
1190
1191 =cut
1192
1193 #----------
1194 sub start {
1195 #----------
1196 my ($self, $seqType) = @_;
1197
1198 $seqType ||= (wantarray ? 'list' : 'query');
1199
1200 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1201
1202 if($seqType =~ /list|array/i) {
1203 return ($self->{'_queryStart'}, $self->{'_sbjctStart'});
1204 } else {
1205 ## Sensitive to member name changes.
1206 $seqType = "_\L$seqType\E";
1207 return $self->{$seqType.'Start'};
1208 }
1209 }
1210
1211 =head2 end
1212
1213 Usage : $hsp->end( [seq_type] );
1214 Purpose : Gets the end coordinate for the query, sbjct, or both sequences
1215 : in the HSP alignment.
1216 Example : $qbeg = $hsp->end('query');
1217 : $sbeg = $hsp->end('sbjct');
1218 : ($qbeg, $sbeg) = $hsp->end();
1219 Returns : scalar context: integer
1220 : array context without args: list of two integers
1221 Argument : In scalar context: seq_type = 'query' or 'sbjct'
1222 : (case insensitive). If not supplied, 'query' is used.
1223 : Array context can be "induced" by providing an argument of 'list' or 'array'.
1224 Throws : n/a
1225
1226 See Also : L<start()|start>, L<range()|range>
1227
1228 =cut
1229
1230 #----------
1231 sub end {
1232 #----------
1233 my ($self, $seqType) = @_;
1234
1235 $seqType ||= (wantarray ? 'list' : 'query');
1236
1237 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1238
1239 if($seqType =~ /list|array/i) {
1240 return ($self->{'_queryStop'}, $self->{'_sbjctStop'});
1241 } else {
1242 ## Sensitive to member name changes.
1243 $seqType = "_\L$seqType\E";
1244 return $self->{$seqType.'Stop'};
1245 }
1246 }
1247
1248
1249
1250 =head2 strand
1251
1252 Usage : $hsp_object->strand( [seq_type] )
1253 Purpose : Get the strand of the query or sbjct sequence.
1254 Example : print $hsp->strand('query');
1255 : ($qstrand, $sstrand) = $hsp->strand();
1256 Returns : -1, 0, or 1
1257 : -1 = Minus strand, +1 = Plus strand
1258 : Returns 0 if strand is not defined, which occurs
1259 : for non-TBLASTN/X reports.
1260 : In scalar context without arguments, returns queryStrand value.
1261 : In array context without arguments, returns a two-element list
1262 : of strings (queryStrand, sbjctStrand).
1263 : Array context can be "induced" by providing an argument of 'list' or 'array'.
1264 Argument : seq_type: 'query' | 'sbjct' or undef
1265 Throws : n/a
1266
1267 See Also : L<_set_seq()|_set_seq>, L<_set_match_stats()|_set_match_stats>
1268
1269 =cut
1270
1271 #-----------
1272 sub strand {
1273 #-----------
1274 my( $self, $seqType ) = @_;
1275 $seqType ||= (wantarray ? 'list' : 'query');
1276
1277 return '' if $seqType eq 'query' and $self->{'_prog'} eq 'TBLASTN';
1278
1279 ## Sensitive to member name format.
1280 $seqType = "_\L$seqType\E";
1281
1282 # $seqType could be '_list'.
1283 $self->{'_queryStrand'} or $self->_set_seq_data() unless $self->{'_set_seq_data'};
1284
1285 if($seqType =~ /list|array/i) {
1286 return ('','') unless defined $self->{'_queryStrand'};
1287 return ($self->{'_queryStrand'}, $self->{'_sbjctStrand'});
1288 }
1289 local $^W = 0;
1290 $STRAND_SYMBOL{$self->{$seqType.'Strand'}} || 0;
1291 }
1292
1293
1294 #####################################################################################
1295 ## INSTANCE METHODS ##
1296 #####################################################################################
1297
1298
1299 =head2 seq
1300
1301 Usage : $hsp->seq( [seq_type] );
1302 Purpose : Get the query or sbjct sequence as a Bio::Seq.pm object.
1303 Example : $seqObj = $hsp->seq('query');
1304 Returns : Object reference for a Bio::Seq.pm object.
1305 Argument : seq_type = 'query' or 'sbjct' (default = 'query').
1306 Throws : Propagates any exception that occurs during construction
1307 : of the Bio::Seq.pm object.
1308 Comments : The sequence is returned in an array of strings corresponding
1309 : to the strings in the original format of the Blast alignment.
1310 : (i.e., same spacing).
1311
1312 See Also : L<seq_str()|seq_str>, L<seq_inds()|seq_inds>, B<Bio::Seq.pm>
1313
1314 =cut
1315
1316 #-------
1317 sub seq {
1318 #-------
1319 my($self,$seqType) = @_;
1320 $seqType ||= 'query';
1321 my $str = $self->seq_str($seqType);
1322 my $num = $self->name;
1323 my $name = $seqType =~ /query/i
1324 ? $self->parent->parent->name
1325 : $self->parent->name;
1326
1327 require Bio::Seq;
1328
1329 new Bio::Seq (-ID => $name,
1330 -SEQ => $str,
1331 -DESC => "Blast HSP #$num, $seqType sequence",
1332 );
1333 }
1334
1335
1336
1337 =head2 seq_str
1338
1339 Usage : $hsp->seq_str( seq_type );
1340 Purpose : Get the full query, sbjct, or 'match' sequence as a string.
1341 : The 'match' sequence is the string of symbols in between the
1342 : query and sbjct sequences.
1343 Example : $str = $hsp->seq_str('query');
1344 Returns : String
1345 Argument : seq_Type = 'query' or 'sbjct' or 'match'
1346 Throws : Exception if the argument does not match an accepted seq_type.
1347 Comments : Calls _set_residues() to set the 'match' sequence if it has
1348 : not been set already.
1349
1350 See Also : L<seq()|seq>, L<seq_inds()|seq_inds>, L<_set_match_seq()|_set_match_seq>
1351
1352 =cut
1353
1354 #------------
1355 sub seq_str {
1356 #------------
1357 my($self,$seqType) = @_;
1358
1359 ## Sensitive to member name changes.
1360 $seqType = "_\L$seqType\E";
1361
1362 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1363
1364 if($seqType =~ /sbjct|query/) {
1365 my $seq = join('',@{$self->{$seqType.'Seq'}});
1366 $seq =~ s/\s+//g;
1367 return $seq;
1368
1369 } elsif( $seqType =~ /match/i) {
1370 # Only need to call _set_match_seq() if the match seq is requested.
1371 my $aref = $self->_set_match_seq() unless ref $self->{'_matchSeq'};
1372 $aref = $self->{'_matchSeq'};
1373
1374 ## DEBUGGING CODE:
1375 # if($self->parent->name eq '1AK5_' and $self->parent->parent->name eq 'YAR073W') {
1376 # print "seq_str():\n @$aref";<STDIN>;
1377 # }
1378
1379 return join('',@$aref);
1380
1381 } else {
1382 $self->throw("Invalid or undefined sequence type: $seqType",
1383 "Valid types: query, sbjct, match");
1384 }
1385 }
1386
1387
1388
1389
1390 =head2 seq_inds
1391
1392 Usage : $hsp->seq_inds( seq_type, class, collapse );
1393 Purpose : Get a list of residue positions (indices) for all identical
1394 : or conserved residues in the query or sbjct sequence.
1395 Example : @ind = $hsp->seq_inds('query', 'identical');
1396 : @ind = $hsp->seq_inds('sbjct', 'conserved');
1397 : @ind = $hsp->seq_inds('sbjct', 'conserved', 1);
1398 Returns : List of integers
1399 : May include ranges if collapse is true.
1400 Argument : seq_type = 'query' or 'sbjct' (default = query)
1401 : class = 'identical' or 'conserved' (default = identical)
1402 : (can be shortened to 'id' or 'cons')
1403 : (actually, anything not 'id' will evaluate to 'conserved').
1404 : collapse = boolean, if true, consecutive positions are merged
1405 : using a range notation, e.g., "1 2 3 4 5 7 9 10 11"
1406 : collapses to "1-5 7 9-11". This is useful for
1407 : consolidating long lists. Default = no collapse.
1408 Throws : n/a.
1409 Comments : Calls _set_residues() to set the 'match' sequence if it has
1410 : not been set already.
1411
1412 See Also : L<seq()|seq>, L<_set_residues()|_set_residues>, L<collapse_nums()|collapse_nums>, B<Bio::Tools::Blast::Sbjct::seq_inds()>
1413
1414 =cut
1415
1416 #---------------
1417 sub seq_inds {
1418 #---------------
1419 my ($self, $seq, $class, $collapse) = @_;
1420
1421 $seq ||= 'query';
1422 $class ||= 'identical';
1423 $collapse ||= 0;
1424
1425 $self->_set_residues() unless defined $self->{'_identicalRes_query'};
1426
1427 $seq = ($seq !~ /^q/i ? 'sbjct' : 'query');
1428 $class = ($class !~ /^id/i ? 'conserved' : 'identical');
1429
1430 ## Sensitive to member name changes.
1431 $seq = "_\L$seq\E";
1432 $class = "_\L$class\E";
1433
1434 my @ary = sort { $a <=> $b } keys %{ $self->{"${class}Res$seq"}};
1435
1436 return $collapse ? &collapse_nums(@ary) : @ary;
1437 }
1438
1439
1440
1441
1442 =head2 get_aln
1443
1444 Usage : $hsp->get_aln()
1445 Purpose : Get a Bio::UnivAln.pm object constructed from the query + sbjct
1446 : sequences of the present HSP object.
1447 Example : $aln_obj = $hsp->get_aln();
1448 Returns : Object reference for a Bio::UnivAln.pm object.
1449 Argument : n/a.
1450 Throws : Propagates any exception ocurring during the construction of
1451 : the Bio::UnivAln object.
1452 Comments : Requires Bio::UnivAln.pm.
1453 : The Bio::UnivAln.pm object is constructed from the query + sbjct
1454 : sequence objects obtained by calling seq().
1455 : Gap residues are included (see $GAP_SYMBOL). It is important that
1456 : Bio::UnivAln.pm recognizes the gaps correctly. A strategy for doing
1457 : this is being considered. Currently it is hard-wired.
1458
1459 See Also : L<seq()|seq>, B<Bio::UnivAln.pm>
1460
1461 =cut
1462
1463 #------------
1464 sub get_aln {
1465 #------------
1466 my $self = shift;
1467
1468 require Bio::UnivAln;
1469
1470 my $qseq = $self->seq('query');
1471 my $sseq = $self->seq('sbjct');
1472
1473 my $desc = sprintf "HSP #%s of query %s vs. sbjct %s",
1474 $self->name, $self->parent->parent->name, $self->parent->name;
1475
1476 my $type = $self->{'_prog'} =~ /P$|^T/ ? 'amino' : 'dna';
1477
1478 Bio::UnivAln->new( -seqs => [$qseq, $sseq],
1479 -desc => $desc,
1480 -type => $type,
1481 );
1482 }
1483
1484
1485 =head2 display
1486
1487 Usage : $sbjct_object->display( %named_parameters );
1488 Purpose : Display information about Bio::Tools::Blast::Sbjct.pm data members
1489 : including: length, gaps, score, significance value,
1490 : sequences and sequence indices.
1491 Example : $object->display(-SHOW=>'stats');
1492 Argument : Named parameters: (TAGS CAN BE UPPER OR LOWER CASE)
1493 : -SHOW => 'hsp',
1494 : -WHERE => filehandle (default = STDOUT)
1495 Returns : n/a
1496 Status : Experimental
1497 Comments : For more control over the display of sequence data,
1498 : use seq(), seq_str(), seq_inds().
1499
1500 See Also : L<_display_seq()|_display_seq>, L<seq()|seq>, L<seq_str()|seq_str>, L<seq_inds()|seq_inds>, L<_display_matches()|_display_matches>, B<Bio::Root::Object::display()>
1501
1502 =cut
1503
1504 #-----------
1505 sub display {
1506 #-----------
1507 my( $self, %param ) = @_;
1508
1509 my $sbjctName = $self->parent->name();
1510 my $queryName = $self->parent->parent->name();
1511 my $layout = $self->parent->parent->_layout();
1512
1513 my $OUT = $self->set_display(%param);
1514
1515 printf( $OUT "%-15s: %d\n", "LENGTH TOTAL", $self->length('total') );
1516 printf( $OUT "%-15s: %d\n", "LENGTH QUERY", $self->length('query') );
1517 printf( $OUT "%-15s: %d\n", "LENGTH SBJCT", $self->length('sbjct') );
1518 printf( $OUT "%-15s: %d\n", "GAPS QUERY", $self->gaps('query') );
1519 printf( $OUT "%-15s: %d\n", "GAPS SBJCT", $self->gaps('sbjct') );
1520 printf( $OUT "%-15s: %d\n", "SCORE", $self->{'_score'} );
1521 printf( $OUT "%-15s: %0.1f\n", "BITS", $self->{'_bits'} );
1522 if($layout == 1) {
1523 printf( $OUT "%-15s: %.1e\n", "P-VAL", $self->{'_p'} );
1524 printf( $OUT "%-15s: %.1e\n", "EXPECT", $self->{'_expect'} );
1525 } else {
1526 printf( $OUT "%-15s: %.1e\n", "EXPECT", $self->{'_expect'} );
1527 }
1528
1529 my $queryLength = $self->length('query');
1530
1531 printf( $OUT "%-15s: %d (%0.0f%%)\n", "IDENTICAL", $self->{'_numIdentical'},
1532 $self->{'_numIdentical'}/$queryLength * 100 );
1533 printf( $OUT "%-15s: %d (%0.0f%%) %s \n", "CONSERVED", $self->{'_numConserved'},
1534 $self->{'_numConserved'}/$queryLength * 100,
1535 "includes identical" );
1536
1537 $self->_display_seq('query', $queryName, $OUT);
1538 $self->_display_seq('sbjct', $sbjctName, $OUT);
1539 $self->_display_matches($queryName, $sbjctName, $OUT);
1540 }
1541
1542
1543
1544
1545 =head2 _display_seq
1546
1547 Usage : n/a; called automatically by display()
1548 Purpose : Display information about query and sbjct HSP sequences.
1549 : Prints the start, stop coordinates and the actual sequence.
1550 Example : n/a
1551 Argument :
1552 Returns : printf call.
1553 Status : Experimental
1554 Comments : For more control, use seq(), seq_str(), or seq_inds().
1555
1556 See Also : L<display()|display>, L<seq()|seq>, L<seq_str()|seq_str>, L<seq_inds()|seq_inds>, L<_display_matches()|_display_matches>
1557
1558 =cut
1559
1560 #------------------
1561 sub _display_seq {
1562 #------------------
1563 my( $self, $seqType, $name, $OUT ) = @_;
1564
1565 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1566
1567 # Sensitive to member name changes.
1568 my $mem = "_\L$seqType\E";
1569 printf( $OUT "\n%10s: %s\n%10s %s\n", "\U$seqType\E", "$name", "-----",
1570 ('-'x ((CORE::length $name) + 2)) );
1571 printf( $OUT "%13s: %d\n", "START", $self->{$mem.'Start'} );
1572 printf( $OUT "%13s: %d\n", "STOP", $self->{$mem.'Stop'} );
1573 printf( $OUT "%13s: \n", "SEQ" );
1574 foreach( @{ $self->{$mem.'Seq'}} ) {
1575 printf( $OUT "%15s%s\n", "", $_);
1576 }
1577 }
1578
1579
1580 =head2 _display_matches
1581
1582 Usage : n/a; called automatically by display()
1583 Purpose : Display information about identical and conserved positions
1584 : within both the query and sbjct sequences.
1585 Example : n/a
1586 Argument :
1587 Returns : printf call.
1588 Status : Experimental
1589 Comments : For more control, use seq_inds().
1590
1591 See Also : L<display()|display>, L<seq_inds()|seq_inds>, L<_display_seq()|_display_seq>,
1592
1593 =cut
1594
1595 #--------------------
1596 sub _display_matches {
1597 #--------------------
1598 my( $self, $queryName, $sbjctName, $OUT) = @_;
1599 my($resNum, $count);
1600
1601 $self->_set_residues() unless defined $self->{'_identicalRes_query'};
1602
1603 printf( $OUT "\n%10s: \n%10s\n", "HITS", "-----" );
1604 foreach( @{ $self->{'_matchSeq'}} ) {
1605 printf( $OUT "%15s%s\n", "", $_ );
1606 }
1607
1608 print $OUT "\n\U$queryName\E\n------------\n";
1609 printf( $OUT "\n%5s%s:\n%5s%s\n\t", "", "IDENTICAL RESIDUES IN $queryName (n=$self->{'_numIdentical'})",
1610 "", "--------------------------------------------" );
1611 $count = 0;
1612 foreach $resNum ( sort keys %{ $self->{'_identicalRes_query' }} ) {
1613 $count++;
1614 print $OUT "$resNum";
1615 $count > 0 and print $OUT +( $count % 15 ? ", " : "\n\t");
1616 }
1617
1618 print $OUT "\n";
1619
1620 my $justConserved = ($self->{'_numConserved'})-($self->{'_numIdentical'});
1621 printf( $OUT "\n%5s%s:\n%5s%s\n\t", "","CONSERVED RESIDUES IN $queryName (n=$justConserved)",
1622 "", "--------------------------------------------" );
1623 $count = 0;
1624 foreach $resNum ( sort keys %{ $self->{'_conservedRes_query' }} ) {
1625 $count++;
1626 print $OUT "$resNum";
1627 $count > 0 and print $OUT +( $count % 15 ? ", " : "\n\t");
1628 }
1629
1630
1631 print $OUT "\n\n\U$sbjctName\E\n------------\n";
1632 printf( $OUT "\n%5s%s:\n%5s%s\n\t", "", "IDENTICAL RESIDUES IN $sbjctName (n=$self->{'_numIdentical'})",
1633 "", "--------------------------------------------" );
1634 $count = 0;
1635 foreach $resNum ( sort keys %{ $self->{'_identicalRes_sbjct' }} ) {
1636 $count++;
1637 print $OUT "$resNum";
1638 $count > 0 and print $OUT +( $count % 15 ? ", " : "\n\t");
1639 }
1640
1641 print $OUT "\n";
1642 $justConserved = ($self->{'_numConserved'})-($self->{'_numIdentical'});
1643 printf( $OUT "\n%5s%s:\n%5s%s\n\t", "","CONSERVED RESIDUES IN $sbjctName (n=$justConserved)",
1644 "", "--------------------------------------------" );
1645 $count = 0;
1646 foreach $resNum ( sort keys %{ $self->{'_conservedRes_sbjct' }} ) {
1647 $count++;
1648 print $OUT "$resNum";
1649 $count > 0 and print $OUT +( $count % 15 ? ", " : "\n\t");
1650 }
1651 }
1652
1653
1654
1655
1656 =head2 homol_data
1657
1658 Usage : $data = $hsp_object->homo_data( %named_params );
1659 Purpose : Gets similarity data for a single HSP.
1660 Returns : String:
1661 : "Homology data" for each HSP is in the format:
1662 : "<integer> <start> <stop>"
1663 : where integer is the value returned by homol_score().
1664 Argument : Named params: (UPPER OR LOWERCASE TAGS)
1665 : currently just one param is used:
1666 : -SEQ =>'query' or 'sbjct'
1667 Throws : n/a
1668 Status : Experimental
1669 Comments : This is a very experimental method used for obtaining a
1670 : coarse indication of:
1671 : 1) how strong the similarity is between the sequences in the HSP,
1672 : 3) the endpoints of the alignment (sequence monomer numbers)
1673
1674 See Also : L<homol_score()|homol_score>, B<Bio::Tools::Blast.::homol_data()>, B<Bio::Tools::Blast::Sbjct::homol_data()>
1675
1676 =cut
1677
1678 #---------------
1679 sub homol_data {
1680 #---------------
1681 my ($self, %param) = @_;
1682 my $seq = $param{-SEQ} || $param{'-seq'} || 'sbjct'; # 'query' or 'sbjct'
1683 my $homolScore = $self->homol_score();
1684 # Sensitive to member name changes.
1685 $seq = "_\L$seq\E";
1686
1687 $self->_set_seq_data() unless $self->{'_set_seq_data'};
1688 return ( $homolScore.' '.$self->{$seq.'Start'}.' '.$self->{$seq.'Stop'});
1689 }
1690
1691
1692 =head2 homol_score
1693
1694 Usage : $self->homol_score();
1695 Purpose : Get a homology score (integer 1 - 3) as a coarse representation of
1696 : the strength of the similarity independent of sequence composition.
1697 : Based on the Blast bit score.
1698 Example : $hscore = $hsp->homol_score();
1699 Returns : Integer
1700 Argument : n/a
1701 Throws : n/a
1702 Status : Experimental
1703 Comments : See @Bio::Tools::Blast::HSP::SCORE_CUTOFFS for the specific values.
1704 : Currently, BIT_SCORE HOMOL_SCORE
1705 : --------- -----------
1706 : >=100 --> 3
1707 : 30-100 --> 2
1708 : < 30 --> 1
1709
1710 See Also : L<homol_data()|homol_data>
1711
1712 =cut
1713
1714 #----------------
1715 sub homol_score {
1716 #----------------
1717 my $self = shift;
1718
1719 if( $self->{'_bits'} >= $SCORE_CUTOFFS[0] ) { 1 }
1720 elsif($self->{'_bits'} < $SCORE_CUTOFFS[0] and
1721 $self->{'_bits'} >= $SCORE_CUTOFFS[1] ) { 2 }
1722 else { 3 }
1723 }
1724
1725
1726 #####################################################################################
1727 ## CLASS METHODS ##
1728 #####################################################################################
1729
1730 =head1 CLASS METHODS
1731
1732 =head2 collapse_nums
1733
1734 Usage : @cnums = collapse_nums( @numbers );
1735 Purpose : Collapses a list of numbers into a set of ranges of consecutive terms:
1736 : Useful for condensing long lists of consecutive numbers.
1737 : EXPANDED:
1738 : 1 2 3 4 5 6 10 12 13 14 15 17 18 20 21 22 24 26 30 31 32
1739 : COLLAPSED:
1740 : 1-6 10 12-15 17 18 20-22 24 26 30-32
1741 Argument : List of numbers and sorted numerically.
1742 Returns : List of numbers mixed with ranges of numbers (see above).
1743 Throws : n/a
1744 Comments : Probably belongs in a more general utility class.
1745
1746 See Also : L<seq_inds()|seq_inds>
1747
1748 =cut
1749
1750 #------------------
1751 sub collapse_nums {
1752 #------------------
1753 # This is not the slickest connectivity algorithm, but will do for now.
1754 my @a = @_;
1755 my ($from, $to, $i, @ca, $consec);
1756
1757 $consec = 0;
1758 for($i=0; $i < @a; $i++) {
1759 not $from and do{ $from = $a[$i]; next; };
1760 if($a[$i] == $a[$i-1]+1) {
1761 $to = $a[$i];
1762 $consec++;
1763 } else {
1764 if($consec == 1) { $from .= ",$to"; }
1765 else { $from .= $consec>1 ? "\-$to" : ""; }
1766 push @ca, split(',', $from);
1767 $from = $a[$i];
1768 $consec = 0;
1769 $to = undef;
1770 }
1771 }
1772 if(defined $to) {
1773 if($consec == 1) { $from .= ",$to"; }
1774 else { $from .= $consec>1 ? "\-$to" : ""; }
1775 }
1776 push @ca, split(',', $from) if $from;
1777
1778 @ca;
1779 }
1780
1781
1782 1;
1783 __END__
1784
1785 #####################################################################################
1786 # END OF CLASS
1787 #####################################################################################
1788
1789 =head1 FOR DEVELOPERS ONLY
1790
1791 =head2 Data Members
1792
1793 Information about the various data members of this module is provided for those
1794 wishing to modify or understand the code. Two things to bear in mind:
1795
1796 =over 4
1797
1798 =item 1 Do NOT rely on these in any code outside of this module.
1799
1800 All data members are prefixed with an underscore to signify that they are private.
1801 Always use accessor methods. If the accessor doesn't exist or is inadequate,
1802 create or modify an accessor (and let me know, too!).
1803
1804 =item 2 This documentation may be incomplete and out of date.
1805
1806 It is easy for these data member descriptions to become obsolete as
1807 this module is still evolving. Always double check this info and search
1808 for members not described here.
1809
1810 =back
1811
1812 An instance of Bio::Tools::Blast::HSP.pm is a blessed reference to a hash containing
1813 all or some of the following fields:
1814
1815 FIELD VALUE
1816 --------------------------------------------------------------
1817 (member names are mostly self-explanatory)
1818
1819 _score :
1820 _bits :
1821 _p :
1822 _n : Integer. The 'N' value listed in parenthesis with P/Expect value:
1823 : e.g., P(3) = 1.2e-30 ---> (N = 3).
1824 : Not defined in NCBI Blast2 with gaps.
1825 : To obtain the number of HSPs, use Bio::Tools::Blast::Sbjct::num_hsps().
1826 _expect :
1827 _queryLength :
1828 _queryGaps :
1829 _queryStart :
1830 _queryStop :
1831 _querySeq :
1832 _sbjctLength :
1833 _sbjctGaps :
1834 _sbjctStart :
1835 _sbjctStop :
1836 _sbjctSeq :
1837 _matchSeq : String. Contains the symbols between the query and sbjct lines
1838 which indicate identical (letter) and conserved ('+') matches
1839 or a mismatch (' ').
1840 _numIdentical :
1841 _numConserved :
1842 _identicalRes_query :
1843 _identicalRes_sbjct :
1844 _conservedRes_query :
1845 _conservedRes_sbjct :
1846 _match_indent : The number of leading space characters on each line containing
1847 the match symbols. _match_indent is 13 in this example:
1848 Query: 285 QNSAPWGLARISHRERLNLGSFNKYLYDDDAG
1849 Q +APWGLARIS G+ + Y YD+ AG
1850 ^^^^^^^^^^^^^
1851
1852 INHERITED DATA MEMBERS
1853
1854 _name : From Bio::Root::Object.pm.
1855 :
1856 _parent : From Bio::Root::Object.pm. This member contains a reference to the
1857 : Bio::Tools::Blast::Sbjct.pm object to which this hit belongs.
1858
1859
1860 =cut
1861
1862 1;
1863