Skip to content
Snippets Groups Projects
Commit 57c35d7e authored by Andreas Kusalananda Kähäri's avatar Andreas Kusalananda Kähäri
Browse files

Change to the method get_base_count():

  Requiring that all of $a, $c, $t, and $g be non-zero is too strict
  when trying to avoid that $a+$c+$t+g be non-zero in division.  Just
  test the sum instead.  This would have generated the wrong results for
  sequences which lacked one of A, C, T, or G...

  tr/Aa// is theoretically quicker than tr/Aa/Aa/

  Formatting


Tested.
parent f12295af
No related branches found
No related tags found
No related merge requests found
......@@ -610,38 +610,46 @@ sub subseq {
sub get_base_count {
my $self = shift;
my $a = 0; my $c = 0; my $t = 0; my $g = 0;
my $a = 0;
my $c = 0;
my $t = 0;
my $g = 0;
my $start = 1;
my $end;
my $RANGE = 100_000;
my $len = $self->length;
my $len = $self->length();
my $seq;
while($start <= $len) {
while ( $start <= $len ) {
$end = $start + $RANGE - 1;
$end = $len if ( $end > $len );
$end = $len if($end > $len);
$seq = $self->subseq($start, $end);
$a += $seq =~ tr/Aa/Aa/;
$c += $seq =~ tr/Cc/Cc/;
$t += $seq =~ tr/Tt/Tt/;
$g += $seq =~ tr/Gg/Gg/;
$seq = $self->subseq( $start, $end );
$a += $seq =~ tr/Aa//;
$c += $seq =~ tr/Cc//;
$t += $seq =~ tr/Tt//;
$g += $seq =~ tr/Gg//;
$start = $end + 1;
}
my $actg = $a + $c + $t + $g;
my $gc_content = 0;
if($a || $g || $c || $t) { #avoid divide by 0
$gc_content = sprintf( "%1.2f", (($g + $c)/($a + $g + $t + $c)) * 100);
if ( $actg > 0 ) { # Avoid dividing by 0
$gc_content = sprintf( "%1.2f", ( ( $g + $c )/$actg )*100 );
}
return {'a' => $a,
'c' => $c,
't' => $t,
'g' => $g,
'n' => $len - $a - $c - $t - $g,
'%gc' => $gc_content};
return { 'a' => $a,
'c' => $c,
't' => $t,
'g' => $g,
'n' => $len - $actg,
'%gc' => $gc_content };
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment