#!/usr/bin/perl -w

use strict;

my $p; #Bos taurus protein

use Cwd qw (cwd);
my $dir = cwd; #program.pl saved in directory 

print "Enter Bos Taurus (Cow) protein name:"; #The program asks for the name of protein
$p = <STDIN>; #p variable will be the inserted protein name in the previous step
chomp ($p);   #s'ha de posar

#First, we create a directory to save the obtained results

system("mkdir -p $dir/$p"); #Then we create the directory 
my $directory; #falta fer que si ja existeix el directory no ens el creï
$directory = "$dir/$p"; 
print("$p directory did not exist and was created\n"); #And print that the directory did not exist and was just created


#PREPROCESSING 
#In the next step, the program changes the U for an X in the Cow genome selenoproteins and will also clean the sequences and eliminate the @, # or % that are found between the sequences.  

system ("sed 's/U/X/g' $dir/$p.raw.fa | sed 's/#//g' | sed 's/@//g' | sed 's/%//g' > $directory/$p.fa");
print "U changed to X and $p sequence is clean to proceed\n";
print "\nTBLASTN in process\n"; 

#Now we perform the tblastn
#We perform a tblastn analysis, which compares and aligns the protein query (Bos taurus) to our subject genome (Nanger dama)

system ("tblastn -query $directory/$p.fa -db /mnt/NFS_UPF/soft/genomes/2021/Nanger_dama/genome.fa -out $directory/$p.tblastn -outfmt 7 -evalue 0.001"); #Discard E-values higher than 0.001 outfmt is used to "posar més maco el resultat"
print "tblastn for $p was performed successfully and it is shown below\n";
system ("cat $directory/$p.tblastn"); #Based on the e-value and sequence identity, we select the best hits (scaffolds) manually. 

#FASTAFETCH
#Extraction of the chosen scaffold from the Nanger dama genome 

print "Introduce the scaffold name: \n"; 
my $s = <STDIN>; #this step is done manually
chomp ($s);
system ("fastafetch /mnt/NFS_UPF/soft/genomes/2021/Nanger_dama/genome.fa /mnt/NFS_UPF/soft/genomes/2021/Nanger_dama/genome.index $s > $directory/$p.$s.fastafetch.fa");
print "fastafetch was performed successfully\n";

#FASTASUBSEQ
#We will extract the region of interest from Nanger dama genome

print "Introduce the position of the first nucleotide: \n";
my $start = <STDIN>; #First nucleotide of the best hit of the scaffold. 
chomp ($start);

print "Introduce the position of the last nucleotide: \n";
my $final = <STDIN>; #Last nucleotide of the best hit of the scaffold. 
chomp ($final);
$start = $start - 50000; #We will add 50000bp to the lenght 
if ($start <= 0){ #If the scaffold ends, then we will start at the 1st nucleotide
    $start = 1;
}

my $slength; #This is the lenght of the selected scaffold
my $lines; #This are the lines in the fasta file
my $nuc; #This are the number of nucleotides in one row of the fasta file

$lines = `cat $directory/$p.$s.fastafetch.fa | wc -l`;
$nuc = `head -n 2 $directory/$p.$s.fastafetch.fa | tail -n 1 | wc -c`;
$slength = $lines*$nuc; #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 $directory/$p.$s.fastafetch.fa $start $length > $directory/$p.$s.subseq");
print "The fastasubseq was performed successfully\n";

#EXONERATE+EGREP
#It will store the exons of the protein of interest

system ("exonerate -m p2g --showtargetgff --exhaustive yes -q $directory/$p.fa -t $directory/$p.$s.subseq | egrep -w exon > $directory/$p.$s.exonerate.gff");
print "The exonerate was performed successfully\n";

#FASTASEQFROMGFF
#The program will predict a cDNA sequence 

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

#FASTATRANSLATE
#It will translate the cDNA predicted sequence from Fastaseqfromgff into a protein sequence

system ("fastatranslate -F 1 $directory/$p.$s.pred.nuc > $directory/$p.$s.pred.fa");
print "The fastatranslate was performed successfully\n";

#Change * to X:

system ("sed 's/*/X/g' $directory/$p.$s.pred.fa > $directory/$p.$s.predx.fa");
print "The * was replaced successfully by an X\n";
system ("cat $directory/$p.$s.predx.fa");

#TCOFFEE
#Finally, we will align the predicted protein from Fastatranslate (Nanger dama) and the query sequence (Bos taurus) by using Tcoffee

system ("t_coffee $directory/$p.fa $directory/$p.$s.predx.fa > $directory/$p.$s.tcoffee");
system ("cat $directory/$p.$s.tcoffee");
print "Tcoffee was performed successfully!\n";
