#!/usr/bin/perl -w

use strict;

open (FH,"<$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;#guarda el nom de cada matriu
my $posicio = 0; #posició dins del hash
my $compta = 0;#variable que conta la llargada del hash
my $i; #files del motiu
my $j; #columnes del motiu
my @k= ("A","C","G","T");#vecator que ens agafarà les variables que fan de claus a la matriu

while (<FH>){ 
    
#primera part: aqui contruim les matrius de ocurències
    
    
    if ($_ =~ m/FA/){
	$nom = $_;
	chomp ($nom);
	print "$nom\n";
    }

    if ($_=~ m/(\d+)\t(\d+)\s+(\d+)\s+(\d+)\s+(\d+)$/){ #aqui agafem per linies de l'arxiu de matrius per tal d'incorporar els valors en el nostre hash
	
	$motiu{"A"}[$posicio]= $2;
	$motiu{"C"}[$posicio]= $3;
	$motiu{"G"}[$posicio]= $4;
	$motiu{"T"}[$posicio]= $5;
	
	$posicio = $posicio +1;
	$compta = $compta +1;
    }
    
    if (($_ =~ m/(\/\/)\Z/) && ($compta > 1)){ #quan acabem de captar totes les variables de la matriu volem seguir fent tots els calculs preestablerts, sabem que les matrius acaben amb //.
	
	
	
#segona part: transformar la nostra matriu d'ocurrències en una matriu de pesos.
	
	open (PROMOTOR,"<$ARGV[1]");
	
	my $dna="";#variable que anira agafant cada linia del arxiu del promotor que incorporem.
	my @v; #vector que contindra la nostra seqüència del promotor.
	my $idna = 0; #per contar posicions de la seq promotora dins del vector v
	my $l;#longitud de la seqüència promotora incorporada
	my $A = 0;
	my $C = 0;
	my $G = 0;
	my $T = 0;
	
	
	while (<PROMOTOR>){
	    
	    chomp ($_);

	    if ($_ =~ m/>/){
	    }
	    else {

	    $dna = $dna.$_; #et guarda cada linia en la variable perque quedi una sola linia de lectura en el vector posterior
	}
	}
	close (PROMOTOR);

	@v = split (//,$dna);
	
	$l = scalar (@v);
	
	
	while ($idna <$l){ #ves llegint tota la seqüencia promotora
	    
	    if ( $v[$idna]eq "A"|| $v[$idna] eq 'a' ){
		
		$A = $A + 1;
	    }
	    if ( $v[$idna]eq "C"|| $v[$idna] eq 'c' ){
		
		$C = $C + 1;
	    }
	    if ( $v[$idna]eq "G" || $v[$idna] eq 'g'){
		
		$G = $G + 1;
	    }
	    if ( $v[$idna]eq "T"|| $v[$idna] eq 't' ){
		
		$T = $T + 1;
	    }
	    $idna = $idna +1;
	}
	
	
	my  $logA = (log($A / $l)/log(10));
	my  $logC = (log($C / $l)/log(10));
	my  $logG = (log($G / $l)/log(10));
	my  $logT = (log($T / $l)/log(10));

	#aqui hem calculat els logaritmes que t'indiquen la freqüència en que trobem cada nucleòtid dina la nostra seqüència promotora.
	
	
	my @log = ($logA,$logC,$logG,$logT); #col·loquem cada un dels logaritmes en un vector per tal de realitzar les operacions següents alhora de construir la matriu de pesos.
	my $x = scalar (@log);#longitud d'aquest nou vector
	
	
	
	$j = 0;
	
	
	while ($j < $compta){#construcció de la matriu de pesos
	    $i = 0;
	    $x = 0;
	    
	    my $suma = 0;
	    
	    while ($i <scalar (@k)){
		
		$suma = $suma + $motiu{$k[$i]}[$j];
		$i = $i + 1;
	    }
	    $i = 0;
	    
	    
	    while ($i < scalar (@k)){
		
		if ($motiu{$k[$i]}[$j] == 0){#si l'ocurrència del nucleòtid que observem és 0, asssignem directament un pes de -999.
		    
		    $motiu{$k[$i]}[$j] = -999
		    }
		
		else {
		    
		    $motiu{$k[$i]}[$j]= (log ($motiu{$k[$i]}[$j]/$suma)/log (10)) - $log[$x]; #valors de la matriu de pesos
	        }
		
		$i = $i + 1;
		$x = $x + 1;
	    }
	    
	    $j = $j + 1;
	}
	
#tercera part: aqui calculem l'score millor de la nostra seqüència de promotor i la seva posició en la seqüència
	
	
	my $score = 0; #el que tu calcules
	
	my $scoremax = -999999999;#el score màxim un cop calculem tots els possibles scores amb la matriu de pesos
	
	my $posiciomillor = 0; #posicio millor dina la seqüència
	
    	my $pmotiu; #comptar posicio de la sequencia nomes dins del motiu
	my $f = 0; #comptar posicio dins del hash
       	
	$idna = 0;#posició dins de tota la seqüència
	    
	    while ($idna < $l-($compta)+1) {
		
		$f = 0;
		$score = 0;
		$pmotiu = $idna;
		
		
		while ($f <  $compta) { 
		    
		    if ($v[$pmotiu] eq 'A' || $v[$pmotiu] eq 'a' ){
			
			$score = $score + ($motiu{"A"}[$f]);
		    }
		    
		    if ($v[$pmotiu] eq 'C' || $v[$pmotiu] eq 'c' ){
			
			$score = $score + ($motiu{"C"}[$f]);
		    } 
		    
		    if ($v[$pmotiu] eq 'G' || $v[$pmotiu] eq 'g' ){
			
			$score = $score + ($motiu{"G"}[$f]);
		    }
		    
		    if ($v[$pmotiu] eq 'T' || $v[$pmotiu] eq 't' ){

			$score = $score + ($motiu{"T"}[$f]);
		    }
		    
		    $pmotiu = $pmotiu + 1;
		    $f = $f + 1; 
		    
		}	
		
		
		if ($score > $scoremax) {#aqui ens anirà guardant els scores més grans que vagi trobant
		    $scoremax = $score;
		    $posiciomillor = $idna;
		}
		
		
		$idna = $idna + 1;
		
	    }
	
	
	print "El factor de transcripció $nom té la millor puntuació, que es $scoremax en la posicio $posiciomillor\n";
	
	
#quarta part: aqui contem els scores que ens dona la nostra seqüència un cop ha estat reordenada a l'atzar, i això ho mirem 100 vegades

	my $voltes = 0; #per repetir l'operació 100 cops.
	my $scorerandom = -999999; #millor score del factor de transcripció a la seqüència quan fem la barreja, o random
	my $c = 0; #vegades que el scorerandom és igual o superior al scoremax
	my $pvalue;
	
	while ($voltes < 100 ){
	    
	    my $n = scalar(@v);
	    my $i = $n - 1;
	    
	    while ($i >= 0) {
		
		my $j = int(rand($i+1));
		
		if ($i != $j) {
		    
		    my $tmp = $v[$i];
		    
		    $v[$i] = $v[$j];
		    $v[$j] = $tmp;
		}
		$i = $i - 1;
	    }
	    
	    $idna = 0;
	    $posiciomillor = 0;
	    
	    while ($idna < $l-($compta)+ 1) {
		
		$f = 0;
		$score = 0;
		$pmotiu = $idna;
		
		while ($f < $compta) { 
		    
		    if ($v[$pmotiu] eq 'A' || $v[$pmotiu] eq 'a' ){
			
			$score = $score + ($motiu{"A"}[$f]);
		    }
		    
		    if ($v[$pmotiu] eq 'C' || $v[$pmotiu] eq 'c' ){
			
			$score = $score + ($motiu{"C"}[$f]);
		    } 
		    
		    if ($v[$pmotiu] eq 'G' || $v[$pmotiu] eq 'g' ){
			
			$score = $score + ($motiu{"G"}[$f]);
		    }
		    
		    if ($v[$pmotiu] eq 'T' || $v[$pmotiu] eq 't' ){
			
			$score = $score + ($motiu{"T"}[$f]);
		    }
		    
		    $pmotiu = $pmotiu + 1;
		    $f = $f + 1; 
		    
		}	
	         
	    if ($score >= $scoremax){
		
		$c = $c + 1;
		$idna = $l;
		
	    }
	    	$idna = $idna + 1;
		
	    }

	    $voltes = $voltes + 1;
	}
		
		
	
	if ($c == 0) {
	    
	    print " el p-value del factor de trasncripció $nom és 0,0\n";
	}
	else {
	    
	    $pvalue = $c/100;
	    
	    print "el p-value del factor de transcripció $nom és igual a $pvalue \n";
	    print "\n";
	}
	
	#tornem a posar el hash a zero perquè agafi la segona matriu i ho torni a calcular tot tantes vagades com matrius hi hagi.

	%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]);
	
	$compta = 0;
	$posicio = 0;
	
    }
}

   

close (FH);
