#!/usr/bin/perl -w

use strict;

my $seq;
my @v; #definim un vector amb la seqüència promotora que hi entrem
my $s = 0; #indica la posició (nucleòtid) que s'està llegint dins de la seqüència promotora
my $ra = 0; #suma les A de la seqüència promotora
my $rc = 0; #suma les C de la seqüència promotora
my $rg = 0; #suma les G de la seqüència promotora
my $rt = 0; #suma les T de la seqüència promotora
my $qa = 0; #proporció de A de la seqüència promotora
my $qc = 0; #proporció de C de la seqüència promotora
my $qg = 0; #proporció de G de la seqüència promotora
my $qt = 0; #proporció de T de la seqüència promotora
my $longvector; #definim la longitud del vector (seqüència promotora)
my $va =0;
my $vc =0;
my $vg =0;
my $vt =0;




open (PROMOTOR,"$ARGV[0]"); #obrim la seqüència promotora

while (<PROMOTOR>){

chomp($_);	

@v = split(//,$_); #definim un vector amb la seqüència promotora

$longvector = scalar(@v); #definim una variable amb la longitud del vector (promotor)
$s=0;

while ($s < $longvector){


    if ($v[$s] eq 'A'){
        $ra = $ra + 1; #suma les A de la seqüència promotora
    }
    if ($v[$s] eq 'C'){
        $rc = $rc + 1; #suma les C de la seqüència promotora
    }
    if ($v[$s] eq 'G'){
        $rg = $rg + 1; #suma les G de la seqüència promotora
    }
    if ($v[$s] eq 'T'){
        $rt = $rt + 1; #suma les T de la seqüència promotora
    }

$s = $s + 1; #incrementem 1 posició del nucleòtid de ls seqüència promotora 


}

}
$qa =($ra / $longvector); #proporció de A de la seqüència promotora

$qc =($rc / $longvector); #proporció de C de la seqüència promotora

$qg =($rg / $longvector); #proporció de G de la seqüència promotora

$qt =($rt / $longvector); #proporció de T de la seqüència promotora


close (PROMOTOR);


my %motiu = ("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]);


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

my $c = 0; #variable per anar contant files
my $p = 0; #variable per contar longitud de la matriu (número de columes)
my $log=0; 
my $z=0;
my $nom; #nom de cada un dels factors de transcripció
my $columnes; #definim una variable que compta columnes
my $i=0; #files
my $j=0; #columnes
my $valor=0; #definim variable per a guardar la suma per cada F.T per poder calcular l'score maxim.
my $FT=""; #t'hi guardes el FT de + score
my $k=0; #variable per poder printar la seqüència que es correspon amb l'score maxim
my $max=-10000; #variable on s'enregistrarà la suma amb el score maxim
my $files; #definim una variable que conta files

while (<FITXERMEU>){
    chomp($_);
    
   if ($_=~ m/^FA(.+)$/){
  
	$nom = $1;
	$p = 0;
    
	} 
	
    elsif ($_ =~ m/\d+\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/){
  
	my $suma_col=$1+$2+$3+$4; #definim una variable per la suma de nucleòtids de cada posició
	my $pa=$1/$suma_col; #definim variable per a la proporció de A en cada posició de la matriu
	my $pc=$2/$suma_col; #definim variable per a la proporció de C en cada posició de la matriu
	my $pg=$3/$suma_col; #definim variable per a la proporció de G en cada posició de la matriu
	my $pt=$4/$suma_col; #definim variable per a la proporció de T en cada posició de la matriu
	
	if($pa==0){
		$motiu{"A"}[$p] =-999;
	}	
	else{
		$motiu{"A"}[$p] = (log($pa)/log(10))-(log($qa)/log(10));
	}
	
	if($pc==0){
		$motiu{"C"}[$p] =-999;
	}	
	else{
		$motiu{"C"}[$p] = (log($pc)/log(10))-(log($qc)/log(10));
	}
	
	if($pg==0){
		$motiu{"G"}[$p] =-999;
	}	
	else{
		$motiu{"G"}[$p] = (log($pg)/log(10))-(log($qg)/log(10));
	}
	
	if($pt==0){
		$motiu{"T"}[$p] =-999;
	}	
	else{
		$motiu{"T"}[$p] = (log($pt)/log(10))-(log($qt)/log(10));
	}
        $p = $p + 1;
		
    } 
    elsif ($_=~ (m/\/\/\Z/) && ($c > 2)){
	
	
	my @fila = ("A","C","G","T");
	$i=0; #definim variable per les files de la matriu
	$j=0; #definim variable per les columnes de la matriu
	$columnes=$p; 
	$files=4;
	print "MATRIU $nom\n\n";
	while($i<$files){
			$j=0;
			print "$fila[$i]";
			while($j<$columnes){
				print "\t$motiu{$fila[$i]}[$j]";
				$j=$j+1;
			}
			print "\n";
	$i=$i+1;
	}
	
	$i=0;
	$FT="";
	$k=0; #variable per poder printar la seqüència que es correspon amb l'score maxim
	$valor=0; #definim variable per les columnes de la matriu
	$max=-1000000; #variable on quedarà enregistrat el valor maxim de $valor.A $max se li assigna inicialment un valor tan negatiu (-1000000) ja que a la matriu de logarítmes hi ha valors molt negatius. 
	
	while ($i < ($longvector-$columnes+1)){
		$j=0;
				
		while ($j < $columnes){
					
			my $lletra = $v[$i + $j];
			#print "$lletra\n";
			#print "$i\n";
			$valor = $valor + $motiu{$lletra}[$j];
								
			
			$j = $j + 1;	
									
 			}

		if($valor > $max){ #perquè es vagi enregistrant el valor d'score mirat sempre i quan sigui major al valor maxim ja enregistrat.
			$FT=""; #guardem el F.T de mes score
			$max=$valor;
			$k=0;

			while($k<$columnes){ 
				$FT=$FT.$v[$i + $k]; #concatenem un nucleòtid rere l'altre de la seqüència que ha donat l'score maxim
				$k=$k+1; 
				}
			}
		$valor=0;
		$i = $i + 1;
			
				
		}
	print "SCORE = $max ($FT)\n\n";
	#print "longvector=$longvector";
		

			
	my $pvalue = 0;  #definim una nova variable que guardarà el valor de p-value 
	my $cc = 0; #l'arxiu dels profes es $pvalue = $c / 100
	my $ran = 0; #num. de repeticions
        
	while ($ran <= 100){ #fem el random fins a 100 vegades
      		my @vec = @v;
		my $n = scalar(@vec); #llargada del vector vec 
		my $ii = $n - 1; #definim nova variable files de la matriu

		while ($ii >= 0) {
          	my $jj = int(rand($ii+1));

          		if ($ii != $jj) {
           	 	my $tmp = $vec[$ii];

           	 	$vec[$ii] = $vec[$jj];
            		$vec[$jj] = $tmp;
          		}
			
         	 $ii = $ii - 1;
		#print "aa = $max\n";
	 	}	

		$ii = 0;
		my $max_ii = -1000000000; #definim nova variable per calcular el valor maxim del random, i que ens permetrà després comparar amb el valor de $max anteriorment enregistrat 
		my $jj; #definim nova variable columnes de la matriu
		my $valor = 0; #suma


		while ($ii < ($longvector-$columnes+1)){
			$jj = 0;
								
			while ($jj < $columnes){
			
				my $lletra = $vec[$ii + $jj];
				
				$valor = $valor + $motiu{$lletra}[$jj];
				$jj = $jj + 1;
					
				}
			
			if($valor > $max_ii){
				$max_ii=$valor;
						
				}
										
			$ii = $ii + 1;
			$valor = 0;
				
		}
				
		if ($max_ii > $max){ #comparem el valor de $max_ii del random amb el $max anterior
			$cc = $cc + 1;
			}
		
		$max_ii = 0;
		$ran = $ran + 1;
	} #tanca fer el ran 100 cops
	
	$pvalue = $cc / 100;
	
	print "El p-value es <= $pvalue\n";
	print "$cc vegades S max a l'atzar es major o igual a S max*\n";
		
	$cc = 0;
	print "\n";
	
			
	%motiu = ("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]);

    }
	$c = $c + 1;
}


close (FITXERMEU);