#!/usr/bin/perl -w

use strict;

#Programa que agafa un fitxer amb matrius de TFs seleccionats i un altre fitxer amb sequencia de regio promotora i mostra per pantalla llista dels TF amb un p-value, que com mes baix sigui més probable es que el TF s'uneixi a la regio.

#agafar la sequencia en format FASTA del fitxer i posar-la en el vector @vector

my @vector;
my $string;
my $letra;

open (PR,"<$ARGV[1]");					#obrir fitxer i donar-li un filehandle
print @vector;
while (<PR>){
	if (! ($_ =~ ">20")){
	chomp($_);
	my $linea = $_; 				#fica una linia de la sequencia en $linia
	my $longitud = length($linea);			#calcula la longitut de la linia

		while( $longitud != "0"){	
			$letra = substr($linea,0,1);  	#fica la 1ª lletra de la linea en $letra
						
			$linea = substr($linea,1,$longitud-1); #$linea es $linea - la lletra usada
			push(@vector,$letra); #fica la lletra dins de @vector
			$longitud = length($linea); #calcula la longitut de la nova linia
		}
	}
}
my $scalar;						
$scalar = scalar(@vector);				#longitut total de la sequencia promotora
close (PR);


#agafar matrius d'ocurrencies dels TF, convertir-les en matrius de pesos

open (FH,"<$ARGV[0]");				#obrir fitxer i donar-li un filehandle

my $w = 0;
my $nom = 0;
my $posicions = 0;

my $nA = 0;
my $nC = 0;
my $nG = 0;
my $nT = 0;
my $total = 0;
my $propA = 0;
my $propC = 0;
my $propG = 0;
my $propT = 0;	
				
				#declaracio del hash per emmagatzemar les matrius que llegeix
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]);

while (<FH>){				#a cada volta llegeix una linia del fitxer
chomp ($_);
 
	if ($_ =~ m/FA\s(.+)\s/){	#captura el nom de la matriu amb l'expressio regular
		$nom = $1;
		$w = 0;
	}	
    		
	if (m/0[0-9]\s(\d+)\s(\d+)\s(\d+)\s(\d+)/){	#captura valors d'ocurrencies de la linia  								#i els assigna a les posicions dels  								#vector del hash i a mes suma les	
							#ocurrencies de cada base 		
			$motiu{"A"}[$w] = $1; 
		    	$nA = $nA + $1;
			$motiu{"C"}[$w] = $2; 
			$nC = $nC + $2;
			$motiu{"G"}[$w] = $3; 
			$nG = $nG + $3;
			$motiu{"T"}[$w] = $4; 
 			$nT = $nT + $4;
				
			$w = $w + 1;
				
						       #per saber el nº de posicions de cada matriu
			if ($1 + $2 + $3 + $4 != 0){   
				$posicions = $posicions + 1;
			}
	}	
	
	if ($_ eq '//' && $nom ne '0'){ 			#ha acabat de llegir la matriu
		print "$nom\t";

		my $seqs = $1 + $2 + $3 + $4;		#calcula el nº de sequencies de la matriu
		$total = $nA + $nC + $nG + $nT;
		$propA = $nA / $total;			#calcula proporcions de cada nucleotid
		$propC = $nC / $total;
		$propT = $nT / $total;
		$propG = $nG / $total;
		
						#defineix hash per emmagatzemar les proporcions
		my %proporcions = ("A", $propA, "C", $propC, "G", $propG, "T", $propT);
		
						
						#converteix matriu d'ocurrencies en matriu de pesos
		my $x = 0;
		my @k = keys(%motiu);		
		while ($x < scalar (@k)){
    			my $j = 0;
    			#print "$k[$x]:";
    	
			while ($j < $posicions){
				
				if ($motiu{$k[$x]}[$j] == 0){
					$motiu{$k[$x]}[$j] = -999;
				}
				else{
				$motiu{$k[$x]}[$j] = log($motiu{$k[$x]}[$j]/$seqs) - log($proporcions{$k[$x]});	
				}
			#print "\t$motiu{$k[$x]}[$j]";
			$j = $j + 1;
			}
    			#print "\n";
    			$x = $x + 1;	
		}    #ja l'ha transformat en matriu pesos

		#ara recorre la sequencia promotora i amb la matriu calcula les puntuacions per 		cada possible lloc d'unió
 							
		my $i = 0; 			#posicio dins la seq promotora			
		my $p = 0;			#posicio dins el TF		
		my $score = 0;
		my $letraA = 0;
		my $letraC = 0;
		my $letraG = 0;
		my $letraT = 0;
		my $imax = 0; 			#posicio (i) on esta score max
		my $smax = -999999999999999999; #score maxim

		while ($i < ($scalar-$posicions+1)){	#recorre la sequencia
			$p= 0;
			$score = 0;
	
			while ($p < $posicions){	#recorre les posicions de la matriu
	
							#puntua cada posicio
				if ($vector[$p+$i] eq 'G'){
					$letraG = $motiu{"G"}[$p];
					$score = $score + $letraG;
				}		
				if ($vector[$p+$i] eq 'A'){
					$letraA = $motiu{"A"}[$p];
					$score = $score + $letraA;
				}			
				if ($vector[$p+$i] eq 'C'){
					$letraC = $motiu{"C"}[$p];
					$score = $score + $letraC;
				}
				if ($vector[$p+$i] eq 'T'){
					$letraT = $motiu{"T"}[$p];
					$score = $score + $letraT;
				}
				$p = $p + 1;
				
			}	
					#quan te la puntuacio del TF en una posicio de la seq, 	 						#compara amb score max i si es mes gran se la queda
			if($score > $smax){
				$smax = $score;
				$imax = $i;
			}
			$i = $i + 1;	
			
		}
		print "Binding site: $imax-".($imax+$posicions-1)."\t"."Score: $smax\t";
		#ha acabat de recorrer tota la seq i ha emmagatzemat la puntuacio més alta	
		
		#per calcular el p-value generem 100 permutacions a l'atzar dels nucleòtids de la 			#regio promotora i en puntuem els possibles llocs d'unio del TF amb el mateix 			#algorisme que l'utilitzat per la regio promotora	
				
				#el seguent algorisme permuta a l'atzar els elements de @vector
		my @v = @vector;
		my $count = 0;
		my $c = 0;
		my $p_value = 0;
		
		while ($count <= 100){
     			
			
			my $n = scalar(@v);
        		my $q = $n - 1;
        		
			while ($q >= 0) {
          			my $j = int(rand($q+1));
          			
				if ($q != $j) {
            				my $tmp = $v[$q];

            				$v[$q] = $v[$j];
            				$v[$j] = $tmp;
         		 	}
          			$q = $q - 1;
			}
			
			my $SCALAR = scalar(@v);
			my $I = 0; 				
			my $P = 0;				
			my $SCORE = 0;
			my $LETRAa = 0;
			my $LETRAg = 0;
			my $LETRAc = 0;
			my $LETRAt = 0;
			my $IMAX_rand = 0; 			     #posicio (i) on està score max
			my $SMAX_rand = -9999999999999999999999999;  #Score maxim

			while ($I < ($SCALAR-$posicions+1)){
				
				$P= 0;
				$SCORE = 0;	
				while ($P < $posicions){
			
					if ($v[$P+$I] eq 'G'){
						$LETRAg = $motiu{"G"}[$P];
						$SCORE = $SCORE + $LETRAg;
					}	
					if ($v[$P+$I] eq 'A'){
						$LETRAa = $motiu{"A"}[$P];
						$SCORE = $SCORE + $LETRAa;
					}			
					if ($v[$P+$I] eq 'C'){
						$LETRAc = $motiu{"C"}[$P];
						$SCORE = $SCORE + $LETRAc;
					}	
					if ($v[$P+$I] eq 'T'){
						$LETRAt = $motiu{"T"}[$P];
						$SCORE = $SCORE + $LETRAt;
					}	
					$P = $P + 1;
				}
				if($SCORE > $SMAX_rand){
					$SMAX_rand = $SCORE;
					$IMAX_rand = $I;
				}
				$I = $I + 1;	
			
			} 		#ha calculat score maxim per la permutacio
					
					#si la puntuacio es > o = que la puntuacio maxima del TF en
					#la regio promotora suma al comptador c
			if ($SMAX_rand >= $smax){
				$c = $c + 1;
			}					
			
			$count = $count + 1;
			
			
		} 
		
		$p_value = $c/100;			#calcula el p-value
			print "p-value: $p_value\n";
		
		$c = 0;
		$posicions = 0;  
		%motiu = ();
	}

	
} 
close (FH);	
