G::Seq::GCskew - Analysis methods related to GC skew and genomic strand bias

    This class is a part of G-language Genome Analysis Environment, 
    collecting sequence analysis methods related to GC skew.

    msg_error("WARNING: cum_gcskew is deprecated since v.1.6.13.\n" .
	      "         This method will be removed in future releases.\n" .
	      "         use gcskew(\$gb, -cumulative=>1) instead!\n\n");

    return gcskew(@_, -cumulative=>1);

    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);
    my $first = shift @args;
    my $second = shift @args;

    $first = $gb->{$first}->{start} if($first =~ /^FEATURE/ || /^CDS/);
    $second = $gb->{$second}->{start} if($second =~ /^FEATURE/ || /^CDS/);

    unless(length($second)){
	my ($ori, $ter) = rep_ori_ter($gb);
	$second = $ori;
    }

    my @dist;

    $dist[0] = abs($first - $second);
    $dist[1] = abs($first + length($gb->{SEQ}) - $second);
    $dist[2] = abs($first - (length($gb->{SEQ}) + $second));
    
    my @new = sort {$a <=> $b} @dist;

    return shift @new;

    require Statistics::Descriptive;
    &opt_default(output=>"stdout", purine=>0, keto=>0, window=>4096);
    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);
    my $len = length($gb->{SEQ});
    my $output = opt_val("output");
    my $purine = opt_val("purine");
    my $keto = opt_val("keto");
    my $lowpass = opt_val("filter");
    my $power = opt_val("window");

    &msg_send("\nfind_ori_ter:\n") if ($output eq 'stdout');

    my ($max, $min);
    if(length $lowpass){
	require Math::FFT;
	my $window = int($len / $power);
	my @gcskew = gcskew($gb, -output=>"/dev/null", -window=>$window, -purine=>$purine, -keto=>$keto);

	while(scalar @gcskew > $power){
	    pop @gcskew;
	}

	my $fft = new Math::FFT([@gcskew]);
	my $coeff = $fft->rdft();
	my $coeff2 = [@{$coeff}];
	my $spctrm = $fft->spctrm();
	my $j;
	my $power2 = int($power/100);
	$power2 ++ if ($power2 %2 == 1);

	for($j = $power2 * (100 - $lowpass) - 1; $j < $power; $j ++){
	    $coeff->[$j] = 0;
	}
    
	my $orig = $fft->invrdft($coeff);

	my (@cum, $tmp);
	foreach my $value (@{$orig}){
	    $tmp += $value;
	    push(@cum, $tmp);
	}

	my $stat = Statistics::Descriptive::Full->new();
	$stat->add_data(@cum);
	my $maxi = $stat->maxdex();
	my $mini = $stat->mindex();

	$max = ($maxi + 1) * $window - int($window/2);
	$min = ($mini + 1) * $window - int($window/2);
    }else{
	local *peak_search = sub{
	    my $seq = shift;

	    if($purine){
		$seq =~ tr/atgcn/02021/;
	    }elsif($keto){
		$seq =~ tr/atgcn/20021/;
	    }else{
		$seq =~ tr/atgcn/11021/;
	    }
	    
	    my (@data, $val, $i);
	    for($i = 0; $i <= length($seq); $i ++){
		if(substr($seq, $i, 1) =~ /^\d$/){
		    $val += substr($seq, $i, 1) - 1;
		}
		push(@data, $val);
	    }

	    my $stat = Statistics::Descriptive::Full->new();
	    $stat->add_data(@data);
	    return ($stat->maxdex(), $min = $stat->mindex());
	};

	if($len > 100000){
	    my @cumgcskew = gcskew($gb, -output=>"/dev/null", -window=>10000, 
				   -purine=>$purine, -keto=>$keto, -cumulative=>1);

	    my $stat = Statistics::Descriptive::Full->new();
	    $stat->add_data(@cumgcskew);
	    my $maxi = $stat->maxdex();
	    my $mini = $stat->mindex();

	    ($max, undef) = peak_search(substr($gb->{SEQ}, $maxi * 10000, 20000));
	    (undef, $min) = peak_search(substr($gb->{SEQ}, $mini * 10000, 20000));
	    $max += $maxi * 10000;
	    $min += $mini * 10000;
	    $max = 0 if ($max < $len / 1000 || abs($max - $len) < $len / 1000);
	}else{
	    ($max, $min) = peak_search($gb->{SEQ});
	}
    }

    &msg_send("   Predicted Origin:   " , $max, "\n")  if ($output eq 'stdout');
    &msg_send("   Predicted Terminus: " , $min, "\n\n") if ($output eq 'stdout');

    return ($max, $min);

    require Math::FFT;

    opt_default(window=>4096, purine=>0, keto=>0, at=>0, version=>2, p=>0);
    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);
    my %opt = opt_val();

    my $window = int(length($gb->{SEQ}) / $opt{'window'});
    if($window < 10){
	die("Error in gcsi: number of windows is too large, or the genome is too small.\n" . 
	    "GCSI would not be accurate for window size less than 100 bp.\n");
    }
    my @gcskew = gcskew($gb, -output=>'/dev/null', -window=>$window, -purine=>$opt{'purine'}, -keto=>$opt{'keto'}, -at=>$opt{'at'});

    while(scalar @gcskew > $opt{'window'}){
	pop @gcskew;
    }

    local *gcsicore = sub{
	my @gcskew = @_;
	my @cumgcskew = cumulative(@gcskew);

	my $fft = new Math::FFT([@gcskew]);
	my @spectrum = @{$fft->spctrm()};
	shift @spectrum;
	my $first = shift @spectrum;

	my $dist = abs(max(@cumgcskew)) + abs(min(@cumgcskew));
	my ($gcsi, $sr);

	if($opt{'version'} == 1){
	    $sr = $first/mean(@spectrum);
	    $gcsi = ($sr/6000 + $dist/600)/2;
	}else{
	    $dist *= 4096/$opt{'window'};
	    $sr = 40 * (($first * 6000 * 100) ** 0.4);
	    $gcsi = sqrt($sr/6000 * $dist/600);
	}

	return ($gcsi, $sr, $dist);
    };

    my ($gcsi, $sr, $dist) = gcsicore(@gcskew);

    if($opt{'p'}){
	require Statistics::Distributions;

	my @random;
	for (1..100){
	    my ($rgcsi, undef, undef) = gcsicore(shuffle(@gcskew));
	    push(@random, $rgcsi);
	}

	my $z = abs(($gcsi - mean(@random))/standard_deviation(@random));
	my $p = Statistics::Distributions::uprob($z) * 2;
#	my $rcmd = new Rcmd();
#	say $rcmd->normtest(@random);
	return ($gcsi, $sr, $dist, $z, $p);
    }

    if(wantarray()){
	return ($gcsi, $sr, $dist);
    }else{
	return $gcsi;
    }

    &opt_default(window=>10000, at=>0, purine=>0, keto=>0, output=>"show", 
		 filename=>"gcskew.png", cumulative=>0, slide=>undef);
    my @args       = opt_get(@_);
    my $gb         = opt_as_gb(shift @args);
    my $ref        =\$ gb->{SEQ};
    my $len        = length($$ref);
    my $window     = opt_val("window");
    my $slide      = opt_val("slide") || $window;
    die('Error at gcskew(): window size too small') if ($window < 10);
    my $output     = opt_val("output");
    my $cumulative = '';
    $cumulative    = "Cumulative " if (opt_val("cumulative"));
    my $filename   = opt_val("filename");
    $filename      =~ s/\.png$/\.csv/ if (opt_val("output") eq 'f');
    my $at         = opt_val("at");
    my $purine     = opt_val("purine");
    my $keto       = opt_val("keto");
    my @gcskew     = ();
    my @location   = ();
    my ($tmp, $pos, $j, $i) = (0,0,0,0);

    while($len - $pos >= $window){
	my ($g, $c);

	if($at){
	    $g =  substr($$ref, $pos, $window) =~ tr/a/a/;
	    $c =  substr($$ref, $pos, $window) =~ tr/t/t/;
	}elsif($purine){
	    $g =  substr($$ref, $pos, $window) =~ tr/a/a/;
	    $g += substr($$ref, $pos, $window) =~ tr/g/g/;
	    $c =  substr($$ref, $pos, $window) =~ tr/t/t/;
	    $c += substr($$ref, $pos, $window) =~ tr/c/c/;
	}elsif($keto){
	    $g =  substr($$ref, $pos, $window) =~ tr/t/t/;
	    $g += substr($$ref, $pos, $window) =~ tr/g/g/;
	    $c =  substr($$ref, $pos, $window) =~ tr/a/a/;
	    $c += substr($$ref, $pos, $window) =~ tr/c/c/;
	}else{
	    $g =  substr($$ref, $pos, $window) =~ tr/g/g/;
	    $c =  substr($$ref, $pos, $window) =~ tr/c/c/;
	}

	if(length($cumulative)){
	    if ($c+$g <= 0){
		$tmp += 0;
	    }else{
		$tmp += sprintf("%.6f",($c-$g)/($c+$g));
	    }
	    $gcskew[$i] = $tmp;
	}else{
	    if ($c+$g <= 0){
		$gcskew[$i] = 0;
	    }else{
		$gcskew[$i] = sprintf("%.6f",($c-$g)/($c+$g));
	    }
	}

	$location[$i] = $pos;
	$pos += $slide;
	$i ++;
    }
    $i --;

    my $title = $cumulative . "GC skew";
    if ($at){
	$title = $cumulative . "AT skew";
    }elsif($purine){
	$title = $cumulative . "Purine skew";
    }elsif($keto){
	$title = $cumulative . "Keto skew";
    }

    if ($output eq 'g' || $output eq 'show'){
	mkdir ("graph", 0777);

	_UniMultiGrapher(\@
			 location,\@gcskew,
			 -x=>"bp", -y=>$title,
			 -filename=>$filename,
			 -title=>$title,
			 -style=>"lines", -type=>"columns",
			 );

	msg_gimv("graph/" . $filename) if ($output eq 'show');
    }elsif ($output eq 'f'){
	my $j = 0;
	mkdir ("data", 0777);

	open(OUT, ">data/" . $filename);
	print OUT "location,$title\n";
	for ($j = 0; $j <= $i; $j++){
	    print OUT $location[$j], ",", $gcskew[$j], "\n";
	}
	close(OUT);
    }
    
    return @gcskew;

    &opt_default(window=>10000, at=>0, purine=>0, keto=>0, output=>"show", 
		 application=>"gimv", filename=>"gcwin.png", skipGap=>0);
    my @args = opt_get(@_);

    my $gb = opt_as_gb(shift @args);
    my $ref =\$ gb->{SEQ};
    my $window = opt_val("window");
    my $at = opt_val("at");
    my $purine = opt_val("purine");
    my $keto = opt_val("keto");
    my $skip = opt_val("skipGap");
    my $application = opt_val("application");
    my $filename = opt_val("filename");
    $filename =~ s/\.png$/\.csv/ if (opt_val("output") eq 'f');
    my $opt = opt_val("output");
    my (@gcwin, @location);
    my $j;
    
    my $i = 0;
    my ($g, $c, $rest);
    while(length($$ref) - ($window * $i) >= $window){
        if($at){
            $g = substr($$ref, $window * $i, $window) =~ tr/a/a/;
            $c = substr($$ref, $window * $i, $window) =~ tr/t/t/;
	    $rest =  substr($$ref, $window * $i, $window) =~ tr/g/g/;
	    $rest += substr($$ref, $window * $i, $window) =~ tr/c/c/;
        }elsif($purine){
            $g = substr($$ref, $window * $i, $window) =~ tr/g/g/;
            $c = substr($$ref, $window * $i, $window) =~ tr/a/a/;
	    $rest =  substr($$ref, $window * $i, $window) =~ tr/t/t/;
	    $rest += substr($$ref, $window * $i, $window) =~ tr/c/c/;
        }elsif($keto){
            $g = substr($$ref, $window * $i, $window) =~ tr/g/g/;
            $c = substr($$ref, $window * $i, $window) =~ tr/t/t/;
	    $rest =  substr($$ref, $window * $i, $window) =~ tr/a/a/;
	    $rest += substr($$ref, $window * $i, $window) =~ tr/c/c/;
	}else{
            $g = substr($$ref, $window * $i, $window) =~ tr/g/g/;
            $c = substr($$ref, $window * $i, $window) =~ tr/c/c/;
	    $rest =  substr($$ref, $window * $i, $window) =~ tr/a/a/;
	    $rest += substr($$ref, $window * $i, $window) =~ tr/t/t/;
        }
	$rest += $g + $c;

	if($rest == 0 && $skip){
	    $gcwin[$i] = $gcwin[$i - 1];
	}else{
	    $gcwin[$i] = $rest >= 1 ? sprintf("%.6f",($g+$c)/$rest) : 0;
	}
	$location[$i] = $i * $window;

	$i ++;
    }
    $i --;

    my $title = "GC content";
    if ($at){
        $title = "AT content";
    }elsif($purine){
        $title = "Purine content";
    }elsif($keto){
        $title = "Keto content";
    }
        
    if ($opt eq 'g' || $opt eq 'show'){
	mkdir ("graph", 0777);
	_UniMultiGrapher(\@
			 location,\@ gcwin,
			 -x=>"bp", -y=>$title, 
			 -filename=>$filename,
			 -title=>$title, -style=>"lines", -type=>"columns"
			 );
	msg_gimv("graph/" . $filename)
	    if ($opt eq 'show');;
    }elsif ($opt eq 'f'){
	my $j = 0;
	mkdir ("data", 0777);
	open(OUT, ">data/" . $filename);
	print OUT "location,$title\n";
	for ($j = 0; $j <= $i; $j++){
	    print OUT $location[$j], ",", $gcwin[$j], "\n";
	}
	close(OUT);
    }
    
    return @gcwin;

    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);
    my $wing = shift @args;

    my ($ori, $ter) = rep_ori_ter($gb);

    my (@sectionA, @sectionB, @sectionC, @left, @right);
    
    if($ori < $ter){
	foreach my $cds ($gb->cds()){

	    if($gb->{$cds}->{start} <= $ori){
		push(@sectionA, $cds);
	    }elsif($gb->{$cds}->{start} >= $ori && $gb->{$cds}->{start} <= $ter){
		push(@sectionB, $cds);
	    }elsif($gb->{$cds}->{start} >= $ter){
		push(@sectionC, $cds);
	    }else{
		warn("Something is wrong at G::Seq::GCskew::genes_from_ori()");
	    }
	}

	@left = (reverse(@sectionA), reverse(@sectionC));
	@right = @sectionB;

    }else{
	foreach my $cds ($gb->cds()){

	    if($gb->{$cds}->{start} < $ter){
		push(@sectionA, $cds);
	    }elsif($gb->{$cds}->{start} >= $ter && $gb->{$cds}->{start} <= $ori){
		push(@sectionB, $cds);
	    }elsif($gb->{$cds}->{start} >= $ori){
		push(@sectionC, $cds);
	    }else{
		warn("Something is wrong at G::Seq::GCskew::genes_from_ori()");
	    }
	}
	
	@left = (reverse(@sectionB));
	@right = (@sectionC, @sectionA);
	
    }

    if(lc($wing) =~ /l/){
	return @left;
    }elsif(lc($wing) =~ /r/){
	return @right;
    }else{
	my %hash;
	foreach my $cds ($gb->cds()){
	    $hash{$cds} = dist_in_cc($gb, $gb->{$cds}->{start});
	}

	my @all = sort{ $hash{$a} <=> $hash{$b}} keys %hash;
	return @all;
    }

    &opt_default(divide=>250, at=>0, output=>"show", application=>"gimv",
		  filename=>"genomicskew.png", intron=>0);
    my @args = opt_get(@_);

    my $filename = opt_val("filename");
    $filename =~ s/\.png$/\.csv/ if (opt_val("output") eq 'f');
	
    my $gb = opt_as_gb(shift @args);
    my $divide = opt_val("divide");
    my $opt = opt_val("output");
    my $application = opt_val("application");
    my $at = opt_val("at");
    my $intron = opt_val("intron");
    my (@gcskew, @betskew, @geneskew, @thirdskew);
    my @location = (0..$divide);
    my ($j, $window, $CDS, $BET, $THIRD);

    foreach my $cds ($gb->cds()){
	next if (length $gb->{$cds}->{join});
	my $seq .= $gb->get_gbkseq($gb->{$cds}->{start}, $gb->{$cds}->{end});
	$CDS .= $seq;
       
        for($j = 2; $j <= length($seq); $j += 3){
            if ($gb->{"$cds"}->{direction} eq 'complement'){
                $THIRD .= substr($seq, $j, 1);
            }else{
                $THIRD .= substr($seq, $j - 2, 1);
            }
        }
    }

    foreach my $cds ($gb->intergenic()){
	$BET .= $gb->get_geneseq($cds);
    }

    my $i = 0;
    
    $window = int(length($gb->{SEQ}) / $divide);
    while($i <= $divide){
	my $g = substr($gb->{SEQ}, $window * $i, $window) =~ tr/g/g/;
	$g = substr($gb->{SEQ}, $window * $i, $window) =~ tr/a/a/ if ($at);
	my $c = substr($gb->{SEQ}, $window * $i, $window) =~ tr/c/c/;
	$c = substr($gb->{SEQ}, $window * $i, $window) =~ tr/t/t/ if ($at);
	$gcskew[$i] = 0;
	$gcskew[$i] = sprintf("%.6f",($c-$g)/($c+$g)) unless ($c+$g<1);
	$i ++;
    }
    $i = 0;
    
    $window = int(length($CDS) / $divide);
    while($i <= $divide){
	my $g = substr($CDS, $window * $i, $window) =~ tr/g/g/;
	$g = substr($CDS, $window * $i, $window) =~ tr/a/a/ if ($at);
	my $c = substr($CDS, $window * $i, $window) =~ tr/c/c/;
	$c = substr($CDS, $window * $i, $window) =~ tr/t/t/ if ($at);
	$geneskew[$i] = 0;
	$geneskew[$i] = sprintf("%.6f",($c-$g)/($c+$g)) unless ($c+$g<1);
	$i ++;
    }
    $i = 0;
    
    $window = int(length($BET) / $divide);
    while($i <= $divide){
	my $g = substr($BET, $window * $i, $window) =~ tr/g/g/;
	$g = substr($BET, $window * $i, $window) =~ tr/a/a/ if ($at);
	my $c = substr($BET, $window * $i, $window) =~ tr/c/c/;
	$c = substr($BET, $window * $i, $window) =~ tr/t/t/ if ($at);
	$betskew[$i] = 0;
	$betskew[$i] = sprintf("%.6f",($c-$g)/($c+$g)) unless ($c+$g<1);
	$i ++;
    }
    $i = 0;
    
    $window = int(length($THIRD) / $divide);
    while($i <= $divide){
	my $g = substr($THIRD, $window * $i, $window) =~ tr/g/g/;
	$g = substr($THIRD, $window * $i, $window) =~ tr/a/a/ if ($at);
	my $c = substr($THIRD, $window * $i, $window) =~ tr/c/c/;
	$c = substr($THIRD, $window * $i, $window) =~ tr/t/t/ if ($at);
	$thirdskew[$i] = 0;
	$thirdskew[$i] = sprintf("%.6f",($c-$g)/($c+$g)) unless ($c+$g<1);
	$i ++;
    }

    my $title = "GC skew";
    $title = "AT skew" if ($at);
    
    if ($opt eq "show" || $opt eq "g"){
	mkdir ("graph", 0777);
	_UniMultiGrapher(\@
			 location, 
			 -x=>"bp", -y=>$title,\@
			 gcskew, -x1=>"whole genome",\@
			 geneskew, -x2=>"coding region",\@
			 betskew, -x3=>"intergenic region",\@
			 thirdskew, -x4=>"codon third position",
			 -style=>"lines", -type=>"columns",
			 -filename=>$filename,
			 -title=>$title
			 );

	msg_gimv("graph/" . $filename) if ($opt eq 'show');
    }elsif ($opt eq 'f'){
	my $j = 0;
	mkdir ("data", 0777);
	open(OUT, ">data/" . $filename);
	print OUT "location,$title,coding,intergenic,third codon\n";
	for ($j = 0; $j <= $divide; $j++){
	    print OUT $location[$j], ",", $gcskew[$j], ",", $geneskew[$j], ",",
	    $betskew[$j], ",", $thirdskew[$j], ",", "\n";
	}
	close(OUT);
    }
    
    return (\@gcskew,\@ geneskew,\@ betskew,\@ thirdskew);

    my @args = opt_get(@_);
    my $gb = shift @args;
    my ($ori, $ter) = rep_ori_ter($gb);
    my ($seq1, $seq2);

    if ($ori > $ter){
	$seq1  = substr($gb->{SEQ}, $ori);
	$seq1 .= substr($gb->{SEQ}, 0, $ter);
	$seq2  = complement(substr($gb->{SEQ}, $ter, $ori - $ter));
    }else{
	$seq1 = substr($gb->{SEQ}, $ori, $ter - $ori);
	$seq2 = complement( substr($gb->{SEQ}, $ter) . substr($gb->{SEQ}, 0, $ori) );
    }
    
    if(wantarray()){
	return ($seq1, $seq2);
    }else{
	return $seq1 . $seq2;
    }

    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);
    my $pos = shift @args;

    if($pos =~ /^FEATURE/ || /^CDS/){
	$pos = $gb->{$pos}->{start};
    }
    
    my ($ori, $ter) = rep_ori_ter($gb);

    if($ori < $ter){
	if($pos <= $ori){
	    return 'left';
	}elsif($pos >= $ori && $pos <= $ter){
	    return 'right';
	}elsif($pos >= $ter){
	    return 'left';
	}
    }else{
	if($pos < $ter){
	    return 'right';
	}elsif($pos >= $ter && $pos <= $ori){
	    return 'left';
	}elsif($pos >= $ori){
	    return 'right';
	}
    }

    opt_default(direction=>'direct');
    my @args = opt_get(@_);
    my $gb = shift @args;
    my $pos = shift @args;
    my $direction = opt_val("direction");

    if($pos =~ /^FEATURE/ || /^CDS/){
	$direction = $gb->{$pos}->{direction};
	$pos = $gb->{$pos}->{start};
    }

    my ($ori, $ter) = rep_ori_ter($gb);

    if ($ori > $ter){
	if ($pos < $ter || $pos > $ori){
	    if ($direction eq 'complement'){
		return ("lagging");
	    }else{
		return ("leading");
	    }
	}else{
	    if ($direction eq 'complement'){
		return ("leading");
	    }else{
		return ("lagging");
	    }
	}
    }else{
	if ($pos > $ori && $pos < $ter){
	    if ($direction eq 'complement'){
		return ("lagging");
	    }else{
		return ("leading");
	    }
	}else{
	    if ($direction eq 'complement'){
		return ("leading");
	    }else{
		return ("lagging");
	    }
	}
    }

    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);
    my ($ori, $ter);
    my $id = $gb->{LOCUS}->{id};


    if(length $gb->{FEATURE0}->{terminus}){
	$ori = $gb->{FEATURE0}->{origin};
	$ter = $gb->{FEATURE0}->{terminus};
    }else{

	if ($id eq 'U00096' || $id eq 'NC_000913'){
	    $ori = 3923500 - 1; 
	    $ter = 1588800 - 1; 
	}elsif ($id eq 'AL009126' || $id eq 'NC_000964'){
	    $ori = 1 - 1;
	    $ter = 2017000 - 1;
	}elsif ($id eq 'L42023' || $id eq 'NC_000907'){
	    $ori = 603000 - 1;
	    $ter = 1518000 - 1;
	}elsif ($id eq 'AL513382' || $id eq 'NC_003198'){
	    $ori = 3765000 - 1;
	    $ter = 1437000 - 1;
	}else{
	    msg_error('Using find_ori_ter() to predict...' . "\n");
	    ($ori, $ter) = &G::Seq::GCskew::find_ori_ter($gb, -output=>"/dev/null");
	}

	my $interface = msg_ask_interface();
	msg_interface('NULL');
	my $cds = ($gb->find(-type=>'rep_origin'))[0];
	if($cds){
	    $ori = int(($gb->{$cds}->{end} + $gb->{$cds}->{start})/2) - 1;
	    msg_error("rep_origin feature found. Using $ori as origin.\n");
	}
	msg_interface($interface);

	$gb->{FEATURE0}->{origin} = $ori;
	$gb->{FEATURE0}->{terminus} = $ter;
    }

    return ($ori, $ter);

    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);

    foreach my $cds ($gb->cds()){
	my $geneseq = $gb->get_geneseq($cds);

	my (%gc3, $tot);

	my $i = 0;
	for($i = 2; $i < length $geneseq; $i += 3){
	    $gc3{substr($geneseq, $i, 1)}++;
	    $tot ++;
	}

	$gb->{$cds}->{gc3} = ($gc3{g} + $gc3{c})/$tot;
    }

    return 1;

    my @args = opt_get(@_);
    my $gb = opt_as_gb(shift @args);

    foreach my $cds ($gb->cds()){
	$gb->{$cds}->{strand} = query_strand($gb, $cds);
	$gb->{$cds}->{arm} = query_arm($gb, $cds);
    }

    return 1;