File: blast.t

package info (click to toggle)
bioperl 1.6.1-2
links: PTS, VCS
area: main
in suites: squeeze
size: 40,768 kB
ctags: 12,005
sloc: perl: 174,299; xml: 13,923; sh: 1,941; lisp: 1,803; asm: 109; makefile: 53
file content (1607 lines) | stat: -rw-r--r-- 65,003 bytes
# -*-Perl-*- Test Harness script for Bioperl
# $Id: SearchIO_blast.t 14995 2008-11-16 06:20:00Z cjfields $

use strict;

BEGIN {
	use lib '.';
    use Bio::Root::Test;
    
    test_begin(-tests => 1093);
	
	use_ok('Bio::SearchIO');
}

my ($searchio, $result,$iter,$hit,$hsp);

$searchio = Bio::SearchIO->new('-format' => 'blast',
		'-file'   => test_input_file('ecolitst.bls'));

$result = $searchio->next_result;

is($result->database_name, 'ecoli.aa', 'database_name()');
is($result->database_entries, 4289);
is($result->database_letters, 1358990);

is($result->algorithm, 'BLASTP');
like($result->algorithm_version, qr/^2\.1\.3/);
like($result->query_name, qr/gi|1786183|gb|AAC73113.1| (AE000111) aspartokinase I,\s+homoserine dehydrogenase I [Escherichia coli]/);
is($result->query_accession, 'AAC73113.1');
is($result->query_gi, 1786183);
is($result->query_length, 820);
is($result->get_statistic('kappa'), '0.135');
is($result->get_statistic('kappa_gapped'), '0.0410');
is($result->get_statistic('lambda'), '0.319');
is($result->get_statistic('lambda_gapped'), '0.267');
is($result->get_statistic('entropy'), '0.383');
is($result->get_statistic('entropy_gapped'), '0.140');

is($result->get_statistic('dbletters'), 1358990);
is($result->get_statistic('dbentries'), 4289);
is($result->get_statistic('effective_hsplength'), 47);
is($result->get_statistic('effectivespace'), 894675611);
is($result->get_parameter('matrix'), 'BLOSUM62');
is($result->get_parameter('gapopen'), 11);
is($result->get_parameter('gapext'), 1);
is($result->get_statistic('S2'), '92');
is($result->get_statistic('S2_bits'), '40.0');
float_is($result->get_parameter('expect'), '1.0e-03');
is($result->get_statistic('num_extensions'), '82424');


my @valid = ( [ 'gb|AAC73113.1|', 820, 'AAC73113', '0', 1567, 4058],
	      [ 'gb|AAC76922.1|', 810, 'AAC76922', '1e-91', 332, 850],
	      [ 'gb|AAC76994.1|', 449, 'AAC76994', '3e-47', 184, 467]);
my $count = 0;
while( $hit = $result->next_hit ) {
    my $d = shift @valid;

    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
	float_is($hit->significance, shift @$d);
    is($hit->bits, shift @$d );
    is($hit->raw_score, shift @$d );

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 1);
            is($hsp->query->end, 820);
            is($hsp->hit->start, 1);
            is($hsp->hit->end, 820);
            is($hsp->length('total'), 820);
            is($hsp->start('hit'), $hsp->hit->start);
            is($hsp->end('query'), $hsp->query->end);
            is($hsp->strand('sbjct'), $hsp->subject->strand);# alias for hit
            float_is($hsp->evalue, 0.0);
            is($hsp->score, 4058);
            is($hsp->bits,1567);	    	    
            is(sprintf("%.2f",$hsp->percent_identity), 98.29);
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.9829);
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.9829);
            is($hsp->gaps, 0);
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is(@valid, 0);

$searchio = Bio::SearchIO->new('-format' => 'blast',
			       '-file'   => test_input_file('ecolitst.wublastp'));

$result = $searchio->next_result;

is($result->database_name, 'ecoli.aa');
is($result->database_letters, 1358990);
is($result->database_entries, 4289);
is($result->algorithm, 'BLASTP');
like($result->algorithm_version, qr/^2\.0MP\-WashU/);
like($result->query_name, qr/gi|1786183|gb|AAC73113.1| (AE000111) aspartokinase I,\s+homoserine dehydrogenase I [Escherichia coli]/);
is($result->query_accession, 'AAC73113.1');

is($result->query_length, 820);
is($result->query_gi, 1786183);
is($result->get_statistic('kappa'), 0.136);
is($result->get_statistic('lambda'), 0.319);
is($result->get_statistic('entropy'), 0.384);
is($result->get_statistic('dbletters'), 1358990);
is($result->get_statistic('dbentries'), 4289);
is($result->get_parameter('matrix'), 'BLOSUM62');
is($result->get_statistic('Frame+0_lambda_used'), '0.319');
is($result->get_statistic('Frame+0_kappa_used'), '0.136');
is($result->get_statistic('Frame+0_entropy_used'), '0.384');

is($result->get_statistic('Frame+0_lambda_computed'), '0.319');
is($result->get_statistic('Frame+0_kappa_computed'), '0.136');
is($result->get_statistic('Frame+0_entropy_computed'), '0.384');

is($result->get_statistic('Frame+0_lambda_gapped'), '0.244');
is($result->get_statistic('Frame+0_kappa_gapped'), '0.0300');
is($result->get_statistic('Frame+0_entropy_gapped'), '0.180');

@valid = ( [ 'gb|AAC73113.1|', 820, 'AAC73113', '0', 4141],
	   [ 'gb|AAC76922.1|', 810, 'AAC76922', '3.1e-86', 844],
	   [ 'gb|AAC76994.1|', 449, 'AAC76994', '2.8e-47', 483]);
$count = 0;
while( $hit = $result->next_hit ) {
    my $d = shift @valid;

    if ($count==1) {
        # Test HSP contig data returned by SearchUtils::tile_hsps()
        # Second hit has two hsps that overlap.
	
	# compare with the contig made by hand for these two contigs
	# in t/data/contig-by-hand.wublastp
	# (in this made-up file, the hsps from ecolitst.wublastp
        #  were aligned and contiged, and Length, Identities, Positives 
	#  were counted, by a human (maj) )
	
	my $hand_hit = Bio::SearchIO->new(
	    -format=>'blast', 
	    -file=>test_input_file('contig-by-hand.wublastp')
	    )->next_result->next_hit;
	my $hand_hsp = $hand_hit->next_hsp;
	my @hand_qrng = $hand_hsp->range('query');
	my @hand_srng = $hand_hsp->range('hit');
	my @hand_matches = $hand_hit->matches;

        my($qcontigs, $scontigs) = Bio::Search::SearchUtils::tile_hsps($hit);
        # Query contigs
        is($qcontigs->[0]->{'start'}, $hand_qrng[0]);
        is($qcontigs->[0]->{'stop'}, $hand_qrng[1]);
        is($qcontigs->[0]->{'iden'}, $hand_matches[0]);
        is($qcontigs->[0]->{'cons'}, $hand_matches[1]);
        # Subject contigs
        is($scontigs->[0]->{'start'}, $hand_srng[0]);
        is($scontigs->[0]->{'stop'}, $hand_srng[1]);
        is($scontigs->[0]->{'iden'}, $hand_matches[0]);
        is($scontigs->[0]->{'cons'}, $hand_matches[1]);
    }

    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->raw_score, shift @$d );

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 1);
            is($hsp->query->end, 820);
            is($hsp->hit->start, 1);
            is($hsp->hit->end, 820);
            is($hsp->length('total'), 820);
            
            float_is($hsp->evalue, 0.0);
            float_is($hsp->pvalue, '0.0');
            is($hsp->score, 4141);
            is($hsp->bits,1462.8);	    	    
            is($hsp->percent_identity, 100);
            is($hsp->frac_identical('query'), 1.00);
            is($hsp->frac_identical('hit'), 1.00);
            is($hsp->gaps, 0);
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is(@valid, 0);

# test that add hit really works properly for BLAST objects
# bug 1611
my @hits = $result->hits;
$result->add_hit($hits[0]);
is($result->num_hits, @hits + 1);

# test WU-BLAST -noseqs option
$searchio = Bio::SearchIO->new('-format' => 'blast',
			       '-file'   => test_input_file('ecolitst.noseqs.wublastp'));

$result = $searchio->next_result;

is($result->database_name, 'ecoli.aa');
is($result->database_letters, 1358990);
is($result->database_entries, 4289);
is($result->algorithm, 'BLASTP');
like($result->algorithm_version, qr/^2\.0MP\-WashU/);
like($result->query_name, qr/gi|1786183|gb|AAC73113.1| (AE000111) aspartokinase I,\s+homoserine dehydrogenase I [Escherichia coli]/);
is($result->query_accession, 'AAC73113.1');
is($result->query_gi, 1786183);

is($result->query_length, 820);
is($result->get_statistic('kappa'), 0.135);
is($result->get_statistic('lambda'), 0.319);
is($result->get_statistic('entropy'), 0.384);
is($result->get_statistic('dbletters'), 1358990);
is($result->get_statistic('dbentries'), 4289);
is($result->get_parameter('matrix'), 'BLOSUM62');
is($result->get_statistic('Frame+0_lambda_used'), '0.319');
is($result->get_statistic('Frame+0_kappa_used'), '0.135');
is($result->get_statistic('Frame+0_entropy_used'), '0.384');

is($result->get_statistic('Frame+0_lambda_computed'), '0.319');
is($result->get_statistic('Frame+0_kappa_computed'), '0.135');
is($result->get_statistic('Frame+0_entropy_computed'), '0.384');

is($result->get_statistic('Frame+0_lambda_gapped'), '0.244');
is($result->get_statistic('Frame+0_kappa_gapped'), '0.0300');
is($result->get_statistic('Frame+0_entropy_gapped'), '0.180');

@valid = ( [ 'gb|AAC73113.1|', 820, 'AAC73113', '0', 4141],
	   [ 'gb|AAC76922.1|', 810, 'AAC76922', '6.6e-93', 907],
	   [ 'gb|AAC76994.1|', 449, 'AAC76994', '2.8e-47', 483]);
$count = 0;
while( $hit = $result->next_hit ) {
    my $d = shift @valid;

    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->raw_score, shift @$d );

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 1);
            is($hsp->query->end, 820);
            is($hsp->hit->start, 1);
            is($hsp->hit->end, 820);
            is($hsp->length('total'), 820);
            
            float_is($hsp->evalue, 0.);
            float_is($hsp->pvalue , '0.');
            is($hsp->score, 4141);
            is($hsp->bits,1462.8);	    	    
            is($hsp->percent_identity, 100);
            is($hsp->frac_identical('query'), 1.00);
            is($hsp->frac_identical('hit'), 1.00);
            is($hsp->gaps, 0);
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is(@valid, 0);

# test tblastx 
$searchio = Bio::SearchIO->new('-format' => 'blast',
			       '-file'   => test_input_file('HUMBETGLOA.tblastx'));

$result = $searchio->next_result;
is($result->database_name, 'ecoli.nt');
is($result->database_letters, 4662239);
is($result->database_entries, 400);
is($result->algorithm, 'TBLASTX');
like($result->algorithm_version, qr/^2\.1\.2/);
is($result->query_name, 'HUMBETGLOA');
is($result->query_description, 'Human haplotype C4 beta-globin gene, complete cds.');
is($result->query_length, 3002);
is($result->get_statistic('kappa'), 0.135);
is($result->get_statistic('lambda'), 0.318);
is($result->get_statistic('entropy'), 0.401);
is($result->get_statistic('dbletters'), 4662239);
is($result->get_statistic('dbentries'), 400);
is($result->get_statistic('T'), 13);
is($result->get_statistic('X1'), 16);
is($result->get_statistic('X1_bits'), 7.3);
is($result->get_statistic('X2'), 0);
is($result->get_statistic('X2_bits'), '0.0');
is($result->get_statistic('S1'), 41);
is($result->get_statistic('S1_bits'), 21.7);
is($result->get_statistic('S2'), 53);
is($result->get_statistic('S2_bits'), 27.2);

is($result->get_statistic('decayconst'), 0.1);

is($result->get_parameter('matrix'), 'BLOSUM62');

@valid = ( [ 'gb|AE000479.1|AE000479', 10934, 'AE000479', '0.13', 33.6, 67],
	   [ 'gb|AE000302.1|AE000302', 10264, 'AE000302', '0.61', 31.3, 62],
	   [ 'gb|AE000277.1|AE000277', 11653, 'AE000277', '0.84', 30.8, 61]);
$count = 0;

while( $hit = $result->next_hit ) {
    my $d = shift @valid;
    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->bits, shift @$d );
    is($hit->raw_score, shift @$d );

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 1057);
            is($hsp->query->end, 1134);
            is($hsp->query->strand, 1);
            is($hsp->strand('query'), $hsp->query->strand);
            is($hsp->hit->end, 5893);
            is($hsp->hit->start, 5816);
            is($hsp->hit->strand, -1);
            is($hsp->strand('sbjct'), $hsp->subject->strand);
            is($hsp->length('total'), 26);
            
            float_is($hsp->evalue, 0.13);
            is($hsp->score, 67);
            is($hsp->bits,33.6);
            is(sprintf("%.2f",$hsp->percent_identity), 42.31);
            is(sprintf("%.4f",$hsp->frac_identical('query')), '0.4231');
            is(sprintf("%.4f",$hsp->frac_identical('hit')), '0.4231');
            is($hsp->query->frame(), 0);
            is($hsp->hit->frame(), 1);
            is($hsp->gaps, 0);	    
            is($hsp->query_string, 'SAYWSIFPPLGCWWSTLGPRGSLSPL');
            is($hsp->hit_string, 'AAVWALFPPVGSQWGCLASQWRTSPL');
            is($hsp->homology_string, '+A W++FPP+G  W  L  +   SPL');
            # changed to reflect positional ambiguities, note extra flag
            is(join(' ', $hsp->seq_inds('query', 'nomatch',1)), '1063-1065 1090-1095 1099-1104 1108-1113 1117-1125');
            is(join(' ', $hsp->seq_inds('hit', 'nomatch',1)), '5825-5833 5837-5842 5846-5851 5855-5860 5885-5887');            
            is($hsp->ambiguous_seq_inds, 'query/subject');
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is(@valid, 0);

# test for MarkW bug in blastN

$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('a_thaliana.blastn'));


$result = $searchio->next_result;
is($result->database_name, 'All GenBank+EMBL+DDBJ+PDB sequences (but no EST, STS, GSS,or phase 0, 1 or 2 HTGS sequences) ');
is($result->database_letters, 4677375331);
is($result->database_entries, 1083200);
is($result->algorithm, 'BLASTN');
like($result->algorithm_version, qr/^2\.2\.1/);
is($result->query_name, '');
is($result->query_length, 60);
is($result->get_parameter('gapopen'), 5);
is($result->get_parameter('gapext'), 2);
is($result->get_parameter('ktup'), undef);

is($result->get_statistic('lambda'), 1.37);
is($result->get_statistic('kappa'), 0.711);
is($result->get_statistic('entropy'),1.31 );
is($result->get_statistic('T'), 0);
is($result->get_statistic('A'), 30);
is($result->get_statistic('X1'), '6');
is($result->get_statistic('X1_bits'), 11.9);
is($result->get_statistic('X2'), 15);
is($result->get_statistic('X2_bits'), 29.7);
is($result->get_statistic('S1'), 12);
is($result->get_statistic('S1_bits'), 24.3);
is($result->get_statistic('S2'), 17);
is($result->get_statistic('S2_bits'), 34.2);

is($result->get_statistic('dbentries'), 1083200);

@valid = ( [ 'gb|AY052359.1|', 2826, 'AY052359', '3e-18', 95.6, 48, 1, 60, 
	     '1.0000'],
	   [ 'gb|AC002329.2|AC002329', 76170, 'AC002329', '3e-18', 95.6, 48, 1, 60, 
	     '1.0000' ],
	   [ 'gb|AF132318.1|AF132318', 5383, 'AF132318', '0.04', 42.1, 21, 35, 55, 
	     '0.3500']);
$count = 0;

while( my $hit = $result->next_hit ) {
    my $d = shift @valid;
    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->bits, shift @$d );
    is($hit->raw_score, shift @$d );
    is($hit->start, shift @$d);
    is($hit->end,shift @$d);    
    is(sprintf("%.4f",$hit->frac_aligned_query), shift @$d);
    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 1);
            is($hsp->query->end, 60);
            is($hsp->query->strand, 1);
            is($hsp->hit->start, 154);
            is($hsp->hit->end, 212);
            is($hsp->hit->strand, 1);
            is($hsp->length('total'), 60);	    
            float_is($hsp->evalue, 3e-18);
            is($hsp->score, 48);
            is($hsp->bits,95.6);
            is(sprintf("%.2f",$hsp->percent_identity), 96.67);
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.9667);
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.9831);
            is($hsp->query->frame(), 0);
            is($hsp->hit->frame(), 0);
            is($hsp->query->seq_id, undef);
            is($hsp->hit->seq_id, 'gb|AY052359.1|');			
            is($hsp->gaps('query'), 0);
            is($hsp->gaps('hit'), 1);
            is($hsp->gaps, 1);	    
            is($hsp->query_string, 'aggaatgctgtttaattggaatcgtacaatggagaatttgacggaaatagaatcaacgat');
            is($hsp->hit_string, 'aggaatgctgtttaattggaatca-acaatggagaatttgacggaaatagaatcaacgat');
            is($hsp->homology_string, '|||||||||||||||||||||||  |||||||||||||||||||||||||||||||||||');
			my $aln = $hsp->get_aln;
            is(sprintf("%.2f", $aln->overall_percentage_identity), 96.67);
            is(sprintf("%.2f",$aln->percentage_identity), 98.31);
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
} 
is(@valid, 0);

#WU-BlastX test

$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('dnaEbsub_ecoli.wublastx'));

$result = $searchio->next_result;
is($result->database_name, 'ecoli.aa');
is($result->database_letters, 1358990);
is($result->database_entries, 4289);
is($result->algorithm, 'BLASTX');
like($result->algorithm_version, qr/^2\.0MP\-WashU/);
is($result->query_name, 'gi|142864|gb|M10040.1|BACDNAE');
is($result->query_description, 'B.subtilis dnaE gene encoding DNA primase, complete cds');
is($result->query_accession, 'M10040.1');
is($result->query_gi, 142864);
is($result->query_length, 2001);
is($result->get_parameter('matrix'), 'blosum62');

is($result->get_statistic('lambda'), 0.318);
is($result->get_statistic('kappa'), 0.135);
is($result->get_statistic('entropy'),0.401 );

is($result->get_statistic('dbentries'), 4289);

@valid = ( [ 'gi|1789447|gb|AAC76102.1|', 581, 'AAC76102', '1.1e-74', 671]);
$count = 0;

while( my $hit = $result->next_hit ) {
    my $d = shift @valid;
    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->raw_score, shift @$d );
    is(sprintf("%.4f",$hit->frac_identical('query')), '0.3640');
    is(sprintf("%.4f",$hit->frac_identical('hit')), '0.3660');
    is(sprintf("%.4f",$hit->frac_conserved('query')), '0.5370');
    is(sprintf("%.4f",$hit->frac_conserved('hit')), '0.5400');
    is(sprintf("%.4f",$hit->frac_aligned_query), '0.6200');
    is(sprintf("%.4f",$hit->frac_aligned_hit), '0.7100');

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 21);
            is($hsp->query->end, 1265);
            is($hsp->query->strand, 1);
            is($hsp->hit->start, 1);
            is($hsp->hit->end, 413);
            is($hsp->hit->strand, 0);
            is($hsp->length('total'), 421);	    
            float_is($hsp->evalue, 1.1e-74);
            float_is($hsp->pvalue, '1.1e-74');
            is($hsp->score,671);
            is($hsp->bits,265.8);
            is(sprintf("%.2f",$hsp->percent_identity), 35.87);
            
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.3639);
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.3656);
            is(sprintf("%.4f",$hsp->frac_conserved('query')), 0.5373);
            is(sprintf("%.2f",$hsp->frac_conserved('hit')), 0.54);
            
            is(sprintf("%.4f",$hsp->frac_identical('hsp')), 0.3587);
            is(sprintf("%.4f",$hsp->frac_conserved('hsp')), 0.5297);
            
            is($hsp->query->frame(), 2);
            is($hsp->hit->frame(), 0);
            is($hsp->gaps('query'), 6);
            is($hsp->gaps('hit'), 8);
            is($hsp->gaps, 14);	    
            is($hsp->query_string, 'MGNRIPDEIVDQVQKSADIVEVIGDYVQLKKQGRNYFGLCPFHGESTPSFSVSPDKQIFHCFGCGAGGNVFSFLRQMEGYSFAESVSHLADKYQIDFPDDITVHSGARP---ESSGEQKMAEAHELLKKFYHHLLINTKEGQEALDYLLSRGFTKELINEFQIGYALDSWDFITKFLVKRGFSEAQMEKAGLLIRREDGSGYFDRFRNRVMFPIHDHHGAVVAFSGRALGSQQPKYMNSPETPLFHKSKLLYNFYKARLHIRKQERAVLFEGFADVYTAVSSDVKESIATMGTSLTDDHVKILRRNVEEIILCYDSDKAGYEATLKASELL---QKKGCKVRVAMIPDGLDPDDYIKKFGGEKFKNDIIDASVTVMAFKMQYFRKGKNLSDEGDRLAYIKDVLKEISTLSGSLEQEVYVKQ');
            is($hsp->hit_string, 'MAGRIPRVFINDLLARTDIVDLIDARVKLKKQGKNFHACCPFHNEKTPSFTVNGEKQFYHCFGCGAHGNAIDFLMNYDKLEFVETVEELAAMHNLEVPFE----AGSGPSQIERHQRQTLYQLMDGLNTFYQQSL-QQPVATSARQYLEKRGLSHEVIARFAIGFAPPGWDNVLKRFGGNPENRQSLIDAGMLVTNDQGRSY-DRFRERVMFPIRDKRGRVIGFGGRVLGNDTPKYLNSPETDIFHKGRQLYGLYEAQQDNAEPNRLLVVEGYMDVVALAQYGINYAVASLGTSTTADHIQLLFRATNNVICCYDGDRAGRDAAWRALETALPYMTDGRQLRFMFLPDGEDPDTLVRKEGKEAFEARM-EQAMPLSAFLFNSLMPQVDLSTPDGRARLSTLALPLISQVPGETLR-IYLRQ');
            is($hsp->homology_string, 'M  RIP   ++ +    DIV++I   V+LKKQG+N+   CPFH E TPSF+V+ +KQ +HCFGCGA GN   FL   +   F E+V  LA  + ++ P +    +G+ P   E    Q + +  + L  FY   L        A  YL  RG + E+I  F IG+A   WD + K       +   +  AG+L+  + G  Y DRFR RVMFPI D  G V+ F GR LG+  PKY+NSPET +FHK + LY  Y+A+    +  R ++ EG+ DV       +  ++A++GTS T DH+++L R    +I CYD D+AG +A  +A E        G ++R   +PDG DPD  ++K G E F+  + + ++ + AF         +LS    R       L  IS + G   + +Y++Q');
            is(join(' ', $hsp->seq_inds('query', 'nomatch',1)), '24-29 39-47 54-56 60-71 90-98 129-137 150-152 156-158 180-182 192-194 219-221 228-236 243-251 255-263 267-269 279-284 291-296 300-302 309-311 315-317 321-332 342-344 351-362 366-368 372-374 378-383 387-389 393-398 405-413 417-440 444-449 456-461 468-470 474-476 486-491 495-497 510-518 525-527 531-533 537-557 561-569 573-578 594-599 603-605 609-614 618-620 633-635 654-656 660-665 669-671 678-680 684-686 693-695 705-710 738-740 753-755 759-761 768-773 786-797 801-806 810-812 819-821 831-833 840-860 864-869 894-896 900-902 921-923 927-938 945-947 957-959 972-974 981-986 993-995 999-1013 1017-1019 1029-1037 1050-1052 1062-1067 1077-1079 1083-1085 1089-1091 1098-1103 1107-1109 1113-1115 1122-1124 1128-1130 1137-1163 1173-1184 1188-1208 1212-1217 1224-1226 1230-1232 1236-1244 1248-1250');
            is(join(' ', $hsp->seq_inds('query', 'mismatch',1)), '24-29 39-47 54-56 60-71 90-98 129-137 150-152 156-158 180-182 192-194 219-221 228-236 243-251 255-263 267-269 279-284 291-296 300-302 309-311 315-317 342-344 351-362 366-368 372-374 378-383 387-389 393-398 405-413 420-440 444-449 456-461 468-470 474-476 486-491 495-497 510-518 525-527 531-533 537-557 561-569 573-578 594-599 603-605 609-614 633-635 654-656 660-665 669-671 678-680 684-686 693-695 705-710 738-740 753-755 759-761 768-773 786-797 801-806 810-812 819-821 831-833 840-860 864-869 894-896 900-902 921-923 927-938 945-947 957-959 972-974 981-986 993-995 999-1013 1017-1019 1029-1037 1050-1052 1062-1067 1077-1079 1083-1085 1089-1091 1098-1103 1113-1115 1122-1124 1128-1130 1137-1163 1173-1184 1188-1208 1212-1217 1224-1226 1230-1232 1236-1244');
            is(join(' ', $hsp->seq_inds('hit', 'nomatch',1)), '2 3 7-9 12 14-17 24-26 37-39 44 46 54 58 67 70-72 75-77 79-81 83 87 88 91 92 94 97 99 104 106-108 110-113 115 117 119 120 122 124 125 128-130 132-138 140 141 144 145 148 150 154 155 157 162-164 167 169 171-177 179-181 183 184 190 191 193 195 196 202 209 211 212 214 217 219 222 226 227 237 242 244 247 248 253-256 258 259 261 264 268 271-277 279 280 289 291 298 300-303 306 310 315 318 319 322 324-331 333 337-339 344 348 349 353 355 357 360 361 364 367 369 372-380 384-387 389-395 397 398 401 403 405-407');
            is(join(' ', $hsp->seq_inds('hit', 'mismatch',1)), '2 3 7-9 12 14-17 24-26 37-39 44 46 54 58 67 70-72 75-77 79-81 83 87 88 91 92 94 97 99 104 110-113 115 117 119 120 122 124 125 128-130 132-138 140 141 144 145 148 150 154 155 157 162-164 167 169 171-177 179-181 183 184 190 191 193 195 196 202 209 211 212 214 217 219 222 226 227 237 242 244 247 248 253-256 258 259 261 264 268 271-277 279 280 289 291 298 300-303 306 310 315 318 319 322 324 325 329-331 333 337-339 344 348 349 353 355 357 360 361 364 367 369 372-380 384-387 389-395 397 398 401 403 405-407');
            is(join(' ', $hsp->seq_inds('query', 'gaps',1)), '347 1004');
            is(join(' ', $hsp->seq_inds('hit', 'gaps',1)), '100 131 197 362 408');
            is($hsp->ambiguous_seq_inds, 'query');
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is(@valid, 0);

#Trickier WU-Blast
$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('tricky.wublast'));
$result = $searchio->next_result;
my $hits_left = 1;
while (my $hit = $result->next_hit) {
	# frac_aligned_hit used to be over 1, frac_identical & frac_conserved are still too wrong
	TODO: {
        local $TODO = 'frac_identical & frac_conserved are still too wrong';
        cmp_ok sprintf("%.3f",$hit->frac_identical), '>', 0.9;
        cmp_ok sprintf("%.3f",$hit->frac_conserved), '<=', 1;
    }
    is(sprintf("%.2f",$hit->frac_aligned_query), '0.92');
    is(sprintf("%.2f",$hit->frac_aligned_hit), '0.91');
    $hits_left--;
}
is($hits_left, 0);

# More frac_ method testing, this time on ncbi blastn
$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('frac_problems.blast'));
my @expected = ("1.000", "0.943");
while (my $result = $searchio->next_result) {
    my $hit = $result->next_hit;
    is($hit->frac_identical, shift @expected);
}
is(@expected, 0);

# And even more: frac_aligned_query should never be over 1!
$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('frac_problems2.blast'));
$result = $searchio->next_result;
$hit = $result->next_hit;
is $hit->frac_aligned_query, 0.97;

# Also, start and end should be sane
$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('frac_problems3.blast'));
$result = $searchio->next_result;
$hit = $result->next_hit;
is $hit->start('sbjct'), 207;
is $hit->end('sbjct'), 1051;

#WU-TBlastN test

$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('dnaEbsub_ecoli.wutblastn'));

$result = $searchio->next_result;
is($result->database_name, 'ecoli.nt');
is($result->database_letters, 4662239);
is($result->database_entries, 400);
is($result->algorithm, 'TBLASTN');
like($result->algorithm_version, qr/^2\.0MP\-WashU/);
is($result->query_name, 'gi|142865|gb|AAA22406.1|');
is($result->query_description, 'DNA primase');
is($result->query_accession, 'AAA22406.1');
is($result->query_gi, 142865);
is($result->query_length, 603);
is($result->get_parameter('matrix'), 'blosum62');

is($result->get_statistic('lambda'), '0.320');
is($result->get_statistic('kappa'), 0.136);
is($result->get_statistic('entropy'),0.387 );

is($result->get_statistic('dbentries'), 400);

@valid = ( [ 'gi|1789441|gb|AE000388.1|AE000388', 10334, 'AE000388', '1.4e-73', 671]);
$count = 0;

while( my $hit = $result->next_hit ) {
    my $d = shift @valid;
    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->raw_score, shift @$d );

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 1);
            is($hsp->query->end, 415);
            is($hsp->query->strand, 0);
            is($hsp->hit->start, 4778);
            is($hsp->hit->end, 6016);
            is($hsp->hit->strand, 1);
            is($hsp->length('total'), 421);	    
            float_is($hsp->evalue, 1.4e-73);
            float_is($hsp->pvalue,1.4e-73);
            is($hsp->score,671);
            is($hsp->bits,265.8);
            is(sprintf("%.2f",$hsp->percent_identity), 35.87);
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.3656);	    
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.3639);
            is(sprintf("%.4f",$hsp->frac_conserved('hsp')), 0.5297);
            is($hsp->query->frame(), 0);
            is($hsp->hit->frame(), 1);
            is($hsp->gaps('query'), 6);
            is($hsp->gaps('hit'), 8);
            is($hsp->gaps, 14);	    
            is($hsp->query_string, 'MGNRIPDEIVDQVQKSADIVEVIGDYVQLKKQGRNYFGLCPFHGESTPSFSVSPDKQIFHCFGCGAGGNVFSFLRQMEGYSFAESVSHLADKYQIDFPDDITVHSGARP---ESSGEQKMAEAHELLKKFYHHLLINTKEGQEALDYLLSRGFTKELINEFQIGYALDSWDFITKFLVKRGFSEAQMEKAGLLIRREDGSGYFDRFRNRVMFPIHDHHGAVVAFSGRALGSQQPKYMNSPETPLFHKSKLLYNFYKARLHIRKQERAVLFEGFADVYTAVSSDVKESIATMGTSLTDDHVKILRRNVEEIILCYDSDKAGYEATLKASELL---QKKGCKVRVAMIPDGLDPDDYIKKFGGEKFKNDIIDASVTVMAFKMQYFRKGKNLSDEGDRLAYIKDVLKEISTLSGSLEQEVYVKQ');
            is($hsp->hit_string, 'MAGRIPRVFINDLLARTDIVDLIDARVKLKKQGKNFHACCPFHNEKTPSFTVNGEKQFYHCFGCGAHGNAIDFLMNYDKLEFVETVEELAAMHNLEVPFE----AGSGPSQIERHQRQTLYQLMDGLNTFYQQSL-QQPVATSARQYLEKRGLSHEVIARFAIGFAPPGWDNVLKRFGGNPENRQSLIDAGMLVTNDQGRSY-DRFRERVMFPIRDKRGRVIGFGGRVLGNDTPKYLNSPETDIFHKGRQLYGLYEAQQDNAEPNRLLVVEGYMDVVALAQYGINYAVASLGTSTTADHIQLLFRATNNVICCYDGDRAGRDAAWRALETALPYMTDGRQLRFMFLPDGEDPDTLVRKEGKEAFEARM-EQAMPLSAFLFNSLMPQVDLSTPDGRARLSTLALPLISQVPGETLR-IYLRQ');
            is($hsp->homology_string, 'M  RIP   ++ +    DIV++I   V+LKKQG+N+   CPFH E TPSF+V+ +KQ +HCFGCGA GN   FL   +   F E+V  LA  + ++ P +    +G+ P   E    Q + +  + L  FY   L        A  YL  RG + E+I  F IG+A   WD + K       +   +  AG+L+  + G  Y DRFR RVMFPI D  G V+ F GR LG+  PKY+NSPET +FHK + LY  Y+A+    +  R ++ EG+ DV       +  ++A++GTS T DH+++L R    +I CYD D+AG +A  +A E        G ++R   +PDG DPD  ++K G E F+  + + ++ + AF         +LS    R       L  IS + G   + +Y++Q');
            is(join(' ', $hsp->seq_inds('query', 'nomatch',1)), '2 3 7-9 12 14-17 24-26 37-39 44 46 54 58 67 70-72 75-77 79-81 83 87 88 91 92 94 97 99 101-104 108 111-114 116 118 120 121 123 125 126 129-131 133-140 142 143 146 147 150 152 156 157 159 164-166 169 171 173-179 181-183 185 186 192 193 195 197 198 200 205 212 214 215 217 220 222 225 229 230 240 245 247 250 251 256-259 261 262 264 267 271 274-280 282 283 292 294 301 303-306 309 313 318 321 322 325 327-331 333 337-339 344 348 349 353 355 357 360 361 363 365 368 370 373-381 385-388 390-396 398 399 402 404 406-408 410');
            is(join(' ', $hsp->seq_inds('hit', 'nomatch',1)), '4781-4786 4796-4804 4811-4813 4817-4828 4847-4855 4886-4894 4907-4909 4913-4915 4937-4939 4949-4951 4976-4978 4985-4993 5000-5008 5012-5020 5024-5026 5036-5041 5048-5053 5057-5059 5066-5068 5072-5074 5087-5089 5093-5101 5105-5116 5120-5122 5126-5128 5132-5137 5141-5143 5147-5152 5159-5167 5171-5191 5195-5200 5207-5212 5219-5221 5225-5227 5237-5242 5246-5248 5261-5269 5276-5278 5282-5284 5288-5308 5312-5320 5324-5329 5345-5350 5354-5356 5360-5365 5381-5383 5402-5404 5408-5413 5417-5419 5426-5428 5432-5434 5441-5443 5453-5458 5486-5488 5501-5503 5507-5509 5516-5521 5534-5545 5549-5554 5558-5560 5567-5569 5579-5581 5588-5608 5612-5617 5642-5644 5648-5650 5669-5671 5675-5686 5693-5695 5705-5707 5720-5722 5729-5734 5741-5743 5747-5770 5774-5776 5786-5794 5807-5809 5819-5824 5834-5836 5840-5842 5846-5848 5855-5860 5867-5869 5876-5878 5882-5884 5891-5917 5927-5938 5942-5962 5966-5971 5978-5980 5984-5986 5990-5998');
            is(join(' ', $hsp->seq_inds('query', 'gaps',1)), '109 328');
            is(join(' ', $hsp->seq_inds('hit', 'gaps',1)), '5077 5170 5368 5863 6001');
            is($hsp->ambiguous_seq_inds, 'subject');
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is($count, 1);

# WU-BLAST TBLASTX
$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('dnaEbsub_ecoli.wutblastx'));

$result = $searchio->next_result;
is($result->database_name, 'ecoli.nt');
is($result->database_letters, 4662239);
is($result->database_entries, 400);
is($result->algorithm, 'TBLASTX');
like($result->algorithm_version, qr/^2\.0MP\-WashU/);
is($result->query_name, 'gi|142864|gb|M10040.1|BACDNAE');
is($result->query_description, 'B.subtilis dnaE gene encoding DNA primase, complete cds');
is($result->query_accession, 'M10040.1');
is($result->query_gi, 142864);
is($result->query_length, 2001);
is($result->get_parameter('matrix'), 'blosum62');

is($result->get_statistic('lambda'), 0.318);
is($result->get_statistic('kappa'), 0.135);
is($result->get_statistic('entropy'),0.401 );
is($result->get_statistic('dbentries'), 400);

@valid = ( [ 'gi|1789441|gb|AE000388.1|AE000388', 10334, 'AE000388', '6.4e-70', 318, 148.6],
	   [ 'gi|2367383|gb|AE000509.1|AE000509', 10589, 'AE000509', 1, 59, 29.9]
	   );
$count = 0;

while( my $hit = $result->next_hit ) {
    my $d = shift @valid;
    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    # using e here to deal with 0.9992 coming out right here as well
    float_is($hit->significance, shift @$d);
    is($hit->raw_score, shift @$d );
    is($hit->bits, shift @$d );
    if( $count == 0 ) {
        my $hspcounter = 0;
        while( my $hsp = $hit->next_hsp ) {
            $hspcounter++;
            if( $hspcounter == 3 ) {
                # let's actually look at the 3rd HSP
                is($hsp->query->start, 441);
                is($hsp->query->end, 617);
                is($hsp->query->strand, 1);
                is($hsp->hit->start, 5192);
                is($hsp->hit->end, 5368);
                is($hsp->hit->strand, 1);
                is($hsp->length('total'), 59);	    
                float_is($hsp->evalue, 6.4e-70);
                float_is($hsp->pvalue,6.4e-70);
                is($hsp->score,85);
                is($hsp->bits,41.8);
                is(sprintf("%.2f",$hsp->percent_identity), '32.20');
                is(sprintf("%.3f",$hsp->frac_identical('hit')), 0.322);
                is(sprintf("%.3f",$hsp->frac_identical('query')), 0.322);
                is(sprintf("%.4f",$hsp->frac_conserved('hsp')), 0.4746);
                is($hsp->query->frame(), 2);
                is($hsp->hit->frame(), 1);
                is($hsp->gaps('query'), 0);
                is($hsp->gaps('hit'), 0);
                is($hsp->gaps, 0);	    
                is($hsp->query_string, 'ALDYLLSRGFTKELINEFQIGYALDSWDFITKFLVKRGFSEAQMEKAGLLIRREDGSGY');
                is($hsp->hit_string, 'ARQYLEKRGLSHEVIARFAIGFAPPGWDNVLKRFGGNPENRQSLIDAGMLVTNDQGRSY');
                is($hsp->homology_string, 'A  YL  RG + E+I  F IG+A   WD + K       +   +  AG+L+  + G  Y');
                is(join(' ', $hsp->seq_inds('query', 'nomatch',1)), '444-449 456-461 468-470 474-476 486-491 495-497 510-518 525-527 531-533 537-557 561-569 573-578 594-599 603-605 609-614');
                is(join(' ', $hsp->seq_inds('hit', 'nomatch',1)), '5195-5200 5207-5212 5219-5221 5225-5227 5237-5242 5246-5248 5261-5269 5276-5278 5282-5284 5288-5308 5312-5320 5324-5329 5345-5350 5354-5356 5360-5365');
                is($hsp->ambiguous_seq_inds, 'query/subject');
                last;
            }
        }
        is($hspcounter, 3);
    }
    elsif( $count == 1 ) {
        my $hsps_to_do = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 587);
            is($hsp->query->end, 706);
            is($hsp->query->strand, -1);
            is($hsp->hit->start, 4108);
            is($hsp->hit->end, 4227);
            is($hsp->hit->strand, -1);
            is($hsp->length('total'), 40);	    
            float_is($hsp->evalue, 7.1);
            float_is($hsp->pvalue , '1.00');
            is($hsp->score,59);
            is($hsp->bits,29.9);
            is(sprintf("%.2f",$hsp->percent_identity), '37.50');
            is(sprintf("%.4f",$hsp->frac_identical('hit')), '0.3750');
            is(sprintf("%.4f",$hsp->frac_identical('query')), '0.3750');
            is(sprintf("%.4f",$hsp->frac_conserved('hsp')), '0.4750');
            is($hsp->query->frame(), 2);
            is($hsp->hit->frame(), 2);
            is($hsp->gaps('query'), 0);
            is($hsp->gaps('hit'), 0);
            is($hsp->gaps, 0);
            is($hsp->query_string, 'WLPRALPEKATTAP**SWIGNMTRFLKRSKYPLPSSRLIR');
            is($hsp->hit_string, 'WLSRTTVGSSTVSPRTFWITRMKVKLSSSKVTLPSTKSTR');
            is($hsp->homology_string, 'WL R     +T +P   WI  M   L  SK  LPS++  R');
            $hsps_to_do--;
            last;
        }
        is($hsps_to_do, 0);
    }
    last if( $count++ > @valid );
}
is($count, 2);

# WU-BLAST -echofilter option test (Bug 2388)
$searchio = Bio::SearchIO->new('-format' => 'blast',
			       '-file'   => test_input_file('echofilter.wublastn'));

$result = $searchio->next_result;

is($result->database_name, 'NM_003201.fa');
is($result->database_letters, 1936);
is($result->database_entries, 1);
is($result->algorithm, 'BLASTN');
like($result->algorithm_version, qr/^2\.0MP\-WashU/);
like($result->query_name, qr/ref|NM_003201.1| Homo sapiens transcription factor A, mitochondrial \(TFAM\), mRNA/);
is($result->query_accession, 'NM_003201.1');

is($result->query_length, 1936);
is($result->get_statistic('lambda'), 0.192);
is($result->get_statistic('kappa'), 0.182);
is($result->get_statistic('entropy'), 0.357);
is($result->get_statistic('dbletters'), 1936);
is($result->get_statistic('dbentries'), 1);
is($result->get_parameter('matrix'), '+5,-4');

@valid = ( [ 'ref|NM_003201.1|', 1936, 'NM_003201', '0', 9680],);
$count = 0;
while( $hit = $result->next_hit ) {
    my $d = shift @valid;

    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->raw_score, shift @$d );

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 1);
            is($hsp->query->end, 1936);
            is($hsp->hit->start, 1);
            is($hsp->hit->end, 1936);
            is($hsp->length('total'), 1936);
            
            float_is($hsp->evalue, 0.);
            float_is($hsp->pvalue , '0.');
            is($hsp->score, 9680);
            is($hsp->bits,1458.4);	    	    
            is($hsp->percent_identity, 100);
            is($hsp->frac_identical('query'), 1.00);
            is($hsp->frac_identical('hit'), 1.00);
            is($hsp->gaps, 0);
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is(@valid, 0);


# Do a multiblast report test
$searchio = Bio::SearchIO->new('-format' => 'blast',
			       '-file'   => test_input_file('multi_blast.bls'));

@expected = qw(CATH_RAT CATL_HUMAN CATL_RAT PAPA_CARPA);
my $results_left = 4;
while( my $result = $searchio->next_result ) {
    is($result->query_name, shift @expected, "Multiblast query test");
    $results_left--;
}
is($results_left, 0);

# Test GCGBlast parsing

$searchio = Bio::SearchIO->new('-format' => 'blast',
			      '-file'   => test_input_file('test.gcgblast'));
$result = $searchio->next_result();

is($result->query_name, '/v0/people/staji002/test.gcg');
is($result->algorithm, 'BLASTP');
is($result->algorithm_version, '2.2.1 [Apr-13-2001]');
is($result->database_name, 'pir');
is($result->database_entries, 274514);
is($result->database_letters, 93460074);

$hit = $result->next_hit;
is($hit->description, 'F22B7.10 protein - Caenorhabditis elegans');
is($hit->name, 'PIR2:S44629');
is($hit->length, 628);
is($hit->accession, 'PIR2:S44629');
float_is($hit->significance, 2e-08);
is($hit->raw_score, 136 );
is($hit->bits, '57.0' );
$hsp = $hit->next_hsp;
float_is($hsp->evalue, 2e-08);
is($hsp->bits, '57.0');
is($hsp->score, 136);
is(int($hsp->percent_identity), 28);
is(sprintf("%.2f",$hsp->frac_identical('query')), 0.29);
is($hsp->frac_conserved('total'), 69/135);
is($hsp->gaps('total'), 8);
is($hsp->gaps('hit'), 6);
is($hsp->gaps('query'), 2);

is($hsp->hit->start, 342);
is($hsp->hit->end, 470);
is($hsp->query->start, 3);
is($hsp->query->end, 135);

is($hsp->query_string, 'CAAEFDFMEKETPLRYTKTXXXXXXXXXXXXXXRKIISDMWGVLAKQQTHVRKHQFDHGELVYHALQLLAYTALGILIMRLKLFLTPYMCVMASLICSRQLFGW--LFCKVHPGAIVFVILAAMSIQGSANLQTQ');
is($hsp->hit_string, 'CSAEFDFIQYSTIEKLCGTLLIPLALISLVTFVFNFVKNT-NLLWRNSEEIG----ENGEILYNVVQLCCSTVMAFLIMRLKLFMTPHLCIVAALFANSKLLGGDRISKTIRVSALVGVI-AILFYRGIPNIRQQ');
is($hsp->homology_string, 'C+AEFDF++  T  +   T                 + +   +L +    +     ++GE++Y+ +QL   T +  LIMRLKLF+TP++C++A+L  + +L G   +   +   A+V VI A +  +G  N++ Q');

#test all the database accession number formats
$searchio = Bio::SearchIO->new(-format => 'blast',
				 -file   => test_input_file('testdbaccnums.out'));
$result = $searchio->next_result;

@valid = ( ['pir||T14789','T14789','T14789','CAB53709','AAH01726'],
	   ['gb|NP_065733.1|CYT19', 'NP_065733','CYT19'],
	   ['emb|XP_053690.4|Cyt19','XP_053690'],
	   ['dbj|NP_056277.2|DKFZP586L0724','NP_056277'],
	   ['prf||XP_064862.2','XP_064862'],
	   ['pdb|BAB13968.1|1','BAB13968'],
	   ['sp|Q16478|GLK5_HUMAN','Q16478'],
	   ['pat|US|NP_002079.2','NP_002079'],
	   ['bbs|NP_079463.2|','NP_079463'],
	   ['gnl|db1|NP_002444.1','NP_002444'],
	   ['ref|XP_051877.1|','XP_051877'],
	   ['lcl|AAH16829.1|','AAH16829'],
	   ['gi|1|gb|NP_065733.1|CYT19','NP_065733'],
	   ['gi|2|emb|XP_053690.4|Cyt19','XP_053690'],
	   ['gi|3|dbj|NP_056277.2|DKFZP586L0724','NP_056277'],
	   ['gi|4|pir||T14789','T14789'],
	   ['gi|5|prf||XP_064862.2','XP_064862'],
	   ['gi|6|pdb|BAB13968.1|1','BAB13968'],
	   ['gi|7|sp|Q16478|GLK5_HUMAN','Q16478'],
	   ['gi|8|pat|US|NP_002079.2','NP_002079'],
	   ['gi|9|bbs|NP_079463.2|','NP_079463'],
	   ['gi|10|gnl|db1|NP_002444.1','NP_002444'],
	   ['gi|11|ref|XP_051877.1|','XP_051877'],
	   ['gi|12|lcl|AAH16829.1|','AAH16829'],
	   ['MY_test_ID','MY_test_ID']
	   );

$hit = $result->next_hit;
my $d = shift @valid;
is($hit->name, shift @$d);
is($hit->accession, shift @$d);
my @accnums = $hit->each_accession_number;
foreach my $a (@accnums) {
	is($a, shift @$d);
}
$d = shift @valid;
$hit = $result->next_hit;
is($hit->name, shift @$d);
is($hit->accession, shift @$d);
is($hit->locus, shift @$d);

$hits_left = 23;
while( $hit = $result->next_hit ) {
    my $d = shift @valid;
    is($hit->name, shift @$d);
    is($hit->accession, shift @$d);
    $hits_left--;
}
is($hits_left, 0);

# Parse MEGABLAST

# parse the BLAST-like output
my $infile = test_input_file('503384.MEGABLAST.2');
my $in = Bio::SearchIO->new(-file => $infile,
			   -format => 'blast'); # this is megablast blast-like output
my $r = $in->next_result;
my @dcompare = ( ['Contig3700', 5631, 396, 785,'0.0', 785, '0.0', 396, 639, 12, 
		  8723,9434, 1, 4083, 4794, -1],
                 ['Contig3997', 12734, 335, 664,'0.0', 664, '0.0', 335, 401, 0, 
		  1282, 1704, 1, 1546, 1968,-1 ],
                 ['Contig634', 858, 245, 486,'1e-136', 486, '1e-136', 245, 304, 3, 
		  7620, 7941, 1, 1, 321, -1],
                 ['Contig1853', 2314, 171, 339,'1e-91',339, '1e-91', 171, 204, 0,
		  6406, 6620, 1, 1691, 1905, 1]
    );

is($r->algorithm, 'MEGABLAST');
is($r->query_name, '503384');
is($r->query_description, '11337 bp 2 contigs');
is($r->query_length, 11337);
is($r->database_name, 'cneoA.nt ');
is($r->database_letters, 17206226);
is($r->database_entries, 4935);

$hits_left = 4;
while( my $hit = $r->next_hit ) {
    my $d = shift @dcompare;
    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->raw_score, shift @$d);
    is($hit->bits, shift @$d);
	float_is($hit->significance, shift @$d);
    
    my $hsp = $hit->next_hsp;
    is($hsp->bits, shift @$d);
    float_is($hsp->evalue, shift @$d);
    is($hsp->score, shift @$d);
    is($hsp->num_identical, shift @$d);
    is($hsp->gaps('total'), shift @$d);
    is($hsp->query->start, shift @$d);
    is($hsp->query->end, shift @$d);
    is($hsp->query->strand, shift @$d);
    is($hsp->hit->start, shift @$d);
    is($hsp->hit->end, shift @$d);
    is($hsp->hit->strand, shift @$d);
    $hits_left--;
}
is($hits_left, 0);

# Let's test RPS-BLAST

my $parser = Bio::SearchIO->new(-format => 'blast',
			       -file   => test_input_file('ecoli_domains.rpsblast'));

$r = $parser->next_result;
is($r->query_name, 'gi|1786183|gb|AAC73113.1|');
is($r->query_gi, 1786183);
is($r->num_hits, 7);
$hit = $r->next_hit;
is($hit->name, 'gnl|CDD|3919');
float_is($hit->significance, 0.064);
is($hit->bits, 28.3);
is($hit->raw_score, 63);
$hsp = $hit->next_hsp;
is($hsp->query->start, 599);
is($hsp->query->end,655);
is($hsp->hit->start,23);
is($hsp->hit->end,76);


# Test PSI-BLAST parsing

$searchio = Bio::SearchIO->new('-format' => 'blast',
			       '-file'   => test_input_file('psiblastreport.out'));

$result = $searchio->next_result;

is($result->database_name, '/home/peter/blast/data/swissprot.pr');
is($result->database_entries, 88780);
is($result->database_letters, 31984247);

is($result->algorithm, 'BLASTP');
like($result->algorithm_version, qr/^2\.0\.14/);
is($result->query_name, 'CYS1_DICDI');
is($result->query_length, 343);
is($result->get_statistic('kappa') , 0.0491);
cmp_ok($result->get_statistic('lambda'), '==', 0.270);
cmp_ok($result->get_statistic('entropy'), '==', 0.230);
is($result->get_statistic('dbletters'), 31984247);
is($result->get_statistic('dbentries'), 88780);
is($result->get_statistic('effective_hsplength'), 49);
is($result->get_statistic('effectivespace'), 8124403938);
is($result->get_parameter('matrix'), 'BLOSUM62');
is($result->get_parameter('gapopen'), 11);
is($result->get_parameter('gapext'), 1);

my @valid_hit_data = ( [ 'sp|P04988|CYS1_DICDI', 343, 'P04988', '0', 721],
		       [ 'sp|P43295|A494_ARATH', 313, 'P43295', '1e-75', 281],
		       [ 'sp|P25804|CYSP_PEA', 363, 'P25804', '1e-74', 278]);
my @valid_iter_data = ( [ 127, 127, 0, 109, 18, 0, 0, 0, 0],
			[ 157, 40, 117, 2, 38, 0, 109, 3, 5]);
my $iter_count = 0;
while( $iter = $result->next_iteration ) {
    $iter_count++;
    my $di = shift @valid_iter_data;
    is($iter->number, $iter_count);

    is($iter->num_hits, shift @$di);
    is($iter->num_hits_new, shift @$di);
    is($iter->num_hits_old, shift @$di);
    is(scalar($iter->newhits_below_threshold), shift @$di);
    is(scalar($iter->newhits_not_below_threshold), shift @$di);
    is(scalar($iter->newhits_unclassified), shift @$di);
    is(scalar($iter->oldhits_below_threshold), shift @$di);
    is(scalar($iter->oldhits_newly_below_threshold), shift @$di);
    is(scalar($iter->oldhits_not_below_threshold), shift @$di);

    my $hit_count = 0;
    if ($iter_count == 1) {
        while( $hit = $result->next_hit ) {
            my $d = shift @valid_hit_data;
            
            is($hit->name, shift @$d);
            is($hit->length, shift @$d);
            is($hit->accession, shift @$d);
            float_is($hit->significance, shift @$d);
            is($hit->bits, shift @$d );
            
            if( $hit_count == 1 ) {
                my $hsps_left = 1;
                while( my $hsp = $hit->next_hsp ){
                    is($hsp->query->start, 32);
                    is($hsp->query->end, 340);
                    is($hsp->hit->start, 3);
                    is($hsp->hit->end, 307);
                    is($hsp->length('total'), 316);
                    is($hsp->start('hit'), $hsp->hit->start);
                    is($hsp->end('query'), $hsp->query->end);
                    is($hsp->strand('sbjct'), $hsp->subject->strand);# alias for hit
                    float_is($hsp->evalue, 1e-75);
                    is($hsp->score, 712);
                    is($hsp->bits, 281);
                    is(sprintf("%.1f",$hsp->percent_identity), 46.5);
                    is(sprintf("%.4f",$hsp->frac_identical('query')), 0.4757);
                    is(sprintf("%.3f",$hsp->frac_identical('hit')), 0.482);
                    is($hsp->gaps, 18);
                    $hsps_left--;
                }
                is($hsps_left, 0);
            }
            last if( $hit_count++ > @valid_hit_data );
        }
    }
}
is(@valid_hit_data, 0);
is(@valid_iter_data, 0);

# Test filtering

$searchio = Bio::SearchIO->new( '-format' => 'blast', 
                                '-file'   => test_input_file('ecolitst.bls'),
                                '-signif' => 1e-100);

@valid = qw(gb|AAC73113.1|);
$r = $searchio->next_result;

while( my $hit = $r->next_hit ) {
    is($hit->name, shift @valid);
}

$searchio = Bio::SearchIO->new( '-format' => 'blast', 
                                '-file'   => test_input_file('ecolitst.bls'),
                                '-score' => 100);

@valid = qw(gb|AAC73113.1| gb|AAC76922.1| gb|AAC76994.1|);
$r = $searchio->next_result;

while( my $hit = $r->next_hit ) {
    is($hit->name, shift @valid);
}
is(@valid, 0);

$searchio = Bio::SearchIO->new( '-format' => 'blast', 
                                '-file'   => test_input_file('ecolitst.bls'),
                                '-bits' => 200);

@valid = qw(gb|AAC73113.1| gb|AAC76922.1|);
$r = $searchio->next_result;

while( my $hit = $r->next_hit ) {
    is($hit->name, shift @valid);
}
is(@valid, 0);


my $filt_func = sub{ my $hit=shift; 
                     $hit->frac_identical('query') >= 0.31
                     };

$searchio = Bio::SearchIO->new( '-format' => 'blast', 
                                '-file'   => test_input_file('ecolitst.bls'),
                                '-hit_filter' => $filt_func);

@valid = qw(gb|AAC73113.1| gb|AAC76994.1|);
$r = $searchio->next_result;

while( my $hit = $r->next_hit ) {
    is($hit->name, shift @valid);
}
is(@valid, 0);




# bl2seq parsing testing

# this is blastp bl2seq
$searchio = Bio::SearchIO->new(-format => 'blast',
			      -file   => test_input_file('bl2seq.out'));
$result = $searchio->next_result;
isa_ok($result,'Bio::Search::Result::ResultI');
is($result->query_name, '');
is($result->algorithm, 'BLASTP');
$hit = $result->next_hit;
is($hit->name, 'ALEU_HORVU');
is($hit->length, 362);
$hsp = $hit->next_hsp;
is($hsp->score, 481);
is($hsp->bits, 191);
is(int $hsp->percent_identity, 34);
float_is($hsp->evalue, 2e-53);
is(int($hsp->frac_conserved*$hsp->length), 167);
is($hsp->query->start, 28);
is($hsp->query->end, 343);
is($hsp->hit->start, 60);
is($hsp->hit->end,360);
is($hsp->gaps, 27);

# this is blastn bl2seq 
$searchio = Bio::SearchIO->new(-format => 'blast',
			      -file   => test_input_file('bl2seq.blastn.rev'));
$result = $searchio->next_result;
isa_ok($result,'Bio::Search::Result::ResultI');
is($result->query_name, '');
is($result->algorithm, 'BLASTN');
is($result->query_length, 180);
$hit = $result->next_hit;
is($hit->length, 179);
is($hit->name, 'human');
$hsp = $hit->next_hsp;
is($hsp->score, 27);
is($hsp->bits, '54.0');
is(int $hsp->percent_identity, 88);
float_is($hsp->evalue, 2e-12);
is(int($hsp->frac_conserved*$hsp->length), 83);
is($hsp->query->start, 94);
is($hsp->query->end, 180);
is($hsp->query->strand, 1);
is($hsp->hit->strand, -1);
is($hsp->hit->start, 1);
is($hsp->hit->end,94);
is($hsp->gaps, 7);

# this is blastn bl2seq 
$searchio = Bio::SearchIO->new(-format => 'blast',
			      -file   => test_input_file('bl2seq.blastn'));
$result = $searchio->next_result;
isa_ok($result,'Bio::Search::Result::ResultI');
is($result->query_name, '');
is($result->query_length, 180);
is($result->algorithm, 'BLASTN');
$hit = $result->next_hit;
is($hit->name, 'human');
is($hit->length, 179);
$hsp = $hit->next_hsp;
is($hsp->score, 27);
is($hsp->bits, '54.0');
is(int $hsp->percent_identity, 88);
float_is($hsp->evalue, 2e-12);
is(int($hsp->frac_conserved*$hsp->length), 83);
is($hsp->query->start, 94);
is($hsp->query->end, 180);
is($hsp->query->strand, 1);
is($hsp->hit->strand, 1);
is($hsp->hit->start, 86);
is($hsp->hit->end,179);
is($hsp->gaps, 7);

# this is blastp bl2seq
$searchio = Bio::SearchIO->new(-format => 'blast',
			      -file   => test_input_file('bl2seq.bug940.out'));
$result = $searchio->next_result;
isa_ok($result,'Bio::Search::Result::ResultI');
is($result->query_name, 'zinc');
is($result->algorithm, 'BLASTP');
is($result->query_description, 'finger protein 135 (clone pHZ-17) [Homo sapiens]. neo_id RS.ctg14243-000000.6.0');
is($result->query_length, 469);
$hit = $result->next_hit;
is($hit->name, 'gi|4507985|');
is($hit->ncbi_gi, 4507985);
is($hit->description,'zinc finger protein 135 (clone pHZ-17) [Homo sapiens]. neo_id RS.ctg14243-000000.6.0');
is($hit->length, 469);
$hsp = $hit->next_hsp;
is($hsp->score, 1626);
is($hsp->bits, 637);
is(int $hsp->percent_identity, 66);
float_is($hsp->evalue, 0.0);
is(int($hsp->frac_conserved*$hsp->length), 330);
is($hsp->query->start, 121);
is($hsp->query->end, 469);
is($hsp->hit->start, 1);
is($hsp->hit->end,469);
is($hsp->gaps, 120);
ok($hit->next_hsp); # there is more than one HSP here, 
                    # make sure it is parsed at least

# cannot distinguish between blastx and tblastn reports
# so we're only testing a blastx report for now

# this is blastx bl2seq
$searchio = Bio::SearchIO->new(-format => 'blast',
			      -file   => test_input_file('bl2seq.blastx.out'));
$result = $searchio->next_result;
isa_ok($result,'Bio::Search::Result::ResultI');
is($result->query_name, 'AE000111.1');
is($result->query_description, 'Escherichia coli K-12 MG1655 section 1 of 400 of the complete genome');
is($result->algorithm, 'BLASTX');
is($result->query_length, 720);
$hit = $result->next_hit;
is($hit->name, 'AK1H_ECOLI');
is($hit->description,'P00561 Bifunctional aspartokinase/homoserine dehydrogenase I (AKI-HDI) [Includes: Aspartokinase I ; Homoserine dehydrogenase I ]');
is($hit->length, 820);
$hsp = $hit->next_hsp;
is($hsp->score, 634);
is($hsp->bits, 248);
is(int $hsp->percent_identity, 100);
float_is($hsp->evalue, 2e-70);
is(int($hsp->frac_conserved*$hsp->length), 128);
is($hsp->query->start, 1);
is($hsp->query->end, 384);
is($hsp->hit->start, 1);
is($hsp->hit->end,128);
is($hsp->gaps, 0);
is($hsp->query->frame,0);
is($hsp->hit->frame,0);
is($hsp->query->strand,-1);
is($hsp->hit->strand,0);

# this is tblastx bl2seq (self against self)
$searchio = Bio::SearchIO->new(-format => 'blast',
			      -file   => test_input_file('bl2seq.tblastx.out'));
$result = $searchio->next_result;
isa_ok($result,'Bio::Search::Result::ResultI');
is($result->query_name, 'Escherichia');
is($result->algorithm, 'TBLASTX');
is($result->query_description, 'coli K-12 MG1655 section 1 of 400 of the complete genome');
is($result->query_length, 720);
$hit = $result->next_hit;
is($hit->name, 'gi|1786181|gb|AE000111.1|AE000111');
is($hit->ncbi_gi, 1786181);
is($hit->description,'Escherichia coli K-12 MG1655 section 1 of 400 of the complete genome');
is($hit->length, 720);
$hsp = $hit->next_hsp;
is($hsp->score, 1118);
is($hsp->bits, 515);
is(int $hsp->percent_identity, 95);
float_is($hsp->evalue, 1e-151);
is(int($hsp->frac_conserved*$hsp->length), 229);
is($hsp->query->start, 1);
is($hsp->query->end, 720);
is($hsp->hit->start, 1);
is($hsp->hit->end,720);
is($hsp->gaps, 0);
is($hsp->query->frame,0);
is($hsp->hit->frame,0);
is($hsp->query->strand,1);
is($hsp->hit->strand,1);

# this is NCBI tblastn
$searchio = Bio::SearchIO->new(-format => 'blast',
										-file   => test_input_file('tblastn.out'));
$result = $searchio->next_result;
isa_ok($result,'Bio::Search::Result::ResultI');
is($result->algorithm, 'TBLASTN');
$hit = $result->next_hit;
is($hit->name,'gi|10040111|emb|AL390796.6|AL390796');

# Test Blast parsing with B=0 (WU-BLAST)
$searchio = Bio::SearchIO->new(-file   => test_input_file('no_hsps.blastp'),
			      -format => 'blast');
$result = $searchio->next_result;
is($result->query_name, 'mgri:MG00189.3');
$hit = $result->next_hit;
is($hit->name, 'mgri:MG00189.3');
is($hit->description, 'hypothetical protein 6892 8867 +');
is($hit->bits, 3098);
float_is($hit->significance, 0.);

$hit = $result->next_hit;
is($hit->name, 'fgram:FG01141.1');
is($hit->description, 'hypothetical protein 47007 48803 -');
is($hit->bits, 2182);
float_is($hit->significance, 4.2e-226);
is($result->num_hits, 415);
# Let's now test if _guess_format is doing its job correctly
my %pair = ( 'filename.blast'  => 'blast',
	     'filename.bls'    => 'blast',
	     'f.blx'           => 'blast',
	     'f.tblx'          => 'blast',
	     'fast.bls'        => 'blast',
	     'f.fasta'         => 'fasta',
	     'f.fa'            => 'fasta',
	     'f.fx'            => 'fasta',
	     'f.fy'            => 'fasta',
	     'f.ssearch'       => 'fasta',
	     'f.SSEARCH.m9'    => 'fasta',
	     'f.m9'            => 'fasta',
	     'f.psearch'       => 'fasta',
	     'f.osearch'       => 'fasta',
	     'f.exon'          => 'exonerate',
	     'f.exonerate'     => 'exonerate',
	     'f.blastxml'      => 'blastxml',
	     'f.xml'           => 'blastxml');
while( my ($file,$expformat) = each %pair ) {
    is(Bio::SearchIO->_guess_format($file),$expformat, "$expformat for $file");
}


# Test Wes Barris's reported bug when parsing blastcl3 output which
# has integer overflow

$searchio = Bio::SearchIO->new(-file => test_input_file('hsinsulin.blastcl3.blastn'),
			      -format => 'blast');
$result = $searchio->next_result;
is($result->query_name, 'human');
is($result->database_letters(), '-24016349'); 
# this is of course not the right length, but is the what blastcl3 
# reports, the correct value is
is($result->get_statistic('dbletters'),'192913178');
is($result->get_statistic('dbentries'),'1867771');


# test for links and groups being parsed out of WU-BLAST properly
$searchio = Bio::SearchIO->new(-format => 'blast',
			       -file   => test_input_file('brassica_ATH.WUBLASTN'));
ok($result = $searchio->next_result);
ok($hit = $result->next_hit);
ok($hsp = $hit->next_hsp);
is($hsp->links,'(1)-3-2');
is($hsp->query->strand, 1);
is($hsp->hit->strand, 1);
is($hsp->hsp_group, '1');

## Web blast result parsing

$searchio = Bio::SearchIO->new(-format => 'blast',
			       -file   => test_input_file('catalase-webblast.BLASTP'));
ok($result = $searchio->next_result);
ok($hit = $result->next_hit);
is($hit->name, 'gi|40747822|gb|EAA66978.1|', 'full hit name');
is($hit->accession, 'EAA66978', 'hit accession');
is($hit->ncbi_gi, 40747822);
ok($hsp = $hit->next_hsp);
is($hsp->query->start, 1, 'query start');
is($hsp->query->end, 528, 'query start');

# tests for new BLAST 2.2.13 output
$searchio = Bio::SearchIO->new(-format => 'blast',
							  -file   => test_input_file('new_blastn.txt'));

$result = $searchio->next_result;
is($result->database_name, 'All GenBank+EMBL+DDBJ+PDB sequences (but no EST, STS,GSS,environmental samples or phase 0, 1 or 2 HTGS sequences)');
is($result->database_entries, 3742891);
is($result->database_letters, 16670205594);
is($result->algorithm, 'BLASTN');
is($result->algorithm_version, '2.2.13 [Nov-27-2005]');
is($result->query_name, 'pyrR,');
is($result->query_length, 558);
is($result->get_statistic('kappa'), '0.711');
is($result->get_statistic('kappa_gapped'), '0.711');
is($result->get_statistic('lambda'), '1.37');
is($result->get_statistic('lambda_gapped'), '1.37');
is($result->get_statistic('entropy'), '1.31');
is($result->get_statistic('entropy_gapped'), '1.31');
is($result->get_statistic('dbletters'), '-509663586');
is($result->get_statistic('dbentries'), 3742891);
is($result->get_statistic('effective_hsplength'), undef);
is($result->get_statistic('effectivespace'), undef);
is($result->get_parameter('matrix'), 'blastn matrix:1 -3');
is($result->get_parameter('gapopen'), 5);
is($result->get_parameter('gapext'), 2);
is($result->get_statistic('S2'), '60');
is($result->get_statistic('S2_bits'), '119.4');
float_is($result->get_parameter('expect'), '1e-23');
is($result->get_statistic('num_extensions'), '117843');

@valid = ( [ 'gi|41400296|gb|AE016958.1|', 4829781, 'AE016958', 41400296, '6e-059', 119, 236],
	      [ 'gi|54013472|dbj|AP006618.1|', 6021225, 'AP006618', 54013472, '4e-026', 64, 127],
	      [ 'gi|57546753|dbj|BA000030.2|', 9025608, 'BA000030', 57546753, '1e-023', 60, 119]);
$count = 0;
while( $hit = $result->next_hit ) {
    my $d = shift @valid;

    is($hit->name, shift @$d);
    is($hit->length, shift @$d);
    is($hit->accession, shift @$d);
    is($hit->ncbi_gi, shift @$d);
    float_is($hit->significance, shift @$d);
    is($hit->raw_score, shift @$d );
    is($hit->bits, shift @$d );

    if( $count == 0 ) {
        my $hsps_left = 1;
        while( my $hsp = $hit->next_hsp ) {
            is($hsp->query->start, 262);
            is($hsp->query->end, 552);
            is($hsp->hit->start, 1166897);
            is($hsp->hit->end, 1167187);
            is($hsp->length('total'), 291);
            is($hsp->hit_features, 'PyrR');
            is($hsp->start('hit'), $hsp->hit->start);
            is($hsp->end('query'), $hsp->query->end);
            is($hsp->strand('sbjct'), $hsp->subject->strand);# alias for hit
            float_is($hsp->evalue, 6e-59);
            is($hsp->score, 119);
            is($hsp->bits,236);	    	    
            is(sprintf("%.2f",$hsp->percent_identity), 85.22);
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.8522);
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.8522);
            is($hsp->gaps, 0);
            $hsps_left--;
        }
        is($hsps_left, 0);
    }
    last if( $count++ > @valid );
}
is(@valid, 0);

# Bug 2189
$searchio = Bio::SearchIO->new(-format => 'blast',
							  -file   => test_input_file('blastp2215.blast'));

$result = $searchio->next_result;
is($result->database_entries, 4460989);
is($result->database_letters, 1533424333);
is($result->algorithm, 'BLASTP');
is($result->algorithm_version, '2.2.15 [Oct-15-2006]');
is($result->query_name, 'gi|15608519|ref|NP_215895.1|');
is($result->query_gi, 15608519);
is($result->query_length, 193);
@hits = $result->hits;
is(scalar(@hits), 10);
is($hits[1]->accession,'1W30');
is($hits[4]->significance,'2e-72');
is($hits[7]->bits,'254');
$result = $searchio->next_result;
is($result->database_entries, 4460989);
is($result->database_letters, 1533424333);
is($result->algorithm, 'BLASTP');
is($result->algorithm_version, '2.2.15 [Oct-15-2006]');
is($result->query_name, 'gi|15595598|ref|NP_249092.1|');
is($result->query_length, 423);
@hits = $result->hits;
is(scalar(@hits), 10);
is($hits[1]->accession,'ZP_00972546');
is($hits[2]->ncbi_gi, 116054132);
is($hits[4]->significance, '0.0');
is($hits[7]->bits, 624);

# Bug 2246
$searchio = Bio::SearchIO->new(-format => 'blast',
                               -verbose => -1,
							  -file   => test_input_file('bug2246.blast'));
$result = $searchio->next_result;
$hit = $result->next_hit;
is $hit->name, 'UniRef50_Q9X0H5';
is $hit->length, 0;
is $hit->accession, 'UniRef50_Q9X0H5';
is $hit->description, 'Cluster: Histidyl-tRNA synthetase; n=4; Thermoto...';
is $hit->bits, 23;
float_is($hit->significance, 650);

# Bug 1986
$searchio = Bio::SearchIO->new(-format => 'blast',
                               -verbose => -1,
							  -file   => test_input_file('bug1986.blastp'));
$result = $searchio->next_result;
$hit = $result->next_hit;
is $hit->name, 'ENSP00000350182';
is $hit->length, 425;
is $hit->accession, 'ENSP00000350182';
is $hit->description, 'pep:novel clone::BX322644.8:4905:15090:-1 gene:ENSG00000137397 transcript:ENST00000357569 ';
is $hit->raw_score, 301;
is $hit->bits, 120;
float_is($hit->significance, 3e-27);
$hit = $result->next_hit;
is $hit->name, 'ENSP00000327738';
is $hit->length, 468;
is $hit->accession, 'ENSP00000327738';
is $hit->description, 'pep:known-ccds chromosome:NCBI36:4:189297592:189305643:1 gene:ENSG00000184108 transcript:ENST00000332517 CCDS3851.1';
is $hit->raw_score, 289;
is $hit->bits, 115;
float_is($hit->significance, 8e-26);

# Bug 1986, pt. 2

# handle at least the first iteration with BLAST searches using databases
# containing non-unique IDs

my $file = test_input_file('bug1986.blast2');
my %unique_accs;
open (my $IN,$file) or die $!;

while (<$IN>) {
  last if (/^Sequences/);
}
$count = 1;
while (<$IN>) {
  chomp;
  next if m{^\s*$};
  next unless ($_);
  last if m{^>};
  my ($accession) = split(/\s+/);
  #print "Real Hit $count = $accession\n";
  $unique_accs{$accession}++;
  #last if ($count == 10);
  ++$count;
}
close ($IN);

is ($count, 495);
is (scalar(keys %unique_accs), 490);

my %search_accs;

$searchio = Bio::SearchIO->new(-format => 'blast',
                               -verbose => -1,
							  -file   => $file);
$result = $searchio->next_result;
$count = 1;
while (my $hit = $result->next_hit) {
	$search_accs{$hit->accession}++;
	$count++;
}

is ($count, 495);
is (scalar(keys %search_accs), 490);

is_deeply(\%unique_accs, \%search_accs);

# bug 2391 - long query names

$file = test_input_file('bug2391.megablast');

$searchio = Bio::SearchIO->new(-format => 'blast',
							  -file   => $file);
$result = $searchio->next_result;

# data is getting munged up with long names
is($result->query_name, 'c6_COX;c6_QBL;6|31508172;31503325;31478402|rs36223351|1|dbSNP|C/G');
is($result->query_description, '');
is($result->algorithm, 'MEGABLAST');

# bug 2399 - catching Expect(n) values

$file = test_input_file('bug2399.tblastn');

$searchio = Bio::SearchIO->new(-format => 'blast',
							  -file   => $file);
my $total_n = 0;
while(my $query = $searchio->next_result) {
    while(my $subject = $query->next_hit) {
        $total_n += grep{$_->n} $subject->hsps;
    }
}
is($total_n, 10);

sub cmp_evalue ($$) {
	my ($tval, $aval) = @_;
	is(sprintf("%g",$tval), sprintf("%g",$aval));
}