#!/usr/bin/perl -w

use strict;
use Math::Libm ':all';

#############################################################################
############################# VARIABLES #####################################
#############################################################################

my $fitxerfasta;                         #nom de l'arxiu que conte les sequencies en format fasta
my $lmotiu;                              #la longitud del motiu que lusuari introdueix
my @input;                               #per treure 2 vectors referenciats de la funcio llegirseq
my @ident;                               #emmagatzema les identificacions de les sequencies
my @seq;                                 #emmagatzema els aminoacids de cada sequencia
my %mat_pesos;                           #freqs absoluts de cada aminoacid en el motiu-bckg
my %mat_relatives;                       #freqs relatives de cada aminoacid en el motiu-bckg
my $counts=0;                            #comptar les voltes que fa la iteracio
my %aastotals;                           #freqs absolutes x aminacid, sns importar la posicio
my @pos_candidats;                       # millors posicions d'inici del motiu per cada seq
my @pesos;                               #pesos tots candidats de cada sequencia
my @out;                                 #per treure una hash referenciat i un valor de la funcio M_step
my $loglkh= -99999999;                   #la versemblan.bŽça del model que estem calculant
my $current_loglkh= -999999999999;       #per comparar la versemb del model anterior-actual
my @aas = ( "D", "E", "A", "R", "N", "C", "F", "Q", "G", "H", "I", "L", "K", "M", "P", "S", "T", "W", "Y", "V" );
my $dif=1;
#############################################################################
########################### COS DEL PROGRAMA ################################
#############################################################################

if (@ARGV < 2) {

    print "EM.pl fitxerfasta longitudmotiu\n";  #informa que falten arguments si no hi hem posat els 2
    exit;
}

$fitxerfasta = $ARGV[0];                      #argument amb el fitxer fasta per processar
$lmotiu = $ARGV[1];                           #argument amb la longitud del motiu q volem buscar
                                             
@input = llegirseq ($fitxerfasta);             #llegeix les sequencies i retorna:

@ident = @{$input[0]};                          #vector amb les identitats

@seq = @{$input[1]};                            #vector amb cada sequencia en una sola linia i en una posicio
                                                #del vector

%aastotals= contaraas (\@seq);                 #comptatge de les ocurrencies de cada aminoacid en totes les sequencies

%mat_pesos = sliding (\@seq);                   # genera amb el metode indicat a la pagina, la millor matriu inicial 
 
%mat_relatives = frec_relatives (\%mat_pesos);           #fa les frequencies relatives de la millor 
                                                         # matriu que ha trovat

$counts =0;

$loglkh = -99999999;                                    # reinicia els valors de loglkh
$current_loglkh = -99999999999999;                      #

while ($dif>(1/100) && $counts < 1000) {                #  mentre el loglkh millori mes de 0,01 en valor absolut 
                                                                #  i no hagi fet mes de mil iteracions

    $current_loglkh = $loglkh;                                  
 
    @pesos = E_step (\@seq, \%mat_relatives);                    # fes E i M step, i avalua el model

    @out = M_step (\@seq, \@pesos,\%aastotals);

    %mat_relatives=%{$out[0]};

    $loglkh= $out[1];

    $counts = $counts +1;
    $dif = diferencia ($loglkh,$current_loglkh);               #calcul de la diferencia en valor absolut

}

@pesos = E_step ( \@seq, \%mat_relatives);                     #la matriu de pesos bona la pases per totes les seq

@pos_candidats = posiciocandidats (\@pesos);                   #et quedes amb la posicio de cada seq que dona mes probabilitat

print "\n\nLa versemblansa del model es: $loglkh\n\n\n";       # Comensa l'output

output (\@seq, \@ident, \@pos_candidats, \%mat_relatives);

##############################################################################
#################################### FUNCIONS ################################
##############################################################################



sub llegirseq {                               #  posar les ident de la sequencia en el vector @ident
                                              #  i les sequencies concatenades en el vector @seq
    my $seq1 = "";
    my @ident;
    my @seq;
    my $i=0;
    my @input;

    if (!open(SEQUENCIES,"<$_[0]")) {         #si no es pot obrir mostra error i tanca el programa

	print "EM.pl impossible obrir $fitxerfasta\n";
	exit;

    }


    while () {

	if (m/(>.+)/) {

	    $ident[$i]=$1;
	    $seq1="";
	    $i=$i+1;

	} else {

	    m/([DEARNCFQGHILKMPSTWYVdearncfqghilkmpstwyv]+)/;
	    $seq1=$seq1.$1;

	}
	$seq[@ident-1]=uc($seq1);
    }

    $input[0] = \@ident;                      #referenciem els vectors per posarlos al vector input
    $input[1] = \@seq;
    return @input;                            #vector amb les ident de les seqs i les seqs concatenades
}
                                              
sub matriuinicial {                           
                                              
    my @seq =@{$_[0]};
    my @posinici =@{$_[1]};
    my $i=0;
    my %mat_pesos = iniciarmatriu ($lmotiu);  #  hash per anotar les freqs absolutes dels nucleotids
                                              #  en el background i en les pos del motiu

    while ($i<@seq) {

	my $inici;
	my $j=0;
	my @linia = ( $seq[$i] =~ m/[ACTG]/g );
	$inici = $posinici [$i];

	while ($j<@linia) {

	    if ($j == $inici) {               #  si estem dins el motiu, suma 1 al
                                              #  aminoacid que correspongui en cada posicio
		my $k = 0;

		while ( $k < $lmotiu ) {

		    $mat_pesos{$linia[$j+$k]}[$k+1] = $mat_pesos{$linia[$j+$k]}[$k+1] + 1;
		    $k = $k + 1;

		}

		$j = $j + $lmotiu;

	    }                                 #   quan no estem dins el motiu, suma 1 al aminoacid que correspongui
                                              #   en la posicio 0 del hash (background)
	    $mat_pesos{$linia[$j]}[0] = $mat_pesos{$linia[$j]}[0] + 1;
	    $j = $j +1;
	}
	$i = $i + 1;
     }

    return %mat_pesos;                        #  hash amb les freqs absolutes de cada
}                                             #  nucleotid en el background
                                              #  i suposant el lloc d'inici del motiu

sub contaraas {                        #  calculem els aminoacids de cada que hi ha
                                              #  en totes les sequencies
    my @seq = @{$_[0]};
    my $i =0;
    my %aastotals;           #comencem inicialitzant tots els valors a 0
    
    while ($i<@aas) {

	$aastotals{$aas[$i]} = 0;
	$i = $i +1;

    }

    $i=0;

    while ($i<@seq) {                         #  per cada nucleotid que trobem de cada,
                                              #  sumarem 1 al valor total d'aquell aminoacid
	my $j=0;
	my @linia = ($seq[$i] =~ m/\w/g);

	while ($j < @linia) {

	    $aastotals{$linia[$j]}= $aastotals{$linia[$j]} + 1;
	    $j =$j+1;
	}
	$i = $i +1;
    }
    return %aastotals;                       #hash amb els valors absoluts de cada aminoacid
}


sub iniciarmatriu {                           #  inicialitzem la matriu amb el background i les
                                              #  posicions del motiu amb 0 a totes les posicions
    my $lmotiu = $_[0];
    my $i=0;
    my %mat_pesos;

    while ($i<$lmotiu+1) {
	
	my $j=0;
	
	while ($j<@aas) {

	    $mat_pesos{$aas[$j]}[$i]=0;
	    $j=$j+1;

	}

	$i = $i +1;

    }
    return %mat_pesos;
}

sub diferencia {

    my $loglkh = $_[0];
    my $current_loglkh = $_[1];
    my $dif;

    if ($loglkh=~ m/\-(\d+.\d+)/) {
	$loglkh = $1;
    }

    if ($current_loglkh=~ m/\-(\d+.\d+)/) {
	$current_loglkh = $1;
    }

    if ($loglkh<$current_loglkh) {

	$dif = ($current_loglkh-$loglkh);

    } else {

	$dif = ($loglkh-$current_loglkh);
    }
    return $dif;
}
                                              #  calcul de les freqs relatives dels aminoacids
sub frec_relatives {                          #  en el background i les posicions del motiu

    my %mat_pesos = %{$_[0]};                 #hash amb les frequencies absolutes
    my %mat_relatives;                        #hash que anem a calcular amb les freqs relatives
    my $i=0;

    while ($i<$lmotiu+1) {                    #  recorrem totes les posicions del %mat_pesos i
                                              #  dividim cada valor per la suma de valors de la columna
	my $j =0;
	my $suma = 1;

	while ($j<@aas) {

	    $suma = $suma + $mat_pesos{$aas[$j]}[$i];
	    $j=$j+1;
	}

	$j=0;

	while ($j<@aas) {

	    $mat_relatives{$aas[$j]}[$i] = ($mat_pesos{$aas[$j]}[$i]+(1/20))/$suma;
	    $j = $j +1;
	}

	$i = $i +1;
   }
    return %mat_relatives;                    #hash de les frequencies relatives
}
                                              #  pas en que puntuem cada candidat en cada
sub E_step {                                  #  sequencia segons el model que estem usant
                                              #  el model el representa el %mat_relatives
    my @seq = @{$_[0]};
    my %mat_relatives = %{$_[1]};
    my @pesos;
    my $i=0;


    while ($i < @seq) {                       # per cada sequencia:

	my @linia = ($seq[$i] =~ m/\w/g); #  cada sequencia la fragmentem en un @linia
	my $inici = 0;                        #  que conte un nucleotid en cada posicio
	my $suma = 0;
        my @sum;

	while ($inici <= @linia-$lmotiu) {   #per cada possible lloc d'inici del motiu:
                                             
	    $pesos[$i][$inici] = 0;

	    my $j=0;

	    while ($j < @linia) {             #per cada aminoacid:


		if ($j == $inici) {           #  quan entrem dins el motiu, puntuem el candidat
                                              #  segons la posicio q ocupen els aminoacids amb els valors de %mat_relatives
		    my $k=0;

		    while ($k < $lmotiu) {

			$pesos[$i][$inici] = $pesos[$i][$inici]+log($mat_relatives{$linia[$j+$k]}[$k+1]);
			$k = $k+1;
		    }

		    $j =$j + $lmotiu;
		}

		if ($j != @linia) {           #  fora del motiu puntuem el candidat segons
                                              #  la freq del background, sns importar la posicio
		    
		    
		    $pesos[$i][$inici] = $pesos[$i][$inici]+log($mat_relatives{$linia[$j]}[0]);
		    $j = $j+1;
		}
	    }
	    
	    $sum[$inici]=$pesos[$i][$inici];
	    
	    $inici = $inici +1;
	}

	$suma=logsum(@sum);
	$inici=0;                             #  reinicialitzem altre cop per normalitzar
                                              #  els pesos de cada candidat en cada sequencia
	while ($inici<@{$pesos[$i]}) {

	    $pesos[$i][$inici] = exp($pesos[$i][$inici] - $suma);
	    $inici = $inici + 1;
	}

	$i=$i+1;
    }

    return @pesos;                            #  retorna una matriu amb les probabilitats 
                                              #  per cada posicio que, donada una  sequencia, 
}                                             #  el motiu comenci en aquella posicio


sub M_step {                                  #  updating de la matriu bckg-motiu per millorar
                                              #  el model segons les puntuacions dels candidats
    my @seq= @{$_[0]};
    my @pesos = @{$_[1]};
    my %aastotals = %{$_[2]};
    my $i=0;
    my %mat_relatives;
    my %mat_pesos = iniciarmatriu ($lmotiu);  #reinicialitzem el hash amb 0 a totes posicions
    my %suma;               
    my %sumalog;
		    
    while ($i<@aas) {

	$sumalog{$aas[$i]}=0;
	$suma{$aas[$i]}=0;
	$i=$i+1;
    }

    $i =0;

    my $versem;                               #expressa la versemblansa del model


    while ($i<@seq) {                         #per cada sequencia:

	my $j=0;

	while ($j<@{$pesos[$i]}){           #per cada lloc d'inici possible:

	    my @motiu = (substr($seq[$i],$j,$lmotiu) =~ m/\w/g); #agafem el motiu en un vector
	    my $k=0;
 
	    while ($k<@motiu){               #recorrem el motiu i  suma el pes d'aquell candidat al aminoacid i la posicio que troba

		$mat_pesos{$motiu[$k]}[$k+1] = $mat_pesos{$motiu[$k]}[$k+1] + $pesos[$i][$j]; 
		$k = $k +1;

	    }
	    $j = $j +1;
	}

	$i=$i+1;
    }

    $i = 1;

    while ($i<$lmotiu+1){
	my $j =0;

	while ($j<@aas) {

	    $suma{$aas[$j]} = $suma{$aas[$j]} + $mat_pesos{$aas[$j]}[$i];     # suma aminoacids de cada dins del motiu
                                                                              #  per fer el bakground restant de total
	    $j=$j+1;
	}  
	                                            
	$i = $i +1;
    }

    $i =0;

    while ($i<@aas) {

	$mat_pesos{$aas[$i]}[0] = $aastotals{$aas[$i]}-$suma{$aas[$i]};      # Fa el background
	$i=$i+1;
    }


    %mat_relatives = frec_relatives (\%mat_pesos);   #  recalculem les freqs relatives tenint
                                                     #  en compte les freq absolutes actualitzades
    $i=0;

    while ($i < $lmotiu+1) {                         #tornem a recorrer la matriu per calcular la versemblansa
 
	my $j =0;

	while ($j<@aas) {

	    $sumalog{$aas[$j]} = $sumalog{$aas[$j]} + ($mat_pesos{$aas[$j]}[$i] * log($mat_relatives{A}[$i]));
	    $j = $j+1;
	}

	$i = $i +1;
    }

    $i =0;
    $versem =0;

    while ($i<@aas) {

	$versem = $versem + $sumalog{$aas[$i]};
	$i = $i +1;

    }

    my @out = ( \%mat_relatives, $versem );         #retorna la matriu ponderada i el loglkh de la matriu

    return @out;

}


sub logsum {
my $len = scalar(@_);
my $x = 0;

if ($len > 1) {

$x = $_[0];
for (my $i = 1; $i < $len; $i++) {
$x = $x + log1p(exp($_[$i]-$x));
}

}

return $x;
}


sub sliding {                                                    

   
    my @seq = @{$_[0]};
    my $i =0;
    my $rand = int(rand(@seq));    #per una sequencia escollida al atzar 
    my @punts;   
    my %mat_pesos;
    my $suma2=0; 
    my @posicions;
   
  

    while ($i $suma2) {     #si la suma es mes gran que la de l'anterior finestra ens quedem les posicions d'inici

	    $suma2 = $suma ;
	    @posicions = @{$punts[0]};
	}

	$i = $i+5;
    }

    
    %mat_pesos = matriuinicial (\@seq,\@posicions);  #fem la matriu inicial amb aquestes posicions inicials calculades
    return %mat_pesos;

}

sub suma {

    my $i=0;
    my $suma =0;
    while ($i<@_) {
	$suma = $suma + $_[$i];
	$i =$i+1;
    }
    return $suma;
}

sub max {
  my $i=0;
  my $max = -99999;
  my $posmax = 0;

  while ($i < scalar(@_)) {

    if ($_[$i] >= $max) {
      $max = $_[$i];
      $posmax = $i;
    }

    $i = $i + 1;

  }

  return ($posmax,$max);

}





sub posiciocandidats {                        #  un cop acceptat el millor model, selecciona el millor candidat
                                              #  de cada sequencia, guardant la posicio que ocupa i la probabilitat que te
    my @pesos = @{$_[0]};
    my $i=0;
    my @pos_candidats;

    while ($i<@pesos) {                       #per cada seq:

	my $j=0;
	my $max = -999999999;

	while ($j<@{$pesos[$i]}) {           #per cada lloc d'inici:

	    if ($pesos[$i][$j] > $max ) {     #si el valor que mirem es millor que el que tenim ens el quedem, i si no el descartem

		$max = $pesos[$i][$j];
		$pos_candidats[$i][0] = $j;
		$pos_candidats[$i][1] = $max;

	    }
	    $j= $j+1;
	}
	$i = $i +1;
    }
    return @pos_candidats;                    #vector que conte el millor candidat de cada sequencia
}


sub output {                                  #disseny de l'output, que mostrara el programa per la pantalla

    my @seq = @{$_[0]};
    my @ident = @{$_[1]};
    my @pos_candidats = @{$_[2]};
    my %mat_relatives = %{$_[3]};
    my $i=0;
    my %mat_loglkh;
    my %mat_info;
    my @bits_info;
 
    print "Identificador \t Motiu \t   Posicio  \t  Pes\n";

    while ($i< @ident) {

	my $motiu = substr($seq[$i],$pos_candidats[$i][0],$lmotiu);
	print "$ident[$i]  \t$motiu  \t",($pos_candidats[$i][0]+1),"  \t";
	printf "%.4f\n",$pos_candidats[$i][1]; 
	$i = $i +1;
    }
    $i = 1;

    while ($i < $lmotiu +1) { #calcul de la matriu de pesos amb el logaritme de la rao de versemblanŽça
                              # i dels bits d'informacio de cada posicio
	my $j =0;
	$bits_info[$i]=0;
	while ($j<@aas) {

	    $mat_loglkh{$aas[$j]}[$i] =log ($mat_relatives{$aas[$j]}[$i]) - log ($mat_relatives{$aas[$j]}[0]);


	    $mat_info{$aas[$j]}[$i] = ($mat_relatives{$aas[$j]}[$i]*log($mat_relatives{$aas[$j]}[$i])/log(2));

	
	    $bits_info[$i]= $bits_info[$i] + $mat_info{$aas[$j]}[$i];

	    $j =$j+1;
	}
	$bits_info[$i]= $bits_info[$i] +(log(20)/log(2));
	$i =$i+1;
    }
    print "\n\n";

    print "Matriu de frequencies relatives del motiu\n\n";
    printarhash (\%mat_relatives);

    $i=1;    
    
    print "Bits Info\nB\t";
    
    while ($i<@bits_info){

	printf "%.2f\t",$bits_info[$i];
	$i=$i+1;
    }
    print "\n\n";
    print "Matriu de pesos del motiu\n\n";
    printarhash (\%mat_loglkh);

   
}

sub printarhash {

    my %hash=%{$_[0]};
    my $i =1;
  
    print "Posicio\t";

    while ($i<$lmotiu+1) {

	print "pos$i\t";
	$i = $i+1;
    }

    $i=0;
    print "\n";

    while ($i<@aas) {
	my $j=1;
	print $aas[$i],"\t";
	while ($j<$lmotiu+1){
	    printf "%.2f\t",$hash{$aas[$i]}[$j];
	    $j=$j+1;
	}
	print "\n";
	$i =$i+1;
    }
    print "\n\n";
}





##############################################################################################
##############################################################################################
########################                                           ###########################
########################  implemented by: EVA MERCE & SERGI REGOT  ###########################
########################                                           ###########################
##############################################################################################
##############################################################################################