#!/user/bin/perl -w
use strict;


my $line_fugu;       # Llegim les línies del proteoma de Fugu
my $line_zebra;      # Llegim les línies del proteoma de Zebrafish
my $marca_protf = 0;    
my $marca_protz = 0;

open (RESULTATS,"> percentatges_vertebrats.txt");

while (<STDIN>){        # Llegim les parelles d'homòlegs

    /(\w+)\s+(\w+)/;

    my $idfugu = $1;    # Guarda l'identificador de la seqüència de Fugu
    my $idzebra = $2;   # Guarda l'identificador del seu homòleg en Zebrafish
    
   


######  OBTENIM LA SEQÜÈNCIA DE LA PROTEÏNA DE FUGU A ALINIAR  ######

    open(FUGU,"</disc8/genomes/T.rubripes/ensembl/Fugu_rubripes.FUGU4.feb.pep.fa") || die "no puc obrir el fitxer!!\n";                           # Obrim proteoma de fugu
   

    open(CLUSIN1,">tempfugu.fa");       # FileHandle per guardar la seqüència
    
    $line_fugu = <FUGU>;

    while($line_fugu){
	
	if(($marca_protf == 1) && ($line_fugu =~ /^>/)) {# Evitem que continue guardant la resta de seq
	    $marca_protf = 0;
	}

	if($line_fugu =~ /$idfugu/) {   # Si ens troba l'identificador al proteoma

	    $marca_protf = 1;
	}
	if($marca_protf == 1) {
	    print CLUSIN1 "$line_fugu"; # Guardar l'identificador al FH 
       	}

    $line_fugu=<FUGU>;

    }

    close(FUGU);
    close(CLUSIN1);
    

   
    
######  OBTENIM LA SEQÜÈNCIA DE LA PROTEÏNA DE ZEBRAFISH A ALINIAR  ######

    open(ZEBRA,"</disc8/genomes/D.rerio/ensembl/Danio_rerio.ZFISH5.feb.pep.fa") || die "no puc obrir el fitxer!!\n";                                 # Obrim proteoma de zebra
    
    open(CLUSIN2,">tempzebra.fa");         # FH per guardar la seqüència
 
    $line_zebra = <ZEBRA>;
 
   while($line_zebra) {
       if(($marca_protz == 1) && ($line_zebra =~ /^>/)) {  # Per evitar guardar, després de la nostra seq
	    $marca_protz = 0;                              # les seq posteriors del proteoma
	}

	if($line_zebra =~ /$idzebra/) {                    # Si troba l'identificador al proteoma
	    $marca_protz = 1;
	}
	if($marca_protz == 1) {
	    print CLUSIN2 "$line_zebra";                   # Guardar l'identificador al FH 
       	}

    $line_zebra=<ZEBRA>;

    }
    close(ZEBRA);
    close(CLUSIN2);
    



######  AGRUPEM LES DUES SEQÜÈNCIES A ALINIAR PER ENVIAR-HO A CLUSTALW  ######

    open CLUSIN1,"<tempfugu.fa";
    open CLUSIN2,"<tempzebra.fa";
    open(CLUSTAL, ">tempclustal.fa");    # Guardem les 2 seqüències seguides

    while(<CLUSIN1>){
	print CLUSTAL "$_";              # Primer guardem la seqüència de Fugu               
    }
    while(<CLUSIN2>){
	print CLUSTAL "$_";              # A continuació, tot seguit, la seqüència de Zebrafish
    }
    close(CLUSTAL);
    close(CLUSIN1);
    close(CLUSIN2);




######  CLUSTALW  ######

# Enviem l'arxiu creat a clustalw

`/disc8/bin/clustalw1.8/clustalw tempclustal.fa`;

# Es pot veure un exemple de fitxer de clustalw clicant a "fitxer" a la pàgina web,
# útil per entendre la  següent part del programa.


######  OBTENCIÓ LONGITUD SEQÜÈNCIA I SÍMBOLS DE L'ALINIAMENT  ######

my @seql;              # Guardem les línies de la seqüència d'aminoàcids a un vector
my $seq;               # Guardem la seqüència d'aa seguida
my @v;                 # Guardem cada aa en una posició
my $longitudvector;    # Comptem el total d'aa
my @estrl;             # Guardem les línies amb els símbols de l'aliniament
my $estr;              # Guardem la seqüència de símbols seguida


if (!open(LONG,"< tempclustal.aln")) {
  print "clustalw.pl: impossible obrir tempclustal.aln del id $idfugu \n";
  exit(1);
}
my $a = 0;    # Comptador
my $b = 0;    # Comptador

while (<LONG>){                   # Llegim resultat de fer Clustalw

    if(m/\ANEWSINFRUP(\w+)\s+([\w\-]+)/) { #només reconeix la seqüència de Fugu

	$seql[$a] = $2;           # Fem que en seql es guardi la seqüència d'aminoàcids, per línies
	$a = $a +1;
    } elsif (/\A([\*\.:\s]+)\Z/) {# Posem \davant símbols per treure'ls el seu significat a les expr. reg.
	$estrl[$b] = $1;          # Fem que en estrl es guardi la seqüència de símbols, per línies
	$b = $b +1; 
    }

}
close (LONG);

# Concatenem la seqüència d'aminoàcids

$seq = "";           # Perquè comenci amb paraula buida. Posem tota la seqüència d'aa seguida

my $i=0;             # Comptador 

while ($i < scalar(@seql)) {
    chomp($seql[$i]);            # Per treure espai final de línia
    $seq = $seq . $seql[$i];     # Concatenació de totes les línies
    $i = $i + 1;
}

@v = split(//,$seq);             
$longitudvector = scalar(@v);    # Total d'aminoàcids de la seqüència més llarga de l'aliniament


# Concatenem la seqüència de símbols de l'aliniament

$estr  = "";        # Perquè comenci amb paraula buida                

my $j =0;           # Comptador

while ($j < scalar(@estrl)) {
    chomp($estrl[$j]);
    $estr = $estr.$estrl[$j];    # Concatenació de totes les línies
    $j = $j +1;
}




######  COMPTAR ELS MATCHS ENTRE SEQÜÈNCIES  ######

my @w;       # Guardem cada símbol en una posició
my $k =0;    # Comptador. Recorre les posicions del vector
my $n;       # Número total de símbols
my $s =0;    # Comptador de les estrelles. Per calcular identitat
my $t =0;    # Comptador dels símbols de similaritat

@w = split(//,$estr);          # Posa cada símbol a una posició
$n = scalar(@w);              

while ($k < $n) {
    if($w[$k] eq "*") {
	$s = $s +1;            # Suma el número de * total a l'aliniament
    }
    

    if(($w[$k] eq "*") || ($w[$k] eq ".") || ($w[$k] eq ":")) {
	   $t = $t +1;         # Suma els símbols de similaritat
       }
    $k = $k +1;             
}



######  CÀLCUL DEL %IDENTITAT I %SIMILARITAT  ######

my $identitat;
my $similaritat;


$identitat = $s/$longitudvector;

$similaritat = $t/$longitudvector;


print "$idfugu\t$idzebra\t$identitat\t$similaritat\t\n"; #xivato

print RESULTATS "$idfugu\t$idzebra\t$identitat\t$similaritat\t\n";
}