#!/usr/bin/perl -w

use strict;

my $nom; #serà perquè s'imprimeixi el nom de cada factor de transcripció
my %hash; #serà on guardarem la matriu del punt 1 de l'exercici, convertida en un hash

my @v2=(); #són les quatre els vectors del hash, A, C, T, G
my @v3=();
my @v4=();
my @v5=(); 

my $num_files = 0; #és per a anar comptant el número de files de la matriu inicial
my $p = 0;
my $x = 0; #servirà per evitar que llegeixi les dues // i pari el programa

my $ap = 0; #a fa referència al nucleòtid adenina, i p a %. Serveix per calcular el % de A a la nostra seqüència promotora. Punt2
my $cp = 0; #Serveix per calcular el % de C a la nostra seqüència promotora. Punt2
my $tp = 0; #Serveix per calcular el % de T a la nostra seqüència promotora. Punt2
my $gp = 0; #Serveix per calcular el % de G a la nostra seqüència promotora. Punt2

my %result = ("A" => [0, 0, 0, 0, 0, 0, 0],
	      "C" => [0, 0, 0, 0, 0, 0, 0],
	      "G" => [0, 0, 0, 0, 0, 0, 0],
	      "T" => [0, 0, 0, 0, 0, 0, 0]); #serà el hash on emmagatzarem els resultats de les freqüències



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

while (<MATRIUS>) {

    if (m/(FA)/){ 
	$x = 1;
	$nom = $_;
    }

    if (m/(\d+)[\t\s]+(\d+)[\t\s]+(\d+)[\t\s]+(\d+)[\t\s]+(\d+)/){

	$v2[$p] = $2;
	$v3[$p] = $3;
	$v4[$p] = $4;
	$v5[$p] = $5;

	$p = $p + 1;
	$num_files = $num_files + 1;

    }
   
    if (m/(\/\/\Z)/ && $x==1){

	%hash = ("A" => \@v2,
		 "C" => \@v3,
		 "G" => \@v4,
		 "T" => \@v5);

	$p = 0;
	

my $suma;
my $col = 0; #serveix per anar contant columnes
	


#--------------------------------------------------càlcul de les freqüències

	while ($col < scalar(@{$hash{"C"}})){ 
	    
	    $suma = 0;
	    
	    $suma = $hash{"A"}[$col] + $hash{"C"}[$col] + $hash{"G"}[$col] + $hash{"T"}[$col];
	    
	    
	    $result{"A"}[$col] = $hash{"A"}[$col] / $suma;
	    
	    $result{"C"}[$col] = $hash{"C"}[$col] / $suma;

	    $result{"G"}[$col] = $hash{"G"}[$col] / $suma;

	    $result{"T"}[$col] = $hash{"T"}[$col] / $suma;


	    $col = $col + 1;

	}

#-----------------------------------------------------------Comença el punt 2--------------------------------------------------------------

my $z = 0; #serveix per anar-nos desplaçant a la seqüència, per anar contant
my $a = 0; #serveix per contar quantes adenines hi ha a la seqüència
my $c = 0; #serveix per contar quantes citosines hi ha a la seqüència
my $t = 0; #serveix per contar quantes timidines hi ha a la seqüència
my $g = 0; #serveix per contar quantes guanines hi ha a la seqüència
	

	open (SEQ, "< $ARGV[1]"); #obrim el fitxer on hi ha la seqüència promotora del gen per calcular el % de cada nucleòtid
	my $log;
	my@v;
	my $SEQ;
	while(<SEQ>){
	    if(/>/){}
	    else{
		chomp($_);
		$SEQ= $SEQ . $_;	
	    }
	}
	@v= split(//,$SEQ);
	$log = scalar(@v);
	close (SEQ);



	$z = 0;
	    while ($z < $log){

		if ($v[$z] eq "g" || $v[$z] eq "G") {
		    
		    $g = $g + 1;
		}
		
		if ($v[$z] eq "a" || $v[$z] eq "A"){
		    
		    $a = $a + 1;
		}
		
		if ($v[$z] eq "c" || $v[$z] eq "C"){
		    
		    $c = $c + 1;
		}
		
		if ($v[$z] eq "t" || $v[$z] eq "T"){
		    
		    $t = $t + 1;
		}
		
		$z = $z + 1
		} 
	    
	    $ap = $a/$log;
	    $cp = $c/$log;
	    $tp = $t/$log;
	    $gp = $g/$log;
	    
	

#-------------------------un cop hem calculat el % de cada nucleòtid a la seqüència promotora apliquem l'equació del logaritme 
    
	my %pesos; #serà on emmagatzarem els resultats de la matriu de pesos
	my $f = 0; #serveix per anar contant files
	$col = 0;

   
	while ($col < scalar(@{$hash{"A"}})){

	    if ($result{"A"}[$col] == 0){

		$pesos{"A"}[$col] = -999;
		}
	    
	    else {

		$pesos{"A"}[$col] = log(($result{"A"}[$col])/(log(10))) - log(($ap)/(log(10)))
		
	    }

	    $col = $col + 1;
	}
	
	$col = 0;


	while ($col < scalar(@{$hash{"C"}})){
	
	    if ($result{"C"}[$col] == 0){

		$pesos{"C"}[$col] = -999;
		} 
	    
	    else {

		$pesos{"C"}[$col] = log(($result{"C"}[$col])/(log(10))) - log(($cp)/(log(10)));
	    }

	    $col = $col + 1;
	}

	$col = 0;


	while ($col < scalar(@{$hash{"G"}})){

	    if ($result{"G"}[$col] == 0){

		$pesos{"G"}[$col] = -999; 
		}
	    
	    else {

		$pesos{"G"}[$col] = log(($result{"G"}[$col])/(log(10))) - log(($gp)/(log(10)));
		}

	    $col = $col + 1;
	}

	$col = 0;


	while ($col < scalar(@{$hash{"T"}})){

	    if ($result{"T"}[$col] == 0){

		$pesos{"T"}[$col] = -999;
		}
	    
	    else {

		$pesos{"T"}[$col] = log(($result{"T"}[$col])/(log(10))) - log(($tp)/(log(10)));
		}

	    $col = $col + 1;
	}

 
    
#--------------------------------------------- s'acaba el punt 2 i comença el punt 3-------------------------------------------------------

	open (SEQ, "< $ARGV[1]");
	
	my $passet = 0; #anirà comptant les posicions dins de la finestra, dins el factor de transcripció
	my $pass = 0; #anirà comptant la primera posició de la finestra i la farà avançar a la seqüència
	my $puntuacio = 0; #és l'score que s'obté a la nostra seqüència
	my $long = 0; #és la longitud de la seqüència
	my $max = 0; 
	my $pos = 0;
	my $millor = -99999999999;
	my $posicio = 0;
       	my @seq;
	my @lin;
	my $LIN = 0;
	my $final;

	while (<SEQ>){
	    if(/>/){}
	    else{
		$LIN = $_;
		
		chomp ($LIN); 
		
		@lin = split (//,$LIN);
		
		push(@seq,@lin);
	    }
 	}

	$long = scalar(@seq);

	my  $finestra = scalar(@{$pesos{"A"}});
	

	while($pass < $long - $finestra){
	   
	    $passet = 0; 
	    $puntuacio = 0;	
	    
	    while ($passet < $finestra){

		if ($seq[$pass + $passet] eq "a" || $seq[$passet + $pass] eq "A"){
		    
		    $puntuacio = $puntuacio + $pesos{"A"}[$passet];
		}
		    
		if ($seq[$passet + $pass] eq "c" || $seq[$passet + $pass] eq "C"){
		    
		   $puntuacio = $puntuacio + $pesos{"C"}[$passet];
	  	}
		
		if ($seq[$passet + $pass] eq "g" || $seq[$passet + $pass] eq "G"){
		    
		   $puntuacio = $puntuacio + $pesos{"G"}[$passet];
		}
		
		if ($seq[$passet + $pass] eq "t" || $seq[$passet +$pass] eq "T"){
		    
		    $puntuacio = $puntuacio + $pesos{"T"}[$passet];
		}
		
		$passet = $passet + 1;
	    }

	        
	    if ($puntuacio > $millor){

		$millor = $puntuacio;
		$posicio = $pass;
		$final = $posicio+ $finestra-1;

	    }

	    $pass = $pass +1;
	
}

print "El factor de transcripció $nom s'uneix a la seqüència donada. La seva puntuació màxima és $millor i se situa entre els nucleòtids $posicio i $final.\n";

#------------------------------------------------s'acaba el punt 3 i comença el punt 4-----------------------------------------------------
 

my $h = 0;   #numero de vegades que mirem la seq random, que al final serà 100
my $vegades = 0; #numero de vegades que la puntuació a la  seqüència random serà superior a $millor del punt 3  	    
my $b;     #es situa a cada posició de la seq random des de la max longitud i anar disminuïnt	    
my $y;
my $millorandom;
my $tmp;
@v = @seq;
my $w = 0;
my $pvalue;
my $puntuacio2;
my $r = 0;

	while ($h < 100) { 
	    
	    $b = $long -1;
	    while ($b >= 0){
		$y = int(rand($b+1));
		if ($b != $w){
		    $tmp = $v[$b];
		    $v[$b] = $v[$y];
		    $v[$y] = $tmp;
		}
		
		$b = $b - 1;
	    }
	   
	    $pass = 0;
	    $vegades=0;

    while($pass < $long-$finestra){
		$passet = 0; 
		$puntuacio2 = 0; 
		
		while ($passet < $finestra){
		  
		    if ($v[$passet + $pass] eq "a" || $v[$passet + $pass] eq "A"){			     
			$puntuacio2 = $puntuacio2 + $pesos{"A"}[$passet];
		    }

		    if ($v[$passet + $pass] eq "c" || $v[$passet + $pass] eq "C"){
			
			$puntuacio2 = $puntuacio2 + $pesos{"C"}[$passet];
		    }
		    
		    if ($v[$passet + $pass] eq "g" || $v[$passet + $pass] eq "G"){
			
			$puntuacio2 = $puntuacio2 + $pesos{"G"}[$passet];
		    }
		    
		    if ($v[$passet + $pass] eq "t" || $v[$passet+ $pass] eq "T"){
			
			$puntuacio2 = $puntuacio2 + $pesos{"T"}[$passet];
			
		    }
		    
		    $passet = $passet + 1;
		}

		if ($puntuacio2 >= $millor){

		    $vegades = $vegades + 1;
		    $puntuacio2 = 0; 
		}
		
		$pass = $pass +1;
		
	    }
 
	    if ($vegades > 0){
		$r = $r + 1;

	    }

	      		
	    $puntuacio2 = 0;
	    $h = $h +1;
	    }
	    
	    $pvalue = $r/100;


print "El p-value associat a aquest factor de transcripció és igual o més petit que $pvalue.\n\n";
	
    
	$num_files = 0;	#tornem a posar la matriu a l'inici per fer tantes voltes com matrius
	@v2 = ();
	@v3 = ();
	@v4 = ();
	@v5 = ();

	} #tanca el while de les 100 vegades

}	

#--------------------------------------------fi punt 4--------------------------------------------------------------

	       

close(MATRIUS);


