#!/usr/bin/perl -w

use strict;

###############################################################################
#                                                                             #
#    IDENTIFICACIó COMPUTACIONAL DE LLOCS D'UNIó A FACTORS DE TRANSCRIPCIó    #
#                                                                             #
#                               Samanta Yubero                                #
#                                Sonia Laguna                                 #
#                                Manuel Gómez                                 #
#                                                                             #
###############################################################################



###### 1.  LECTURA DELS ARXIUS D'ENTRADA I ENREGISTRAMENT DE LES MATRIUS EN HASHES  ######


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


my %motiu = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0],   #Es crea un hash inicial buit.
	     "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; #Enregistra els noms del FT's.

my $posicio = 0; #Avança la posició dins el vector de cada nucleòtid del hash.

my $comptar = 0; #Compta quantes posicions té el FT.

my $i;

while () {

    if ($_ =~ m/FA./) { #Això reconeix el nom de cada FT.

	$nom = $_;

        chomp($nom);

	print "$nom\n";

    }

    if ($_=~ m/(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\Z/) { #S'enregistren els nombres de la matriu de cada FT al hash.
 
	$motiu{"A"}[$posicio]= $2;

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

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

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

	$posicio = $posicio + 1;

	$comptar = $comptar + 1;

    }
    
    if (($_ =~ m/(\/\/)\Z/) && ($comptar > 1)) { #Reconeix l'inici de cada FT.


######  2. CONSTRUCCIó DE LA MATRIU DE PESOS  ######

   ## Càlcul de la proporció d'ocurr&egrove;ncies de cada nucleòtid a la seqüència promotora ##

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

	my $longitut; 

	my $g = 0;
	
	my $c = 0;

	my $t = 0;

	my $a = 0;

	my $pos = 0; #Recorre el vector de posició en posició (nucleòtid a nucleòtid).


	while () {
    
	    chomp ($_); #Es treuen els canvis de linea.

	    my @v = split(//,$_); #Es transforma la seqüència promotora en un vector.

	    $longitut = scalar(@v); #Calcula la longitut del vector.

	    while ($pos < $longitut) {

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

		    $g = $g + 1;

		}

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

		    $c = $c + 1;

		}

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

		    $t = $t + 1;

		}

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

		    $a = $a + 1;

		}    

		$pos= $pos + 1;
	    }

	    $a = $a/$longitut;
	    $c = $c/$longitut;
	    $t = $t/$longitut;
	    $g = $g/$longitut;

	}

	close (PROMOTOR); 


  ## Construcció de la matriu de pesos ##

	$posicio = 0; #Es reinicia la variable a 0, per anar posició per posició dins el hash de vectors.

	my %pmi = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0],  #Es crea un nou hash per enregistrar la matriu de pesos.
	           "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]);
	
	while ($posicio < $comptar) {

	    if ($motiu{"A"}[$posicio] == 0) { #Assigna una puntuació negativa molt baixa si a la matriu d'ocurrències apareix un 0.

		$pmi{"A"}[$posicio] = -999 - ((log $a)/(log 10));  
	    }

	    else { #Calcula la proporció d'ocurrències de cada nucleòtid al motiu.

		$pmi{"A"}[$posicio]= (log (($motiu{"A"}[$posicio])/(($motiu{"A"}[$posicio])+($motiu{"C"}[$posicio])+($motiu{"G"}[$posicio])+($motiu{"T"}[$posicio])))/(log 10)) - (log ($a)/log(10));
	    }

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

		$pmi{"C"}[$posicio] = -999 - ((log $c)/(log 10));  
	    }

	    else {

		$pmi{"C"}[$posicio]= (log (($motiu{"C"}[$posicio])/(($motiu{"A"}[$posicio])+($motiu{"C"}[$posicio])+($motiu{"G"}[$posicio])+($motiu{"T"}[$posicio])))/(log 10)) - (log ($c)/log(10));

	    }

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

		$pmi{"G"}[$posicio] = -999 - ((log $g)/(log 10)); 
	    }

	    else {
		$pmi{"G"}[$posicio]= (log (($motiu{"G"}[$posicio])/(($motiu{"A"}[$posicio])+($motiu{"C"}[$posicio])+($motiu{"G"}[$posicio])+($motiu{"T"}[$posicio])))/(log 10)) - (log ($g)/log(10));

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

		$pmi{"T"}[$posicio] = -999 - ((log $t)/(log 10));  
	    }

	    else {

		$pmi{"T"}[$posicio]= (log (($motiu{"T"}[$posicio])/(($motiu{"A"}[$posicio])+($motiu{"C"}[$posicio])+($motiu{"G"}[$posicio])+($motiu{"T"}[$posicio])))/(log 10)) - (log ($t)/log(10));

	    }
	
	    $posicio = $posicio + 1;
	}


###### 3.  CàLCUL DE LA PUNTUACIó MàXIMA PER A CADA FT  ######

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

	my $punts = 0; #Enregistrarà la puntuació dins cada finestra.

	my $lloc; #S'utilitza per a iniciar el recompte de punts dins cada finestra.

	my $posmax = 0; #és la posició dins el promotor on la puntuació de la finestra és màxima.

	my $finestra = $comptar; #La finestra serà igual al numero de posicions de cada FT. 

	my $inici = 0; #L'utilitzem per a recorrer el vector.

	my $sestel = -999; #és la puntuació màxima per a cada FT. La igualem a -999 ja que correspondrà amb la unió amb una posició incorrecta.

	my $nucl = 0; #Ens mou dins cada finestra.

	my $y = ; #Enregistra el promotor dins una variable.
	
	chomp ($y);

	my @w = split(//,$y); #Crea un vector a partir de la variable promotor.

	$longitut = scalar(@w);

	while ($inici <= ($longitut - $finestra)) {

	    $punts = 0; #S'inicien a 0 per a cada nova finestra.
	    $lloc = $inici; #El lloc ha de ser igual a l'inici per què coincideixi dins el vector.
	    $nucl = 0; #Ha de ser 0 a cada volta per què es mogui dins la finestra i per què coincideixi la posició de cada nucleòtid de la finestra mamb la posició dins el hash.
 
	    while ($nucl <= ($finestra - 1)) { #La finestra - 1 és per què la finestra no tengui posicions buides al final del promotor.

		if (($w[$lloc] eq 'G') || ($w[$lloc] eq 'g')) {

		    $punts = $punts + ($pmi{"G"}[$nucl]);

		}
		if (($w[$lloc] eq 'C') || ($w[$lloc] eq 'c')) {

		    $punts = $punts + ($pmi{"C"}[$nucl]);

		}
		if (($w[$lloc] eq 'T') || ($w[$lloc] eq 't')) {
		    
		    $punts = $punts + ($pmi{"T"}[$nucl]);
	   
		}
		if (($w[$lloc] eq 'A') || ($w[$lloc] eq 'a')) {

		    $punts = $punts + ($pmi{"A"}[$nucl]); 
		}  
  
		$lloc = $lloc + 1;
		$nucl = $nucl + 1;

	    }

	    if ($punts > $sestel) {

		$sestel = $punts;

		$posmax = $inici;
				  
	    }

	    $inici = $inici + 1;
	    	    	    
	}

	print "s* = $sestel\n";
	print "la posició és $posmax\n";

	close (PROM);

	
###### 4.  CàLCUL DEL P-VALUE  ######

	my $cent = 0; #Per fer 100 voltes.
	my $score = -999; #és la puntuació màxima random
	my $cops = 0; #Cada cop que la puntuació màxima random supera a la puntuació màxima del promotor original.
	my $pvalue; #Enregistra el p-value.
	$i = 0; #es reinicia a 0 la variable, per si ja l'haviem utilitzat.
	my $j = 0;

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

	my $f = ; #Es torna a enregistrar el promotor dins una altre variable.

	chomp ($f);

	while ($cent < 100) { #Ens donarà 100 repeticions.

	    my @z = split(//,$f);
	    my $n = scalar(@z); 
	    my $i = $n - 1; 


  ## Crea seqüències aleatòries del promotor ##

	    while ($i >= 0) { 

		$j = int(rand($i+1)); 

		if ($i != $j) { 

		    my $tmp = $z[$i]; 
		    $z[$i] = $z[$j]; 
		    $z[$j] = $tmp; 
		} 
		$i = $i - 1; 
	    }
   
	    $finestra = $comptar; #Reiniciem les variables.
	    $inici = 0;
	    $nucl = 0;


  ## Recompte de la puntuació per a cada seqüència random ##
	      
	    while ($inici <= ($n - $finestra)) {
		$punts = 0;
		$lloc = $inici;
		$nucl = 0;
		
		while ($nucl <= ($finestra - 1)) {
		    
		    if (($z[$lloc] eq 'G') || ($z[$lloc] eq 'g')) {

			$punts = $punts + ($pmi{"G"}[$nucl]);

		    }

		    if (($z[$lloc] eq 'C') || ($z[$lloc] eq 'c')) {

			$punts = $punts + ($pmi{"C"}[$nucl]);

		    }

		    if (($z[$lloc] eq 'T') || ($z[$lloc] eq 't')) {

	    		$punts = $punts + ($pmi{"T"}[$nucl]);

	   	    }

		    if (($z[$lloc] eq 'A') || ($z[$lloc] eq 'a')) {

			$punts = $punts + ($pmi{"A"}[$nucl]); 

		    }  

  		    $lloc = $lloc + 1;
		    $nucl = $nucl + 1;
		}

		if ($punts > $sestel) {  #Si la puntuació d'una finestra random és major que la puntuació màxima del promotor, feim que es sumi un cop i que ja no miri més en aquella seqüència aleatòria.
		    
		    $cops = $cops + 1;
		    $inici = $n;

		}

		$inici = $inici + 1;
	       
	    }
	   
	    $cent = $cent + 1;

 	}

	$pvalue = $cops/100;
	print "el p-value associat al FT és <= $pvalue\n";
	print "\n";

  ## Si volguessim imprimir les matrius (la de ocurrències ò la de pesos) ##

####################################################
#	my @k= ("A","C","G","T");                  #
#                                                  #
#	$i = 0;                                    #
#                                                  #
#	while ($i < scalar (@k)) {                 #
#                                                  #
#	    my $j = 0;                             #
#                                                  #
#	    print $k[$i];                          #
#                                                  #
#	    while ($j < $compta) {                 #
#                                                  #
#		print "\t $motiu{$k[$i]}[$j]";     # Es substiuiria el $motiu per $pmi per mostrar la matriu de pesos en lloc de la de 
#                                                  #   ocurrències.
#		$j = $j + 1;                       #
#                                                  #
#	    }                                      #
#                                                  #
#	    print "\n";                            #
#                                                  #
#	    $i = $i + 1;                           #
#	}                                          #
####################################################    

	%motiu = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0], #Es reinicien les variables inicials a 0 per poder iniciar el procés del següent FT.
	          "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]);

	$comptar = 0;

	$posicio = 0;
	
    }


}

close (PROMOT);

close (FACTORS);