[Bioperl-l] Error parsing BLAST report

Freimuth, Robert freimuth at pathology.wustl.edu
Wed Aug 16 15:56:37 EDT 2006


Hello,

I'm trying to parse a BLAST report using the following code:

use warnings;
use strict;

use Bio::SearchIO;

my $file = 'NP_006065_blast.out';

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

while( my $result = $searchio->next_result() )
{
    while( my $hit = $result->next_hit )
    {
            my $hit_acc_num = $hit->accession();

            # get the total length of the aligned region for query or
sbjct seq
            # (includes all HSPs, calculated after tiling)

            my $align_len = $hit->length_aln( 'query' );

            print "Alignment length for $hit_acc_num is $align_len\n";
    }
}

There are 104 one-line descriptions in the report, and alignments for
each one of them (the blast report was created using
b_num_alignments_shown => 500 and v_num_descriptions_shown => 500).
However, when I run the above code I get 14 errors like the following:

-------------------- WARNING ---------------------
MSG: There is no HSP data for hit 'ENSP00000327738'.
You have called a method (Bio::Search::Hit::GenericHit::length_aln)
that requires HSP data and there was no HSP data for this hit,
most likely because it was absent from the BLAST report.
Note that by default, BLAST lists alignments for the first 250 hits,
but it lists descriptions for 500 hits. If this is the case,
and you care about these hits, you should re-run BLAST using the
-b option (or equivalent if not using blastall) to increase the number
of alignments.
 
---------------------------------------------------

There is an alignment for this (and the other 13 sequences) in the
report.  In fact, if I edit the report and delete all but the
description and the alignment for ENSP00000327738, it parses fine (no
error).

I continued editing the report and produced the following minimal test
case that reproduces the error.  Note that the description for
ENSP00000350182 appears twice, BUT THE ERROR IS FOR ENSP00000327738.

*********** BLAST REPORT FOR TEST CASE ***********

BLASTP 2.2.11 [Jun-05-2005]


Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A.
Schaffer, 
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), 
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Query= NP_006065
         (442 letters)

Database: Homo_sapiens.NCBI36.apr.pep.fa 
           48,851 sequences; 23,910,368 total letters

Searching..................................................done

                                                                 Score
E
Sequences producing significant alignments:                      (bits)
Value

ENSP00000350182 pep:novel clone::BX322644.8:4905:15090:-1 gene:E...
120   3e-27
ENSP00000350182 pep:novel clone::BX322644.8:4905:15090:-1 gene:E...
120   3e-27
ENSP00000327738 pep:known-ccds chromosome:NCBI36:4:189297592:189...
115   8e-26

>ENSP00000350182 pep:novel clone::BX322644.8:4905:15090:-1
gene:ENSG00000137397
           transcript:ENST00000357569 
          Length = 425

 Score =  120 bits (301), Expect = 3e-27
 Identities = 76/261 (29%), Positives = 140/261 (53%), Gaps = 21/261
(8%)

Query: 9   IEKEVTCPICLELLTEPLSLDCGHSFCQACITAKIKESVIISRGESSCPVCQTRFQPGNL
68
           +++EV CPICL++L +P+++DCGH+FC  CIT +I E+   S G   CP+C+T  +   +
Sbjct: 10  LQEEVICPICLDILQKPVTIDCGHNFCLKCIT-QIGET---SCGFFKCPLCKTSVRKNAI
65

Query: 69  RPNRHLANIVERVKEVKMSP-QEGQKRDVCEHHGKKLQIFCKEDGKVICWVCELSQEHQG
127
           R N  L N+VE+++ ++ S  Q  +K   C  H +    FC++DGK +C+VC  S++H+ 
Sbjct: 66  RFNSLLRNLVEKIQALQASEVQSKRKEATCPRHQEMFHYFCEDDGKFLCFVCCESKDHKS
125

Query: 128 HQTFRINEVVKECQEKLQVALQRLIKEDQEAEKLED------DIRQERTAWKIERQKILK
181
           H    I E  +  Q ++Q  +Q L ++++E  +++       D+  ++   + E+Q+IL 
Sbjct: 126 HNVSLIEEAAQNYQGQIQEQIQVLQQKEKETVQVKAQGVHRVDVFTDQV--EHEKQRILT
183

Query: 182 GFNEMRVILDNEEQRELQKL----EEGEVNVLDNLAAATDQLVQQRQDASTLISDLQRRL
237
            F  +  +L+ E+   L ++     EG       +A+   QL     D   L+  L+ + 
Sbjct: 184 EFELLHQVLEEEKNFLLSRIYWLGHEGTEAGKHYVASTEPQL----NDLKKLVDSLKTKQ
239

Query: 238 TGSSVEMLQDVIDVMKRSESW 258
                ++L+     + RSE +
Sbjct: 240 NMPPRQLLEVTQPHLPRSEEF 260


>ENSP00000327738 pep:known-ccds
chromosome:NCBI36:4:189297592:189305643:1
           gene:ENSG00000184108 transcript:ENST00000332517
           CCDS3851.1
          Length = 468

 Score =  115 bits (289), Expect = 8e-26
 Identities = 101/410 (24%), Positives = 180/410 (43%), Gaps = 39/410
(9%)

Query: 8   DIEKEVTCPICLELLTEPLSLDCGHSFCQACITAKIKESVIISRGESSCPVCQTRFQPGN
67
           ++ +E+TC ICL+  + P++ +CGHSFC  C+    +E         SCP C    +  +
Sbjct: 9   NLREELTCFICLDYFSSPVTTECGHSFCLVCLLRSWEE----HNTPLSCPECWRTLEGPH
64

Query: 68  LRPNRHLANIVERVKEVKMSPQEGQKRDVCEHHGK-----KLQIFCKEDGKVICWVCELS
122
            + N  L  +    ++++   Q  Q  D    +G+     K     ++ G        ++
Sbjct: 65  FQSNERLGRLASIARQLR--SQVLQSEDEQGSYGRMPTTAKALSDDEQGGSAF-----VA
117

Query: 123 QEHQGHQTFRINEVVKECQEKLQVALQRLIKEDQEA------EKLEDDIRQERTAWKIER
176
           Q H  ++    +E  +  +EKLQ  L  L    +EA      EK    + QE T  K  +
Sbjct: 118 QSHGANRVHLSSEAEEHHREKLQEILNLLRVRRKEAQAVLTHEKERVKLCQEET--KTCK
175

Query: 177 QKILKGFNEMRVILDNEEQRELQKLEEGEVNVLDNLAAATDQLVQQRQDASTLISDLQRR
236
           Q ++  + +M   L  EEQ +LQ LE+ E   +  L     +L QQ +  S +I+ ++  
Sbjct: 176 QVVVSEYMKMHQFLKEEEQLQLQLLEQEEKENMRKLRNNEIKLTQQIRSLSKMIAQIESS
235

Query: 237 LTGSSVEMLQDVIDVMKRSESWTXXXXXXXXXXXXXXFRVPDLSGMLQVLKELTDVQYYW
296
              S+ E L++V   ++RSE                   +  ++GM ++L++ +      
Sbjct: 236 SQSSAFESLEEVRGALERSE----PLLLQCPEATTTELSLCRITGMKEMLRKFS------
285

Query: 297 VDVMLNPGSATSNVAISVDQRQVKTVRTCTFKNSNPCDF-SAFGVFGCQYFSSGKYYWEV
355
            ++ L+P +A + + +S D + VK   +      NP  F  +  V G Q F+SG++YWEV
Sbjct: 286 TEITLDPATANAYLVLSEDLKSVKYGGSRQQLPDNPERFDQSATVLGTQIFTSGRHYWEV
345

Query: 356 DVSGKIAWILGVHSKISSLNKRKSSGFAFDPSVNYSKVYSRYRPQYGYWV 405
           +V  K  W +G+     S    +       P   +S +  +    Y  WV
Sbjct: 346 EVGNKTEWEVGICKDSVS----RKGNLPKPPGDLFSLIGLKIGDDYSLWV 391


  Database: Homo_sapiens.NCBI36.apr.pep.fa
    Posted date:  Jun 15, 2006  8:56 PM
  Number of letters in database: 23,910,368
  Number of sequences in database:  48,851
  
Lambda     K      H
   0.319    0.133    0.398 

Gapped
Lambda     K      H
   0.267   0.0410    0.140 


Matrix: BLOSUM62
Gap Penalties: Existence: 11, Extension: 1
Number of Hits to DB: 20,900,506
Number of Sequences: 48851
Number of extensions: 899179
Number of successful extensions: 6075
Number of sequences better than 1.0e-25: 105
Number of HSP's better than  0.0 without gapping: 18
Number of HSP's successfully gapped in prelim test: 87
Number of HSP's that attempted gapping in prelim test: 5632
Number of HSP's gapped (non-prelim): 157
length of query: 442
length of database: 23,910,368
effective HSP length: 107
effective length of query: 335
effective length of database: 18,683,311
effective search space: 6258909185
effective search space used: 6258909185
T: 11
A: 40
X1: 16 ( 7.4 bits)
X2: 38 (14.6 bits)
X3: 64 (24.7 bits)
S1: 41 (21.8 bits)
S2: 289 (115.9 bits)

*********** END BLAST REPORT FOR TEST CASE ***********

Any ideas?

Thanks,

Bob




More information about the Bioperl-l mailing list