#!/usr/bin/perl -w

use strict;


### VARIABLES ###

my $seqprob;           #se li associa el nom del fitxer amb seqüència problema
my $out;               #se li associa el nom del fitxer de sortida
my $ref;               #llegeix la referència de la seqüència en fasta
my @fasta;             #llegeix les altres linies de la seqüència en fasta
my $fastap;            #posició del vector @fasta
my $seq;               #conté la seqüència llegida a @fasta sense salts de línia
my $seq1;              #seqüència en la segona pauta de lectura
my $seq2;              #seqüència en la tercera pauta de lectura
my $seq3;              #seqüència en la quarta pauta de lectura
my $seq4;              #seqüència en la cinquena pauta de lectura
my $seq5;              #seqüència en la sisena pauta de lectura
my $long;              #llargada mínima del vector
my $complem;           #si es té en compte la cadena complementaria
my @orf1;              #conté tots els orf pel primer marc de lectura
my @orf2;              #conté tots els orf pel segon marc de lectura
my @orf3;              #conté tots els orf pel tercer marc de lectura
my @orf4;              #conté tots els orf pel quart marc de lectura
my @orf5;              #conté tots els orf pel cinquè marc de lectura
my @orf6;              #conté tots els orf pel sisè marc de lectura
my @seqv;              #@seqv conté $seq codó a codó; $inici conté el codó d'inici de l'orf
my (@orftot,@atgPos, @stopPos);  #@orftot conte tots els orf d'una seqüència; @atgPos conté la posició on comença cada orf
my ($orftot,$atgPos, $stopPos);
my @tradv;             #@tradv conte els orf traduïts
my $tradv;
my $longseq;           #longitud de la seqüència
my $inici;             #guarda la opció que l'orf comenci o no amb Met
my $char;              #caràcters per línia a l'imprimir els orf
my %aminoacid=("TTT" => "F",  #conté la relació codó-aminoàcid (codi d'una lletra)
"TTC" => "F",
"TTA" => "L",
"TTG" => "L",
"CTT" => "L",
"CTC" => "L",
"CTA" => "L",
"CTG" => "L",
"ATT" => "I",
"ATC" => "I",
"ATA" => "I",
"ATG" => "M",
"GTT" => "V",
"GTC" => "V",
"GTA" => "V",
"GTG" => "V",
"TCT" => "S",
"TCC" => "S",
"TCA" => "S",
"TCG" => "S",
"CCT" => "P",
"CCC" => "P",
"CCA" => "P",
"CCG" => "P",
"ACT" => "T",
"ACC" => "T",
"ACA" => "T",
"ACG" => "T",
"GCT" => "A",
"GCC" => "A",
"GCA" => "A",
"GCG" => "A",
"TAT" => "Y",
"TAC" => "Y",
"TAA" => "*",
"TAG" => "*",
"CAT" => "H",
"CAC" => "H",
"CAA" => "Q",
"CAG" => "Q",
"AAT" => "N",
"AAC" => "N",
"AAA" => "K",
"AAG" => "K",
"GAT" => "D",
"GAC" => "D",
"GAA" => "E",
"GAG" => "E",
"TGT" => "C",
"TGC" => "C",
"TGA" => "*",
"TGG" => "W",
"CGT" => "R",
"CGC" => "R",
"CGA" => "R",
"CGG" => "R",
"AGT" => "S",
"AGC" => "S",
"AGA" => "R",
"AGG" => "R",
"GGT" => "G",
"GGC" => "G",
"GGA" => "G",
"GGG" => "G");


### FUNCIONS ####

#busca els orf que comencen per ATG dins una seqüència
#rep la seqüència ("seqv)
#retorna el vector que conté els orf i els que conten les posicions d'inici i final de de cada orf
sub orfExtract
{
    my $gravar=0;               #variable que determina si es grava (=1) o no (=0) el vector
    my ($args) = @_;            #conté els arguments que se li passen a la funció
    my @args = @$args;
    my @seqv=@args;             #conté la seqüència a analitzar
    my $k=0;                    #posició del vector @orftot
    my $x=0;                    #posició del vector @seqv
    my $i=0;                    #posició del vector @atgPos
    my $j=0;                    #posició del vector @orf
    my $l=0;                    #posició del vector @stopPos
    my @orf=();                 #conté un orf
    my @orftot = ();            #conté el conjunt d'orf d'una pauta de lectura
    my @atgPos =();             #conté la posició d'inici de cada orf
    my @stopPos = ();           #conté la posició de final de cada orf
    
    for($x = 0; $x < scalar(@seqv); $x++)    #recorre la seqüència codó a codó
    {
	
	if ($seqv [$x] eq 'ATG')             #si el codó és ATG (Met, inici)
	{
	    
	    if ($gravar == 0)                #i si $gravar és =0 (per tant, no es té un compte un ATG a meitat de l'orf)
	    {
		$atgPos[$i] = ($x+1)*3-2;    #grava la posició d'inici de l'orf
		$i++;                        #avança una posició a @atgPos
		$gravar = 1;                 #iguala $gravar a 1
		$j = 0;
		@orf = ();                   #buida el vector orf
	    }
	
	}
	
	if ( ($seqv[$x] eq ('TAA') || $seqv[$x] eq ('TAG') || $seqv[$x] eq ('TGA')) && ($gravar == 1))    #si el codó llegit és TAA, TAG o ATG (stop) i $gravar és 1
	{ 
	    $orf[$j] = $seqv[$x];            #copia la posició al vector orf
	    $stopPos[$l] = ($x+1)*3;         #grava la posició de final de l'orf
	    $l++;                            #avança una posició a @stopPos
	    $gravar = 0;                     #iguala $gravar a 0
	    
	    if( scalar(@orf) >= $long)       #si la longitus de @orf és major que $long
	    {
		$orftot[$k] = [@orf];        #grava el vector @orf a la posició pertinent del vector @orftotal
		$k++; 
	    }	    
	    else                             #si la longitus de @orf és menor que $long 
	    {  
		pop(@atgPos);                #elimina la darrera posició del vector @atgPos
		$i--;                        #tira el comptador enrera una posició
		pop(@stopPos);               #elimina la darrera posició del vector @stopPos
		$l--;                        #tira el comptador enrera una posició
	    } 
	
	} 
	
	if($gravar == 1)                     #si el valor de $gravar és 1
	{
	    $orf[$j] = $seqv [$x];           #va copiant les posicions de $seqv al vector @orf
	    $j++;
	}
 
    } 
    
    if ($gravar == 1)                        #si $gravar és 1, és a dir, s'ha trobat una Met però no un Stop
    {
	pop(@atgPos);                        #elimina la darrera posició del vector @atgPos
	$i--;                                #tira el comptador enrera una posició	
    }
    return (\@{orftot},\@atgPos,\@stopPos);  #retorna el valors @orftot, @atgPos i @stopPos

}


#busca els orf sense que comencin per ATG (és a dir, d'inici a Stop o d'Stop a Stop) dins una seqüència
#rep la seqüència ("seqv)
#retorna el vector que conté els orf i els que conten les posicions d'inici i final de de cada orf
sub orfExtractNoATG
{

 my ($args, $pauta) = @_;               #conté els arguments que se li passen a la funció
    my @args = @$args;
    my @seqv=@args;             #conté la seqüència a analitzar
    my $x=0;                    #posició del vector @seqv
    my $i=0;                    #posició del vector @atgPos
    my $j=0;                    #posició del vector @orf
    my $l=0;                    #posició del vector @stopPos
    my $c=1;                    #comptado pel salt de línia
    my @orf=();                 #conté un orf
    my $atgPos;                 #conté la posició d'inici de cada orf
    my $stopPos;                #conté la posició de final de cada orf

    $atgPos = 1;                             #grava la posició de la primera 'a' de 'atg' en la seqüència
    $i++;                                    #avança una posició a @atgPos

    for($x = 0; $x < scalar(@seqv); $x++)    #recorre la seqüència codó a codó
    {
        $orf[$j] = $aminoacid{$seqv[$x]};    #va copiant les posicions de $seqv al vector @orf
        $j++;
       
	if ($seqv[$x] eq ('TAA') || $seqv[$x] eq ('TAG') || $seqv[$x] eq ('TGA'))   #si el codó llegit és TAA, TAG o TGA (stop)     
        {
	    if (scalar (@orf) >= $long)      #si l'orf té la longitud mínima
	    {
		$stopPos = ($x+1)*3;         #grava la posició de final de l'orf 

		print OUT "$ref$pauta$atgPos-$stopPos\n";   #imprimeix en fitxer la referència completa
		$c = 1;
            
		for ($i=0; $i< scalar (@orf);$i++)           #mentre $j=0 i fins que finalitzi el vector @tradv
		{
		    print OUT "$orf[$i]";    #imprimeix en fitxer la seqüència d'aminoàcids
		    
		    if ($c % $char == 0 )    #si $c és múltiple de 80
		    {
			print OUT "\n";      #salt de línia
		    }
		    
		    $c++;
		}
		print OUT "\n\n";            #salt de línia doble
	    }

	    @orf=();                         #buida el vector @orf per gravar-hi el següent orf
	    $j = 0;                          #j torna a ser la priemra posició del vector
	    $atgPos = ($x+1+1)*3-2;          #grava la posició d'inici de l'orf	    
            
        }

    }

}


# tradueix la seqüència de DNA a proteïna
# rep el vector de vectors que conté tots els orf ("orftot")
# retorna un vector de vectors amb els orf traduïts (@tradv)
sub traduccio
{
    my ($orfv) = @_;            #conté els arguments que se li passen a la funció
    my @orfv = @$orfv;
    my @tradv = ();             #vector que conté tots els orf traduïts
    my @trad = ();              #vector on s'emmagatzema cada orf mentre es tradueix
    my $i=0;                    #conté les posicions de la primera dimensió del vector @orfv
    my $j=0;                    #conté les posicions de la segona dimensió del vector @orfv
    
    for $i ( 0 .. $#orfv )                   #mentre $i=0 i fins que finalitzi el vector @orfv
	{
	    @trad = ();                      #iguala el vector #@trad a 0

	    for $j ( 0 .. $#{$orfv[$i]} )                #mentre $j=0 i fins que finalitzi el vector @orfVec
	    {
		$trad[$j] = $aminoacid{$orfv[$i][$j]};   #aplica el hash a cada codó  de @orfv i el guarda a @trad
	    }
	    
	    $tradv[$i] = [@trad];            #grava el vector @trad al vector de vectors @tradv

	}

    return \@{tradv};                        #retorna @tradv

}


#imprimeix els orf, incloent la referència, la pauta de lectura i la posició del codó d'inici
#rep el vector amb l'orf traduït (tradv), els vectors amb la posicions d'inici i final de l'orf (@atgPos i @stopPos) i la pauta de lectura
sub mostrar
{
    my ($tradv,$atgPos,$stopPos,$pauta) = @_;  #conté els arguments que se li passen a la funció
    my @atgPos = @$atgPos;
    my @stopPos = @$stopPos; 
    my @tradv = @$tradv;              
    my $i=0;                          #conté les posicions de la segona dimensió del vector @atgPos
    my $l=0;                          #conté les posicions de la segona dimensió del vector @stopPos
    my $j=0;                          #conté les posicions de la primera dimensió del vector @orfv
    my $c=1;                          #comptador pel salt de línia
	
    for $i ( 0 .. $#tradv )                   #mentre $i=0 i fins que finalitzi el vector @tradv
    {
	print OUT "$ref$pauta $atgPos[$i]-$stopPos[$l]\n"; #imprimeix en fitxer la referència completa
	$l++;
	$c = 1;
	
	for $j ( 0 .. $#{$tradv[$i]} )        #mentre $j=0 i fins que finalitzi el vector @tradv
	{

	    print OUT $tradv[$i][$j] ;        #imprimeix en fitxer la seqüència d'aminoàcids
	    if ($c % $char == 0 )                #si $c és múltiple de 80
	    {
		print OUT "\n";               #salt de línia
	    }
	    $c++;
	}

	print OUT "\n\n";                     #salt de línia doble

    }
	
}

#ajusta les posicions, ja que les posicions preses considera el primer nucleòtid de la pauta com el primer de la seqüència. Tambe fa la correccio per la cadena complementaria
#rep els vectors amb les posicions d'inici i final, el corrector específic de pauta ($num) i la opció de cadena complementaria ($comp)
#retorna  els vectors amb les posicions d'inici i final, ja modificats
sub AjustaPos
{
    my ($atgPos,$stopPos,$num,$comp) = @_;  #contŽé els arguments que se li passen a la funciŽó#
    my @atgPos = @$atgPos;
    my @stopPos = @$stopPos;               
    my $i=0;                          #contŽé les posicions de la segona dimensiŽó del vector @atgPos

    #print  scalar(@atgPos). " \n";
    #print  scalar(@stopPos). " \n";
    for $i ( 0 .. $#atgPos )
    {

	if ($comp eq 'y')
	{
	    $atgPos[$i] = $longseq - $atgPos[$i] + $num;
	    $stopPos[$i] = $longseq - $stopPos[$i] + $num;
	}
	else
	{
	    $atgPos[$i] = $atgPos[$i] + $num;
	    $stopPos[$i] = $stopPos[$i] + $num;
	}
    }
    
    return (\@atgPos,\@stopPos);
}


### COS PRINCIPAL DEL PROGRAMA ###

#presentació

print "\n\t#########################################################\n";
print "\t#                                                       #\n";
print "\t#                                                       #\n";
print "\t#                                                       #\n";
print "\t#                      FINDORF 2003                     #\n";
print "\t#                                                       #\n";
print "\t#                       versió 1.0                      #\n";
print "\t#                                                       #\n";
print "\t#                                                       #\n";
print "\t#          Ramírez-Soriano A - Molina-Tomàs MC          #\n";
print "\t#                                                       #\n";
print "\t#                                                       #\n";
print "\t#                                                       #\n";
print "\t#########################################################\n\n";


#introduir el nom del fitxer de la seqüència problema
 
print "introduzca el nombre del archivo que contiene la sequencia en la que desea buscar los orf:\n"; 

$seqprob = <STDIN>; 

chomp ($seqprob); 


#introduir el nom del fitxer de sortida 

print "introduzca el nombre del arxivo de salida:\n";

$out = <STDIN>; 

chomp ($out);


#definir codó d'inici

print "desea que el orf empieze por Met? (y/n):\n";

$inici = <STDIN>;

chomp ($inici);


#fixar una llargada mínima de l'orf

print "indique el número mínimo de codones que debe tener un orf para ser considerado significativo:\n";
 
$long = <STDIN>;


#tenir o no en compte la cadena complementaria

print "desea buscar orf tambien en la cadena complementaria? (y/n):\n";

$complem = <STDIN>;

chomp ($complem);


#caràcters per línia en fasta

print "quantos aminoacidos desea que aparezcan por línea al fichero de salida para los orf?\n";

$char = <STDIN>;

chomp ($char);


#obrir el fitxer amb la seqüència problema i associar-la a la variable $seq

open(FITXER,"<$seqprob") || die "no es pot obrir fitxer d'entrada\n";


#obrir el fitxer amb la seqüència problema

$ref = <FITXER>;        #llegeix la referència del fitxer (l�ínia 1)
chomp($ref);
 
@fasta = <FITXER>;      #llegeix la resta de línies del fitxer 

$seq = "";              #s'associa la variable a paraula buida 


#obté el contingut de la seqüència en fasta (ja se li eliminat la referència) en una sola línia

foreach $fastap (@fasta)        #per cada posició del vector @fasta 
{ 
    chomp($fastap );            #elimina el salt de línia 
    $seq = $seq . $fastap;      #concatena una línia al contingut de la variable seqüència
}

 
close(FITXER);                  #tancar el fitxer

$longseq = length($seq);

$seq = "\U$seq";                #passar totes les bases a majuscules 

$seq1 = substr($seq,1);         #elimina la primera posició de seq, obtenint la segona pauta de lectura

$seq2 = substr($seq,2);         #elimina les dues primeres posicions de seq, obtenint la tercera pauta de lectura




open(OUT, ">$out") || die "no es pot obrir el fitxer de sortida\n";    #obrir el fitxer de sortida


#busca els orf en cada pauta de lectura

@seqv = ($seq =~ m/.../g);                                 #guarda $seq al vector @seqv codó a codó
if($inici eq 'y')                                          #si volem que la seqüència comenci amb ATG
{
    ($orftot,$atgPos,$stopPos) = &orfExtract(\@{seqv});    #crida a la funció orfExtract, que extreu els orf
    @orftot = @$orftot;
    @atgPos = @$atgPos;
    @stopPos = @$stopPos;
    $tradv = &traduccio(\@{orftot});                       #crida a la funció traduccio, que tradueix a proteïna
    @tradv = @$tradv;
    &mostrar((\@{tradv},\@atgPos,\@stopPos,''));           #crida a la funció mostrar, que imprimeix els resultats en format fasta
}
else                                                       #si no
{
    &orfExtractNoATG ((\@{seqv}), '');                     #crida a la funció orfExtractNoATG, que executa tot el procés d'extracció i mostrat dels orf
}

@seqv = ($seq1 =~ m/.../g);
if($inici eq 'y')
{
    ($orftot,$atgPos,$stopPos) = &orfExtract(\@{seqv});
    @orftot = @$orftot;
    @atgPos = @$atgPos;
    @stopPos = @$stopPos;
    ($atgPos,$stopPos) = &AjustaPos(\@atgPos,\@stopPos,1,'n');    #crida la funció AjustaPos, que ajusta la posició segons la pauta
    @atgPos = @$atgPos;
    @stopPos = @$stopPos;
    $tradv = &traduccio(\@{orftot});
    @tradv = @$tradv;
    &mostrar((\@{tradv},\@atgPos,\@stopPos,''));
}
else
{
    &orfExtractNoATG ((\@{seqv}), '');
}

@seqv = ($seq2 =~ m/.../g);
if($inici eq 'y')
{
    ($orftot,$atgPos,$stopPos) = &orfExtract(\@{seqv});
    @orftot = @$orftot;
    @atgPos = @$atgPos;
    @stopPos = @$stopPos;
    ($atgPos,$stopPos) = &AjustaPos(\@atgPos,\@stopPos, 2,'n');
    @stopPos = @$stopPos;
    @atgPos = @$atgPos;
    $tradv = &traduccio(\@{orftot});
    @tradv = @$tradv;
    &mostrar((\@{tradv},\@atgPos,\@stopPos,''));
    @orf3 = @tradv;
}
else
{
    &orfExtractNoATG ((\@{seqv}), '');
}

if ($complem eq 'y')
{
    $seq3 = reverse($seq);               #gira la seqüència
    $seq3 =~ tr/acgtACGT/tgcaTGCA/;      #fa el complementari de la seqüència revertida: quarta pauta de lectura

    $seq4 = substr($seq3,1);             #elimina la primera posició de seq3, obtenint la cinquena pauta de lectura

    $seq5 = substr($seq3,2);             #elimina les dues primeres posicions de seq3, obtenint la sisena pauta de lectura


    @seqv = ($seq3 =~ m/.../g);
    if($inici eq 'y')
    {
	($orftot,$atgPos,$stopPos) = &orfExtract(\@{seqv});
	@orftot = @$orftot;
	@atgPos = @$atgPos;
	@stopPos = @$stopPos;
	($atgPos,$stopPos) = &AjustaPos(\@atgPos,\@stopPos, 1,'y');
	@stopPos = @$stopPos;
	@atgPos = @$atgPos;
	$tradv = &traduccio(\@{orftot});
	@tradv = @$tradv;
	&mostrar((\@{tradv},\@atgPos,\@stopPos,'complem'));
	@orf4 = @tradv;
    }
    else
    { 
	&orfExtractNoATG ((\@{seqv}), 'complem');
    }
	
    @seqv = ($seq4 =~ m/.../g);
    if($inici eq 'y')
    {
	($orftot,$atgPos,$stopPos) = &orfExtract(\@{seqv});
	@orftot = @$orftot;
	@atgPos = @$atgPos;
	@stopPos = @$stopPos;
	($atgPos,$stopPos) = &AjustaPos(\@atgPos,\@stopPos, 0,'y');
	@stopPos = @$stopPos;
	@atgPos = @$atgPos;
	$tradv = &traduccio(\@{orftot});
	@tradv = @$tradv;
	&mostrar((\@{tradv},\@atgPos,\@stopPos,'complem'));
	@orf5 = @tradv;
    }
    else
    {
	&orfExtractNoATG ((\@{seqv}), 'complem');
    }

    @seqv = ($seq5 =~ m/.../g);
    if($inici eq 'y')
    {
	($orftot,$atgPos,$stopPos) = &orfExtract(\@{seqv});
	@orftot = @$orftot;
	@atgPos = @$atgPos;
	@stopPos = @$stopPos;
	($atgPos,$stopPos) = &AjustaPos(\@atgPos,\@stopPos, -1,'y');
	@stopPos = @$stopPos;
	@atgPos = @$atgPos;
	$tradv = &traduccio(\@{orftot});
	@tradv = @$tradv;
	&mostrar((\@{tradv},\@atgPos,\@stopPos,'complem'));
	@orf6 = @tradv;
    }
    else
    {
	&orfExtractNoATG ((\@{seqv}), 'complem');
    }
}

close(OUT);