← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Bio/Tools/BPlite.pm

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Tools/BPlite.pm
1
 
# $Id: BPlite.pm,v 1.42.4.1 2006/10/02 23:10:31 sendu Exp $
2
 
##############################################################################
3
 
# Bioperl module Bio::Tools::BPlite
4
 
##############################################################################
5
 
#
6
 
# The original BPlite.pm module has been written by Ian Korf !
7
 
# see http://sapiens.wustl.edu/~ikorf
8
 
#
9
 
# You may distribute this module under the same terms as perl itself
10
 
 
11
 
=head1 NAME
12
 
 
13
 
Bio::Tools::BPlite - Lightweight BLAST parser
14
 
 
15
 
=head1 SYNOPSIS
16
 
 
17
 
 use Bio::Tools::BPlite;
18
 
 my $report = new Bio::Tools::BPlite(-fh=>\*STDIN);
19
 
 
20
 
  {
21
 
    $report->query;
22
 
    $report->database;
23
 
    while(my $sbjct = $report->nextSbjct) {
24
 
        $sbjct->name;
25
 
        while (my $hsp = $sbjct->nextHSP) {
26
 
            $hsp->score;
27
 
            $hsp->bits;
28
 
            $hsp->percent;
29
 
            $hsp->P;
30
 
            $hsp->EXP;
31
 
            $hsp->match;
32
 
            $hsp->positive;
33
 
            $hsp->length;
34
 
            $hsp->querySeq;
35
 
            $hsp->sbjctSeq;
36
 
            $hsp->homologySeq;
37
 
            $hsp->query->start;
38
 
            $hsp->query->end;
39
 
            $hsp->hit->start;
40
 
            $hsp->hit->end;
41
 
            $hsp->hit->seq_id;
42
 
            $hsp->hit->overlaps($exon);
43
 
        }
44
 
    }
45
 
 
46
 
    # the following line takes you to the next report in the stream/file
47
 
    # it will return 0 if that report is empty,
48
 
    # but that is valid for an empty blast report.
49
 
    # Returns -1 for EOF.
50
 
 
51
 
    last if ($report->_parseHeader == -1);
52
 
    redo;
53
 
  }
54
 
 
55
 
 
56
 
=head1 DESCRIPTION
57
 
 
58
 
BPlite is a package for parsing BLAST reports. The BLAST programs are a family
59
 
of widely used algorithms for sequence database searches. The reports are
60
 
non-trivial to parse, and there are differences in the formats of the various
61
 
flavors of BLAST. BPlite parses BLASTN, BLASTP, BLASTX, TBLASTN, and TBLASTX
62
 
reports from both the high performance WU-BLAST, and the more generic
63
 
NCBI-BLAST.
64
 
 
65
 
Many people have developed BLAST parsers (I myself have made at least three).
66
 
BPlite is for those people who would rather not have a giant object
67
 
specification, but rather a simple handle to a BLAST report that works well
68
 
in pipes.
69
 
 
70
 
=head2 Object
71
 
 
72
 
BPlite has three kinds of objects, the report, the subject, and the HSP. To
73
 
create a new report, you pass a filehandle reference to the BPlite constructor.
74
 
 
75
 
 my $report = new Bio::Tools::BPlite(-fh=>\*STDIN); # or any other filehandle
76
 
 
77
 
The report has two attributes (query and database), and one method (nextSbjct).
78
 
 
79
 
 $report->query;     # access to the query name
80
 
 $report->database;  # access to the database name
81
 
 $report->nextSbjct; # gets the next subject
82
 
 while(my $sbjct = $report->nextSbjct) {
83
 
     # canonical form of use is in a while loop
84
 
 }
85
 
 
86
 
A subject is a BLAST hit, which should not be confused with an HSP (below). A
87
 
BLAST hit may have several alignments associated with it. A useful way of
88
 
thinking about it is that a subject is a gene and HSPs are the exons. Subjects
89
 
have one attribute (name) and one method (nextHSP).
90
 
 
91
 
 $sbjct->name;    # access to the subject name
92
 
 $sbjct->nextHSP; # gets the next HSP from the sbjct
93
 
 while(my $hsp = $sbjct->nextHSP) {
94
 
     # canonical form is again a while loop
95
 
 }
96
 
 
97
 
An HSP is a high scoring pair, or simply an alignment.  HSP objects
98
 
inherit all the useful methods from RangeI/SeqFeatureI/FeaturePair,
99
 
but provide an additional set of attributes (score, bits, percent, P,
100
 
match, EXP, positive, length, querySeq, sbjctSeq, homologySeq) that
101
 
should be familiar to anyone who has seen a blast report.
102
 
 
103
 
For lazy/efficient coders, two-letter abbreviations are available for the 
104
 
attributes with long names (qs, ss, hs). Ranges of the aligned sequences in
105
 
query/subject and other information (like seqname) are stored
106
 
in SeqFeature objects (i.e.: $hsp-E<gt>query, $hsp-E<gt>subject which is equal to
107
 
$hsp-E<gt>feature1, $hsp-E<gt>feature2). querySeq, sbjctSeq and homologySeq do only
108
 
contain the alignment sequences from the blast report.
109
 
 
110
 
 $hsp->score;
111
 
 $hsp->bits;
112
 
 $hsp->percent;
113
 
 $hsp->P;
114
 
 $hsp->match;
115
 
 $hsp->positive;
116
 
 $hsp->length;
117
 
 $hsp->querySeq;      $hsp->qs;
118
 
 $hsp->sbjctSeq;      $hsp->ss;
119
 
 $hsp->homologySeq;   $hsp->hs;
120
 
 $hsp->query->start;
121
 
 $hsp->query->end;
122
 
 $hsp->query->seq_id;
123
 
 $hsp->hit->primary_tag; # "similarity"
124
 
 $hsp->hit->source_tag;  # "BLAST"
125
 
 $hsp->hit->start;
126
 
 $hsp->hit->end;
127
 
 ...
128
 
 
129
 
So a very simple look into a BLAST report might look like this.
130
 
 
131
 
 my $report = new Bio::Tools::BPlite(-fh=>\*STDIN);
132
 
 while(my $sbjct = $report->nextSbjct) {
133
 
     print ">",$sbjct->name,"\n";
134
 
     while(my $hsp = $sbjct->nextHSP) {
135
 
                print "\t",$hsp->start,"..",$hsp->end," ",$hsp->bits,"\n";
136
 
     }
137
 
 }
138
 
 
139
 
The output of such code might look like this:
140
 
 
141
 
 >foo
142
 
     100..155 29.5
143
 
     268..300 20.1
144
 
 >bar
145
 
     100..153 28.5
146
 
     265..290 22.1
147
 
 
148
 
 
149
 
=head1 AUTHORS
150
 
 
151
 
Ian Korf (ikorf@sapiens.wustl.edu, http://sapiens.wustl.edu/~ikorf), 
152
 
Lorenz Pollak (lorenz@ist.org, bioperl port)
153
 
 
154
 
=head1 ACKNOWLEDGEMENTS
155
 
 
156
 
This software was developed at the Genome Sequencing Center at Washington
157
 
Univeristy, St. Louis, MO.
158
 
 
159
 
=head1 CONTRIBUTORS
160
 
 
161
 
Jason Stajich, jason@cgt.mc.duke.edu
162
 
 
163
 
=head1 COPYRIGHT
164
 
 
165
 
Copyright (C) 1999 Ian Korf. All Rights Reserved.
166
 
 
167
 
=head1 DISCLAIMER
168
 
 
169
 
This software is provided "as is" without warranty of any kind.
170
 
 
171
 
=cut
172
 
 
173
 
package Bio::Tools::BPlite;
174
 
 
175
 
use strict;
176
 
 
177
 
use Bio::Tools::BPlite::Sbjct; # we want to use Sbjct
178
 
use Symbol;
179
 
 
180
 
use base qw(Bio::Root::Root Bio::SeqAnalysisParserI Bio::Root::IO);
181
 
 
182
 
# new comes from a RootI now
183
 
 
184
 
=head2 new
185
 
 
186
 
 Title   : new
187
 
 Function: Create a new Bio::Tools::BPlite object
188
 
 Returns : Bio::Tools::BPlite
189
 
 Args    : -file     input file (alternative to -fh)
190
 
           -fh       input stream (alternative to -file)
191
 
 
192
 
=cut
193
 
 
194
 
sub new {
195
 
  my ($class, @args) = @_; 
196
 
  my $self = $class->SUPER::new(@args);
197
 
    $self->warn("Use of Bio::Tools::BPlite is deprecated".
198
 
                   "Use Bio::SearchIO classes instead");
199
 
  # initialize IO
200
 
  $self->_initialize_io(@args);
201
 
 
202
 
  $self->{'QPATLOCATION'} = [];  # Anonymous array of query pattern locations for PHIBLAST
203
 
 
204
 
  if ($self->_parseHeader) {$self->{'REPORT_DONE'} = 0} # there are alignments
205
 
  else                     {$self->{'REPORT_DONE'} = 1} # empty report
206
 
  
207
 
  return $self; # success - we hope!
208
 
}
209
 
 
210
 
# for SeqAnalysisParserI compliance
211
 
 
212
 
=head2 next_feature
213
 
 
214
 
 Title   : next_feature
215
 
 Usage   : while( my $feat = $res->next_feature ) { # do something }
216
 
 Function: SeqAnalysisParserI implementing function. This implementation
217
 
           iterates over all HSPs. If the HSPs of the current subject match
218
 
           are exhausted, it will automatically call nextSbjct().
219
 
 Example :
220
 
 Returns : A Bio::SeqFeatureI compliant object, in this case a
221
 
           Bio::Tools::BPlite::HSP object, and FALSE if there are no more
222
 
           HSPs.
223
 
 Args    : None
224
 
 
225
 
=cut
226
 
 
227
 
sub next_feature{
228
 
   my ($self) = @_;
229
 
   my ($sbjct, $hsp);
230
 
   $sbjct = $self->{'_current_sbjct'};
231
 
   unless( defined $sbjct ) {
232
 
       $sbjct = $self->{'_current_sbjct'} = $self->nextSbjct;
233
 
       return  unless defined $sbjct;
234
 
   }   
235
 
   $hsp = $sbjct->nextHSP;
236
 
   unless( defined $hsp ) {
237
 
       $self->{'_current_sbjct'} = undef;
238
 
       return $self->next_feature;
239
 
   }
240
 
   return $hsp || undef;
241
 
}
242
 
 
243
 
=head2 query
244
 
 
245
 
 Title    : query
246
 
 Usage    : $query = $obj->query();
247
 
 Function : returns the query object
248
 
 Example  :
249
 
 Returns  : query object
250
 
 Args     :
251
 
 
252
 
=cut
253
 
 
254
 
sub query    {shift->{'QUERY'}}
255
 
 
256
 
=head2 qlength
257
 
 
258
 
 Title    : qlength
259
 
 Usage    : $len = $obj->qlength();
260
 
 Function : returns the length of the query 
261
 
 Example  :
262
 
 Returns  : length of query
263
 
 Args     :
264
 
 
265
 
=cut
266
 
 
267
 
sub qlength  {shift->{'LENGTH'}}
268
 
 
269
 
=head2 pattern
270
 
 
271
 
 Title    : pattern
272
 
 Usage    : $pattern = $obj->pattern();
273
 
 Function : returns the pattern used in a PHIBLAST search
274
 
 
275
 
=cut
276
 
 
277
 
sub pattern {shift->{'PATTERN'}}
278
 
 
279
 
=head2 query_pattern_location
280
 
 
281
 
 Title    : query_pattern_location
282
 
 Usage    : $qpl = $obj->query_pattern_location();
283
 
 Function : returns reference to array of locations in the query sequence
284
 
            of pattern used in a PHIBLAST search
285
 
 
286
 
=cut
287
 
 
288
 
sub query_pattern_location {shift->{'QPATLOCATION'}}
289
 
 
290
 
=head2 database
291
 
 
292
 
 Title    : database
293
 
 Usage    : $db = $obj->database();
294
 
 Function : returns the database used in this search
295
 
 Example  :
296
 
 Returns  : database used for search
297
 
 Args     :
298
 
 
299
 
=cut
300
 
 
301
 
sub database {shift->{'DATABASE'}}
302
 
 
303
 
=head2 nextSbjct
304
 
 
305
 
 Title    : nextSbjct
306
 
 Usage    : $sbjct = $obj->nextSbjct();
307
 
 Function : Method of iterating through all the Sbjct retrieved 
308
 
            from parsing the report 
309
 
 Example  : while ( my $sbjct = $obj->nextSbjct ) {}
310
 
 Returns  : next Sbjct object or null if finished
311
 
 Args     :
312
 
 
313
 
=cut
314
 
 
315
 
sub nextSbjct {
316
 
  my ($self) = @_;
317
 
  
318
 
  $self->_fastForward or return;
319
 
  local $_;
320
 
  #######################
321
 
  # get all sbjct lines #
322
 
  #######################
323
 
  my $def = $self->_readline();  
324
 
  while(defined ($_ = $self->_readline() ) ) {
325
 
    if    (! /\w/)           {next}
326
 
    elsif (/Strand HSP/o)    {next} # WU-BLAST non-data
327
 
    elsif (/^\s{0,2}Score/o) {$self->_pushback($_); last}
328
 
    elsif (/^Histogram|^Searching|^Parameters|
329
 
            ^\s+Database:|
330
 
            ^\s+Posted date:/ox) {
331
 
        $self->_pushback($_); 
332
 
        last;
333
 
    } else {
334
 
        $def .= $_;
335
 
    }
336
 
  }
337
 
  if( ! $def ) { 
338
 
      return;
339
 
  }
340
 
  $def =~ s/\s+/ /g;
341
 
  $def =~ s/\s+$//g;
342
 
  
343
 
  my $length;
344
 
  if( $def =~ s/Length = ([\d,]+)$//g ) {
345
 
      $length = $1;
346
 
  }
347
 
  return unless $def =~ /^>/;
348
 
  $def =~ s/^>//;
349
 
 
350
 
  ####################
351
 
  # the Sbjct object #
352
 
  ####################
353
 
  my $sbjct = new Bio::Tools::BPlite::Sbjct('-name'=>$def,
354
 
                                            '-length'=>$length,
355
 
                                            '-parent'=>$self);
356
 
  return $sbjct;
357
 
}
358
 
 
359
 
# begin private routines
360
 
 
361
 
sub _parseHeader {
362
 
  my ($self) = @_;
363
 
 
364
 
  # normally, _parseHeader will break out of the parse as soon as it
365
 
  # reaches a new Subject (i.e. the first one after the header) if you
366
 
  # call _parseHeader twice in a row, with nothing in between, all you
367
 
  # accomplish is a ->nextSubject call..  so we need a flag to
368
 
  # indicate that we have *entered* a header, before we are allowed to
369
 
  # leave it!
370
 
 
371
 
  my $header_flag = 0; # here is the flag/ It is "false" at first, and
372
 
                       # is set to "true" when any valid header element
373
 
                       # is encountered
374
 
  local $_;
375
 
  $self->{'REPORT_DONE'} = 0;  # reset this bit for a new report
376
 
  while(defined($_ = $self->_readline() ) ) {
377
 
      s/\(\s*\)//;      
378
 
      if (/^Query=(?:\s+(.+))?/s) {
379
 
          $header_flag = 1;     # valid header element found
380
 
          my $query = $1;
381
 
          while( defined($_ = $self->_readline() ) ) {
382
 
              # Continue reading query name until encountering either
383
 
              # a line that starts with "Database" or a blank line.
384
 
              # The latter condition is needed in order to be able to
385
 
              # parse megablast output correctly, since Database comes
386
 
              # before (not after) the query.
387
 
              if( ($_ =~ /^Database/) || ($_ =~ /^$/) ) {
388
 
                  $self->_pushback($_); last;
389
 
              }       
390
 
              $query .= $_;
391
 
          }
392
 
          $query =~ s/\s+/ /g;
393
 
          $query =~ s/\s+$//;
394
 
          $query =~ s/^>//;
395
 
 
396
 
          my $length = 0;
397
 
          if( $query =~ /\(([\d,]+)\s+\S+\)\s*$/ ) {      
398
 
              $length = $1;
399
 
              $length =~ s/,//g;
400
 
          } else { 
401
 
              $self->debug("length is 0 for '$query'\n");
402
 
          }
403
 
          $self->{'QUERY'} = $query;
404
 
          $self->{'LENGTH'} = $length;
405
 
      }
406
 
      elsif (/^(<b>)?(T?BLAST[NPX])\s+([\w\.-]+)\s+(\[[\w-]*\])/o) { 
407
 
          $self->{'BLAST_TYPE'} = $2; 
408
 
          $self->{'BLAST_VERSION'} = $3;
409
 
      }                         # BLAST report type - not a valid header element # JB949
410
 
      
411
 
      # Support Paracel BTK output
412
 
      elsif ( $_ =~ /(^[A-Z0-9_]+)\s+BTK\s+/ ) { 
413
 
          $self->{'BLAST_TYPE'} = $1;
414
 
          $self->{'BTK'} = 1;
415
 
      } 
416
 
      elsif ($_ =~ /^Database:\s+(.+)/) {$header_flag = 1;$self->{'DATABASE'} = $1} # valid header element found
417
 
      elsif ($_ =~ /^\s*pattern\s+(\S+).*position\s+(\d+)\D/) {   
418
 
          # For PHIBLAST reports
419
 
          $header_flag = 1;     # valid header element found
420
 
          $self->{'PATTERN'} = $1;
421
 
          push (@{$self->{'QPATLOCATION'}}, $2);
422
 
      } 
423
 
      elsif (($_ =~ /^>/) && ($header_flag==1)) {$self->_pushback($_); return 1} # only leave if we have actually parsed a valid header!
424
 
      elsif (($_ =~ /^Parameters|^\s+Database:/) && ($header_flag==1)) { 
425
 
      # if we entered a header, and saw nothing before the stats at the end, 
426
 
      # then it was empty
427
 
          $self->_pushback($_);
428
 
          return 0;             # there's nothing in the report
429
 
      } elsif( /Reference:\s+Aaron E\. Darling/ ) {
430
 
          $self->{'BTK'} = 1;
431
 
      }  
432
 
      # bug fix suggested by MI Sadowski via Martin Lomas
433
 
      # see bug report #1118
434
 
      if( ref($self->_fh()) !~ /GLOB/ && 
435
 
          $self->_fh()->can('EOF') && eof($self->_fh()) ) {
436
 
          $self->warn("unexpected EOF in file\n");
437
 
          return -1;
438
 
      }
439
 
  }
440
 
  return -1; # EOF
441
 
}
442
 
 
443
 
sub _fastForward {
444
 
    my ($self) = @_;
445
 
    return 0 if $self->{'REPORT_DONE'}; # empty report
446
 
    local $_;
447
 
    while(defined( $_ = $self->_readline() ) ) {
448
 
        if (/^Histogram|^Searching|^Parameters|^\s+Database:|
449
 
             ^\s+Posted date:/xo) {
450
 
            return 0;
451
 
        } elsif( $self->{'BTK'} && /^BLAST/o ) {
452
 
            return 0;
453
 
        } elsif( /^>/ ) {
454
 
            $self->_pushback($_);       
455
 
            return 1;
456
 
        }
457
 
    }
458
 
    unless( $self->{'BTK'} ) { # Paracel BTK reports have no footer
459
 
        $self->warn("Possible error (1) while parsing BLAST report!");
460
 
    }
461
 
}
462
 
 
463
 
1;
464
 
__END__