#!/usr/bin/perl -w

use strict;


############################    ENTRADA DE DADES    ############################



# si el nombre d'arguments al programa es mes petit que 7
# no l'hem pogut cridar amb els noms de fitxers que hem d'obrir
# despres

if (scalar(@ARGV) < 7) {
	print "buscoexons.pl fitxercodons fitxeracc fitxerdonn fitxerseq cutoffacc cutoffdonn llargminim\n";
	exit(1);
};


# assigno els noms del vector d'arguments (parametres) a 7
# variables amb noms mes significatius

my $fitxercod = $ARGV[0];
my $fitxeracc = $ARGV[1];
my $fitxerdonn = $ARGV[2];
my $fitxerseq = $ARGV[3];
my $cutoffacc = $ARGV[4];
my $cutoffdonn = $ARGV[5];
my $llargminim = $ARGV[6];

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


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

my %pcodons;

while (<PROPCODONS>) {

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

	m/(\w+)\s+(\d+\.\d+)/o;
	
	# 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
	
	$pcodons{$1}=$2;
};


close(PROPCODONS);


# obrim el fitxer especificat al segon argument passat al
# programa des del shell

if (!open(SEQUENCIA,"< $fitxerseq")) {
    print "buscoexons.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 diferents trossos de la sequencia
# canviarem el marcador de registre de '\n' a 


my $primeralinia = <SEQUENCIA>;

$/ = '>';

my $x = <SEQUENCIA>;

close(SEQUENCIA);

$/ = "\n";


# voldriem tenir la sequencia en una sola variable sense la primera linia
# '>identificador'. per tal de fer aixo ficarem el contingut llegit
# (el registre sencer) en un vector on cada posicio contindra una linia

my @v = split(/\n/og,$x);

my $seq="";
my $i=0;

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

# Ara poso la sequencia en un vector @nt (de nucleotids), on cada posicio
# correspon a un nucleotid

my @nt = split(//og,$seq);


# encaixo la variable $seq amb una expressio regular que ens donara els 3
# vectors de codons @c0, @c1 i @c2, un per cada pauta d lectura.

my @c0;
my @c1;
my @c2;

@c0 = ($seq=~ m/.../og);
$seq = substr($seq,1);
@c1 = ($seq=~ m/.../og);
$seq = substr($seq,1);
@c2 = ($seq=~ m/.../og);


# Crido la funcio de lectura de les matrius de pesos dels donnors
# i acceptors i les emmagatzemo en els respectius hashos de vectors

my %hvacc;
my %hvdonn;

%hvacc = &llegeix_matriu($fitxeracc);
%hvdonn = &llegeix_matriu($fitxerdonn);

# Averiguo llargada de la sequencia de la matriu per tal d'utilitzar-la


			############################
##########################    BUSCO EXONS!!!!!    #####################
			############################

# Crido la funcio STOPs que em donarà la posició dels stop codons
# per a cada pauta de lectura.

my @vstop0 = &STOPs(\@c0);
my @vstop1 = &STOPs(\@c1);
my @vstop2 = &STOPs(\@c2);


# Crido la funció ORFs que em donarà les sequencies entre codons stop
# en cada posició del vector resultant, per cada pauta de lectura.

my @orfs0 = &ORFs(\@vstop0,\@c0);
my @orfs1 = &ORFs(\@vstop1,\@c1);
my @orfs2 = &ORFs(\@vstop2,\@c2);


# Crido la funcio predir_da, que em permet predecir donnors i acceptors amb les
# matrius de pesos corresponents i emmagatzemar la seva posicio (nucleotid
# apartir del qual ENTRA dins l'exo) en la sequencia original en una de les
# files de la matriu de sortida i la seva puntuacio, que es mes alta o igual 
# al cutoff (variable que prove del 5e parametre, recollit des del shell)en 
#l'altra fila.


my @predaccept0 = &predir_da(\@orfs0,$cutoffacc,\%hvacc,\@vstop0,0);
my @predaccept1 = &predir_da(\@orfs1,$cutoffacc,\%hvacc,\@vstop1,1);
my @predaccept2 = &predir_da(\@orfs2,$cutoffacc,\%hvacc,\@vstop2,2);

my @preddonors0 = &predir_da(\@orfs0,$cutoffdonn,\%hvdonn,\@vstop0,0);
my @preddonors1 = &predir_da(\@orfs1,$cutoffdonn,\%hvdonn,\@vstop1,1);
my @preddonors2 = &predir_da(\@orfs2,$cutoffdonn,\%hvdonn,\@vstop2,2);



# Crido la funcio buscoexons, que em permet recorrer la matriu d'acceptors de
# determinada pauta de lectura i
# aparellar-la amb TOTS els donnors del propi orf, tornant-nos una matriu on les
# files son parells donnor-acceptor i les columnes corresponen a posico d'inici 
# la primera fila, posicio de fi la segona i l'score la tercera.


my @exons0 = &BuscoExons(\@predaccept0,\@preddonors0,0,\@vstop0,$llargminim);
my @exons1 = &BuscoExons(\@predaccept1,\@preddonors1,1,\@vstop1,$llargminim);
my @exons2 = &BuscoExons(\@predaccept2,\@preddonors2,2,\@vstop2,$llargminim);




# Crido la funcio exoscore, q em retorna una matriu identica a
# la d'entrada, excepte la 3a fila, on ha sumat la versemblanca 
# de la sequencia de l'exo, mesurada segons els codons d'aquesta.


 
my @definitiva0 = &exoscore(\@exons0,\@c0,\%pcodons,0);
my @definitiva1 = &exoscore(\@exons1,\@c1,\%pcodons,1);
my @definitiva2 = &exoscore(\@exons2,\@c2,\%pcodons,2);


#Crido la funcio matriunica, que uneix en una sola matriu els exons de les 3 pautes de lectura

my @unica = &matriunica(\@definitiva0,\@definitiva1,\@definitiva2);



# Crido la funcio ordre  que m'ordena per ordre de scores la matriu d'entrada.

my @redefinitiva = &ordre(\@unica);



# Printo els resultats en format GFF

my $j = 0;


if (scalar(@redefinitiva) > 0) {
    while ($j < scalar(@{$redefinitiva[0]})) {
	print "$ARGV[3]\tbuscoexon\tinternal\t$redefinitiva[0][$j]\t$redefinitiva[1][$j]\t";
	printf "%.2f",$redefinitiva[2][$j];
	print "\t+\t$redefinitiva[3][$j]\n";
	$j = $j + 1;
    }
}






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

# NOM: llegeix_matriu
# PROPOSIT: llegir la matriu de pesos i emmagatzemar-la en un hash
#           de vectors.
# PARAMETRES: primer -  nom del fitxer de matriu de pesos
#
# RETORNA: un hash de vectors on trobem les lletres A,C,G,T com a
# claus i un vector amb els pesos de cada lletra a les diferents
# posicions de la sequencia com a valors.

sub llegeix_matriu() {
 
    my %hvpesos;
    my $matriu = $_[0];
    
# Ara obriré el fitxer de la matriu de pesos

    if (!open(MATRIPESOS,"< $matriu")) {
	print "buscoexons.pl: impossible obrir $matriu\n";
	exit(1);
    };
  
# I la anire modificant per tal de llegir la informacio rellevant i treure la palla.
  
    my $flagvalor=0;
    my $linia=0;
	
    while (<MATRIPESOS>) {
		
# ha de començar a llegir en el simbol P0

	if (m/P0/o) {
	    $flagvalor=1;
	    next;
	};
#Quan arribem al simbol XX aturem l'entrada de dades.

	if ( m/XX\s+(\d+)/o ) {
	    $flagvalor=0;
	    $hvpesos{P}=$1;
	    next;
	};
		

	if ($flagvalor==1) {
			
	    m/\d+\s+([+-]?\d*\.?\d+)\s+([+-]?\d*\.?\d+)\s+([+-]?\d*\.?\d+)\s+([+-]?\d*\.?\d+)/o;

# Ho poso les dades en un hash de vectors que tindra el nom %hvpesos
	    
	    $hvpesos{A}[$linia]=$1;
	    $hvpesos{C}[$linia]=$2;
	    $hvpesos{G}[$linia]=$3;
	    $hvpesos{T}[$linia]=$4;

	    $linia=$linia+1;
	}; 
	
    };
    close(MATRIPESOS);
    
    return (%hvpesos);
	
} # llegeix_matriu








# NOM: STOPs
# PROPOSIT: troba els stop codons i els emmagatzema en un vector la posicio
#           del nucleòtid (no del codo) en la que es troba.
# PARAMETRES: primer -  vector de condons
#
# RETORNA: Un vector on hi son cada component del vector es una posicio dins
#		d'una sequència on s'hauria eliminat 0, 1 o 2 nucleotids de la
#		sequència original, segons la seva pauta de lectura


sub STOPs() {
    my @vc = @{$_[0]};
    my @vstop;
    my $pos=0;
    my $i=0;
    while ($i < scalar(@vc)) {
	if ($vc[$i] eq "TAA" || $vc[$i] eq "TGA" || $vc[$i] eq "TAG") {
	    $vstop[$pos] = $i * 3;
	    $pos = $pos + 1;         
	};
	$i = $i + 1;
    }; # while $i
	
    return @vstop;
} # STOPs






# NOM: ORFs
# PROPOSIT: Agafa els stop codons i els utilitza per trobar les sequencies
#		intermitges, és a dir, els ORFs, de cada pauta de lectura.
#           del nucleotid (no del codo) en la que es troba.
# PARAMETRES: primer -  vector de stops
#             segon - vector de codons de la pauta de lectura propia
# RETORNA: Un vector amb cada posició corresponent a la sequencia d'un orf.


sub ORFs() {
    my @vstop = @{$_[0]};
    my @vc = @{$_[1]};
    my @orfs;
    my $seq = "";

    my $i = 0;
    while ($i < scalar(@vc)) {
	$seq = $seq . $vc[$i];
	$i = $i + 1;
    };
    $i=0;
    if (scalar(@vstop) > 0) {
	$orfs[$i] = substr($seq, 0, $vstop[0]);
	$i = $i + 1;
    }
    else {
        $orfs[0] = $seq;
	return @orfs;
    }
    
    $i = 1;
        while ($i < scalar(@vstop)) {
	$orfs[$i] = substr($seq , $vstop[$i-1] + 3 , $vstop[$i] - $vstop[$i-1]);
	$i = $i + 1;
    };
    
      $orfs[$i] = substr($seq , $vstop[$i-1] + 3);
    
	
    return (@orfs);
} # ORFs











# NOM: predir_da
# PROPOSIT: llegeix un vector d'ORFs, un cutoff i una matriu de pesos
#           on a $matriupesos{P} hi ha la posicio d'inici/fi de l'exo.
#
# PARAMETRES: primer -  vector d'ORFs
#             segon - cutoff apartir del qual s'accepta el donnor/acceptor
#             tercer - hash de vectors  amb els pesos dels nucleotids
#                      segons la seva posicio en el donnor/acceptor
#             quart - vector dels codons STOPs
#             cinque - pauta de lectura per corretgir la posicio
#             
# RETORNA: una matriu amb un columna per cada donor/acceptor i dues
#          files: una per posicio respecte la sequencia original, i una altra per
#          puntuacio (pes).


sub predir_da() {
    my @orfs = @{$_[0]};
    my $cutoff = $_[1];
    my %matriuda = %{$_[2]};
    my @vstop = @{$_[3]};
    my $pauta = $_[4];

    my $h=0;
    my $i;
    my $j;
    my $k;
    my @unposibleda;
    my $llargada;
    my $unposibleda;
    my $unorf;
    my @unorf;
    my $puntuac;
    my @matrpospunt;
    my $lletra;
    

    $llargada = scalar(@{$matriuda{A}});
    
    $j = 0;
    while ($j < scalar(@orfs)) {
	$unorf = $orfs[$j];
	@unorf = split(//og,$unorf);
	
	$k=0;
#Obro una finestra
	while ($k <= scalar(@unorf) - $llargada) {
	    $unposibleda = substr($unorf , $k , $llargada);
	    @unposibleda = split(//og,$unposibleda);
	    my $puntuac;
# Puntuo la finestra	
	    $i = 0;
	    $puntuac = 0.0;
	    while ($i < $llargada) {
		$lletra = $unposibleda[$i];
		$puntuac = $puntuac + $matriuda{$lletra}[$i];
		$i = $i + 1;
	    };
# Selecciono segons cutoff
	    if ($puntuac >= $cutoff) {

		if($j == 0){
		    $matrpospunt[0][$h] = $k + $matriuda{P} - 1 + $pauta;
		    $matrpospunt[1][$h] = $puntuac;
		    $h = $h + 1;
		}else{
		    $matrpospunt[0][$h] = $vstop[$j-1] + 3 + $k + $matriuda{P} - 1 + $pauta;
		    $matrpospunt[1][$h] = $puntuac;
		    $h = $h + 1;
		};
	    };
	   
	    $k = $k + 1;
	};
	$j = $j + 1;
    };
    return @matrpospunt;
} # predir_da










# NOM: BuscoExons
# PROPOSIT: recorrer la matriu d'acceptors de
#           determinada pauta de lectura i
#           aparellar-la amb TOTS els donnors del propi orf.
#
# PARAMETRES: primer - matriu d'acceptors
#             segon - matriu de donnors
#             tercer - pauta de lectura
#             quart - vector de codons, el qual no esta  corretgit
#                     segons la pauta de lectura, motiu pel que
#                     aquesta s'inclou com a parametre
#             cinque - llargada minima de l'exó (es un argument que no ha de baixar de 18)
#
# RETORNA: una matriu on les files son parells donnor-acceptor 
#         i les columnes corresponen a posicio d'inici 
#         la primera fila, posicio de fi la segona i l'score la tercera.


sub BuscoExons() {
    my @matracc = @{$_[0]};
    my @matrdonn = @{$_[1]};
    my $pauta = $_[2];
    my @vstop = @{$_[3]};
    my $llargmin = $_[4];

    
    my $h = 0; # itera per les posicions dels codons de stop
    my $i = 0; # itera pels acceptors
    my $j = 0; # itera pels donors
    my $k = 0; # va afegint exons
    my $flagvalor; # es 1 quan entre l'acceptor i el donnor hi ha un stop
    my @matrexons = ();


    if (scalar(@matracc) == 0 || scalar(@matrdonn) == 0) {
	return @matrexons;
    };

    my $numdonn = @{$matrdonn[0]};
    my $numacc = @{$matracc[0]};

    while($i < $numacc){
	$j = 0;
	while($j < $numdonn){
	    if($matrdonn[0][$j] > $matracc[0][$i] + $llargmin){ 
		$flagvalor = 0;
		$h = 0;
		while ($h < scalar(@vstop)){
		    if($matrdonn[0][$j] > ($vstop[$h] + $pauta) && $vstop[$h] + $pauta > $matracc[0][$i]){
			$flagvalor = 1;
			last;
		    }else{
			$h = $h + 1;
		    };
		};
		if($flagvalor == 0){
		    $matrexons[0][$k] = $matracc[0][$i];  
		    $matrexons[1][$k] = $matrdonn[0][$j];  
		    $matrexons[2][$k] = $matrdonn[1][$j] + $matracc[1][$i];
		    $k = $k + 1;
		} else {
		    $i = $i + 1;
		    last;
		};
	    } else {
		$j = $j + 1;
	    };
	    $j = $j + 1;
	};
	$i = $i + 1;
    };
    return @matrexons;
}# BuscoExons
    







# NOM: Exoscore
# PROPOSIT: calcula la rao de versemblanca de les sequencies compreses entre
#           acceptors i donnors, amb una pauta de lectura donada
# PARAMETRES: primer - matriu on tenim en files acceptor, donnor, score dels dos sumat.
#             segon  - vector de codons de la pauta de lectura corresponent
#             tercer - hash on tenim per a cada codo, la seva versemblanca
#             quart - pauta de lectura
# RETORNA: una matriu on tenim els elements de la matriu d'entrada, amb la      
#          3a fila modificada, ja que sumem la rao de versemblanca en escala
#          logaritmica de l'exo en concret al score d'acceptor i donnor.
        
    
    sub exoscore () {
	my @matrexons = @{$_[0]};
	my @codons = @{$_[1]};
	my %pcodons = %{$_[2]};
	my $pauta = $_[3];
	
	if (scalar(@matrexons) == 0){
	    return @matrexons;
	};
	
	my @versacum; # presenta a cada posicio la suma de la 
 # versemblanca del codo, mes la dels codons anteriors
	
 # calculo la rao de versemblanca en escala logaritmica
	
	my $r = 0.0;
	my $i = 0;
	my $v;

        $versacum[0]=0.0;

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

	    $v = $codons[$i];
	    $r = log($pcodons{$v}) - log(1/64);

	    $versacum[$i+1] = $r + $versacum[$i]; 

	    $i = $i + 1;
	}; # while

	my $j; # recorre les posicions dels acceptors
	my $k; # recorre les posicions dels donnors 
	my $h = 0; # indica la fila de la matriu 
	
	while ($h < scalar(@{$matrexons[0]})){
	    $j = $matrexons[0][$h] - $pauta;
	    $k = $matrexons[1][$h] - $pauta;
	    
	    my $l = 0; # ens diu el codo "sencer" inicial, on comenca a contar.
	    my $m = 0; # ens diu el codo "sencer" final, on acaba de contar.
	    
	    if ($j % 3 == 0){
		$l = ($j / 3) - 1;
	    }else{
		$l = $j / 3;
	    };
	    if ($k % 3 == 2){
		$m = $k / 3;
	    }else{
		$m = ($k / 3) - 1;
	    };
	    my $versexo; # indica la versemblanca de l'exo
	    
	    $versexo = $versacum[$m+1] - $versacum[$l+1];

	    $matrexons[2][$h] = $matrexons[2][$h]+ $versexo;
	    
	    $h = $h + 1;
	}; # while
	
	return @matrexons;
    }# exoscore
    
    
 
    
    

# NOM: matriunica
# PROPOSIT: unificar les matrius d cada pauta de lectura
#
# PARAMETRES: primer - matriu d'exons de la pauta 0
#             segon - matriu d'exons de la pauta 1
#             tercer - matriu d'exons de la pauta 2
#
# RETORNA: Una matriu com les d'entrada, amb tots els exons, però amb   
#         una fila mes, on  apareix la pauta a la que corresponen els exons.
    

    sub matriunica () {
	my @matriu0 = @{$_[0]};
	my @matriu1 = @{$_[1]};
	my @matriu2 = @{$_[2]};
	
	my @matriunica;
	
	my $i = 0;
	
	my $j = 0;
	unless (scalar(@matriu0) == 0){
	while ($j < scalar(@{$matriu0[0]})){
	    $matriunica[0][$i] = $matriu0[0][$j];
	    $matriunica[1][$i] = $matriu0[1][$j];
	    $matriunica[2][$i] = $matriu0[2][$j];
	    $matriunica[3][$i] = 0;
	    $i = $i + 1;
	    $j = $j + 1;
	}
    }
	my $k = 0;
	unless (scalar(@matriu1) == 0){
	    while ($k < scalar(@{$matriu1[0]})){
		$matriunica[0][$i] = $matriu1[0][$k];
		$matriunica[1][$i] = $matriu1[1][$k];
		$matriunica[2][$i] = $matriu1[2][$k];
		$matriunica[3][$i] = 1;
		$i = $i + 1;
		$k = $k + 1;
	    }
	}
	
	my $l = 0;
	unless (scalar(@matriu2) == 0){
	    while ($l < scalar(@{$matriu2[0]})){
		$matriunica[0][$i] = $matriu2[0][$l];
		$matriunica[1][$i] = $matriu2[1][$l];
		$matriunica[2][$i] = $matriu2[2][$l];
		$matriunica[3][$i] = 2;
		$i = $i + 1;
		$l = $l + 1;
	    }
	}
	return @matriunica;
    }







    
# NOM: ordre
# PROPOSIT: ordenar la matriu d'exons segons el score
#
# PARAMETRES: primer - matriu d'exons amb el seu score
#             
# RETORNA: Una matriu com la d'entrada, però amb les columnes (exons) 
#          ordenades segons l'score conjunt (acceptor/donnor/sequencia)
    
    sub ordre (){
	
	my @matrescores = @{$_[0]};
    
	my @ordrescores;
	my @scores;
	
	my $h = 0;


	if (scalar(@matrescores) == 0){
	    return @matrescores;
	};

	while($h < scalar(@{$matrescores[0]})){
	    $scores[$h] = $matrescores[2][$h];
	    $h = $h + 1;
	};

	my $k = 0;
	my $pos;
	my $max;
	while($k < scalar(@scores)){

	    for($i = 0, $max = -99999999; $i < scalar(@scores); $i++){
		if($scores[$i] > $max){
		    $max = $scores[$i];
		    $pos = $i;
		};
	    }; #for

	    $scores[$pos] = -99999999;
	    $ordrescores[3][$k] = $matrescores[3][$pos];
	    $ordrescores[2][$k] = $matrescores[2][$pos];
	    $ordrescores[1][$k] = $matrescores[1][$pos];
	    $ordrescores[0][$k] = $matrescores[0][$pos];
 
	    $k = $k + 1;
	};#while
	    
	return @ordrescores;
    }# ordre;