#!/usr/bin/perl -w

use strict;


#########################################################
#		Treball de Bioinformàtica		#
#							#
#	Arnau Biayna Hierro NIA:30739			#
#	Daniela Stephanie Rivera Vallvé NIA: 33064	#
#							#
#########################################################

use constant ITERACIONS    => 100;
#

print "Carregant";
my %humanFactorsHash = llegeixHumanFactors($ARGV[0]);
my @v = llegeixSeqPromotora($ARGV[1]);
my ($totalNucleotids,%numPromo) = calculaProporcioNuncleotidsSP(\@v);
my %matriuPesos = calculaMatriuDePesos(\%humanFactorsHash,\%numPromo);
my %scoreMAX = calculaScore($totalNucleotids,\%humanFactorsHash,\@v,\%numPromo,\%matriuPesos);
my %pValues = calculaPValues($totalNucleotids,\%humanFactorsHash,\@v,\%scoreMAX,\%matriuPesos);
imprimeixPValues(\%pValues);

sub llegeixHumanFactors{
#llegim el fitxer de matrius, matriu a matriu, 
#enregistrant-la en memoria com a un hash de hashos
	my %humanFactorsHash;
	my $nomMatriu;
	my $posMotiu;
	my $nucleotid;
	my $fitxerHumanFactors = shift;
	#print "$fitxerHumanFactors fitxerHumanFactors\n";
	open(HUMANFACTORS,"< ".$fitxerHumanFactors);
	my $linia = <HUMANFACTORS>;
	while (<HUMANFACTORS>) {
	
		if (m/FA[\s]+(.+)/) {
			$nomMatriu = $1;
			#print $nomMatriu," nomMatriu\n";
  		}
		if (m/(XX)/){};
		if (m/(\d+)\s(\d+)\s(\d+)\s(\d+)\s(\d+)/) {
			$posMotiu = $1;
			#print $posMotiu," posMotiu\n $2 $3 $4 $5\n";
			$humanFactorsHash{$nomMatriu}{$posMotiu}{"A"}=$2;
			$humanFactorsHash{$nomMatriu}{$posMotiu}{"C"}=$3;
			$humanFactorsHash{$nomMatriu}{$posMotiu}{"G"}=$4;
			$humanFactorsHash{$nomMatriu}{$posMotiu}{"T"}=$5;
  		}
	}
	close(HUMANFACTORS);
	#print $humanFactorsHash{"AR [T00040]"}{"00"}{"A"};
	#print $humanFactorsHash{"AR [T00040]"}{"03"}{"T"};
	#print $humanFactorsHash{"AP-1 [T00029]"}{"04"}{"A"};
	#print $humanFactorsHash{"NRSF [T06124]"}{"00"}{"A"};
	#print %humanFactorsHash;
	#passem per referencia.
	#imprimeixDades(\%humanFactorsHash);
	return %humanFactorsHash;
}

sub llegeixSeqPromotora{
	my $fitxerSeqPromotora = shift;
	my $linia;
	open(SEQUNCIAPROMOTORA,"< $fitxerSeqPromotora");
		$linia = <SEQUNCIAPROMOTORA>;
		my @v;
		my $string = "";
		while($linia=<SEQUNCIAPROMOTORA>)
		{
			$string =~ s/\s//;
			chomp($linia);
			$string = $string.$linia;;
		}
		#pasamos la secuancia a mayusculas
		@v = split(//,uc($string));

		#print "V @v\n";
	close(SEQUNCIAPROMOTORA);
	return @v;
}


sub calculaProporcioNuncleotidsSP{
#Calculem la proporció d'ocurrències del 
#nucleòtid x al llarg de tota la seqüència promotora
	my $param1 = shift(@_); 
	my @v = @$param1;

	#my $totalNucleotids = ${$param2};
	my %numPromo;
	$numPromo{"A"} = 0;
	$numPromo{"C"} = 0;
	$numPromo{"G"} = 0;
	$numPromo{"T"} = 0;

	my $l;
	my $i;

	$l=scalar(@v);
	#print "$l\n";
	$i=$l-1;
	while ($i >= 0) {
		$numPromo{$v[$i]} = $numPromo{$v[$i]} + 1;
		$i = $i - 1;
	}

	my $totalNucleotids = $numPromo{"A"} + $numPromo{"C"} + $numPromo{"G"} + $numPromo{"T"};
	$numPromo{"A"} = $numPromo{"A"} / $totalNucleotids;
	$numPromo{"C"} = $numPromo{"C"} / $totalNucleotids;
	$numPromo{"G"} = $numPromo{"G"} / $totalNucleotids;
	$numPromo{"T"} = $numPromo{"T"} / $totalNucleotids;

	#print $numPromo{"A"}," ",$numPromo{"C"}," ",$numPromo{"G"}," ",$numPromo{"T"};
	return ($totalNucleotids,%numPromo);
}

sub calculaMatriuDePesos{
#Calcula matriu de pesos wi(x) = log(piM(x)) - log(pB(x)).
	my ( $param1, $param2 ) = @_; 	

	my %humanFactorsHash = %$param1;
	my %numPromo = %$param2;	
	my %matriuPesos;

	my $nomMatriu;
	my $posMotiu;
	my $nucleotid;

	for $nomMatriu ( keys %humanFactorsHash ) {
		for $posMotiu ( keys %{ $humanFactorsHash{$nomMatriu} } ) {
			for $nucleotid ( keys %{ $humanFactorsHash{$nomMatriu}{$posMotiu} } ) {
			  $matriuPesos{$nomMatriu}{$posMotiu}{$nucleotid} = -999;
				#print "PM $humanFactorsHash{$nomMatriu}{$posMotiu}{$nucleotid}\n";
			  if ($humanFactorsHash{$nomMatriu}{$posMotiu}{$nucleotid} != 0)
			  {
				if($numPromo{$nucleotid} != 0)
				{
			  my $logPM = log($humanFactorsHash{$nomMatriu}{$posMotiu}{$nucleotid});
				my $logPB = log($numPromo{$nucleotid});
			  	$matriuPesos{$nomMatriu}{$posMotiu}{$nucleotid} = $logPM - $logPB;
			  	}
			  } 
			

			  #print "$matriuPesos{$nomMatriu}{$posMotiu}{$nucleotid} ";
        		}
			#print "\n";
		}
		#print "\n";
	}
	return %matriuPesos;
}

sub calculaScore{
	#Calculem el score
	my $totalNucleotids = shift;
	my ( $param1, $param2, $param3, $param4) = @_; 	
	
	
	my %humanFactorsHash = %$param1;
	my @v = @$param2;
	my %numPromo = %$param3;
	my %matriuPesos = %$param4;

	my $nomMatriu;
	my $posMotiu;
	my $score = 0;
	my %scoreMAX;
	my $f = 0;
	my $j;
	my $i;
	
	#Per a cada posible combinacio de nucleotids
	for $nomMatriu ( keys %humanFactorsHash ) {
		#f son el numero de posicions de nucleotids dintre de cada motiu
		my @k = keys %{ $humanFactorsHash{$nomMatriu}};
		$f = scalar(@k);
		#print "$f f\n";
		$scoreMAX{$nomMatriu}{"Score"} = -99999;
		$scoreMAX{$nomMatriu}{"Posicio"} = -1;
		#print "totalNucleotids: $totalNucleotids\n";
		for ($i = 0; $i<=$totalNucleotids-$f; $i++) {
			$posMotiu = 0;
			$score = 0;
			for ($j = $i; $j<=$i+$f-1; $j++){
				
				if ($posMotiu<10) {$posMotiu = "0".$posMotiu;}
				#print "$nomMatriu $posMotiu $v[$j]\n";
				$score = $score + $matriuPesos{$nomMatriu}{$posMotiu}{$v[$j]};
				#print "$score + $matriuPesos{$nomMatriu}{$posMotiu}{$v[$j]}\n";
				$posMotiu = $posMotiu + 1;
			}
			#si el score calculat supera el parcial fins 
			#ara substituim i guardem la seva posició
			if ($score>$scoreMAX{$nomMatriu}{"Score"}) 
			{
				$scoreMAX{$nomMatriu}{"Score"}=$score;
				$scoreMAX{$nomMatriu}{"Posicio"} = $i;
				#print "Score superat $score $i\n";	
			}
		}
		
	}
	print "AP-1 [T00029]:$scoreMAX{\"AP-1 [T00029]\"}{\"Score\"}\n";
	return %scoreMAX;
}

sub comparaScores{

	my ( $param1, $param2, $param3) = @_; 	

	my %c = %$param1;
	my %scoreMAX = %$param2;
	my %scoreMAXPermutacio = %$param3;

	my $nomMatriu;

	for $nomMatriu ( keys %scoreMAX ) {
		#comptem en quants d'aquests casos la puntuació màxima 
		#observada és més gran o igual que s* ($scoreMAX). 
		#Denotem aquest nombre de cops per c.
		if ($scoreMAXPermutacio{$nomMatriu}{"Score"}>=$scoreMAX{$nomMatriu}{"Score"}){
			$c{$nomMatriu} = $c{$nomMatriu} + 1;
		}	
	}
	return %c;
}


sub calculaPermutacioFisherYates{
#assumim tenim la sequencia promotora en un vector @v, el seguent
#algorisme permutara a l'atzar els elements de @v
        my $param1 = shift(@_);
	my @v = @$param1;

	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;
        }
	#print "\nvector\n\n--------------------------\n";
	#print @v,"\n";
	return @v;
}

sub calculaPValues{
	my @permutacio;
	my @scorePermutacio;
	my $totalNucleotids = shift;
	my ( $param1, $param2, $param3, $param4 ) = @_; 	

	my %humanFactorsHash = %$param1;
	my @v = @$param2;
	my %scoreMAX = %$param3;
	my %matriuPesos = %$param4;

	my %c;
	my %pValues;
	my $i;
	my $nomMatriu;
	
	#Inicialitzo el hash de cs
	for $nomMatriu ( keys %scoreMAX ) {
		$c{$nomMatriu} = 0;
	}

	#repetim aquests dos passos 100 cops.
	#hem definit una constant que es diu ITERACIONS inicialitzada a 100.
	for ($i = 0; $i<ITERACIONS; $i++) {
	  @permutacio = calculaPermutacioFisherYates(\@v);
	  my %scoreMAXPermutacio = calculaScore($totalNucleotids,\%humanFactorsHash,\@permutacio,\%numPromo,\%matriuPesos);
	  %c = comparaScores(\%c,\%scoreMAX,\%scoreMAXPermutacio);
	}

	#direm que el p-value associat al FT és igual o més petit que c/100.
	for $nomMatriu ( keys %c ) {
		$pValues{$nomMatriu} = $c{$nomMatriu}/ITERACIONS; 
		#print "$nomMatriu: $c{$nomMatriu} \n"
	}
	return %pValues;
}

sub imprimeixPValues{
	my $param1 = shift(@_);
	my %pValues = %$param1;

	shift;
	my $nomMatriu;
	print "Els resultats són: \n";
	print "Nom Motiu \t\tp-value \t Posició d'inici \t Score\n";
	

	foreach $nomMatriu (sort {$pValues{$a} cmp $pValues{$b}} sort {$pValues{$a} cmp $scoreMAX{$b}} keys %pValues)
	{
    		print "$nomMatriu: \t";
		if (length($nomMatriu) < 14 ) {print "\t";}
		print "$pValues{$nomMatriu} \t\t $scoreMAX{$nomMatriu}{\"Posicio\"} \t\t\t";
		print" $scoreMAX{$nomMatriu}{\"Score\"}\n";
	}
}
sub imprimeixDades{

	my $param1 = shift(@_);
	my %humanFactorsHash=%$param1;
	my $nomMatriu;
	my $posMotiu;
	my $nucleotid;

	for $nomMatriu ( keys %humanFactorsHash ) {
		print "$nomMatriu: \n";
		for $posMotiu ( keys %{ $humanFactorsHash{$nomMatriu} } ) {
			for $nucleotid ( keys %{ $humanFactorsHash{$nomMatriu}{$posMotiu} } ) {
				print "$humanFactorsHash{$nomMatriu}{$posMotiu}{$nucleotid} ";
        		}
			print "\n";
		}
		print "\n";
	}

}
