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

my $p; #proteïna
my $seq; #fasta sequence


#Primer, el programa ens pregunta quina és la proteïna amb la que volem treballar, per accedir al query d'interès

print "Amb quina proteïna vols treballar?";
$p = <STDIN>;
chomp ($p);

#Es canvien les U per les X en el fitxer de text

system ("sed 's/U/X/g' $p.fa.txt > $p.fa");

print "S'han canviat les U per X\n";

print "Running blast…\n";

system ("tblastn -query $p.fa -db /cursos/20428/BI/genomes/2017/Astyanax_mexicanus/genome.fa -outfmt 7 -evalue 0.001 -out $p.blast");

print "Blast completat\n";

#El programa ens mostra per pantalla el resultat del blast. L'argument -outfmt ens donarà el blast simplificat (taula dels hits). És important tenir en compte, que el programa només ens mostra aquells hits que són etsadisticament significatius, és a dir, aquells amb un e-value <= 0.001 (argument -evalue 0.001).

system ("more $p.blast");

print "\n";

print "Inici de l'anàlisis:\n";


############## PROGRAMA ANALISIS ###################

my $s; #scaffold
my $i; #primer nucleòtid
my $f; #últim nucleòtid
my $l; #llargada del hit

print "Entrar l'scaffold: ";
$s = <STDIN>;
chomp ($s);

print "Primer nucleòtid (el més petit): ";
$i = <STDIN>;
chomp ($i);

print "Últim nucleòtid (el més gran): ";
$f = <STDIN>;
chomp ($f);


#El programa extreu l'scaffold d'interés amb el Fastafetch

print "Inici Fastafetch\n";

system ("fastafetch /cursos/20428/BI/genomes/2017/Astyanax_mexicanus/genome.fa  /cursos/20428/BI/genomes/2017/Astyanax_mexicanus/genome.index '$s' > $s.fetch.fa");

print  "Fastafetch completat\n";

#Ara el programa extendrà 50000 nucleòtids per davant i per darrere del nucleòtid d'inici i final. Després, calcularà la llargada de l'scaffold.

if ($i >= 50000)
{$i=$i-50000;
 $f=$f+50000;}
else
{$i=0;};

$l=$f-$i;

#El programa realitza el Fastasubseq per tallar l'scaffold en les posicions d'interès

print "Inici Fastasubseq\n";

system ("fastasubseq $s.fetch.fa $i $l > $p.$s.subseq.fa"); 

print "Fastasubseq completat\n";


#El programa realitza l'Exonerate i seguidament l'Egrep, compara el query amb el fastasubseq i selecciona només els exons i els fusiona en un sol arxiu

print "Inici de l'Exonerate i Egrep\n";

system ("exonerate --exhaustive yes -m p2g --showtargetgff -q $p.fa -t $p.$s.subseq.fa > $p.$s.exonerate"); 

system ("egrep -w exon $p.$s.exonerate > $p.$s.exons.exonerate");

print  "Exonerate i Egrep completats\n" ;


#El programa obté la sequencia de cDNA generada anteriorment amb formar gff. Utilitzem el FastaseqfromGFF

print "Inici FastaseqfromGFF\n";

system ("fastaseqfromGFF.pl $p.$s.subseq.fa $p.$s.exons.exonerate > $p.$s.cDNA.fa");

print "FastaseqfromGFF completat\n" ;


#El programa tradueix el cDNA obtingut a proteïna amb el Fastatranslate

system ("fastatranslate -f $p.$s.cDNA.fa -F 1 > $p.$s.prot.fa");

print "Fastatranslate completat\n" ;


#El programa canvia els * per X del fitxer .prot.fa obtingut de el Fastatranslate.

system ("sed 's/*/X/g' $p.$s.prot.fa > $p.$s.protX.fa");

print "S'han canviat els * per X\n";


#El programa compara la predicció de proteina obtinguda amb la proteïna inicial de Danio rerio extreta del SelenoDB. Utilitza T-Coffee

print "Inici T-Coffee\n";

system ("t_coffee $p.fa $p.$s.protX.fa > $p.$s.tcoffee");

print "T-Coffee completat\n";


#El programa ens mostra el resultat del T-Coffee

system ("more $p.$s.tcoffee");


#El programa crea una carpeta anomenada $p.$s i mou tots els arxius creats en l'anàlisis dins d'aquesta carpeta.

system ("mkdir $p.$s");

system ("mv $s.fetch.fa $p.$s");
system ("mv $p.$s.subseq.fa $p.$s");
system ("mv $p.$s.exonerate $p.$s");
system ("mv $p.$s.exons.exonerate $p.$s");
system ("mv $p.$s.cDNA.fa $p.$s");
system ("mv $p.$s.prot.fa $p.$s");
system ("mv $p.$s.protX.fa $p.$s"); 
system ("mv $p.$s.tcoffee $p.$s");

print "Enhorabona super! Has acabat l'anàlisis!\n";