Bio/SearchIO/blasttable.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::SearchIO::blasttable
   4 #
   5 # Cared for by Jason Stajich <jason-at-bioperl-dot-org>
   6 #
   7 # Copyright Jason Stajich
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::SearchIO::blasttable - Driver module for SearchIO for parsing NCBI -m 8/9 format
  16
  17 =head1 SYNOPSIS
  18
  19   # do not use this module directly
  20   use Bio::SearchIO;
  21   my $parser = Bio::SearchIO->new(-file   => $file,
  22                                  -format => 'blasttable');
  23
  24   while( my $result = $parser->next_result ) {
  25   }
  26
  27 =head1 DESCRIPTION
  28
  29 This module will support parsing NCBI -m 8 or -m 9 tabular output
  30 and WU-BLAST -mformat 2 or -mformat 3 tabular output.
  31
  32 =head1 FEEDBACK
  33
  34 =head2 Mailing Lists
  35
  36 User feedback is an integral part of the evolution of this and other
  37 Bioperl modules. Send your comments and suggestions preferably to
  38 the Bioperl mailing list.  Your participation is much appreciated.
  39
  40   bioperl-l@bioperl.org                  - General discussion
  41   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  42
  43 =head2 Reporting Bugs
  44
  45 Report bugs to the Bioperl bug tracking system to help us keep track
  46 of the bugs and their resolution. Bug reports can be submitted via
  47 the web:
  48
  49   http://bugzilla.open-bio.org/
  50
  51 =head1 AUTHOR - Jason Stajich
  52
  53 Email jason-at-bioperl-dot-org
  54
  55 =head1 APPENDIX
  56
  57 The rest of the documentation details each of the object methods.
  58 Internal methods are usually preceded with a _
  59
  60 =cut
  61
  62 # Let the code begin...
  63
  64
  65 package Bio::SearchIO::blasttable;
  66 use vars qw(%MAPPING %MODEMAP $DEFAULT_WRITER_CLASS $DefaultProgramName);
  67 use strict;
  68 use Bio::Search::Result::ResultFactory;
  69 use Bio::Search::Hit::HitFactory;
  70 use Bio::Search::HSP::HSPFactory;
  71
  72 $DefaultProgramName = 'BLASTN';
  73 $DEFAULT_WRITER_CLASS = 'Bio::Search::Writer::HitTableWriter';
  74
  75 # mapping of terms to Bioperl hash keys
  76 %MODEMAP = (
  77             'Result'             => 'result',
  78             'Hit'                => 'hit',
  79             'Hsp'                => 'hsp'
  80             );
  81
  82 %MAPPING = (
  83              'Hsp_bit-score'  => 'HSP-bits',
  84              'Hsp_score'      => 'HSP-score',
  85              'Hsp_evalue'     => 'HSP-evalue',
  86              'Hsp_query-from' => 'HSP-query_start',
  87              'Hsp_query-to'   => 'HSP-query_end',
  88              'Hsp_hit-from'   => 'HSP-hit_start',
  89              'Hsp_hit-to'     => 'HSP-hit_end',
  90              'Hsp_positive'   => 'HSP-conserved',
  91              'Hsp_identity'   => 'HSP-identical',
  92              'Hsp_mismatches' => 'HSP-mismatches',
  93              'Hsp_qgapblocks' => 'HSP-query_gapblocks',
  94              'Hsp_hgapblocks' => 'HSP-hit_gapblocks',
  95              'Hsp_gaps'       => 'HSP-hsp_gaps',
  96              'Hsp_hitgaps'    => 'HSP-hit_gaps',
  97              'Hsp_querygaps'  => 'HSP-query_gaps',
  98              'Hsp_align-len'  => 'HSP-hsp_length',
  99              'Hsp_query-frame'=> 'HSP-query_frame',
 100              'Hsp_hit-frame'  => 'HSP-hit_frame',
 101
 102              'Hit_id'        => 'HIT-name',
 103              'Hit_len'       => 'HIT-length',
 104              'Hit_accession' => 'HIT-accession',
 105              'Hit_def'       => 'HIT-description',
 106              'Hit_signif'    => 'HIT-significance',
 107              'Hit_score'     => 'HIT-score',
 108              'Hit_bits'      => 'HIT-bits',
 109
 110              'Result_program'  => 'RESULT-algorithm_name',
 111              'Result_version'  => 'RESULT-algorithm_version',
 112              'Result_query-def'=> 'RESULT-query_name',
 113              'Result_query-len'=> 'RESULT-query_length',
 114              'Result_query-acc'=> 'RESULT-query_accession',
 115              'Result_querydesc'=> 'RESULT-query_description',
 116              'Result_db'       => 'RESULT-database_name',
 117              'Result_db-len'   => 'RESULT-database_entries',
 118              'Result_db-let'   => 'RESULT-database_letters',
 119              );
 120
 121 use base qw(Bio::SearchIO);
 122
 123 =head2 new
 124
 125  Title   : new
 126  Usage   : my $obj = Bio::SearchIO::blasttable->new();
 127  Function: Builds a new Bio::SearchIO::blasttable object
 128  Returns : an instance of Bio::SearchIO::blasttable
 129  Args    :
 130
 131
 132 =cut
 133
 134 sub _initialize {
 135     my ($self,@args) = @_;
 136     $self->SUPER::_initialize(@args);
 137
 138     my ($pname) = $self->_rearrange([qw(PROGRAM_NAME)],
 139                                     @args);
 140     $self->program_name($pname || $DefaultProgramName);
 141     $self->_eventHandler->register_factory('result', Bio::Search::Result::ResultFactory->new(-type => 'Bio::Search::Result::GenericResult'));
 142     $self->_eventHandler->register_factory('hit', Bio::Search::Hit::HitFactory->new(-type => 'Bio::Search::Hit::GenericHit'));
 143     $self->_eventHandler->register_factory('hsp', Bio::Search::HSP::HSPFactory->new(-type => 'Bio::Search::HSP::GenericHSP'));
 144 }
 145
 146
 147 =head2 next_result
 148
 149  Title   : next_result
 150  Usage   : my $result = $parser->next_result
 151  Function: Parse the next result from the data stream
 152  Returns : L<Bio::Search::Result::ResultI>
 153  Args    : none
 154
 155
 156 =cut
 157
 158 sub next_result{
 159    my ($self) = @_;
 160    my ($lastquery,$lasthit);
 161    local $/ = "\n";
 162    local $_;
 163    my ($alg, $ver);
 164    while( defined ($_ = $self->_readline) ) {
 165           # WU-BLAST -mformat 3 only
 166           if(m{^#\s((?:\S+?)?BLAST[NPX])\s(\d+\.\d+.+\d{4}\])}) {
 167             ($alg, $ver) = ($1, $2);
 168                         # only one header for whole file with WU-BLAST
 169                         # so $alg and $ver won't get set properly for
 170                         # each result
 171                         $self->program_name($alg) if $alg;
 172                         $self->element({'Name' => 'Result_version',
 173                                                         'Data' => $ver}) if $ver;
 174             next;
 175           }
 176       # -m 9 only
 177       elsif(m{^#\s+((?:\S+?)?BLAST[NPX])\s+(.+)}) {
 178             ($alg, $ver) = ($1, $2);
 179             next;
 180        }
 181        next if /^\#/ || /^\s+$/;
 182
 183           my @fields = split;
 184           my ($qname,$hname, $percent_id, $hsp_len, $mismatches,$gapsm,
 185               $qstart,$qend,$hstart,$hend,$evalue,$bits);
 186           # WU-BLAST-specific
 187           my ($num_scores, $raw_score, $identities, $positives, $percent_pos,
 188               $qgap_blocks,$qgaps, $sgap_blocks, $sgaps, $qframe,
 189               $sframe);
 190           # NCBI -m8 and -m9
 191           if (@fields == 12) {
 192               ($qname,$hname, $percent_id, $hsp_len, $mismatches,$gapsm,
 193                $qstart,$qend,$hstart,$hend,$evalue,$bits) = @fields;
 194           }
 195           # WU-BLAST -mformat 2 and 3
 196           elsif ((@fields == 22) or (@fields == 24)) {
 197               ($qname,$hname,$evalue,$num_scores, $bits, $raw_score, $hsp_len,
 198                $identities, $positives,$mismatches, $percent_id, $percent_pos,
 199                $qgap_blocks, $qgaps, $sgap_blocks, $sgaps, $qframe, $qstart,
 200                $qend, $sframe, $hstart,$hend,) = @fields;
 201               # we need total gaps in the alignment
 202               $gapsm=$qgaps+$sgaps;
 203           }
 204           else {}
 205
 206        # Remember Jim's code is 0 based
 207        if( defined $lastquery &&
 208            $lastquery ne $qname ) {
 209            $self->end_element({'Name' => 'Hit'});
 210            $self->end_element({'Name' => 'Result'});
 211            $self->_pushback($_);
 212            return $self->end_document;
 213        } elsif( ! defined $lastquery ) {
 214            $self->{'_result_count'}++;
 215            $self->start_element({'Name' => 'Result'});
 216            $self->element({'Name' => 'Result_program',
 217                            'Data' => $alg || $self->program_name});
 218        $self->element({'Name' => 'Result_version',
 219                            'Data' => $ver}) if $ver;
 220            $self->element({'Name' => 'Result_query-def',
 221                            'Data' => $qname});
 222            $self->start_element({'Name' => 'Hit'});
 223            $self->element({'Name' => 'Hit_id',
 224                            'Data' => $hname});
 225            # we'll store the 1st hsp bits as the hit bits
 226            $self->element({'Name' => 'Hit_bits',
 227                            'Data' => $bits});
 228            # we'll store the 1st hsp value as the hit evalue
 229            $self->element({'Name' => 'Hit_signif',
 230                            'Data' => $evalue});
 231
 232        } elsif( $lasthit ne $hname ) {
 233            if( $self->in_element('hit') ) {
 234                $self->end_element({'Name' => 'Hit'});
 235            }
 236            $self->start_element({'Name' => 'Hit'});
 237            $self->element({'Name' => 'Hit_id',
 238                            'Data' => $hname});
 239            # we'll store the 1st hsp bits as the hit bits
 240            $self->element({'Name' => 'Hit_bits',
 241                            'Data' => $bits});
 242            # we'll store the 1st hsp value as the hit evalue
 243            $self->element({'Name' => 'Hit_signif',
 244                            'Data' => $evalue});
 245        }
 246        my $identical = $hsp_len - $mismatches - $gapsm;
 247        $self->start_element({'Name' => 'Hsp'});
 248        $self->element({'Name' => 'Hsp_evalue',
 249                        'Data' => $evalue});
 250        $self->element({'Name' => 'Hsp_bit-score',
 251                        'Data' => $bits});
 252        $self->element({'Name' => 'Hsp_identity',
 253                        'Data' => $identical});
 254        $self->element({'Name' => 'Hsp_positive',
 255                        'Data' => $identical});
 256        $self->element({'Name' => 'Hsp_gaps',
 257                        'Data' => $gapsm});
 258        $self->element({'Name' => 'Hsp_query-from',
 259                        'Data' => $qstart});
 260        $self->element({'Name' => 'Hsp_query-to',
 261                        'Data' => $qend});
 262
 263        $self->element({'Name' => 'Hsp_hit-from',
 264                        'Data' => $hstart });
 265        $self->element({'Name' => 'Hsp_hit-to',
 266                        'Data' => $hend });
 267        $self->element({'Name' => 'Hsp_align-len',
 268                        'Data' => $hsp_len});
 269        $self->end_element({'Name' => 'Hsp'});
 270        $lastquery = $qname;
 271        $lasthit   = $hname;
 272    }
 273    # fencepost
 274    if( defined $lasthit && defined $lastquery ) {
 275        if( $self->in_element('hit') ) {
 276            $self->end_element({'Name' => 'Hit'});
 277        }
 278        $self->end_element({'Name' => 'Result'});
 279        return $self->end_document;
 280    }
 281 }
 282
 283 =head2 start_element
 284
 285  Title   : start_element
 286  Usage   : $eventgenerator->start_element
 287  Function: Handles a start element event
 288  Returns : none
 289  Args    : hashref with at least 2 keys 'Data' and 'Name'
 290
 291
 292 =cut
 293
 294 sub start_element{
 295    my ($self,$data) = @_;
 296     # we currently don't care about attributes
 297     my $nm = $data->{'Name'};
 298    if( my $type = $MODEMAP{$nm} ) {
 299         $self->_mode($type);
 300         if( $self->_will_handle($type) ) {
 301             my $func = sprintf("start_%s",lc $type);
 302             $self->_eventHandler->$func($data->{'Attributes'});
 303         }
 304         unshift @{$self->{'_elements'}}, $type;
 305     }
 306     if($nm eq 'Result') {
 307         $self->{'_values'} = {};
 308         $self->{'_result'}= undef;
 309         $self->{'_mode'} = '';
 310     }
 311
 312 }
 313
 314 =head2 end_element
 315
 316  Title   : start_element
 317  Usage   : $eventgenerator->end_element
 318  Function: Handles an end element event
 319  Returns : none
 320  Args    : hashref with at least 2 keys 'Data' and 'Name'
 321
 322
 323 =cut
 324
 325 sub end_element {
 326     my ($self,$data) = @_;
 327     my $nm = $data->{'Name'};
 328     my $rc;
 329     # Hsp are sort of weird, in that they end when another
 330     # object begins so have to detect this in end_element for now
 331
 332     if( my $type = $MODEMAP{$nm} ) {
 333         if( $self->_will_handle($type) ) {
 334             my $func = sprintf("end_%s",lc $type);
 335             $rc = $self->_eventHandler->$func($self->{'_reporttype'},
 336                                               $self->{'_values'});
 337         }
 338         shift @{$self->{'_elements'}};
 339
 340     } elsif( $MAPPING{$nm} ) {
 341         if ( ref($MAPPING{$nm}) =~ /hash/i ) {
 342             my $key = (keys %{$MAPPING{$nm}})[0];
 343             $self->{'_values'}->{$key}->{$MAPPING{$nm}->{$key}} = $self->{'_last_data'};
 344         } else {
 345             $self->{'_values'}->{$MAPPING{$nm}} = $self->{'_last_data'};
 346         }
 347     } else {
 348         $self->warn( "unknown nm $nm ignoring\n");
 349     }
 350     $self->{'_last_data'} = ''; # remove read data if we are at
 351                                 # end of an element
 352     $self->{'_result'} = $rc if( $nm eq 'Result' );
 353     return $rc;
 354
 355 }
 356
 357 =head2 element
 358
 359  Title   : element
 360  Usage   : $eventhandler->element({'Name' => $name, 'Data' => $str});
 361  Function: Convience method that calls start_element, characters, end_element
 362  Returns : none
 363  Args    : Hash ref with the keys 'Name' and 'Data'
 364
 365
 366 =cut
 367
 368 sub element{
 369    my ($self,$data) = @_;
 370    $self->start_element($data);
 371    $self->characters($data);
 372    $self->end_element($data);
 373 }
 374
 375
 376 =head2 characters
 377
 378  Title   : characters
 379  Usage   : $eventgenerator->characters($str)
 380  Function: Send a character events
 381  Returns : none
 382  Args    : string
 383
 384
 385 =cut
 386
 387 sub characters{
 388    my ($self,$data) = @_;
 389
 390    return unless ( defined $data->{'Data'} );
 391    if( $data->{'Data'} =~ /^\s+$/ ) {
 392        return unless $data->{'Name'} =~ /Hsp\_(midline|qseq|hseq)/;
 393    }
 394
 395    if( $self->in_element('hsp') &&
 396        $data->{'Name'} =~ /Hsp\_(qseq|hseq|midline)/ ) {
 397
 398        $self->{'_last_hspdata'}->{$data->{'Name'}} .= $data->{'Data'};
 399    }
 400
 401    $self->{'_last_data'} = $data->{'Data'};
 402 }
 403
 404 =head2 _mode
 405
 406  Title   : _mode
 407  Usage   : $obj->_mode($newval)
 408  Function:
 409  Example :
 410  Returns : value of _mode
 411  Args    : newvalue (optional)
 412
 413
 414 =cut
 415
 416 sub _mode{
 417     my ($self,$value) = @_;
 418     if( defined $value) {
 419         $self->{'_mode'} = $value;
 420     }
 421     return $self->{'_mode'};
 422 }
 423
 424 =head2 within_element
 425
 426  Title   : within_element
 427  Usage   : if( $eventgenerator->within_element($element) ) {}
 428  Function: Test if we are within a particular element
 429            This is different than 'in' because within can be tested
 430            for a whole block.
 431  Returns : boolean
 432  Args    : string element name
 433
 434
 435 =cut
 436
 437 sub within_element{
 438    my ($self,$name) = @_;
 439    return 0 if ( ! defined $name &&
 440                  ! defined  $self->{'_elements'} ||
 441                  scalar @{$self->{'_elements'}} == 0) ;
 442    foreach (  @{$self->{'_elements'}} ) {
 443        if( $_ eq $name  ) {
 444            return 1;
 445        }
 446    }
 447    return 0;
 448 }
 449
 450 =head2 in_element
 451
 452  Title   : in_element
 453  Usage   : if( $eventgenerator->in_element($element) ) {}
 454  Function: Test if we are in a particular element
 455            This is different than 'in' because within can be tested
 456            for a whole block.
 457  Returns : boolean
 458  Args    : string element name
 459
 460
 461 =cut
 462
 463 sub in_element{
 464    my ($self,$name) = @_;
 465    return 0 if ! defined $self->{'_elements'}->[0];
 466    return ( $self->{'_elements'}->[0] eq $name)
 467 }
 468
 469
 470 =head2 start_document
 471
 472  Title   : start_document
 473  Usage   : $eventgenerator->start_document
 474  Function: Handles a start document event
 475  Returns : none
 476  Args    : none
 477
 478
 479 =cut
 480
 481 sub start_document{
 482     my ($self) = @_;
 483     $self->{'_lasttype'} = '';
 484     $self->{'_values'} = {};
 485     $self->{'_result'}= undef;
 486     $self->{'_mode'} = '';
 487     $self->{'_elements'} = [];
 488 }
 489
 490
 491 =head2 end_document
 492
 493  Title   : end_document
 494  Usage   : $eventgenerator->end_document
 495  Function: Handles an end document event
 496  Returns : Bio::Search::Result::ResultI object
 497  Args    : none
 498
 499
 500 =cut
 501
 502 sub end_document{
 503    my ($self,@args) = @_;
 504    return $self->{'_result'};
 505 }
 506
 507 =head2 result_count
 508
 509  Title   : result_count
 510  Usage   : my $count = $searchio->result_count
 511  Function: Returns the number of results we have processed
 512  Returns : integer
 513  Args    : none
 514
 515
 516 =cut
 517
 518 sub result_count {
 519     my $self = shift;
 520     return $self->{'_result_count'};
 521 }
 522
 523 sub report_count { shift->result_count }
 524
 525
 526 =head2 program_name
 527
 528  Title   : program_name
 529  Usage   : $obj->program_name($newval)
 530  Function: Get/Set the program name
 531  Returns : value of program_name (a scalar)
 532  Args    : on set, new value (a scalar or undef, optional)
 533
 534
 535 =cut
 536
 537 sub program_name{
 538     my $self = shift;
 539
 540     $self->{'program_name'} = shift if @_;
 541     return $self->{'program_name'} || $DefaultProgramName;
 542 }
 543
 544
 545 =head2 _will_handle
 546
 547  Title   : _will_handle
 548  Usage   : Private method. For internal use only.
 549               if( $self->_will_handle($type) ) { ... }
 550  Function: Provides an optimized way to check whether or not an element of a
 551            given type is to be handled.
 552  Returns : Reference to EventHandler object if the element type is to be handled.
 553            undef if the element type is not to be handled.
 554  Args    : string containing type of element.
 555
 556 Optimizations:
 557
 558 =over 2
 559
 560 =item 1
 561
 562 Using the cached pointer to the EventHandler to minimize repeated
 563 lookups.
 564
 565 =item 2
 566
 567 Caching the will_handle status for each type that is encountered so
 568 that it only need be checked by calling
 569 handler-E<gt>will_handle($type) once.
 570
 571 =back
 572
 573 This does not lead to a major savings by itself (only 5-10%).  In
 574 combination with other optimizations, or for large parse jobs, the
 575 savings good be significant.
 576
 577 To test against the unoptimized version, remove the parentheses from
 578 around the third term in the ternary " ? : " operator and add two
 579 calls to $self-E<gt>_eventHandler().
 580
 581 =cut
 582
 583 sub _will_handle {
 584     my ($self,$type) = @_;
 585     my $handler = $self->{'_handler'};
 586     my $will_handle = defined($self->{'_will_handle_cache'}->{$type})
 587                              ? $self->{'_will_handle_cache'}->{$type}
 588                              : ($self->{'_will_handle_cache'}->{$type} =
 589                                $handler->will_handle($type));
 590
 591     return $will_handle ? $handler : undef;
 592 }
 593
 594 1;