#!/usr/bin/perl -w

use strict;

# -------------------- subprograma 1 ---------------------------------------------------------

open (FH,"<$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 $i;					#files del hash
my $j;					#columnes del hash
my $posicio = 0;			#portara el compte de quina és la longitud de la matriu 					horitzontalment (anirà avançant per cada fila, així, 						podrem determinar si les matrius són de 7, de 8 o de 9). 						Aquesta variable l'arrossegarem al llarg dels 4 					subprogrames.
my $compta = 0;				#portara el compte de les files que hi ha a la matriu 						(seran 4, una pel nt A, una pel C, una pel G i una pel T)
my $name;				# nom de les matrius


while (<FH>){

					# inicialment farem l'expresio regular per tal de que 						reconegui el nom de la matriu, seguidament ho farem perquè 						ens reconegui els números, per útlim escriurem l'expresio 						regular perquè reconegui el pas d'una matriu a l'altre
    if ($_ =~ m/FA/){

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

    if ($_=~ m/(\d+)\t(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/){

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

	$posicio = $posicio +1;
	$compta = $compta +1;

    }

    if (($_ =~ m/(\/\/)\Z/) && ( $compta > 1)){
							#seguidament escriurem les instruccions 							per printar la matriu que ens dona l'arxiu 								de matrius a la terminal
	my @k= ("A","C","G","T");
	my $i = 0;

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

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

	    while ($j < $compta){

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

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

	$compta = 0;


# ---------------------------------- subprograma 2 ---------------------------------------------


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]);
						# intorduïm un nou hash (%pesos) per tal 							d'emplenar-lo amb el resultat obtingut de la 							matriu de pesos

						#en la primera part del programa calculem la 							proporció d'ocurrències de cada nt al llarg de 							tota la seq promotora.
my $ARGV;
open (PROMOTOR,"<$ARGV[1]") || die();
						#a través del ARGV, li diem que ens obri el fitxer 							que es troba en la segona posició després d'haver 							cridat el programa desde la terminal,per tant, 							ens obrirà la seqüència promotora, que és el que 							ens interessa en aquest moment.
my $a = 0;
my $c = 0;
my $g = 0;
my $t = 0;

my $pa;
my $pc;
my $pg;
my $pt;

my @nt;
my $seq;

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

    else {
	chomp $_;
	$seq = $seq . $_;
	}
}

@nt = split(//,$seq);

  	for (my $ii=0;$ii<=$#nt;$ii++){
	    if ($nt[$ii] eq "a") {$a = $a+1;}
	    if ($nt[$ii] eq "c") {$c = $c+1;}
	    if ($nt[$ii] eq "g") {$g = $g+1;}
	    if ($nt[$ii] eq "t") {$t = $t+1;}
	}					
						# dins el "for"hi posem les condicions invariants: 							$i és la variable que anira recorrent les 							diferents posicions del vector (per tant, dels nt 							de la seqüència promotora), li diem que $i vagi 						sumant consecutivament ($i++) fins que sigui menor 							o igual que la longitud de la seqüència promotora, 							del vector "nt". Posem $# per indicar-li que la 						posició "i" ha d'anar recorrent tota la seqüència 							fins al final.
						# introduïm un if per cada nt diferent, així, si 							troba aquest nt a la seqüència promotora, sumara 1 							a la suma total d'aquests
						# seguidament, fem la divisió de la quantitat de 							cada nucleòtid per la longitud total de la 							seqüència promotora, per poder-ho aplicar a la 							fórmula que hem de resoldre per aconseguir el 							resultat de la matriu de pesos.
	$pa = $a/scalar(@nt);
	$pc = $c/scalar(@nt);
	$pg = $g/scalar(@nt);
	$pt = $t/scalar(@nt);
					# en la segona part del programa calculem la proporció de 						cada nt respecte el total que veiem a la matriu

					# per a cada nt (A,T,G,C) fem el mateix càlcul, per tant, 						tot ens pot anar dins el mateix "for": dins la condició 					del for posem el contador $i a 0, i li diem que	ens vagi 						recorrent les diferents posicions del hash de matrius
					# delcarem una nova variable, $opera (d'operacio), que 						prendrà el valor de $i (per tant, a cada volta anirà 						variant) i ens servirà per realitzar l'operació que hem de 						fer per cada cas.

					# finalment, introduïm el resultat obtingut, en la 						variable $pesos, que s'enregistrarà a una determinada 						posició de la matriu
for(my $ii=0;$ii<$posicio;$ii++){

    my $divident = $hash{"A"}[$ii] + $hash{"C"}[$ii] + $hash{"G"}[$ii] + $hash{"T"}[$ii];
    
    my $opera = $hash{"A"}[$ii];
    my $pes;
    
    if ($opera == 0){
	$pes = -999;	
    }
    else{
	$pes = log($opera/$divident) - log($pa);
    }
    $pesos{"A"}[$ii] = $pes;



    $opera = $hash{"C"}[$ii];
    
    if ($opera == 0){
	$pes = -999;
    }
    else{	
	$pes = log($opera/$divident) - log($pc);
    }
    $pesos{"C"}[$ii] = $pes;	



    $opera = $hash{"G"}[$ii];
    

    if ($opera == 0){	
	$pes = -999;
    }
    else{	
	$pes = log($opera/$divident) - log($pg);
    }
    $pesos{"G"}[$ii] = $pes;

    
    $opera = $hash{"T"}[$ii];
    
    if ($opera == 0){	
	$pes = -999;
    }
    else{	
	$pes = log($opera/$divident) - log($pt);
    }
    $pesos{"T"}[$ii] = $pes;	
}
						# seguidament, li diem que ens printi la matriu 						obtingua per pantalla:
my @keys = keys(%pesos);
my $iii = 0;

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

    while ($j < $posicio){
	printf("\t%.2f", $pesos{$keys[$iii]}[$j] );
	
	$j = $j + 1;
    }

    print "\n";

    $iii = $iii + 1;
}
   		
# ------------aplicacio de la matriu de pesos a la sequencia promotora-----------------------

# ---------------------------subprograma 3---------------------------------------------------


my $s = 0;
my $pos =0;

for (my $i=0;$i<=$#nt-$posicio;$i++){

	my $score = 0;

	for (my $j=0;$j<$posicio;$j++){
		if (($nt[$i+$j])eq "a"){
			$score = $score + $pesos{"A"}[$j];
		} 
		if (($nt[$i+$j])eq "c"){
			$score = $score + $pesos{"C"}[$j];
		} 
		if (($nt[$i+$j])eq "g"){
			$score = $score + $pesos{"G"}[$j];
		} 
		if (($nt[$i+$j])eq "t"){
			$score = $score + $pesos{"T"}[$j];
		} 	
	}

	if ($score>$s){
		$s = $score;
		$pos = $i;
	}
}

print "score maxim: $s\n";
print "posicio: $pos\n";



#-----------------------------------subprograma 4---------------------------------------------

my $pvalue;
my $count = 0; 			#variable que anirà comptant les vegades que l'score obtingut es 					superior o igual a l'score trobat de la nostra seqÜència promotora
 
for (my $x=0;$x<100;$x++){  
	my $n = scalar(@nt);
        my $i = $n - 1;
        while ($i >= 0) {
          my $j = int(rand($i+1));
          if ($i != $j) {
            my $tmp = $nt[$i];
            $nt[$i] = $nt[$j];
            $nt[$j] = $tmp;
          }
          $i = $i - 1;
	  
        }
	
	my $maxim = 0;
	my $score;

	for (my $i=0;$i<=$#nt-$posicio;$i++){
	
		$score = 0;

		for (my $j=0;$j<$posicio;$j++){
		    
			if (($nt[$i+$j])eq "a"){
				$score = $score + $pesos{"A"}[$j];
			} 
			if (($nt[$i+$j])eq "c"){
				$score = $score + $pesos{"C"}[$j];
			} 
			if (($nt[$i+$j])eq "g"){
				$score = $score + $pesos{"G"}[$j];
			} 
			if (($nt[$i+$j])eq "t"){
				$score = $score + $pesos{"T"}[$j];
			} 	
		}
		if ($score>$maxim){
			$maxim = $score;
		}
	}
	if ($maxim>=$s){
		$count = $count +1;
	}			

}

$pvalue = $count/100;


print "El p-value associat al FT és $pvalue\n";
						#com mes petit sigui el valor de p-value vol dir 							que el FT té més probabilitat d'unir-se a la 							seqüència promotora, i que aquesta unió NO sigui a 							l'atzar.

# ------------------------------------------------------------------------------------------------# reinicialitzacio de variables per tal de poder fer el càlcul amb les 13 matrius (els 13 factors de transcripció que tenim inicialment)

%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]);

	$posicio=0;


	
    }
}

close (PROMOTOR);
close(FH);
			# finalment, tanquem els FH oberts a l'inici del programa	
