#!/usr/bin/perl -w

use strict;


# si el nombre d'arguments al programa es mes petit que 2
# no l'hem pogut cridar amb els noms de fitxers que hem d'obrir
# despres

if (scalar(@ARGV) < 2) {
  print "dnaloglkh.pl fitxercodons fitxerseq\n";
  exit(1);
}


# assignem els noms del vector d'arguments (parametres) a dues
# variables amb noms mes significatius

my $fitxercod = $ARGV[0];
my $fitxerseq = $ARGV[1];


# obrim el fitxer especificat al primer argument passat al
# programa des del shell

if (!open(PROPCODONS,"< $fitxercod")) {
  print "dnaloglkh.pl: impossible obrir $fitxercod\n";
  exit(1);
}


# llegim les proporcions de codons del fitxer i les enregistrem
# al hash %pcodons

my %pcodons;

while (<PROPCODONS>) {

  # a la variable anonima tenim la linia llegida i l'encaixem amb
  # la seguent expressio regular

  m/(\w+) (\d+\.\d+)/;

  # l'expressio regular anterior encaixa amb dues paraules la primera
  # de les quals es alfabetica i la segona un numero amb decimals.
  # les dues paraules estan agrupades amb parentesis lo qual significa
  # que ens podem referir a la primera com $1 i a la segona com $2

 $pcodons{$1} = $2 ;

}

close(PROPCODONS);


# obrim el fitxer especificat al segon argument passat al
# programa des del shell

if (!open(SEQUENCIA,"< $fitxerseq")) {
  print "dnaloglkh.pl: impossible obrir $fitxerseq\n";
  exit(1);
}


# el fitxer amb la sequencia esta en format FASTA. aixo implica que te
# una primera linia amb el simbol '>' seguit d'un identificador, i despres
# segueix un seguit de linies amb differents trossos de la sequencia
# canviarem el marcador de registre de '\n' a '>' per tal de llegir la
# sequencia tota d'un cop.

my $primeralinia = <SEQUENCIA>;

$/ = '>';

my $x = <SEQUENCIA>;

close(SEQUENCIA);


# voldriem tenir la sequencia en una sola variable sense la primera linia
# '>identificador'. per tal de fer aixo ficarem el contingut llegit
# (el registre sencer) en un vector on cada posicio contindra una linia

my @v = split(/\n/,$x);


# en la variable $seq ficarem la sequencia sencera. aixo ho farem enganxant
# els trossos de sequencia que tenim a cada posicio del vector @v

my $seq="";
my $i=0;

while ($i < scalar(@v)) {
  $seq = $seq.$v[$i];
  $i = $i + 1;
}


# utilitzarem el mateix vector @v per tal d'emmagatzemar els codons de la
# sequencia. per tal de trocejar la sequencia en codons, l'encaixarem amb
# l'expressio regular que representa 3 simbols qualsevols utilitzant el
# modificador 'g' per tal de que faci totes les correspondencies possibles
# al llarg de la sequencia i aixi obtenir el vector de codons

$seq = substr($seq,$ARGV[2]);

@v = ($seq =~ m/.../g);


# ara solament ens queda ja calcular la rao de versemblanc,a en escala
# logaritmica. la rao de versemblanc,a es la probabilitat conjunta dels
# codons segons la taula donada, respecte (dividit per) la probabilitat
# conjunta uniforme dels codons (la probabilitat uniforme d'un unic codo
# es 1/64). aixo ho calcularem en escala logaritmica per tal de poder
# manipular els nombres molt petits que es poden produir al multiplicar
# moltes xifres que tenen valors entre 0 i 1. el logaritme del producte
# es la suma de logaritmes i el logaritme de la divisio es la resta de
# logaritmes.

my $r = 0.0;

$i=0;
while ($i < scalar(@v)) {
    if ($pcodons{$v[$i]} == 0){
	$pcodons{$v[$i]} = 0.00000000000000000000000000000000001;
    }
  $r= $r + log($pcodons{$v[$i]}) - log(1/64)  ;

  $i = $i + 1;

}

printf "%.2f\n",$r;
