#!/usr/bin/perl -w

use strict;

#############################################################
###################                   #######################                   
###################       DADES       #######################
###################                   #######################
#############################################################
##                                                         ##
##                                                         ##
##      AUTORES: Natàlia Pérez i Laia Pibernat             ##
##                                                         ##
##                                                         ##
#############################################################



#############################################################
##El programa següent prediu les posicions d'interacció de ##
## factors de transcrició en una seqüència promotora així  ##
##  com mostra el valor de l'score i el p-value de cada    ##
##                     interacció                          ##
#############################################################



############# DECLARACIÓ DE VARIABLES ########################

my %motiu = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0],   #emmagatzema els valors de les matrius dels factors de transcripció proporcionades
             "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 = 0;                                        # 
my @v;                                            #  
my $posicio = 0;                                  #
my $score = 0;                                    #



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


while (<FITXERMEU>){           #anem llegint les línies del document, d'una en una.
   
   chomp($_);


   if (m/(FA)/) {      #expressió regular per a que trobi la primera línia amb l'identificador del factor de transcripció
   
       print "$_\n";   #mostra la primera línia d'informació de cada FT

  }
   

   if (m/\d+\t(\d+)\t(\d+)\t(\d+)\t(\d+)/) { #expressió regular per trobar la fila de números
      

       $motiu{"A"}[$i] = $1;   # assignem els valors de les variables anònimes de les files a la posició corresponent dels vectors del hash
       $motiu{"C"}[$i] = $2;   # primera fila que llegirà, posició 0 dels vectors.
       $motiu{"G"}[$i] = $3;
       $motiu{"T"}[$i] = $4;
       $i = $i +1; 

    }

   if ((m/\/\//) && ($i > 2)){   #condició un cop tenim la hash plena


       ###conversió en matriu de pesos



       ### càlcul de les proporcions de nucleòtids en la seqüència promotora (pB(x))

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

       my $x = <SEQPROM>;
       chomp($x);
       
       my @vector = split(//,$x);
       
       my $recorre = 0;
       my $pa = 0;
       my $pc = 0 ;
       my $pt = 0;
       my $pg = 0;
       my $logpa = 0;
       my $logpc = 0;
       my $logpt = 0;
       my $logpg = 0;

       while ($recorre < scalar(@vector)){


	   if ($vector[$recorre] eq "A" || $vector[$recorre] eq "a") {
	       $pa = $pa + 1;
	   }

	   if( $vector[$recorre] eq "C" || $vector[$recorre] eq "c") {
	       $pc = $pc + 1;
	   }

	   if ($vector[$recorre] eq "T" || $vector[$recorre] eq "t") {
	       $pt = $pt + 1;
	   }
	   if ($vector[$recorre] eq "G" || $vector[$recorre] eq "g") {
	       $pg = $pg + 1;
	   }

	   $recorre = $recorre + 1;
       }
       

       $logpa = log($pa/scalar(@vector));
       $logpc = log($pc/scalar(@vector)); 
       $logpt = log($pt/scalar(@vector));
       $logpg = log($pg/scalar(@vector));



       close (SEQPROM);

       ### càlcul de piM a partir de la hash actual--> obtenim hash de log de freqüències relatives
      
       my %frequencia = ("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 $h =0;
       

       while ($h < $i){              #recorrem taula hash
	 
	
	   if ($motiu{"A"}[$h] == 0){
	       $frequencia{"A"}[$h]= -999;
	      }
	   else {
	      $frequencia{"A"}[$h] = log ($motiu{"A"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
	   }
	   if ($motiu{"C"}[$h] == 0){
	       $frequencia{"C"}[$h]= -999;
	      }
	   else{
	       $frequencia{"C"}[$h] = log ($motiu{"C"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
	    }
	   if($motiu{"T"}[$h] == 0){
	       $frequencia{"T"}[$h]= -999;
           }
	   else{
	      $frequencia{"T"}[$h] = log ($motiu{"T"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
	    }
	   if($motiu{"G"}[$h] == 0){
	       $frequencia{"G"}[$h]= -999;
           }
	   else{
	       $frequencia{"G"}[$h] = log ($motiu{"G"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
	   }

	    $h = $h +1;
	 

       }

       #### resta de log, construcció de la matriu de 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 $w =0;
       

       while ($w < $i){              #recorrem taula hash
	 
	
	   $pesos{"A"}[$w]=  $frequencia{"A"}[$w]- $logpa;
	   $pesos{"C"}[$w]=  $frequencia{"C"}[$w]- $logpc;
	   $pesos{"T"}[$w]=  $frequencia{"T"}[$w]- $logpt;
	   $pesos{"G"}[$w]=  $frequencia{"G"}[$w]- $logpg;
	   
	   $w = $w +1;


       }

       ##  printar la taula hash    
      
       my @k = ("A", "C", "G", "T");
       my $p = 0;

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

	   my $j = 0;

	   print "$k[$p] ";

	   while ($j < $i ) {

	       print "\t$pesos{$k[$p]}[$j]";

	       $j = $j + 1;

	   }
    
       print "\n";

       $p = $p + 1; 
       }
       
       ### anàlisi seq promotora
       

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

       my $y = <PROMOTER>;
       
      
       @v=split(//,$y);


       my $pv = 0;
       my @van;
       my $millor = -999999; 
       my $posicio = 0;

       while ($pv < (scalar(@v)-$i+1)) { 
	    	

	   my $pva = 0;
	   my $voltes= 0;
	   

	   while ($voltes < $i ){   ## omplim vector d'anàlisi, serà tan llarg com la longitud del motiu ($i)

	       $van[$pva] = $v[$pv];

	       $pv = $pv +1;
	       $pva = $pva +1;
	       $voltes = $voltes +1;
	   } 
	 
	 
    
	   ## càlcul d'score del vector anàlisi

	  
	   $score = 0; 
	   $pva = 0;
  
	   while ($pva < scalar(@van)) {
	      
	      
               if ($van[$pva] eq 'A' || $van[$pva] eq 'a' ){
    
		   $score = $score + $pesos{"A"}[$pva];

	       }
	       if ($van[$pva] eq 'C' || $van[$pva] eq 'c'){
		   
		   $score = $score + $pesos{"C"}[$pva];

	       }
	
	       if ($van[$pva] eq 'T'|| $van[$pva] eq 't'){
    
		   $score = $score + $pesos{"T"}[$pva];

	       }

	       if ($van[$pva] eq 'G' || $van[$pva] eq 'g'){
    
		   $score = $score + $pesos{"G"}[$pva];

		   
	       }
	      
	       $pva = $pva +1 ;
	
	   }
	   
  
           if ($score > $millor){

	       $millor = $score;
	       $posicio = $pv - $i;
	   }
	
	   $pv = $pv - $i +1;

      }

       print "posició: $posicio\n";
       print "score: $millor\n";
     
      
      close (PROMOTER);


###### PART 4: estimar la freqüència d'observar a l'atzar aquest score:

##Algoritme per agafar una seqüència random:

       my $finsa100 = 0;
       my $c = 0;
       
       while ($finsa100 < 100){

	   my $n = scalar(@v);
	   my $ii = $n - 1;
      
	   while ($ii >= 0) {
	       
	       my $jj = int(rand($ii+1));
        
	       if ($ii != $jj) {
           
		   my $tmp = $v[$ii];

		   $v[$ii] = $v[$jj];
		   $v[$jj] = $tmp;
	       }

	       $ii = $ii - 1;
	   }


##Realitzem el score de la seqüència random actual (es farà per les 100 vegades!):	   
	   

	   my $pvrandom = 0;
	   my @vanrandom;
	   my $millorrandom = -999999;
	   	

	   while ($pvrandom < scalar(@v)-$i+1) { 
	      
	       my $pvarandom=0;
	       my $voltesrandom= 0;
	       

	       while ($voltesrandom < $i ){   ## omplim vector d'analisi, serà tan llarg com la 						longitud del motiu ($i)
		   
		   $vanrandom[$pvarandom] = $v[$pvrandom];

		   $pvrandom = $pvrandom +1;
		   $pvarandom = $pvarandom +1;
		   $voltesrandom = $voltesrandom +1;
	       } 
	      
	 
    
	   ## càlcul d'score del vector anàlisi
	      
	 
	       my $scorerandom = 0; 
	       $pvarandom = 0;
	       
	       	   
	       while ($pvarandom < scalar(@vanrandom)) {
		        
		   if ($vanrandom[$pvarandom] eq 'A' || $vanrandom[$pvarandom] eq 'a' ){
		       $scorerandom = $scorerandom + $pesos{"A"}[$pvarandom];
		   }
		   if ($vanrandom[$pvarandom] eq 'C' || $vanrandom[$pvarandom] eq 'c'){
		       $scorerandom = $scorerandom + $pesos{"C"}[$pvarandom];
		   }
		   if ($vanrandom[$pvarandom] eq 'T'|| $vanrandom[$pvarandom] eq 't'){
    		       $scorerandom = $scorerandom + $pesos{"T"}[$pvarandom];
		   }
		   if ($vanrandom[$pvarandom] eq 'G' || $vanrandom[$pvarandom] eq 'g'){
    		       $scorerandom = $scorerandom + $pesos{"G"}[$pvarandom];
		   }
 		   $pvarandom = $pvarandom +1 ;
	       }
	       
	
	       if ($scorerandom > $millorrandom){
     
		   $millorrandom = $scorerandom;
		            
	       }
	   
	      $pvrandom = $pvrandom - $i +1;

	  
	   }

	      
	   if ($millorrandom > $millor){
		           
		      $c = $c  + 1;
	   } 

	   $finsa100 = $finsa100 + 1;
       }

  
      
       my $pvalue = $c / 100;

       print "p-value $pvalue\n";
       print "\n";







      
       # deixar la hash  a 0.

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


       $i = 0;         #comencem a omplir la nova hash


   }


}


close (FITXERMEU);
       

