#!/usr/bin/perl -w

use strict;


##################### 1ª parte del programa #################################




my %matriz;  #hash#
my $linea;   #variable que recorre las lineas del archivo#
my $a = 0;   #variable que determinara las posiciones dentro de los vectores donde guardaremos las lineas#
my @v;       #vector donde guardamos temporalmente los elementos de una linea#
my @v2;      #vectores de los valores que corresponderan a cada palabra clave del hash#
my @v3;
my @v4;
my @v5;
my $contador = 0;  #variable que cuenta las lineas#
my $nom;     #variable donde guardaremos el identificador de cada factor de transcripcion#


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

while (<FILE>){     
 
    chomp;
    $linea = $_;
    
###Guardamos el nombre de los identificadores de los factores de transcripcion###

    if($linea =~ m/FA\s+(.*)/)
    {
	$nom = $1;
	print "NOM $nom\n";
	@v2 = ();
	@v3 = ();
	@v4 = ();
	@v5 = ();
    }
    
###Cogemos solo las lineas que empiecen por 0-digito###

    if ($linea =~ m/^(0)\d+\s+\d*\s+\d*\s+\d*\s+\d*/){
	
	@v = split(/\t/,$linea);
	
	#Hacemos corresponder los valores introducidos en el vector v (contenido de la linea) en una posicion determinada de cada vector#

	$v2[$a] = $v[1];
	$v3[$a] = $v[2];
	$v4[$a] = $v[3];
	$v5[$a] = $v[4];
	

	$a = $a + 1; 
	
		
        #construccion del hash

	my @k = keys(%matriz);
	
	
	$contador = $contador + 1;
	
    }
    
        #Hacemos que pare de llenar vectores cuando encuentre //, sin contar las de las primeras lineas#
    
    
    if (($linea =~ m/(\/\/)\Z/) && ($contador > 1)){
	
	#Llenamos el hash con los vectores, asociandolos a la palabra clave correspondiente#

	%matriz = ("A" => \@v2,
		   "C" => \@v3,
		   "G" => \@v4,
		   "T" => \@v5);
	$a = 0;

	my @k = keys(%matriz);   
  
	#programa para printar el hash#

	my $x = 0;

	while ($x < scalar(@k)){

	    close(FILE2);
	    	       	    
	    $x = $x + 1;
	}

	print "\n";
	
	$x=0;
	while($x < 4)

	{
	    my $y = 0;
	    while ($y < scalar(@v2)){
		
			
		$y = $y + 1;
	    }

	    $x = $x+1;
	    
	}

		

#####################2ª parte programa#######################
	

#Cálculo de pb(x):
	
	open(FILE2,"<$ARGV[1]");
	
	my $secuencia; #variable donde introduciremos la secuencia promotora#
	my @a;         #vectores donde incluiremos cada nucleotido#
	my @t;
	my @c;
	my @g;
	
	my $cuentaa;  #variables que cuentan el numero de ocurrencias de cada nucleotido#
	my $cuentat; 
	my $cuentac;
	my $cuentag;
	
	my $sumaa = 0; #variable en la que acumularemos la suma de las ocurrencias de cada nucleotido#
	my $sumat = 0;
	my $sumac = 0;
	my $sumag = 0;
	
	my $sumatotal = 0; #variable de la suma total#
	
	my $propa = 0;   #variables en la que almacenaremos el resultado de la proporcion de cada uno de los nucleotidos#
	my $propt = 0;
	my $propc = 0;
	my $propg = 0;
	

	while (<FILE2>){
	    
	    #quitamos el identificador de la secuencia del archivo FASTA#

	    if (/\>/){
		next;  }
	    
	    chomp;

	    #cambiamos los nucleótidos escritos en minuscula por mayuscula#

	    s/a/A/g;
	    s/c/C/g;
	    s/g/G/g;
	    s/t/T/g;
	



	    #concatenamos toda la secuencia promotora para evitar los saltos de linea#

	    $secuencia = $secuencia.$_;
	}
	
	
	   #mediante una expresion regular buscamos cada nucleotido en la secuencia y los introducimos en su vector correspondiente#


	@a = ($secuencia =~ m/(A)/g);
	@t = ($secuencia =~ m/(T)/g);
	@c = ($secuencia =~ m/(C)/g);
	@g = ($secuencia =~ m/(G)/g);
	
	
	  #contamos el numero de correspondencias de cada nucleotido contando el numero de posiciones de cada vector#

	$cuentaa = scalar(@a);
	$cuentat = scalar(@t);
	$cuentac = scalar(@c);
	$cuentag = scalar(@g);
	
	 #contamos ocurrencias de cada nucleotido para toda la secuencia#

	$sumaa = $sumaa + $cuentaa;
	$sumat = $sumat + $cuentat;
	$sumac = $sumac + $cuentac;
	$sumag = $sumag + $cuentag;
	
	 #sumamos las ocurrencias de todos los nucleotidos para poder hacer la proporcion de cada uno#

	$sumatotal = $sumaa + $sumac + $sumag + $sumag;
	
	$propa = $sumaa / $sumatotal;
	$propt = $sumat / $sumatotal;
	$propc = $sumac / $sumatotal;
	$propg = $sumag / $sumatotal;
	

	close(FILE2);
	

#calculamos pim(x)#
	
	
	my $total;
	my $sumapim = 0;
	my $divpim = 0;
	my $pos1=0;

	#recorremos las posiciones del hash cambiando los valores iguales a 0 por -999#
	
	while ($pos1 < scalar(@v2)){
	    
	    if ($matriz{"A"}[$pos1] == 0){
		$matriz{"A"}[$pos1] = -999;
		
	    }
	    
	    #en el caso de que el valor no sea 0, calculamos la proporción de cada nucleotido en una posicion concreta y calculamos el valorde esa posicion en la matriz de pesos, de forma que sustituimos los hashes anteriores por matrices de pesos#


	    else{
		
		$sumapim = (${$matriz{"A"}}[$pos1]) + (${$matriz{"T"}}[$pos1]) + (${$matriz{"G"}}[$pos1]) + (${$matriz{"C"}}[$pos1]);
		$divpim = (${$matriz{"A"}}[$pos1]) / $sumapim;
		
		${$matriz{"A"}}[$pos1] = log $divpim - log $propa;
		
		
	    }
	    
	    	    
	    
	    if (${$matriz{"C"}}[$pos1] == 0){
		
		${$matriz{"C"}}[$pos1] = -999;
		
	    }
	    
	    else {
		
		$sumapim = (${$matriz{"A"}}[$pos1]) + (${$matriz{"T"}}[$pos1]) + (${$matriz{"G"}}[$pos1]) + (${$matriz{"C"}}[$pos1]);
		$divpim = (${$matriz{"C"}}[$pos1]) / $sumapim;
		
		${$matriz{"C"}}[$pos1] = log $divpim - log $propc;
		
		
	    }
	    
	    
	    if (${$matriz{"T"}}[$pos1] == 0){
		
		${$matriz{"T"}}[$pos1] = -999;
		
	    }
	    
	    else {
		
		$sumapim = (${$matriz{"A"}}[$pos1]) + (${$matriz{"T"}}[$pos1]) + (${$matriz{"G"}}[$pos1]) + (${$matriz{"C"}}[$pos1]);
		$divpim = (${$matriz{"T"}}[$pos1]) / $sumapim;
		
		${$matriz{"T"}}[$pos1] = log $divpim - log $propt;
		
		
	    }
	    
	    
	    
	    if (${$matriz{"G"}}[$pos1] == 0){
		
		${$matriz{"G"}}[$pos1] = -999;
		
	    }
	    
	    else {
		
		$sumapim = (${$matriz{"A"}}[$pos1]) + (${$matriz{"T"}}[$pos1]) + (${$matriz{"G"}}[$pos1]) + (${$matriz{"C"}}[$pos1]);
		$divpim = (${$matriz{"G"}}[$pos1]) / $sumapim;
		
		${$matriz{"G"}}[$pos1] = log $divpim  - log $propg;
		
		
	    }
	    
	    $pos1 = $pos1 + 1;
	    
	}

	#programa para printar la matriz de pesos#
	
	$x = 0;

	while ($x < scalar(@k)){
	    
	    
	    $x = $x + 1;

	}
	
#	print "\n";
	
	$x=0;
	
	while($x < 4)

	{
	    my $y = 0;

	    while ($y < scalar(@v2)){
		
		
	#	 print "\t$matriz{$k[$x]}[$y]";

		$y = $y + 1;

	    }
	    
	    $x = $x+1;
	    
	  #  print "\n";
	    
	}



############################3ª parte programa####################
	

	my $noi = 0;           #variable que recorre las posiciones de la secuencia promotora#
	my $window;            #variable tamaño de "sliding window"#
       	my $scoremax = -999999999;   #variable de "score" máximo#
	my $posmax = 0;              #posicion donde encontramos el mayor "score"#
	
	$window = scalar(@v2);
	
	#introducimos la secuencia promotora en un vector#

	my @seq = split(//,$secuencia);
	

	#definimos la ventana#

	my $sor;    #variable que recorre las posiciones de "sliding window"#
	
	#recorremos la secuencia#

	while ($noi <= (scalar(@seq) - $window)){
	    
	    $sor= 0;
	    
	    my $s = 0;   #variable "score"#

	    #Miramos la ventana#
	    
	    while ($sor < scalar(@v2)){

		$s = $s+$matriz{$seq[$sor+$noi]}[$sor];

		$sor = $sor+1;		    
	    }
	    
	    #comparamos el "score" obtenido en cada ventana y registramos el maximo y su posicion#
	     	    
	    if ($s > $scoremax){
		
		$scoremax = $s;

		$posmax = $noi;
	    }

	    $noi = $noi +1;
	}
	
	print "la puntuacion maxima es: $scoremax en la posicion: $posmax\n";
	
	
########################4ºparte##############################
	
	#realizamos el mismo procedimiento que en la parte tres pero con una secuencia aleatoria#
	
	my $posmaxr;
       	my $sorr;
	my $cien = 0;    #variable que cuenta el numero de veces que hacemos el bucle#
	my $yasta = 0;   #variable que cuenta el numero de veces que el "score" maximo obtenido aleatoriamente es mayor al maximo obtenido anteriormente con nuestra secuencia promotora#
	
	
	while ($cien < 100){

	    my $scoremaxr = -999999999;
	    
	    my $noir = 0;
	    
	    
	    #permutacion para crear secuencia aleatoria a partir de nuestra secuencia promotora#
	    
	    my @seqr=@seq;
	    
	    my $nr = scalar(@seqr);
	    
	    my $ir = $nr - 1;
	    
	    while ($ir >= 0) {
		
		my $jr = int(rand($ir+1));
		
		if ($ir != $jr) {
		    
		    my $tmp = $seqr[$ir];
		    
		    $seqr[$ir] = $seqr[$jr];
		    
		    $seqr[$jr] = $tmp;
		    
		}
		
		$ir = $ir - 1;
		
	    }
	    
	    #procedimiento parte tres adaptado a las nuevas variables (iguales que las anteriores pero acabadas en "r"#

	    while ($noir <= (scalar(@seqr) - $window)){

		$sorr= 0;
		
		my $sr = 0;
					
		
		while ($sorr < scalar(@v2)){

		    $sr = $sr+$matriz{$seqr[$sorr+$noir]}[$sorr];

		    $sorr = $sorr+1;		    
		}
		
		
		if ($sr > $scoremaxr){
		    
		    $scoremaxr = $sr;

		    $posmaxr = $noir;

		}

		$noir = $noir +1;
		
	    }
	    
	    $cien = $cien + 1;
	    
	    
	    if ($scoremaxr > $scoremax){
		
		$yasta = $yasta + 1;
		
		
	    }
	}
	
	#calculamos el p-value#
	
	my $end = $yasta/100;
	
	print "$end\n";
	
    }
    
    
}