#!/usr/bin/perl -w

use strict;

my $contigstart;
my $start;
my $contiglength;
my $genome = "/mnt/NFS_UPF/bioinfo/BI/genomes/2018/Notechis_scutatus/genome.fa";
my $contig;
my $genomelength = "/mnt/NFS_UPF/bioinfo/BI/genomes/2018/Notechis_scutatus/genome.lengths";
my $contigend;
my $end;
my $output;
my $index = "/mnt/NFS_UPF/bioinfo/BI/genomes/2018/Notechis_scutatus/genome.index";
my $predict;
my $contigvell = "res"; #per diferenciar quan estiguem parlant del mateix contig que abans i per tant son exons o un contig amb un nom diferent#
my $length;

my $proteina = $ARGV [0]; #El primer argument que posem és el nom de la proteina que volem veure si està al genome de la nostra espècie#

my $directoryprot = "./proteines_Notechis_scutatus/$proteina";

chdir ("$directoryprot");

system("sed 's/\#//g' $proteina.aa.fa | sed 's/@//g' | sed 's/%//g' | sed 's/U/X/g' > $proteina.mod.aa.fa");#canviem tots els caracters que no ens interessen#
       
system("tblastn -query $proteina.mod.aa.fa -db $genome -outfmt 7 > $proteina.blast");

system ("cat $proteina.blast");

#Un cop tens això tu tries quins contigs vols utilitzar en funció del E-value, la llargada de l'alineament i el % de coincidència#

system ("mkdir tcoffee");
system ("mkdir exonerate");

print "enter contig name: ";

$contig = <STDIN>; 
chomp($contig);

while ($contig ne "fi"){

print "enter contig start point: ";

$contigstart = <STDIN>;
chomp($contigstart);

print "enter contig end point: ";

$contigend = <STDIN>;
chomp($contigend);

$contiglength = `grep $contig $genomelength | cut -f1 -d ' '`; 



if ($contigstart < $contigend){

$start = $contigstart - 50000;
$end = $contigend + 50000; 

if ($start < 0){

	$start = 1;
}
if ($end > $contiglength){

	$end = $contiglength;

}

$length = $end - $start;

}
else {
$start = $contigstart + 50000;
$end = $contigend - 50000; 
if ($end < 0){

	$end = 1;
}
if ($start > $contiglength){

	$start = $contiglength;
}

$length = $start - $end;

}

system("fastafetch $genome $index $contig > $contig.fa");

if ($contigstart < $contigend){

system("fastasubseq  $contig.fa $start $length > genomic.$contig.fa");
}

else {
system("fastasubseq  $contig.fa $end $length > genomic.$contig.fa");
}

system ("exonerate -m p2g --showtargetgff -q $proteina.mod.aa.fa -t genomic.$contig.fa -n 1 > exonerate.$proteina.$contig.gff 2>error_exonerate.$contig.txt");

system("egrep -w exon exonerate.$proteina.$contig.gff > $proteina.pred.$contig.gff");

if (-z "$proteina.pred.$contig.gff"){
print("Exonerate shows no results\n"); 
}

system("fastaseqfromGFF.pl genomic.$contig.fa $proteina.pred.$contig.gff > $contig.nuc.fa");

system ("fastatranslate $contig.nuc.fa -F 1 > $contig.pep.fa");

system ("sed s/*/X/g $contig.pep.fa > $contig.pepX.fa");

system ("t_coffee $contig.pepX.fa $proteina.mod.aa.fa > $contig.tcoffee.txt 2>error_tcoffee.$contig.txt");

system("mv $contig.tcoffee.txt ./tcoffee");
system("mv $proteina.pred.$contig.gff ./exonerate");


    
print "enter another contig name: ";
   
$contig = <STDIN>; 
chomp($contig);

}
chdir ("../..");

