#!/usr/bin/perl -w

use strict;

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

my $linea=0;
my $sumaA=0;
my $sumaC=0;
my $sumaG=0;
my $sumaT=0;
my $largo=0;
my $sumatot=0;
my $scoremax;
my $propA=0;
my $propC=0;
my $propG=0;
my $propT=0;

while (<SEQPROM>){
    chomp($_);
    my @v=split(//,$_);
    my $longitud=scalar(@v);
    my $i=0;
    if ($linea>0) {
	while ($i<$longitud) {
	    if ($v[$i] eq "A" || $v[$i] eq "a") {
		$sumaA=$sumaA+1;
	    }
	    if ($v[$i] eq "C" || $v[$i] eq "c") {
		$sumaC=$sumaC+1;
	    }

	    if ($v[$i] eq "G" || $v[$i] eq "g") {
		$sumaG=$sumaG+1;

	    }

	    if ($v[$i] eq "T" || $v[$i] eq "t") {
		$sumaT=$sumaT+1;

	    }


	    $i=$i+1;
	}

    }
    $linea=$linea+1;
}

$sumatot = $sumaA + $sumaC + $sumaG + $sumaT;


$propA = $sumaA/$sumatot;
$propC = $sumaC/$sumatot;
$propG = $sumaG/$sumatot;
$propT = $sumaT/$sumatot;

close (SEQPROM);



my %motivo = ("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 $i=0;
my $cuenta=0;


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

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 $nom="";
while (<FH>){
    chomp($_);

    if (m/(FA)\s(.+)/) {
	$nom=$2;
	next ;

    }

    elsif ($_=~ m/(\d+)\t(\d+)\s+(\d+)\s+(\d+)\s+(\d+)$/){
	$motivo{"A"}[$i]=$2;
	$motivo{"C"}[$i]=$3;
	$motivo{"G"}[$i]=$4;
	$motivo{"T"}[$i]=$5;
	$largo = $largo+1;
	$i=$i+1;
	$cuenta=$cuenta+1;
	next;
    }



    if (($_=~m/(\/\/)\Z/) && ($i>1)) {

	my @k = keys(%motivo);
	my $lugar = 0;


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

	while ($posicion<$cuenta){
	    my $j=0;
	    my $sumaA=0;
	    my $sumaC=0;
	    my $sumaG=0;
	    my $sumaT=0;
	    $sumaA = 
$motivo{"A"}[$posicion]/($motivo{"A"}[$posicion]+$motivo{"C"}[$posicion]+$motivo{"G"}[$posicion]+$motivo{"T"}[$posicion]);
	    $sumaC = 
$motivo{"C"}[$posicion]/($motivo{"A"}[$posicion]+$motivo{"C"}[$posicion]+$motivo{"G"}[$posicion]+$motivo{"T"}[$posicion]);
	    $sumaG = 
$motivo{"G"}[$posicion]/($motivo{"A"}[$posicion]+$motivo{"C"}[$posicion]+$motivo{"G"}[$posicion]+$motivo{"T"}[$posicion]);
	    $sumaT = 
$motivo{"T"}[$posicion]/($motivo{"A"}[$posicion]+$motivo{"C"}[$posicion]+$motivo{"G"}[$posicion]+$motivo{"T"}[$posicion]);

	    if($motivo{"A"}[$posicion]==0){

		$pesos{"A"}[$posicion]=-999;

	    }
	    elsif (($motivo{"A"}[$posicion])>0) {

		$pesos{"A"}[$posicion] = log($sumaA)-log $propA;
	    }

	    if($motivo{"C"}[$posicion]==0){


		$pesos{"C"}[$posicion]=-999;
	    }
	    elsif (($motivo{"C"}[$posicion])>0){

		$pesos{"C"}[$posicion] = log($sumaC)-log$propC;
	    }

	    if($motivo{"G"}[$posicion]==0){
		$pesos{"G"}[$posicion]=-999;
	    }
	    elsif (($motivo{"G"}[$posicion])>0) {

		$pesos{"G"}[$posicion] = log($sumaG)-log$propG;
	    }

	    if($motivo{"T"}[$posicion]==0){

		$pesos{"T"}[$posicion]=-999;
	    }
	    elsif (($motivo{"T"}[$posicion])>0){

		$pesos{"T"}[$posicion] = log($sumaT)-log$propT;
	    }


	    $posicion=$posicion+1;


	}
	@k = keys(%pesos);

	$lugar = 0;

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

	    while ($j < $cuenta) {
		$j = $j + 1;
	    }

	    $lugar = $lugar + 1;
	}



	my @v;
	my $secuencia;



	open (SEQPROM,"<$ARGV[0]");
	while (<SEQPROM>){
	    if(/>/){}
	    else{
		chomp($_);
		$secuencia = $secuencia . $_;
	    }
	}

	@v = split(//,$secuencia);
	my $totalnt=scalar(@v);

	my $inicial=0;
	my $ventana=0;
	my $scoreposicion=0;
	my $scoref=0;
	my $posiciomax=0;
	my $padromax;
	$scoremax=-9999999999;
	while ($inicial<$totalnt-$largo) {
	    $ventana=0;
	    $scoreposicion=0;
	    my $padro;
	    while ($ventana<$largo) {

		if ($v[$inicial+$ventana] eq "A" || $v[$inicial+$ventana] eq "a") {
		    $scoreposicion = $scoreposicion + $pesos{"A"}[$ventana];
		}
		if ($v[$inicial+$ventana] eq "C" || $v[$inicial+$ventana] eq "c") {
		    $scoreposicion = $scoreposicion + $pesos{"C"}[$ventana];
		}
		if ($v[$inicial+$ventana] eq "G" || $v[$inicial+$ventana] eq "g") {
		    $scoreposicion = $scoreposicion + $pesos{"G"}[$ventana];
		}
		if ($v[$inicial+$ventana] eq "T" || $v[$inicial+$ventana] eq "t") {
		    $scoreposicion = $scoreposicion + $pesos{"T"}[$ventana];
		}
		$ventana=$ventana+1;
		$padro = $padro . $v[$inicial+$ventana];
	    }
	    if ($scoreposicion>$scoremax) {
		$posiciomax= $inicial+2;
		$scoremax=$scoreposicion;
		$padromax = $padro;
	    }
	    $inicial = $inicial +1;
	}
	my $vez=0;
	my $padro;
	my $contar=0;

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

		    $secuenciarandom[$ii] = $secuenciarandom[$jj];
		    $secuenciarandom[$jj] = $tmp;
		}
		$ii = $ii - 1;
	    }
	    my $scorerand=-9999999;
	    $inicial=0;
	   my $pos = 0;
	    while ($inicial<$totalnt-$largo) {
		$ventana=0;
		$scoreposicion=0;

		while ($ventana<$largo) {

		    if ($secuenciarandom[$inicial+$ventana] eq "A" || 
$secuenciarandom[$inicial+$ventana] eq "a") {
			$scoreposicion = $scoreposicion + $pesos{"A"}[$ventana];
		    }
		    if ($secuenciarandom[$inicial+$ventana] eq "C" || 
$secuenciarandom[$inicial+$ventana] eq "c") {
			$scoreposicion = $scoreposicion + $pesos{"C"}[$ventana];
		    }

		    if ($secuenciarandom[$inicial+$ventana] eq "G" || 
$secuenciarandom[$inicial+$ventana] eq "g") {
			$scoreposicion = $scoreposicion + $pesos{"G"}[$ventana];
		    }
		    if ($secuenciarandom[$inicial+$ventana] eq "T" || 
$secuenciarandom[$inicial+$ventana] eq "t") {
			$scoreposicion = $scoreposicion + $pesos{"T"}[$ventana];
		    }
		    $ventana=$ventana+1;
		}


		if ($scoreposicion>$scorerand) {

		    $scorerand=$scoreposicion;
		    $pos=$inicial;

		}
		$inicial=$inicial +1;
	    }

	    if ($scorerand>=$scoremax) {

		$contar=$contar+1;
	    }

	    $vez = $vez+1;
	}


	my $c=$contar/100;
	print "The maximum score for $nom matrix was $scoremax (pvalue <= $c), at 
position $posiciomax [$padromax]\n";


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

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

	$i=0;
	$largo = 0;
	$cuenta=0;

    }

    close (SEQPROM);
}
close (FH);