#!/usr/bin/perl -w

use strict;

print "Enter the name of the protein: \n";
my $prot = <STDIN>;
chomp ($prot);

#We begin replacing the U of the protein to an X

system ("sed s/U/X/g $prot > $prot.x");
print "The U has been replaced successfully\n";

#Now we perform the tblastn

system ("tblastn -query $prot.x -db /mnt/NFS_UPF/soft/genomes/2019/Colinus_virginianus/genome.fa -out $prot.tblastn -outfmt 7 -evalue 0.001");
print "The tblastn was performed successfully and it is shown below\n";
system ("cat $prot.tblastn");

#Next we proceed to perform the fastafetch

print "Introduce the scaffold name: \n";
my $scaffold = <STDIN>;
chomp ($scaffold);
system ("fastafetch /mnt/NFS_UPF/soft/genomes/2019/Colinus_virginianus/genome.fa /mnt/NFS_UPF/soft/genomes/2019/Colinus_virginianus/genome.index $scaffold > $prot.fastafetch");
print "fastafetch was performed successfully\n";

#Now we do the fastasubseq

print "Introduce the position of the first nucleotide (in other words, the smallest number): \n";
my $start = <STDIN>;
chomp ($start);

print "Introduce the position of the last nucleotide (in other words, the biggest number): \n";
my $final = <STDIN>;
chomp ($final);
$start = $start - 50000;
if ($start <= 0){
    $start = 1;
}

my $slength; #This is the vector length
my $lines; #This are the lines in the fasta file
my $nucleotides; #This are the number of nucleotides in one row of the fasta file

$lines = `cat $prot.fastafetch | wc -l`;
$lines -=2;
$nucleotides = `head -n 2 $prot.fastafetch | tail -n 1 | wc -c`;
$nucleotides -=1;
$slength = $lines*$nucleotides; #The vector length will be the number of nucleotides in one row multiplyied by the number of rows
$final = $final + 50000;
if ($final > $slength){
    $final = $slength;
}

my $length = $final - $start; #This is the length we will use for the fastasubseq
system ("fastasubseq $prot.fastafetch $start $length > $prot.subseq");
print "The fastasubseq was performed successfully\n";

#Next we will perform the exonerate

system ("exonerate -m p2g --showtargetgff -q $prot.x -t $prot.subseq | egrep -w exon > $prot.exonerate.gff");
print "The exonerate was performed successfully\n";

#At this point we will perform the fastaseqfromgff

system ("fastaseqfromGFF.pl $prot.subseq $prot.exonerate.gff > $prot.pred.nuc");
print "The fastaseqfromGFF was performed successfully\n";

#Now we will perform the fastatranslate

system ("fastatranslate -F 1 $prot.pred.nuc > $prot.pred.aa");
print "The fastatranslate was performed successfully\n";

#We are almost there! We are going to replace the * to an X

system ("sed s/*/X/g $prot.pred.aa > $prot.pred.aa.x");
print "The * was replaced successfully\n";

#Finally, we do the tcoffee

system ("t_coffee $prot.x $prot.pred.aa.x > tcoffee_$prot");
system ("cat tcoffee_$prot");
print "The tcoffee was performed successfully. Congratulations, the analysis is finished!\n";
