← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Bio/Seq/SequenceTrace.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/Seq/SequenceTrace.pm
1
 
# $Id: SequenceTrace.pm,v 1.6 2003/06/06 12:25:03 heikki Exp $
 
1
# $Id: SequenceTrace.pm,v 1.14.4.1 2006/10/02 23:10:27 sendu Exp $
2
2
#
3
3
# BioPerl module for Bio::Seq::SequenceTrace
4
4
#
30
30
Bioperl modules. Send your comments and suggestions preferably to one
31
31
of the Bioperl mailing lists.  Your participation is much appreciated.
32
32
 
33
 
  bioperl-l@bioperl.org             - General discussion
34
 
  http://bio.perl.org/MailList.html - About the mailing lists
 
33
  bioperl-l@bioperl.org                  - General discussion
 
34
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
35
35
 
36
36
=head2 Reporting Bugs
37
37
 
38
38
Report bugs to the Bioperl bug tracking system to help us keep track
39
 
the bugs and their resolution.  Bug reports can be submitted via email
40
 
or the web:
 
39
the bugs and their resolution.  Bug reports can be submitted via the
 
40
web:
41
41
 
42
 
  bioperl-bugs@bio.perl.org
43
 
  http://bugzilla.bioperl.org/
 
42
  http://bugzilla.open-bio.org/
44
43
 
45
44
=head1 AUTHOR - Chad Matsalla
46
45
 
55
54
 
56
55
package Bio::Seq::SequenceTrace;
57
56
 
58
 
use vars qw(@ISA);
59
57
 
60
58
use strict;
61
 
use Bio::Root::Root;
62
59
use Bio::Seq::QualI;
63
60
use Bio::PrimarySeqI;
64
61
use Bio::PrimarySeq;
65
62
use Bio::Seq::PrimaryQual;
66
 
use Bio::Seq::TraceI;
67
 
use Bio::Seq::SeqWithQuality;
68
63
use Dumpvalue();
69
64
 
70
65
my $dumper = new Dumpvalue();
71
66
 
72
 
@ISA = qw(Bio::Root::Root Bio::Seq::SeqWithQuality Bio::Seq::TraceI);
 
67
use base qw(Bio::Root::Root Bio::Seq::Quality Bio::Seq::TraceI);
73
68
 
74
69
=head2 new()
75
70
 
112
107
               ACCURACY_G
113
108
               ACCURACY_C )], @args);
114
109
          # first, deal with the sequence and quality information
115
 
     if ($swq && ref($swq) eq "Bio::Seq::SeqWithQuality") {
 
110
     if ($swq && ref($swq) eq "Bio::Seq::Quality") {
116
111
          $self->{swq} = $swq;
117
112
     }
118
113
     else {
119
114
          $self->throw("A Bio::Seq::SequenceTrace object must be created with a
120
 
               Bio::Seq::SeqWithQuality object. You provided this type of object: "
 
115
               Bio::Seq::Quality object. You provided this type of object: "
121
116
               .ref($swq));
122
117
     }
123
118
     if (!$acc_a) {
174
169
          $self->throw('You must provide a valid base channel (atgc) to use trace()');
175
170
     }
176
171
     $base_channel =~ tr/A-Z/a-z/;
177
 
     if ($base_channel !~ /a|t|g|c/) {
 
172
     if ($base_channel !~ /[acgt]/) {
178
173
          $self->throw('You must provide a valid base channel (atgc) to use trace()');
179
174
     }
180
175
     if ($values) {
190
185
          return $self->{trace}->{$base_channel};
191
186
     }
192
187
     else {
193
 
          return undef;
 
188
          return;
194
189
     }
195
190
}
196
191
 
293
288
 
294
289
sub alphabet {
295
290
        my $self = shift;
296
 
        return $self->{swq}->{seq_ref}->alphabet();     
 
291
        return $self->{swq}->{seq_ref}->alphabet(@_);
297
292
}
298
293
 
299
294
=head2 display_id()
313
308
        field for example, or extending the sequence object) to solve
314
309
        this. Notice that $seq->id() maps to this function, mainly for
315
310
        legacy/convience issues.
316
 
        This method sets the display_id for the SeqWithQuality object.
 
311
        This method sets the display_id for the Quality object.
317
312
 Returns : A string
318
313
 Args    : If a scalar is provided, it is set as the new display_id for
319
 
        the SeqWithQuality object.
 
314
        the Quality object.
320
315
 Status  : Virtual
321
316
 
322
317
=cut
341
336
        for the implemetation, allowing multiple objects to have the same
342
337
        accession number in a particular implementation. For sequences
343
338
        with no accession number, this method should return "unknown".
344
 
        This method sets the accession_number for the SeqWithQuality
 
339
        This method sets the accession_number for the Quality
345
340
        object. 
346
341
 Returns : A string (the value of accession_number)
347
342
 Args    : If a scalar is provided, it is set as the new accession_number
348
 
        for the SeqWithQuality object.
 
343
        for the Quality object.
349
344
 Status  : Virtual
350
345
 
351
346
 
371
366
        way the implementaiton can control clients can expect one id to
372
367
        map to one object. For sequences with no accession number, this
373
368
        method should return a stringified memory location.
374
 
        This method sets the primary_id for the SeqWithQuality
 
369
        This method sets the primary_id for the Quality
375
370
        object.
376
371
 Returns : A string. (the value of primary_id)
377
372
 Args    : If a scalar is provided, it is set as the new primary_id for
378
 
        the SeqWithQuality object.
 
373
        the Quality object.
379
374
 
380
375
=cut
381
376
 
393
388
 Title   : desc()
394
389
 Usage   : $qual->desc($newval); _or_ 
395
390
           $description = $qual->desc();
396
 
 Function: Get/set description text for this SeqWithQuality object.
 
391
 Function: Get/set description text for this Quality object.
397
392
 Returns : A string. (the value of desc)
398
393
 Args    : If a scalar is provided, it is set as the new desc for the
399
 
        SeqWithQuality object.
 
394
           Quality object.
400
395
 
401
396
=cut
402
397
 
403
398
sub desc {
404
 
        # a mechanism to set the desc for the SeqWithQuality object.
 
399
        # a mechanism to set the desc for the Quality object.
405
400
        # probably will be used most often by set_common_features()
406
401
   my ($self,$value) = @_;
407
402
   if( defined $value) {
419
414
        for display_id().
420
415
 Returns : A string. (the value of id)
421
416
 Args    : If a scalar is provided, it is set as the new id for the
422
 
        SeqWithQuality object.
 
417
           Quality object.
423
418
 
424
419
=cut
425
420
 
438
433
 Usage   : $string    = $obj->seq(); _or_
439
434
        $obj->seq("atctatcatca");
440
435
 Function: Returns the sequence that is contained in the imbedded in the
441
 
        PrimarySeq object within the SeqWithQuality object
 
436
        PrimarySeq object within the Quality object
442
437
 Returns : A scalar (the seq() value for the imbedded PrimarySeq object.)
443
 
 Args    : If a scalar is provided, the SeqWithQuality object will
 
438
 Args    : If a scalar is provided, the Quality object will
444
439
        attempt to set that as the sequence for the imbedded PrimarySeq
445
440
        object. Otherwise, the value of seq() for the PrimarySeq object
446
441
        is returned.
456
451
sub seq {
457
452
        my ($self,$value) = @_;
458
453
        if( defined $value) {
459
 
                $self->{swq}->{seq_ref}->seq($value);
 
454
                $self->{swq}->seq($value);
460
455
        }
461
 
        return $self->{swq}->{seq_ref}->seq();
 
456
        return $self->{swq}->seq();
462
457
}
463
458
 
464
459
=head2 qual()
467
462
 Usage   : @quality_values  = @{$obj->qual()}; _or_
468
463
        $obj->qual("10 10 20 40 50");
469
464
 Function: Returns the quality as imbedded in the PrimaryQual object
470
 
        within the SeqWithQuality object.
 
465
        within the Quality object.
471
466
 Returns : A reference to an array containing the quality values in the 
472
467
        PrimaryQual object.
473
 
 Args    : If a scalar is provided, the SeqWithQuality object will
 
468
 Args    : If a scalar is provided, the Quality object will
474
469
        attempt to set that as the quality for the imbedded PrimaryQual
475
470
        object. Otherwise, the value of qual() for the PrimaryQual
476
471
        object is returned.
486
481
        my ($self,$value) = @_;
487
482
 
488
483
        if( defined $value) {
489
 
                $self->{swq}->{qual_ref}->qual($value);
 
484
                $self->{swq}->qual($value);
490
485
        }
491
 
        return $self->{swq}->{qual_ref}->qual();
 
486
        return $self->{swq}->qual();
492
487
}
493
488
 
494
489
 
498
493
 
499
494
 Title   : length()
500
495
 Usage   : $length = $seqWqual->length();
501
 
 Function: Get the length of the SeqWithQuality sequence/quality.
502
 
 Returns : Returns the length of the sequence and quality if they are
503
 
        both the same. Returns "DIFFERENT" if they differ.
 
496
 Function: Get the length of the Quality sequence/quality.
 
497
 Returns : Returns the length of the sequence and quality
504
498
 Args    : None.
505
499
 
506
500
=cut
518
512
 Usage   : $qualobj = $seqWqual->qual_obj(); _or_
519
513
        $qualobj = $seqWqual->qual_obj($ref_to_primaryqual_obj);
520
514
 Function: Get the PrimaryQual object that is imbedded in the
521
 
        SeqWithQuality object or if a reference to a PrimaryQual object
 
515
        Quality object or if a reference to a PrimaryQual object
522
516
        is provided, set this as the PrimaryQual object imbedded in the
523
 
        SeqWithQuality object.
524
 
 Returns : A reference to a Bio::Seq::SeqWithQuality object.
 
517
        Quality object.
 
518
 Returns : A reference to a Bio::Seq::Quality object.
525
519
 
526
520
=cut
527
521
 
528
522
sub qual_obj {
529
523
    my ($self,$value) = @_;
530
 
    return $self->{swq}->qual_obj($value);
 
524
#    return $self->{swq}->qual_obj($value);
 
525
    return $self->{swq};
531
526
}
532
527
 
533
528
 
534
529
=head2 seq_obj
535
530
 
536
531
 Title   : seq_obj()
537
 
 Usage   : $seqobj = $seqWqual->qual_obj(); _or_
 
532
 Usage   : $seqobj = $seqWqual->seq_obj(); _or_
538
533
        $seqobj = $seqWqual->seq_obj($ref_to_primary_seq_obj);
539
534
 Function: Get the PrimarySeq object that is imbedded in the
540
 
        SeqWithQuality object or if a reference to a PrimarySeq object is
 
535
        Quality object or if a reference to a PrimarySeq object is
541
536
        provided, set this as the PrimarySeq object imbedded in the
542
 
        SeqWithQuality object.
 
537
        Quality object.
543
538
 Returns : A reference to a Bio::PrimarySeq object.
544
539
 
545
540
=cut
546
541
 
547
542
sub seq_obj {
548
543
    my ($self,$value) = @_;
549
 
    return $self->{swq}->seq_obj($value);
 
544
#    return $self->{swq}->seq_obj($value);
 
545
    return $self->{swq};
550
546
}
551
547
 
552
548
=head2 _set_descriptors
554
550
 Title   : _set_descriptors()
555
551
 Usage   : $seqWqual->_qual_obj($qual,$seq,$id,$acc,$pid,$desc,$given_id,
556
552
        $alphabet);
557
 
 Function: Set the descriptors for the SeqWithQuality object. Try to
 
553
 Function: Set the descriptors for the Quality object. Try to
558
554
        match the descriptors in the PrimarySeq object and in the
559
555
        PrimaryQual object if descriptors were not provided with
560
556
        construction.
766
762
        $start2 = shift(@subs);
767
763
        $end2 =  pop(@subs);
768
764
     my $new_object =  new Bio::Seq::SequenceTrace(
769
 
               -swq =>   new Bio::Seq::SeqWithQuality(
 
765
               -swq =>   new Bio::Seq::Quality(
770
766
                             -seq => $self->subseq($start,$end),
771
767
                             -qual     =>   $self->subqual($start,$end),
772
768
                             -id    =>   $self->id()
800
796
sub _synthesize_traces {
801
797
     my ($self) = shift;
802
798
     $self->peak_indices(qw());
803
 
     my $version = 2;
 
799
#ml     my $version = 2;
804
800
          # the user should be warned if traces already exist
805
801
          #
806
802
          #
807
 
     ( my $sequence = $self->seq() ) =~ tr/a-z/A-Z/;
808
 
     my @quals = @{$self->qual()};
809
 
     my $info;
 
803
#ml     ( my $sequence = $self->seq() ) =~ tr/a-z/A-Z/;
 
804
#ml     my @quals = @{$self->qual()};
 
805
#ml     my $info;
810
806
         # build the ramp for the first base.
811
807
         # a ramp looks like this "1 4 13 29 51 71 80 71 51 29 13 4 1" times the quality score.
812
808
         # REMEMBER: A C G T
827
823
     $self->initialize_traces("0",$total_length+2);
828
824
         # now populate them
829
825
    my ($current_base,$place_base_at,$peak_quality,$ramp_counter,$current_ramp,$ramp_position);
830
 
    my $sequence_length = $self->length();
 
826
#ml    my $sequence_length = $self->length();
831
827
    my $half_ramp = int($ramp_data->{'ramp_width'}/2);
832
828
    for ($pos = 0; $pos<$self->length();$pos++) {
833
 
          $current_base = $self->seq_obj()->subseq($pos+1,$pos+1);
 
829
          $current_base = uc $self->seq_obj()->subseq($pos+1,$pos+1);
834
830
               # print("Synthesizing the ramp for $current_base\n");
835
831
          my $all_bases = "ATGC";
836
832
          $peak_quality = $self->qual_obj()->qualat($pos+1);
851
847
          $self->peak_index_at($pos+1,
852
848
              $place_base_at+1
853
849
          );
854
 
          my $other_bases = $self->_get_other_bases($current_base);
 
850
#ml          my $other_bases = $self->_get_other_bases($current_base);
855
851
          # foreach ( split('',$other_bases) ) {
856
852
          #          push @{$self->{'text'}->{"v3_base_accuracy"}->{$_}},0;
857
853
          #}