#! /usr/bin/perl -w

use strict;

# Este programa localiza los ORF's situados entre codones stop, sin tener en cuenta el primer y ultimo ORF's incompletos ni aquellos ORF's "vacios" que se
# generan cuando el programa encuentra dos codones stop seguidos. Ademas, esto lo hace en direccion forward unicamente y para las tres pautas de lectura posibles.
# Este programa solo sera util para la prediccion de zonas codificantes en secuencias procariotas porque no utiliza matrices de pesos para detectar acceptor y
# donor sites, imprescindibles para la determinacion de exones e intrones en secuencias eucariotas.

# Si el numero de argumentos escritos despues del nombre del programa es inferior a 2 no lo podremos ejecutar, el programa nos devolvera un texto de ayuda 
# indicando que argumentos necesita para poder ser ejecutado correctamente.

if (scalar(@ARGV) < 2) {
    print STDERR "debe introducir: orffinder.pl archivo_secuencia_fasta.fa archivo_tabla_uso_codones.txt\n";
    exit(1);
}

# Si el programa no encuentra el archivo fasta especificado nos dara el siguiente mensaje de error, si lo encuentra cargara todo su contenido en un filehandle.

if (!open(FASTAFILE, "< $ARGV[0]")) {
    print STDERR "imposible abrir $ARGV[0]\n";
    exit(1);
}

# Si el programa no encuentra el archivo con la tabla de uso de codones nos dara el siguiente mensaje de error, si lo encuentra cargara todo su contenido en un
# filehandle.

if (!open(CODONUSAGE, "< $ARGV[1]")) {
    print STDERR "imposible abrir $ARGV[1]\n";
    exit(1);
}

# Declaramos el hash que contendra todos los ORF's encontrados para cada pauta de lectura.

my %orfs = ();

# Asignamos el contenido del filehandle FASTAFILE a un vector con el mismo nombre para tener cada base de la secuencia de DNA en una posicion del vector, despues
# cerramos el filehandle.

my @fastafile = <FASTAFILE>;

close(FASTAFILE);

# Declaramos un hash que contendra la tabla de uso de codones porque nos interesa tener los tripletes como las keys del hash y las probabilidades de los diferentes 
# codones como los values del hash.

my %codonusage = ();

# Pasamos el contenido del filehandle CODONUSAGE al hash mediante una expresion regular, de modo que asignamos la categoria de key a lo que esta en la primera 
# columna de la tabla ($1) y la categoria de value a la segunda columna de la tabla ($2). Finalmente cerramos el filehandle.

while (<CODONUSAGE>) {
    m/(\w+)\s+(\d+\.\d+)/o;
    $codonusage{$1}=$2;
}

close(CODONUSAGE);

# Declaramos una nueva variable ($seq) que contendra la secuencia de DNA unicamente y otra variable ($id) que contendra el identificador (primera linea del
# archivo fasta). A partir de ahora trabajaremos unicamente sobre $seq, al no tratarse de un vector ya no tenemos cada nucleotido en una posicion sino toda la
# secuencia como un unico bloque.

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

while ($i < scalar(@fastafile)) {
    chomp($fastafile[$i]);
    if ($fastafile[$i]=~m/^>/o) {
	$fastafile[$i]=~s/^>//o;
	$id = $fastafile[$i];
    }	
    else {
	$seq = $seq.$fastafile[$i];
	}
    $i = $i+1;
}

# Declaramos un vector (@pautas) que contiene las 3 pautas de lectura posibles y llamamos a la funcion "orfs" (definida mas abajo) para cada pauta de lectura.

my @pautas = (0,1,2);

foreach my $k (@pautas) {
    &orfs($k);
}



# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # #  F U N C I O N E S  # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #



# Definimos la funcion "orfs". Esta funcion busca ORF's en la secuencia situados entre dos codones stop, ademas elimina el primer y ultimo ORF's incompletos y
# aquellos ORF's "vacios" que se generan cuando la funcion encuentra dos codones stop seguidos. Dentro de esta funcion hemos hecho una llamada a la siguiente
# funcion ("codonusage") que definimos mas abajo. Los resultados se muestran por pantalla en forma de tabla (una para cada pauta de lectura) y tambien son
# redireccionados a un archivo de salida (resultados.out).

sub orfs() {

    my $pauta = $_[0];                                       # Variable que contiene la pauta de lectura correspondiente.
    my @vecseq = split(//o,$seq);                            # Hacemos este split para conseguir un vector que tendra cada base de la secuencia en una posicion.
    my $b=$pauta;                                            # Posicion nucleotidica en cada momento dentro de la secuencia, inicialmente es igual a la pauta.
    my $ini=$b;                                              # Variable que indica el inicio de cada ORF.
    my $orfstr='';                                           # Variable que contiene los ORF's temporalmente, hasta que son almacenados en el hash.
    my $codon;                                               # Variable que lee la secuencia de DNA de triplete en triplete.
    my $outformat = "%-5s %8d  %8d  %8d  %8.2f       +\n";   # Variable que contiene el formato del "printf" de las lineas 147 y 148.
                  
                   #          1         2         3         4         5         6         7         8         9        10
                   # 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
    
    open (SALIDA,">> resultados.out");                                                                  # Redirigimos los resultados del programa almacenados
                                                                                                        # en un filehandle a un archivo de salida .out.
    print STDOUT    "Pauta   Inicio     Final  Longitud     Score    Strand\n";
    print STDOUT    "======================================================\n";
    print SALIDA    "Pauta   Inicio     Final  Longitud     Score    Strand\n";
    print SALIDA    "======================================================\n";

    while ($b < (scalar(@vecseq) - (3 + $pauta))) {
	
	$codon = $vecseq[$b].$vecseq[$b+1].$vecseq[$b+2];
	
	if ($codon ne "TGA" && $codon ne "TAA" && $codon ne "TAG") {                                    # Va concatenando los codones leidos en 
	    $orfstr = $orfstr.$codon;                                                                   # $orfstr hasta que encuentra un codon stop.
	    $b = $b+3;                                                                                  # Corre 3 posiciones dentro del ORF actual.
	}
	
	else {                                                                                          # Si encuentra un codon stop en la secuencia.
	    
	    if (length($orfstr) > 0) {                                                                  # Si el ORF actual no esta vacio lo almacena en
		push @{ $orfs{$pauta} }, [ $orfstr, $ini+1, $b, &codonusage($orfstr) ];                 # su posicion correspondiente dentro del hash "%orfs".
		$ini = $b + 3;                                                                          # Para saltar el codon stop y no incluirlo en el orf.
	    }
	    
	    $orfstr = '';                                                                               # Vaciamos $orfstr para que acumule el ORF siguiente.
	    $b = $b+3;                                                                                  # Corre 3 posiciones dentro de la secuencia completa.
	}
    }
    
    if (scalar(@{ $orfs{$pauta} }) > 0) {                                                               # Eliminamos el primer ORF incompleto por la izquierda
	shift @{ $orfs{$pauta} };                                                                       # haciendo un shift del hash %orfs para cada pauta.
	
	foreach my $s (@{ $orfs{$pauta} }) {                                                            # Imprimimos el contenido del hash, mostrando para cada
	    printf STDOUT $outformat, $pauta, $s->[1], $s->[2], length($s->[0]), $s->[3];               # ORF: pauta, origen, final, longitud y score obtenido en el
	    printf SALIDA $outformat, $pauta, $s->[1], $s->[2], length($s->[0]), $s->[3];               # calculo de la razon de verosimilitud.
	}
    }
    
    else {                                                                                              # Si alguna de las pautas no genera ningun ORF, imprime
	                                                                                                # este mensaje de aviso.
	printf STDOUT "%s : NO ENCONTRAMOS ORF's PARA ESTA PAUTA\n", $pauta;
	printf SALIDA "%s : NO ENCONTRAMOS ORF's PARA ESTA PAUTA\n", $pauta;
    }
    
    print STDOUT "\n";
    print SALIDA "\n";

    close(SALIDA);
}


# Definimos la funcion "codonusage". Esta funcion calcula la razon de verosimilitud en escala logaritmica para cada codon. Dicha razon es la probabilidad conjunta
# de los codones segun la tabla de uso de codones especificada al ejecutar el programa, respecto (dividido por) la probabilidad conjunta uniforme de los codones
# (la probabilidad uniforme de un unico codon es 1/64). La calculamos en escala logaritmica para poder manipular los numeros muy pequenos que se pueden generar
# al multiplicar muchas cifras con valores entre 0 y 1. El logaritmo del producto es la suma de logaritmos y el logaritmo de la division es la resta de logaritmos.
# Esta funcion esta incluida con una llamada dentro de la funcion "orfs", de modo que se ejecutara automaticamente para cada pauta de lectura.

sub codonusage() {

    my $seq = $_[0];                                                # Contiene la secuencia de DNA completa en un unico bloque.
    my @vcodons = ($seq =~ m/.../g);                                # Almacena los tripletes generados a partir de la secuencia en un vector (@vcodons).
    my $rao = 0;                                                    # Contendra el resultado del calculo de la razon de verosimilitud para cada codon.
    my $i = 0;
    
    while ($i < scalar(@vcodons)) {

	$rao = $rao + log($codonusage{$vcodons[$i]}) - log(1/64);   # Calculo de la razon de verosimilitud para cada codon.
	$i = $i + 1;
    }

    return sprintf("%.2f",$rao);                                    # Nos devuelve el resultado del calculo con un formato de 2 decimales. Al estar integrada
}                                                                   # esta funcion con la anterior ("orfs"), el score calculado para cada codon aparecera
                                                                    # integrado en la tabla de resultados.
