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

### Remember! First of all you must introduce in the shell the following commands:
### sudo mkdir /mnt/public
### sudo mount.cifs -o username=u124040 //fs-aules.ds.upf.edu/PUBLIC/20428 /mnt/public
### export PATH=/mnt/NFS_UPF/bioinfo/BI/bin:$PATH
### export PATH=/mnt/NFS_UPF/bioinfo/BI/soft/genewise/x86_64/bin:$PATH
### export WISECONFIGDIR=/mnt/NFS_UPF/bioinfo/BI/soft/genewise/x86_64/wise2.2.0/wisecfg/

my $q; #protein we will use as query, obtained from SelenoDB 1.0
my $u; #the content of the file that contains the query with U that we need to convert into an X
my $s; #scaffold chosen from the Blast we already performed, where we found the best hits
my $start; #first nucleotide of the hits in the same scaffold (lower number)
my $final; #last nucleotide of the hits in the same scaffold (higher number)
my $length; #hit lenght
my $k; #the content of the file ending with .p.fa
my $rev; #the answer to the question about if our hit is reverse
my $n; #the content of the file ending with .genomic.fa, obtained with fastasubseq
my $runagain; #the answer to the question about if we want to run the program again

### First the program asks for the name of the protein we want to work with, obtained from SelenoDB 1.0 from the human genome, the genome closer to the Pseudonaja textilis. We will use this protein as the query, and the program will name the files accordingly ### 

print "Remember! First of all you must introduce in the shell the proper commands!\n";

print "Enter the query name you downloaded from SelenoDB 1.0: \n";
$q = <STDIN>;
chomp ($q);

### First of all, the program changes the U into an X in the query file ###

open(QUERY, "$q");
undef $/;
$u = <QUERY>;
close QUERY;
$/ ="\n";

$u =~ s/U/X/g;

open(XQUERY, '>', "$q"); #we save it in the same file
print XQUERY $u;
close XQUERY;

print "In the query file, the U has been changed into an X! ;)\n";

print "Running Blast! \n";

### The program will find the region where the gene of this protein is located using tblastn. This blast will automatically filter the hits to show only those with an e-value under 0.01 ###

system ("tblastn -query $q -db /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Pseudonaja_textilis/genome.fa -outfmt 7 -out $q.blast -evalue 0.001");

print "Blast done! :D Choose the proper hits: \n";

### The program shows on the screen the results from the Blast. Now we can choose which hits we want to analyze ###

system ("cat $q.blast");

print "\n";

### Now we will perform the following steps to complete the whole analysis ###

print "Enter the scaffold we are interested in: ";
$s = <STDIN>;
chomp ($s);

print "Enter the position of the first nucleotide of the hit/s: ";
$start = <STDIN>;
chomp ($start);

print "Enter the position of the last nucleotide of the hit/s: ";
$final = <STDIN>;
chomp ($final);

### The program makes a directory in which every file will be saved. The name will have the following structure: query.scaffold ###

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

### Now the program will extract the sequence (scaffold) that contains the region we are interested in from the index of the genome of P. textiles that was already given ###

system ("fastafetch /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Pseudonaja_textilis/genome.fa /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Pseudonaja_textilis/genome.index $s > $s.fa");

print "Fastafetch done! This is the length of your scaffold:\n";

system ("fastalength $s.fa");

print "Enter the scaffold length: ";
my $fastalength=<STDIN>;
chomp ($fastalength);

### The program extends the sequence 50 000 nucleotides on both ends ###

$start=$start-50000;

if ($start<0){
    $start=0;
}
print "We will consider this starting point: $start\n";

$length=($final+50000)-$start;

if ($fastalength< ($start+$length)){
    $length=$fastalength-$start-1;
}

print "...and this length: $length\n";

### The program will extract the genomic region of the scaffold that potentially contains the gene we are looking for ###

system ("fastasubseq $s.fa $start $length > $q.$s.genomic.fa"); 

print "Fastasubseq done! Now we will run the Exonerate...\n";

### Exonerate will give us the annotation of the gene that codes for the query protein. We won't use the "exhaustive" because it works too slow. With egrep, we obtain a file with the GFF annotation of the exonerate program we run ###

system ("exonerate -m p2g --showtargetgff -q $q -t $q.$s.genomic.fa | egrep -w exon > $q.$s.exonerate.gff");

print "Exonerate done! ;) But wait... let's do the fastaseqfromGFF to get the cDNA!\n";

### Now the program uses the GFF file to extract the cDNA sequence, without introns ###

system ("fastaseqfromGFF.pl $q.$s.genomic.fa $q.$s.exonerate.gff > $q.$s.nt.fa");

print "FastaseqfromGFF done! Let's translate it to proteins! \n";

### Now the program translates the cDNA sequence to proteins from all the possible reading frames ###

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

print "Fastatranslate done!\n";

### Now the program changes the * for X in the file ending with .p.fa and saves this new content in a new file ending with .pX.fa ###

open(INPUT, "$q.$s.p.fa");
undef $/;
$k = <INPUT>;
close INPUT;
$/ ="\n";

$k =~ s/\*/X/g;

open(my $fh, '>', "$q.$s.pX.fa");
print $fh $k;
close $fh;

### Now the program will do the global alignment of the resulting protein with the original sequence from the human genome by using TCOFFEE ###

print "Running T-Coffee…\n";

system ("t_coffee $q $q.$s.pX.fa > $q.$s.tcoffe");

### And now we will use the Genewise! ###

print "Is your hit in the reverse chain? (Check it in the .exonerate.gff)  ";
$rev = <STDIN>;
chomp ($rev);

print "Running Genewise…\n";

if ($rev eq "yes"){
system ("genewise -pep -pretty -cdna -gff -trev $q $q.$s.genomic.fa > $q.$s.genewiserev.fa"); #reverse
system ("mv $q.$s.genewiserev.fa $q.$s"); #moves to the directory
}

else {
system ("genewise -pep -pretty -cdna -gff $q $q.$s.genomic.fa > $q.$s.genewisefor.fa"); #forward
system ("mv $q.$s.genewisefor.fa $q.$s"); #moves to the directory

}

### Open the input file
open (SUBSEQ, "$q.$s.genomic.fa");

### Open the output file
open(NSUBSEQ, ">", "$q.$s.genomicNSeblastian.fa");
while(<SUBSEQ>) {
	if (/^>/){ # if line starts with '>'
		print NSUBSEQ $_;
	} elsif ($_ !~/^>/){ # if line does not starts with '>'
		s/[^atcgATGCN"\n"]/N/g;	# substitute all bases different than the canonical nucleotid bases (a, t, c, g, A, T, G, C, N) by an N.
		print NSUBSEQ $_;
	}
}
close (NSUBSEQ);
close (SUBSEQ);


### Finally all files are moved to the new directory ###

system ("mv $s.fa $q.$s");
system ("mv $q.$s.genomic.fa $q.$s");
system ("mv $q.$s.exonerate.gff $q.$s");
system ("mv $q.$s.nt.fa $q.$s");
system ("mv $q.$s.p.fa $q.$s");
system ("mv $q.blast $q.$s");
system ("mv $q.$s.pX.fa $q.$s"); 
system ("mv $q.$s.tcoffe $q.$s");
system ("mv $q.aln $q.$s");
system ("mv $q.html $q.$s"); 
system ("mv $q.$s.genomicNSeblastian.fa $q.$s");

print "THE END! GOOD JOB! All your files have been saved in the directory $q.$s\n";

### The program asks us if we want to analyze another scaffold ###

print "Do you want to analyze another scaffold? ";

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

if ($runagain eq "yes") {
    system ("perl perlPtextilis.pl");
}

else {
    print "Okay! Bye! :) \n";
}
