← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty-proposed

« back to all changes in this revision

Viewing changes to t/SearchIO/blast_pull.t

  • 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:
t/SearchIO/blast_pull.t
 
1
# -*-Perl-*- Test Harness script for Bioperl
 
2
# $Id: SearchIO_blast_pull.t 11525 2007-06-27 10:16:38Z sendu $
 
3
 
 
4
use strict;
 
5
 
 
6
BEGIN {
 
7
        use lib '.';
 
8
    use Bio::Root::Test;
 
9
    
 
10
    test_begin(-tests => 289);
 
11
        
 
12
        use_ok('Bio::SearchIO');
 
13
}
 
14
 
 
15
 
 
16
my $searchio = Bio::SearchIO->new(-format => 'blast_pull', -file => test_input_file('new_blastn.txt'));
 
17
 
 
18
my $result = $searchio->next_result;
 
19
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)';
 
20
is $result->database_entries, 3742891;
 
21
is $result->database_letters, 16670205594;
 
22
is $result->algorithm, 'BLASTN';
 
23
is $result->algorithm_version, '2.2.13 [Nov-27-2005]';
 
24
is $result->query_name, 'pyrR,';
 
25
is $result->query_length, 558;
 
26
is $result->get_statistic('kappa'), 0.711;
 
27
is $result->get_statistic('kappa_gapped'), 0.711;
 
28
is $result->get_statistic('lambda'), 1.37;
 
29
is $result->get_statistic('lambda_gapped'), 1.37;
 
30
is $result->get_statistic('entropy'), 1.31;
 
31
is $result->get_statistic('entropy_gapped'), 1.31;
 
32
is $result->get_statistic('dbletters'), -509663586;
 
33
is $result->get_statistic('dbentries'), 3742891;
 
34
is $result->get_statistic('effectivespace'), 8935230198384;
 
35
is $result->get_parameter('matrix'), 'blastn matrix:1 -3';
 
36
is $result->get_parameter('gapopen'), 5;
 
37
is $result->get_parameter('gapext'), 2;
 
38
is $result->get_statistic('S2'), 60;
 
39
is $result->get_statistic('S2_bits'), 119.4;
 
40
float_is $result->get_parameter('expect'), 1e-23;
 
41
is $result->get_statistic('num_extensions'), 117843;
 
42
 
 
43
my @valid = ( [ 'gi|41400296|gb|AE016958.1|', 4829781, 'AE016958', '6e-059', 236],
 
44
              [ 'gi|54013472|dbj|AP006618.1|', 6021225, 'AP006618', '4e-026', 127],
 
45
              [ 'gi|57546753|dbj|BA000030.2|', 9025608, 'BA000030', '1e-023', 119]);
 
46
my $count = 0;
 
47
while (my $hit = $result->next_hit ) {
 
48
    my $d = shift @valid;
 
49
 
 
50
    is $hit->name, shift @$d;
 
51
    is $hit->length, shift @$d;
 
52
    is $hit->accession, shift @$d;
 
53
    float_is($hit->significance, shift @$d);
 
54
    is $hit->raw_score, shift @$d;
 
55
 
 
56
    if( $count == 0 ) {
 
57
        my $hsps_left = 1;
 
58
        while( my $hsp = $hit->next_hsp ) {
 
59
            is $hsp->query->start, 262;
 
60
            is $hsp->query->end, 552;
 
61
            is $hsp->hit->start, 1166897;
 
62
            is $hsp->hit->end, 1167187;
 
63
            is $hsp->length('hsp'), 291;
 
64
            is $hsp->start('hit'), $hsp->hit->start;
 
65
            is $hsp->end('query'), $hsp->query->end;
 
66
            is $hsp->strand('sbjct'), $hsp->subject->strand;# alias for hit
 
67
            float_is($hsp->evalue, 6e-59);
 
68
            is $hsp->score, 119;
 
69
            is $hsp->bits,236;              
 
70
            is sprintf("%.1f",$hsp->percent_identity), 85.2;
 
71
            is sprintf("%.3f",$hsp->frac_identical('query')), 0.852;
 
72
            is sprintf("%.3f",$hsp->frac_identical('hit')), 0.852;
 
73
            is $hsp->gaps, 0;
 
74
            $hsps_left--;
 
75
        }
 
76
        is $hsps_left, 0;
 
77
    }
 
78
    last if( $count++ > @valid );
 
79
}
 
80
is @valid, 0;
 
81
 
 
82
# descriptionless hit
 
83
$searchio = Bio::SearchIO->new('-format' => 'blast_pull',
 
84
                                                           '-file' => test_input_file('blast_no_hit_desc.txt'));
 
85
 
 
86
$result = $searchio->next_result;
 
87
my $hit = $result->next_hit;
 
88
is $hit->name, 'chr1';
 
89
is $hit->description, '';
 
90
 
 
91
# further (NCBI blastn/p) tests taken from SearchIO.t
 
92
$searchio = Bio::SearchIO->new('-format' => 'blast_pull',
 
93
                                                           '-file' => test_input_file('ecolitst.bls'));
 
94
 
 
95
$result = $searchio->next_result;
 
96
 
 
97
is($result->database_name, 'ecoli.aa', 'database_name()');
 
98
is($result->database_entries, 4289);
 
99
is($result->database_letters, 1358990);
 
100
 
 
101
is($result->algorithm, 'BLASTP');
 
102
like($result->algorithm_version, qr/^2\.1\.3/);
 
103
like($result->query_name, qr/gi|1786183|gb|AAC73113.1| (AE000111) aspartokinase I,\s+homoserine dehydrogenase I [Escherichia coli]/);
 
104
is($result->query_length, 820);
 
105
is($result->get_statistic('kappa'), '0.135');
 
106
is($result->get_statistic('kappa_gapped'), '0.0410');
 
107
is($result->get_statistic('lambda'), '0.319');
 
108
is($result->get_statistic('lambda_gapped'), '0.267');
 
109
is($result->get_statistic('entropy'), '0.383');
 
110
is($result->get_statistic('entropy_gapped'), '0.140');
 
111
 
 
112
is($result->get_statistic('dbletters'), 1358990);
 
113
is($result->get_statistic('dbentries'), 4289);
 
114
is($result->get_statistic('effective_hsplength'), 47);
 
115
is($result->get_statistic('effectivespace'), 894675611);
 
116
is($result->get_parameter('matrix'), 'BLOSUM62');
 
117
is($result->get_parameter('gapopen'), 11);
 
118
is($result->get_parameter('gapext'), 1);
 
119
is($result->get_statistic('S2'), '92');
 
120
is($result->get_statistic('S2_bits'), '40.0');
 
121
float_is($result->get_parameter('expect'), '1.0e-03');
 
122
is($result->get_statistic('num_extensions'), '82424');
 
123
 
 
124
 
 
125
@valid = ( [ 'gb|AAC73113.1|', 820, 'AAC73113', '0', 1567],
 
126
              [ 'gb|AAC76922.1|', 810, 'AAC76922', '1e-91', 332],
 
127
              [ 'gb|AAC76994.1|', 449, 'AAC76994', '3e-47', 184]);
 
128
$count = 0;
 
129
while (my $hit = $result->next_hit) {
 
130
    my $d = shift @valid;
 
131
 
 
132
    is($hit->name, shift @$d);
 
133
    is($hit->length, shift @$d);
 
134
    is($hit->accession, shift @$d);
 
135
    float_is($hit->significance, shift @$d);
 
136
    is($hit->raw_score, shift @$d );
 
137
 
 
138
    if( $count == 0 ) {
 
139
        my $hsps_left = 1;
 
140
        while (my $hsp = $hit->next_hsp) {
 
141
            is($hsp->query->start, 1);
 
142
            is($hsp->query->end, 820);
 
143
            is($hsp->hit->start, 1);
 
144
            is($hsp->hit->end, 820);
 
145
            is($hsp->length('hsp'), 820);
 
146
            is($hsp->start('hit'), $hsp->hit->start);
 
147
            is($hsp->end('query'), $hsp->query->end);
 
148
            is($hsp->strand('sbjct'), $hsp->subject->strand);# alias for hit
 
149
            float_is($hsp->evalue, 0.0);
 
150
            is($hsp->score, 4058);
 
151
            is($hsp->bits,1567);                    
 
152
            is(sprintf("%.2f",$hsp->percent_identity), 98.29);
 
153
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.9829);
 
154
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.9829);
 
155
            is($hsp->gaps, 0);
 
156
            $hsps_left--;
 
157
        }
 
158
        is($hsps_left, 0);
 
159
    }
 
160
    last if( $count++ > @valid );
 
161
}
 
162
is(@valid, 0);
 
163
 
 
164
$searchio = Bio::SearchIO->new('-format' => 'blast_pull',
 
165
                                                          '-file'   => test_input_file('a_thaliana.blastn'));
 
166
 
 
167
$result = $searchio->next_result;
 
168
is($result->database_name, 'All GenBank+EMBL+DDBJ+PDB sequences (but no EST, STS, GSS,or phase 0, 1 or 2 HTGS sequences)');
 
169
is($result->database_letters, 4677375331);
 
170
is($result->database_entries, 1083200);
 
171
is($result->algorithm, 'BLASTN');
 
172
like($result->algorithm_version, qr/^2\.2\.1/);
 
173
is($result->query_name, '');
 
174
is($result->query_length, 60);
 
175
is($result->get_parameter('gapopen'), 5);
 
176
is($result->get_parameter('gapext'), 2);
 
177
is($result->get_parameter('ktup'), undef);
 
178
 
 
179
is($result->get_statistic('lambda'), 1.37);
 
180
is($result->get_statistic('kappa'), 0.711);
 
181
is($result->get_statistic('entropy'),1.31 );
 
182
is($result->get_statistic('T'), 0);
 
183
is($result->get_statistic('A'), 30);
 
184
is($result->get_statistic('X1'), '6');
 
185
is($result->get_statistic('X1_bits'), 11.9);
 
186
is($result->get_statistic('X2'), 15);
 
187
is($result->get_statistic('X2_bits'), 29.7);
 
188
is($result->get_statistic('S1'), 12);
 
189
is($result->get_statistic('S1_bits'), 24.3);
 
190
is($result->get_statistic('S2'), 17);
 
191
is($result->get_statistic('S2_bits'), 34.2);
 
192
 
 
193
is($result->get_statistic('dbentries'), 1083200);
 
194
 
 
195
@valid = ( [ 'gb|AY052359.1|', 2826, 'AY052359', '3e-18', 96, 1, 60, 
 
196
             '1.0000'],
 
197
           [ 'gb|AC002329.2|AC002329', 76170, 'AC002329', '3e-18', 96, 1, 60, 
 
198
             '1.0000' ],
 
199
           [ 'gb|AF132318.1|AF132318', 5383, 'AF132318', '0.04', 42, 35, 55, 
 
200
             '0.3500']);
 
201
$count = 0;
 
202
 
 
203
while( my $hit = $result->next_hit ) {
 
204
    my $d = shift @valid;
 
205
    is($hit->name, shift @$d);
 
206
    is($hit->length, shift @$d);
 
207
    is($hit->accession, shift @$d);
 
208
    float_is($hit->significance, shift @$d);
 
209
    is($hit->raw_score, shift @$d );
 
210
    is($hit->start, shift @$d);
 
211
    is($hit->end,shift @$d);    
 
212
    is(sprintf("%.4f",$hit->frac_aligned_query), shift @$d);
 
213
    if( $count == 0 ) {
 
214
        my $hsps_left = 1;
 
215
        while( my $hsp = $hit->next_hsp ) {
 
216
            is($hsp->query->start, 1);
 
217
            is($hsp->query->end, 60);
 
218
            is($hsp->query->strand, 1);
 
219
            is($hsp->hit->start, 154);
 
220
            is($hsp->hit->end, 212);
 
221
            is($hsp->hit->strand, 1);
 
222
            is($hsp->length('hsp'), 60);            
 
223
            float_is($hsp->evalue, 3e-18);
 
224
            is($hsp->score, 48);
 
225
            is($hsp->bits,95.6);
 
226
            is(sprintf("%.2f",$hsp->percent_identity), 96.67);
 
227
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.9667);
 
228
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.9831);
 
229
            is($hsp->query->frame(), 0);
 
230
            is($hsp->hit->frame(), 0);
 
231
            is($hsp->query->seq_id, undef);
 
232
            is($hsp->hit->seq_id, 'gb|AY052359.1|');
 
233
            is($hsp->gaps('query'), 0);
 
234
            is($hsp->gaps('hit'), 1);
 
235
            is($hsp->gaps, 1);
 
236
            is($hsp->query_string, 'aggaatgctgtttaattggaatcgtacaatggagaatttgacggaaatagaatcaacgat');
 
237
            is($hsp->hit_string, 'aggaatgctgtttaattggaatca-acaatggagaatttgacggaaatagaatcaacgat');
 
238
            is($hsp->homology_string, '|||||||||||||||||||||||  |||||||||||||||||||||||||||||||||||');
 
239
            is(sprintf("%.2f",$hsp->get_aln->overall_percentage_identity), 96.67);
 
240
            is(sprintf("%.2f",$hsp->get_aln->percentage_identity), 98.31);
 
241
            $hsps_left--;
 
242
        }
 
243
        is($hsps_left, 0);
 
244
    }
 
245
    last if( $count++ > @valid );
 
246
} 
 
247
is(@valid, 0);
 
248
 
 
249
$searchio = Bio::SearchIO->new('-format' => 'blast_pull',
 
250
                                                          '-file'   => test_input_file('frac_problems.blast'));
 
251
my @expected = ("1.000", "0.943");
 
252
while (my $result = $searchio->next_result) {
 
253
    my $hit = $result->next_hit;
 
254
        if (@expected == 2) {
 
255
                is($hit->frac_identical, shift @expected);
 
256
        }
 
257
        else {
 
258
                TODO: {
 
259
                        local $TODO = 'frac_identical failing!';
 
260
                        is($hit->frac_identical, shift @expected);
 
261
                }
 
262
        }
 
263
}
 
264
is(@expected, 0);
 
265
 
 
266
# And even more: frac_aligned_query should never be over 1!
 
267
$searchio = Bio::SearchIO->new('-format' => 'blast_pull',
 
268
                                                          '-file'   => test_input_file('frac_problems2.blast'));
 
269
$result = $searchio->next_result;
 
270
$hit = $result->next_hit;
 
271
is $hit->frac_aligned_query, 0.97;
 
272
 
 
273
# Also, start and end should be sane
 
274
$searchio = Bio::SearchIO->new('-format' => 'blast_pull',
 
275
                                                          '-file'   => test_input_file('frac_problems3.blast'));
 
276
$result = $searchio->next_result;
 
277
$hit = $result->next_hit;
 
278
is $hit->start('sbjct'), 207;
 
279
is $hit->end('sbjct'), 1051;
 
280
 
 
281
# Do a multiblast report test
 
282
$searchio = Bio::SearchIO->new('-format' => 'blast_pull',
 
283
                                                           '-file'   => test_input_file('multi_blast.bls'));
 
284
 
 
285
@expected = qw(CATH_RAT CATL_HUMAN CATL_RAT PAPA_CARPA);
 
286
my $results_left = 4;
 
287
while( my $result = $searchio->next_result ) {
 
288
    is($result->query_name, shift @expected, "Multiblast query test");
 
289
    $results_left--;
 
290
}
 
291
is($results_left, 0);
 
292
 
 
293
# Web blast result parsing
 
294
$searchio = Bio::SearchIO->new(-format => 'blast_pull',
 
295
                                                           -file   => test_input_file('catalase-webblast.BLASTP'));
 
296
ok($result = $searchio->next_result);
 
297
ok($hit = $result->next_hit);
 
298
is($hit->name, 'gi|40747822|gb|EAA66978.1|', 'full hit name');
 
299
is($hit->accession, 'EAA66978', 'hit accession');
 
300
ok(my $hsp = $hit->next_hsp);
 
301
is($hsp->query->start, 1, 'query start');
 
302
is($hsp->query->end, 528, 'query start');
 
303
 
 
304
# tests for new BLAST 2.2.13 output
 
305
$searchio = Bio::SearchIO->new(-format => 'blast_pull',
 
306
                                                          -file   => test_input_file('new_blastn.txt'));
 
307
 
 
308
$result = $searchio->next_result;
 
309
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)');
 
310
is($result->database_entries, 3742891);
 
311
is($result->database_letters, 16670205594);
 
312
is($result->algorithm, 'BLASTN');
 
313
is($result->algorithm_version, '2.2.13 [Nov-27-2005]');
 
314
is($result->query_name, 'pyrR,');
 
315
is($result->query_length, 558);
 
316
is($result->get_statistic('kappa'), '0.711');
 
317
is($result->get_statistic('kappa_gapped'), '0.711');
 
318
is($result->get_statistic('lambda'), '1.37');
 
319
is($result->get_statistic('lambda_gapped'), '1.37');
 
320
is($result->get_statistic('entropy'), '1.31');
 
321
is($result->get_statistic('entropy_gapped'), '1.31');
 
322
is($result->get_statistic('dbletters'), '-509663586');
 
323
is($result->get_statistic('dbentries'), 3742891);
 
324
is($result->get_statistic('effective_hsplength'), undef);
 
325
is($result->get_statistic('effectivespace'), 8935230198384);
 
326
is($result->get_parameter('matrix'), 'blastn matrix:1 -3');
 
327
is($result->get_parameter('gapopen'), 5);
 
328
is($result->get_parameter('gapext'), 2);
 
329
is($result->get_statistic('S2'), '60');
 
330
is($result->get_statistic('S2_bits'), '119.4');
 
331
float_is($result->get_parameter('expect'), '1e-23');
 
332
is($result->get_statistic('num_extensions'), '117843');
 
333
 
 
334
@valid = ( [ 'gi|41400296|gb|AE016958.1|', 4829781, 'AE016958', '6e-059', 236],
 
335
              [ 'gi|54013472|dbj|AP006618.1|', 6021225, 'AP006618', '4e-026', 127],
 
336
              [ 'gi|57546753|dbj|BA000030.2|', 9025608, 'BA000030', '1e-023', 119]);
 
337
$count = 0;
 
338
while( $hit = $result->next_hit ) {
 
339
    my $d = shift @valid;
 
340
 
 
341
    is($hit->name, shift @$d);
 
342
    is($hit->length, shift @$d);
 
343
    is($hit->accession, shift @$d);
 
344
    float_is($hit->significance, shift @$d);
 
345
    is($hit->raw_score, shift @$d );
 
346
 
 
347
    if( $count == 0 ) {
 
348
        my $hsps_left = 1;
 
349
        while( my $hsp = $hit->next_hsp ) {
 
350
            is($hsp->query->start, 262);
 
351
            is($hsp->query->end, 552);
 
352
            is($hsp->hit->start, 1166897);
 
353
            is($hsp->hit->end, 1167187);
 
354
            is($hsp->length('hsp'), 291);
 
355
            is($hsp->start('hit'), $hsp->hit->start);
 
356
            is($hsp->end('query'), $hsp->query->end);
 
357
            is($hsp->strand('sbjct'), $hsp->subject->strand);# alias for hit
 
358
            float_is($hsp->evalue, 6e-59);
 
359
            is($hsp->score, 119);
 
360
            is($hsp->bits,236);             
 
361
            is(sprintf("%.2f",$hsp->percent_identity), 85.22);
 
362
            is(sprintf("%.4f",$hsp->frac_identical('query')), 0.8522);
 
363
            is(sprintf("%.4f",$hsp->frac_identical('hit')), 0.8522);
 
364
            is($hsp->gaps, 0);
 
365
            $hsps_left--;
 
366
        }
 
367
        is($hsps_left, 0);
 
368
    }
 
369
    last if( $count++ > @valid );
 
370
}
 
371
is(@valid, 0);
 
372
 
 
373
# Bug 2189
 
374
$searchio = Bio::SearchIO->new(-format => 'blast_pull',
 
375
                                                          -file   => test_input_file('blastp2215.blast'));
 
376
 
 
377
$result = $searchio->next_result;
 
378
is($result->database_entries, 4460989);
 
379
is($result->database_letters, 1533424333);
 
380
is($result->algorithm, 'BLASTP');
 
381
is($result->algorithm_version, '2.2.15 [Oct-15-2006]');
 
382
is($result->query_name, 'gi|15608519|ref|NP_215895.1|');
 
383
is($result->query_length, 193);
 
384
my @hits = $result->hits;
 
385
is(scalar(@hits), 10);
 
386
is($hits[1]->accession,'1W30');
 
387
float_is($hits[4]->significance,'2e-72');
 
388
is($hits[7]->score,'254');
 
389
$result = $searchio->next_result;
 
390
is($result->database_entries, 4460989);
 
391
is($result->database_letters, 1533424333);
 
392
is($result->algorithm, 'BLASTP');
 
393
is($result->algorithm_version, '2.2.15 [Oct-15-2006]');
 
394
is($result->query_name, 'gi|15595598|ref|NP_249092.1|');
 
395
is($result->query_length, 423);
 
396
@hits = $result->hits;
 
397
is(scalar(@hits), 10);
 
398
is($hits[1]->accession,'ZP_00972546');
 
399
float_is($hits[4]->significance, '0.0');
 
400
is($hits[7]->score, 624);
 
401
 
 
402
# Bug 2246
 
403
$searchio = Bio::SearchIO->new(-format => 'blast_pull',
 
404
                               -verbose => -1,
 
405
                                                          -file   => test_input_file('bug2246.blast'));
 
406
$result = $searchio->next_result;
 
407
$hit = $result->next_hit;
 
408
is $hit->name, 'UniRef50_Q9X0H5';
 
409
is $hit->length, undef;
 
410
is $hit->accession, 'UniRef50_Q9X0H5';
 
411
is $hit->description, 'Cluster: Histidyl-tRNA synthetase; n=4; Thermoto...';
 
412
is $hit->raw_score, 23;
 
413
float_is($hit->significance, 650);
 
414
 
 
415
#*** need to /fully/ test a multi-result, multi-hit, multi-hsp report!