Bio/SearchIO/blasttable.pm

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