#!/usr/bin/perl -w

use strict;

#####################################################################################################
###############################  PROCEDIMENT 0: INICIALITZAR VARIABLES   ############################
#####################################################################################################

my @seq1;              #vector per la sequencia 1
my @seq2;              #vector per la sequencia 2
my $ID1;               #Identificador per sequencia 1
my $ID2;               #Identificador per sequencia 2
my $file1;             #Nom del fitxer on es troba la primera sequencia
my $file2;             #Nom del fitxer on es troba la segona sequencia
my $file3;             #Nom del fitxer de la matriu de pesos
my %pesos;             #Hash per guardar la matriu de pesos
my @elements;          #Llista dels elements identificatius del hash
my ($e1,$e2,$e3,$e4);  #Parametres d'entrada a partir de funcions 
                         #(per varios retorns no individuals)
my @matriu;            #matriu d'alineament
my @matriuLEFT;        #matriu de valors des de l'esquerra a [X][Y][0]
                         #i moviments a [X][Y][1]
my @matriuUP;          #matriu de valors des de dalta [X][Y][0]
                         #i moviments a [X][Y][1]
my @alineats;          #sequencia de moviments durant alineament a partir
                         #de matriuLEFT i matriuUP en posicio [X][Y][1] mes codi U,D o L
my $Igap;              #valor d'iniciar un gap
my $Cgap;              #valor de continuar un gap
my $score;             #valor maxim de la matriu
my @seq1AL;            #sequencia 1 alineada
my @seq2AL;            #sequencia 2 alineada
my @seqAL;             #modul de simbols d'alineament
my $identity;          #numero d'Aa coincidnets en l'alineament
my $similarity;        #numero d'Aa amb valor superior a 0 en l'alineament
my $gaps;              #numero de gaps afegits a l'alineament

#####################################################################################################
################################   PROCEDIMENT 1: INTRODUIR DADES   #################################
#####################################################################################################

#######################################
###   ASSIGNAR PARAMETRES ENTRADA   ###
#######################################

$file1 = $ARGV[0];                       #cada entrada d'arxiu correspon a una sequencia
$file2 = $ARGV[1];
$file3 = $ARGV[2];                       #aquesta ultima correspon a la matriu
$Igap = $ARGV[3];                        #inicial gap penalty
$Cgap = $ARGV[4];                        #continue gap penalty

################################
###   INTRODUIR SEQUENCIES   ###
################################

($ID1,@seq1) = info_sequencia($file1);   #adquirim les sequencies amb la funcio
($ID2,@seq2) = info_sequencia($file2);

unshift (@seq1,"-");                     #afegim un "-" (gap) al principi de la sequencia
unshift (@seq2,"-");         

#####################################
###   INTRODUIR MATRIU DE PESOS   ###
#####################################

($e1,$e2) = obtenir_pesos($file3);      #adquirim la matriu amb la funcio

%pesos = %{$e1};                         #Restaurem els valors   
@elements = @{$e2};   

#####################################################################################################
#################################   PROCEDIMENT 2: OMPLIR MATRIU   ##################################
#####################################################################################################

($e1,$score,$e2) = omplir_matriu(\@seq1,\@seq2,\%pesos,$Igap,$Cgap);

@matriu = @{$e1};         #per la matriu, la seq2 representa les files
                          #per la matriu, la seq1 representa les columnes
@alineats = @{$e2};

#imprimir_matriu(\@alineats,11,11);
#imprimir_matriu(\@matriu,11,11);
#####################################################################################################
################################   PROCEDIMENT 3: RECUPERAR SEQUENCIES   ############################
#####################################################################################################

($e1,$e2,$e3,$identity,$similarity,$gaps) = recuperar_cami(\@seq1,\@seq2,\@alineats,\@matriu,\%pesos);

@seq1AL = @{$e1};      #Rescatem resultats
@seq2AL = @{$e2};
@seqAL = @{$e3};

#####################################################################################################
#################################   PROCEDIMENT 4: IMPRIMIR RESULTATS   #############################
#####################################################################################################

imprimir_resultats(\@seq1AL,\@seq2AL,\@seqAL,$ID1,$ID2,$score,$file3,$Igap,$Cgap,$identity,$similarity,$gaps);











#####################################################################################################
############################################              ###########################################
############################################   FUNCIONS   ###########################################
############################################              ###########################################
#####################################################################################################

#########
###FUNCIO ESCANEJAR_SEQUENCIA
###
###OBJECTIU: Escaneja la sequencia i la assigna al vector
###
###PARAMETRES: [0] --> arxiu on es troba la sequencia
###
###RETORNA: l'identificador i la sequencia
#########

sub info_sequencia{

    my $seq = "";          #sobre aquesta variable anirem unificant les linees de la sequencia (inicialitza amb conjunt buit)
    my @sequencia;         #aqui guardem la sequencia amb un element per posicio
    my $identificador;     #enmagatzema el nom de la sequencia
    my $arxiu = $_[0];

    if (!open (FILE, $arxiu)){                        #obrim l'arxiu sobre el que volem treballar    
	die "No s'ha pogut obrir l'arxiu $arxiu.\n";  #si no el trobes pararia el programa
    }
    else{
 
	$identificador = <FILE>;              #escanejem la primera linea, que fa referencia al nom
	chomp ($identificador);

	while (<FILE>){                       #per cada linea que trobem i fins la ultima...
	    chomp ($_);                       #...
	    $seq = $seq . $_;                 #les anem unificant sobre $seq
	}

	$seq = uc($seq);                      #les fiquem en majuscules per facilitar la comparacio

	@sequencia = split (//,$seq);         #dividim la sequencia per elements sobre un array

	close (FILE);
    }

    return ($identificador,@sequencia);
}

#########
###FUNCIO OBTENIR PESOS
###
###OBJECTIU: Escanejar la matriu de pesos
###
###PARAMETRES: [0] --> fitxer de la matriu
###
###RETORNA: La matriu de pesos i els identificadors. 
###         Ambdos casos desreferenciats 
#########

sub obtenir_pesos{

    my $arxiu = $_[0];        #arxiu d'entrada
    my %matriu;               #hash on guardarem l'info
    my @aminoacids;           #guarda el nom dels elements de la primera fila per etiquetar els valors
    my $posicio = 0;          #diu la posicio d'scan util on ens trobem dins la matriu
    my $trobat = 0;           #per scanejar la linea d'identificadors d'elements
    my $i;                    #contador
    my @valors;               #guarda per cada fila els valors d'aquesta

    if (!open (FILE, $arxiu)){                        #obrim l'arxiu sobre el que volem treballar    
	die "No s'ha pogut obrir l'arxiu $arxiu.\n";  #si no el trobes pararia el programa
    }
    else{
	while (<FILE>){
	    if (!(m/\#/)){                            #si no te la marca d'info adicional
		if (!$trobat){                        #el primer cop que pasem ho fem per aqui per agafar l'identificador
		    chomp ($_);
		    $_ = substr($_,1);                #hem d'eliminar el primer caracter que es un espai
		    $_ = uc($_);                      #UPER-CASE: ho passem a majuscules
		    @aminoacids = split (/\s+/,$_);   #ho repartim en un vector
		    $trobat++;
		}	
		else{                                 #per els demes casos ho considerem valors i ho asignem al hash
		    chomp ($_);
		    $_ = substr($_,1);
		    @valors = split (/\s+/,$_);
		    for ($i=0;$i<scalar(@aminoacids);$i++){     #per cada posicio avancem per tots els Aa
			$matriu{$aminoacids[$posicio]}{$aminoacids[$i]} = $valors[$i];
		    }
		    $posicio++;
		}
	    }
	}
	close (FILE);
    }
    return (\%matriu,\@aminoacids);;
}

#########
###FUNCIO OMPLIR MATRIU
###
###OBJECTIU: Plenar la matriu d'alineament
###
###PARAMETRES: [0] --> sequencia1 desreferenciada
###            [1] --> sequencia2 desreferenciada
###            [2] --> matriu de pesos desreferenciada
###            [3] --> inici gap penalty
###            [4] --> continuar gap penalty
###
###RETORNA: La matriu plena, l'score maxim i la sequencia de moviments. 
###         La matriu i la sequencia desreferenciades. 
#########

sub omplir_matriu{

    my @seq1 = @{$_[0]};
    my @seq2 = @{$_[1]};
    my %pesos = %{$_[2]};
    my $Igap = $_[3];
    my $Cgap = $_[4];
    my @matriuLEFT;                        #matriu de valors des de l'esquerra a [X][Y][0]
                                             #i moviments a [X][Y][1]
    my @matriuUP;                          #matriu de valors des de dalt a [X][Y][0]
                                             #i moviments a [X][Y][1]
    my ($i,$j);                            #$i = files ; $j = columnes
    my @matriu;
    my $score;                             #valor maxim d'alineament
    my $Agap;                              #ACTUAL GAP: mira des de quin gap venim
    my @alineats;                          #matriu de "cami anterior"


                                                                            #####
    $matriu[0][0] = 0;                                                         ##
    $matriuUP[0][0][0] = 0;                                                    ##
    $matriuLEFT[0][0][0] = 0;                                                  ##
    $matriuUP[0][0][1] = 0;                                                    ##
    $matriuLEFT[0][0][1] = 0;                                                  ##
    for ($i=1;$i<scalar(@seq2);$i++){      #Replenem de 0 la primera columna   ##
	$matriu[$i][0] = $matriu[$i-1][0] + $Cgap;                             ##
	$matriuUP[$i][0][0] = $matriuUP[$i-1][0][0] + $Cgap;                   ##
   	$matriuLEFT[$i][0][0] = 0;                                             ##
	$matriuLEFT[$i][0][1] = 0;                                             ## 
	$matriuUP[$i][0][1] = $i;                                              ##
	$alineats[$i][0] = "U$i";                                              ##
    }                                                                          ####inicialitzem la matriu
    for ($j=1;$j<scalar(@seq1);$j++){      #Replenem de 0 la primera fila      ####progresio nomes conta
	$matriu[0][$j] = $matriu[0][$j-1] + $Cgap;                             ####estensio de gap
	$matriuLEFT[0][$j][0] = $matriuLEFT[0][$j-1][0] +$Cgap;                  ##
	$matriuUP[0][$j][0] = 0;                                               ##
	$matriuUP[0][$j][1] = 0;                                               ##
	$matriuLEFT[0][$j][1] = $j;                                            ##
	$alineats[0][$j] = "L$j";                                              ##
    }                                                                          ##
    $alineats[0][0] = 'S'; #START                                              ##
    $alineats[scalar(@seq2)-1][scalar(@seq1)-1] = 'E'; #END                 #####

    for ($i=1;$i<scalar(@seq2);$i++){
	for ($j=1;$j<scalar(@seq1);$j++){
            
      	    #si ja hem acabat una de les dues sequencies el valor dels gaps seguents sera 0
	    #serveix per ajustar per darrere compensant ajust d'inici de sequencia

#	    ($Igap,$Cgap) = avaluar_gap($i,$j,scalar(@seq1)-1,scalar(@seq2)-1,$Igap,$Cgap); 

	    #avaluacio del valor optim

	    ($matriuLEFT[$i][$j][0],$matriuLEFT[$i][$j][1]) = millor_opcio ($i,$j-1,$Igap,$Cgap,
									    $matriu[$i][$j-1],$matriuLEFT[$i][$j-1][0],
									    $matriuLEFT[$i][$j-1][1]);
	    ($matriuUP[$i][$j][0],$matriuUP[$i][$j][1]) = millor_opcio ($i-1,$j,$Igap,$Cgap,
									$matriu[$i-1][$j],$matriuUP[$i-1][$j][0],
									$matriuUP[$i-1][$j][1]);

	    ($matriu[$i][$j],$Agap) = max (($matriuUP[$i][$j][0]),                                 #venim des de dalt
					   ($matriuLEFT[$i][$j][0]),                               #venim des de esquerra
					   ($matriu[$i-1][$j-1]+$pesos{$seq1[$j]}{$seq2[$i]}));    #aliniem

	    if (!$Agap){
		$alineats[$i][$j] = 'D';
	    }
	    else{
		if($Agap > 0){
		    $alineats[$i][$j] = "U$matriuUP[$i][$j][1]";
		}
		else{
		    $alineats[$i][$j] = "L$matriuLEFT[$i][$j][1]";
		}
	    }

#	    ($Igap,$Cgap) = restaurar_gap($_[3],$_[4]);
	}
    }
    $score = $matriu[scalar(@seq2)-1][scalar(@seq1)-1]; 
    
    return (\@matriu,$score,\@alineats);
}

#########
###FUNCIO MILLOR OPCIO
###
###OBJECTIU: Decideix si es millor iniciar gap o continuar-ne un
###
###PARAMETRES: [0] --> moviment per linees
###            [1] --> moviment per columnes
###            [2] --> inici gap
###            [3] --> continuacio gap
###            [4] --> anterior total
###            [5] --> anterior direccional
###            [6] --> valor de salt anterior
###
###RETORNA: Valor maxim. Grandaria del salt
#########

sub millor_opcio

{
    
    my $i = $_[0]; 
    my $j = $_[1]; 
    my $Igap = $_[2];
    my $Cgap = $_[3];
    my $valor1=$_[5]+$Cgap;
    my $valor2=$_[4]+$Igap;
    my $lastsalt = $_[6];

    if ($valor1>=$valor2) 
    {
	return ($valor1,$lastsalt+1);
    }
    else
    {
	return ($valor2,1);
    }

}

#########
###FUNCIO MAX
###
###OBJECTIU: Donats tres valors trobar el maxim
###
###PARAMETRES: [0] --> valor1    (venim de dalt)
###            [1] --> valor2    (venim d'esquerra)
###            [2] --> valor3    (aliniem)
###
###RETORNA: Valor maxim. Cami per on s'ha vingut
###            -1 --> venim des de l'esquerra (afegim gap a seq2)
###             0 --> venim en diagonal (aliniem)
###             1 --> venim des de dalt (afegim gap a seq1)
#########

sub max{

    my $valor1 = $_[0];
    my $valor2 = $_[1];
    my $valor3 = $_[2];

    if (($valor1!=$valor2)&&($valor1!=$valor3)&&($valor2!=$valor3)){   #busca el major de tres valors diferents
	if (($valor1 > $valor2)&&($valor1 > $valor3)){
	    return ($valor1,1);
	}
	elsif (($valor2 > $valor1)&&($valor2 > $valor3)){
	    return ($valor2,-1);
	}
	elsif (($valor3 > $valor1)&&($valor3 > $valor2)){
	    return ($valor3,0);
	}
    }
    else{                                     #busca el major de tres valors on dos o mes son iguals
	if ($valor1==$valor2){
	    if ($valor1 > $valor3){
		return ($valor1,1);
	    }
	    else{
		return ($valor3,0);
	    }
	}
	elsif ($valor1==$valor3){
	    if ($valor1 > $valor2){
		return ($valor3,0);           #com que 3 i 1 son iguals preferim tornar 3 per ser mes restrictius
	    }                                 #amb els gaps
	    else{
		return ($valor2,-1);
	    }
	}
	elsif ($valor2==$valor3){
	    if ($valor2 > $valor1){
		return ($valor3,0);           #com que 3 i 1 son iguals preferim tornar 3 per ser mes restrictius
	    }                                 #amb els gaps
	    else{
		return ($valor1,1);
	    }
	}
    }
}

#########
###FUNCIO AVALUAR GAP
###
###OBJECTIU: Decideix el valor del gap a partir dels pasos anteriors
###
###PARAMETRES: [0] --> posicio i
###            [1] --> posicio j
###            [2] --> longitud seq1
###            [3] --> longitud seq2
###            [4] --> inici gap penalty
###            [5] --> continua gap penalty
###
###RETORNA: El valor de gap ajustat pel final. 
#########

sub avaluar_gap{

    my $i = $_[0];
    my $j = $_[1];
    my $seq1 = $_[2];
    my $seq2 = $_[3];
    my $Igap = $_[4];
    my $Cgap = $_[5];

    if (($i == $seq2)||($j == $seq1)){
	return (0,0);
    }
    else{
	return ($Igap,$Cgap);
    }
}

#########
###FUNCIO RESTAURAR GAP
###
###OBJECTIU: Restaura els valors inicials del gap penalty
###
###PARAMETRES: [0] --> inici gap penalty
###            [1] --> continue gap penalty
###
###RETORNA: Els valors de gap. 
#########

sub restaurar_gap{

    my $Igap = $_[0];
    my $Cgap = $_[1];

    return ($Igap,$Cgap);
}

#########
###FUNCIO RECUPERAR CAMI
###
###OBJECTIU: A partir d'una matriu omplerta recuperem l'alineament optim
###
###PARAMETRES: [0] --> sequencia1 desreferenciada
###            [1] --> sequencia2 desreferenciada
###            [2] --> matriu de cami desreferenciada
###            [3] --> matriu de valors d'alineament
###            [4] --> hash de pesos desreferenciat
###
###RETORNA: tres arrays corresponents a les tres matrius i els valors de identitat, similaritat i gaps
###         Els tres arrays desreferenciats. 
#########

sub recuperar_cami{

    my @seq1 = @{$_[0]};
    my @seq2 = @{$_[1]};
    my @cami = @{$_[2]};
    my @matriu = @{$_[3]};
    my %pesos = %{$_[4]};
    my @seq1AL;
    my @seq2AL;
    my @seqAL;
    my $i = scalar(@seq2)-1;
    my $j = scalar(@seq1)-1;
    my $k;
    my $identity=0;
    my $similarity=0;
    my $gaps=0;

    while ($cami[$i][$j] ne 'S'){                    #ens movem per la matriu
	if ($cami[$i][$j] eq 'D'){                   #ens movem diagonal (0)--> aparellament
	    unshift (@seq1AL,$seq1[$j]);       #Aparellar implica asignar la posicio de cada sequencia 
	    unshift (@seq2AL,$seq2[$i]);       #a la sequencia modificada per alinear
	    if ($seq1[$j] eq $seq2[$i]){
		unshift (@seqAL,"*");          #si es tracta del mateix Aa marcarem *
		$identity++;
	    }
	    else{
		if ($pesos{$seq1[$j]}{$seq2[$i]} > 0){
		    unshift(@seqAL,":");       #si els dos Aa es substitueixen mes que per atzar (>0) usem :
		    $similarity++;
		}
		else{
		    unshift(@seqAL,".");       #en qualsevol altre cas no posem res
		}
	    }
	    $i--;
	    $j--;
	}
	else{
	    if ($cami[$i][$j] =~ m/U(\d+)/){             #ens movem cap dalt (1) --> gap a seq1
		for ($k=0;$k<$1;$k++){
		    unshift (@seq1AL,"-");          #asignem gap a seq1 i afegim Aa de seq2
		    unshift (@seq2AL,$seq2[$i]);
		    unshift(@seqAL," ");            #desaparellament no mostra signe
		    $gaps++;
		    $i--;
		}
	    }
	    elsif ($cami[$i][$j] =~ m/L(\d+)/){              #ens movem cap a esquerra (-1) --> gap a seq2
		for ($k=0;$k<$1;$k++){
		    unshift (@seq1AL,$seq1[$j]);    #asignem gap a seq2 i afegim AA de seq1
		    unshift (@seq2AL,"-");
		    unshift(@seqAL," ");            #desaparellament no mostra signe
		    $gaps++;
		    $j--;
		}
	    }
	}
    }
    return (\@seq1AL,\@seq2AL,\@seqAL,$identity,$similarity,$gaps);
}

#########
###FUNCIO IMPRIMIR RESULTATS
###
###OBJECTIU: Imprimir l'alineament optim d'una manera entenedora i visualment atractiva
###
###PARAMETRES: [0] --> sequencia1 desreferenciada
###            [1] --> sequencia2 desreferenciada
###            [2] --> sequencia de simbols desreferenciada
###            [3] --> identificador seq1
###            [4] --> identificador seq2
###            [5] --> score maxim de l'alineament
###            [6] --> nom de la matriu de substitucio
###            [7] --> gap penalty
###            [8] --> extended penalty
###            [9] --> identitat
###            [10]--> similaritat
###            [11]--> numero de gaps
###
###RETORNA: No retorna, imprimeix resultats
#########

sub imprimir_resultats{

    my @seq1 = @{$_[0]};
    my @seq2 = @{$_[1]};
    my @seqAL = @{$_[2]};
    my $ID1 = $_[3];
    my $ID2 = $_[4];
    my $score = $_[5];
    my $pesos = $_[6];
    my $Igap = $_[7];
    my $Cgap = $_[8];
    my $identity = $_[9];
    my $similarity = $_[10];
    my $gaps = $_[11];
    my @pesos = split (//,$pesos);
    my $i=0;
    my $j=0;
    my $k=60;
    my $Lseq1=scalar(@seq1);
    my $IDpercent=$identity/$Lseq1 * 100;
    my $SIpercent=$similarity/$Lseq1 * 100;
    my $GApercent=$gaps/$Lseq1 * 100;

    print "#SEQUENCE ALIGNMENT v 0.1\n#\n";
    print "#SEQUENCE 1 : $ID1\n#SEQUENCE 2 : $ID2\n#\n";
    print "#SCORING MATRIX : ";
	for ($i=0;$pesos[$i] ne ".";$i++){
	    print "$pesos[$i]";
	}
    print "\n#\n";
    print "#Gap_penalty : $Igap\n#Extended_penalty : $Cgap\n#\n";
    print "#Identity : $identity\/$Lseq1\t";
    printf "(%.2f", $IDpercent;
    print "\%)\n";
    print "#Similarity : $similarity\/$Lseq1\t";
    printf "(%.2f",$SIpercent;
    print "\%)\n";
    print "#Gaps : $gaps\/$Lseq1\t";
    printf "(%.2f",$GApercent;
    print "\%)\n";
    print "#\n#SCORE : $score\n";

    $ID1 = substr($ID1,0,10);
    $ID2 = substr($ID2,0,10);
    while ($j<scalar(@seq1)){
	print "$ID1\t";
	for ($i=$j;($i<$k)&&($i<scalar(@seq1));$i++){ 
	    print "$seq1[$i]";
	}
	print "\t$i";
        print "\n";
	print "         \t";
	for ($i=$j;($i<$k)&&($i<scalar(@seqAL));$i++){
	    print "$seqAL[$i]";
	}
	print "\n";
	print"$ID2\t";
	for ($i=$j;($i<$k)&&($i<scalar(@seq2));$i++){
	    print "$seq2[$i]";
	}
	print "\t$i";
	print "\n";
	$j+=60;
	$k+=60;

	print "\n\n";
    }
}