#!/usr/bin/perl

use strict;

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

my $posicio=0;
my $comptar=0;
my $suma=0;
my $divisioA=0;
my $divisioC=0;
my $divisioG=0;
my $divisioT=0;

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

while (<FITXERMEU>) {
   
    if (m/FA/){
	print;
	my $nom= $_; 
    }

    elsif ($_=~m/(\d+)\t(\d+)\s+(\d+)\s+(\d+)\s+(\d+)$/){
	   
	    $motiu{"A"}[$posicio]=$2;
	    $motiu{"C"}[$posicio]=$3;
	    $motiu{"G"}[$posicio]=$4;
	    $motiu{"T"}[$posicio]=$5;

	    $posicio=$posicio+1;
	    $comptar=$comptar+1;
	    
    }

    ####### Printem la matriu de cada factor de transcripcio #######
  
    if (($_=~m/(\/\/)\Z/) && ($posicio>1)) { 
     
	my @k = ("A","C","G","T");
	my $i = 0;

	while ($i < scalar(@k)) {
	    my $j = 0;

	    print $k[$i];
	    while ($j < $comptar) { 
		print "\t$motiu{$k[$i]}[$j]";
		$j = $j + 1;
	    }
	    print "\n";
	    $i = $i + 1;
     
	}

    ####### Deixem de printar #######
 	
	######## Comença l'arxiu adjunt per tal de veure la proporció de cada nucleòtid en la seqüència ########

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

	my $A=0;
	my $C=0;
	my $T=0;
	my $G=0;
	my @v;
	my $x=0;
	my $l=0;

	while (<FITXERMEU2>) {
	    chomp ($_);
	     
	    if (m/\>/){
	    }

	    else {
		@v= split (//, $_);
		$l= $l + scalar(@v);

	    while ($x<scalar(@v)){
	
		if ($v[$x]eq 'A' || $v[$x] eq 'a' ){

		    $A=$A+1;
		}
		if ($v[$x]eq 'C' || $v[$x] eq 'c'){

		    $C=$C+1;
		}
		if ($v[$x]eq 'G' || $v[$x] eq 'g'){

		    $G=$G+1;
		}
		if ($v[$x]eq 'T' || $v[$x] eq 't'){

		    $T=$T+1;
		}
		$x=$x+1;
	    }
	    }
	    
	    $x=0;
	  
  
	}

	$A= $A/$l;
	$C= $C/$l;
	$T= $T/$l;
	$G= $G/$l;





	close (FITXERMEU2);

        ####### Tanquem el fitxer de la sequencia del promotor #######

	####### Calcul de la matriu de pesos #######
	
	$posicio=0;
	
	while ($posicio<$comptar){
	

	    $suma= $motiu{"A"}[$posicio]+$motiu{"C"}[$posicio]+$motiu{"G"}[$posicio]+$motiu{"T"}[$posicio];
	
	    $divisioA= $motiu{"A"}[$posicio]/$suma;
	    $divisioC= $motiu{"C"}[$posicio]/$suma;
	    $divisioG= $motiu{"G"}[$posicio]/$suma;
	    $divisioT= $motiu{"T"}[$posicio]/$suma;

	    if ($divisioA==0){

		$motiu{"A"}[$posicio]= -999; 
	    }
	    else{

		$motiu{"A"}[$posicio]=(log ($divisioA)/ log (10)) - (log ($A)/ log (10));  
	    }

	    if ($divisioC==0){

		$motiu{"C"}[$posicio]= -999;
	    }
	    else{
		$motiu{"C"}[$posicio]= (log ($divisioC)/ log (10)) -(log ($C)/ log (10));
	    }
	    if ($divisioG==0){
		
		$motiu{"G"}[$posicio]= -999; 
	    }
	    else {
		$motiu{"G"}[$posicio]= (log ($divisioG)/log (10)) - (log ($G)/log (10)); 
	    }
	    if ($divisioT==0){

		$motiu{"T"}[$posicio]= -999; 	
	    }

	    else { 
		$motiu{"T"}[$posicio]= (log ($divisioT)/log (10))- (log ($T)/ log (10)); 
	    }

	    $posicio= $posicio+1;
	    $suma=0;
        }

	####### Printem la matriu de pesos #######

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

	while ($i < scalar(@k)) {
	    my $j = 0;

	    print $k[$i];
	    while ($j < $comptar) { 
		print "\t$motiu{$k[$i]}[$j]";
		$j = $j + 1;
	    }
	    print "\n";
	    $i = $i + 1;
	}
 
	####### Deixem de printar #######
       
	open (FITXERMEU2,"<$ARGV[1]");
	    my $seq;

	while (<FITXERMEU2>) {
	    chomp ($_);
	    

	    if (m/\>/){
	    }

	    else {
	    $seq= $seq.$_;
	    @v= split (//, $seq);
	    $l= scalar(@v);
	    
	    }
	}
	    
	    my $nt=1;
	    my $nt2=0;
	    my $valor=0;
	    $posicio=0;
	    my $score=-1000000;
	    my $inici=0;
	    my $final=0;
	    my $compara=0;

	    while ($nt <= (($l- $comptar) + 1)) {      

		$compara=0;

		while ($compara < $comptar){
		
		    if ($v[$nt2] eq 'A' || $v[$nt2] eq 'a' ){
			
			$valor= $valor + $motiu{"A"}[$posicio];
		    }
		    
		    if ($v[$nt2]eq 'C' || $v[$nt2] eq 'c'){

			$valor= $valor + $motiu{"C"}[$posicio];
		    }
		
		    if ($v[$nt2]eq 'G' || $v[$nt2] eq 'g'){

			$valor= $valor + $motiu{"G"}[$posicio];
		    }
	
		    if ($v[$nt2]eq 'T' || $v[$nt2] eq 't'){

			$valor= $valor + $motiu{"T"}[$posicio];
		    }
		
		    $posicio= $posicio+1;
		    $nt2=$nt2+1;
		    $compara= $compara +1; 
		   
		}

		if ($valor > $score){

		    $score = $valor;
		    $inici = $nt; 
		    $final = $inici + ($comptar-1);
		   
		}
	
		$nt= $nt +1; 
		$nt2= $nt - 1;
		$posicio=0;
		$valor=0;
	

	    }
	 
	    print "La puntuacio maxima obtinguda es $score\n";
	    print "La posicio d'inici es $inici i la posicio final es $final\n";
	 
	close (FITXERMEU2);

        ####### Programa per permutar la seqüència i càlcul del p value associat#######

	my $cops = 0;
	my $c=0;
	
	open (FITXERMEU2,"<$ARGV[1]");
	    my $seq2;

	while ($cops < 100){
	  
	    
	    while (<FITXERMEU2>) {
	    chomp ($_);
	    

	    if (m/\>/){
	    }

	    else {
	    $seq2= $seq2.$_;
	    @v= split (//, $seq2);
	    $l= scalar(@v);
	    
	}
	}


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

		    $v[$i] = $v[$j];
		    $v[$j] = $tmp;
		}
		$i = $i - 1;
	    }
	
	    my $nt1=0;
	    my $nt22=0;
	    my $valor2=0;
	    my $posicio2=0;
	    my $compara2=0;
	  
	
	    while ($nt1 <= (($l- $comptar) +1)) {      

		$compara2=0;

		while ($compara2 < $comptar){
		
		    if ($v[$nt22] eq 'A' || $v[$nt22] eq 'a' ){
			
			$valor2= $valor2 + $motiu{"A"}[$posicio2];
		    }
		    
		    if ($v[$nt22]eq 'C' || $v[$nt22] eq 'c'){

			$valor2= $valor2 + $motiu{"C"}[$posicio2];
		    }
		
		    if ($v[$nt22]eq 'G' || $v[$nt22] eq 'g'){

			$valor2= $valor2 + $motiu{"G"}[$posicio2];
		    }
	
		    if ($v[$nt22]eq 'T' || $v[$nt22] eq 't'){

			$valor2= $valor2 + $motiu{"T"}[$posicio2];
		    }
		
		    $posicio2= $posicio2+1;
		    $nt22=$nt22+1;
		    $compara2= $compara2 +1; 
		   
		}
	

		if ($valor2 >= $score){
		    $nt1 = $l- $comptar +2;
		    $c= $c+1;
		   
		}
	
		$nt1= $nt1 +1; 
		$nt22= $nt1;
		$posicio2=0;
		$valor2=0;
	
	

	    }
	 $cops= $cops+1;
	}
	
	my $pvalue=$c/100;
        print "El p value associat es $pvalue\n";


	close (FITXERMEU2);

	####### Tornem a començar amb el seguent FT #########

	%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]);
	
	$comptar=0;
	$posicio=0;


    }
  
}  

close (FITXERMEU);