#!/usr/bin/perl -w
package CodonBias;
require Exporter;

use strict;
use lib "../stat";
use lib "../../stat";
use Vector;
use lib "../Perl/modules";
use GenomicsUtils;

our (@ISA, @EXPORT);
@ISA = qw(Exporter);
@EXPORT = qw();

{
    sub aaCodons($) {
	my $codonTableR = shift;
	my %aaC = ();
	map { addToHashList(\%aaC, $codonTableR->{$_}, $_) } keys %$codonTableR;
	return %aaC;
    }

    sub sumAA($$$) {
	my ($inRef, $outRef, $aaCodonsRef) = @_;
	foreach my $aa (keys %$aaCodonsRef) {
	    $outRef->{$aa} = Vector::sum( map { $inRef->{$_} } @{$aaCodonsRef->{$aa}} );
	}
    }

    sub CBI($$) {
	my ($codonR, $aaCodonsR) = @_;
	my $cbi = 0;
	my $totalCodon = 0;
	foreach my $aa (keys %$aaCodonsR) {
	    if ($aa ne "*" && scalar @{$aaCodonsR->{$aa}} > 1) {
		my $sum = 0;
		map { $sum += $codonR->{$_}==0 ? .5 : $codonR->{$_} } @{ $aaCodonsR->{$aa} };
		$totalCodon += $sum;
		foreach my $codon (@{$aaCodonsR->{$aa}}) {
		    my $count = $codonR->{$codon}==0 ? .5 : $codonR->{$codon};
		    $cbi += $count*log($sum/$count);
		}
	    }
	}
	return($cbi/$totalCodon);
    }

    sub BestCodonCount($$) {
	my ($codonsR, $aaToCodonsR) = @_;
	my %bestCodonCount = ();
	foreach my $aa (keys %$aaToCodonsR) {
	    $bestCodonCount{$aa} = 1;
	    foreach (@{ $aaToCodonsR->{$aa} }) {
		$bestCodonCount{$aa} = $codonsR->{$_} if $codonsR->{$_} > $bestCodonCount{$aa};
	    }
	}
	return %bestCodonCount;
    }

    # codonR is the codon usage of the gene, and aaCodons is the list
    sub CAI($$$$) {
	my ($codonR, $aaCodonsR, $refCodonR, $aaRefBestR) = @_;
	my $cai = 0;
	my $tCodons = 0;
	foreach my $aa (keys %$aaCodonsR) {
	    if ($aa ne "*" && scalar @{$aaCodonsR->{$aa}} > 1) {
		foreach my $codon (@{ $aaCodonsR->{$aa} }) {
		    if ($codonR->{$codon} > 0) {
			my $refCount = $refCodonR->{$codon};
			$refCount = .5 if $refCount == 0;
			my $weight = $aaRefBestR->{$aa}/$refCount;
			$cai += $codonR->{$codon} * log($weight) if $weight > 1 && $weight != 0;
			$tCodons += $codonR->{$codon};
		    }
		}
	    }
	}
	return exp($cai/$tCodons);
    }
}
return 1;
