#!/usr/bin/perl -w

use strict;

####################### PREDICCIO D'EXONS: PREDICTOR ####################

# AUTORES: Ana Moreno, Roser Zaurin
# DATA: 2 abril 2003
# CURS: 4art 
# ESTUDIS: Biologia Humana, Facultat de Ciencies de la Salut i de la Vida

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

# OBJECTIU: Prediccio dels possibles exons dins d'una sequencia d'ADN 
# donada en format FASTA, mostrant, en un document de sortida en format
# GFF, la posicio de donnors i acceptors sites que limiten l'exo; la 
# puntuacio tenint en compte el biaix codificant; i la pauta de lectura 
# utilitzada.  




####################### PRE-TRACTAMENT DE DADES #########################

# si el nombre d'arguments al programa es mes petit que 4
# el programa no es pot executar perque no te tots els fitxers
# que necessita

if (scalar(@ARGV) < 6) {
    print "predictor.pl sequencia matriu_d_donnors matriu_d_acceptors proporcionscodons llindar largada_min_exo\n";
    exit(1);
}


# assignem els noms del vector d'arguments a variables
# amb noms mes significatius

my $fitxerseq = $ARGV[0];
my $matriu_d = $ARGV[1];
my $matriu_a = $ARGV[2];
my $fitxercod = $ARGV[3];
my $llindar = $ARGV[4];
my $llargada_min_exo = $ARGV[5];

# obrim el fitxer especificat al quart argument, que es on hi
#  ha les proporcions de cada codo en regio codificant, passat al
# programa des del shell

if (!open(PROPORCIONSCODONS,"< $fitxercod")) {
    print "predictor.pl: impossible obrir $fitxercod\n";
    exit(1);
}


# llegim les proporcions de codons del fitxer i les enregistrem
# al hash %proporcions_codons

my %proporcions_codons;

while (<PROPORCIONSCODONS>) {

    # a la variable anonima tenim la linia llegida i l'encaixem amb
    # la seguent expressio regular

    m/(\w+) (\d+\.\d+)/;

    # l'expressio regular anterior encaixa amb dues paraules la primera
    # de les quals es alfabetica i la segona un numero amb decimals.
    # les dues paraules estan agrupades amb parentesis lo qual significa
    # que ens podem referir a la primera com $1 i a la segona com $2

    $proporcions_codons{$1}=$2;

}

close(PROPORCIONSCODONS);

# obrim el fitxer especificat al segon argument, on hi ha la 
# sequencia problema,  passat al programa des del shell

if (!open(SEQUENCIA,"< $fitxerseq")) {
    print "predictor.pl: impossible obrir $fitxerseq\n";
    exit(1);
}


# el fitxer amb la sequencia esta en format FASTA. aixo implica que te
# una primera linia amb el simbol '>' seguit d'un identificador, i despres
# segueix un seguit de linies amb differents trossos de la sequencia

my $iden = <SEQUENCIA>;  # llegim la primera linia del format FASTA ">XXXX"
my @seql = <SEQUENCIA>;  # llegim les linies de la sequencia
my $seq  = "";           # declarem una variable on anira a parar tota la
                         # sequencia seguida, i la inicialitzem amb la
                         # paraula buida

close(SEQUENCIA);

# en la variable $seq posarem la sequencia sencera. aixo ho farem enganxant
# els trossos de sequencia que tenim a cada posicio del vector @seql

my $i=0;

while ($i < scalar(@seql)) {
    chomp($seql[$i]);
    $seq = $seq . $seql[$i];
    $i = $i + 1;
}

my $seq0 = "\U$seq";

# definim la mateixa variable per les altres dues pautes de lectura

my $seq1=substr ($seq0, 1);
my $seq2=substr ($seq1, 1);


##################### FI PRETRACTAMENT DE DADES  ###############################

#cridem la funcio per calcular stops

my @matriu0 = &calcula_stops($seq0);
my @matriu1 = &calcula_stops($seq1);
my @matriu2 = &calcula_stops($seq2);



#cridem la funcio que ens llegeix les matrius de donors i acceptors 
#per poder executar les dues funcions seguents 

my %matriu_1= &llegeix_matriu_pesos($matriu_d); 
my %matriu_2= &llegeix_matriu_pesos($matriu_a);

#cridem una funcio que ens busca la posicio dels donors i els acceptors 
#i ho posem en dues matrius de dues dimensions i tres columnes -orf, 
#posicions donors/acceptors, puntuacio de cada donor/acceptor- 


my @matriu_puntuacio_donnors_0 = &calcula_donor_acceptor(\%matriu_1,\@matriu0,$llindar); 
my @matriu_puntuacio_donnors_1 = &calcula_donor_acceptor(\%matriu_1,\@matriu1,$llindar);
my @matriu_puntuacio_donnors_2 = &calcula_donor_acceptor(\%matriu_1,\@matriu2,$llindar);

my @matriu_puntuacio_acceptors_0 = &calcula_donor_acceptor(\%matriu_2,\@matriu0,$llindar);
my @matriu_puntuacio_acceptors_1 = &calcula_donor_acceptor(\%matriu_2,\@matriu1,$llindar);
my @matriu_puntuacio_acceptors_2 = &calcula_donor_acceptor(\%matriu_2,\@matriu2,$llindar);


#afegim una columna que sera el frame 

@matriu_puntuacio_acceptors_0 = &calcula_frame(\@matriu_puntuacio_acceptors_0,\@matriu0);
@matriu_puntuacio_acceptors_1 = &calcula_frame(\@matriu_puntuacio_acceptors_1,\@matriu1);
@matriu_puntuacio_acceptors_2 = &calcula_frame(\@matriu_puntuacio_acceptors_2,\@matriu2);


   

#cridem la funcio que ens faci totes les combinacions possibles en cada 
#orf entre donor i acceptor sites per aixi construir tots els exons 
#putatius 

my @donnor_acc_0 = &combi(\@matriu_puntuacio_donnors_0,
                          \@matriu_puntuacio_acceptors_0,\@matriu0);
my @donnor_acc_1 = &combi(\@matriu_puntuacio_donnors_1,
                          \@matriu_puntuacio_acceptors_1,\@matriu1);
my @donnor_acc_2 = &combi(\@matriu_puntuacio_donnors_2,
                          \@matriu_puntuacio_acceptors_2,\@matriu2);

#corretgim les posicions absolutes i relatives de les pautes 1 i 2 sumant
#un i dos nucleotids respectivament

my $da;
for($da=0; $da<scalar(@donnor_acc_1); $da++){
    $donnor_acc_1[$da][1]=$donnor_acc_1[$da][1]+1;
    $donnor_acc_1[$da][2]=$donnor_acc_1[$da][2]+1;
    $donnor_acc_1[$da][3]=$donnor_acc_1[$da][3]+1;
    $donnor_acc_1[$da][4]=$donnor_acc_1[$da][4]+1;
}

for($da=0; $da<scalar(@donnor_acc_2); $da++){
    $donnor_acc_2[$da][1]=$donnor_acc_2[$da][1]+2;
    $donnor_acc_2[$da][2]=$donnor_acc_2[$da][2]+2;
    $donnor_acc_2[$da][3]=$donnor_acc_2[$da][3]+2;
    $donnor_acc_2[$da][4]=$donnor_acc_2[$da][4]+2;
}

#cridem la funcio que ens sumara el biaix codificant per a cada exo putatiu
#que ha trobat per a cada pauta de lectura


my $is;
for ($is=0; $is < scalar(@donnor_acc_0) ; $is++) {
    my $punt_biaix_0 = &biaix_codificant (\%proporcions_codons,substr($donnor_acc_0[$is][5],                                                                                                                    $matriu_puntuacio_acceptors_0[$is][3]));
    $donnor_acc_0[$is][6] = $punt_biaix_0+$donnor_acc_0[$is][0];
    $donnor_acc_0[$is][7] = $matriu_puntuacio_acceptors_0[$is][3];
}



for ($is=0; $is < scalar(@donnor_acc_1) ; $is++) {
    my $punt_biaix_1 = &biaix_codificant (\%proporcions_codons,substr($donnor_acc_1[$is][5],
                                           $matriu_puntuacio_acceptors_1[$is][3]));
    $donnor_acc_1[$is][6] = $punt_biaix_1+$donnor_acc_1[$is][0];
    $donnor_acc_1[$is][7] = $matriu_puntuacio_acceptors_1[$is][3];
}


for ($is=0; $is < scalar(@donnor_acc_2) ; $is++) {
    my $punt_biaix_2 = &biaix_codificant (\%proporcions_codons,substr($donnor_acc_2[$is][5],
                                           $matriu_puntuacio_acceptors_2[$is][3]));
    $donnor_acc_2[$is][6] = $punt_biaix_2+$donnor_acc_2[$is][0];
    $donnor_acc_2[$is][7] = $matriu_puntuacio_acceptors_2[$is][3];
}



#cridem la funcio que busca ATG

my @initials0 = &busca_ATG(\@matriu0);
my @initials1 = &busca_ATG(\@matriu1);
my @initials2 = &busca_ATG(\@matriu2);


#aprofitem la funcio combi per aparellar els ATG amb els donnors

my @exo_inicial0 = &combi(\@matriu_puntuacio_donnors_0,\@initials0,\@matriu0);
my @exo_inicial1 = &combi(\@matriu_puntuacio_donnors_1,\@initials1,\@matriu1);
my @exo_inicial2 = &combi(\@matriu_puntuacio_donnors_2,\@initials2,\@matriu2);


#modifiquem les columnes corresponents a les posicions absolutes i relatives
#segons els nucleotids trets per fer les pautes de lectura

my $init;
for($init=0; $init<scalar(@exo_inicial1); $init++) {

    $exo_inicial1[$init][1]=$exo_inicial1[$init][1]+1;
    $exo_inicial1[$init][2]=$exo_inicial1[$init][2]+1;
    $exo_inicial1[$init][3]=$exo_inicial1[$init][3]+1;
    $exo_inicial1[$init][4]=$exo_inicial1[$init][4]+1;
}

for($init=0; $init<scalar(@exo_inicial2); $init++) {

    $exo_inicial2[$init][1]=$exo_inicial2[$init][1]+2;
    $exo_inicial2[$init][2]=$exo_inicial2[$init][2]+2;
    $exo_inicial2[$init][3]=$exo_inicial2[$init][3]+2;
    $exo_inicial2[$init][4]=$exo_inicial2[$init][4]+2;
}


#cridem la matriu biaix_codificant per atribuir score a cada exo inicial
#que se li sumara a la puntuacio del donnor

my $r;
for ($r=0; $r < scalar(@exo_inicial0); $r++) {
    my $inicial_biaix_0 = &biaix_codificant(\%proporcions_codons,$exo_inicial0[$r][5]);
    $exo_inicial0[$r][6] = $inicial_biaix_0 + $exo_inicial0[$r][0];
    $exo_inicial0[$r][7] = "0"; #el frame dels exons inicials sempre sera 0
}


for ($r=0; $r < scalar(@exo_inicial1); $r++) {
    my $inicial_biaix_1 = &biaix_codificant(\%proporcions_codons,$exo_inicial1[$r][5]);
    $exo_inicial1[$r][6] = $inicial_biaix_1 + $exo_inicial1[$r][0];
    $exo_inicial1[$r][7] = "0"; #el frame dels exons inicials sempre sera 0

}


for ($r=0; $r < scalar(@exo_inicial2); $r++) {
    my $inicial_biaix_2 = &biaix_codificant(\%proporcions_codons,$exo_inicial2[$r][5]);
    $exo_inicial2[$r][6] = $inicial_biaix_2 + $exo_inicial2[$r][0];
    $exo_inicial2[$r][7] = "0"; #el frame dels exons inicials sempre sera 0
}


#ara comencem a predir els exons finals. Necessitem cridar la funcio combi
#i per aixo crearem una matriu d'stops adequada per aquesta funcio a partir
#de l'@matriu0, @matriu1, @matriu2.

my @terminals0;
my @terminals1;
my @terminals2;

my $istop;

for ($istop=0; $istop<scalar (@matriu0); $istop++) {
$terminals0[$istop][0]=$istop;
$terminals0[$istop][1]=length($matriu0[$istop][0]);
$terminals0[$istop][2]=0;
}

for ($istop=0; $istop<scalar (@matriu1); $istop++) {
$terminals1[$istop][0]=$istop;
$terminals1[$istop][1]=length($matriu1[$istop][0]);
$terminals1[$istop][2]=0;
}

for ($istop=0; $istop<scalar (@matriu2); $istop++) {
$terminals2[$istop][0]=$istop;
$terminals2[$istop][1]=length($matriu2[$istop][0]);
$terminals2[$istop][2]=0;
}


#cridem la funcio combi per aparellar acceptors amb stops

my @exo_terminal0 = &combi(\@terminals0,\@matriu_puntuacio_acceptors_0,\@matriu0);
my @exo_terminal1 = &combi(\@terminals1,\@matriu_puntuacio_acceptors_1,\@matriu1);
my @exo_terminal2 = &combi(\@terminals2,\@matriu_puntuacio_acceptors_2,\@matriu2);


my $fin;
for($fin=0; $fin<scalar(@exo_terminal1); $fin++) {

    $exo_terminal1[$fin][1]=$exo_terminal1[$fin][1]+1;
    $exo_terminal1[$fin][2]=$exo_terminal1[$fin][2]+1;
    $exo_terminal1[$fin][3]=$exo_terminal1[$fin][3]+1;
    $exo_terminal1[$fin][4]=$exo_terminal1[$fin][4]+1;
}

for($fin=0; $fin<scalar(@exo_terminal2); $fin++) {

    $exo_terminal2[$fin][1]=$exo_terminal2[$fin][1]+2;
    $exo_terminal2[$fin][2]=$exo_terminal2[$fin][2]+2;
    $exo_terminal2[$fin][3]=$exo_terminal2[$fin][3]+2;
    $exo_terminal2[$fin][4]=$exo_terminal2[$fin][4]+2;
}


#cridem la matriu biaix_codificant per atribuir score a cada exo terminal
#que se li sumara a la puntuacio de l'acceptor (a mes tenim en compte les diferents
#pautes de lectura i ens quedem la de puntuacio maxima)

my $m;
for ($m=0; $m < scalar(@exo_terminal0); $m++) {
    my $terminal_biaix_0 = &biaix_codificant (\%proporcions_codons,substr($exo_terminal0[$m][5],
                                              $matriu_puntuacio_acceptors_0[$m][3]));
    $exo_terminal0[$m][6] = $terminal_biaix_0+$exo_terminal0[$m][0];
    $exo_terminal0[$m][7] = $matriu_puntuacio_acceptors_0[$m][3];

}

for ($m=0; $m < scalar(@exo_terminal1); $m++) {
    my $terminal_biaix_1 = &biaix_codificant (\%proporcions_codons,substr($exo_terminal1[$m][5],
                                              $matriu_puntuacio_acceptors_1[$m][3])); 
    
    $exo_terminal1[$m][6] = $terminal_biaix_1+$exo_terminal1[$m][0];
    $exo_terminal1[$m][7] = $matriu_puntuacio_acceptors_1[$m][3];
}

for ($m=0; $m < scalar(@exo_terminal2); $m++) {
    my $terminal_biaix_2 = &biaix_codificant (\%proporcions_codons,substr($exo_terminal2[$m][5],
                                              $matriu_puntuacio_acceptors_2[$m][3]));

    $exo_terminal2[$m][6] = $terminal_biaix_2+$exo_terminal2[$m][0];
    $exo_terminal2[$m][7] = $matriu_puntuacio_acceptors_2[$m][3];
}



#a continuacio, predim possibles single exons aprofitant que ja tenim
#les matrius @initials i @terminals preparades per passar-les per la
#funcio combi que ens aparellara ATG amb STOP

my @single0 = &combi(\@terminals0,\@initials0,\@matriu0);
my @single1 = &combi(\@terminals1,\@initials1,\@matriu1);
my @single2 = &combi(\@terminals2,\@initials2,\@matriu2);


#el resultat el passem per la funcio biaix_codificant que ens donara la 
#puntuacio d'aquests exons ja que fins ara es 0 perque no hi ha puntuacio
#per inici i per terminacio.


my $t;
for ($t=0; $t < scalar(@single0); $t++) {
    my $single_biaix0 = &biaix_codificant (\%proporcions_codons,$single0[$t][5]);
    $single0[$t][6] = $single_biaix0 + $single0[$t][0];
    $single0[$t][7] = 0; #el frame dels single exons sempre sera 0
}

for ($t=0; $t < scalar(@single1); $t++) {
    my $single_biaix1 = &biaix_codificant (\%proporcions_codons,$single1[$t][5]);
    $single1[$t][6] = $single_biaix1 + $single1[$t][0];
    $single1[$t][7] = 0; #el frame dels single exons sempre sera 0
}

for ($t=0; $t < scalar(@single2); $t++) {
    my $single_biaix2 = &biaix_codificant (\%proporcions_codons,$single2[$t][5]);
    $single2[$t][6] = $single_biaix2 + $single2[$t][0];
    $single2[$t][7] = 0; #el frame dels single exons sempre sera 0
}



#unim les matrius per preparar les dades per imprir en format GFF i afegim dues noves columnes
#una que digui si l'exo es intern, terminal, inicial o single


my @unica_intern= &unir_tres_matrius(\@donnor_acc_0,\@donnor_acc_1,\@donnor_acc_2);
my $re;
for($re=0; $re<scalar(@unica_intern); $re++){
    $unica_intern[$re][8]="internal";
}



my @unica_inicial= &unir_tres_matrius(\@exo_inicial0,\@exo_inicial1,\@exo_inicial2);
for($re=0; $re<scalar(@unica_inicial); $re++){
    $unica_inicial[$re][8]="initial ";
}




my @unica_terminal= &unir_tres_matrius(\@exo_terminal0,\@exo_terminal1,\@exo_terminal2);
for($re=0; $re<scalar(@unica_terminal); $re++){
    $unica_terminal[$re][8]="terminal";
   
}


my @unica_single= &unir_tres_matrius(\@single0,\@single1,\@single2);
for($re=0; $re<scalar(@unica_single); $re++){
    $unica_single[$re][8]="single  ";
    
}


#finalment, nomes ens queda unir aquestes 4 matrius resultants en una
#sola que ens servira per mostrar el resultat final en format GFF.
#(podriem haver fet una funcio que unis matrius de 2 en 2, en lloc de 
#una que les uneix de tres en tres i una altra de quatre en quatre, pero
#vem pensar que aixi no l'hauriem de cridar tants cops)

my @matriu_final= &unir_quatre_matrius(\@unica_intern,\@unica_inicial,\@unica_terminal,\@unica_single);



#ja nomes queda ordenar la matriu final segons les puntuacions totals

@matriu_final= sort{${$a}[3] <=> ${$b}[3]}(@matriu_final);


#print!!!!
my $g;

chomp ($iden);
print "Resultados del análisis de la secuencia con identificador $iden:\n";
print "\n";
for ($g=0; $g<scalar(@matriu_final); $g++) {

    if($matriu_final[$g][4]-$matriu_final[$g][3]>=$ARGV[5]){ #per establir el 
                                                             #tamany minim de l'exo
    
    print "predictor\t$matriu_final[$g][8]\t$matriu_final[$g][3]\t$matriu_final[$g][4]\t";
    printf "\t%.3f",$matriu_final[$g][6];
    print "\t\t+\t$matriu_final[$g][7]\n";
    
    }

} 





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

# NOM: calcula_stops
# PROPOSIT: aquesta funcio ens busca els codons stop -TAA, TGA o TAG-
#          
# PARAMETRES: primer - vector codons on tenim la sequencia sencera i ens
#                      la llegeix de tres en tres nucleotids
# RETORNA: el resultat s'enregistra en un matriu on la primera columna
#          son els marcs de lectura oberts i la segona la posicio absoluta
#          de l'stop que el limita



sub calcula_stops { 

    my $sequencia=$_[0]; #sequencia sencera on volem trobar els orfs
    
    my $iorfs=0;         #numero d'orf
    my $posicio_codo=0;  #numero de codo en cada orf
    my @codons = ($sequencia=~m/.../g); #dividim la seq en codons
    my $resultat = "";   #l'inicialitzem amb la paraula buida
    my @orfs;            #vector de sortida amb tots els orfs


    while ($posicio_codo < scalar(@codons) ) { 
   
        #busquem els codons stop que poden ser tres-TAA, TAG, TGA-

	if ($codons[$posicio_codo] ne "TAA" && $codons[$posicio_codo] ne "TGA" && 
	    $codons[$posicio_codo] ne "TAG") { 
	    $resultat=$resultat.$codons[$posicio_codo]; 
	} 
	else  { 
	    $resultat=$resultat.$codons[$posicio_codo]; 
	    $orfs[$iorfs][0]=$resultat; 
       	    $orfs[$iorfs][1]=$posicio_codo*3+3;
	    $iorfs=$iorfs+1; 
	    $resultat=""; 
	} 

	$posicio_codo=$posicio_codo+1; 

    } 
    return @orfs; 
} 



#NOM: llegeix_matriu 
#PROPOSIT: aconseguir passar per totes les posicions de la matriu,
#          llegeix la matriu
#PARAMETRES: nom del fitxer de matriu de pesos (de donnors o d'acceptors 
#            respectivament)
#RETORNA: un hash de vectors on les claus son els nucleotids A, C, T,
#         G i un vector amb els pesos de cada lletra a les diferents 
#         posicions de la sequencia


sub llegeix_matriu_pesos { 

    
    my $nom_fitxer_matriu_pesos=$_[0]; 
    my %matriu_pesos;
    
    if (!open(MATRIU_DA,"<$nom_fitxer_matriu_pesos")) { 
	print "predictor.pl: impossible obrir $nom_fitxer_matriu_pesos\n"; 
	exit (1); 
    } 

    my $nucleotids;      
    my @nucleotids;      
    my $llegir_matriu=0; 
    my $linies=0; 
    my $posicio;


    while (<MATRIU_DA>) { 

	if ($llegir_matriu==0 && m/\AP0\s+(\w+)\s+(\w+)\s+(\w+)\s+(\w+)/) { 

	    $nucleotids[0] = $1; 
	    $nucleotids[1] = $2; 
	    $nucleotids[2] = $3; 
	    $nucleotids[3] = $4; 

	    $llegir_matriu = 1; 
	} 

	if ($llegir_matriu == 1 && 
	    m/\A\d+\s+([\-\d\.]+)\s+([\-\d\.]+)\s+([\-\d\.]+)\s+([\-\d\.]+)/) { 

	    $matriu_pesos{$nucleotids[0]}[$linies]=$1; 
	    $matriu_pesos{$nucleotids[1]}[$linies]=$2; 
	    $matriu_pesos{$nucleotids[2]}[$linies]=$3; 
	    $matriu_pesos{$nucleotids[3]}[$linies]=$4; 

	    $linies = $linies + 1; 
	} 

	if ($llegir_matriu == 1 && m/\AXX\s+(\d+)/) { 
	    $matriu_pesos{punt_de_tall} = $1; 
	    $llegir_matriu = 0; 
	} 
    } 

    close (MATRIU_DA); 

    return %matriu_pesos; 
} 


#NOM: calcula_donor_acceptor 
#PROPOSIT: calcular la posicio i puntuacio de donnors i accepors
#PARAMETRES: primer:matriu de pesos de donnors o d'acceptors 
#            segon:matriu que compta els diferents orfs i la posicio
#                  absoluta d'on acaben els STOPS
#RETORNA: una matriu amb tres columnes:orf; posicio del donnor/acceptor;
#         pes del donnor/acceptor

sub calcula_donor_acceptor { 

    my %matriu_calcula = %{$_[0]};  #matriu de pesos 
    my @matriu_alafuncio= @{$_[1]}; #matriu amb cada orf trobat i 
                                    #posicio absoluta de cada stop
    my $llindar =$_[2];                #llindar

    my $llargada_finestra=scalar(@{$matriu_calcula{A}}); 
    my $i;             #em serveix per desplacar la finestra
    my $line=0;        #desplaca la linia de la matriu de sortida on 
                       #guardo el resultat
    my $iorfs=0;       #numero d'orfs
    my $nt_finestra;   #subsequencia de nucleotids-->una finestra
    my @donor_acceptor_putative; 
    my @matriu_puntuacio;   
    

#definim una finestra; i definim un vector on posarem totes les finestres 
#possibles, on cada posicio d'aquest vector sera un nucleotid 
    
        
    while ($iorfs < scalar(@matriu_alafuncio)) {

        $i=0;
	while ($i < length($matriu_alafuncio[$iorfs][0])-$llargada_finestra) { 
	    
	    $nt_finestra=substr($matriu_alafuncio[$iorfs][0], $i, $llargada_finestra);
	    @donor_acceptor_putative=split(//,$nt_finestra);

	    
	    my $h; #recorrem els nucleotids dins d'un mateix parell 
                   #donnor-acceptor
	    my $r=0; #pes de cada nucleotid en la finestra on treballem

	    
	    for ($h=0; $h < $llargada_finestra; $h++) {
		$r=$r+$matriu_calcula{$donor_acceptor_putative[$h]}[$h]; 
	    }   

	    if ($r > $llindar) {  
		$matriu_puntuacio[$line][0]=$iorfs;
		$matriu_puntuacio[$line][1]=$i + $matriu_calcula{punt_de_tall}-1;
		$matriu_puntuacio[$line][2]=$r;

		$line=$line+1;
	    }

	    $i=$i+1;
            $r=0; 
	} 

	$iorfs=$iorfs + 1; 

    }
    
    return @matriu_puntuacio;

} 

#NOM: calcula_frame
#PROPOSIT: calcular el frame real al que pertenyen els exons que predim 
#PARAMETRES: primer: matriu de puntuacio d'acceptors que hem construit anteriorment 
#            segon: matriu on tenim els orfs i la posicio absoluta d'on acaben els STOPS
#            
#RETORNA: la matriu de puntuacio d'acceptors pero amb una columna mes (el frame)
      



sub calcula_frame{

    my @matriu_puntuacio_acceptors = @{$_[0]};
    my @matriu = @{$_[1]};
    my $w; 

for ($w=0; $w<scalar(@matriu_puntuacio_acceptors); $w++) {
    my $orfes= $matriu_puntuacio_acceptors[$w][0];
    if ($orfes==0) {
	my $residu;
        $residu=$matriu_puntuacio_acceptors[$w][1]%3;
	if ($residu==0) {
	    $matriu_puntuacio_acceptors[$w][3]=0;
	}
        if ($residu==1) {
	    $matriu_puntuacio_acceptors[$w][3]=2;
	}
      	if ($residu==2) {
       	    $matriu_puntuacio_acceptors[$w][3]=1;
       	}
    }
    else {
        my $residu=($matriu_puntuacio_acceptors[$w][1]-$matriu[$orfes-1][1])%3;
	if ($residu==0) {
	    $matriu_puntuacio_acceptors[$w][3]=0;
	}
        if ($residu==1) {
	    $matriu_puntuacio_acceptors[$w][3]=2;
	}
      	if ($residu==2) {
       	    $matriu_puntuacio_acceptors[$w][3]=1;
       	}
    }
}
    return @matriu_puntuacio_acceptors;
}

#NOM: combi
#PROPOSIT: combinar les dues matrius obtingudes a la funcio anterior
#          fent parelles acceptor-donnor
#PARAMETRES: primer:matriu amb posicio i puntuacio de donnor; segon:
#            matriu amb posicio i puntuacio d'acceptor i tercer:
#            matriu amb tots els orfs i la posicio absoluta de cada
#            stop
#RETORNA: matriu amb sis columnes: suma de puntuacions; posicio d'acceptor
#         relativa; posicio de donnor relativa; posicio d'acceptor absoluta;
#         posicio de donnor absoluta i sequencia entre acceptor i donnor.

sub combi {

    my @punt_donnors = @{$_[0]}; #matriu amb puntuacions i posicions de donnors
    my @punt_acceptors = @{$_[1]}; #matriu amb puntuacions i posicions 
                                   #d'acceptors
    my @matriu_combi= @{$_[2]}; #matriu amb orfs i posicions absolutes 
                                #de cada stop
   
    my $i=0;    #recorrem les posicions dels acceptors
    my $j=0;    #recorrem les posicions dels donnors per cada acceptor
    my $idx=0;  #fila de la matriu de sortida
    my $suma=0; #suma del pes de l'acceptor i el del donnor quan s'aparellin
                #formant un possible exo
    my $donnors_size=scalar(@punt_donnors);
    my $acceptor_size=scalar(@punt_acceptors);
    my @donnor_acc;


    for ($i=0 ; $i < $acceptor_size ; $i++) {
	
        #recoloca la j a lloc cada vegada que la $i incrementa una unitat
	
	$j=0;
	while ($j < $donnors_size && $punt_acceptors[$i][0] != $punt_donnors[$j][0]) { 
	    $j++;
	}

        #el que ve a continuacio nomes s'executa si el numero d'orf coincidex entre la
        #matriu dels donnors i la dels acceptors
 
	while ($j < $donnors_size && $punt_acceptors[$i][0] == $punt_donnors[$j][0]) { 
	    
	    if ($punt_acceptors[$i][1] < $punt_donnors[$j][1]) {
		$suma = $punt_acceptors[$i][2] + $punt_donnors[$j][2];
		$donnor_acc[$idx][0]=$suma;
		$donnor_acc[$idx][1]=$punt_acceptors[$i][1];  #posicions relatives acc
		$donnor_acc[$idx][2]=$punt_donnors[$j][1];   #posicions relatives donnors
		
		if($punt_acceptors[$i][0]==0) {
		    $donnor_acc[$idx][3]=$punt_acceptors[$i][1];
		    $donnor_acc[$idx][4]=$punt_donnors[$j][1];
		}
		else {
		    $donnor_acc[$idx][3]=$punt_acceptors[$i][1]+
                                         $matriu_combi[$punt_acceptors[$i][0]-1][1];
		    $donnor_acc[$idx][4]=$punt_donnors[$j][1]+
                                         $matriu_combi[$punt_acceptors[$i][0]-1][1];
		}
		
		my $seq=$matriu_combi[$punt_acceptors[$i][0]][0];
		my $first=$donnor_acc[$idx][1];
		my $llarg=$donnor_acc[$idx][2]-$donnor_acc[$idx][1];
		$donnor_acc[$idx][5] = substr($seq,$first,$llarg);

		$idx++;
	    }
	    $j++;
	}
    }
    
    return @donnor_acc;
    
}



#NOM: biaix_codificant
#PROPOSIT: calcula la rao de versemblanca de les sequencies compreses
#          entre acceptors i donnors
#PARAMETRES: primer: fitxer on hi ha cada codo i la seva proporcio
#             segon: fitxer de la sequencia limitada per un acceptor 
#                    i un donnor
#             tercer: matriu que ens ha retornat la funcio anterior
#RETORNA: una matriu on tenim els elements de la matriu d'entrada amb
#         una nova columna on hi ha l'score de cada exo putatiu predit


sub biaix_codificant {

    my %propcod = %{$_[0]}; #matriu biaix codificant 
    my $sequencia = $_[1];  #correspon a la columna de la matriu 
                            #@donnor_acc on hi hem guardat el tros de
                            #sequencia entre l'acceptor i el donnor    

    my @per_codons = ($sequencia =~ m/.../g);

    my $r=0; #resultat de la suma de pesos de cada codo al passar la matriu
    my $n;   #numero del codo 
    

    for ($n = 0; $n < scalar(@per_codons); $n++) {
	$r=$r+log($propcod{$per_codons[$n]})-log(1/64);
    }
    return $r;
}




#NOM: busca_ATG
#PROPOSIT: buscar tots els codons ATG per predir els exons inicials
#PARAMETRES: sequencia sencera en cada pauta de lectura
#RETORNA: matriu on la primera columna sera el numero de orf
#         i la segona la posicio relativa de cada ATG  


sub busca_ATG {

    my @matriu_busca= @{$_[0]}; 
    my $nt=0;    
    my $fila;
    my $i=0;
    my @initials;

for ($fila=0; $fila < scalar(@matriu_busca); $fila++) {
    
    my @codons=($matriu_busca[$fila][0] =~ m/.../g);
   
    while ($nt < scalar(@codons)) {
    
   
    if($codons[$nt] ne "ATG") {
	$nt=$nt+1;
    }

   

    else  {
	$initials[$i][0]=$fila;
	$initials[$i][1]=$nt*3;
	$initials[$i][2]=0; #igualem aquesta columna a 0 per poder
                            #introduir-la a la funcio combi 
                            #(es la seva puntuacio)

	$i=$i+1;
    }
    
    $nt=$nt+1;

    }
    
    $nt=0;

}
    return @initials;
}


#NOM: unir_tres_matrius
#PROPOSIT: unir les matrius que introduim una sota de l'altre
#PARAMETRES: les tres matrius a unir
#RETORNA:una sola matriu amb tantes files com la suma de les
#        files de les tres matrius

sub unir_tres_matrius {

    my @matriu_primera = @{$_[0]};
    my @matriu_segona = @{$_[1]};
    my @matriu_tercera = @{$_[2]};
    my $fila=0; #recorrera totes les files de la matriu final
    my @matriu_unica;


  my $j;
  if (scalar(@matriu_primera) != 0) {
      for ($j=0; $j<scalar(@matriu_primera); $j++) {
	   $matriu_unica[$fila][0]= $matriu_primera[$j][0];
	   $matriu_unica[$fila][1]= $matriu_primera[$j][1];
	   $matriu_unica[$fila][2]= $matriu_primera[$j][2];
	   $matriu_unica[$fila][3]= $matriu_primera[$j][3];
	   $matriu_unica[$fila][4]= $matriu_primera[$j][4];
	   $matriu_unica[$fila][5]= $matriu_primera[$j][5];
	   $matriu_unica[$fila][6]= $matriu_primera[$j][6];
           $matriu_unica[$fila][7]= $matriu_primera[$j][7];
	   $fila=$fila+1;
       }
  }

  my $q;
  if (scalar(@matriu_segona) != 0) {
      for ($q=0; $q<scalar(@matriu_segona); $q++) {
	   $matriu_unica[$fila][0]= $matriu_segona[$q][0];
	   $matriu_unica[$fila][1]= $matriu_segona[$q][1];
	   $matriu_unica[$fila][2]= $matriu_segona[$q][2];
	   $matriu_unica[$fila][3]= $matriu_segona[$q][3];
	   $matriu_unica[$fila][4]= $matriu_segona[$q][4];
	   $matriu_unica[$fila][5]= $matriu_segona[$q][5];
	   $matriu_unica[$fila][6]= $matriu_segona[$q][6];
           $matriu_unica[$fila][7]= $matriu_segona[$q][7];
	   $fila=$fila+1;
       }
  }

   my $g;
   if (scalar(@matriu_tercera) != 0) {
       for ($g=0; $g<scalar(@matriu_tercera); $g++) {
   	   $matriu_unica[$fila][0]= $matriu_tercera[$g][0];
	   $matriu_unica[$fila][1]= $matriu_tercera[$g][1];
	   $matriu_unica[$fila][2]= $matriu_tercera[$g][2];
	   $matriu_unica[$fila][3]= $matriu_tercera[$g][3];
	   $matriu_unica[$fila][4]= $matriu_tercera[$g][4];
	   $matriu_unica[$fila][5]= $matriu_tercera[$g][5];
	   $matriu_unica[$fila][6]= $matriu_tercera[$g][6];
           $matriu_unica[$fila][7]= $matriu_tercera[$g][7];
	   $fila=$fila+1;
       }
   }
return @matriu_unica;
}



#NOM: unir_quatre_matrius
#PROPOSIT: uneix 4 matrius donades
#PARAMETRES: les 4 matrius
#RETORNA: una matriu unica suma de les quatre

sub unir_quatre_matrius {

    my @matriu_primera = @{$_[0]};
    my @matriu_segona = @{$_[1]};
    my @matriu_tercera = @{$_[2]};
    my @matriu_quarta = @{$_[3]};
    my $fila=0; #recorrera totes les files de la matriu final
    my @matriu_unica;


  my $j;
  if (scalar(@matriu_primera) != 0) {
      for ($j=0; $j<scalar(@matriu_primera); $j++) {
	   $matriu_unica[$fila][0]= $matriu_primera[$j][0];
	   $matriu_unica[$fila][1]= $matriu_primera[$j][1];
	   $matriu_unica[$fila][2]= $matriu_primera[$j][2];
	   $matriu_unica[$fila][3]= $matriu_primera[$j][3];
	   $matriu_unica[$fila][4]= $matriu_primera[$j][4];
	   $matriu_unica[$fila][5]= $matriu_primera[$j][5];
	   $matriu_unica[$fila][6]= $matriu_primera[$j][6];
	   $matriu_unica[$fila][7]= $matriu_primera[$j][7];
	   $matriu_unica[$fila][8]= $matriu_primera[$j][8];

	   $fila=$fila+1;
       }
  }

  my $q;
  if (scalar(@matriu_segona) != 0) {
      for ($q=0; $q<scalar(@matriu_segona); $q++) {
	   $matriu_unica[$fila][0]= $matriu_segona[$q][0];
	   $matriu_unica[$fila][1]= $matriu_segona[$q][1];
	   $matriu_unica[$fila][2]= $matriu_segona[$q][2];
	   $matriu_unica[$fila][3]= $matriu_segona[$q][3];
	   $matriu_unica[$fila][4]= $matriu_segona[$q][4];
	   $matriu_unica[$fila][5]= $matriu_segona[$q][5];
	   $matriu_unica[$fila][6]= $matriu_segona[$q][6];
	   $matriu_unica[$fila][7]= $matriu_segona[$q][7];
           $matriu_unica[$fila][8]= $matriu_segona[$q][8];

	   $fila=$fila+1;
       }
  }

   my $g;
   if (scalar(@matriu_tercera) != 0) {
       for ($g=0; $g<scalar(@matriu_tercera); $g++) {
	   $matriu_unica[$fila][0]= $matriu_tercera[$g][0];
	   $matriu_unica[$fila][1]= $matriu_tercera[$g][1];
	   $matriu_unica[$fila][2]= $matriu_tercera[$g][2];
	   $matriu_unica[$fila][3]= $matriu_tercera[$g][3];
	   $matriu_unica[$fila][4]= $matriu_tercera[$g][4];
	   $matriu_unica[$fila][5]= $matriu_tercera[$g][5];
	   $matriu_unica[$fila][6]= $matriu_tercera[$g][6];
	   $matriu_unica[$fila][7]= $matriu_tercera[$g][7];
	   $matriu_unica[$fila][8]= $matriu_tercera[$g][8];
	      
	   $fila=$fila+1;
       }
   }

   my $p;
   if (scalar(@matriu_quarta) != 0) {
       for ($p=0; $p<scalar(@matriu_quarta); $p++) {
	   $matriu_unica[$fila][0]= $matriu_quarta[$p][0];
	   $matriu_unica[$fila][1]= $matriu_quarta[$p][1];
	   $matriu_unica[$fila][2]= $matriu_quarta[$p][2];
	   $matriu_unica[$fila][3]= $matriu_quarta[$p][3];
	   $matriu_unica[$fila][4]= $matriu_quarta[$p][4];
	   $matriu_unica[$fila][5]= $matriu_quarta[$p][5];
	   $matriu_unica[$fila][6]= $matriu_quarta[$p][6];
	   $matriu_unica[$fila][7]= $matriu_quarta[$p][7];
	   $matriu_unica[$fila][8]= $matriu_quarta[$p][8];
	  
	   $fila=$fila+1;
       }
   }
return @matriu_unica;
}