#!/usr/bin/perl -w

use strict;

########### Vamos a abrir la secuencia mediante filehandle que va a asignar a un vector el exon 1 de la secuencia, a otro el 2 y a otro el 3. De esta manera obtenemos tres vectores donde cada uno de ellos contiene los exones 1,2 y 3 de todas las secuencias respectivamente.



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

my @secuencias = <MIFICHERO>;

close(MIFICHERO);
 

my $exon1;
my $exon2;
my $exon3;

my $contador = 0;

my (@exons1, @exons2, @exons3);


while ($contador < scalar@secuencias) {
   
    $secuencias[$contador] =~ m/^.*?([ATCGN]+)\s([ATCGN]+)\s([ATCGN]+)/;
   
    
    $exon1 = $1;
    $exon2 = $2;
    $exon3 = $3;

    push @exons1, $exon1;
    push @exons2, $exon2;
    push @exons3, $exon3;
    
    $contador = $contador + 1;


}


##########  Cogemos el vector que contiene los exones 1, leemos cada elemento del vector (1 exon) que dividimos en ventanas (substr) y lo guardamos y contamos medienta un hash (motivos). La ventana corre de uno en uno y solamente cogemos aquellos motivos en los que no aparece N.

my $contadorel1 = 0;
my $contadorfr1 = 0;
my $longitud = 5;
my $ventana;
my %motivos;


while ($contadorel1 < scalar@exons1) {

   
    $contadorfr1 = 0;

    while ($contadorfr1 <= length($exons1[$contadorel1]) - $longitud) {

	       
	$ventana = substr($exons1[$contadorel1],$contadorfr1,$longitud);

	
	if ($ventana =~ /^([ACGTacgt]+)$/) {

	    if (!exists($motivos{$ventana})) {
		$motivos{$ventana} = 0;
	    }
    

	    $motivos{$ventana} = $motivos{$ventana} + 1;
	}

	$contadorfr1 = $contadorfr1 + 1;

    }

   
    $contadorel1 = $contadorel1 + 1;

}

      


my $k =  keys(%motivos);

#print keys(%motivos),"\n";

#print "$k\n";

#######  Esto nos permite imprimir el hash para comprobar que el código anterior funcionaba correctamente.

my @k = keys(%motivos);
my @valor = values(%motivos);
my $i = 0;

while ($i < keys(%motivos)) {

  #print "($k[$i],$valor[$i])\n";

  $i = $i + 1;

}


##### Hacemos lo mismo para los exones alternativos (@exons2)

my $contadorel2 = 0;
my $contadorfr2 = 0;
my $longitud2 = 5;
my $ventana2;
my %motivos2;


while ($contadorel2 < scalar@exons2) {

  
    $contadorfr2 = 0;

    while ($contadorfr2 <= length($exons2[$contadorel2]) - $longitud2) {

       
	$ventana2 = substr($exons2[$contadorel2],$contadorfr2,$longitud2);

	if ($ventana2 =~ /^([ACGTacgt]+)$/) {

		if (!exists($motivos2{$ventana2})) {
		    $motivos2{$ventana2} = 0;
		}

		$motivos2{$ventana2} = $motivos2{$ventana2} + 1;
	    }

	
	$contadorfr2 = $contadorfr2 + 1;

    }

   
    $contadorel2 = $contadorel2 + 1;

}
  

my @k2 = keys(%motivos2);
my @valor2 = values(%motivos2);
my $i2 = 0;

while ($i2 < keys(%motivos2)) {

    #print "($k2[$i2],$valor2[$i2])\n";

    $i2 = $i2 + 1;

}


############   En este apartado se realiza el mismo proceso que en los dos apartados anteriores pero con la diferencia que no creamos un nuevos hash sino que cada uno de los nuevos motivos encontrados se añaden al hash %motivos, para formar el conjunto de motivos de los exones constitutivos. 


my $contadorel3 = 0;
my $contadorfr3 = 0;
my $longitud3 = 5;
my $ventana3;



while ($contadorel3 < scalar@exons3) {

  
    $contadorfr3 = 0;

    while ($contadorfr3 <= length($exons3[$contadorel3]) - $longitud3) {

       
	$ventana3 = substr($exons3[$contadorel3],$contadorfr3,$longitud3);


	if ($ventana3 =~ /^([ACGTacgt])+$/) {

	if (!exists($motivos{$ventana3})) {
	    $motivos{$ventana3} = 0;
	}

	$motivos{$ventana3} = $motivos{$ventana3} + 1;
	
}


	
	$contadorfr3 = $contadorfr3 + 1;

    }

   
    $contadorel3 = $contadorel3 + 1;

}


    
@k = keys(%motivos);
@valor = values(%motivos);
$i = 0;

while ($i < keys(%motivos)) {

  #print "($k[$i],$valor[$i])\n";

  $i = $i + 1;

}




########### Esto nos da el número de elementos que tiene el hash ,es decir, el número de motivos diferentes que existen  en los exones uno y tres. Además nos dice el número de apariciones de cada motivo. 


#my $k =  keys(%motivos);

#print keys(%motivos),"\n";

#print "$k\n";


########### Esto nos da el número de elementos que tiene el hash, es decir, el número de motivos diferentes del exon dos. Además nos dice el número de apariciones de cada motivo. 

#my $k2 =  keys(%motivos2);

#print keys(%motivos2),"\n";

#print "$k2\n";



###########  A partir de aquí se calcularan los valores relativos del hash respecto al valor total de ventanas, para ello se debe sumar el número de veces que se repite cada motivo, es decir, las frecuencias absolutas para luego dividir cada frecuencia absoluta entre el valor total. Este primera parte es para exones constitutivos.

my $j = 0;
my $sumas = 0;

while ($j < keys(%motivos)) {

    $sumas = $sumas + $valor[$j]; 

    $j = $j + 1;

}


$j = 0;


while ($j < keys(%motivos)) {

    $valor[$j] = $valor[$j] / $sumas; 

    #print "($k[$j],$valor[$j])\n";

    $j = $j + 1;

}


########## Esta parte para los exones alternativos.

my $j2 = 0;
my $sumas2 = 0;

while ($j2 < keys(%motivos2)) {

     $sumas2 = $sumas2 + $valor2[$j2]; 

     $j2 = $j2 + 1;

}

#print "$sumas2\n";

$j2 = 0;

while ($j2 < keys(%motivos2)) {

    $valor2[$j2] = $valor2[$j2] / $sumas2; 

   # print "($k2[$j2],$valor2[$j2])\n"; 

    $j2 = $j2 + 1;

}



##################  En el siguiente apartado, se van a comparar los dos hashes (constitutivos y alternativos) para determinar aquellos motivos que son comunes. El value asignado a cada motivo, es la resta de las frecuencias relativas multiplicada por 100000 para no perder información de los decimales relevantes.

##################  También permite agrupar en un hash los motivos que aparecen exclusivamente en los exones constitutivos.

  
my %comunes;
my %nocomunes;
my $a = 0;
my $b = 0;
my $palabra;



my @keys = keys(%motivos);

foreach my $key (@keys)
{
    if ($motivos2{$key}) {
    
#	 $valor2[$j2] = $valor2[$j2] / $sumas2; 
	my $a = $motivos{$key} / $sumas;
	my $b = $motivos2{$key} / $sumas2;
	$comunes{$key} = ($a - $b) * 100000;

   
    }
else {
	$nocomunes{$key} = $motivos{$key} / $sumas;
    }
    
}

my @comun_keys = keys(%comunes);
my @noncomun = keys(%nocomunes);


map {print "$_ : $comunes{$_}\n" } @comun_keys;



#map {print "$comunes{$_}\n" } @comun_keys;


#map {print "$_ \n" } @comun_keys;


print " Motivos exclusivos de exones constitutivos\n";

map {print "$_ : $nocomunes{$_}\n" } @noncomun;

##########  Esta parte permite separar los motivos que aparecen exclusivamnete en exones alternativos.

my %comunes2;
my %nocomunes2;
my $a2 = 0;
my $b2 = 0;



my @keys2 = keys(%motivos2);

foreach my $key2 (@keys2)
{
    if ($motivos{$key2}) {
    
#	 $valor2[$j2] = $valor2[$j2] / $sumas2; 
        $a = $motivos{$key2} / $sumas;
	$b = $motivos2{$key2} / $sumas2;
#	$comunes{$key} = ($a - $b) * 100000;

   
    }
else {
	$nocomunes2{$key2} = $motivos2{$key2} / $sumas2;
    }
    
}

#my @comun_keys = keys(%comunes);
my @noncomun2 = keys(%nocomunes2);

#map {print "$_ : $comunes{$_}\n" } @comun_keys;



#map {print "$comunes{$_}\n" } @comun_keys;


#map {print "$_ \n" } @comun_keys;

print " Motivos exclusivos de exones alternativos\n";

map {print "$_ : $nocomunes2{$_}\n" } @noncomun2;
