#!/usr/bin/perl -w

use strict;

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


my %hash = ("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 $pos =0; # posicion
my $line = 0;# lineas de la matriz
my $name; # nombre de la matriz
my $i;# variable que recorrera las filas del motivo
my $j;#variable que recorrera las columnas del motivo

#abrimos fichero

while (<FITXER>){

#PRIMERA PARTE


#expresion regular que nos reconoce el nombre de la matriz para luego imprimir el nombre por pantalla
    if ($_ =~ m/FA/){

	$name = $_;
	chomp ($name);
	print "$name\n";
    }

# expresion regular que nos asocia las posiciones del motivo a una variable determinada
    if ($_=~ m/(\d+)\t(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/){

	
	$hash{"A"}[$pos]= $2;
	$hash{"C"}[$pos]= $3;
	$hash{"G"}[$pos]= $4;
	$hash{"T"}[$pos]= $5;

	$pos = $pos +1;
	$line = $line +1;

    }

#expresion regular que nos determina el inicio y final que cada matriz
    if (($_ =~ m/(\/\/)\Z/) && ( $line > 1)){

#imprimimos los hashes de vectores por la terminal

	my @k= ("A","C","G","T");
	my $i = 0;

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

	    my $j = 0;
	    print $k[$i];

	    while ($j < $line){

		print "\t $hash{$k[$i]}[$j]";
		$j = $j + 1;
	    }

	    print "\n";
	    $i = $i + 1;
	}




	$line = 0;



#SEGUNDA PARTE, introducimos un hash nuevo  vacío para luego meter los valores de los pesos


my %pesos = ("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 $ARGV;
open (FITXER2,"<$ARGV[1]") || die();

#abrimos nuestro segundo archivo, la secuencia del promotor
#ponemos los contadores para cada nucleotido a 0

my $a = 0;
my $c = 0;
my $g = 0;
my $t = 0;

#variables para la proporcion de cada nucleotido respecto a la secuencia entera, necesario para luego realizar el logaritmo y el calculo de los pesos.
my $pa;
my $pc;
my $pg;
my $pt;

my @nucleotide;
my $secuence;

while (<FITXER2>) {
   
    if(/>/){}

    else {
	chomp $_;
	$secuence = $secuence . $_;
	}
}
# convertimos la secuencia en un vector para utilizar los nucleotidos como variables independientes

@nucleotide = split(//,$secuence);

#el for nos sirve igual que el while, con las condiciones iniciales determinadas.La variable x nos recorre las posiciones del vector de la secuencia del promotor, no puede ser mas grande que la secuencia total , ponemos el contador a cero, y que avance una a una, posicion por posicion.  
# primero contamos por sepaardo cuantos nucleotidos hay a lo largo de la secuenciaue dicen,  si encuentras una posicion del vector q corresponde a alguno de los cuatro nucleotidos suma mas 1 al contador respectivo que inicialmente estaba a 0
.
  	for (my $x=0;$x<=$#nucleotide;$x++){
	    if ($nucleotide[$x] eq "a") {$a = $a+1;}
	    if ($nucleotide[$x] eq "c") {$c = $c+1;}
	    if ($nucleotide[$x] eq "g") {$g = $g+1;}
	    if ($nucleotide[$x] eq "t") {$t = $t+1;}
	}

# hacemos el calculo de la proporción, dividido entre el numero total de nucleotidos de la secuencia

	$pa = $a/scalar(@nucleotide);
	$pc = $c/scalar(@nucleotide);
	$pg = $g/scalar(@nucleotide);
	$pt = $t/scalar(@nucleotide);



#Ootro bucle , donde x recorrera las posiciones dentro del hash inicial. Como condiciones iniciales, ponemos el contador a 0, determinamos que no sobrepase el numero total de posiciones, que va entre 6 y 9, y que vaya sumandose mas uno cada vez que entre en el bucle


for(my $x=0;$x<$pos;$x++){
    my $peso;

#$suma, suma las ocurrencias de cada nucleotido  para una misma posicion dentro del hash de vectores inicial
    my $suma = $hash{"A"}[$x] + $hash{"C"}[$x] + $hash{"G"}[$x] + $hash{"T"}[$x];
  
    my $value = $hash{"A"}[$x];
# en el enunciado, para facilitar l atarea, nos recomendaron poner -999 en el caso que en una posicion de la matriz hubiera un 0   
    
    if ($value == 0){
	$peso = -999;	
    }
#hacemos el calculo, la resta de logaritmos, como nos indica el enunciado
    else{
	$peso = log($value/$suma) - log($pa);
    }
# el valor obtenido lo pasamos a otra variable, que es la que se introducirá dentro de la matriz de pesos
    $pesos{"A"}[$x] = $peso;


# hacemos lo mismo para C, G y T
    $value = $hash{"C"}[$x];
    
    if ($value == 0){
	$peso = -999;
    }
    else{	
	$peso = log($value/$suma ) - log($pc);
    }
    $pesos{"C"}[$x] = $peso;	



    $value = $hash{"G"}[$x];
    

    if ($value == 0){	
	$peso = -999;
    }
    else{	
	$peso = log($value/$suma) - log($pg);
    }
    $pesos{"G"}[$x] = $peso;

    
    $value = $hash{"T"}[$x];
    
    if ($value == 0){	
	$peso = -999;
    }
    else{	
	$peso = log($value/$suma) - log($pt);
    }
    $pesos{"T"}[$x] = $peso;	
}




#ahora queremos que nos lo muestre por pantalla, para comprobar que todo sale bien, pero no lo imprimiremos cuando tengamos el programa final, dejaremos el print como un comentario 

my @v = keys(%pesos);
my $y = 0;

while ($y < scalar(@v)){
    
    my $j = 0,
    print "$v[$y]:" ;

    while ($j < $pos){
	printf("\t%.2f", $pesos{$v[$y]}[$j] );
	#print "\t$pesos{$k[$iii]}[$j]";
	$j = $j + 1;
    }

    print "\n";

    $y = $y + 1;
}
   

#TERCERA PARTE

# En esta parte queremos obtener el score maximo, y para ello compararemos los resultados que iremos obteniendo succesivamente. 
# Consiste en ir rastreando la secuencia con ventanas de varias posiciones, como en el dibujo del enunciado. La primera ventana tendra los  primeros nucleotidos de la secuencia, calcularemos el score, segun la formula que nos dan en el enunciado.
#Luego corremos la ventana una posicion, es decir la segunda ventana empieza por el segundo nucleotido de la secuencia . La tercera ventana, empieza por el tercero ...etc
#/i/ ira recorriendo cada nucleotido cada vez que tengas que calcular el score y /j/ cambiara cada vez que corras la ventana
#el $max es nuestro score maximo ,tras haber rastreado la secuencia, lo que pide la pregunta.  
 #inicializamos  los contadores a 0

my $max = 0;
my $position =0;
# como condiciones iniciales queremos que la i no llegue hasta la ultima posicion de la secuencia, sino hasta la ultima posicion desde donde se pueda crear una ventana.
for (my $i=0;$i<=$#nucleotide-$pos;$i++){

	my $score = 0;

	for (my $j=0;$j<$pos;$j++){
		if (($nucleotide[$i+$j])eq "a"){
			$score = $score + $pesos{"A"}[$j];
		} 
		if (($nucleotide[$i+$j])eq "c"){
			$score = $score + $pesos{"C"}[$j];
		} 
		if (($nucleotide[$i+$j])eq "g"){
			$score = $score + $pesos{"G"}[$j];
		} 
		if (($nucleotide[$i+$j])eq "t"){
			$score = $score + $pesos{"T"}[$j];
		} 	
	}
#Comparamos el score obtenido en cada caso  con el maximo obtenido,  y obtendremos un maximo y la posicion desde donde se calculo 
#si uno de los scores que obtienes es superior al maximo score que tenias hasta ahora, entonces denominamelo $max. Y denominame $position el numero de la posicion donde has obtenido este score que se convirtio en maximo.

	if ($score>$max){
		$max = $score;
		$position = $i;
	}
}

print "score maxim: $max\n";
print "posicio: $position\n";



#CUARTA PARTE

# Esta pregunta consiste en obtener el maximo del maximo, en este caso lo llamamos $superior, que lo obtendremos al realizar una permutacion de la secuencia 100 veces, es decir, obtener 100 secuencias al azar
#$superior, es un marcador que tiene que empezar en 0, como antes con el max, porque habra que ir comparando con el resultado anterior, si es inferior, no pasara nada, pero si es superior pasara de ser $maximum a $superior


my $pvalue;
my $superior = 0;


# Trozo de programa cedido por el enunciado para conseguir permutar la secuencia 
for (my $x=0;$x<100;$x++){  
	my $n = scalar(@nucleotide);
        my $i = $n - 1;
        while ($i >= 0) {
          my $j = int(rand($i+1));
          if ($i != $j) {
            my $tmp = $nucleotide[$i];
            $nucleotide[$i] = $nucleotide[$j];
            $nucleotide[$j] = $tmp;
          }
          $i = $i - 1;
	  
        }
#Aqui hacemos lo mismo que antes, calcular el score de cada permutacion.
	
	my $maximum = 0;
	my $score;

	for (my $i=0;$i<=$#nucleotide-$pos;$i++){
	
		$score = 0;

		for (my $j=0;$j<$pos;$j++){
		    
			if (($nucleotide[$i+$j])eq "a"){
				$score = $score + $pesos{"A"}[$j];
			} 
			if (($nucleotide[$i+$j])eq "c"){
				$score = $score + $pesos{"C"}[$j];
			} 
			if (($nucleotide[$i+$j])eq "g"){
				$score = $score + $pesos{"G"}[$j];
			} 
			if (($nucleotide[$i+$j])eq "t"){
				$score = $score + $pesos{"T"}[$j];
			} 	
		}
#Primero comparamos los diferentes scores de las 100 permutaciones,cada una obtenida por el azar. El maximo score obtenido tras ir comparando uno a uno,se convertira en $maximum

		if ($score>$maximum){
			$maximum = $score;
		}
	}

#Finalmente comparamos a ver si $maximum es superior o igual al score maximo obtenido con nuestra secuencia del paso anterior, $max
	if ($maximum>=$max){
		$superior = $superior +1;
	}			

}

$pvalue = $superior/100;
# Cuanto mas pequeño sea p-value,más probabilidad de que no sea al azar.

# Porfin pedimos que nos printe el resultado final....el ultimo paso conseguido.

print "El p-value associado al factor de transcripcion es  $pvalue\n";










# ponemos a 0 las variables.

%pesos = ("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;


	
    }
}
#cerramos nuestros ficheros para finalizar nuestro programa y de esta forma ejecutarlo.
close (FITXER2);
close(FITXER);