File Coverage

Bio/Tools/Alignment/Trim.pm
Criterion Covered Total %
statement 26 196 13.2
branch 0 64 0.0
condition 2 30 6.6
subroutine 9 21 42.8
pod 11 11 100.0
total 48 322 14.9


line stmt bran cond sub pod time code
1             # Bio::Tools::Alignment::Trim.pm
2             #
3             # Please direct questions and support issues to
4             #
5             # Cared for by Chad Matsalla
6             #
7             # Copyright Chad Matsalla
8             #
9             # You may distribute this module under the same terms as perl itself
10              
11             # POD documentation - main docs before the code
12              
13             =head1 NAME
14              
15             Bio::Tools::Alignment::Trim - A kludge to do specialized trimming of
16             sequence based on quality.
17              
18             =head1 SYNOPSIS
19              
20             use Bio::Tools::Alignment::Trim;
21             $o_trim = Bio::Tools::Alignment::Trim->new();
22             $o_trim->set_reverse_designator("R");
23             $o_trim->set_forward_designator("F");
24              
25              
26             =head1 DESCRIPTION
27              
28             This is a specialized module designed by Chad for Chad to trim sequences
29             based on a highly specialized list of requirements. In other words, write
30             something that will trim sequences 'just like the people in the lab would
31             do manually'.
32              
33             I settled on a sliding-window-average style of search which is ugly and
34             slow but does _exactly_ what I want it to do.
35              
36             Mental note: rewrite this.
37              
38             It is very important to keep in mind the context in which this module was
39             written: strictly to support the projects for which Consed.pm was
40             designed.
41              
42             =head1 FEEDBACK
43              
44             =head2 Mailing Lists
45              
46             User feedback is an integral part of the evolution of this and other
47             Bioperl modules. Send your comments and suggestions preferably to one
48             of the Bioperl mailing lists. Your participation is much appreciated.
49              
50             bioperl-l@bioperl.org - General discussion
51             http://bioperl.org/wiki/Mailing_lists - About the mailing
52             lists
53              
54             =head2 Support
55              
56             Please direct usage questions or support issues to the mailing list:
57              
58             I
59              
60             rather than to the module maintainer directly. Many experienced and
61             reponsive experts will be able look at the problem and quickly
62             address it. Please include a thorough description of the problem
63             with code and data examples if at all possible.
64              
65             =head2 Reporting Bugs
66              
67             Report bugs to the Bioperl bug tracking system to help us keep track
68             the bugs and their resolution. Bug reports can be submitted via the
69             web:
70              
71             https://github.com/bioperl/bioperl-live/issues
72              
73             =head1 AUTHOR - Chad Matsalla
74              
75             Email bioinformatics-at-dieselwurks.com
76              
77             =head1 APPENDIX
78              
79             The rest of the documentation details each of the object methods.
80             Internal methods are usually preceded with a _
81              
82             =cut
83              
84             package Bio::Tools::Alignment::Trim;
85              
86 1     1   3 use strict;
  1         1  
  1         22  
87 1     1   3 use Dumpvalue;
  1         1  
  1         15  
88              
89 1     1   3 use vars qw(%DEFAULTS);
  1         1  
  1         58  
90              
91 1     1   6 use base qw(Bio::Root::Root);
  1         1  
  1         310  
92              
93             BEGIN {
94 1     1   1261 %DEFAULTS = ( 'f_designator' => 'f',
95             'r_designator' => 'r',
96             'windowsize' => '10',
97             'phreds' => '20');
98             }
99              
100             =head2 new()
101              
102             Title : new()
103             Usage : $o_trim = Bio::Tools::Alignment::Trim->new();
104             Function: Construct the Bio::Tools::Alignment::Trim object. No parameters
105             are required to create this object. It is strictly a bundle of
106             functions, as far as I am concerned.
107             Returns : A reference to a Bio::Tools::Alignment::Trim object.
108             Args : (optional)
109             -windowsize (default 10)
110             -phreds (default 20)
111              
112              
113             =cut
114              
115             sub new {
116 1     1 1 1 my ($class,@args) = @_;
117 1         4 my $self = $class->SUPER::new(@args);
118 1         5 my($windowsize,$phreds) =
119             $self->_rearrange([qw(
120             WINDOWSIZE
121             PHREDS
122             )],
123             @args);
124 1   33     6 $self->{windowsize} = $windowsize || $DEFAULTS{'windowsize'};
125 1   33     5 $self->{phreds} = $phreds || $DEFAULTS{'phreds'};
126             # print("Constructor set phreds to ".$self->{phreds}."\n") if $self->verbose > 0;
127             $self->set_designators($DEFAULTS{'f_designator'},
128 1         3 $DEFAULTS{'r_designator'});
129 1         3 return $self;
130             }
131              
132             =head2 set_designators($forward_designator,$reverse_designator)
133              
134             Title : set_designators(,)
135             Usage : $o_trim->set_designators("F","R")
136             Function: Set the string by which the system determines whether a given
137             sequence represents a forward or a reverse read.
138             Returns : Nothing.
139             Args : two scalars: one representing the forward designator and one
140             representing the reverse designator
141              
142             =cut
143              
144             sub set_designators {
145 2     2 1 4 my $self = shift;
146 2         6 ($self->{'f_designator'},$self->{'r_designator'}) = @_;
147             }
148              
149             =head2 set_forward_designator($designator)
150              
151             Title : set_forward_designator($designator)
152             Usage : $o_trim->set_forward_designator("F")
153             Function: Set the string by which the system determines if a given
154             sequence is a forward read.
155             Returns : Nothing.
156             Args : A string representing the forward designator of this project.
157              
158             =cut
159              
160             sub set_forward_designator {
161 1     1 1 1 my ($self,$desig) = @_;
162 1         2 $self->{'f_designator'} = $desig;
163             }
164              
165             =head2 set_reverse_designator($reverse_designator)
166              
167             Title : set_reverse_designator($reverse_designator)
168             Function: Set the string by which the system determines if a given
169             sequence is a reverse read.
170             Usage : $o_trim->set_reverse_designator("R")
171             Returns : Nothing.
172             Args : A string representing the forward designator of this project.
173              
174             =cut
175              
176             sub set_reverse_designator {
177 1     1 1 1 my ($self,$desig) = @_;
178 1         2 $self->{'r_designator'} = $desig;
179             }
180              
181             =head2 get_designators()
182              
183             Title : get_designators()
184             Usage : $o_trim->get_designators()
185             Returns : A string describing the current designators.
186             Args : None
187             Notes : Really for informational purposes only. Duh.
188              
189             =cut
190              
191             sub get_designators {
192 0     0 1   my $self = shift;
193 0           return("forward: ".$self->{'f_designator'}." reverse: ".$self->{'r_designator'});
194             }
195              
196             =head2 trim_leading_polys()
197              
198             Title : trim_leading_polys()
199             Usage : $o_trim->trim_leading_polys()
200             Function: Not implemented. Does nothing.
201             Returns : Nothing.
202             Args : None.
203             Notes : This function is not implemented. Part of something I wanted to
204             do but never got around to doing.
205              
206             =cut
207              
208             sub trim_leading_polys {
209 0     0 1   my ($self, $sequence) = @_;
210             }
211              
212             =head2 dump_hash()
213              
214             Title : dump_hash()
215             Usage : $o_trim->dump_hash()
216             Function: Unimplemented.
217             Returns : Nothing.
218             Args : None.
219             Notes : Does nothing.
220              
221             =cut
222              
223             sub dump_hash {
224 0     0 1   my $self = shift;
225 0           my %hash = %{$self->{'qualities'}};
  0            
226             } # end dump_hash
227              
228             =head2 trim_singlet($sequence,$quality,$name,$class)
229              
230             Title : trim_singlet($sequence,$quality,$name,$class)
231             Usage : ($r_trim_points,$trimmed_sequence) =
232             @{$o_trim->trim_singlet($sequence,$quality,$name,$class)};
233             Function: Trim a singlet based on its quality.
234             Returns : a reference to an array containing the forward and reverse
235             trim points and the trimmed sequence.
236             Args : $sequence : A sequence (SCALAR, please)
237             $quality : A _scalar_ of space-delimited quality values.
238             $name : the name of the sequence
239             $class : The class of the sequence. One of qw(singlet
240             singleton doublet pair multiplet)
241             Notes : At the time this was written the bioperl objects SeqWithQuality
242             and PrimaryQual did not exist. This is what is with the clumsy
243             passing of references and so on. I will rewrite this next time I
244             have to work with it. I also wasn't sure whether this function
245             should return just the trim points or the points and the sequence.
246             I decided that I always wanted both so that's how I implemented
247             it.
248             - Note that the size of the sliding windows is set during construction of
249             the Bio::Tools::Alignment::Trim object.
250              
251             =cut
252              
253             sub trim_singlet {
254 0     0 1   my ($self,$sequence,$quality,$name,$class) = @_;
255             # this split is done because I normally store quality values in a
256             # space-delimited scalar rather then in an array.
257             # I do this because serialization of the arrays is tough.
258 0           my @qual = split(' ',$quality);
259 0           my @points;
260 0           my $sequence_length = length($sequence);
261 0           my ($returnstring,$processed_sequence);
262             # smooth out the qualities
263 0           my $r_windows = &_sliding_window(\@qual,$self->{windowsize});
264             # find out the leading and trailing trimpoints
265 0           my $start_base = $self->_get_start($r_windows,$self->{windowsize},$self->{phreds});
266 0           my (@new_points,$trimmed_sequence);
267             # do you think that any sequence shorter then 100 should be
268             # discarded? I don't think that this should be the decision of this
269             # module.
270             # removed, 020926
271 0           $points[0] = $start_base;
272             # whew! now for the end base
273             # required parameters: reference_to_windows,windowsize,$phredvalue,start_base
274             my $end_base = &_get_end($r_windows,$self->{windowsize},
275 0           $self->{phreds},$start_base);
276 0           $points[1] = $end_base;
277             # now do the actual trimming
278             # CHAD : I don't think that it is a good idea to call chop_sequence here
279             # because chop_sequence also removes X's and N's and things
280             # and that is not always what is wanted
281 0           return \@points;
282             }
283              
284             =head2 trim_doublet($sequence,$quality,$name,$class)
285              
286             Title : trim_doublet($sequence,$quality,$name,$class)
287             Usage : ($r_trim_points,$trimmed_sequence) =
288             @{$o_trim->trim_singlet($sequence,$quality,$name,$class)};
289             Function: Trim a singlet based on its quality.
290             Returns : a reference to an array containing the forward and reverse
291             Args : $sequence : A sequence
292             $quality : A _scalar_ of space-delimited quality values.
293             $name : the name of the sequence
294             $class : The class of the sequence. One of qw(singlet
295             singleton doublet pair multiplet)
296             Notes : At the time this was written the bioperl objects SeqWithQuality
297             and PrimaryQual did not exist. This is what is with the clumsy
298             passing of references and so on. I will rewrite this next time I
299             have to work with it. I also wasn't sure whether this function
300             should return just the trim points or the points and the sequence.
301             I decided that I always wanted both so that's how I implemented
302             it.
303              
304             =cut
305              
306             #'
307             sub trim_doublet {
308 0     0 1   my ($self,$sequence,$quality,$name,$class) = @_;
309 0           my @qual = split(' ',$quality);
310 0           my @points;
311 0           my $sequence_length = length($sequence);
312 0           my ($returnstring,$processed_sequence);
313             # smooth out the qualities
314 0           my $r_windows = &_sliding_window(\@qual,$self->{windowsize});
315             # determine where the consensus sequence starts
316 0           my $offset = 0;
317 0           for (my $current = 0; $current<$sequence_length;$current++) {
318 0 0         if ($qual[$current] != 0) {
319 0           $offset = $current;
320 0           last;
321             }
322             }
323             # start_base required: r_quality,$windowsize,$phredvalue
324 0           my $start_base = $self->_get_start($r_windows,$self->{windowsize},$self->{phreds},$offset);
325 0 0         if ($start_base > ($sequence_length - 100)) {
326 0           $points[0] = ("FAILED");
327 0           $points[1] = ("FAILED");
328 0           return @points;
329             }
330 0           $points[0] = $start_base;
331             #
332             # whew! now for the end base
333             #
334             # required parameters: reference_to_windows,windowsize,$phredvalue,start_base
335             # |
336             # 010420 NOTE: We will no longer get the end base to avoid the Q/--\___/-- syndrome
337 0           my $end_base = $sequence_length;
338 0           my $start_of_trailing_zeros = &count_doublet_trailing_zeros(\@qual);
339 0           $points[1] = $end_base;
340             # CHAD : I don't think that it is a good idea to call chop_sequence here
341             # because chop_sequence also removes X's and N's and things
342             # and that is not always what is wanted
343 0           return @points;
344             } # end trim_doublet
345              
346             =head2 chop_sequence($name,$class,$sequence,@points)
347              
348             Title : chop_sequence($name,$class,$sequence,@points)
349             Usage : ($start_point,$end_point,$chopped_sequence) =
350             $o_trim->chop_sequence($name,$class,$sequence,@points);
351             Function: Chop a sequence based on its name, class, and sequence.
352             Returns : an array containing three scalars:
353             1- the start trim point
354             2- the end trim point
355             3- the chopped sequence
356             Args :
357             $name : the name of the sequence
358             $class : The class of the sequence. One of qw(singlet
359             singleton doublet pair multiplet)
360             $sequence : A sequence
361             @points : An array containing two elements- the first contains
362             the start trim point and the second conatines the end trim
363             point.
364              
365             =cut
366              
367             sub chop_sequence {
368 0     0 1   my ($self,$name,$class,$sequence,@points) = @_;
369 0           print("Coming into chop_sequence, \@points are @points\n");
370 0           my $fdesig = $self->{'f_designator'};
371 0           my $rdesig = $self->{'r_designator'};
372 0 0 0       if (!$points[0] && !$points[1]) {
373 0           $sequence = "junk";
374 0           return $sequence;
375             }
376 0 0 0       if ($class eq "singlet" && $name =~ /$fdesig$/) {
    0 0        
    0 0        
    0 0        
    0          
377 0           $sequence = substr($sequence,$points[0],$points[1]-$points[0]);
378             }
379             elsif ($class eq "singlet" && $name =~ /$rdesig$/) {
380 0           $sequence = substr($sequence,$points[0],$points[1]-$points[0]);
381             }
382             elsif ($class eq "singleton" && $name =~ /$fdesig$/) {
383 0           $sequence = substr($sequence,$points[0],$points[1]-$points[0]);
384             }
385             elsif ($class eq "singleton" && $name =~ /$rdesig$/) {
386 0           $sequence = substr($sequence,$points[0],$points[1]-$points[0]);
387             }
388             elsif ($class eq "doublet") {
389 0           $sequence = substr($sequence,$points[0],$points[1]-$points[0]);
390             }
391             # this is a _terrible_ to do this! i couldn't seem to find a better way
392             # i thought something like s/(^.*[Xx]{5,})//g; might work, but no go
393             # no time to find a fix!
394 0           my $length_before_trimming = length($sequence);
395 0           my $subs_Xs = $sequence =~ s/^.*[Xx]{5,}//g;
396 0 0         if ($subs_Xs) {
397 0           my $length_after_trimming = length($sequence);
398 0           my $number_Xs_trimmed = $length_before_trimming - $length_after_trimming;
399 0           $points[0] += $number_Xs_trimmed;
400             }
401 0           $length_before_trimming = length($sequence);
402 0           my $subs_Ns = $sequence =~ s/[Nn]{1,}$//g;
403 0 0         if ($subs_Ns) {
404 0           my $length_after_trimming = length($sequence);
405 0           my $number_Ns_trimmed = $length_before_trimming - $length_after_trimming;
406 0           $points[1] -= $number_Ns_trimmed;
407 0           $points[1] -= 1;
408             }
409 0           push @points,$sequence;
410 0           print("chop_sequence \@points are @points\n");
411 0           return @points;
412             }
413              
414             =head2 _get_start($r_quals,$windowsize,$phreds,$offset)
415              
416             Title : _get_start($r_quals,$windowsize,$phreds,$offset)
417             Usage : $start_base = $self->_get_start($r_windows,5,20);
418             Function: Provide the start trim point for this sequence.
419             Returns : a scalar representing the start of the sequence
420             Args :
421             $r_quals : A reference to an array containing quality values. In
422             context, this array of values has been smoothed by then
423             sliding window-look ahead algorithm.
424             $windowsize : The size of the window used when the sliding window
425             look-ahead average was calculated.
426             $phreds :
427             $offset :
428              
429             =cut
430              
431             sub _get_start {
432 0     0     my ($self,$r_quals,$windowsize,$phreds,$offset) = @_;
433 0 0         print("Using $phreds phreds\n") if $self->verbose > 0;
434             # this is to help determine whether the sequence is good at all
435 0           my @quals = @$r_quals;
436 0           my ($count,$count2,$qualsum);
437 0 0         if ($offset) { $count = $offset; } else { $count = 0; }
  0            
  0            
438             # search along the length of the sequence
439 0           for (; ($count+$windowsize) <= scalar(@quals); $count++) {
440             # sum all of the quality values in this window.
441 0           my $cumulative=0;
442 0           for($count2 = $count; $count2 < $count+$windowsize; $count2++) {
443 0 0         if (!$quals[$count2]) {
444             # print("Quals don't exist here!\n");
445             }
446             else {
447 0           $qualsum += $quals[$count2];
448             # print("Incremented qualsum to ($qualsum)\n");
449             }
450 0           $cumulative++;
451             }
452             # print("The sum of this window (starting at $count) is $qualsum. I counted $cumulative bases.\n");
453             # if the total of windowsize * phreds is
454 0 0 0       if ($qualsum && $qualsum >= $windowsize*$phreds) { return $count; }
  0            
455 0           $qualsum = 0;
456             }
457             # if ($count > scalar(@quals)-$windowsize) { return; }
458 0           return $count;
459             }
460              
461             =head2 _get_end($r_qual,$windowsize,$phreds,$count)
462              
463             Title : _get_end($r_qual,$windowsize,$phreds,$count)
464             Usage : my $end_base = &_get_end($r_windows,20,20,$start_base);
465             Function: Get the end trim point for this sequence.
466             Returns : A scalar representing the end trim point for this sequence.
467             Args :
468             $r_qual : A reference to an array containing quality values. In
469             context, this array of values has been smoothed by then
470             sliding window-look ahead algorithm.
471             $windowsize : The size of the window used when the sliding window
472             look-ahead average was calculated.
473             $phreds :
474             $count : Start looking for the end of the sequence here.
475              
476             =cut
477              
478             sub _get_end {
479 0     0     my ($r_qual,$windowsize,$phreds,$count) = @_;
480 0           my @quals = @$r_qual;
481 0           my $total_bases = scalar(@quals);
482 0           my ($count2,$qualsum,$end_of_quals,$bases_counted);
483 0 0         if (!$count) { $count=0; }
  0            
484 0           BASE: for (; $count < $total_bases; $count++) {
485 0           $bases_counted = 0;
486 0           $qualsum = 0;
487 0           POSITION: for($count2 = $count; $count2 < $total_bases; $count2++) {
488 0           $bases_counted++;
489              
490 0 0         if ($count2 == $total_bases-1) {
    0          
491 0           $qualsum += $quals[$count2];
492 0           $bases_counted++;
493 0           last BASE;
494             }
495             elsif ($bases_counted == $windowsize) {
496 0           $qualsum += $quals[$count2];
497 0 0         if ($qualsum < $bases_counted*$phreds) {
498 0           return $count+$bases_counted+$windowsize;
499             }
500 0           next BASE;
501             }
502             else {
503 0           $qualsum += $quals[$count2];
504             }
505             }
506 0 0         if ($qualsum < $bases_counted*$phreds) {
507 0           return $count+$bases_counted+$windowsize;
508             }
509             else { }
510 0           $qualsum = 0;
511             } # end for
512 0 0         if ($end_of_quals) {
513 0           my $bases_for_average = $total_bases-$count2;
514 0           return $count2;
515             }
516             else { }
517 0 0         if ($qualsum) { } # print ("$qualsum\n");
518 0           return $total_bases;
519             } # end get_end
520              
521             =head2 count_doublet_trailing_zeros($r_qual)
522              
523             Title : count_doublet_trailing_zeros($r_qual)
524             Usage : my $start_of_trailing_zeros = &count_doublet_trailing_zeros(\@qual);
525             Function: Find out when the trailing zero qualities start.
526             Returns : A scalar representing where the zeros start.
527             Args : A reference to an array of quality values.
528             Notes : Again, this should be rewritten to use PrimaryQual objects.
529             A more detailed explanation of why phrap puts these zeros here should
530             be written and placed here. Please email and hassle the author.
531              
532              
533             =cut
534              
535             sub count_doublet_trailing_zeros {
536 0     0 1   my ($r_qual) = shift;
537 0           my $number_of_trailing_zeros = 0;
538 0           my @qualities = @$r_qual;
539 0           for (my $current=scalar(@qualities);$current>0;$current--) {
540 0 0 0       if ($qualities[$current] && $qualities[$current] != 0) {
541 0           $number_of_trailing_zeros = scalar(@qualities)-$current;
542 0           return $current+1;
543             }
544             }
545 0           return scalar(@qualities);
546             } # end count_doublet_trailing_zeros
547              
548             =head2 _sliding_window($r_quals,$windowsize)
549              
550             Title : _sliding_window($r_quals,$windowsize)
551             Usage : my $r_windows = &_sliding_window(\@qual,$windowsize);
552             Function: Create a sliding window, look-forward-average on an array
553             of quality values. Used to smooth out differences in qualities.
554             Returns : A reference to an array containing the smoothed values.
555             Args : $r_quals: A reference to an array containing quality values.
556             $windowsize : The size of the sliding window.
557             Notes : This was written before PrimaryQual objects existed. They
558             should use that object but I haven't rewritten this yet.
559              
560             =cut
561              
562             #'
563             sub _sliding_window {
564 0     0     my ($r_quals,$windowsize) = @_;
565 0           my (@window,@quals,$qualsum,$count,$count2,$average,@averages,$bases_counted);
566 0           @quals = @$r_quals;
567 0           my $size_of_quality = scalar(@quals);
568             # do this loop for all of the qualities
569 0           for ($count=0; $count <= $size_of_quality; $count++) {
570 0           $bases_counted = 0;
571 0           BASE: for($count2 = $count; $count2 < $size_of_quality; $count2++) {
572 0           $bases_counted++;
573             # if the search hits the end of the averages, stop
574             # this is for the case near the end where bases remaining < windowsize
575 0 0         if ($count2 == $size_of_quality) {
    0          
576 0           $qualsum += $quals[$count2];
577 0           last BASE;
578             }
579             # if the search hits the size of the window
580             elsif ($bases_counted == $windowsize) {
581 0           $qualsum += $quals[$count2];
582 0           last BASE;
583             }
584             # otherwise add the quality value
585 0 0         unless (!$quals[$count2]) {
586 0           $qualsum += $quals[$count2];
587             }
588             }
589 0 0 0       unless (!$qualsum || !$windowsize) {
590 0           $average = $qualsum / $bases_counted;
591 0 0         if (!$average) { $average = "0"; }
  0            
592 0           push @averages,$average;
593             }
594 0           $qualsum = 0;
595             }
596             # 02101 Yes, I repaired the mismatching numbers between averages and windows.
597             # print("There are ".scalar(@$r_quals)." quality values. They are @$r_quals\n");
598             # print("There are ".scalar(@averages)." average values. They are @averages\n");
599 0           return \@averages;
600            
601             }
602              
603             =head2 _print_formatted_qualities
604              
605             Title : _print_formatted_qualities(\@quals)
606             Usage : &_print_formatted_qualities(\@quals);
607             Returns : Nothing. Prints.
608             Args : A reference to an array containing quality values.
609             Notes : An internal procedure used in debugging. Prints out an array nicely.
610              
611             =cut
612              
613             sub _print_formatted_qualities {
614 0     0     my $rquals = shift;
615 0           my @qual = @$rquals;
616 0           for (my $count=0; $count
617 0 0         if (($count%10)==0) { print("\n$count\t"); }
  0            
618 0 0         if ($qual[$count]) { print ("$qual[$count]\t");}
  0            
619 0           else { print("0\t"); }
620             }
621 0           print("\n");
622             }
623              
624             =head2 _get_end_old($r_qual,$windowsize,$phreds,$count)
625              
626             Title : _get_end_old($r_qual,$windowsize,$phreds,$count)
627             Usage : Deprecated. Don't use this!
628             Returns : Deprecated. Don't use this!
629             Args : Deprecated. Don't use this!
630              
631             =cut
632              
633             #'
634             sub _get_end_old {
635 0     0     my ($r_qual,$windowsize,$phreds,$count) = @_;
636 0           warn("Do Not Use this function (_get_end_old)");
637 0           my $target = $windowsize*$phreds;
638 0           my @quals = @$r_qual;
639 0           my $total_bases = scalar(@quals);
640 0           my ($count2,$qualsum,$end_of_quals);
641 0 0         if (!$count) { $count=0; }
  0            
642 0           BASE: for (; $count < $total_bases; $count++) {
643 0           for($count2 = $count; $count2 < $count+$windowsize; $count2++) {
644 0 0         if ($count2 == scalar(@quals)-1) {
645 0           $qualsum += $quals[$count2];
646 0           $end_of_quals = 1;
647 0           last BASE;
648              
649             }
650 0           $qualsum += $quals[$count2];
651             }
652 0 0         if ($qualsum < $windowsize*$phreds) {
653 0           return $count+$windowsize;
654             }
655 0           $qualsum = 0;
656             } # end for
657             } # end get_end_old
658              
659              
660             # Autoload methods go after =cut, and are processed by the autosplit program.
661              
662             1;
663             __END__