Bio/SearchIO/waba.pm

   1 #
   2 # BioPerl module for Bio::SearchIO::waba
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason@bioperl.org>
   7 #
   8 # Copyright Jason Stajich
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::SearchIO::waba - SearchIO parser for Jim Kent WABA program
  17 alignment output
  18
  19 =head1 SYNOPSIS
  20
  21     # do not use this object directly, rather through Bio::SearchIO
  22
  23     use Bio::SearchIO;
  24     my $in = Bio::SearchIO->new(-format => 'waba',
  25                                -file   => 'output.wab');
  26     while( my $result = $in->next_result ) {
  27         while( my $hit = $result->next_hit ) {
  28             while( my $hsp = $result->next_hsp ) {
  29
  30             }
  31         }
  32     }
  33
  34 =head1 DESCRIPTION
  35
  36 This parser will process the waba output (NOT the human readable format).
  37
  38 =head1 FEEDBACK
  39
  40 =head2 Mailing Lists
  41
  42 User feedback is an integral part of the evolution of this and other
  43 Bioperl modules. Send your comments and suggestions preferably to
  44 the Bioperl mailing list.  Your participation is much appreciated.
  45
  46   bioperl-l@bioperl.org                  - General discussion
  47   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  48
  49 =head2 Support
  50
  51 Please direct usage questions or support issues to the mailing list:
  52
  53 I<bioperl-l@bioperl.org>
  54
  55 rather than to the module maintainer directly. Many experienced and
  56 reponsive experts will be able look at the problem and quickly
  57 address it. Please include a thorough description of the problem
  58 with code and data examples if at all possible.
  59
  60 =head2 Reporting Bugs
  61
  62 Report bugs to the Bioperl bug tracking system to help us keep track
  63 of the bugs and their resolution. Bug reports can be submitted via the
  64 web:
  65
  66   https://github.com/bioperl/bioperl-live/issues
  67
  68 =head1 AUTHOR - Jason Stajich
  69
  70 Email jason-at-bioperl.org
  71
  72 =head1 APPENDIX
  73
  74 The rest of the documentation details each of the object methods.
  75 Internal methods are usually preceded with a _
  76
  77 =cut
  78
  79
  80 # Let the code begin...
  81
  82
  83 package Bio::SearchIO::waba;
  84 use vars qw(%MODEMAP %MAPPING @STATES);
  85 use strict;
  86
  87 # Object preamble - inherits from Bio::Root::Root
  88
  89 use Bio::Search::Result::ResultFactory;
  90 use Bio::Search::HSP::HSPFactory;
  91
  92 use POSIX;
  93
  94 BEGIN {
  95     # mapping of NCBI Blast terms to Bioperl hash keys
  96     %MODEMAP = ('WABAOutput' => 'result',
  97                 'Hit'         => 'hit',
  98                 'Hsp'         => 'hsp'
  99                 );
 100     @STATES = qw(Hsp_qseq Hsp_hseq Hsp_stateseq);
 101     %MAPPING =
 102         (
 103           'Hsp_query-from'=>  'HSP-query_start',
 104           'Hsp_query-to'  =>  'HSP-query_end',
 105           'Hsp_hit-from'  =>  'HSP-hit_start',
 106           'Hsp_hit-to'    =>  'HSP-hit_end',
 107           'Hsp_qseq'      =>  'HSP-query_seq',
 108           'Hsp_hseq'      =>  'HSP-hit_seq',
 109           'Hsp_midline'   =>  'HSP-homology_seq',
 110           'Hsp_stateseq'  =>  'HSP-hmmstate_seq',
 111           'Hsp_align-len' =>  'HSP-hsp_length',
 112
 113           'Hit_id'        => 'HIT-name',
 114           'Hit_accession' => 'HIT-accession',
 115
 116           'WABAOutput_program'  => 'RESULT-algorithm_name',
 117           'WABAOutput_version'  => 'RESULT-algorithm_version',
 118           'WABAOutput_query-def'=> 'RESULT-query_name',
 119           'WABAOutput_query-db' => 'RESULT-query_database',
 120           'WABAOutput_db'       => 'RESULT-database_name',
 121           );
 122 }
 123
 124
 125 use base qw(Bio::SearchIO);
 126
 127 =head2 new
 128
 129  Title   : new
 130  Usage   : my $obj = Bio::SearchIO::waba->new();
 131  Function: Builds a new Bio::SearchIO::waba object
 132  Returns : Bio::SearchIO::waba
 133  Args    : see Bio::SearchIO
 134
 135 =cut
 136
 137 sub _initialize {
 138     my ($self,@args) = @_;
 139     $self->SUPER::_initialize(@args);
 140     $self->_eventHandler->register_factory('result', Bio::Search::Result::ResultFactory->new(-type => 'Bio::Search::Result::WABAResult'));
 141
 142     $self->_eventHandler->register_factory('hsp', Bio::Search::HSP::HSPFactory->new(-type => 'Bio::Search::HSP::WABAHSP'));
 143 }
 144
 145
 146 =head2 next_result
 147
 148  Title   : next_result
 149  Usage   : my $hit = $searchio->next_result;
 150  Function: Returns the next Result from a search
 151  Returns : Bio::Search::Result::ResultI object
 152  Args    : none
 153
 154 =cut
 155
 156 sub next_result{
 157     my ($self) = @_;
 158     local $/ = "\n";
 159     local $_;
 160
 161     my ($curquery,$curhit);
 162     my $state = -1;
 163     $self->start_document();
 164     my @hit_signifs;
 165     while( defined ($_ = $self->_readline )) {
 166
 167         if( $state == -1 ) {
 168             my ($qid, $qhspid,$qpercent, $junk,
 169                 $alnlen,$qdb,$qacc,$qstart,$qend,$qstrand,
 170                 $hitdb,$hacc,$hstart,$hend,
 171                 $hstrand) =
 172                     ( /^(\S+)\.(\S+)\s+align\s+ # get the queryid
 173                       (\d+(\.\d+)?)\%\s+     # get the percentage
 174                       of\s+(\d+)\s+  # get the length of the alignment
 175                       (\S+)\s+           # this is the query database
 176                       (\S+):(\-?\d+)\-(\-?\d+) # The accession:start-end for query
 177                       \s+([\-\+])        # query strand
 178                       \s+(\S+)\.         # hit db
 179                       (\S+):(\-?\d+)\-(\-?\d+) # The accession:start-end for hit
 180                       \s+([\-\+])\s*$    # hit strand
 181                       /ox );
 182
 183             # Curses.  Jim's code is 0 based, the following is to readjust
 184             if( $hstart < 0 ) { $hstart *= -1}
 185             if( $hend   < 0 ) { $hend   *= -1}
 186             if( $qstart < 0 ) { $qstart *= -1}
 187             if( $qend   < 0 ) { $qend   *= -1}
 188             $hstart++; $hend++; $qstart++; $qend++;
 189             if( ! defined $alnlen ) {
 190                 $self->warn("Unable to parse the rest of the WABA alignment info for: '$_'");
 191                 last;
 192             }
 193             $self->{'_reporttype'} = 'WABA'; # hardcoded - only
 194                                              # one type of WABA AFAIK
 195             if( defined $curquery &&
 196                 $curquery ne $qid ) {
 197                 $self->end_element({'Name' => 'Hit'});
 198                 $self->_pushback($_);
 199                 $self->end_element({'Name' => 'WABAOutput'});
 200                 return $self->end_document();
 201             }
 202
 203             if( defined $curhit &&
 204                 $curhit ne $hacc) {
 205                 # slight duplication here -- keep these in SYNC
 206                 $self->end_element({'Name' => 'Hit'});
 207                 $self->start_element({'Name' => 'Hit'});
 208                 $self->element({'Name' => 'Hit_id',
 209                                 'Data' => $hacc});
 210                 $self->element({'Name' => 'Hit_accession',
 211                                 'Data' => $hacc});
 212
 213             } elsif ( ! defined $curquery ) {
 214                 $self->start_element({'Name' => 'WABAOutput'});
 215                 $self->{'_result_count'}++;
 216                 $self->element({'Name' => 'WABAOutput_query-def',
 217                                 'Data' => $qid });
 218                 $self->element({'Name' => 'WABAOutput_program',
 219                                 'Data' => 'WABA'});
 220                 $self->element({'Name' => 'WABAOutput_query-db',
 221                                 'Data' => $qdb});
 222                 $self->element({'Name' => 'WABAOutput_db',
 223                                 'Data' => $hitdb});
 224
 225                 # slight duplication here -- keep these N'SYNC ;-)
 226                 $self->start_element({'Name' => 'Hit'});
 227                 $self->element({'Name' => 'Hit_id',
 228                                 'Data' => $hacc});
 229                 $self->element({'Name' => 'Hit_accession',
 230                                 'Data' => $hacc});
 231             }
 232
 233
 234             # strand is inferred by start,end values
 235             # in the Result Builder
 236             if( $qstrand eq '-' ) {
 237                 ($qstart,$qend) = ($qend,$qstart);
 238             }
 239             if( $hstrand eq '-' ) {
 240                 ($hstart,$hend) = ($hend,$hstart);
 241             }
 242
 243             $self->start_element({'Name' => 'Hsp'});
 244             $self->element({'Name' => 'Hsp_query-from',
 245                             'Data' => $qstart});
 246             $self->element({'Name' => 'Hsp_query-to',
 247                             'Data' => $qend});
 248             $self->element({'Name' => 'Hsp_hit-from',
 249                             'Data' => $hstart});
 250             $self->element({'Name' => 'Hsp_hit-to',
 251                             'Data' => $hend});
 252             $self->element({'Name' => 'Hsp_align-len',
 253                             'Data' => $alnlen});
 254
 255             $curquery = $qid;
 256             $curhit   = $hacc;
 257             $state = 0;
 258         } elsif( ! defined $curquery ) {
 259             $self->warn("skipping because no Hit begin line was recognized\n$_") if( $_ !~ /^\s+$/ );
 260             next;
 261         } else {
 262             chomp;
 263             $self->element({'Name' => $STATES[$state++],
 264                             'Data' => $_});
 265             if( $state >= scalar @STATES ) {
 266                 $state = -1;
 267                 $self->end_element({'Name' => 'Hsp'});
 268             }
 269         }
 270     }
 271     if( defined $curquery  ) {
 272         $self->end_element({'Name' => 'Hit'});
 273         $self->end_element({'Name' => 'WABAOutput'});
 274         return $self->end_document();
 275     }
 276     return;
 277 }
 278
 279 =head2 start_element
 280
 281  Title   : start_element
 282  Usage   : $eventgenerator->start_element
 283  Function: Handles a start element event
 284  Returns : none
 285  Args    : hashref with at least 2 keys 'Data' and 'Name'
 286
 287
 288 =cut
 289
 290 sub start_element{
 291    my ($self,$data) = @_;
 292     # we currently don't care about attributes
 293     my $nm = $data->{'Name'};
 294    if( my $type = $MODEMAP{$nm} ) {
 295         $self->_mode($type);
 296         if( $self->_eventHandler->will_handle($type) ) {
 297             my $func = sprintf("start_%s",lc $type);
 298             $self->_eventHandler->$func($data->{'Attributes'});
 299         }
 300         unshift @{$self->{'_elements'}}, $type;
 301     }
 302     if($nm eq 'WABAOutput') {
 303         $self->{'_values'} = {};
 304         $self->{'_result'}= undef;
 305         $self->{'_mode'} = '';
 306     }
 307
 308 }
 309
 310 =head2 end_element
 311
 312  Title   : start_element
 313  Usage   : $eventgenerator->end_element
 314  Function: Handles an end element event
 315  Returns : none
 316  Args    : hashref with at least 2 keys 'Data' and 'Name'
 317
 318
 319 =cut
 320
 321 sub end_element {
 322     my ($self,$data) = @_;
 323     my $nm = $data->{'Name'};
 324     my $rc;
 325     # Hsp are sort of weird, in that they end when another
 326     # object begins so have to detect this in end_element for now
 327     if( $nm eq 'Hsp' ) {
 328         foreach ( qw(Hsp_qseq Hsp_midline Hsp_hseq) ) {
 329             $self->element({'Name' => $_,
 330                             'Data' => $self->{'_last_hspdata'}->{$_}});
 331         }
 332         $self->{'_last_hspdata'} = {}
 333     }
 334
 335     if( my $type = $MODEMAP{$nm} ) {
 336         if( $self->_eventHandler->will_handle($type) ) {
 337             my $func = sprintf("end_%s",lc $type);
 338             $rc = $self->_eventHandler->$func($self->{'_reporttype'},
 339                                               $self->{'_values'});
 340         }
 341         shift @{$self->{'_elements'}};
 342
 343     } elsif( $MAPPING{$nm} ) {
 344         if ( ref($MAPPING{$nm}) =~ /hash/i ) {
 345             my $key = (keys %{$MAPPING{$nm}})[0];
 346             $self->{'_values'}->{$key}->{$MAPPING{$nm}->{$key}} = $self->{'_last_data'};
 347         } else {
 348             $self->{'_values'}->{$MAPPING{$nm}} = $self->{'_last_data'};
 349         }
 350     } else {
 351         $self->warn( "unknown nm $nm ignoring\n");
 352     }
 353     $self->{'_last_data'} = ''; # remove read data if we are at
 354                                 # end of an element
 355     $self->{'_result'} = $rc if( $nm eq 'WABAOutput' );
 356     return $rc;
 357
 358 }
 359
 360 =head2 element
 361
 362  Title   : element
 363  Usage   : $eventhandler->element({'Name' => $name, 'Data' => $str});
 364  Function: Convience method that calls start_element, characters, end_element
 365  Returns : none
 366  Args    : Hash ref with the keys 'Name' and 'Data'
 367
 368
 369 =cut
 370
 371 sub element{
 372    my ($self,$data) = @_;
 373    $self->start_element($data);
 374    $self->characters($data);
 375    $self->end_element($data);
 376 }
 377
 378
 379 =head2 characters
 380
 381  Title   : characters
 382  Usage   : $eventgenerator->characters($str)
 383  Function: Send a character events
 384  Returns : none
 385  Args    : string
 386
 387
 388 =cut
 389
 390 sub characters{
 391    my ($self,$data) = @_;
 392
 393    return unless ( defined $data->{'Data'} );
 394    if( $data->{'Data'} =~ /^\s+$/ ) {
 395        return unless $data->{'Name'} =~ /Hsp\_(midline|qseq|hseq)/;
 396    }
 397
 398    if( $self->in_element('hsp') &&
 399        $data->{'Name'} =~ /Hsp\_(qseq|hseq|midline)/ ) {
 400
 401        $self->{'_last_hspdata'}->{$data->{'Name'}} .= $data->{'Data'};
 402    }
 403
 404    $self->{'_last_data'} = $data->{'Data'};
 405 }
 406
 407 =head2 _mode
 408
 409  Title   : _mode
 410  Usage   : $obj->_mode($newval)
 411  Function:
 412  Example :
 413  Returns : value of _mode
 414  Args    : newvalue (optional)
 415
 416
 417 =cut
 418
 419 sub _mode{
 420     my ($self,$value) = @_;
 421     if( defined $value) {
 422         $self->{'_mode'} = $value;
 423     }
 424     return $self->{'_mode'};
 425 }
 426
 427 =head2 within_element
 428
 429  Title   : within_element
 430  Usage   : if( $eventgenerator->within_element($element) ) {}
 431  Function: Test if we are within a particular element
 432            This is different than 'in' because within can be tested
 433            for a whole block.
 434  Returns : boolean
 435  Args    : string element name
 436
 437
 438 =cut
 439
 440 sub within_element{
 441    my ($self,$name) = @_;
 442    return 0 if ( ! defined $name &&
 443                  ! defined  $self->{'_elements'} ||
 444                  scalar @{$self->{'_elements'}} == 0) ;
 445    foreach (  @{$self->{'_elements'}} ) {
 446        if( $_ eq $name  ) {
 447            return 1;
 448        }
 449    }
 450    return 0;
 451 }
 452
 453 =head2 in_element
 454
 455  Title   : in_element
 456  Usage   : if( $eventgenerator->in_element($element) ) {}
 457  Function: Test if we are in 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 in_element{
 467    my ($self,$name) = @_;
 468    return 0 if ! defined $self->{'_elements'}->[0];
 469    return ( $self->{'_elements'}->[0] eq $name)
 470 }
 471
 472
 473 =head2 start_document
 474
 475  Title   : start_document
 476  Usage   : $eventgenerator->start_document
 477  Function: Handles a start document event
 478  Returns : none
 479  Args    : none
 480
 481
 482 =cut
 483
 484 sub start_document{
 485     my ($self) = @_;
 486     $self->{'_lasttype'} = '';
 487     $self->{'_values'} = {};
 488     $self->{'_result'}= undef;
 489     $self->{'_mode'} = '';
 490     $self->{'_elements'} = [];
 491 }
 492
 493
 494 =head2 end_document
 495
 496  Title   : end_document
 497  Usage   : $eventgenerator->end_document
 498  Function: Handles an end document event
 499  Returns : Bio::Search::Result::ResultI object
 500  Args    : none
 501
 502
 503 =cut
 504
 505 sub end_document{
 506    my ($self,@args) = @_;
 507    return $self->{'_result'};
 508 }
 509
 510 =head2 result_count
 511
 512  Title   : result_count
 513  Usage   : my $count = $searchio->result_count
 514  Function: Returns the number of results we have processed
 515  Returns : integer
 516  Args    : none
 517
 518
 519 =cut
 520
 521 sub result_count {
 522     my $self = shift;
 523     return $self->{'_result_count'};
 524 }
 525
 526 sub report_count { shift->result_count }
 527
 528 1;