#!/usr/bin/perl -w

use strict;

open (FITXERMEU,"<$ARGV[0]");


my %motiu = ("A" => [0,0,0,0,0,0,0,0,0],
	     "C" => [0,0,0,0,0,0,0,0,0],
	     "G" => [0,0,0,0,0,0,0,0,0],
	     "T" => [0,0,0,0,0,0,0,0,0]);
my $nom;

my $posicio = 0;
 


my $i;





while (<FITXERMEU>){


    if ($_ =~ m/FA\s+(.+)/){

	$nom = $1;

	print "\nNom: $nom\n\n"; #printar el nom del FT 

    }

    if ($_=~ m/\d+\t(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/){
 
	$motiu{"A"}[$posicio]= $1;

	$motiu{"C"}[$posicio]= $2;

	$motiu{"G"}[$posicio]= $3;

	$motiu{"T"}[$posicio]= $4;

	$posicio = $posicio +1;


    }


    
 
    if (($_ =~ m/(\/\/)/) && ($posicio > 1)) { 


   
	open(SEQUENCIA, "<$ARGV[1]");

	my $dna = " ";
	my @seq;
	my $longitud = 0;

	while (<SEQUENCIA>){
	    chomp ($_);
	    if ($_ =~ m/>/){
	    }
	    else {
		$dna = $dna.$_;
	    }
	}

	close (SEQUENCIA); #una sola variable amb totes les dades de la seq.



	@seq = split(//,$dna);

	$longitud = scalar(@seq); #longitud de la seq.






	my $pos = 0;

	my $A = 0;

	my $C = 0;

	my $G = 0;

	my $T = 0;


	while ($pos < $longitud) {

	    if ($seq[$pos] eq 'A' || $seq[$pos] eq 'a' ) {

		$A = $A + 1;

	    }

	    if ($seq[$pos] eq 'C' || $seq[$pos] eq 'c') {

		$C = $C + 1;

	    }

	    if ($seq[$pos] eq 'G' || $seq[$pos] eq 'g') {

		$G = $G + 1;

	    }

	    if ($seq[$pos] eq 'T' || $seq[$pos] eq 't') {

		$T = $T + 1;

	    }

	    $pos = $pos + 1;

	}


	my $probA;

	my $probC;

	my $probG;

	my $probT;


	if ($A == 0) {

	    $probA = -999;

	}

	else { $probA = log($A / $longitud) / log(10);

	   }

	if ($C == 0) {

	    $probC = -999;

	}

	else { $probC = log($C / $longitud) / log(10);

	   }

	if ($G == 0) {

	    $probG = -999;

	}

	else { $probG = log($G / $longitud) / log(10);

	   }

	if ($T == 0) {

	    $probT = -999;

	}

	else { $probT = log($T / $longitud) / log(10);

	   }

	

# Declarem una nova matriu on ens incorpori els logaritmes de probabilitat en cada posició. A més hi restarem els logaritmes de la seqüència promotora.

	my %logaritme = ("A" => [0,0,0,0,0,0,0,0,0],
			 "C" => [0,0,0,0,0,0,0,0,0],
			 "G" => [0,0,0,0,0,0,0,0,0],
			 "T" => [0,0,0,0,0,0,0,0,0]);

	$pos = 0;

	while ($pos < $posicio) { 

	    if ($motiu{"A"}[$pos] == 0) {

		$logaritme{"A"}[$pos] = -999;

	    }
	    
	    else {

		$logaritme{"A"}[$pos] = (log(($motiu{"A"}[$pos]) / ($motiu{"A"}[$pos] + $motiu{"C"}[$pos] +  $motiu{"G"}[$pos] + $motiu{"T"}[$pos])) / log(10)) - $probA;

	    }

	    if ($motiu{"C"}[$pos] == 0) {

		$logaritme{"C"}[$pos] = -999;

	    }

	    else {

		$logaritme{"C"}[$pos] = (log(($motiu{"C"}[$pos]) / ($motiu{"A"}[$pos] + $motiu{"C"}[$pos] +  $motiu{"G"}[$pos] + $motiu{"T"}[$pos])) / log(10)) - $probC;

	    }

	    if ($motiu{"G"}[$pos] == 0) {

		$logaritme{"G"}[$pos] = -999;

	    }
	    
	    else {

		$logaritme{"G"}[$pos] = (log(($motiu{"G"}[$pos]) / ($motiu{"A"}[$pos] + $motiu{"C"}[$pos] +  $motiu{"G"}[$pos] + $motiu{"T"}[$pos])) / log(10)) - $probG;

	    }

	    if ($motiu{"T"}[$pos] == 0) {

		$logaritme{"T"}[$pos] = -999;

	    }

	    else {

		$logaritme{"T"}[$pos] = (log(($motiu{"T"}[$pos]) / ($motiu{"A"}[$pos] + $motiu{"C"}[$pos] +  $motiu{"G"}[$pos] + $motiu{"T"}[$pos])) / log(10)) - $probT;

	    }

	    $pos = $pos + 1;

	}



#calculem la millor puntuació i guardem la posició que correspon a la millor puntuació


    

	my $magatzem = 0; #guarda la puntuació de cada motiu
	
	my $millor = -9999999999999 ; #guarda la millor puntuació
	
	my $motiu = $posicio; #representa la llargada de cada matriu de pesos
			
	my $xocolata = 0; #guarda la millor posició
	
	my $final; 
		
	my $p = 0; #comptar posició de la seqüència només dins del motiu

	my $l = 0; #comptar posició dins del motiu

	my $r; #anar comptant posisions al llarg de la seqüència
	
	$r = 0;

	while ($r < $longitud-($motiu)+1) {
	    
	 $l = 0;
	 $magatzem = 0;
	 $p=$r;

	    while ($l < $motiu) { 
		
		if ($seq[$p] eq 'G' || $seq[$p] eq 'g' ){
		    $magatzem = $logaritme{"G"}[$l] + $magatzem;
		}

		if ($seq[$p] eq 'A' || $seq[$p] eq 'a' ){
		    $magatzem = $logaritme{"A"}[$l] + $magatzem;
		}

		if ($seq[$p] eq 'C' || $seq[$p] eq 'c' ){
		    $magatzem = $logaritme{"C"}[$l] + $magatzem;
		}

		if ($seq[$p] eq 'T' || $seq[$p] eq 't' ){
		    $magatzem = $logaritme{"T"}[$l] + $magatzem;
		}
		$p = $p + 1;
		$l = $l + 1; 
			
	    }	
	    if ($magatzem > $millor) {
		$millor = $magatzem;
		$xocolata = $r;
	
		$final=$xocolata+$motiu-1;
	    }

	    $r = $r +1;
	      
	}

	

	print "La millor puntuació s'obté quan aquest FT s'uneix a la seqüència promotora en els nucleòtids $xocolata - $final\nLa puntuació és: $millor\n";
	

#calcular pvalue


	my $voltes = 0; #repetir l'operació 100 cops.
	
	my $c = 0; #vegades que el millorrandom és igual o superior al millor

	my $pvalue;
	



	while ($voltes < 100 ){
   
	    my $i = $longitud - 1;
	    
	    while ($i >= 0) {
		
		my $j = int(rand($i+1));
		
		if ($i != $j) {
		    
		    my $tmp = $seq[$i];
		    
		    $seq[$i] = $seq[$j];
		    $seq[$j] = $tmp;
		}
		$i = $i - 1;
	    }
	    
	my $millorrandom = -999999; #millor puntuació del factor de transcripció en la seqüència random
   

	$r = 0;


	while ($r < $longitud-($motiu)+1) {
	    
	 $l = 0;
	 $magatzem = 0;
	 $p=$r;

	    while ($l < $motiu) { 
		
		if ($seq[$p] eq 'G' || $seq[$p] eq 'g' ){
		    $magatzem = $logaritme{"G"}[$l] + $magatzem;
		}

		if ($seq[$p] eq 'A' || $seq[$p] eq 'a' ){
		    $magatzem = $logaritme{"A"}[$l] + $magatzem;
		}

		if ($seq[$p] eq 'C' || $seq[$p] eq 'c' ){
		    $magatzem = $logaritme{"C"}[$l] + $magatzem;
		}

		if ($seq[$p] eq 'T' || $seq[$p] eq 't' ){
		    $magatzem = $logaritme{"T"}[$l] + $magatzem;
		}
		$p = $p + 1;
		$l = $l + 1; 
			
	    }	
	    if ($magatzem > $millorrandom) {
		$millorrandom = $magatzem;
	
	    }
	    
	    $r = $r +1;
	      
     }
  
	
	    
	    if ($millorrandom >= $millor){
		
		$c = $c + 1;
		
	    }
	    
	    $voltes = $voltes + 1;
	}
	
	
	    
	$pvalue = $c/100;
	    
	print "El p-value és igual o mes petit que $pvalue \n";
	print "\n";





	%motiu = ("A" => [0,0,0,0,0,0,0,0,0],
	  	  "C" => [0,0,0,0,0,0,0,0,0],
	    	  "G" => [0,0,0,0,0,0,0,0,0],
	     	  "T" => [0,0,0,0,0,0,0,0,0]);

	
	$posicio = 0;
	
    } #tanquem el tercer if per indicar que ja tenim la primera matriu de pesos feta i reinicialitzem el hash inicial per llegir el segon factor de transcripció.


}

close (FITXERMEU);
