#!/usr/bin/perl -w

use strict;

use IPC::System::Simple qw(capture);

my $i;
my $l;
my @array;
my $scaffold;
my $start;
my $length;

@array = ("sel15", "eEFsec", "SELENOE","GPx1a","GPx1b","GPx2","GPx3a","GPx3b","GPx4a", "GPx4b", "GPx7", "GPx8", "DIO1", "DIO2","DIO3a", "DIO3b", "MsrA", "PSTK", "SBP2", "SecS", "SEPHS2", "SELENOH", "SELENOI", "SELENOJ", "SELENOK", "SELENOL", "SELENOM", "SELENON", "SELENOO1", "SELENOO2", "SELENOP1", "SELENOP2", "MSRB1a", "MSRB1b", "MSRB2", "MSRB3", "SELENOS", "SELENOT1", "SELENOT1b", "SELENOT2", "SELENOU1a", "SELENOU2", "SELENOU3", "SELENOW2", "TXNRD2", "TXNRD3", "SECp43");


$l = scalar @array;
$i = 0;

# Creem totes les carpetes

system ("mkdir -p output.blast output.fastafetch output.fastasubseq output.exonerate output.fastaseqfromgff output.fastatranslate output.tcoffee output.genewise");
print "Carpetes creades amb èxit\n";
#
#
#

while ($i<$l) {
print "\n";
print "\n";
print "COMENCEM AMB LA PROTEINA $array[$i]\n";
print "Ara farem el tblastn\n";
# print "$array[$i]\n";
system ("tblastn -evalue 0.005 -query input.aa/$array[$i].aa.fa -db /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Xiphophorus_hellerii/genome.fa -out output.blast/$array[$i].blast -outfmt 6");
print "tblastn de la selenoproteina $array[$i] DONE\n";
#
#
#
print "Ara farem el fastafetch\n";
$scaffold = `cut -f 2 output.blast/$array[$i].blast | head -n 1 | tail -n 1`;
chomp $scaffold;
print "Scaffold triat: $scaffold\n";
    
system ("fastafetch /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Xiphophorus_hellerii/genome.fa /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Xiphophorus_hellerii/genome.index $scaffold > output.fastafetch/$array[$i]xh.contig.fa");
print "fastafetch de l'array $array[$i] DONE\n";
#
#
#
print "Ara farem el fastasubseq\n";


# Definim la variable $start: la petita -50.000 (amb un if columna x < columna y; else {agafa columna y i resta-li 50.000).
# Definim una variable $length que serà 100000 sempre.


my $columna9;
my $columna10;

$columna9 = `cut -f 9 output.blast/$array[$i].blast | head -n 1 | tail -n 1`;
chomp $columna9;
$columna10 = `cut -f 10 output.blast/$array[$i].blast | head -n 1 | tail -n 1`;
chomp $columna10;
print "columna 9: $columna9\n";
print "columna 10: $columna10\n";

# Definim quina columna és la més petita

if ($columna9 lt $columna10) { $start = $columna9 - 50000;}
else {$start = $columna10 - 50000;}
print "start: $start\n";

# Definim la llargada de la variable

$length = 100000;

#
#
#

# Fem el fastasubseq

system ("fastasubseq output.fastafetch/$array[$i]xh.contig.fa $start $length > output.fastasubseq/$array[$i]xh.hit.fa");
print "fastasubseq de la selenoproteina $array[$i] DONE\n";
#
#
# Hem provat de fer el genewise i ens ha donat el mateix que l'exonerate. Donat que l'exonerate és més acurat, ens basarem en els resultats de l'exonerate.
# print "Ara fem el genewise\n";
# system ("genewise -pep -pretty -cdna -gff -both input.aa/$array[$i].aa.fa output.fastasubseq/$array[$i]xh.hit.fa > output.genewise/$array[$i]xh.genewise.fa");
# print "Genewise de la selenoproteina $array[$i] DONE\n";
#
#
print "Extraiem els exons de la nostra selenoproteïna amb l’exonerate\n";
system ("exonerate -m p2g --showtargetgff -q input.aa/$array[$i].aa.fa -t output.fastasubseq/$array[$i]xh.hit.fa -n 1 | egrep -w exon > output.exonerate/$array[$i]xh.exonerate.gff");
print "exonerate de la selenoproteina $array[$i] DONE\n";
#
#
#
print "Fem la predicció amb el fastaseqfromGFF\n";
system ("fastaseqfromGFF.pl output.fastasubseq/$array[$i]xh.hit.fa output.exonerate/$array[$i]xh.exonerate.gff > output.fastaseqfromgff/$array[$i]xh.pred.nuc.fa");
print "fastasubseqfromGFF de la selenoproteina $array[$i] DONE\n";
#
#
#
print "Traduim la sequencia a proteina amb el fastatranslate\n";
system ("fastatranslate -f output.fastaseqfromgff/$array[$i]xh.pred.nuc.fa -F 1 > output.fastatranslate/$array[$i]xh.pred.aa.fa");
print "fastatranslate de la selenoproteina $array[$i] DONE\n";
#
#
# Substituim els * per X
system ("sed -i 's/*/X/g' output.fastatranslate/$array[$i]xh.pred.aa.fa");
#

print "Comparem les proteïnes amb el Tcoffee\n";
system ("t_coffee input.aa/$array[$i].aa.fa output.fastatranslate/$array[$i]xh.pred.aa.fa > output.tcoffee/$array[$i]xh.result.fa 2>/dev/null");
print "Tcoffee de la selenoproteina $array[$i] DONE\n";
#
#
#
# Passem a la seguent proteina 
$i = $i+1;
}
