3 # BioPerl module for Bio::Seq::SequenceTrace
5 # Cared for by Chad Matsalla <bioinformatics@dieselwurks.com
7 # Copyright Chad Matsalla
9 # You may distribute this module under the same terms as perl itself
11 # POD documentation - main docs before the code
15 Bio::Seq::SequenceTrace - Bioperl object packaging a sequence with its trace
23 This object stores a sequence with its trace.
29 User feedback is an integral part of the evolution of this and other
30 Bioperl modules. Send your comments and suggestions preferably to one
31 of the Bioperl mailing lists. Your participation is much appreciated.
33 bioperl-l@bioperl.org - General discussion
34 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
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 the
42 http://bugzilla.open-bio.org/
44 =head1 AUTHOR - Chad Matsalla
46 Email bioinformatics@dieselwurks.com
49 The rest of the documentation details each of the object methods.
50 Internal methods are usually preceded with a _
55 package Bio
::Seq
::SequenceTrace
;
62 use Bio
::Seq
::PrimaryQual
;
65 my $dumper = new Dumpvalue
();
67 use base
qw(Bio::Root::Root Bio::Seq::Quality Bio::Seq::TraceI);
72 Usage : $st = Bio::Seq::SequenceTrace->new
73 ( -swq => Bio::Seq::SequenceWithQuality,
74 -trace_a => \@trace_values_for_a_channel,
75 -trace_t => \@trace_values_for_t_channel,
76 -trace_g => \@trace_values_for_g_channel,
77 -trace_c => \@trace_values_for_c_channel,
78 -accuracy_a => \@a_accuracies,
79 -accuracy_t => \@t_accuracies,
80 -accuracy_g => \@g_accuracies,
81 -accuracy_c => \@c_accuracies,
82 -peak_indices => '0 5 10 15 20 25 30 35'
84 Function: Returns a new Bio::Seq::SequenceTrace object from basic
86 Returns : a new Bio::Seq::SequenceTrace object
87 Arguments: I think that these are all describes in the usage above.
92 my ($class, @args) = @_;
93 my $self = $class->SUPER::new
(@args);
94 # default: turn OFF the warnings
95 $self->{supress_warnings
} = 1;
96 my($swq,$peak_indices,$trace_a,$trace_t,
97 $trace_g,$trace_c,$acc_a,$acc_t,$acc_g,$acc_c) =
98 $self->_rearrange([qw(
108 ACCURACY_C )], @args);
109 # first, deal with the sequence and quality information
110 if ($swq && ref($swq) eq "Bio::Seq::Quality") {
114 $self->throw("A Bio::Seq::SequenceTrace object must be created with a
115 Bio::Seq::Quality object. You provided this type of object: "
119 # this means that you probably did not provide traces and accuracies
120 # and that they need to be synthesized
121 $self->set_accuracies();
124 $self->accuracies('a',$acc_a);
125 $self->accuracies('t',$acc_t);
126 $self->accuracies('g',$acc_g);
127 $self->accuracies('c',$acc_c);
130 $self->_synthesize_traces();
133 $self->trace('a',$trace_a);
134 $self->trace('t',$trace_t);
135 $self->trace('g',$trace_g);
136 $self->trace('c',$trace_c);
137 $self->peak_indices($peak_indices);
139 $self->id($self->swq_obj()->id());
150 =head2 trace($base,\@new_values)
152 Title : trace($base,\@new_values)
153 Usage : @trace_Values = @{$obj->trace($base,\@new_values)};
154 Function: Returns the trace values as a reference to an array containing the
155 trace values. The individual elements of the trace array are not validated
156 and can be any numeric value.
157 Returns : A reference to an array.
159 Arguments: $base : which color channel would you like the trace values for?
160 - $base must be one of "A","T","G","C"
161 \@new_values : a reference to an array of values containing trace
167 my ($self,$base_channel,$values) = @_;
168 if (!$base_channel) {
169 $self->throw('You must provide a valid base channel (atgc) to use trace()');
171 $base_channel =~ tr/A-Z/a-z/;
172 if ($base_channel !~ /[acgt]/) {
173 $self->throw('You must provide a valid base channel (atgc) to use trace()');
176 if (ref($values) eq "ARRAY") {
177 $self->{trace
}->{$base_channel} = $values;
180 my @trace = split(' ',$values);
181 $self->{trace
}->{$base_channel} = \
@trace;
184 if ($self->{trace
}->{$base_channel}) {
185 return $self->{trace
}->{$base_channel};
193 =head2 peak_indices($new_indices)
195 Title : peak_indices($new_indices)
196 Usage : $indices = $obj->peak_indices($new_indices);
197 Function: Return the trace index points for this object.
199 Args : If used, the trace indices will be set to the provided value.
204 my ($self,$peak_indices)= @_;
206 if (ref($peak_indices) eq "ARRAY") {
207 $self->{peak_indices
} = $peak_indices;
210 my @indices = split(' ',$peak_indices);
211 $self->{peak_indices
} = \
@indices;
214 if (!$self->{peak_indices
}) {
216 $self->{peak_indices
} = \
@temp;
218 return $self->{peak_indices
};
222 =head2 _reset_peak_indices()
224 Title : _rest_peak_indices()
225 Usage : $obj->_reset_peak_indices();
226 Function: Reset the peak indices.
229 Notes : When you create a sub_trace_object, the peak indices
230 will still be pointing to the apporpriate location _in the
231 original trace_. In order to fix this, the initial value must
232 be subtracted from each value here. ie. The first peak index
237 sub _reset_peak_indices
{
239 my $length = $self->length();
240 my $subtractive = $self->peak_index_at(1);
242 $self->peak_index_at(1,"null");
243 for (my $counter=2; $counter<= $length; $counter++) {
244 my $original = $self->peak_index_at($counter);
245 $new = $original - $subtractive;
246 $self->peak_index_at($counter,$new);
255 =head2 peak_index_at($position)
257 Title : peak_index_at($position)
258 Usage : $peak_index = $obj->peak_index_at($postition);
259 Function: Return the trace iindex point at this position
261 Args : If used, the trace index at this position will be
262 set to the provided value.
267 my ($self,$position,$value)= @_;
269 if ($value eq "null") {
270 $self->peak_indices->[$position-1] = "0";
273 $self->peak_indices->[$position-1] = $value;
276 return $self->peak_indices()->[$position-1];
282 Usage : $molecule_type = $obj->alphabet();
283 Function: Get the molecule type from the PrimarySeq object.
284 Returns : What what PrimarySeq says the type of the sequence is.
291 return $self->{swq
}->{seq_ref
}->alphabet(@_);
297 Usage : $id_string = $obj->display_id();
298 Function: Returns the display id, aka the common name of the Quality
300 The semantics of this is that it is the most likely string to be
301 used as an identifier of the quality sequence, and likely to have
302 "human" readability. The id is equivalent to the ID field of the
303 GenBank/EMBL databanks and the id field of the Swissprot/sptrembl
304 database. In fasta format, the >(\S+) is presumed to be the id,
305 though some people overload the id to embed other information.
306 Bioperl does not use any embedded information in the ID field,
307 and people are encouraged to use other mechanisms (accession
308 field for example, or extending the sequence object) to solve
309 this. Notice that $seq->id() maps to this function, mainly for
310 legacy/convience issues.
311 This method sets the display_id for the Quality object.
313 Args : If a scalar is provided, it is set as the new display_id for
320 my ($self,$value) = @_;
321 if( defined $value) {
322 $self->{swq
}->display_id($value);
324 return $self->{swq
}->display_id();
328 =head2 accession_number()
330 Title : accession_number()
331 Usage : $unique_biological_key = $obj->accession_number();
332 Function: Returns the unique biological id for a sequence, commonly
333 called the accession_number. For sequences from established
334 databases, the implementors should try to use the correct
335 accession number. Notice that primary_id() provides the unique id
336 for the implemetation, allowing multiple objects to have the same
337 accession number in a particular implementation. For sequences
338 with no accession number, this method should return "unknown".
339 This method sets the accession_number for the Quality
341 Returns : A string (the value of accession_number)
342 Args : If a scalar is provided, it is set as the new accession_number
343 for the Quality object.
349 sub accession_number
{
350 my( $self, $acc ) = @_;
352 $self->{swq
}->accession_number($acc);
354 $acc = $self->{swq
}->accession_number();
355 $acc = 'unknown' unless defined $acc;
363 Usage : $unique_implementation_key = $obj->primary_id();
364 Function: Returns the unique id for this object in this implementation.
365 This allows implementations to manage their own object ids in a
366 way the implementaiton can control clients can expect one id to
367 map to one object. For sequences with no accession number, this
368 method should return a stringified memory location.
369 This method sets the primary_id for the Quality
371 Returns : A string. (the value of primary_id)
372 Args : If a scalar is provided, it is set as the new primary_id for
378 my ($self,$value) = @_;
380 $self->{swq
}->primary_id($value);
382 return $self->{swq
}->primary_id();
389 Usage : $qual->desc($newval); _or_
390 $description = $qual->desc();
391 Function: Get/set description text for this Quality object.
392 Returns : A string. (the value of desc)
393 Args : If a scalar is provided, it is set as the new desc for the
399 # a mechanism to set the desc for the Quality object.
400 # probably will be used most often by set_common_features()
401 my ($self,$value) = @_;
402 if( defined $value) {
403 $self->{swq
}->desc($value);
405 return $self->{swq
}->desc();
411 Usage : $id = $qual->id();
412 Function: Return the ID of the quality. This should normally be (and
413 actually is in the implementation provided here) just a synonym
415 Returns : A string. (the value of id)
416 Args : If a scalar is provided, it is set as the new id for the
422 my ($self,$value) = @_;
423 if (!$self) { $self->throw("no value for self in $value"); }
424 if( defined $value ) {
425 $self->{swq
}->display_id($value);
427 return $self->{swq
}->display_id();
433 Usage : $string = $obj->seq(); _or_
434 $obj->seq("atctatcatca");
435 Function: Returns the sequence that is contained in the imbedded in the
436 PrimarySeq object within the Quality object
437 Returns : A scalar (the seq() value for the imbedded PrimarySeq object.)
438 Args : If a scalar is provided, the Quality object will
439 attempt to set that as the sequence for the imbedded PrimarySeq
440 object. Otherwise, the value of seq() for the PrimarySeq object
442 Notes : This is probably not a good idea because you then should call
443 length() to make sure that the sequence and quality are of the
444 same length. Even then, how can you make sure that this sequence
445 belongs with that quality? I provided this to give you rope to
446 hang yourself with. Tie it to a strong device and use a good
452 my ($self,$value) = @_;
453 if( defined $value) {
454 $self->{swq
}->seq($value);
456 return $self->{swq
}->seq();
462 Usage : @quality_values = @{$obj->qual()}; _or_
463 $obj->qual("10 10 20 40 50");
464 Function: Returns the quality as imbedded in the PrimaryQual object
465 within the Quality object.
466 Returns : A reference to an array containing the quality values in the
468 Args : If a scalar is provided, the Quality object will
469 attempt to set that as the quality for the imbedded PrimaryQual
470 object. Otherwise, the value of qual() for the PrimaryQual
472 Notes : This is probably not a good idea because you then should call
473 length() to make sure that the sequence and quality are of the
474 same length. Even then, how can you make sure that this sequence
475 belongs with that quality? I provided this to give you a strong
476 board with which to flagellate yourself.
481 my ($self,$value) = @_;
483 if( defined $value) {
484 $self->{swq
}->qual($value);
486 return $self->{swq
}->qual();
495 Usage : $length = $seqWqual->length();
496 Function: Get the length of the Quality sequence/quality.
497 Returns : Returns the length of the sequence and quality
504 return $self->seq_obj()->length();
511 Title : qual_obj($different_obj)
512 Usage : $qualobj = $seqWqual->qual_obj(); _or_
513 $qualobj = $seqWqual->qual_obj($ref_to_primaryqual_obj);
514 Function: Get the PrimaryQual object that is imbedded in the
515 Quality object or if a reference to a PrimaryQual object
516 is provided, set this as the PrimaryQual object imbedded in the
518 Returns : A reference to a Bio::Seq::Quality object.
523 my ($self,$value) = @_;
524 # return $self->{swq}->qual_obj($value);
532 Usage : $seqobj = $seqWqual->seq_obj(); _or_
533 $seqobj = $seqWqual->seq_obj($ref_to_primary_seq_obj);
534 Function: Get the PrimarySeq object that is imbedded in the
535 Quality object or if a reference to a PrimarySeq object is
536 provided, set this as the PrimarySeq object imbedded in the
538 Returns : A reference to a Bio::PrimarySeq object.
543 my ($self,$value) = @_;
544 # return $self->{swq}->seq_obj($value);
548 =head2 _set_descriptors
550 Title : _set_descriptors()
551 Usage : $seqWqual->_qual_obj($qual,$seq,$id,$acc,$pid,$desc,$given_id,
553 Function: Set the descriptors for the Quality object. Try to
554 match the descriptors in the PrimarySeq object and in the
555 PrimaryQual object if descriptors were not provided with
558 Args : $qual,$seq,$id,$acc,$pid,$desc,$given_id,$alphabet as found
560 Notes : Really only intended to be called by the new() method. If
561 you want to invoke a similar function try
562 set_common_descriptors().
567 sub _set_descriptors
{
568 my ($self,$qual,$seq,$id,$acc,$pid,$desc,$given_id,$alphabet) = @_;
569 $self->{swq
}->_seq_descriptors($qual,$seq,$id,$acc,$pid,$desc,$given_id,$alphabet);
572 =head2 subseq($start,$end)
574 Title : subseq($start,$end)
575 Usage : $subsequence = $obj->subseq($start,$end);
576 Function: Returns the subseq from start to end, where the first base
577 is 1 and the number is inclusive, ie 1-2 are the first two
578 bases of the sequence.
580 Args : Two positions.
585 my ($self,@args) = @_;
586 # does a single value work?
587 return $self->{swq
}->subseq(@args);
590 =head2 baseat($position)
592 Title : baseat($position)
593 Usage : $base_at_position_6 = $obj->baseat("6");
594 Function: Returns a single base at the given position, where the first
595 base is 1 and the number is inclusive, ie 1-2 are the first two
596 bases of the sequence.
603 my ($self,$val) = @_;
604 return $self->{swq
}->subseq($val,$val);
607 =head2 subqual($start,$end)
609 Title : subqual($start,$end)
610 Usage : @qualities = @{$obj->subqual(10,20);
611 Function: returns the quality values from $start to $end, where the
612 first value is 1 and the number is inclusive, ie 1-2 are the
613 first two bases of the sequence. Start cannot be larger than
614 end but can be equal.
615 Returns : A reference to an array.
616 Args : a start position and an end position
621 my ($self,@args) = @_;
622 return $self->{swq
}->subqual(@args);
625 =head2 qualat($position)
627 Title : qualat($position)
628 Usage : $quality = $obj->qualat(10);
629 Function: Return the quality value at the given location, where the
630 first value is 1 and the number is inclusive, ie 1-2 are the
631 first two bases of the sequence. Start cannot be larger than
632 end but can be equal.
639 my ($self,$val) = @_;
640 return $self->{swq
}->qualat($val);
643 =head2 sub_peak_index($start,$end)
645 Title : sub_peak_index($start,$end)
646 Usage : @peak_indices = @{$obj->sub_peak_index(10,20);
647 Function: returns the trace index values from $start to $end, where the
648 first value is 1 and the number is inclusive, ie 1-2 are the
649 first two trace indices for this channel.
650 Returns : A reference to an array.
651 Args : a start position and an end position
656 my ($self,$start,$end) = @_;
658 $self->throw("in sub_peak_index, start [$start] has to be greater than end [$end]");
661 if( $start <= 0 || $end > $self->length ) {
662 $self->throw("You have to have start positive and length less than the total length of sequence [$start:$end] Total ".$self->length."");
665 # remove one from start, and then length is end-start
669 my @sub_peak_index_array = @
{$self->{peak_indices
}}[$start..$end];
671 # return substr $self->seq(), $start, ($end-$start);
672 return \
@sub_peak_index_array;
676 =head2 sub_trace($start,$end)
678 Title : sub_trace($base_channel,$start,$end)
679 Usage : @trace_values = @{$obj->sub_trace('a',10,20)};
680 Function: returns the trace values from $start to $end, where the
681 first value is 1 and the number is inclusive, ie 1-2 are the
682 first two bases of the sequence. Start cannot be larger than
683 end but can be e_peak_index.
684 Returns : A reference to an array.
685 Args : a start position and an end position
690 my ($self,$base_channel,$start,$end) = @_;
692 $self->throw("in sub_trace, start [$start] has to be greater than end [$end]");
695 if( $start <= 0 || $end > $self->trace_length() ) {
696 $self->throw("You have to have start positive and length less than the total length of traces [$start:$end] Total ".$self->trace_length."");
699 # remove one from start, and then length is end-start
703 my @sub_peak_index_array = @
{$self->trace($base_channel)}[$start..$end];
705 # return substr $self->seq(), $start, ($end-$start);
706 return \
@sub_peak_index_array;
710 =head2 trace_length()
712 Title : trace_length()
713 Usage : $trace_length = $obj->trace_length();
714 Function: Return the length of the trace if all four traces (atgc)
715 are the same. Otherwise, throw an error.
723 if ( !$self->trace('a') || !$self->trace('t') || !$self->trace('g') || !$self->trace('c') ) {
724 $self->warn("One or more of the trace channels are missing. Cannot give you a length.");
727 my $lengtha = scalar(@
{$self->trace('a')});
728 my $lengtht = scalar(@
{$self->trace('t')});
729 my $lengthg = scalar(@
{$self->trace('g')});
730 my $lengthc = scalar(@
{$self->trace('c')});
731 if (($lengtha == $lengtht) && ($lengtha == $lengthg) && ($lengtha == $lengthc) ) {
734 $self->warn("Not all of the trace indices are the same length".
735 " Here are their lengths: a: $lengtha t:$lengtht ".
736 " g: $lengthg c: $lengthc");
741 =head2 sub_trace_object($start,$end)
743 Title : sub_trace_object($start,$end)
744 Usage : $smaller_object = $object->sub_trace_object('1','100');
745 Function: Get a subset of the sequence, its quality, and its trace.
746 Returns : A reference to a Bio::Seq::SequenceTrace object
747 Args : a start position and an end position
749 - the start and end position refer to the positions of _bases_.
750 - for example, to get a sub SequenceTrace for bases 5-10,
752 - you will get the bases, qualities, and the trace values
753 - you can then use this object to synthesize a new scf
758 sub sub_trace_object
{
759 my ($self,$start,$end) = @_;
761 my @subs = @
{$self->sub_peak_index($start,$end)};
762 $start2 = shift(@subs);
764 my $new_object = Bio
::Seq
::SequenceTrace
->new(
765 -swq
=> Bio
::Seq
::Quality
->new(
766 -seq
=> $self->subseq($start,$end),
767 -qual
=> $self->subqual($start,$end),
770 -trace_a
=> $self->sub_trace('a',$start2,$end2),
771 -trace_t
=> $self->sub_trace('t',$start2,$end2),
772 -trace_g
=> $self->sub_trace('g',$start2,$end2),
773 -trace_c
=> $self->sub_trace('c',$start2,$end2),
774 -peak_indices
=> $self->sub_peak_index($start,$end)
777 $new_object->set_accuracies();
778 $new_object->_reset_peak_indices();
782 =head2 _synthesize_traces()
784 Title : _synthesize_traces()
785 Usage : $obj->_synthesize_traces();
786 Function: Synthesize false traces for this object.
789 Notes : This method is intended to be invoked when this
790 object is created with a SWQ object- that is to say that
791 there is a sequence and a set of qualities but there was
792 no actual trace data.
796 sub _synthesize_traces
{
798 $self->peak_indices(qw());
800 # the user should be warned if traces already exist
803 #ml ( my $sequence = $self->seq() ) =~ tr/a-z/A-Z/;
804 #ml my @quals = @{$self->qual()};
806 # build the ramp for the first base.
807 # a ramp looks like this "1 4 13 29 51 71 80 71 51 29 13 4 1" times the quality score.
809 # note to self-> smooth this thing out a bit later
811 @
{$ramp_data->{'ramp'}} = qw( 1 4 13 29 51 75 80 75 51 29 13 4 1 );
812 # the width of the ramp
813 $ramp_data->{'ramp_width'} = scalar(@
{$ramp_data->{'ramp'}});
814 # how far should the peaks overlap?
815 $ramp_data->{'ramp_overlap'} = 1;
816 # where should the peaks be located?
817 $ramp_data->{'peak_at'} = 7;
818 $ramp_data->{'ramp_total_length'} =
819 $self->seq_obj()->length() * $ramp_data->{'ramp_width'}
820 - $self->seq_obj()->length() * $ramp_data->{'ramp_overlap'};
822 my $total_length = $ramp_data->{ramp_total_length
};
823 $self->initialize_traces("0",$total_length+2);
825 my ($current_base,$place_base_at,$peak_quality,$ramp_counter,$current_ramp,$ramp_position);
826 #ml my $sequence_length = $self->length();
827 my $half_ramp = int($ramp_data->{'ramp_width'}/2);
828 for ($pos = 0; $pos<$self->length();$pos++) {
829 $current_base = uc $self->seq_obj()->subseq($pos+1,$pos+1);
830 # print("Synthesizing the ramp for $current_base\n");
831 my $all_bases = "ATGC";
832 $peak_quality = $self->qual_obj()->qualat($pos+1);
833 # where should the peak for this base be placed? Modeled after a mktrace scf
834 $place_base_at = ($pos * $ramp_data->{'ramp_width'}) -
835 ($pos * $ramp_data->{'ramp_overlap'}) -
836 $half_ramp + $ramp_data->{'ramp_width'} - 1;
837 # print("Placing this base at this position: $place_base_at\n");
838 push @
{$self->peak_indices()},$place_base_at;
839 $ramp_position = $place_base_at - $half_ramp;
840 if ($current_base =~ "N" ) {
843 for ($current_ramp = 0; $current_ramp < $ramp_data->{'ramp_width'}; $current_ramp++) {
844 # print("Placing a trace value here: $current_base ".($ramp_position+$current_ramp+1)." ".$peak_quality*$ramp_data->{'ramp'}->[$current_ramp]."\n");
845 $self->trace_value_at($current_base,$ramp_position+$current_ramp+1,$peak_quality*$ramp_data->{'ramp'}->[$current_ramp]);
847 $self->peak_index_at($pos+1,
850 #ml my $other_bases = $self->_get_other_bases($current_base);
851 # foreach ( split('',$other_bases) ) {
852 # push @{$self->{'text'}->{"v3_base_accuracy"}->{$_}},0;
861 =head2 _dump_traces($transformed)
863 Title : _dump_traces("transformed")
864 Usage : &_dump_traces($ra,$rc,$rg,$rt);
865 Function: Used in debugging. Prints all traces one beside each other.
867 Args : References to the arrays containing the traces for A,C,G,T.
868 Notes : Beats using dumpValue, I'll tell ya. Much better then using
870 - if a scalar is included as an argument (any scalar), this
871 procedure will dump the _delta'd trace. If you don't know what
872 that means you should not be using this.
879 my (@sA,@sT,@sG,@sC);
880 print ("Count\ta\tc\tg\tt\n");
881 my $length = $self->trace_length();
882 for (my $curr=1; $curr <= $length; $curr++) {
883 print(($curr-1)."\t".$self->trace_value_at('a',$curr).
884 "\t".$self->trace_value_at('c',$curr).
885 "\t".$self->trace_value_at('g',$curr).
886 "\t".$self->trace_value_at('t',$curr)."\n");
891 =head2 _initialize_traces()
893 Title : _initialize_traces()
894 Usage : $trace_object->_initialize_traces();
895 Function: Creates empty arrays to hold synthetic trace values.
901 sub initialize_traces
{
902 my ($self,$value,$length) = @_;
903 foreach (qw(a t g c)) {
905 for (my $count=0; $count<$length; $count++) {
906 $temp[$count] = $value;
908 $self->trace($_,\
@temp);
912 =head2 trace_value_at($channel,$position)
914 Title : trace_value_at($channel,$position)
915 Usage : $value = $trace_object->trace_value_at($channel,$position);
916 Function: What is the value of the trace for this base at this position?
917 Returns : A scalar represnting the trace value here.
918 Args : a base channel (a,t,g,c)
919 a position ( < $trace_object->trace_length() )
924 my ($self,$channel,$position,$value) = @_;
926 $self->trace($channel)->[$position] = $value;
928 return $self->sub_trace($channel,($position),($position))->[0];
931 sub _deprecated_get_scf_version_2_base_structure
{
932 # this sub is deprecated- check inside SeqIO::scf
934 my (@structure,$current);
935 my $length = $self->length();
936 for ($current=1; $current <= $self->length() ; $current++) {
937 my $base_here = $self->seq_obj()->subseq($current,$current);
938 $base_here = lc($base_here);
940 $probabilities->{$base_here} = $self->qual_obj()->qualat($current);
941 my $other_bases = "atgc";
943 $other_bases =~ s/$base_here/$empty/e;
944 foreach ( split('',$other_bases) ) {
945 $probabilities->{$_} = "0";
949 $self->peak_index_at($current),
950 $probabilities->{'a'},
951 $probabilities->{'t'},
952 $probabilities->{'g'},
953 $probabilities->{'c'}
960 sub _deprecated_get_scf_version_3_base_structure
{
963 $structure = join('',$self->peak_indices());
968 =head2 accuracies($channel,$position)
970 Title : trace_value_at($channel,$position)
971 Usage : $value = $trace_object->trace_value_at($channel,$position);
972 Function: What is the value of the trace for this base at this position?
973 Returns : A scalar represnting the trace value here.
974 Args : a base channel (a,t,g,c)
975 a position ( < $trace_object->trace_length() )
981 my ($self,$channel,$value) = @_;
983 if (ref($value) eq "ARRAY") {
984 $self->{accuracies
}->{$channel} = $value;
987 my @acc = split(' ',$value);
988 $self->{accuracies
}->{$channel} = \
@acc;
991 return $self->{accuracies
}->{$channel};
995 =head2 set_accuracies()
997 Title : set_sccuracies()
998 Usage : $trace_object->set_accuracies();
999 Function: Take a sequence's quality and synthesize proper scf-style
1000 base accuracies that can then be accessed with
1001 accuracies("a") or something like it.
1007 sub set_accuracies
{
1010 my $length = $self->length();
1011 for ($count=1; $count <= $length; $count++) {
1012 my $base_here = $self->seq_obj()->subseq($count,$count);
1013 my $qual_here = $self->qual_obj()->qualat($count);
1014 $self->accuracy_at($base_here,$count,$qual_here);
1015 my $other_bases = $self->_get_other_bases($base_here);
1016 foreach (split('',$other_bases)) {
1017 $self->accuracy_at($_,$count,"null");
1026 Usage : $trace_object->scf_dump();
1027 Function: Prints out the contents of the structures representing
1028 the SequenceTrace in a manner similar to io_lib's scf_dump.
1029 Returns : Nothing. Prints out the contents of the structures
1030 used to represent the sequence and its trace.
1032 Notes : Used in debugging, obviously.
1039 for ($count=1;$count<=$self->length();$count++) {
1040 my $base_here = lc($self->seq_obj()->subseq($count,$count));
1041 print($base_here." ".sprintf("%05d",$self->peak_index_at($count))."\t");
1042 foreach (sort qw(a c g t)) {
1043 print(sprintf("%03d",$self->accuracy_at($_,$count))."\t");
1047 $self->_dump_traces();
1050 =head2 _get_other_bases($this_base)
1052 Title : _get_other_bases($this_base)
1053 Usage : $other_bases = $trace_object->_get_other_bases($this_base);
1054 Function: A utility routine to return bases other then the one provided.
1055 I was doing this over and over so I put it here.
1056 Returns : Three of a,t,g and c.
1057 Args : A base (atgc)
1058 Notes : $obj->_get_other_bases("a") returns "tgc"
1062 sub _get_other_bases
{
1063 my ($self,$this_base) = @_;
1064 $this_base = lc($this_base);
1065 my $all_bases = "atgc";
1067 $all_bases =~ s/$this_base/$empty/e;
1072 =head2 accuracy_at($base,$position)
1074 Title : accuracy_at($base,$position)
1075 Usage : $accuracy = $trace_object->accuracy_at($base,$position);
1077 Returns : Returns the accuracy of finding $base at $position.
1078 Args : 1. a base channel (atgc) 2. a value to _set_ the accuracy
1079 Notes : $obj->_get_other_bases("a") returns "tgc"
1085 my ($self,$base,$position,$value) = @_;
1088 if ($value eq "null") {
1089 $self->{accuracies
}->{$base}->[$position-1] = "0";
1092 $self->{accuracies
}->{$base}->[$position-1] = $value;
1095 return $self->{accuracies
}->{$base}->[$position-1];