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

#This program was developed by: Sara Congost, Julia Fabà, Gemma Nomdedéu-Sancho, Clara Sabiote Cartes, Cèlia Ventura i Gabarró. 

#Download data, unzip datafiles and place program inside Grup2 folder before running the program. Remember to give chmod permissions (chmod 777 ./code.pl).

my $p; #Danio_rerio protein

use Cwd qw (cwd);
my $dir = cwd; #code.pl direction

print "Enter Danio_rerio protein name:"; #The program asks for the name of protein
$p = <STDIN>;
chomp ($p);



#::::::::::::::CREATE DIRECTORY:::::::::::::::::::::::
#The program creates a directory where all the results will be saved

my $directory;
$directory = "$dir/results/$p";
if (! -e $directory) {
system("mkdir $dir/results/$p"); 
print("$p directory was not existant\n");
}



#::::::::::::::U SUBSTITUTION:::::::::::::::::::::::
#The program changes the U residue in Zebrafish genome for an X

system ("sed 's/U/X/g' $dir/data/$p/$p.aa.fa > $directory/$p.fa");
print "U changed to X\n";
print "\nTBLASTN in process\n"; 


#::::::::::::::::::::::::TBLASTN::::::::::::::::::::::::::
#It will compare the protein query to our subject

my $blast;
$blast = "$directory/$p.tblastn";
if (! -e $blast) { 
print("Blast was not existant\n");
system ("tblastn -query $directory/$p.fa -db $dir/data/genome.fa -evalue 0.01 -out $blast");
}

my $blast_tab;
$blast_tab = "$directory/$p.tblastn.tab";
if (! -e $blast_tab) {
print("Blast_tab was not existant\n");
system ("tblastn -query $directory/$p.fa -db $dir/data/genome.fa -outfmt 6 -evalue 0.01 -out $blast_tab");
}

system ("cat $blast_tab"); #The program gives us all the scaffolds with significant hits
system ("cut -f2 $blast_tab > $directory/$p.potential.scaffolds.fa"); 
system ("sort $directory/$p.potential.scaffolds.fa > $directory/$p.potential.sorted.scaffolds.fa");
system ("uniq $directory/$p.potential.sorted.scaffolds.fa > $directory/$p.potential.group.scaffolds.fa");

print "\nThe possible scaffolds are:\n";
system ("cat $directory/$p.potential.group.scaffolds.fa");


#::::::::::::SCAFFOLDS LOOP:::::::::::
#The program will create a loop to run all the following program for every scaffold

my $s; #scaffold

my $potentials;
$potentials = "$directory/$p.potential.group.scaffolds.fa";
open (INFO, '<', $potentials)or die $!;
chomp $potentials;

my $line;
while (<INFO>) {
$line=$_;
chomp $line;
$s=$line;
print "****$line****";

print "Running $s";


#::::::::::::::::FASTAFETCH::::::::::::::::
#It will extract the scaffold of interest from Anarrhichtys ocellatus' genome

system ("grep $s $blast_tab > $directory/$s.hits");
system ("cat $directory/$s.hits");

my $min;
$min = "$directory/$p.$s.min.fa";
my $max;
$max = "$directory/$p.$s.max.fa";

#hits positions extraction:
system ("cut -f9 $directory/$s.hits > $directory/$p.$s.positions.fa | cut -f10 $directory/$s.hits >> $directory/$p.$s.positions.fa");
system ("sort -n $directory/$p.$s.positions.fa > $directory/$p.$s.positions.sorted.fa");

#minimum and maximum positions extraction:
system ("head -1 $directory/$p.$s.positions.sorted.fa > $min");
system ("tail -1 $directory/$p.$s.positions.sorted.fa > $max"); 

print ("\n"); 
print ("Min=");
system ("cat $min");
print ("\n");
print ("Max=");
system ("cat $max");
print ("\n");
open FILE, "$min";
my $minimum;
$minimum =<FILE>; 
close FILE; 


#Start position definition:
my $start;
if ($minimum-50000>0){
$start=$minimum-50000}
else {$start=1};
#print "\nHit start is:$start\n";


#Running fastafetch:
system ("fastafetch $dir/data/genome.fa $dir/data/genome.index $s > $directory/$p.$s.fa");
print ("*******$s*********");
print "Fastafetch succesful\n";
#system ("cat $dir/results/$p/$p.$s.fa");


#Scaffold length definition:
system ("fastalength $directory/$p.$s.fa > $directory/$p.$s.length.fa");
system ("cat $directory/$p.$s.length.fa");
system ("cut -d ' ' -f1 $directory/$p.$s.length.fa > $directory/$p.$s.numlength.fa");

my $length;
$length = "$directory/$p.$s.numlength.fa";

my $numlength;
open (LEN, '<', $length) or die $!;
while (<LEN>){
$numlength=$_;
chomp $numlength;
print ("******$numlength*********");
}
close (LEN);
print ("************** $numlength ********");


#HIT LENGTH:
my $maximum;
open (FILE, '<', $max) or die $!;
while (<FILE>){
$maximum = $_;
chomp $maximum;
}
close FILE;
my $totlength;
if (($maximum+50000) < $numlength){
$totlength=($maximum-$start)+50000}
else {$totlength=$numlength-$start};


print ("\nThe start position is $start\n");
print("\n The scaffold length is $numlength\n");
print ("\n The hit length will be $totlength\n");



#::::::::::::::::::FASTASUBSEQ::::::::::::
#It will extract the region of interest from Anarrhichtys ocellatus' genome

my $subseq;
$subseq = "$directory/$p.$s.subseq.fa";
if (! -e $subseq) {
print("$p subseq was not existant\n"); 
system ("fastasubseq $directory/$p.$s.fa $start $totlength > $subseq");
}

#system ("cat $subseq"); 



#::::::::::::::::EXONERATE+egrep::::::::::::::::::
#It will store the exons of the protein of interest

my $exonerate;
$exonerate = "$directory/$p.$s.exonerate.gff";
if (! -e $exonerate) {
print("$p exonerate was not existant\n"); 
system ("exonerate -m p2g --showtargetgff -q $directory/$p.fa -t $subseq | egrep -w exon > $exonerate");
}

system ("cat $exonerate"); 




#:::::::::::::::FASTASEQFROMGFF:::::::::::::
#The program will predict a cDNA sequence 

my $prednuc;
$prednuc = "$directory/$p.$s.pred.nuc.fa";
if (! -e $prednuc) {
print("Predicted nucleotide sequence was not existant\n");
system ("fastaseqfromGFF.pl $directory/$p.$s.subseq.fa $directory/$p.$s.exonerate.gff > $prednuc");
}

print "\nThe predicted nucleotide sequence is:\n"; 
system ("cat $prednuc");



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

system ("fastatranslate -f $directory/$p.$s.pred.nuc.fa -F 1 > $directory/$p.$s.pred.prot.fa"); 

print "\nThe predicted aminoacid sequence is:\n"; 
system ("cat $directory/$p.$s.pred.prot.fa"); 

#Change * to X:
system ("sed 's/*/X/g' $directory/$p.$s.pred.prot.fa > $directory/$p.$s.pred.x.prot.fa");
system ("cat $directory/$p.$s.pred.x.prot.fa");



#:::::::::::::::::::TCOFFEE::::::::::::::::::::
#It will align the predicted protein from Fastatranslate (Anarrhichtys ocellatus) and the query sequence (Danio Rerio)

#my $tcoffee;
#$tcoffee = "$dir/results/$p/$p.$s.tcoffee";
#if (! -e $tcoffee) {
#print("$p $s tcoffee was not existant\n");
system ("t_coffee $directory/$p.$s.pred.x.prot.fa $directory/$p.fa > $directory/$p.$s.tcoffee"); 
#}
system ("cat $directory/$p.$s.tcoffee");
}


#:::::::::::::CLOSE SCAFFOLDS LOOP::::::::::::::::

close (INFO);



  




 
