← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Bio/Tools/Sigcleave.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2007-09-21 22:52:22 UTC
  • mfrom: (1.2.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20070921225222-tt20m2yy6ycuy2d8
Tags: 1.5.2.102-1
* Developer release.
* Upgraded source package to debhelper 5 and standards-version 3.7.2.
* Added libmodule-build-perl and libtest-harness-perl to
  build-depends-indep.
* Disabled automatic CRAN download.
* Using quilt instead of .diff.gz to manage modifications.
* Updated Recommends list for the binary package.
* Moved the "production-quality" scripts to /usr/bin/.
* New maintainer: Debian-Med packaging team mailing list.
* New uploaders: Charles Plessy and Steffen Moeller.
* Updated Depends, Recommends and Suggests.
* Imported in Debian-Med's SVN repository on Alioth.
* Executing the regression tests during package building.
* Moved the Homepage: field out from the package's description.
* Updated watch file.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Tools/Sigcleave.pm
2
2
# PACKAGE : Bio::Tools::Sigcleave
3
3
# AUTHOR  : Chris Dagdigian, dag@sonsorol.org
4
4
# CREATED : Jan 28 1999
5
 
# REVISION: $Id: Sigcleave.pm,v 1.18 2003/06/04 08:36:43 heikki Exp $
 
5
# REVISION: $Id: Sigcleave.pm,v 1.22.4.1 2006/10/02 23:10:32 sendu Exp $
6
6
#
7
7
# Copyright (c) 1997-9 bioperl, Chris Dagdigian and others. All Rights Reserved.
8
8
#           This module is free software; you can redistribute it and/or 
49
49
                                                -threshold=>'3.5',
50
50
                                                );
51
51
 
52
 
 
53
52
  # now you can detect procaryotic signal sequences as well as eucaryotic
54
53
  $sig->matrix('eucaryotic'); # or 'procaryotic'
55
54
 
130
129
=head2 Reporting Bugs
131
130
 
132
131
Report bugs to the Bioperl bug tracking system to help us keep track
133
 
the bugs and their resolution. Bug reports can be submitted via email
134
 
or the web:
 
132
the bugs and their resolution. Bug reports can be submitted via the
 
133
web:
135
134
 
136
 
    bioperl-bugs@bio.perl.org
137
 
    http://bugzilla.bioperl.org/
 
135
  http://bugzilla.open-bio.org/
138
136
 
139
137
=head1 AUTHOR
140
138
 
141
 
Chris Dagdigian, dag@sonsorol.org  & others
 
139
Chris Dagdigian, dag-at-sonsorol.org  & others
142
140
 
143
141
=head1 CONTRIBUTORS
144
142
 
145
 
Heikki Lehvaslaiho, heikki@ebi.ac.uk
 
143
Heikki Lehvaslaiho, heikki-at-bioperl-dot-org
146
144
 
147
145
=head1 VERSION
148
146
 
149
 
Bio::Tools::Sigcleave.pm, $Id: Sigcleave.pm,v 1.18 2003/06/04 08:36:43 heikki Exp $
 
147
Bio::Tools::Sigcleave, $Id: Sigcleave.pm,v 1.22.4.1 2006/10/02 23:10:32 sendu Exp $
150
148
 
151
149
=head1 COPYRIGHT
152
150
 
170
168
use and are not meant to be called by the user; they are 
171
169
preceded by an underscore ("_").
172
170
 
173
 
 
174
171
=cut
175
172
 
176
173
#
181
178
##
182
179
#
183
180
 
184
 
 
185
181
package Bio::Tools::Sigcleave;
186
182
 
187
 
use Bio::Root::Root;
188
183
use Bio::PrimarySeq;
189
184
 
190
 
@ISA = qw(Bio::Root::Root);
 
185
use base qw(Bio::Root::Root);
191
186
use strict;
192
187
use vars qw ($ID %WeightTable_euc  %WeightTable_pro );
193
188
$ID  = 'Bio::Tools::Sigcleave';
257
252
 
258
253
 
259
254
foreach my $i (keys %WeightTable_euc) {
260
 
    my $expected = $WeightTable_euc{$i}[15];
261
 
    if ($expected > 0) {
262
 
        for (my $j=0; $j<16; $j++) {
263
 
            if ($WeightTable_euc{$i}[$j] == 0) {
264
 
                $WeightTable_euc{$i}[$j] = 1; 
265
 
                if ($j == 10 || $j == 12) {
266
 
                    $WeightTable_euc{$i}[$j] = 1.e-10;
 
255
        my $expected = $WeightTable_euc{$i}[15];
 
256
        if ($expected > 0) {
 
257
                for (my $j=0; $j<16; $j++) {
 
258
                        if ($WeightTable_euc{$i}[$j] == 0) {
 
259
                                $WeightTable_euc{$i}[$j] = 1; 
 
260
                                if ($j == 10 || $j == 12) {
 
261
                                        $WeightTable_euc{$i}[$j] = 1.e-10;
 
262
                                }
 
263
                        }
 
264
                        $WeightTable_euc{$i}[$j] = log($WeightTable_euc{$i}[$j]/$expected);
267
265
                }
268
 
            }
269
 
            $WeightTable_euc{$i}[$j] = log($WeightTable_euc{$i}[$j]/$expected);
270
266
        }
271
 
    }
272
267
}
273
268
 
274
269
 
275
270
foreach my $i (keys %WeightTable_pro) {
276
 
    my $expected = $WeightTable_pro{$i}[15];
277
 
    if ($expected > 0) {
278
 
        for (my $j=0; $j<16; $j++) {
279
 
            if ($WeightTable_pro{$i}[$j] == 0) {
280
 
                $WeightTable_pro{$i}[$j] = 1; 
281
 
                if ($j == 10 || $j == 12) {
282
 
                    $WeightTable_pro{$i}[$j] = 1.e-10;
 
271
        my $expected = $WeightTable_pro{$i}[15];
 
272
        if ($expected > 0) {
 
273
                for (my $j=0; $j<16; $j++) {
 
274
                        if ($WeightTable_pro{$i}[$j] == 0) {
 
275
                                $WeightTable_pro{$i}[$j] = 1; 
 
276
                                if ($j == 10 || $j == 12) {
 
277
                                        $WeightTable_pro{$i}[$j] = 1.e-10;
 
278
                                }
 
279
                        }
 
280
                        $WeightTable_pro{$i}[$j] = log($WeightTable_pro{$i}[$j]/$expected);
283
281
                }
284
 
            }
285
 
            $WeightTable_pro{$i}[$j] = log($WeightTable_pro{$i}[$j]/$expected);
286
282
        }
287
 
    }
288
283
}
289
284
 
290
 
 
291
 
 
292
285
#####################################################################################
293
286
##                                 CONSTRUCTOR                                     ##
294
287
#####################################################################################
327
320
#----------------
328
321
sub threshold {
329
322
#----------------
330
 
    my ($self, $value) = @_;
331
 
    if( defined $value) {
332
 
        $self->throw("I need a number, not [$value]")
333
 
            if  $value !~ /^[+-]?[\d\.]+$/;
334
 
        $self->{'_threshold'} = $value;
335
 
    }
336
 
    return $self->{'_threshold'} || 3.5 ;
 
323
        my ($self, $value) = @_;
 
324
        if( defined $value) {
 
325
                $self->throw("I need a number, not [$value]")
 
326
                  if  $value !~ /^[+-]?[\d\.]+$/;
 
327
                $self->{'_threshold'} = $value;
 
328
        }
 
329
        return $self->{'_threshold'} || 3.5 ;
337
330
}
338
331
 
339
332
=head1 matrix
352
345
#----------------
353
346
sub matrix {
354
347
#----------------
355
 
    my ($self, $value) = @_;
356
 
    if( defined $value) {
357
 
        $self->throw("I need 'eucaryotic' or 'procaryotic', not [$value]")
358
 
            unless  $value eq 'eucaryotic' or $value eq 'procaryotic';
359
 
        $self->{'_matrix'} = $value;
360
 
    }
361
 
    return $self->{'_matrix'} || 'eucaryotic' ;
362
 
 
 
348
        my ($self, $value) = @_;
 
349
        if( defined $value) {
 
350
                $self->throw("I need 'eucaryotic' or 'procaryotic', not [$value]")
 
351
                  unless  $value eq 'eucaryotic' or $value eq 'procaryotic';
 
352
                $self->{'_matrix'} = $value;
 
353
        }
 
354
        return $self->{'_matrix'} || 'eucaryotic' ;
363
355
}
364
356
 
365
357
=head1 seq
366
358
 
367
359
 Title     : seq
368
 
 Usage     : $value = $self->seq('procaryotic')
369
 
 Purpose   : Read/write method sigcleave seq.
370
 
 Returns   : float.
371
 
 Argument  : new value: 'eucaryotic' or 'procaryotic'
372
 
 Throws    : on non-number argument
373
 
 Comments  : defaults to 3.5
 
360
 Usage     : $value = $self->seq($seq_object)
 
361
 Purpose   : set the Seq object to be used
 
362
 Returns   : Seq object
 
363
 Argument  : protein sequence or Seq object
374
364
 See Also   : n/a
375
365
 
376
366
=cut
378
368
#----------------
379
369
sub seq {
380
370
#----------------
381
 
    my ($self, $value) = @_;
382
 
    if( defined $value) {
383
 
        if ($value->isa('Bio::PrimarySeqI')) {
384
 
            $self->{'_seq'} = $value;
385
 
        } else {
386
 
            $self->{'_seq'} = Bio::PrimarySeq->new(-seq=>$value, 
387
 
                                                   -alphabet=>'protein');
 
371
        my ($self, $value) = @_;
 
372
        if( defined $value) {
 
373
                if ($value->isa('Bio::PrimarySeqI')) {
 
374
                        $self->{'_seq'} = $value;
 
375
                } else {
 
376
                        $self->{'_seq'} = Bio::PrimarySeq->new(-seq => $value, 
 
377
                                                                                                                                -alphabet => 'protein');
 
378
                }
388
379
        }
389
 
    }
390
 
    return $self->{'_seq'};
 
380
        return $self->{'_seq'};
391
381
}
392
382
 
393
 
 
394
 
 
395
383
=head1 _Analyze
396
384
 
397
385
 Title     : _Analyze
446
434
    $pep =~ tr/a-z/A-Z/;
447
435
 
448
436
    for ($seqPos = $seqBegin; $seqPos < $seqEnd; $seqPos++) {
449
 
        $istart = (0 > $seqPos + $pVal)? 0 : $seqPos + $pVal;
450
 
        $iend = ($seqPos + $nVal - 1 < $seqEnd)? $seqPos + $nVal - 1 : $seqEnd;
451
 
        $icol= $iend - $istart + 1;
452
 
        $weight = 0.00;
453
 
        for ($k=0; $k<$icol; $k++) {
454
 
            $c = substr($pep, $istart + $k, 1);
 
437
                 $istart = (0 > $seqPos + $pVal)? 0 : $seqPos + $pVal;
 
438
                 $iend = ($seqPos + $nVal - 1 < $seqEnd)? $seqPos + $nVal - 1 : $seqEnd;
 
439
                 $icol= $iend - $istart + 1;
 
440
                 $weight = 0.00;
 
441
                 for ($k=0; $k<$icol; $k++) {
 
442
                         $c = substr($pep, $istart + $k, 1);
455
443
 
456
 
            ## CD: The if(defined) stuff was put in here because Sigcleave.pm
457
 
            ## CD: kept getting warnings about undefined vals during 'make test' ...
458
 
            if ($matrix eq 'eucaryotic') {
459
 
                $weight += $WeightTable_euc{$c}[$k] if defined $WeightTable_euc{$c}[$k];
460
 
            } else {
461
 
                $weight += $WeightTable_pro{$c}[$k] if defined $WeightTable_pro{$c}[$k];
462
 
            }
463
 
        }
464
 
        $signals{$seqPos+1} = sprintf ("%.1f", $weight) if $weight >= $minWeight;
 
444
                         ## CD: The if(defined) stuff was put in here because Sigcleave.pm
 
445
                         ## CD: kept getting warnings about undefined vals during 'make test' ...
 
446
                         if ($matrix eq 'eucaryotic') {
 
447
                                 $weight += $WeightTable_euc{$c}[$k] if defined $WeightTable_euc{$c}[$k];
 
448
                         } else {
 
449
                                 $weight += $WeightTable_pro{$c}[$k] if defined $WeightTable_pro{$c}[$k];
 
450
                         }
 
451
                 }
 
452
                 $signals{$seqPos+1} = sprintf ("%.1f", $weight)        if $weight >= $minWeight;
465
453
    }
466
 
 
467
454
    $self->{"_signal_scores"} = { %signals };
468
455
}
469
456
 
489
476
#----------------
490
477
sub signals {
491
478
#----------------
492
 
    my $self = shift;
493
 
    my %results;
494
 
    my $position;
495
 
 
496
 
    # do the calculations
497
 
    $self->_Analyze;
498
 
 
499
 
    foreach $position ( sort keys %{ $self->{'_signal_scores'} } ) {
500
 
        $results{$position} = $self->{'_signal_scores'}{$position};
501
 
    }
502
 
    return %results;
 
479
        my $self = shift;
 
480
        my %results;
 
481
        my $position;
 
482
 
 
483
        # do the calculations
 
484
        $self->_Analyze;
 
485
 
 
486
        foreach $position ( sort keys %{ $self->{'_signal_scores'} } ) {
 
487
                $results{$position} = $self->{'_signal_scores'}{$position};
 
488
        }
 
489
        return %results;
503
490
}
504
491
 
505
492
 
523
510
#----------------
524
511
sub result_count {
525
512
#----------------
526
 
    my $self = shift;
527
 
    $self->_Analyze;
528
 
    return keys %{ $self->{'_signal_scores'} };
 
513
        my $self = shift;
 
514
        $self->_Analyze;
 
515
        return keys %{ $self->{'_signal_scores'} };
529
516
}
530
517
 
531
518
 
565
552
    $output = "SIGCLEAVE of $name from: 1 to $seqlen\n\n";
566
553
 
567
554
    if ($hitcount > 0) {
568
 
        $output .= "Report scores over $thresh\n";
569
 
        foreach $pos ((sort { $results{$b} cmp $results{$a} } keys %results)) {
570
 
            my $start = $pos - 15;
571
 
            $start = 1 if $start < 1;
572
 
            my $sig = substr($pep,$start -1,$pos-$start );
573
 
 
574
 
            $output .= sprintf ("Maximum score %1.1f at residue %3d\n",$results{$pos},$pos);
575
 
            $output .= "\n";
576
 
            $output .= " Sequence:  ";
577
 
            $output .= $sig;
578
 
            $output .= "-" x (15- length($sig));
579
 
            $output .= "-";
580
 
            $output .= substr($pep,$pos-1,50);
581
 
            $output .= "\n";
582
 
            $output .= " " x 12;
583
 
            $output .= "| \(signal\)      | \(mature peptide\)\n";
584
 
            $output .= sprintf("          %3d             %3d\n\n",$start,$pos);
585
 
 
586
 
            if (($hitcount > 1) && ($cnt == 1)) {
587
 
                $output .= " Other entries above $thresh\n\n";
588
 
            }
589
 
            $cnt++;
590
 
        }
 
555
                 $output .= "Report scores over $thresh\n";
 
556
                 foreach $pos ((sort { $results{$b} cmp $results{$a} } keys %results)) {
 
557
                         my $start = $pos - 15;
 
558
                         $start = 1 if $start < 1;
 
559
                         my $sig = substr($pep,$start -1,$pos-$start );
 
560
 
 
561
                         $output .= sprintf ("Maximum score %1.1f at residue %3d\n",$results{$pos},$pos);
 
562
                         $output .= "\n";
 
563
                         $output .= " Sequence:  ";
 
564
                         $output .= $sig;
 
565
                         $output .= "-" x (15- length($sig));
 
566
                         $output .= "-";
 
567
                         $output .= substr($pep,$pos-1,50);
 
568
                         $output .= "\n";
 
569
                         $output .= " " x 12;
 
570
                         $output .= "| \(signal\)      | \(mature peptide\)\n";
 
571
                         $output .= sprintf("          %3d             %3d\n\n",$start,$pos);
 
572
 
 
573
                         if (($hitcount > 1) && ($cnt == 1)) {
 
574
                                 $output .= " Other entries above $thresh\n\n";
 
575
                         }
 
576
                         $cnt++;
 
577
                 }
591
578
    }
592
579
    $output;
593
580
}