Bio/Search/SearchUtils.pm

   1 =head1 NAME
   2
   3 Bio::Search::SearchUtils - Utility functions for Bio::Search:: objects
   4
   5 =head1 SYNOPSIS
   6
   7   # This module is just a collection of subroutines, not an object.
   8
   9 =head1 DESCRIPTION
  10
  11 The SearchUtils.pm module is a collection of subroutines used
  12 primarily by Bio::Search::Hit::HitI objects for some of the additional
  13 functionality, such as HSP tiling. Right now, the SearchUtils is just
  14 a collection of methods, not an object.
  15
  16 =head1 AUTHOR
  17
  18 Steve Chervitz E<lt>sac@bioperl.orgE<gt>
  19
  20 =head1 CONTRIBUTORS
  21
  22 Sendu Bala, bix@sendu.me.uk
  23
  24 =cut
  25
  26 package Bio::Search::SearchUtils;
  27 use Bio::Root::Version;
  28
  29 use strict;
  30
  31 =head2 tile_hsps
  32
  33  Usage     : tile_hsps( $sbjct );
  34            : This is called automatically by methods in Bio::Search::Hit::GenericHit
  35            : that rely on having tiled data.
  36            :
  37            : If you are interested in getting data about the constructed HSP contigs:
  38            : my ($qcontigs, $scontigs) = Bio::Search::SearchUtils::tile_hsps($hit);
  39            : if (ref $qcontigs) {
  40            :    print STDERR "Query contigs:\n";
  41            :    foreach (@{$qcontigs}) {
  42            :         print "contig start is $_->{'start'}\n";
  43            :         print "contig stop is $_->{'stop'}\n";
  44            :    }
  45            : }
  46            : See below for more information about the contig data structure.
  47            :
  48  Purpose   : Collect statistics about the aligned sequences in a set of HSPs.
  49            : Calculates the following data across all HSPs:
  50            :    -- total alignment length
  51            :    -- total identical residues
  52            :    -- total conserved residues
  53  Returns   : If there was only a single HSP (so no tiling was necessary)
  54                tile_hsps() returns a list of two non-zero integers.
  55              If there were multiple HSP,
  56                tile_hsps() returns a list of two array references containin HSP contig data.
  57              The first array ref contains a list of HSP contigs on the query sequence.
  58              The second array ref contains a list of HSP contigs on the subject sequence.
  59              Each contig is a hash reference with the following data fields:
  60                'start' => start coordinate of the contig
  61                'stop'  => start coordinate of the contig
  62                'iden'  => number of identical residues in the contig
  63                'cons'  => number of conserved residues in the contig
  64                'strand'=> strand of the contig
  65                'frame' => frame of the contig
  66  Argument  : A Bio::Search::Hit::HitI object
  67  Throws    : n/a
  68  Comments  :
  69            : This method performs more careful summing of data across
  70            : all HSPs in the Sbjct object. Only HSPs that are in the same strand
  71            : and frame are tiled. Simply summing the data from all HSPs
  72            : in the same strand and frame will overestimate the actual
  73            : length of the alignment if there is overlap between different HSPs
  74            : (often the case).
  75            :
  76            : The strategy is to tile the HSPs and sum over the
  77            : contigs, collecting data separately from overlapping and
  78            : non-overlapping regions of each HSP. To facilitate this, the
  79            : HSP.pm object now permits extraction of data from sub-sections
  80            : of an HSP.
  81            :
  82            : Additional useful information is collected from the results
  83            : of the tiling. It is possible that sub-sequences in
  84            : different HSPs will overlap significantly. In this case, it
  85            : is impossible to create a single unambiguous alignment by
  86            : concatenating the HSPs. The ambiguity may indicate the
  87            : presence of multiple, similar domains in one or both of the
  88            : aligned sequences. This ambiguity is recorded using the
  89            : ambiguous_aln() method.
  90            :
  91            : This method does not attempt to discern biologically
  92            : significant vs. insignificant overlaps. The allowable amount of
  93            : overlap can be set with the overlap() method or with the -OVERLAP
  94            : parameter used when constructing the Hit object.
  95            :
  96            : For a given hit, both the query and the sbjct sequences are
  97            : tiled independently.
  98            :
  99            :    -- If only query sequence HSPs overlap,
 100            :          this may suggest multiple domains in the sbjct.
 101            :    -- If only sbjct sequence HSPs overlap,
 102            :          this may suggest multiple domains in the query.
 103            :    -- If both query & sbjct sequence HSPs overlap,
 104            :          this suggests multiple domains in both.
 105            :    -- If neither query & sbjct sequence HSPs overlap,
 106            :          this suggests either no multiple domains in either
 107            :          sequence OR that both sequences have the same
 108            :          distribution of multiple similar domains.
 109            :
 110            : This method can deal with the special case of when multiple
 111            : HSPs exactly overlap.
 112            :
 113            : Efficiency concerns:
 114            :  Speed will be an issue for sequences with numerous HSPs.
 115            :
 116  Bugs      : Currently, tile_hsps() does not properly account for
 117            : the number of non-tiled but overlapping HSPs, which becomes a problem
 118            : as overlap() grows. Large values overlap() may thus lead to
 119            : incorrect statistics for some hits. For best results, keep overlap()
 120            : below 5 (DEFAULT IS 2). For more about this, see the "HSP Tiling and
 121            : Ambiguous Alignments" section in L<Bio::Search::Hit::GenericHit>.
 122
 123 See Also   : L<_adjust_contigs>(), L<Bio::Search::Hit::GenericHit|Bio::Search::Hit::GenericHit>
 124
 125 =cut
 126
 127 #--------------
 128 sub tile_hsps {
 129 #--------------
 130     my $sbjct = shift;
 131
 132     #print STDERR "Calling tile_hsps(): $sbjct\n";
 133     #$sbjct->verbose(1);  # to activate debugging
 134     $sbjct->tiled_hsps(1);
 135
 136     if( $sbjct->num_hsps == 0 || $sbjct->n == 0 ) {
 137         #print STDERR "_tile_hsps(): no hsps, nothing to tile! (", $sbjct->num_hsps, ")\n";
 138         _warn_about_no_hsps($sbjct);
 139         return (undef, undef);
 140
 141     } elsif( $sbjct->n == 1 or $sbjct->num_hsps == 1) {
 142         ## Simple summation scheme. Valid if there is only one HSP.
 143         #print STDERR "_tile_hsps(): single HSP, easy stats.\n";
 144         my $hsp = $sbjct->hsp;
 145         $sbjct->length_aln('query', $hsp->length('query'));
 146         $sbjct->length_aln('hit', $hsp->length('sbjct'));
 147         $sbjct->length_aln('total', $hsp->length('total'));
 148         $sbjct->matches( $hsp->matches() );
 149         $sbjct->gaps('query', $hsp->gaps('query'));
 150         $sbjct->gaps('sbjct', $hsp->gaps('sbjct'));
 151
 152         _adjust_length_aln($sbjct);
 153         return (1, 1);
 154     } else {
 155         #print STDERR "Sbjct: _tile_hsps: summing multiple HSPs\n";
 156         $sbjct->length_aln('query', 0);
 157         $sbjct->length_aln('sbjct', 0);
 158         $sbjct->length_aln('total', 0);
 159         $sbjct->matches( 0, 0);
 160         $sbjct->gaps('query', 0);
 161         $sbjct->gaps('hit', 0);
 162     }
 163
 164     ## More than one HSP. Must tile HSPs.
 165 #    print "\nTiling HSPs for $sbjct\n";
 166     my($hsp, $qstart, $qstop, $sstart, $sstop);
 167     my($frame, $strand, $qstrand, $sstrand);
 168     my(@qcontigs, @scontigs);
 169     my $qoverlap = 0;
 170     my $soverlap = 0;
 171     my $max_overlap = $sbjct->overlap;
 172     my $hit_qgaps = 0;
 173     my $hit_sgaps = 0;
 174     my $hit_len_aln = 0;
 175     my %start_stop;
 176     my $v = $sbjct->verbose;
 177     foreach $hsp ( $sbjct->hsps() ) {
 178         #$sbjct->debug( sprintf("  HSP: %s %d..%d\n",$hsp->query->seq_id, $hsp->query->start, $hsp->hit->end)) if $v > 0; #$hsp->str('query');
 179 #       printf "  Length = %d; Identical = %d; Conserved = %d; Conserved(1-10): %d",$hsp->length, $hsp->length(-TYPE=>'iden'),
 180 #       $hsp->length(-TYPE=>'cons'),
 181 #       $hsp->length(-TYPE=>'cons',
 182 #                    -START=>0,-STOP=>10);
 183
 184         ($qstart, $qstop) = $hsp->range('query');
 185         ($sstart, $sstop) = $hsp->range('sbjct');
 186         $frame = $hsp->frame('hit');
 187         $frame = -1 unless defined $frame;
 188
 189         ($qstrand, $sstrand) = ($hsp->query->strand,
 190                                 $hsp->hit->strand);
 191
 192         # Note: No correction for overlap.
 193
 194         my ($qgaps, $sgaps)  = ($hsp->gaps('query'), $hsp->gaps('hit'));
 195         $hit_qgaps += $qgaps;
 196         $hit_sgaps += $sgaps;
 197         $hit_len_aln += $hsp->length;
 198
 199         ## Collect contigs in the query sequence.
 200         $qoverlap = &_adjust_contigs('query', $hsp, $qstart, $qstop,
 201                                      \@qcontigs, $max_overlap, $frame,
 202                                      $qstrand);
 203
 204         ## Collect contigs in the sbjct sequence
 205         #  (needed for domain data and gapped Blast).
 206         $soverlap = &_adjust_contigs('sbjct', $hsp, $sstart, $sstop,
 207                                      \@scontigs, $max_overlap, $frame,
 208                                      $sstrand);
 209
 210         ## Collect overall start and stop data for query and
 211         #  sbjct over all HSPs.
 212         unless ( defined $start_stop{'qstart'} ) {
 213             $start_stop{'qstart'}  = $qstart;
 214             $start_stop{'qstop'}   = $qstop;
 215             $start_stop{'sstart'}  = $sstart;
 216             $start_stop{'sstop'}   = $sstop;
 217         } else {
 218             $start_stop{'qstart'} = ($qstart < $start_stop{'qstart'} ?
 219                                      $qstart : $start_stop{'qstart'} );
 220             $start_stop{'qstop'}  = ($qstop  > $start_stop{'qstop'} ?
 221                                      $qstop  : $start_stop{'qstop'} );
 222             $start_stop{'sstart'} = ($sstart < $start_stop{'sstart'} ?
 223                                      $sstart : $start_stop{'sstart'} );
 224             $start_stop{'sstop'}  = ($sstop  > $start_stop{'sstop'} ?
 225                                      $sstop  : $start_stop{'sstop'} );
 226         }
 227     }
 228
 229     # Store the collected data in the Hit object
 230     $sbjct->gaps('query', $hit_qgaps);
 231     $sbjct->gaps('hit', $hit_sgaps);
 232     $sbjct->length_aln('total', $hit_len_aln);
 233
 234     $sbjct->start('query',$start_stop{'qstart'});
 235     $sbjct->end('query', $start_stop{'qstop'});
 236     $sbjct->start('hit', $start_stop{'sstart'});
 237     $sbjct->end('hit', $start_stop{'sstop'});
 238     ## Collect data across the collected contigs.
 239
 240     #$sbjct->debug( "\nQUERY CONTIGS:\n"."  gaps = $sbjct->{'_gaps_query'}\n");
 241
 242     # Account for strand/frame.
 243     # Strategy: collect data on a per strand+frame basis and
 244     #           save the most significant one.
 245     my (%qctg_dat);
 246     foreach (@qcontigs) {
 247         ($frame, $strand) = ($_->{'frame'}, $_->{'strand'});
 248
 249         if( $v > 0 ) {
 250             #$sbjct->debug(sprintf( "$frame/$strand len is getting %d for %d..%d\n",
 251             #                      ($_->{'stop'} - $_->{'start'} + 1), $_->{'start'}, $_->{'stop'}));
 252         }
 253
 254         $qctg_dat{ "$frame$strand" }->{'length_aln_query'} += $_->{'stop'} - $_->{'start'} + 1;
 255         $qctg_dat{ "$frame$strand" }->{'totalIdentical'}   += $_->{'iden'};
 256         $qctg_dat{ "$frame$strand" }->{'totalConserved'}   += $_->{'cons'};
 257         $qctg_dat{ "$frame$strand" }->{'qstrand'}   = $strand;
 258     }
 259
 260     # Find longest contig.
 261     my @sortedkeys = sort { $qctg_dat{$b}->{'length_aln_query'}
 262                             <=> $qctg_dat{$a}->{'length_aln_query'} }
 263     keys %qctg_dat;
 264
 265     # Save the largest to the sbjct:
 266     my $longest = $sortedkeys[0];
 267     #$sbjct->debug( "longest is ". $qctg_dat{ $longest }->{'length_aln_query'}. "\n");
 268     $sbjct->length_aln('query', $qctg_dat{ $longest }->{'length_aln_query'});
 269     $sbjct->matches($qctg_dat{ $longest }->{'totalIdentical'},
 270                     $qctg_dat{ $longest }->{'totalConserved'});
 271     $sbjct->strand('query', $qctg_dat{ $longest }->{'qstrand'});
 272
 273     ## Collect data for sbjct contigs. Important for gapped Blast.
 274     ## The totalIdentical and totalConserved numbers will be the same
 275     ## as determined for the query contigs.
 276
 277     #$sbjct->debug( "\nSBJCT CONTIGS:\n"."  gaps = ". $sbjct->gaps('sbjct'). "\n");
 278     my (%sctg_dat);
 279     foreach(@scontigs) {
 280         #$sbjct->debug("  sbjct contig: $_->{'start'} - $_->{'stop'}\n".
 281         #            "         iden = $_->{'iden'}; cons = $_->{'cons'}\n");
 282         ($frame, $strand) = ($_->{'frame'}, $_->{'strand'});
 283         $sctg_dat{ "$frame$strand" }->{'length_aln_sbjct'}   += $_->{'stop'} - $_->{'start'} + 1;
 284         $sctg_dat{ "$frame$strand" }->{'frame'}  = $frame;
 285         $sctg_dat{ "$frame$strand" }->{'sstrand'}  = $strand;
 286     }
 287
 288     @sortedkeys = sort { $sctg_dat{ $b }->{'length_aln_sbjct'}
 289                          <=> $sctg_dat{ $a }->{'length_aln_sbjct'}
 290                      } keys %sctg_dat;
 291
 292     # Save the largest to the sbjct:
 293     $longest = $sortedkeys[0];
 294
 295     $sbjct->length_aln('sbjct', $sctg_dat{ $longest }->{'length_aln_sbjct'});
 296     $sbjct->frame( $sctg_dat{ $longest }->{'frame'} );
 297     $sbjct->strand('hit', $sctg_dat{ $longest }->{'sstrand'});
 298
 299     if($qoverlap) {
 300         if($soverlap) { $sbjct->ambiguous_aln('qs');
 301                         #$sbjct->debug("\n*** AMBIGUOUS ALIGNMENT: Query and Sbjct\n\n");
 302                     }
 303         else { $sbjct->ambiguous_aln('q');
 304                #$sbjct->debug( "\n*** AMBIGUOUS ALIGNMENT: Query\n\n");
 305            }
 306     } elsif($soverlap) {
 307         $sbjct->ambiguous_aln('s');
 308         #$sbjct->debug( "\n*** AMBIGUOUS ALIGNMENT: Sbjct\n\n");
 309     }
 310
 311     _adjust_length_aln($sbjct);
 312
 313     return ( [@qcontigs], [@scontigs] );
 314 }
 315
 316
 317
 318 # Title    : _adjust_length_aln
 319 # Usage    : n/a; internal use only; called by tile_hsps.
 320 # Purpose  : Adjust length of aligment based on BLAST flavor.
 321 # Comments : See comments in logica_length()
 322 sub _adjust_length_aln {
 323     my $sbjct = shift;
 324     my $algo = $sbjct->algorithm;
 325     my $hlen = $sbjct->length_aln('sbjct');
 326     my $qlen = $sbjct->length_aln('query');
 327
 328     $sbjct->length_aln('sbjct', logical_length($algo, 'sbjct', $hlen));
 329     $sbjct->length_aln('query', logical_length($algo, 'query', $qlen));
 330 }
 331
 332 =head2 logical_length
 333
 334  Usage     : logical_length( $alg_name, $seq_type, $length );
 335  Purpose   : Determine the logical length of an aligned sequence based on
 336            : algorithm name and sequence type.
 337  Returns   : integer representing the logical aligned length.
 338  Argument  : $alg_name = name of algorigthm (e.g., blastx, tblastn)
 339            : $seq_type = type of sequence (e.g., query or hit)
 340            : $length = physical length of the sequence in the alignment.
 341  Throws    : n/a
 342  Comments  : This function is used to account for the fact that number of identities
 343              and conserved residues is reported in peptide space while the query
 344              length (in the case of BLASTX and TBLASTX) and/or the hit length
 345              (in the case of TBLASTN and TBLASTX) are in nucleotide space.
 346              The adjustment affects the values reported by the various frac_XXX
 347              methods in GenericHit and GenericHSP.
 348
 349 =cut
 350
 351 sub logical_length {
 352     my ($algo, $type, $len) = @_;
 353     my $logical = $len;
 354     if($algo =~ /^(?:PSI)?T(?:BLASTN|FAST(?:X|Y|XY))/oi ) {
 355         $logical = $len/3 if $type =~ /sbjct|hit|tot/i;
 356     } elsif($algo =~ /^(?:BLASTX|FAST(?:X|Y|XY))/oi ) {
 357         $logical = $len/3 if $type =~ /query|tot/i;
 358     } elsif($algo =~ /^TBLASTX/oi ) {
 359         $logical = $len/3;
 360     }
 361     return $logical;
 362 }
 363
 364
 365 #=head2 _adjust_contigs
 366 #
 367 # Usage     : n/a; internal function called by tile_hsps
 368 # Purpose   : Builds HSP contigs for a given BLAST hit.
 369 #           : Utility method called by _tile_hsps()
 370 # Returns   :
 371 # Argument  :
 372 # Throws    : Exceptions propagated from Bio::Search::Hit::BlastHSP::matches()
 373 #           : for invalid sub-sequence ranges.
 374 # Status    : Experimental
 375 # Comments  : This method does not currently support gapped alignments.
 376 #           : Also, it does not keep track of the number of HSPs that
 377 #           : overlap within the amount specified by overlap().
 378 #           : This will lead to significant tracking errors for large
 379 #           : overlap values.
 380 #
 381 #See Also   : L<tile_hsps>(), L<Bio::Search::Hit::BlastHSP::matches|Bio::Search::Hit::BlastHSP>
 382 #
 383 #=cut
 384
 385 sub _adjust_contigs {
 386     my ($seqType, $hsp, $start, $stop, $contigs_ref,
 387         $max_overlap, $frame, $strand) = @_;
 388     my $overlap = 0;
 389     my ($numID, $numCons);
 390
 391     foreach (@$contigs_ref) {
 392         # Don't merge things unless they have matching strand/frame.
 393         next unless ($_->{'frame'} == $frame && $_->{'strand'} == $strand);
 394
 395         # Test special case of a nested HSP. Skip it.
 396         if ($start >= $_->{'start'} && $stop <= $_->{'stop'}) {
 397             $overlap = 1;
 398             next;
 399         }
 400
 401         # Test for overlap at beginning of contig, or precedes consecutively
 402         if ($start < $_->{'start'} && $stop >= ($_->{'start'} + $max_overlap - 1)) {
 403             eval {
 404                 ($numID, $numCons) = $hsp->matches(-SEQ   =>$seqType,
 405                                -START => $start,
 406                                -STOP  => $_->{'start'} - 1);
 407             };
 408             if($@) { warn "\a\n$@\n"; }
 409             else {
 410                 $_->{'start'} = $start; # Assign a new start coordinate to the contig
 411                 $_->{'iden'} += $numID; # and add new data to #identical, #conserved.
 412                 $_->{'cons'} += $numCons;
 413                 push(@{$_->{hsps}}, $hsp);
 414                 $overlap     = 1;
 415             }
 416         }
 417
 418         # Test for overlap at end of contig, or follows consecutively
 419         if ($stop > $_->{'stop'} and $start <= ($_->{'stop'} - $max_overlap + 1)) {
 420             eval {
 421                 ($numID,$numCons) = $hsp->matches(-SEQ   =>$seqType,
 422                               -START => $_->{'stop'} + 1,
 423                               -STOP  => $stop);
 424             };
 425             if($@) { warn "\a\n$@\n"; }
 426             else {
 427                 $_->{'stop'}  = $stop; # Assign a new stop coordinate to the contig
 428                 $_->{'iden'} += $numID; # and add new data to #identical, #conserved.
 429                 $_->{'cons'} += $numCons;
 430                 push(@{$_->{hsps}}, $hsp);
 431                 $overlap    = 1;
 432             }
 433         }
 434
 435         last if $overlap;
 436     }
 437
 438     if ($overlap && @$contigs_ref > 1) {
 439         ## Merge any contigs that now overlap
 440         my $max = $#{$contigs_ref};
 441         for my $i (0..$max) {
 442             ${$contigs_ref}[$i] || next;
 443             my ($i_start, $i_stop) = (${$contigs_ref}[$i]->{start}, ${$contigs_ref}[$i]->{stop});
 444
 445             for my $u ($i+1..$max) {
 446                 ${$contigs_ref}[$u] || next;
 447                 my ($u_start, $u_stop) = (${$contigs_ref}[$u]->{start}, ${$contigs_ref}[$u]->{stop});
 448
 449                 if ($u_start < $i_start && $u_stop >= ($i_start + $max_overlap - 1)) {
 450                     # find the hsps within the contig that have sequence
 451                     # extending before $i_start
 452                     my ($ids, $cons) = (0, 0);
 453                     my $use_start = $i_start;
 454                     foreach my $hsp (sort { $b->end($seqType) <=> $a->end($seqType) } @{${$contigs_ref}[$u]->{hsps}}) {
 455                         my $hsp_start = $hsp->start($seqType);
 456                         $hsp_start < $use_start || next;
 457
 458                         my ($these_ids, $these_cons);
 459                         eval {
 460                             ($these_ids, $these_cons) = $hsp->matches(-SEQ => $seqType, -START => $hsp_start, -STOP => $use_start - 1);
 461                         };
 462                         if($@) { warn "\a\n$@\n"; }
 463                         else {
 464                             $ids  += $these_ids;
 465                             $cons += $these_cons;
 466                         }
 467
 468                         last if $hsp_start == $u_start;
 469                         $use_start = $hsp_start;
 470                     }
 471                     ${$contigs_ref}[$i]->{start} = $u_start;
 472                     ${$contigs_ref}[$i]->{'iden'} += $ids;
 473                     ${$contigs_ref}[$i]->{'cons'} += $cons;
 474                     push(@{${$contigs_ref}[$i]->{hsps}}, @{${$contigs_ref}[$u]->{hsps}});
 475
 476                     ${$contigs_ref}[$u] = undef;
 477                 }
 478                 elsif ($u_stop > $i_stop && $u_start <= ($i_stop - $max_overlap + 1)) {
 479                     # find the hsps within the contig that have sequence
 480                     # extending beyond $i_stop
 481                     my ($ids, $cons) = (0, 0);
 482                     my $use_stop = $i_stop;
 483                     foreach my $hsp (sort { $a->start($seqType) <=> $b->start($seqType) } @{${$contigs_ref}[$u]->{hsps}}) {
 484                         my $hsp_end = $hsp->end($seqType);
 485                         $hsp_end > $use_stop || next;
 486
 487                         my ($these_ids, $these_cons);
 488                         eval {
 489                             ($these_ids, $these_cons) = $hsp->matches(-SEQ => $seqType, -START => $use_stop + 1, -STOP => $hsp_end);
 490                         };
 491                         if($@) { warn "\a\n$@\n"; }
 492                         else {
 493                             $ids  += $these_ids;
 494                             $cons += $these_cons;
 495                         }
 496
 497                         last if $hsp_end == $u_stop;
 498                         $use_stop = $hsp_end;
 499                     }
 500                     ${$contigs_ref}[$i]->{'stop'}  = $u_stop;
 501                     ${$contigs_ref}[$i]->{'iden'} += $ids;
 502                     ${$contigs_ref}[$i]->{'cons'} += $cons;
 503                     push(@{${$contigs_ref}[$i]->{hsps}}, @{${$contigs_ref}[$u]->{hsps}});
 504
 505                     ${$contigs_ref}[$u] = undef;
 506                 }
 507                 elsif ($u_start >= $i_start && $u_stop <= $i_stop) {
 508                     # nested, drop this contig
 509                     #*** ideally we might do some magic to keep the stats of the
 510                     #    better hsp...
 511                     ${$contigs_ref}[$u] = undef;
 512                 }
 513             }
 514         }
 515
 516         my @merged;
 517         foreach (@$contigs_ref) {
 518             push(@merged, $_ || next);
 519         }
 520         @{$contigs_ref} = @merged;
 521     }
 522     elsif (! $overlap) {
 523         ## If there is no overlap, add the complete HSP data.
 524         ($numID,$numCons) = $hsp->matches(-SEQ=>$seqType);
 525         push @$contigs_ref, {'start' =>$start, 'stop' =>$stop,
 526                              'iden'  =>$numID, 'cons' =>$numCons,
 527                              'strand'=>$strand,'frame'=>$frame,'hsps'=>[$hsp]};
 528     }
 529
 530     return $overlap;
 531 }
 532
 533 =head2 get_exponent
 534
 535  Usage     : &get_exponent( number );
 536  Purpose   : Determines the power of 10 exponent of an integer, float,
 537            : or scientific notation number.
 538  Example   : &get_exponent("4.0e-206");
 539            : &get_exponent("0.00032");
 540            : &get_exponent("10.");
 541            : &get_exponent("1000.0");
 542            : &get_exponent("e+83");
 543  Argument  : Float, Integer, or scientific notation number
 544  Returns   : Integer representing the exponent part of the number (+ or -).
 545            : If argument == 0 (zero), return value is "-999".
 546  Comments  : Exponents are rounded up (less negative) if the mantissa is >= 5.
 547            : Exponents are rounded down (more negative) if the mantissa is <= -5.
 548
 549 =cut
 550
 551 sub get_exponent {
 552     my $data = shift;
 553
 554     my($num, $exp) = split /[eE]/, $data;
 555
 556     if( defined $exp) {
 557         $num = 1 if not $num;
 558         $num >= 5 and $exp++;
 559         $num <= -5 and $exp--;
 560     } elsif( $num == 0) {
 561         $exp = -999;
 562     } elsif( not $num =~ /\./) {
 563         $exp = CORE::length($num) -1;
 564     } else {
 565         $exp = 0;
 566         $num .= '0' if $num =~ /\.$/;
 567         my ($c);
 568         my $rev = 0;
 569         if($num !~ /^0/) {
 570             $num = reverse($num);
 571             $rev = 1;
 572         }
 573         do { $c = chop($num);
 574              $c == 0 && $exp++;
 575          } while( $c ne '.');
 576
 577         $exp = -$exp if $num == 0 and not $rev;
 578         $exp -= 1 if $rev;
 579     }
 580     return $exp;
 581 }
 582
 583 =head2 collapse_nums
 584
 585  Usage     : @cnums = collapse_nums( @numbers );
 586  Purpose   : Collapses a list of numbers into a set of ranges of consecutive terms:
 587            : Useful for condensing long lists of consecutive numbers.
 588            :  EXPANDED:
 589            :     1 2 3 4 5 6 10 12 13 14 15 17 18 20 21 22 24 26 30 31 32
 590            :  COLLAPSED:
 591            :     1-6 10 12-15 17 18 20-22 24 26 30-32
 592  Argument  : List of numbers sorted numerically.
 593  Returns   : List of numbers mixed with ranges of numbers (see above).
 594  Throws    : n/a
 595
 596 See Also   : L<Bio::Search::Hit::BlastHit::seq_inds()|Bio::Search::Hit::BlastHit>
 597
 598 =cut
 599
 600 sub collapse_nums {
 601 # This is probably not the slickest connectivity algorithm, but will do for now.
 602     my @a = @_;
 603     my ($from, $to, $i, @ca, $consec);
 604
 605     $consec = 0;
 606     for($i=0; $i < @a; $i++) {
 607         not $from and do{ $from = $a[$i]; next; };
 608     # pass repeated positions (gap inserts)
 609     next if $a[$i] == $a[$i-1];
 610         if($a[$i] == $a[$i-1]+1) {
 611             $to = $a[$i];
 612             $consec++;
 613         } else {
 614             if($consec == 1) { $from .= ",$to"; }
 615             else { $from .= $consec>1 ? "\-$to" : ""; }
 616             push @ca, split(',', $from);
 617             $from =  $a[$i];
 618             $consec = 0;
 619             $to = undef;
 620         }
 621     }
 622     if(defined $to) {
 623         if($consec == 1) { $from .= ",$to"; }
 624         else { $from .= $consec>1 ? "\-$to" : ""; }
 625     }
 626     push @ca, split(',', $from) if $from;
 627
 628     @ca;
 629 }
 630
 631
 632 =head2 strip_blast_html
 633
 634  Usage     : $boolean = &strip_blast_html( string_ref );
 635            : This method is exported.
 636  Purpose   : Removes HTML formatting from a supplied string.
 637            : Attempts to restore the Blast report to enable
 638            : parsing by Bio::SearchIO::blast.pm
 639  Returns   : Boolean: true if string was stripped, false if not.
 640  Argument  : string_ref = reference to a string containing the whole Blast
 641            :              report containing HTML formatting.
 642  Throws    : Croaks if the argument is not a scalar reference.
 643  Comments  : Based on code originally written by Alex Dong Li
 644            : (ali@genet.sickkids.on.ca).
 645            : This method does some Blast-specific stripping
 646            : (adds back a '>' character in front of each HSP
 647            : alignment listing).
 648            :
 649            : THIS METHOD IS VERY SENSITIVE TO BLAST FORMATTING CHANGES!
 650            :
 651            : Removal of the HTML tags and accurate reconstitution of the
 652            : non-HTML-formatted report is highly dependent on structure of
 653            : the HTML-formatted version. For example, it assumes that first
 654            : line of each alignment section (HSP listing) starts with a
 655            : <a name=..> anchor tag. This permits the reconstruction of the
 656            : original report in which these lines begin with a ">".
 657            : This is required for parsing.
 658            :
 659            : If the structure of the Blast report itself is not intended to
 660            : be a standard, the structure of the HTML-formatted version
 661            : is even less so. Therefore, the use of this method to
 662            : reconstitute parsable Blast reports from HTML-format versions
 663            : should be considered a temorary solution.
 664
 665 =cut
 666
 667 sub strip_blast_html {
 668       # This may not best way to remove html tags. However, it is simple.
 669       # it won't work under following conditions:
 670       #    1) if quoted > appears in a tag  (does this ever happen?)
 671       #    2) if a tag is split over multiple lines and this method is
 672       #       used to process one line at a time.
 673
 674     my ($string_ref) = shift;
 675
 676     ref $string_ref eq 'SCALAR' or
 677         croak ("Can't strip HTML: ".
 678                "Argument is should be a SCALAR reference not a ${\ref $string_ref}\n");
 679
 680     my $str = $$string_ref;
 681     my $stripped = 0;
 682
 683     # Removing "<a name =...>" and adding the '>' character for
 684     # HSP alignment listings.
 685     $str =~ s/(\A|\n)<a name ?=[^>]+> ?/>/sgi and $stripped = 1;
 686
 687     # Removing all "<>" tags.
 688     $str =~ s/<[^>]+>|&nbsp//sgi and $stripped = 1;
 689
 690     # Re-uniting any lone '>' characters.
 691     $str =~ s/(\A|\n)>\s+/\n\n>/sgi and $stripped = 1;
 692
 693     $$string_ref = $str;
 694     $stripped;
 695 }
 696
 697 =head2 result2hash
 698
 699  Title    : result2hash
 700  Usage    : my %data = &Bio::Search::SearchUtils($result)
 701  Function : converts ResultI data to simple hash
 702  Returns  : hash
 703  Args     : ResultI
 704  Note     : used mainly as a utility for running SearchIO tests
 705
 706 =cut
 707
 708 sub result2hash {
 709     my ($result) = @_;
 710     my %hash;
 711     $hash{'query_name'} = $result->query_name;
 712     my $hitcount = 1;
 713     my $hspcount = 1;
 714     foreach my $hit ( $result->hits ) {
 715         $hash{"hit$hitcount\_name"}   =  $hit->name;
 716         # only going to test order of magnitude
 717         # too hard as these don't always match
 718 #       $hash{"hit$hitcount\_signif"} =
 719 #           ( sprintf("%.0e",$hit->significance) =~ /e\-?(\d+)/ );
 720         $hash{"hit$hitcount\_bits"}   =  sprintf("%d",$hit->bits);
 721
 722         foreach my $hsp ( $hit->hsps ) {
 723             $hash{"hsp$hspcount\_bits"}   = sprintf("%d",$hsp->bits);
 724             # only going to test order of magnitude
 725             # too hard as these don't always match
 726 #           $hash{"hsp$hspcount\_evalue"} =
 727 #               ( sprintf("%.0e",$hsp->evalue) =~ /e\-?(\d+)/ );
 728             $hash{"hsp$hspcount\_qs"}     = $hsp->query->start;
 729             $hash{"hsp$hspcount\_qe"}     = $hsp->query->end;
 730             $hash{"hsp$hspcount\_qstr"}   = $hsp->query->strand;
 731             $hash{"hsp$hspcount\_hs"}     = $hsp->hit->start;
 732             $hash{"hsp$hspcount\_he"}     = $hsp->hit->end;
 733             $hash{"hsp$hspcount\_hstr"}   = $hsp->hit->strand;
 734
 735             #$hash{"hsp$hspcount\_pid"}     = sprintf("%d",$hsp->percent_identity);
 736             #$hash{"hsp$hspcount\_fid"}     = sprintf("%.2f",$hsp->frac_identical);
 737             $hash{"hsp$hspcount\_gaps"}    = $hsp->gaps('total');
 738             $hspcount++;
 739         }
 740         $hitcount++;
 741     }
 742     return %hash;
 743 }
 744
 745 sub _warn_about_no_hsps {
 746     my $hit = shift;
 747     my $prev_func=(caller(1))[3];
 748     $hit->warn("There is no HSP data for hit '".$hit->name."'.\n".
 749                "You have called a method ($prev_func)\n".
 750                "that requires HSP data and there was no HSP data for this hit,\n".
 751                "most likely because it was absent from the BLAST report.\n".
 752                "Note that by default, BLAST lists alignments for the first 250 hits,\n".
 753                "but it lists descriptions for 500 hits. If this is the case,\n".
 754                "and you care about these hits, you should re-run BLAST using the\n".
 755                "-b option (or equivalent if not using blastall) to increase the number\n".
 756                "of alignments.\n"
 757               );
 758 }
 759
 760 1;