← Back to branch summary

~ubuntu-branches/ubuntu/saucy/bioperl/saucy-proposed

« back to all changes in this revision

Viewing changes to Bio/SearchIO/blast.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2009-03-10 07:19:11 UTC
  • mfrom: (1.2.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20090310071911-fukqzw54pyb1f0bd
Tags: 1.6.0-2
https://launchpad.net/bugs/324001
* Removed patch system (not used):
  - removed instuctions in debian/rules;
  - removed quilt from Build-Depends in debian/control.
* Re-enabled tests:
  - uncommented test command in debian/rules;
  - uncommented previously missing build-dependencies in debian/control.
  - Re-enabled tests and uncommented build-dependencies accordingly.
* Removed libmodule-build-perl and libtest-harness-perl from
  Build-Depends-Indep (provided by perl-modules).
* Better cleaning of empty directories using find -type d -empty -delete
  instead of rmdir in debian/rules (LP: #324001).

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/SearchIO/blast.pm
1
 
# $Id: blast.pm,v 1.99.4.4 2006/10/11 19:55:17 sendu Exp $
 
1
# $Id: blast.pm 15383 2009-01-16 16:45:08Z cjfields $
2
2
#
3
3
# BioPerl module for Bio::SearchIO::blast
4
4
#
21
21
#          Parse more blast statistics, lambda, entropy, etc
22
22
#          from WU-BLAST in frame-specific manner
23
23
# 20060216 - cjf - fixed blast parsing for BLAST v2.2.13 output
 
24
# 20071104 - dmessina - added support for WUBLAST -echofilter
 
25
# 20071121 - cjf - fixed several bugs (bugs 2391, 2399, 2409)
24
26
 
25
27
=head1 NAME
26
28
 
33
35
   # Bio::SearchIO system.
34
36
 
35
37
    use Bio::SearchIO;
36
 
    my $searchio = new Bio::SearchIO(-format => 'blast',
 
38
    my $searchio = Bio::SearchIO->new(-format => 'blast',
37
39
                                     -file   => 't/data/ecolitst.bls');
38
40
    while( my $result = $searchio->next_result ) {
39
41
        while( my $hit = $result->next_hit ) {
80
82
you can supply the program type with -report_type in the SearchIO
81
83
constructor i.e.
82
84
 
83
 
  my $parser = new Bio::SearchIO(-format => 'blast',
 
85
  my $parser = Bio::SearchIO->new(-format => 'blast',
84
86
                                 -file => 'bl2seq.tblastn.report',
85
87
                                 -report_type => 'tblastn');
86
88
 
139
141
 
140
142
 
141
143
use base qw(Bio::SearchIO);
 
144
use Data::Dumper;
142
145
 
143
146
BEGIN {
144
147
 
157
160
        'Hsp_bit-score'   => 'HSP-bits',
158
161
        'Hsp_score'       => 'HSP-score',
159
162
        'Hsp_evalue'      => 'HSP-evalue',
 
163
        'Hsp_n',          => 'HSP-n',
160
164
        'Hsp_pvalue'      => 'HSP-pvalue',
161
165
        'Hsp_query-from'  => 'HSP-query_start',
162
166
        'Hsp_query-to'    => 'HSP-query_end',
175
179
        'Hsp_hit-frame'   => 'HSP-hit_frame',
176
180
        'Hsp_links'       => 'HSP-links',
177
181
        'Hsp_group'       => 'HSP-hsp_group',
 
182
        'Hsp_features'    => 'HSP-hit_features',
178
183
 
179
184
        'Hit_id'        => 'HIT-name',
180
185
        'Hit_len'       => 'HIT-length',
181
186
        'Hit_accession' => 'HIT-accession',
182
187
        'Hit_def'       => 'HIT-description',
183
188
        'Hit_signif'    => 'HIT-significance',
184
 
 
185
189
        # For NCBI blast, the description line contains bits.
186
190
        # For WU-blast, the  description line contains score.
187
191
        'Hit_score' => 'HIT-score',
195
199
        'BlastOutput_query-def'           => 'RESULT-query_name',
196
200
        'BlastOutput_query-len'           => 'RESULT-query_length',
197
201
        'BlastOutput_query-acc'           => 'RESULT-query_accession',
 
202
        'BlastOutput_query-gi'            => 'RESULT-query_gi',
198
203
        'BlastOutput_querydesc'           => 'RESULT-query_description',
199
204
        'BlastOutput_db'                  => 'RESULT-database_name',
200
205
        'BlastOutput_db-len'              => 'RESULT-database_entries',
311
316
=head2 new
312
317
 
313
318
 Title   : new
314
 
 Usage   : my $obj = new Bio::SearchIO::blast(%args);
 
319
 Usage   : my $obj = Bio::SearchIO::blast->new(%args);
315
320
 Function: Builds a new Bio::SearchIO::blast object 
316
321
 Returns : Bio::SearchIO::blast
317
322
 Args    : Key-value pairs:
370
375
    my ( $self, @args ) = @_;
371
376
    $self->SUPER::_initialize(@args);
372
377
 
373
 
 # Blast reports require a specialized version of the SREB due to the
374
 
 # possibility of iterations (PSI-BLAST). Forwarding all arguments to it.
375
 
 # An issue here is that we want to set new default object factories if none are
376
 
 # supplied.
 
378
    # Blast reports require a specialized version of the SREB due to the
 
379
    # possibility of iterations (PSI-BLAST). Forwarding all arguments to it. An
 
380
    # issue here is that we want to set new default object factories if none are
 
381
    # supplied.
377
382
 
378
 
    my $handler = new Bio::SearchIO::IteratedSearchResultEventBuilder(@args);
 
383
    my $handler = Bio::SearchIO::IteratedSearchResultEventBuilder->new(@args);
379
384
    $self->attach_EventHandler($handler);
380
385
    
381
 
    # 2006-04-26 move this to the attach_handler function in this
382
 
    # module so we can really reset the handler 
 
386
    # 2006-04-26 move this to the attach_handler function in this module so we
 
387
    # can really reset the handler 
383
388
    # Optimization: caching
384
389
    # the EventHandler since it is used a lot during the parse.
 
390
    
385
391
    # $self->{'_handler_cache'} = $handler;
386
392
 
387
393
    my ( $min_qlen, $check_all, $overlap, $best, $rpttype ) = $self->_rearrange(
427
433
    my $data   = '';
428
434
    my $flavor = '';
429
435
    $self->{'_seentop'} = 0;     # start next report at top
 
436
    $self->{'_seentop'} = 0;
430
437
    my ( $reporttype, $seenquery, $reportline );  
431
438
    my ( $seeniteration, $found_again );
432
439
    my $incl_threshold = $self->inclusion_threshold;
436
443
    my $gapped_stats = 0;    # for switching between gapped/ungapped
437
444
                             # lambda, K, H
438
445
    local $_ = "\n";   #consistency
 
446
    PARSER:
439
447
    while ( defined( $_ = $self->_readline ) ) {
440
448
        next if (/^\s+$/);       # skip empty lines
441
449
        next if (/CPU time:/);
442
450
        next if (/^>\s*$/);
443
451
        if (
444
 
               /^([T]?BLAST[NPX])\s*(.+)$/i               # NCBI BLAST
445
 
            || /^(PSITBLASTN)\s+(.+)$/i                   # PSIBLAST
446
 
            || /^(RPS-BLAST)\s*(.+)$/i                    # RPSBLAST
447
 
            || /^(MEGABLAST)\s*(.+)$/i                    # MEGABLAST
 
452
               /^((?:\S+?)?BLAST[NPX]?)\s+(.+)$/i  # NCBI BLAST, PSIBLAST
 
453
                                                   # RPSBLAST, MEGABLAST
448
454
            || /^(P?GENEWISE|HFRAME|SWN|TSWN)\s+(.+)/i    #Paracel BTK
449
455
          )
450
456
        {
451
 
            $self->debug("blast.pm: Start of new report: $1 $2\n");
 
457
            ($reporttype, my $reportversion) = ($1, $2);
 
458
            # need to keep track of whether this is WU-BLAST
 
459
            if ($reportversion && $reportversion =~ m{WashU$}) {
 
460
                $self->{'_wublast'}++;
 
461
            }
 
462
            $self->debug("blast.pm: Start of new report: $reporttype, $reportversion\n");
452
463
            if ( $self->{'_seentop'} ) { 
453
464
                # This handles multi-result input streams
454
465
                $self->_pushback($_);
455
 
                $self->in_element('hsp')
456
 
                  && $self->end_element( { 'Name' => 'Hsp' } );
457
 
                $self->in_element('hit')
458
 
                  && $self->end_element( { 'Name' => 'Hit' } );
459
 
                $self->within_element('iteration')
460
 
                  && $self->end_element( { 'Name' => 'Iteration' } );
461
 
                $self->end_element( { 'Name' => 'BlastOutput' } );
462
 
                return $self->end_document();
 
466
                last PARSER;
463
467
            }
464
468
            $self->_start_blastoutput;
465
 
            $reporttype = $1;
466
469
            if ($reporttype =~ /RPS-BLAST/) {
467
470
                $reporttype .= '(BLASTP)'; # default RPS-BLAST type
468
471
            }
477
480
            $self->element(
478
481
                {
479
482
                    'Name' => 'BlastOutput_version',
480
 
                    'Data' => $2
 
483
                    'Data' => $reportversion
481
484
                }
482
485
            );
483
486
            $self->element(
499
502
            if ( defined $seeniteration ) {
500
503
                $self->within_element('iteration')
501
504
                  && $self->end_element( { 'Name' => 'Iteration' } );
502
 
                $self->_start_iteration;
 
505
                $self->start_element( { 'Name' => 'Iteration' } );
503
506
            }
504
507
            else {
505
 
                $self->_start_iteration;
 
508
                $self->start_element( { 'Name' => 'Iteration' } );
506
509
            }
507
510
            $seeniteration = 1;
508
511
        }
509
512
        elsif (/^Query=\s*(.*)$/) {
 
513
            my $q    = $1;
510
514
            $self->debug("blast.pm: Query= found...$_\n");
511
 
            my $q    = $1;
512
515
            my $size = 0;
513
516
            if ( defined $seenquery ) {
514
517
                $self->_pushback($reportline) if $reportline;
515
518
                $self->_pushback($_);
516
 
                $self->in_element('hsp')
517
 
                  && $self->end_element( { 'Name' => 'Hsp' } );
518
 
                $self->in_element('hit')
519
 
                  && $self->end_element( { 'Name' => 'Hit' } );
520
 
                $self->within_element('iteration')
521
 
                  && $self->end_element( { 'Name' => 'Iteration' } );
522
 
                if ($bl2seq_fix) {
523
 
                    $self->element(
524
 
                        {
525
 
                            'Name' => 'BlastOutput_program',
526
 
                            'Data' => $reporttype
527
 
                        }
528
 
                    );
529
 
                }
530
 
                $self->end_element( { 'Name' => 'BlastOutput' } );
531
 
                return $self->end_document();
 
519
                last PARSER;
532
520
            }
533
521
            else {
534
522
                if ( !defined $reporttype ) {
536
524
                    if ( defined $seeniteration ) {
537
525
                        $self->in_element('iteration')
538
526
                          && $self->end_element( { 'Name' => 'Iteration' } );
539
 
                        $self->_start_iteration;
 
527
                        $self->start_element( { 'Name' => 'Iteration' } );
540
528
                    }
541
529
                    else {
542
 
                        $self->_start_iteration;
 
530
                        $self->start_element( { 'Name' => 'Iteration' } );
543
531
                    }
544
532
                    $seeniteration = 1;
545
533
                }
551
539
                    $self->_pushback($_);
552
540
                    last;
553
541
                }
554
 
                chomp;
555
542
                # below line fixes length issue with BLAST v2.2.13; still works 
556
543
                # with BLAST v2.2.12
557
544
                if ( /\((\-?[\d,]+)\s+letters.*\)/ || /^Length=(\-?[\d,]+)/ ) {
560
547
                    last;
561
548
                }
562
549
                else {
563
 
                    $q .= " $_";
564
 
                    $q =~ s/ +/ /g;
 
550
                    # bug 2391
 
551
                    $q .= ($q =~ /\w$/ && $_ =~ /^\w/) ? " $_" : $_;
 
552
                    $q =~ s/\s+/ /g; # this catches the newline as well
565
553
                    $q =~ s/^ | $//g;
566
554
                }
567
555
 
574
562
                    'Name' => 'BlastOutput_query-def',
575
563
                    'Data' => $nm
576
564
                }
577
 
            );
 
565
            ) if $nm;
578
566
            $self->element(
579
567
                {
580
568
                    'Name' => 'BlastOutput_query-len',
588
576
                    'Data' => $desc
589
577
                }
590
578
            );
591
 
            my ( $acc, $version ) = &_get_accession_version($nm);
 
579
            my ( $gi, $acc, $version ) = $self->_get_seq_identifiers($nm);
592
580
            $version = defined($version) && length($version) ? ".$version" : "";
593
 
            $acc = '' unless defined($acc);
594
581
            $self->element(
595
582
                {
596
583
                    'Name' => 'BlastOutput_query-acc',
597
584
                    'Data' => "$acc$version"
598
585
                }
599
 
            );
 
586
            ) if $acc;
600
587
        }
 
588
                # added check for WU-BLAST -echofilter option (bug 2388)
 
589
                elsif (/^>Unfiltered[+-]1$/) {
 
590
                # skip all of the lines of unfiltered sequence
 
591
                        while($_ !~ /^Database:/) {
 
592
                                $self->debug("Bypassing features line: $_");
 
593
                        $_ = $self->_readline;
 
594
                }
 
595
                        $self->_pushback($_);
 
596
                }
601
597
        elsif (/Sequences producing significant alignments:/) {
602
598
            $self->debug("blast.pm: Processing NCBI-BLAST descriptions\n");
603
599
            $flavor = 'ncbi';
611
607
            # for a line with a leading >. Blank lines occur with this section
612
608
            # for psiblast.
613
609
            if ( !$self->in_element('iteration') ) {
614
 
                $self->_start_iteration;
 
610
                $self->start_element( { 'Name' => 'Iteration' } );
615
611
            }
616
 
          descline:
617
 
            while ( defined( $_ = $self->_readline() ) ) {
618
 
                if (/^>/ 
619
 
                    || /^\s+Database:\s+?/
620
 
                    || /^Parameters:/
621
 
                    || /^\s+Subset/
622
 
                    || /^\s*Lambda/
623
 
                    || /^\s*Histogram/
624
 
                    ) {
625
 
                    $self->_pushback($_); # Catch leading > (end of section)
626
 
                    last descline;
627
 
                }
628
 
                elsif (/([\d\.\+\-eE]+)\s+([\d\.\+\-eE]+)(\s+\d+)?\s*$/) {
629
 
 
630
 
                    # the last match is for gapped BLAST output
631
 
                    # which will report the number of HSPs for the Hit
632
 
                    my ( $score, $evalue ) = ( $1, $2 );
633
 
 
 
612
            # these elements are dropped with some multiquery reports; add
 
613
            # back here
 
614
            $self->element(
 
615
                {
 
616
                    'Name' => 'BlastOutput_db-len',
 
617
                    'Data' => $self->{'_blsdb_length'}
 
618
                } 
 
619
            ) if $self->{'_blsdb_length'};
 
620
            $self->element(
 
621
                {
 
622
                    'Name' => 'BlastOutput_db-let',
 
623
                    'Data' => $self->{'_blsdb_letters'}
 
624
                }
 
625
            ) if $self->{'_blsdb_letters'};
 
626
            $self->element(
 
627
                {
 
628
                    'Name' => 'BlastOutput_db',
 
629
                    'Data' => $self->{'_blsdb'}
 
630
                }
 
631
            ) if $self->{'_blsdb_letters'};
 
632
            
 
633
            # changed 8/28/2008 to exit hit table if blank line is found after an
 
634
            # appropriate line
 
635
            my $h_regex;
 
636
            my $seen_block;
 
637
          DESCLINE:
 
638
            while ( defined( my $descline = $self->_readline() ) ) {
 
639
                if ($descline =~ m{^\s*$}) {
 
640
                    last DESCLINE if $seen_block;
 
641
                    next DESCLINE;
 
642
                }
 
643
                # any text match is part of block...
 
644
                $seen_block++;
 
645
                # GCG multiline oddness...
 
646
                if ($descline =~ /^(\S+)\s+Begin:\s\d+\s+End:\s+\d+/xms) {
 
647
                    my ($id, $nextline) = ($1, $self->_readline);
 
648
                    $nextline =~ s{^!}{};
 
649
                    $descline = "$id $nextline";
 
650
                }
 
651
                # NCBI style hit table (no N)
 
652
                if ($descline =~ /(?<!cor)          # negative lookahead
 
653
                    (\d*\.?(?:[\+\-eE]+)?\d+)       # number (float or scientific notation)
 
654
                    \s+                             # space
 
655
                    (\d*\.?(?:[\+\-eE]+)?\d+)       # number (float or scientific notation)
 
656
                    \s*$/xms) {
 
657
                    
 
658
                    my ( $score, $evalue ) = ($1, $2);
 
659
                    
634
660
                    # Some data clean-up so e-value will appear numeric to perl
635
661
                    $evalue =~ s/^e/1e/i;
636
662
 
637
663
                    # This to handle no-HSP case
638
 
                    my @line = split;
 
664
                    my @line = split ' ',$descline;
639
665
                    
640
666
                    # we want to throw away the score, evalue
641
667
                    pop @line, pop @line;
646
672
 
647
673
                    # add the last 2 entries s.t. we can reconstruct
648
674
                    # a minimal Hit object at the end of the day
649
 
                    push @hit_signifs,
650
 
                      [ $evalue, $score, shift @line, join( ' ', @line ) ];
651
 
                    
652
 
                }
653
 
                elsif (/^CONVERGED/i) {
 
675
                    push @hit_signifs, [ $evalue, $score, shift @line, join( ' ', @line ) ];
 
676
                } elsif ($descline =~ /^CONVERGED/i) {
654
677
                    $self->element(
655
678
                        {
656
679
                            'Name' => 'Iteration_converged',
657
680
                            'Data' => 1
658
681
                        }
659
682
                    );
 
683
                } else {
 
684
                    $self->_pushback($descline); # Catch leading > (end of section)
 
685
                    last DESCLINE;
660
686
                }
661
 
                @hit_signifs = sort {$a->[0] <=> $b->[0]} @hit_signifs;
662
687
            }
663
688
        }
664
689
        elsif (/Sequences producing High-scoring Segment Pairs:/) {
670
695
            $flavor = 'wu';
671
696
 
672
697
            if ( !$self->in_element('iteration') ) {
673
 
                $self->_start_iteration;
 
698
                $self->start_element( { 'Name' => 'Iteration' } );
674
699
            }
675
700
 
676
701
            while ( defined( $_ = $self->_readline() )
712
737
                            'Data' => $l
713
738
                        }
714
739
                    );
 
740
                    # cache for next round in cases with multiple queries
 
741
                    $self->{'_blsdb'} = $db;
 
742
                    $self->{'_blsdb_length'} = $s;
 
743
                    $self->{'_blsdb_letters'} = $l;                    
715
744
                    last;
716
745
                }
717
746
                else {
728
757
        }
729
758
        # bypasses this NCBI blast 2.2.13 extra output for now...
730
759
                # Features in/flanking this part of subject sequence:
731
 
        elsif (/^\sFeatures\s\w+\sthis\spart\sof\ssubject\ssequence:/) {
732
 
                # junk following lines up to start of HSP
733
 
                        while($_ !~ /^\sScore\s=/) {
734
 
                                $self->debug("Bypassing features line: $_");
 
760
        elsif (/^\sFeatures\s\w+\sthis\spart/xmso) {
 
761
            my $featline;
 
762
            $_ = $self->_readline;
 
763
                        while($_ !~ /^\s*$/) {
 
764
                chomp;
 
765
                $featline .= $_;
735
766
                        $_ = $self->_readline;
736
767
                }
737
768
                        $self->_pushback($_);
 
769
            $featline =~ s{(?:^\s+|\s+^)}{}g;
 
770
            $featline =~ s{\n}{;}g;
 
771
            $self->{'_last_hspdata'}->{'Hsp_features'} = $featline;
738
772
        }
739
773
        
740
774
        # move inside of a hit
751
785
            # Query=
752
786
            if ( !$self->within_element('result') ) {
753
787
                $self->_start_blastoutput;
754
 
                $self->_start_iteration;
 
788
                $self->start_element( { 'Name' => 'Iteration' } );
755
789
            }
756
790
            elsif ( !$self->within_element('iteration') ) {
757
 
                $self->_start_iteration;
 
791
                $self->start_element( { 'Name' => 'Iteration' } );
758
792
            }
759
793
            $self->start_element( { 'Name' => 'Hit' } );
760
794
            my $id         = $1;
767
801
                    'Data' => $id
768
802
                }
769
803
            );
770
 
            my ( $acc, $version ) = &_get_accession_version($id);
 
804
            my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
771
805
            $self->element(
772
806
                {
773
807
                    'Name' => 'Hit_accession',
774
808
                    'Data' => $acc
775
809
                }
776
810
            );
777
 
 
778
811
            # add hit significance (from the hit table)
779
 
            # this is where Bug 1986 goes awry
780
 
            
781
 
            my $v = shift @hit_signifs;
782
 
            if ( defined $v ) {
783
 
                $self->element(
784
 
                    {
785
 
                        'Name' => 'Hit_signif',
786
 
                        'Data' => $v->[0]
787
 
                    }
788
 
                );
789
 
                $self->element(
790
 
                    {
791
 
                        'Name' => 'Hit_score',
792
 
                        'Data' => $v->[1]
793
 
                    }
794
 
                );
 
812
            # this is where Bug 1986 went awry
 
813
            
 
814
            # Changed for Bug2409; hit->significance and hit->score/bits derived
 
815
            # from HSPs, not hit table unless necessary
 
816
            
 
817
            HITTABLE:
 
818
            while (my $v = shift @hit_signifs) {        
 
819
                my $tableid = $v->[2];
 
820
                if ($tableid !~ m{\Q$id\E}) {
 
821
                    $self->debug("Hit table ID $tableid doesn't match current hit id $id, checking next hit table entry...\n");
 
822
                    next HITTABLE;
 
823
                } else {
 
824
                    last HITTABLE;
 
825
                }
795
826
            }
796
827
            while ( defined( $_ = $self->_readline() ) ) {
797
828
                next if (/^\s+$/);
808
839
                    last;
809
840
                }
810
841
                else {
 
842
                    s/^\s(?!\s)/\x01/; #new line to concatenate desc lines with <soh>
811
843
                    $restofline .= $_;
812
844
                }
813
845
            }
957
989
        }
958
990
        elsif (
959
991
            ( $self->in_element('hit') || $self->in_element('hsp') )
960
 
            &&    # ncbi blast
 
992
            &&    # ncbi blast, works with 2.2.17
961
993
            m/Score\s*=\s*(\S+)\s*bits\s* # Bit score
962
994
                (?:\((\d+)\))?,            # Missing for BLAT pseudo-BLAST fmt 
963
 
                \s*Expect(?:\(\d+\+?\))?\s*=\s*(\S+) # E-value
 
995
                \s*Expect(?:\((\d+\+?)\))?\s*=\s*([^,\s]+) # E-value
964
996
                /ox
965
997
          )
966
998
        {         # parse NCBI blast HSP
968
1000
              && $self->end_element( { 'Name' => 'Hsp' } );
969
1001
 
970
1002
            # Some data clean-up so e-value will appear numeric to perl
971
 
            my ( $bits, $score, $evalue ) = ( $1, $2, $3 );
 
1003
            my ( $bits, $score, $n, $evalue ) = ( $1, $2, $3, $4 );
972
1004
            $evalue =~ s/^e/1e/i;
973
 
 
974
1005
            $self->start_element( { 'Name' => 'Hsp' } );
975
1006
            $self->element(
976
1007
                {
990
1021
                    'Data' => $evalue
991
1022
                }
992
1023
            );
 
1024
            $self->element(
 
1025
                {
 
1026
                    'Name' => 'Hsp_n',
 
1027
                    'Data' => $n
 
1028
                }
 
1029
            ) if defined $n;
993
1030
            $score = '' unless defined $score;    # deal with BLAT which
994
1031
                                                  # has no score only bits
995
1032
            $self->debug("Got NCBI HSP score=$score, evalue $evalue\n");
1140
1177
 
1141
1178
            # This is for the case when we specify -b 0 (or B=0 for WU-BLAST)
1142
1179
            # and still want to construct minimal Hit objects
1143
 
            while ( my $v = shift @hit_signifs ) {
1144
 
                next unless defined $v;
1145
 
                $self->start_element( { 'Name' => 'Hit' } );
1146
 
                my $id   = $v->[2];
1147
 
                my $desc = $v->[3];
1148
 
                $self->element(
1149
 
                    {
1150
 
                        'Name' => 'Hit_id',
1151
 
                        'Data' => $id
1152
 
                    }
1153
 
                );
1154
 
                my ( $acc, $version ) = &_get_accession_version($id);
1155
 
                $self->element(
1156
 
                    {
1157
 
                        'Name' => 'Hit_accession',
1158
 
                        'Data' => $acc
1159
 
                    }
1160
 
                );
1161
 
 
1162
 
                if ( defined $v ) {
1163
 
                    $self->element(
1164
 
                        {
1165
 
                            'Name' => 'Hit_signif',
1166
 
                            'Data' => $v->[0]
1167
 
                        }
1168
 
                    );
1169
 
                    $self->element(
1170
 
                        {
1171
 
                            'Name' => 'Hit_score',
1172
 
                            'Data' => $v->[1]
1173
 
                        }
1174
 
                    );
1175
 
                }
1176
 
                $self->element(
1177
 
                    {
1178
 
                        'Name' => 'Hit_def',
1179
 
                        'Data' => $desc
1180
 
                    }
1181
 
                );
1182
 
                $self->end_element( { 'Name' => 'Hit' } );
1183
 
            }
1184
 
 
 
1180
            $self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1185
1181
            $self->within_element('iteration')
1186
1182
              && $self->end_element( { 'Name' => 'Iteration' } );
1187
1183
 
1202
1198
            );
1203
1199
            while ( defined( $_ = $self->_readline ) ) {
1204
1200
                if (
1205
 
                       /^(PSI)?([T]?BLAST[NPX])\s*(.+)/i
1206
 
                    || /^MEGABLAST\s*(.+)/i
 
1201
                    /^((?:\S+)?BLAST[NPX]?)\s+(.+)$/i  # NCBI BLAST, PSIBLAST
 
1202
                                                      # RPSBLAST, MEGABLAST
1207
1203
                    || /^(P?GENEWISE|HFRAME|SWN|TSWN)\s+(.+)/i    #Paracel BTK
1208
1204
                  )
1209
1205
                {
1215
1211
                elsif (/^Query=/) {
1216
1212
                    $self->_pushback($reportline) if $reportline;
1217
1213
                    $self->_pushback($_);
1218
 
 
1219
 
                    # -- Superfluous I think, but adding nonetheless
1220
 
                    $self->in_element('hsp')
1221
 
                      && $self->end_element( { 'Name' => 'Hsp' } );
1222
 
                    $self->in_element('hit')
1223
 
                      && $self->end_element( { 'Name' => 'Hit' } );
1224
 
 
1225
 
                    # --
1226
 
                    if ($bl2seq_fix) {
1227
 
                        $self->element(
1228
 
                            {
1229
 
                                'Name' => 'BlastOutput_program',
1230
 
                                'Data' => $reporttype
1231
 
                            }
1232
 
                        );
1233
 
                    }
1234
 
                    $self->end_element( { 'Name' => 'BlastOutput' } );
1235
 
                    return $self->end_document();
 
1214
                    last PARSER;
1236
1215
                }
1237
1216
 
1238
1217
                # here is where difference between wublast and ncbiblast
1846
1825
          && $self->end_element( { 'Name' => 'Hsp' } );
1847
1826
        $self->within_element('hit')
1848
1827
          && $self->end_element( { 'Name' => 'Hit' } );
 
1828
        # cleanup extra hits
 
1829
        $self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1849
1830
        $self->within_element('iteration')
1850
1831
          && $self->end_element( { 'Name' => 'Iteration' } );
1851
1832
        if ($bl2seq_fix) {
1870
1851
    $self->{'_handler_rc'} = undef;
1871
1852
}
1872
1853
 
1873
 
sub _start_iteration {
1874
 
    my $self = shift;
1875
 
    $self->start_element( { 'Name' => 'Iteration' } );
1876
 
 
1877
 
    #   $self->{'_hit_info'} = undef;
1878
 
}
1879
 
 
1880
1854
=head2 _will_handle
1881
1855
 
1882
1856
 Title   : _will_handle
1949
1923
            my $func = sprintf( "start_%s", lc $type );
1950
1924
            $self->{'_handler_rc'} = $handler->$func( $data->{'Attributes'} );
1951
1925
        }
1952
 
        else {
 
1926
        #else {
1953
1927
            #$self->debug( # changed 4/29/2006 to play nice with other event handlers
1954
1928
            #    "Bio::SearchIO::InternalParserError ".
1955
1929
            #    "\nCan't handle elements of type \'$type.\'"
1956
1930
            #);
1957
 
        }
 
1931
        #}
1958
1932
        unshift @{ $self->{'_elements'} }, $type;
1959
1933
        if ( $type eq 'result' ) {
1960
1934
            $self->{'_values'} = {};
1961
1935
            $self->{'_result'} = undef;
1962
 
        }
1963
 
        else {
1964
 
 
 
1936
        } else {
1965
1937
            # cleanup some things
1966
1938
            if ( defined $self->{'_values'} ) {
1967
1939
                foreach my $k (
1978
1950
 
1979
1951
=head2 end_element
1980
1952
 
1981
 
 Title   : start_element
 
1953
 Title   : end_element
1982
1954
 Usage   : $eventgenerator->end_element
1983
1955
 Function: Handles an end element event
1984
 
 Returns : none
 
1956
 Returns : hashref with an element's worth of data
1985
1957
 Args    : hashref with at least 2 keys 'Data' and 'Name'
1986
1958
 
1987
1959
 
1991
1963
    my ( $self, $data ) = @_;
1992
1964
    
1993
1965
    my $nm   = $data->{'Name'};
1994
 
    my $type = $MODEMAP{$nm};
 
1966
    my $type;
1995
1967
    my $rc;
1996
1968
    if ( $nm eq 'BlastOutput_program' ) {
1997
1969
        if ( $self->{'_last_data'} =~ /(t?blast[npx])/i ) {
2003
1975
    # Hsps are sort of weird, in that they end when another
2004
1976
    # object begins so have to detect this in end_element for now
2005
1977
    if ( $nm eq 'Hsp' ) {
2006
 
        foreach (qw(Hsp_qseq Hsp_midline Hsp_hseq)) {
 
1978
        foreach (qw(Hsp_qseq Hsp_midline Hsp_hseq Hsp_features)) {
2007
1979
            $self->element(
2008
1980
                {
2009
1981
                    'Name' => $_,
2010
1982
                    'Data' => $self->{'_last_hspdata'}->{$_}
2011
1983
                }
2012
 
            );
 
1984
            ) if defined $self->{'_last_hspdata'}->{$_};
2013
1985
        }
2014
1986
        $self->{'_last_hspdata'} = {};
2015
1987
        $self->element(
2051
2023
 
2052
2024
    }
2053
2025
    elsif ( $MAPPING{$nm} ) {
2054
 
 
2055
2026
        if ( ref( $MAPPING{$nm} ) =~ /hash/i ) {
2056
2027
 
2057
2028
            # this is where we shove in the data from the
2086
2057
 
2087
2058
sub element {
2088
2059
    my ( $self, $data ) = @_;
2089
 
    $self->start_element($data);
 
2060
    # Note that start element isn't needed for character data
 
2061
    # Not too SAX-y, though
 
2062
    #$self->start_element($data);
2090
2063
    $self->characters($data);
2091
2064
    $self->end_element($data);
2092
2065
}
2339
2312
    $self->{'_check_all'};
2340
2313
}
2341
2314
 
2342
 
=head2 _get_accession_version
2343
 
 
2344
 
 Title   : _get_accession_version
2345
 
 Usage   : my ($acc,$ver) = &_get_accession_version($id)
2346
 
 Function:Private function to get an accession,version pair
2347
 
           for an ID (if it is in NCBI format)
2348
 
 Returns : 2-pule of accession, version
2349
 
 Args    : ID string to process
2350
 
 
2351
 
 
2352
 
=cut
2353
 
 
2354
 
sub _get_accession_version {
2355
 
    my $id = shift;
2356
 
 
2357
 
    # handle case when this is accidently called as a class method
2358
 
    if ( ref($id) && $id->isa('Bio::SearchIO') ) {
2359
 
        $id = shift;
2360
 
    }
2361
 
    return unless defined $id;
2362
 
    my ( $acc, $version );
2363
 
    if ( $id =~ /(gb|emb|dbj|sp|pdb|bbs|ref|lcl)\|(.*)\|(.*)/ ) {
2364
 
        ( $acc, $version ) = split /\./, $2;
2365
 
    }
2366
 
    elsif ( $id =~ /(pir|prf|pat|gnl)\|(.*)\|(.*)/ ) {
2367
 
        ( $acc, $version ) = split /\./, $3;
2368
 
    }
2369
 
    else {
2370
 
 
2371
 
        #punt, not matching the db's at ftp://ftp.ncbi.nih.gov/blast/db/README
2372
 
        #Database Name                     Identifier Syntax
2373
 
        #============================      ========================
2374
 
        #GenBank                           gb|accession|locus
2375
 
        #EMBL Data Library                 emb|accession|locus
2376
 
        #DDBJ, DNA Database of Japan       dbj|accession|locus
2377
 
        #NBRF PIR                          pir||entry
2378
 
        #Protein Research Foundation       prf||name
2379
 
        #SWISS-PROT                        sp|accession|entry name
2380
 
        #Brookhaven Protein Data Bank      pdb|entry|chain
2381
 
        #Patents                           pat|country|number
2382
 
        #GenInfo Backbone Id               bbs|number
2383
 
        #General database identifier           gnl|database|identifier
2384
 
        #NCBI Reference Sequence           ref|accession|locus
2385
 
        #Local Sequence identifier         lcl|identifier
2386
 
        $acc = $id;
2387
 
    }
2388
 
    return ( $acc, $version );
 
2315
# commented out, using common base class util method
 
2316
#=head2 _get_accession_version
 
2317
#
 
2318
# Title   : _get_accession_version
 
2319
# Usage   : my ($acc,$ver) = &_get_accession_version($id)
 
2320
# Function:Private function to get an accession,version pair
 
2321
#           for an ID (if it is in NCBI format)
 
2322
# Returns : 2-pule of accession, version
 
2323
# Args    : ID string to process
 
2324
#
 
2325
#
 
2326
#=cut
 
2327
#
 
2328
#sub _get_accession_version {
 
2329
#    my $id = shift;
 
2330
#
 
2331
#    # handle case when this is accidently called as a class method
 
2332
#    if ( ref($id) && $id->isa('Bio::SearchIO') ) {
 
2333
#        $id = shift;
 
2334
#    }
 
2335
#    return unless defined $id;
 
2336
#    my ( $acc, $version );
 
2337
#    if ( $id =~ /(gb|emb|dbj|sp|pdb|bbs|ref|lcl)\|(.*)\|(.*)/ ) {
 
2338
#        ( $acc, $version ) = split /\./, $2;
 
2339
#    }
 
2340
#    elsif ( $id =~ /(pir|prf|pat|gnl)\|(.*)\|(.*)/ ) {
 
2341
#        ( $acc, $version ) = split /\./, $3;
 
2342
#    }
 
2343
#    else {
 
2344
#
 
2345
#        #punt, not matching the db's at ftp://ftp.ncbi.nih.gov/blast/db/README
 
2346
#        #Database Name                     Identifier Syntax
 
2347
#        #============================      ========================
 
2348
#        #GenBank                           gb|accession|locus
 
2349
#        #EMBL Data Library                 emb|accession|locus
 
2350
#        #DDBJ, DNA Database of Japan       dbj|accession|locus
 
2351
#        #NBRF PIR                          pir||entry
 
2352
#        #Protein Research Foundation       prf||name
 
2353
#        #SWISS-PROT                        sp|accession|entry name
 
2354
#        #Brookhaven Protein Data Bank      pdb|entry|chain
 
2355
#        #Patents                           pat|country|number
 
2356
#        #GenInfo Backbone Id               bbs|number
 
2357
#        #General database identifier           gnl|database|identifier
 
2358
#        #NCBI Reference Sequence           ref|accession|locus
 
2359
#        #Local Sequence identifier         lcl|identifier
 
2360
#        $acc = $id;
 
2361
#    }
 
2362
#    return ( $acc, $version );
 
2363
#}
 
2364
 
 
2365
# general private method used to make minimal hits from leftover
 
2366
# data in the hit table
 
2367
 
 
2368
sub _cleanup_hits {
 
2369
    my ($self, $hits) = @_;
 
2370
    while ( my $v = shift @{ $hits }) {
 
2371
        next unless defined $v;
 
2372
        $self->start_element( { 'Name' => 'Hit' } );
 
2373
        my $id   = $v->[2];
 
2374
        my $desc = $v->[3];
 
2375
        $self->element(
 
2376
            {
 
2377
                'Name' => 'Hit_id',
 
2378
                'Data' => $id
 
2379
            }
 
2380
        );
 
2381
        my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
 
2382
        $self->element(
 
2383
            {
 
2384
                'Name' => 'Hit_accession',
 
2385
                'Data' => $acc
 
2386
            }
 
2387
        );
 
2388
        if ( defined $v ) {
 
2389
            $self->element(
 
2390
                {
 
2391
                    'Name' => 'Hit_signif',
 
2392
                    'Data' => $v->[0]
 
2393
                }
 
2394
            );
 
2395
            if (exists $self->{'_wublast'}) {
 
2396
                $self->element(
 
2397
                    {
 
2398
                        'Name' => 'Hit_score',
 
2399
                        'Data' => $v->[1]
 
2400
                    }
 
2401
                );
 
2402
            } else {
 
2403
                $self->element(
 
2404
                    {
 
2405
                        'Name' => 'Hit_bits',
 
2406
                        'Data' => $v->[1]
 
2407
                    }
 
2408
                );
 
2409
            }
 
2410
            
 
2411
        }
 
2412
        $self->element(
 
2413
            {
 
2414
                'Name' => 'Hit_def',
 
2415
                'Data' => $desc
 
2416
            }
 
2417
        );
 
2418
        $self->end_element( { 'Name' => 'Hit' } );
 
2419
    }
2389
2420
}
2390
2421
 
 
2422
 
2391
2423
1;
2392
2424
 
2393
2425
__END__