#!/usr/bin/perl -w

use strict;


######################### Relacion de especies y sus respectivos codigos numericos#######################

my %especies=("Arabidopsis thaliana"=>"3702",
	      "Bos taurus"=>"9913",
	      "Caenorhabditis elegans"=>"6239",
	      "Danio rerio"=>"7955",
	      "Drosophila melanogaster"=>"7227",
	      "Homo sapiens"=>"9606",
	      "Hordeum vulgare"=>"4513",
	      "Lycopersicon esculentum"=>"4081",
	      "Medicago truncatula"=>"3880",
	      "Mus musculus"=>"10090",
	      "Oryza sativa"=>"4530",
	      "Rattus norvegicus"=>"10116",
	      "Sus scrofa"=>"9823",
	      "Triticum aestivum"=>"4565",
	      "Xenopus laevis"=>"8355",
	      "Zea mays"=>"4577");
#########################################################################################################



#### fingerprint.pl ####

#### Este programa permite convertir los datos de una base de datos tipo UCSC(tipo RefGene), en un ####
#### archivo que contenga el EXON FINGERPRINT de cada uno de los exones para cada gen.             ####
######################################################################################################


### Primer fingerprint ###

my $refGene1=$ARGV[0];

my $Especie1=$ARGV[1];                 #variable que guardara temporalmente el nombre de la 1a especie
my $especie1=$especies{$Especie1};  #variable que almacena el nombre de la especie

open (REFGENE,"< $refGene1");

open (STDOUT,"> fingerprint1.txt");

while (<REFGENE>){
    my $id;       #variable con los identificadores
    my $nex;      #variable con el numero de exones
    my $i=0;      #contador de lineas
    my @input;    #vector con todas las posiciones separadas por linea

    @input= split (/\s+/, $_);
    
    $id=$input[0];
    
    print $id, "\t";       #imprimimos ID del gen que esta leyendo

    my $y;                 #contador de exones
    my @principio;         #vector que recoge las coordenadas de INICIO de los exones
    my @final;             #vector que recoge las coordenadas de FINAL de los exones
 
    
#### Almacenamos las posiciones de inicio y final de exon en dos vectores ####

    @principio= ($input[8]=~ m /(\d+)/g);
    @final= ($input[9]=~ m /(\d+)/g);    


#### Contamos el numero real de exones codificantes ####

    $nex=$input[7];          #Numero real de exones.
    my $enc1=0;              #numero de exones no codificantes desde 5'.
    my $enc2=0;              #numero de exones no codificantes desde 3'.

    #extremo5'
    $y=0;
    while ($y<$input[7]){
	
	if (($principio[$y]<=$input[5])&&($input[5]<=$final[$y])){    #Si el principio del CDS está 
                                                                      #entre el principio y final de exon.
	    
	    $y=$input[7];                                             #para salir del  bucle.
	}
	else {
	    $enc1=$enc1+1;
	    $y=$y+1;
	}
    }

    #extremo3'
    $y=$input[7]-1;
    while ($y>$enc1){

	if (($principio[$y]<=$input[6])&&($input[6]<=$final[$y])){   #Si el principio del CDS está 
                                                                     #entre el principio y final de exon

	    $y=$enc1;                                                #para salir del bucle.
	}
	else {
	    
	    $enc2=$enc2+1;
	    $y=$y-1;
	}
    }
    
    $nex=$nex-$enc1-$enc2;                                          #Número real de exones codificantes.
    
    print $nex, "\t";

    my $tamanyo;      #tamanyo del exon
    my $tammod;       #tamaño real de cada exon al quitar lo que no codifica
    my $fase5=3;      #fase 5'
    my $fase3;        #fase 3'
    my $a=0;          #Es 3-fase3,lo que (nucleotidos que le faltan al ultimo codon del anterior exon)
                      # le restamos al tamanyo del siguiente exon antes de especificar su fase.

   
    $y=$enc1;               #Empezamos a calcular la fase a partir del primer exon que forma el CDS.
    while ($y<=($nex+$enc1-1)){
	$tamanyo=$final[$y]-$principio[$y];
	
	
	if ($nex==1){                      #si solo hay un exon CODIFICANTE.
	    $tammod=$input[6]-$input[5];
	}
	elsif ($y==$enc1){
	    $tammod=$final[$y]-$input[5];   #$input[5] especifica el primer nucleotido 
                                             #codificante.(donde empieza el CDS)
	}
	elsif ($y==($nex+$enc1-1)){
	    $tammod=$input[6]-$principio[$y];  #$input[6] especifica el ultimo nucleotido 
                                                #codificante.(donde acaba el CDS)
	}
	else {
	    $tammod=$tamanyo;
	}

	#modificamos el tamanyo para no tener en cuenta los nucleotidos del anterior exon.

	$tammod=$tammod-$a;

	#calculamos la fase.
	
	if ($tammod% 3==0){
	    
	    $fase3= 0;
	}
	elsif ($tammod% 3==1){
	    
	    $fase3= 1;
	}
	elsif ($tammod% 3==2){
	    
	    $fase3= 2;
	}
	if ($y==($nex+$enc1-1)){
	    $fase3= 3;
	}

	#print el fingerprint de cada exon;
	print $fase5, ":", $fase3, ":", ($tammod+$a), "\t";

	$fase5= $fase3;       # la fase 5' de un exon es la misma que la 3' del anterior.
	$y=$y+1;              # pasamos al exon siguiente.
	$a= 3-$fase3;         # preparamos el modificador de tamaño del exon $a.
		
    }
    
    print "\n";               #cambio de linea (porque cambiamos de gen);

}

close (REFGENE);


############################################################################################################################

### Segundo fingerprint ###

my $refGene2=$ARGV[2];

open (REFGENE,"< $refGene2");

open (STDOUT,"> fingerprint2.txt");

while (<REFGENE>){
    my $id;       #variable con los identificadores
    my $nex;      #variable con el numero de exones
    my $i=0;      #contador de lineas
    my @input;    #vector con todas las posiciones separadas por linea

    @input= split (/\s+/, $_);
    
    $id=$input[0];
    
    print $id, "\t";       #imprimimos ID del gen que esta leyendo

    my $y;                 #contador de exones
    my @principio;         #vector que recoge las coordenadas de INICIO de los exones
    my @final;             #vector que recoge las coordenadas de FINAL de los exones
 
    
#### Almacenamos las posiciones de inicio y final de exon en dos vectores ####

    @principio= ($input[8]=~ m /(\d+)/g);
    @final= ($input[9]=~ m /(\d+)/g);    


#### Contamos el numero real de exones codificantes ####

    $nex=$input[7];          #Numero real de exones.
    my $enc1=0;              #numero de exones no codificantes desde 5'.
    my $enc2=0;              #numero de exones no codificantes desde 3'.

    #extremo5'
    $y=0;
    while ($y<$input[7]){
	
	if (($principio[$y]<=$input[5])&&($input[5]<=$final[$y])){    #Si el principio del CDS está 
                                                                      #entre el principio y final de exon.
	    
	    $y=$input[7];                                             #para salir del  bucle.
	}
	else {
	    $enc1=$enc1+1;
	    $y=$y+1;
	}
    }

    #extremo3'
    $y=$input[7]-1;
    while ($y>$enc1){

	if (($principio[$y]<=$input[6])&&($input[6]<=$final[$y])){   #Si el principio del CDS está 
                                                                     #entre el principio y final de exon

	    $y=$enc1;                                                #para salir del bucle.
	}
	else {
	    
	    $enc2=$enc2+1;
	    $y=$y-1;
	}
    }
    
    $nex=$nex-$enc1-$enc2;                                          #Número real de exones codificantes.
    
    print $nex, "\t";

    my $tamanyo;      #tamanyo del exon
    my $tammod;       #tamaño real de cada exon al quitar lo que no codifica
    my $fase5=3;      #fase 5'
    my $fase3;        #fase 3'
    my $a=0;          #Es 3-fase3,lo que (nucleotidos que le faltan al ultimo codon del anterior exon)
                      # le restamos al tamanyo del siguiente exon antes de especificar su fase.

   
    $y=$enc1;               #Empezamos a calcular la fase a partir del primer exon que forma el CDS.
    while ($y<=($nex+$enc1-1)){
	$tamanyo=$final[$y]-$principio[$y];
	
	
	if ($nex==1){                      #si solo hay un exon CODIFICANTE.
	    $tammod=$input[6]-$input[5];
	}
	elsif ($y==$enc1){
	    $tammod=$final[$y]-$input[5];   #$input[5] especifica el primer nucleotido 
                                             #codificante.(donde empieza el CDS)
	}
	elsif ($y==($nex+$enc1-1)){
	    $tammod=$input[6]-$principio[$y];  #$input[6] especifica el ultimo nucleotido 
                                                #codificante.(donde acaba el CDS)
	}
	else {
	    $tammod=$tamanyo;
	}

	#modificamos el tamanyo para no tener en cuenta los nucleotidos del anterior exon.

	$tammod=$tammod-$a;

	#calculamos la fase.
	
	if ($tammod% 3==0){
	    
	    $fase3= 0;
	}
	elsif ($tammod% 3==1){
	    
	    $fase3= 1;
	}
	elsif ($tammod% 3==2){
	    
	    $fase3= 2;
	}
	if ($y==($nex+$enc1-1)){
	    $fase3= 3;
	}

	#print el fingerprint de cada exon;
	print $fase5, ":", $fase3, ":", ($tammod+$a), "\t";

	$fase5= $fase3;       # la fase 5' de un exon es la misma que la 3' del anterior.
	$y=$y+1;              # pasamos al exon siguiente.
	$a= 3-$fase3;         # preparamos el modificador de tamaño del exon $a.
		
    }
    
    print "\n";               #cambio de linea (porque cambiamos de gen);

}

close (REFGENE);



####################################################################################################################
####################################################################################################################


#### ortofingerprint.pl ####

#### Este programa permite aparear los exon fingerprint de dos genes ortólogos, a partir de la  ####
#### lectura de los archivos donde se hayan almacenado los exon fingerprint de las dos especies ####
#### y el archivo con la relación de ortología entre estas dos especies.                        ####



#hash q guarde los fingerprint de la primera especie.

open (FINGER,"< fingerprint1.txt");


my %fesp1;            ##hash que contiene los exonfingerprint, usando el identificador como clave.
     
while (<FINGER>){
	
    m/(\w+)\s+([^\n]+)/;   ##expresion regular que captura en primera posicion el identificador 
	                        # y en segunda el numero de exones con el fingerprint

    $fesp1{$1}=$2;
}

close (FINGER);


#hash q guarde el fingerprint de la segunda especie

open (FINGER,"< fingerprint2.txt");

my %fesp2;            ##hash que contiene los exonfingerprint, usando el identificador como clave.
     
while (<FINGER>){
	
    m/(\w+)\s+([^\n]+)/;   ##expresion regular que captura en primera posicion el identificador 
	                        # y en segunda el numero de exones con el fingerprint

    $fesp2{$1}=$2;
}

close (FINGER);

#

#intruducir la lista de genes ortologos 

my $ortologos=$ARGV[3];


#nombre del archivo donde iran los fingerprint apareados

open (STDOUT, "> fingers_apareados.txt");


## Aparear los genes ortólogos.

open (ARCHIVO, "<$ortologos");     ##abre el archivo donde NOSOTROS hemos guardado la con la lista de ortologos

while (<ARCHIVO>){
    chomp $_;
    my @par= split (/\|/,$_);                         #cada gen sera un elemento del vector

    my $ID1;                                          #almacena el identificador del gen de la primera especie
    my $ID2;                                          #almacena el identificador del gen de la segunda especie
            

        #hay que tener en cuenta el orden de los genes en el archivo de ortologos

    if ($par[0] eq $especie1){                     
	$par[5]=~ m/(NM_\d+)\./;                           #busca mediante un expresión regular el identificador
	$ID1 = $1;                                          
    }
    else{
	$par[5]=~ m/(NM_\d+)\./;                           #busca mediante un expresión regular el identificador
	$ID2 = $1;
    }


    if ($par[1] eq $especie1){
	$par[8]=~ m/(NM_\d+)\./;                           #busca mediante un expresión regular el identificador
	$ID1 = $1;                                          
    }
    else{
	$par[8]=~ m/(NM_\d+)\./;                           #busca mediante un expresión regular el identificador
        $ID2 = $1;
    }


    if (($fesp1{$ID1})&&($fesp2{$ID2})){                #solo mostrara los fingerprint si se conoce el de los 
	                                                     #dos genes ortologos
	print "$ID1\t$fesp1{$ID1}\t|\t$ID2\t$fesp2{$ID2}\n";
    }
}

close (ARCHIVO);


####################################################################################################################
####################################################################################################################



#### alinear_exon_fingerprint.pl ####

#### Este programa le va a permitir alinear exon fingerprint de genes ortologos, proporcionando ####
#### una puntuación basada en el macht o mismatch de los elementos del triplete que conforma el ####
#### exon fingerprint, penalizando los gaps.                                                    ####
#### Este programa utiliza programación dinámica.

####################################################################################################


open (FINAL,"< fingers_apareados.txt");
open (STDOUT, ">salida.txt");

while (<FINAL>){
    
    my @score;    #objetivo ultimo final e inconmensurable. 
                  #Contiene el score del mejor alineamiento de los 2 fingerprints.

    # Para cada par de genes vamos a querer:
    # - Una variable singular que contenga el número de exones
    # - otra variable singular que guarde el ID.
    # - Un vector que contenga los fingerprints.

    
    my $nex1;     #contiene el número de exones del gen de la primera especie
    my $nex2;     #contiene el número de exones del gen de la segunda especie

    my $ID1;      #contiene el identificador del gen de la primera especie
    my $ID2;      #contiene el identificador del gen de la segunda especie

    my @finger1;  #fingerprint del gen de la primera especie
    my @finger2;  #fingerprint del gen de la segunda especie


    my @vectemp;  #vector temporal


    @vectemp=split (/\|/, $_);  #separamos los dos fingerprints de los genes ortologos,

     
    $vectemp[0]=~ m/(NM_\d+)\t(\d+)\t([\d\:\t]+)/;  # dividimos el PRIMER GEN según sus elementos 
                                                     # mediante expresiones regulares
    $ID1=$1;
    $nex1=$2;
    @finger1=split(/\t/,$3);
    

    $vectemp[1]=~ m/\t(NM_\d+)\t(\d+)\t([\d\:\t]+)/;  # dividimos el SEGUNDO GEN según sus elementos 
                                                     # mediante expresiones regulares
    $ID2=$1;
    $nex2=$2;
    @finger2=split(/\t/,$3);


#### ahora ya tenemos todos los datos de cada gen en las distintas variables


## dado q necesitaremos el tamaño real de los genes(sólo el CDS) para compararlos, vamos a calcularlo ahora  
## tomaremos el tamanyo de cada exon y los iremos sumando, hasta tener el tamanyo total del CDS.

    my $tam1;   #tamanyo del 1er gen.
    my $tam2;   #tamanyo del 2o gen.
    
    $tam1=tamanyo(\@finger1);
    $tam2=tamanyo(\@finger2);
    
   


## CONDICIONES INICIALES ##
    # estos dos genes pueden hallarse en una serie de situaciones distintas:
           # 1) mismo numero de exones, mismo tamanyo de CDS, los exones tienen tamanyos parecidos.
           # 2) mismo numero de exones, mismo tamanyo de CDS, diferente tamanyo de exones entre si.
           # 3) mismo numero de exones, diferente tamanyo del CDS.
           # 4) diferente numero de exones, mismo tamanyo del CDS.
           # 5) diferente numero de exones, diferente tamanyo del CDS.
     
    #definimos que el tamanyo de dos exones es parecido cuando difieren en menos de 6 nucleotidos (2 codones) por exon
    #para comparar el tamaño del CDS completo, multiplicamos por 6 el numero medio de exones entre los dos genes.
    #segun en que situacion se encuentre, recibira o no una penalizacion en el score final del alineamiento de los dos fingerprints.

    my $bonus;                #sera la penalizacion segun el caso.



    if ($nex1==$nex2) {
	if ((($tam1-$tam2)<=(6*($nex1+$nex2)/2))&&(($tam1-$tam2)>=(-6*($nex1+$nex2)/2))) {
	  
	   
	    my $difexon;                        #determina si los exones de cada fingerprint son parecidos en tamanyo o no.
	    
	    $difexon=diferenciaexones(\@finger1, \@finger2);
	    
############## CASO 1 ######################################################################################
	    
	    if ($difexon==0) {                  #si los exones tienen tamanyos parecidos (no son diferentes)
		                                #si se cumplen estas condiciones estamos en el caso 1).
		$bonus=1;
	    }

############## CASO 2 ######################################################################################

	    else {                              #si los exones tienen tamanyos diferentes

		$bonus=0.95;

	    }
	}
        
############## CASO 3 ######################################################################################
	
	else { 
	     $bonus=1-(abs($tam1-$tam2)/($tam1+$tam2));
	   
	}

    }


############## CASO 4 ######################################################################################    

    else {
    
	if ((($tam1-$tam2)<=(6*($nex1+$nex2)/2))&&(($tam1-$tam2)>=(-6*($nex1+$nex2)/2))) {
	   
	    $bonus=1-(abs($nex1-$nex2)/($nex1+$nex2));
	 	   
	}

############## CASO 5 ########################################################################################
	else { 
	     $bonus=1-((abs($nex1-$nex2)/($nex1+$nex2))*(abs($tam1-$tam2)/($tam1+$tam2)));
	 
	 }
    }
    
###################### ALINEAMIENTO DE FASES #################################################################

##### PASO 1 ####
    my @score1;                      #resultado de alinear los dos fingerprints sin realizarles ningun cambio.
    my @score2;                      #resultado de alinear los dos fingerprints intentando combinar exones.

    @score1=alinea(\@finger1,\@finger2);   #calculamos @score1.

######################


##### PASO 2 #########

    my @refinger=refingerprint(\@finger1,\@finger2);
    my @nuevfinger1=@{$refinger[0]};
    my @nuevfinger2=@{$refinger[1]};
    
    @score2=alinea(\@nuevfinger1,\@nuevfinger2);   #calculamos @score2.

######################

## SCORE + BONUS #####
### NORMALIZACION ####

    $score1[0]= (($score1[0]/((scalar(@score1)-1)* 2))*$bonus)* 100;
    $score2[0]= (($score2[0]/((scalar(@score2)-1)* 2))*$bonus* 0.9)* 100;

######################
    
##### SCORE OPTIMO:###
   
    if ($score1[0]>=$score2[0]) {
	@score=@score1;

    }
    else {
	@score=@score2;
    }

######################

### SALIDA DE DATOS ##
print "$ID1\t$ID2\t\t$score[0]\n";   


###############################################################################################################
###########################################  FUNCIONES  #######################################################
###############################################################################################################


###############################################################################################################

#NOMBRE: tamanyo
#OBJETIVO: calcular el tamanyo del CDS sumando los tamanyos de los exones codificantes.
#RETORNA: el tamanyo del CDS ($tamanyoCDS).

    sub tamanyo{
	my $i=0;                  #contador de exones.
	my @fingerprint=@{$_[0]}; #desreferenciamos el vector que contiene los fingerprints.
	my $tamanyoCDS=0;         #suma de los tamanyos de los exones (tamanyo del CDS).
	
	while ($i<scalar(@fingerprint)){
	    $fingerprint[$i]=~ m/\d\:\d\:(\d+)/;
	    $tamanyoCDS=$tamanyoCDS+$1;
	    $i=$i+1;
	}
	return $tamanyoCDS;
    }

################################################################################################################

###############################################################################################################

#NOMBRE: diferenciaexones
#OBJETIVO: determinar si entre si los exones de un par de genes con el mismo numero de exones difieren o no en tamanyo.
#RETORNA: SI(1) o NO(0), es decir, si los genes tienen o no algun exon de tamanyo significativamente diferente.($difex).

sub diferenciaexones{
    my $i=0;          #contador de exones.
    my @f1=@{$_[0]};  #primer fingerprint
    my @f2=@{$_[1]};  #segundo fingerprint
    my $difex;        #indica si los exones son diferentes o no.
    my $a;            #almacena la resta en valor absoluto
    my $t1;           #tamanyo exon primera secuencia
    my $t2;           #tamanyo exon segunda secuencia
    
    while ($i<scalar(@f1)) {

	$f1[$i]=~ m/\d\:\d\:(\d+)/;
	$t1=$1;
	$f2[$i]=~ m/\d\:\d\:(\d+)/;
	$t2=$1;
	
	if ($t1-$t2>=0) {
	    $a=$t1-$t2;
	}
	else {
	    $a=-1*($t1-$t2);
	}
	if ($a<=6) {
	    $difex=0;
	    $i=$i+1;
	}
	else{
	    $difex=$i;
	    $i=scalar(@f1);
	}
    }

    return $difex;
}
###############################################################################################################

###############################################################################################################

#NOMBRE: alinea
#OBJETIVO: calcula el score del mejor alineamiento posible entre dos fingerprints.
#RETORNA: retorna el score del mejor alineamiento ($score);

sub alinea{

    my $i=0;          #contador de exones.
    my @f1=@{$_[0]};  #primer fingerprint.
    my @f2=@{$_[1]};  #segundo fingerprint.
    my @fases1;       #fases 5':3' de todos los exones del primer gen
    my @fases2;       #fases 5':3' de todos los exones del segundo gen
    my @matpunt;      #matriz de puntuaciones, resultado de alinear los dos fingerprints.
    my @score;        #puntuacion del alineamiento optimo y los emparejamientos.
    my $penGAP=-1;    #penalizacion por alinear con un GAP.

    while ($i<scalar(@f1)) {
	$f1[$i]=~ m/(\d\:\d)\:\d+/;   #capturem les fases del 1er gen.
	$fases1[$i]=$1;
	$i=$i+1;
    }

    $i=0;
    while ($i<scalar(@f2)) {
	$f2[$i]=~ m/(\d\:\d)\:\d+/;   #capturem les fases del 2on gen.
	$fases2[$i]=$1;
	$i=$i+1;
    }

    my $n1=1;       #contador de posiciones en la matriz (filas), la aprovecharemos para contabilizar los exones.
    my $n2=1;       #contador de posiciones en la matriz (columnas), la aprovecharemos para contabilizar los exones.


    $matpunt[0][0]=0;                  #definimos la posicion (0,0) de la matriz

    while ($n2<=scalar(@f2)){           #definimos la primera fila de la matriz como -1.
	$matpunt[0][$n2]=$penGAP;
	$n2=$n2+1;
    }
    while ($n1<=scalar(@f1)){           #definimos la primera columna de la matriz como -1;
	$matpunt[$n1][0]=$penGAP;
	$n1=$n1+1;
    }

    #reinicializamos los contadores de exones.
    
    $n1=1; 
    $n2=1;


    my $a;
    my $b;
    my $c;

    ## para la recuperacion del alineamento.
   
    my @retroceso;     #matriz que contiene un codigo que define el mejor movimiento desde cualquier posicion.
    my @abcmax;        #indica el valor maxiomo de a,b y c y a cual corresponde.
    ##

 $retroceso[0][0]='x';                  #definimos la posicion (0,0) de la matriz de retroceso.

    while ($n2<=scalar(@f2)){           #definimos la primera fila de la matriz de retroceso como c.
	$retroceso[0][$n2]='c';
	$n2=$n2+1;
    }
    while ($n1<=scalar(@f1)){           #definimos la primera columna de la matriz como b;
	$retroceso[$n1][0]='b';
	$n1=$n1+1;
    }

    $n1=1; 
    $n2=1;


    
    while ($n1<=scalar(@f1)){

	$n2=1;                     #reinicializamos los valores de las columnas.


	while ($n2<=scalar(@f2)){
	    
	   
	    $a= $matpunt[$n1-1][$n2-1]+ puntuacion ($fases1[$n1-1],$fases2[$n2-1]);
                                            #$a es el resultado de alinear las fases de dos exones en las posiciones n1 y n2.
	    $b= ($matpunt[$n1-1][$n2])+$penGAP;   #$b es el resultado de alinear el exon en posicion n1 del primer gen con un gap.
	  
	    $c= ($matpunt[$n1][$n2-1])+$penGAP;   #$c es el resultado de alinear el exon en posicion n2 del segundo gen con un gap.



	    @abcmax=maxabc($a,$b,$c);         #calculamos el valor maximo de a,b y c.
	    $matpunt [$n1][$n2]=$abcmax[0];   #escogemos el alineamiento que obtiene la puntuacion maxima
	    $retroceso[$n1][$n2]=$abcmax[1];
	
	    
	    $n2=$n2+1;                    #avanzamos en las columnas.
	}

	$n1=$n1+1;                        #avanzamos en las filas.
    }
    
    $score[0]=$matpunt[$n1-1][$n2-1];        #tomamos el valor final del alineamiento como puntuacion maxima.

    my @temp= retroceso(\@retroceso,$n1, $n2,\@fases1, \@fases2);
    push(@score,@temp);

    return @score;
}




###############################################################################################################

###############################################################################################################

#NOMBRE: puntuacion
#OBJETIVO: especificar la puntuacion de alinear las fases de dos exones.
#RETORNA: el resultado de alinear estas fases.($result)

sub puntuacion{

    my $fase1=$_[0];  #fases del primer exon
    my $fase2=$_[1];  #fases del segundo exon

    my $result;       #puntuacion del alineamiento.
    
    $fase1=~ m/(\d)\:(\d)/;
    my $x=$1;         #fase 5' del primer exon
    my $y=$2;         #fase 3' del primer exon
    
    $fase2=~ m/(\d)\:(\d)/;
    my $w=$1;         #fase 5' del segundo exon
    my $z=$2;         #fase 3' del segundo exon

    if (($x==$w)&&($y==$z)){
	$result=2;
    }
    elsif (($x!=$w)&&($y!=$z)){
	$result=0;
    }
    else{
	$result=1;
    }
    return $result;
}

###############################################################################################################

###############################################################################################################

#NOMBRE: maxabc
#OBJETIVO: calcula el valor de $a,$b y $c, e indica cual es.
#RETORNA: un hash de una posicion indicando la letra con puntuacion maxima y que valor tiene.

sub maxabc{

    my $x=-99999999;   #variable que almacenara el valor maximo.
    my @abcmax;        #indica el valor maximo y en la segunda posicion, la letra a la que corresponde.

    if ($_[0]>$x){
	$x=$_[0];
	$abcmax[0]=$x;
	$abcmax[1]='a';
    }
    if ($_[1]>$x){
	$x=$_[1];
	$abcmax[0]=$x;
	$abcmax[1]='b';
    }
    if ($_[2]>$x){
	$x=$_[2];
	$abcmax[0]=$x;
	$abcmax[1]='c';
    }
    
    return @abcmax;
    
}


###############################################################################################################


###############################################################################################################

#NOMBRE: retroceso
#OBJETIVO: retornar los emparejamientos de exones.
#RETORNA: un vector con los exones emparejados con otros exones o con GAPs(@apareamiento)
#input:  matriz de punteros, longitud fingerprint1 + 2, long fingerprint2 +2, fases1, fases2
#uso:     my @temp= retroceso(\@retroceso,$n1, $n2,\@fases1, \@fases2);

sub retroceso{
    my @recorrido=@{$_[0]};             #matriz con los movimientos optimos.
    my @f1=@{$_[3]};                    #vector con las fases del 1er gen.
    my @f2=@{$_[4]};                    #vector con las fases del 2o gen.
    my $n1=$_[1]-1;                       #indica en que exon de 1 nos encontramos.
    my $n2=$_[2]-1;                       #indica en que exon de 2 nos encontramos.
    my @apareamiento;                   #contiene todos los apareamientos.
   

    
    while ($recorrido[$n1][$n2]ne 'x') {

	
        #movimiento en diagonal.
	if ($recorrido[$n1][$n2] eq 'a'){
	    unshift(@apareamiento,"$f1[$n1-1] | $f2[$n2-1]");

	    $n1=$n1-1;
	    $n2=$n2-1;

	 }
	
	#movimiento hacia la izquierda (alineamos n2 con un gap);
	if ($recorrido[$n1][$n2] eq 'c'){
	    unshift(@apareamiento,"**GAP** | $f2[$n2-1]");

	    $n2=$n2-1;    #corremos hacia la izquierda.

	 }

	#movimiento hacia arriba (alineamos n1 con un gap);
	if ($recorrido[$n1][$n2] eq 'b'){
	    unshift(@apareamiento,"$f1[$n1-1] | **GAP**");

	    $n1=$n1-1;    #corremos hacia arriba.

	   }
	
    }

    return @apareamiento;
}
###############################################################################################################

###############################################################################################################

#NOMBRE: refingerprint
#OBJETIVO: toma los fingerprint de dos genes,compara el tamaño de los exones, y si es posible los fusiona en 
         # un solo exon, cambiando el fingerprint del gen. Fusionara todos los exones posibles.
         # se consideraran tamaños distintos si es superior a 6.
#RETORNA: dos vectores, cada uno con el nuevo fingerprint y además el nuevo numero de exones.


sub refingerprint{
    my @finger1=@{$_[0]};
    my @finger2=@{$_[1]};
      
    my @f1;                  
    my @f2;
  
    my $i1=0;          #contador de exones de primer gen
    my $i2=0;          #contador de exones del segundo gen

    my @tamanyos1;     #vector q guardo los tamanyos de los exones originales del primer gen
    my @tamanyos2;     #vector q guardo los tamanyos de los exones originales del segundo gen

    my $tammod=0;
    my $tammod2=0;
    my $tammodb=0;
    my $tammod2b=0;
    my $refing;
    my $refingb;

    while ($i1<scalar(@finger1)) {
	$finger1[$i1]=~ m/\d\:\d\:(\d+)/;
	$tamanyos1[$i1]=$1;
	$i1=$i1+1;
    }
    while ($i2<scalar(@finger2)){
	$finger2[$i2]=~ m/\d\:\d\:(\d+)/;
	$tamanyos2[$i2]=$1;
	$i2=$i2+1;
    }

    $i1=0;
    $i2=0;
   
  GRAN_LOOP:
    while (($i1<scalar(@finger1))&&($i2<scalar(@finger2)))  {

	if((compara($tamanyos1[$i1],$tamanyos2[$i2]))==1) {    #los dos exones tienen el mismo tamaño
	    push(@f1,$finger1[$i1]);                           #add el fingerprint de ese exon al vector definitivo
	    push(@f2,$finger2[$i2]);                           #add el fingerprint de ese exon al vector definitivo

	    $i1=$i1+1;
	    $i2=$i2+1;
	}
	elsif(($tamanyos1[$i1]>$tamanyos2[$i2])&&($tamanyos2[$i2+1])) {                #los exones tienen distinto tamanyo
	   
	    $tammod=(($tamanyos2[$i2])+($tamanyos2[$i2+1])); 
	    
	    if(compara($tamanyos1[$i1],$tammod)==1) {            #el nuevo exon tiene tamanyo similar al q era +grande
		
		$refingb=refing($finger2[$i2],$finger2[$i2+1]); #calcula las nuevas fases e incluye el nuevo tamanyo
		$refing="$refingb:$tammod";
	
		push(@f1,$finger1[$i1]);
		push(@f2,$refing);

		$i1=$i1+1;
		$i2=$i2+2;                                     #tiene q avanzar dos posiciones en el vector original
	    }
	    elsif($tamanyos1[$i1]>$tammod && ($tamanyos2[$i2+2]) ){   #el nuevo tamanyo sigue siendo menor, tratamos de meter un tercer exon
		
		$tammod2=($tammod+$tamanyos2[$i2+2]);
		
		if($tamanyos1[$i1]+6<$tammod2) {                  #no entra un tercer exon
		    
		    $refingb=refing($finger2[$i2],$finger2[$i2+1]);
		    $refing="$refingb:$tammod";
		    
		    push(@f1,$finger1[$i1]);
		    push(@f2,$refing);
		    
		    $i1=$i1+1;
		    $i2=$i2+2;
                }
		else{
		    $refingb=refing($finger2[$i2],$finger2[$i2+2]);
		    $refing="$refingb:$tammod2";
		    
		    push(@f1,$finger1[$i1]);
		    push(@f2,$refing);

		    $i1=$i1+1;
		    $i2=$i2+3;
		}

	    }
	    else {                                              #el nuevo tamanyo es más grande asi q lo descartamos
		 push(@f1,$finger1[$i1]);                           #add el fingerprint de ese exon al vector definitivo
		 push(@f2,$finger2[$i2]);                           #add el fingerprint de ese exon al vector definitivo
	                                              
		 $i1=$i1+1;
		 $i2=$i2+1;

	     }
	}
	elsif($tamanyos1[$i1+1]) {
	    
	    $tammodb=(($tamanyos1[$i1])+($tamanyos1[$i1+1])); 
	    
	    if(compara($tamanyos2[$i2],$tammodb)==1) {            #el nuevo exon tiene tamanyo similar al q era +grande
		
		$refingb=refing($finger1[$i1],$finger1[$i1+1]); #calcula las nuevas fases e incluye el nuevo tamanyo
		$refing="$refingb:$tammodb";
	
		push(@f2,$finger2[$i2]);
		push(@f1,$refing);

		$i1=$i1+2;
		$i2=$i2+1;                                     #tiene q avanzar dos posiciones en el vector original
	    }

	    elsif( ($tamanyos2[$i2]>$tammodb) &&($tamanyos1[$i1+2]) ){                     #el nuevo tamanyo sigue siendo menor, tratamos de meter un tercer exon
	
		$tammod2b=($tammodb+$tamanyos1[$i1+2]);
		
		if($tamanyos2[$i2]+6<$tammod2b) {                  #no entra un tercer exon
		    
		    $refingb=refing($finger1[$i1],$finger1[$i1+1]);
		    $refing="$refingb:$tammodb";
		   		    
		    push(@f2,$finger2[$i2]);
		    push(@f1,$refing);

		    $i1=$i1+2;
		    $i2=$i2+1;
                }
		else{
		    $refingb=refing($finger1[$i1],$finger1[$i1+2]);
		    $refing="$refingb:$tammod2b";
		    
		    push(@f2,$finger2[$i2]);
		    push(@f1,$refing);

		    $i1=$i1+3;
		    $i2=$i2+1;
		}
		
	    }

	    else {                                              #el nuevo tamanyo es más grande asi q lo descartamos
		 push(@f1,$finger1[$i1]);                           #add el fingerprint de ese exon al vector definitivo
		 push(@f2,$finger2[$i2]);                           #add el fingerprint de ese exon al vector definitivo
	                                              
		 $i1=$i1+1;
		 $i2=$i2+1;

	     }
	
	 	    
	}
	else {                                              #el nuevo tamanyo es más grande asi q lo descartamos
	    push(@f1,$finger1[$i1]);                           #add el fingerprint de ese exon al vector definitivo
	    push(@f2,$finger2[$i2]);                           #add el fingerprint de ese exon al vector definitivo
	    
	    $i1=$i1+1;
	    $i2=$i2+1;
	}
	
    }
    
    
    if ($i1<scalar(@finger1)) {
	while ($i1<scalar(@finger1)) {
	    push(@f1,$finger1[$i1]);
	    $i1=$i1+1;
	}
    }

    if ($i2<scalar(@finger2)) {
	while ($i2<scalar(@finger2)) {
	    push(@f2,$finger2[$i2]);
	    $i2=$i2+1;
	}
    }

    
    return (\@f1,\@f2);
 
                                 
}


###############################################################################################################

#NOMBRE:compara
#OBJETIVO:dados dos numeros (correspondientes a los tamanyos de dos exones) compara si son iguales
          #(igual significa que no difieran en mas de unidades).
#RETORNA: son iguales ($x=1) o diferentes ($x=0)

sub compara{
    my $a=$_[0];
    my $b=$_[1];

    my $x;

    if (abs($a-$b)>6){
	$x=0;
    }
    else{
	$x=1;
    }
   
    return $x;
}


###############################################################################################################


###############################################################################################################


#NOMBRE:refing
#OBJETIVO: introduciendo los fingerprint de dos exones, capturar la nueva fase tras fusionar los exones
#RETORNA: la nueva fase ($x)


sub refing{

    my $exon1=$_[0];
    my $exon2=$_[1];
    my $cinco;                 #fase en 5'
    my $tres;                  #fase en 3'

    my $x;

    $exon1=~ m/(\d)\:\d\:\d+/;
    $cinco=$1;
    $exon2=~ m/\d\:(\d)\:\d+/;
    $tres=$1;

    $x="$cinco:$tres";

    return $x;
}

###############################################################################################################

###############################################################################################################

#NOMBRE: max
#OBJETIVO: calcula el valor maximo dados n valores.
#RETORNA: el valor maximo ($x)

sub max{

    my $x=-99999999;   #variable que almacenara el valor maximo.
    my $i=0;           #contador de posiciones.

    while ($i<scalar(@_)){

	if ($_[$i]>=$x){
	    $x=$_[$i];
	}
	$i=$i+1;
    }
    return $x;
}


###############################################################################################################
###############################################################################################################


###############################################################################################################
}
close (FINAL);

