← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/bioperl/oneiric

« back to all changes in this revision

Viewing changes to Bio/Tools/BPlite/HSP.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2009-03-10 07:19:11 UTC
  • mfrom: (3.1.2 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090310071911-ever3si2bbzx1iks
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/HSP.pm
1
 
###############################################################################
2
 
# Bio::Tools::BPlite::HSP
3
 
###############################################################################
4
 
# HSP = High Scoring Pair (to all non-experts as I am)
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
 
 
12
 
#
13
 
# BioPerl module for Bio::Tools::BPlite::HSP
14
 
#
15
 
# Cared for by Peter Schattner <schattner@alum.mit.edu>
16
 
#
17
 
# Copyright Peter Schattner
18
 
#
19
 
# You may distribute this module under the same terms as perl itself
20
 
 
21
 
# POD documentation - main docs before the code
22
 
 
23
 
=head1 NAME
24
 
 
25
 
Bio::Tools::BPlite::HSP - Blast report High Scoring Pair (HSP)
26
 
 
27
 
=head1 SYNOPSIS
28
 
 
29
 
 use Bio::Tools::BPlite;
30
 
 my $report = new Bio::Tools::BPlite(-fh=>\*STDIN);
31
 
 {
32
 
    while(my $sbjct = $report->nextSbjct) {
33
 
        while (my $hsp = $sbjct->nextHSP) {
34
 
            $hsp->score;
35
 
            $hsp->bits;
36
 
            $hsp->percent;
37
 
            $hsp->P;
38
 
            $hsp->match;
39
 
            $hsp->positive;
40
 
            $hsp->length;
41
 
            $hsp->querySeq;
42
 
            $hsp->sbjctSeq;
43
 
            $hsp->homologySeq;
44
 
            $hsp->query->start;
45
 
            $hsp->query->end;
46
 
            $hsp->hit->start;
47
 
            $hsp->hit->end;
48
 
            $hsp->hit->seq_id;
49
 
            $hsp->hit->overlaps($exon);
50
 
        }
51
 
    }
52
 
 
53
 
    # the following line takes you to the next report in the stream/file
54
 
    # it will return 0 if that report is empty,
55
 
    # but that is valid for an empty blast report.
56
 
    # Returns -1 for EOF.
57
 
 
58
 
    last if $report->_parseHeader == -1;
59
 
 
60
 
 redo
61
 
 }
62
 
 
63
 
=head1 DESCRIPTION
64
 
 
65
 
This object handles the High Scoring Pair data for a Blast report.
66
 
This is where the percent identity, query and hit sequence length,
67
 
P value, etc are stored and where most of the necessary information is located when building logic around parsing a Blast report.
68
 
 
69
 
See L<Bio::Tools::BPlite> for more detailed information on the entire
70
 
BPlite Blast parsing system.
71
 
 
72
 
=head1 FEEDBACK
73
 
 
74
 
=head2 Mailing Lists
75
 
 
76
 
User feedback is an integral part of the evolution of this and other
77
 
Bioperl modules. Send your comments and suggestions preferably to
78
 
the Bioperl mailing list.  Your participation is much appreciated.
79
 
 
80
 
  bioperl-l@bioperl.org                  - General discussion
81
 
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
82
 
 
83
 
=head2 Reporting Bugs
84
 
 
85
 
Report bugs to the Bioperl bug tracking system to help us keep track
86
 
of the bugs and their resolution. Bug reports can be submitted via
87
 
email the web:
88
 
 
89
 
  http://bugzilla.open-bio.org/
90
 
 
91
 
=head1 AUTHOR - Peter Schattner
92
 
 
93
 
Email: schattner@alum.mit.edu
94
 
 
95
 
=head1 CONTRIBUTORS
96
 
 
97
 
Jason Stajich, jason@bioperl.org
98
 
 
99
 
=head1 APPENDIX
100
 
 
101
 
The rest of the documentation details each of the object methods.
102
 
Internal methods are usually preceded with a _
103
 
 
104
 
=cut
105
 
 
106
 
# Let the code begin...
107
 
 
108
 
package Bio::Tools::BPlite::HSP;
109
 
 
110
 
use strict;
111
 
 
112
 
# to disable overloading comment this out:
113
 
#use overload '""' => '_overload';
114
 
 
115
 
# Object preamble - inheriets from Bio::SeqFeature::SimilarityPair
116
 
 
117
 
use Bio::SeqFeature::Similarity;
118
 
 
119
 
use base qw(Bio::SeqFeature::SimilarityPair);
120
 
 
121
 
sub new {
122
 
    my ($class, @args) = @_;
123
 
 
124
 
    # workaround to make sure frame is not set before strand is
125
 
    # interpreted from query/hit info 
126
 
    # this workaround removes the key from the hash
127
 
    # so the superclass does not try and work with it
128
 
    # we'll take care of setting it in this module later on
129
 
 
130
 
    my %newargs = @args;
131
 
    foreach ( keys %newargs ) {
132
 
        if( /frame$/i ) {
133
 
            delete $newargs{$_};
134
 
        } 
135
 
    }
136
 
    # done with workaround
137
 
 
138
 
    my $self = $class->SUPER::new(%newargs);
139
 
    
140
 
    my ($score,$bits,$match,$hsplength,$positive,$gaps,$p,$exp,$qb,$qe,$sb,
141
 
        $se,$qs,$ss,$hs,$qname,$sname,$qlength,$slength,$qframe,$sframe,
142
 
        $blasttype) = 
143
 
            $self->_rearrange([qw(SCORE
144
 
                                  BITS
145
 
                                  MATCH
146
 
                                  HSPLENGTH
147
 
                                  POSITIVE
148
 
                                  GAPS                            
149
 
                                  P
150
 
                                  EXP
151
 
                                  QUERYBEGIN
152
 
                                  QUERYEND
153
 
                                  SBJCTBEGIN
154
 
                                  SBJCTEND
155
 
                                  QUERYSEQ
156
 
                                  SBJCTSEQ
157
 
                                  HOMOLOGYSEQ
158
 
                                  QUERYNAME
159
 
                                  SBJCTNAME
160
 
                                  QUERYLENGTH
161
 
                                  SBJCTLENGTH
162
 
                                  QUERYFRAME
163
 
                                  SBJCTFRAME
164
 
                                  BLASTTYPE
165
 
                                  )],@args);
166
 
    
167
 
    $blasttype = 'UNKNOWN' unless $blasttype;
168
 
    $self->report_type($blasttype);
169
 
    # Determine strand meanings
170
 
    my ($queryfactor, $sbjctfactor) = (1,0); # default
171
 
    if ($blasttype eq 'BLASTP' || $blasttype eq 'TBLASTN' ) {
172
 
        $queryfactor = 0;
173
 
    }
174
 
    if ($blasttype eq 'TBLASTN' || $blasttype eq 'TBLASTX' || 
175
 
        $blasttype eq 'BLASTN' )  {
176
 
        $sbjctfactor = 1;
177
 
    }
178
 
    
179
 
    # Set BLAST type
180
 
    $self->{'BLAST_TYPE'} = $blasttype;
181
 
        
182
 
    # Store the aligned query as sequence feature
183
 
    my $strand;
184
 
    if ($qe > $qb) {            # normal query: start < end
185
 
                if ($queryfactor) { $strand = 1; } else { $strand = undef; }
186
 
                $self->query( Bio::SeqFeature::Similarity->new
187
 
                      (-start=>$qb, -end=>$qe, -strand=>$strand, 
188
 
                       -source=>"BLAST" ) ) }
189
 
    else {                      # reverse query (i dont know if this is possible, but feel free to correct)     
190
 
                if ($queryfactor) { $strand = -1; } else { $strand = undef; }
191
 
                $self->query( Bio::SeqFeature::Similarity->new
192
 
                      (-start=>$qe, -end=>$qb, -strand=>$strand,
193
 
                       -source=>"BLAST" ) ) }
194
 
 
195
 
    # store the aligned hit as sequence feature
196
 
    if ($se > $sb) {            # normal hit
197
 
        if ($sbjctfactor) { $strand = 1; } else { $strand = undef; }
198
 
        $self->hit( Bio::SeqFeature::Similarity->new
199
 
                        (-start=>$sb, -end=>$se, -strand=>$strand,
200
 
                         -source=>"BLAST" ) ) }
201
 
    else { # reverse hit: start bigger than end
202
 
        if ($sbjctfactor) { $strand = -1; } else { $strand = undef; }
203
 
        $self->hit( Bio::SeqFeature::Similarity->new
204
 
                        (-start=>$se, -end=>$sb, -strand=>$strand,
205
 
                         -source=>"BLAST" ) ) }
206
 
    
207
 
    # name the sequences
208
 
    $self->query->seq_id($qname); # query name
209
 
    $self->hit->seq_id($sname);   # hit name
210
 
 
211
 
    # set lengths
212
 
    $self->query->seqlength($qlength); # query length
213
 
    $self->hit->seqlength($slength);   # hit length
214
 
 
215
 
    # set object vars
216
 
    $self->score($score);
217
 
    $self->bits($bits);
218
 
 
219
 
    $self->significance($p);
220
 
    $self->{'EXP'} = $exp;
221
 
    
222
 
    $self->query->frac_identical($match);
223
 
    $self->hit->frac_identical($match);
224
 
    $self->{'HSPLENGTH'} = $hsplength;
225
 
    $self->{'PERCENT'} = int((1000 * $match)/$hsplength)/10;
226
 
    $self->{'POSITIVE'} = $positive;
227
 
    $self->{'GAPS'} = $gaps;
228
 
    $self->{'QS'} = $qs;
229
 
    $self->{'SS'} = $ss;
230
 
    $self->{'HS'} = $hs;
231
 
    
232
 
    $self->frame($qframe, $sframe);
233
 
    return $self;               # success - we hope!
234
 
}
235
 
 
236
 
# to disable overloading comment this out:
237
 
sub _overload {
238
 
        my $self = shift;
239
 
        return $self->start."..".$self->end." ".$self->bits;
240
 
}
241
 
 
242
 
=head2 report_type
243
 
 
244
 
 Title    : report_type
245
 
 Usage    : $type = $sbjct->report_type()
246
 
 Function : Returns the type of report from which this hit was obtained.
247
 
            This usually pertains only to BLAST and friends reports, for which
248
 
            the report type denotes what type of sequence was aligned against
249
 
            what (BLASTN: dna-dna, BLASTP prt-prt, BLASTX translated dna-prt, 
250
 
            TBLASTN prt-translated dna, TBLASTX translated dna-translated dna).
251
 
 Example  : 
252
 
 Returns  : A string (BLASTN, BLASTP, BLASTX, TBLASTN, TBLASTX, UNKNOWN)
253
 
 Args     : a string on set (you should know what you are doing)
254
 
 
255
 
=cut
256
 
 
257
 
sub report_type {
258
 
    my ($self, $rpt) = @_;
259
 
    if($rpt) {
260
 
        $self->{'_report_type'} = $rpt;
261
 
    }
262
 
    return $self->{'_report_type'};
263
 
}
264
 
 
265
 
=head2 EXP
266
 
 
267
 
 Title   : EXP
268
 
 Usage   : my $exp = $hsp->EXP;
269
 
 Function: returns the EXP value for the HSP
270
 
 Returns : string value
271
 
 Args    : none
272
 
 Note    : Patch provided by Sami Ashour for BTK parsing
273
 
 
274
 
 
275
 
=cut
276
 
 
277
 
sub EXP{
278
 
    return $_[0]->{'EXP'};
279
 
}
280
 
 
281
 
 
282
 
=head2 P
283
 
 
284
 
 Title    : P
285
 
 Usage    : $hsp->P();
286
 
 Function : returns the P (significance) value for a HSP 
287
 
 Returns  : (double) significance value
288
 
 Args     :
289
 
 
290
 
=cut
291
 
 
292
 
sub P {
293
 
        my ($self, @args) = @_;
294
 
        my $float = $self->significance(@args);
295
 
        my $match = '([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?'; # Perl Cookbook 2.1
296
 
        if ($float =~ /^$match$/) {
297
 
            # Is a C float
298
 
            return $float;
299
 
        } elsif ("1$float" =~ /^$match$/) {
300
 
            # Almost C float, Jitterbug 974
301
 
            return "1$float";
302
 
        } else {
303
 
                $self->warn("[HSP::P()] '$float' is not a known number format. Returning zero (0) instead.");
304
 
                return 0;
305
 
        }
306
 
}
307
 
 
308
 
=head2 percent
309
 
 
310
 
 Title    : percent
311
 
 Usage    : $hsp->percent();
312
 
 Function : returns the percent matching 
313
 
 Returns  : (double) percent matching
314
 
 Args     : none
315
 
 
316
 
=cut
317
 
 
318
 
sub percent         {shift->{'PERCENT'}}
319
 
 
320
 
 
321
 
=head2 match
322
 
 
323
 
 Title    : match
324
 
 Usage    : $hsp->match();
325
 
 Function : returns the match
326
 
 Example  : 
327
 
 Returns  : (double) frac_identical 
328
 
 Args     :
329
 
 
330
 
=cut
331
 
 
332
 
sub match           {shift->query->frac_identical(@_)}
333
 
 
334
 
=head2 hsplength
335
 
 
336
 
 Title    : hsplength
337
 
 Usage    : $hsp->hsplength();
338
 
 Function : returns the HSP length (including gaps)
339
 
 Returns  : (integer) HSP length
340
 
 Args     : none
341
 
 
342
 
=cut
343
 
 
344
 
sub hsplength              {shift->{'HSPLENGTH'}}
345
 
 
346
 
=head2 positive
347
 
 
348
 
 Title    : positive
349
 
 Usage    : $hsp->positive();
350
 
 Function : returns the number of positive matches (symbols in the alignment
351
 
            with a positive score)
352
 
 Returns  : (int) number of positive matches in the alignment
353
 
 Args     : none
354
 
 
355
 
=cut
356
 
 
357
 
sub positive        {shift->{'POSITIVE'}}
358
 
 
359
 
=head2 gaps
360
 
 
361
 
 Title    : gaps
362
 
 Usage    : $hsp->gaps();
363
 
 Function : returns the number of gaps or 0 if none
364
 
 Returns  : (int) number of gaps or 0 if none
365
 
 Args     : none
366
 
 
367
 
=cut
368
 
 
369
 
sub gaps        {shift->{'GAPS'}}
370
 
 
371
 
=head2 querySeq
372
 
 
373
 
 Title    : querySeq
374
 
 Usage    : $hsp->querySeq();
375
 
 Function : returns the query sequence
376
 
 Returns  : (string) the Query Sequence 
377
 
 Args     : none
378
 
 
379
 
=cut
380
 
 
381
 
sub querySeq        {shift->{'QS'}}
382
 
 
383
 
=head2 sbjctSeq
384
 
 
385
 
 Title    : sbjctSeq
386
 
 Usage    : $hsp->sbjctSeq();
387
 
 Function : returns the Sbjct sequence 
388
 
 Returns  : (string) the Sbjct Sequence 
389
 
 Args     : none
390
 
 
391
 
=cut
392
 
 
393
 
sub sbjctSeq        {shift->{'SS'}}
394
 
 
395
 
=head2 homologySeq
396
 
 
397
 
 Title    : homologySeq
398
 
 Usage    : $hsp->homologySeq();
399
 
 Function : returns the homologous sequence 
400
 
 Returns  : (string) homologous sequence 
401
 
 Args     : none
402
 
 
403
 
=cut
404
 
 
405
 
sub homologySeq     {shift->{'HS'}}
406
 
 
407
 
=head2 qs
408
 
 
409
 
 Title    : qs
410
 
 Usage    : $hsp->qs();
411
 
 Function : returns the Query Sequence (same as querySeq)
412
 
 Returns  : (string) query Sequence 
413
 
 Args     : none
414
 
 
415
 
=cut
416
 
 
417
 
sub qs              {shift->{'QS'}}
418
 
 
419
 
=head2 ss
420
 
 
421
 
 Title    : ss
422
 
 Usage    : $hsp->ss();
423
 
 Function : returns the subject sequence ( same as sbjctSeq) 
424
 
 Returns  : (string) Sbjct Sequence
425
 
 Args     : none
426
 
 
427
 
=cut
428
 
 
429
 
sub ss              {shift->{'SS'}}
430
 
 
431
 
=head2 hs
432
 
 
433
 
 Title    : hs
434
 
 Usage    : $hsp->hs();
435
 
 Function : returns the Homologous Sequence (same as homologySeq ) 
436
 
 Returns  : (string) Homologous Sequence
437
 
 Args     : none
438
 
 
439
 
=cut
440
 
 
441
 
sub hs              {shift->{'HS'}}
442
 
 
443
 
sub frame {
444
 
    my ($self, $qframe, $sframe) = @_;
445
 
    if( defined $qframe ) {
446
 
        if( $qframe == 0 ) {
447
 
            $qframe = undef;
448
 
        } elsif( $qframe !~ /^([+-])?([1-3])/ ) {           
449
 
            $self->warn("Specifying an invalid query frame ($qframe)");
450
 
            $qframe = undef;
451
 
        } else { 
452
 
            if( ($1 eq '-' && $self->query->strand >= 0) || 
453
 
                ($1 eq '+' && $self->query->strand <= 0) ) {
454
 
                $self->warn("Query frame ($qframe) did not match strand of query (". $self->query->strand() . ")");
455
 
            }
456
 
            # Set frame to GFF [0-2]
457
 
            $qframe = $2 - 1;
458
 
        }
459
 
        $self->{'QFRAME'} = $qframe;
460
 
    }
461
 
    if( defined $sframe ) {
462
 
          if( $sframe == 0 ) {
463
 
            $sframe = undef;
464
 
          } elsif( $sframe !~ /^([+-])?([1-3])/ ) {         
465
 
            $self->warn("Specifying an invalid hit frame ($sframe)");
466
 
            $sframe = undef;
467
 
          } else { 
468
 
              if( ($1 eq '-' && $self->hit->strand >= 0) || 
469
 
                  ($1 eq '+' && $self->hit->strand <= 0) ) 
470
 
              {
471
 
                  $self->warn("Hit frame ($sframe) did not match strand of hit (". $self->hit->strand() . ")");
472
 
              }
473
 
              
474
 
              # Set frame to GFF [0-2]
475
 
              $sframe = $2 - 1;
476
 
          }
477
 
          $self->{'SFRAME'} = $sframe;
478
 
      }
479
 
 
480
 
    (defined $qframe && $self->SUPER::frame($qframe) && 
481
 
     ($self->{'FRAME'} = $qframe)) || 
482
 
    (defined $sframe && $self->SUPER::frame($sframe) && 
483
 
     ($self->{'FRAME'} = $sframe));
484
 
 
485
 
    if (wantarray() && 
486
 
        $self->{'BLAST_TYPE'} eq 'TBLASTX') 
487
 
    { 
488
 
        return ($self->{'QFRAME'}, $self->{'SFRAME'}); 
489
 
    } elsif (wantarray())  { 
490
 
        (defined $self->{'QFRAME'} && 
491
 
         return ($self->{'QFRAME'}, undef)) || 
492
 
             (defined $self->{'SFRAME'} && 
493
 
              return (undef, $self->{'SFRAME'})); 
494
 
    } else { 
495
 
        (defined $self->{'QFRAME'} && 
496
 
         return $self->{'QFRAME'}) || 
497
 
        (defined $self->{'SFRAME'} && 
498
 
         return $self->{'SFRAME'}); 
499
 
    }
500
 
}
501
 
 
502
 
1;