branch-1-6/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 I<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 # deep bug fix: set $self->{'_last_data'} to undef if $$data{Data} is
 408 # a valid slot, whose value is undef --
 409 # allows an undef to be propagated to object constructors and
 410 # handled there as desired; in particular, when Hsp_postive => -conserved
 411 # is not defined (in BLASTN, e.g.), the value of hsp's {CONSERVED} property is
 412 # set to the value of {IDENTICAL}.
 413 #/maj
 414 #   return unless ( defined $data->{'Data'} );
 415    return unless ( grep /Data/, keys %$data );
 416    if ( !defined $data->{'Data'} ) {
 417        $self->{'_last_data'} = undef;
 418        return;
 419    }
 420    if( $data->{'Data'} =~ /^\s+$/ ) {
 421        return unless $data->{'Name'} =~ /Hsp\_(midline|qseq|hseq)/;
 422    }
 423
 424    if( $self->in_element('hsp') &&
 425        $data->{'Name'} =~ /Hsp\_(qseq|hseq|midline)/ ) {
 426
 427        $self->{'_last_hspdata'}->{$data->{'Name'}} .= $data->{'Data'};
 428    }
 429
 430    $self->{'_last_data'} = $data->{'Data'};
 431 }
 432
 433 =head2 _mode
 434
 435  Title   : _mode
 436  Usage   : $obj->_mode($newval)
 437  Function:
 438  Example :
 439  Returns : value of _mode
 440  Args    : newvalue (optional)
 441
 442
 443 =cut
 444
 445 sub _mode{
 446     my ($self,$value) = @_;
 447     if( defined $value) {
 448         $self->{'_mode'} = $value;
 449     }
 450     return $self->{'_mode'};
 451 }
 452
 453 =head2 within_element
 454
 455  Title   : within_element
 456  Usage   : if( $eventgenerator->within_element($element) ) {}
 457  Function: Test if we are within a particular element
 458            This is different than 'in' because within can be tested
 459            for a whole block.
 460  Returns : boolean
 461  Args    : string element name
 462
 463
 464 =cut
 465
 466 sub within_element{
 467    my ($self,$name) = @_;
 468    return 0 if ( ! defined $name &&
 469                  ! defined  $self->{'_elements'} ||
 470                  scalar @{$self->{'_elements'}} == 0) ;
 471    foreach (  @{$self->{'_elements'}} ) {
 472        if( $_ eq $name  ) {
 473            return 1;
 474        }
 475    }
 476    return 0;
 477 }
 478
 479 =head2 in_element
 480
 481  Title   : in_element
 482  Usage   : if( $eventgenerator->in_element($element) ) {}
 483  Function: Test if we are in a particular element
 484            This is different than 'in' because within can be tested
 485            for a whole block.
 486  Returns : boolean
 487  Args    : string element name
 488
 489
 490 =cut
 491
 492 sub in_element{
 493    my ($self,$name) = @_;
 494    return 0 if ! defined $self->{'_elements'}->[0];
 495    return ( $self->{'_elements'}->[0] eq $name)
 496 }
 497
 498
 499 =head2 start_document
 500
 501  Title   : start_document
 502  Usage   : $eventgenerator->start_document
 503  Function: Handles a start document event
 504  Returns : none
 505  Args    : none
 506
 507
 508 =cut
 509
 510 sub start_document{
 511     my ($self) = @_;
 512     $self->{'_lasttype'} = '';
 513     $self->{'_values'} = {};
 514     $self->{'_result'}= undef;
 515     $self->{'_mode'} = '';
 516     $self->{'_elements'} = [];
 517 }
 518
 519
 520 =head2 end_document
 521
 522  Title   : end_document
 523  Usage   : $eventgenerator->end_document
 524  Function: Handles an end document event
 525  Returns : Bio::Search::Result::ResultI object
 526  Args    : none
 527
 528
 529 =cut
 530
 531 sub end_document{
 532    my ($self,@args) = @_;
 533    return $self->{'_result'};
 534 }
 535
 536 =head2 result_count
 537
 538  Title   : result_count
 539  Usage   : my $count = $searchio->result_count
 540  Function: Returns the number of results we have processed
 541  Returns : integer
 542  Args    : none
 543
 544
 545 =cut
 546
 547 sub result_count {
 548     my $self = shift;
 549     return $self->{'_result_count'};
 550 }
 551
 552 sub report_count { shift->result_count }
 553
 554
 555 =head2 program_name
 556
 557  Title   : program_name
 558  Usage   : $obj->program_name($newval)
 559  Function: Get/Set the program name
 560  Returns : value of program_name (a scalar)
 561  Args    : on set, new value (a scalar or undef, optional)
 562
 563
 564 =cut
 565
 566 sub program_name{
 567     my $self = shift;
 568
 569     $self->{'program_name'} = shift if @_;
 570     return $self->{'program_name'} || $DefaultProgramName;
 571 }
 572
 573
 574 =head2 _will_handle
 575
 576  Title   : _will_handle
 577  Usage   : Private method. For internal use only.
 578               if( $self->_will_handle($type) ) { ... }
 579  Function: Provides an optimized way to check whether or not an element of a
 580            given type is to be handled.
 581  Returns : Reference to EventHandler object if the element type is to be handled.
 582            undef if the element type is not to be handled.
 583  Args    : string containing type of element.
 584
 585 Optimizations:
 586
 587 =over 2
 588
 589 =item 1
 590
 591 Using the cached pointer to the EventHandler to minimize repeated
 592 lookups.
 593
 594 =item 2
 595
 596 Caching the will_handle status for each type that is encountered so
 597 that it only need be checked by calling
 598 handler-E<gt>will_handle($type) once.
 599
 600 =back
 601
 602 This does not lead to a major savings by itself (only 5-10%).  In
 603 combination with other optimizations, or for large parse jobs, the
 604 savings good be significant.
 605
 606 To test against the unoptimized version, remove the parentheses from
 607 around the third term in the ternary " ? : " operator and add two
 608 calls to $self-E<gt>_eventHandler().
 609
 610 =cut
 611
 612 sub _will_handle {
 613     my ($self,$type) = @_;
 614     my $handler = $self->{'_handler'};
 615     my $will_handle = defined($self->{'_will_handle_cache'}->{$type})
 616                              ? $self->{'_will_handle_cache'}->{$type}
 617                              : ($self->{'_will_handle_cache'}->{$type} =
 618                                $handler->will_handle($type));
 619
 620     return $will_handle ? $handler : undef;
 621 }
 622
 623 1;