#!/usr/bin/perl -w
#Este programa ha sido creado por Miren Aizpiri, Julia Alonso, Jaime Cano y Berta Martinez. Para contacto, comunicarse con jaime.cano01@estudiant.upf.edu
#Es condición necesaria tener una query con las U cambiadas por X


#Esta operación crea un archivo bash necesario para poder cambiar los * por X para realizar el tcoffee.

open(my $fh, '>', 'replace.sh');
print $fh "#!/bin/sh\nperl -pi -e 's/\\*/X/g' \$1;\n";
close $fh;
#print "Done!!\n";
system ("chmod u+x replace.sh");

#Introducción de las variables

my $my_index = "/cursos/BI/genomes/vertebrates/2014/Eidolon_helvum/genome.index";
my $my_genome = "/cursos/BI/genomes/vertebrates/2014/Eidolon_helvum/genome.fa";
my $query;
print "Introduzca la query sin terminación .fa:";
$query = <STDIN>;
chomp $query;
use strict;

#Creamos un archivo con el tamaño de cada cromosoma

print "Preparando los ficheros necesarios...(puede tardar, manténgase a la espera)\n";
unless (-e "genome_len"){
  system ("fastalength /cursos/BI/genomes/vertebrates/2014/Eidolon_helvum/genome.fa > genome_len");
}
print "Archivos generados!\n";
print "Dale a Enter para proceder a analizar el genoma\n";
<>;



#BLAST
#Creamos un BLAST simplificado (m8), utilizando solo aquellos scaffolds con un e-value menor de 10-4
if (-e "$query.fa") 
{
   print "Verificando la existencia del BLAST...!\n";
   if (-e "./$query/BLAST/${query}_blast.fa")
   {
	print "Perfecto!! El BLAST ya estaba creado!!\n"
   }
   else 
   {
	system ("mkdir -p ./$query");
	system ("mkdir -p ./$query/BLAST");
	print "Verificado! No hay archivo BLAST. Procedemos a crearlo para $query, manténgase a la espera...\n";
	system ("blastall -p tblastn -i $query.fa -d $my_genome -o ./$query/BLAST/${query}_blast.fa -e 0.0001 -m 8");
	print "Magnífico!! Ya tiene el archivo ${query}_blast.fa en la carpeta $query, subcarpeta BLAST\n";
   }
} 
else 
{   
   print "Error!! La query no existe...\n";
   print "Revisad la localización y nombre del archivo!\n";
   die;
   
}

#fastafetch
#Para realizar el fastafetch debemos leer los nombres de los scaffolds (hits.txt).

system ("cut ./$query/BLAST/${query}_blast.fa -s -f 2 > ./$query/hits.txt");
my $hits = 0;
open (NOMBLAST,"<./$query/hits.txt");
my @v = <NOMBLAST> ;
close (NOMBLAST);
for (@v)
{
  chomp;
  unless (-e "./$query/fastafetch/${_}.fa")
     { 
        system ("mkdir -p ./$query/fastafetch");
	system ("fastafetch $my_genome $my_index '$_' > './$query/fastafetch/${_}.fa' ");
  	#print "$_ creado!\n";
  	$hits = $hits + 1;
     }
}

print "$hits hits encontrados!\n";
print "Prepararos para hacer el fastasubseq, esperen...\n";

#Determinamos la longitud de los scaffolds relevantes y los guardamos en el hash (%leng)

my %leng;
open (lengths_file,"genome_len");
while(<lengths_file>) {
  my $len_line = $_;
  my @values = split(' ', $len_line);
  my $num_len = $values[0];
  my $key = $values[1];
  $leng{$key}=$num_len;
  #print "Key = ¡$key¡ Value = ¡$num_len¡ \n";
}
close(lengths_file);

#fastasubseq
#Leeremos el archivo generado por blast y leeremos el valor de inicio y final del hit (columnas 8 y 9)

my %hits;
my $a;
my $start_blast;
my $stop_blast;
my $start;
my $len;
system ("mkdir -p ./$query/fastasubseq");
open (BLASTALL, "<./$query/BLAST/${query}_blast.fa");
while(<BLASTALL>) 
{
   
  chomp;
  next if ($_ =~ m/^#/);
  my @blast = split("\t",$_);
  if ($blast[8] < $blast[9]) {
	$start_blast = $blast[8];
	$stop_blast = $blast[9];
    }
    else {
	$start_blast = $blast[9];
	$stop_blast = $blast[8];
    }

  if (exists($hits{$blast[1]})) 
  {
    if ($hits{$blast[1]}{"START"} > $start_blast)
	{
	$hits{$blast[1]}{"START"} = $start_blast;
	}
    if ($hits{$blast[1]}{"STOP"} < $stop_blast)
	{
	$hits{$blast[1]}{"STOP"} = $stop_blast;
	}
  }
  else 
  {
    $hits{$blast[1]}{"START"} = $start_blast;
    $hits{$blast[1]}{"STOP"} = $stop_blast;
  }
#Establecemos tamaño de 50000 nucleótidos para ampliar en nuestro scaffold
 if ($start_blast < 25000) {
    $start = 0;
  }
  else {
$start = $start_blast - 25000;
  }
  my $tot_len = $leng{$blast[1]};
  
  $len = $stop_blast - $start_blast + 50000;
  if ( ($len + $start) > $tot_len){
     $len = $tot_len-$start;
  }
 
#print "start = $start\nlen= $len\ntot = $tot_len\n";

system ("fastasubseq './$query/fastafetch/${blast[1]}.fa' $start $len > './$query/fastasubseq/${blast[1]}_fastasubseq.fa' ");

}

close (BLASTALL);

print "Fastasubseq creado y guardo en la carpeta correspondiente!\n";
print "Procediendo a realizar exonerate...\n";
print " \n";

#exonerate
system ("mkdir -p ./$query/exonerate");
system ("mkdir -p ./$query/exonerate/gff");
system ("mkdir -p ./$query/exonerate/index_gff");
system ("mkdir -p ./$query/exonerate/cDNA");
system ("mkdir -p ./$query/exonerate/r_frame");
system ("mkdir -p ./$query/exonerate/final_proteins");
system ("mkdir -p ./$query/genewise");
system ("mkdir -p ./$query/tcoffee");

opendir (SUBSEQDIR, "./$query/fastasubseq/");
while (my $file = readdir(SUBSEQDIR))
  {
	next if ($file =~ m/^\./);	
	chomp $file;	
	my $protein_file = $file;
	$protein_file =~ s/\_fastasubseq.fa//g;	
	chomp $protein_file;
	system ("exonerate -m p2g --showtargetgff -q './$query.fa' -t './$query/fastasubseq/$file' > './$query/exonerate/gff/$protein_file.gff' ");
	system ("egrep -w exon './$query/exonerate/gff/${protein_file}.gff' > './$query/exonerate/index_gff/${protein_file}_index.gff' ");
	system ("fastaseqfromGFF.pl './$query/fastasubseq/$file' './$query/exonerate/index_gff/${protein_file}_index.gff' > './$query/exonerate/cDNA/${protein_file}_cDNA.fa' ");
	system ("fastatranslate './$query/exonerate/cDNA/${protein_file}_cDNA.fa' > './$query/exonerate/r_frame/${protein_file}_r_frame.fa' ");	
	#system ("./replace.sh ./$query/exonerate/r_frame/${protein_file}_r_frame.fa"); 
  }
print "Exonerate creado y guardado en su carpeta!\n";
print "ENTER para continuar con el siguiente\n";
<>;

#tcoffee
#Fastatranslate produce un archivo con los 6 reading frames posibles. Hay que elegir el más adecuado.

my $j = 1;
opendir (TRANSLATE, "./$query/exonerate/r_frame/");
while (my $fail = readdir(TRANSLATE))
  {
	next if ($fail =~ m/^\./);	
	chomp $fail;
	system ("cat './$query/exonerate/r_frame/$fail' ");
	my $protein_fail = $fail;
	$protein_fail =~ s/\_r_frame.fa//g;
	chomp $protein_fail;
	print " \n";
	print "Esta es nuestra query de interés para comparar:\n";
        system ("more $query.fa\n");
        print " \n";	
	print "Para $query escribe cual de los alineamientos se parece más (-3,-2,-1,1,2,3):";	
	my $transcript = <STDIN>;	
	chomp $transcript;
	if ($transcript eq "-3")
	{
	system ("fastarevcomp './$query/exonerate/cDNA/${protein_fail}_cDNA.fa' > './$query/exonerate/final_proteins/${query}_revcomp.fa' ");
	system ("fastatranslate -F 3 './$query/exonerate/final_proteins/${query}_revcomp.fa'  > './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' ");
	system ("./replace.sh ./$query/exonerate/final_proteins/${protein_fail}_${j}.fa"); 
	system("t_coffee $query.fa './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' > ./$query/tcoffee/${protein_fail}_${j}.tcoffee.fa ");
	}	
	if ($transcript eq "-2")
	{
	system ("fastarevcomp './$query/exonerate/cDNA/${protein_fail}_cDNA.fa' > './$query/exonerate/final_proteins/${query}_revcomp.fa' ");
	system ("fastatranslate -F 2 './$query/exonerate/final_proteins/${query}_revcomp.fa'  > './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' ");
	system ("./replace.sh ./$query/exonerate/final_proteins/${protein_fail}_${j}.fa");
	system("t_coffee $query.fa './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' > ./$query/tcoffee/${protein_fail}_${j}.tcoffee.fa ");
	}
	if ($transcript eq "-1")
	{
	system ("fastarevcomp './$query/exonerate/cDNA/${protein_fail}_cDNA.fa' > './$query/exonerate/final_proteins/${query}_revcomp.fa' ");
	system ("fastatranslate -F 1 './$query/exonerate/final_proteins/${query}_revcomp.fa'  > './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' ");
	system ("./replace.sh ./$query/exonerate/final_proteins/${protein_fail}_${j}.fa");
	system("t_coffee $query.fa './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' > ./$query/tcoffee/${protein_fail}_${j}.tcoffee.fa ");
	}
	if ($transcript eq "1")
	{
	system ("fastatranslate -F 1 './$query/exonerate/cDNA/${protein_fail}_cDNA.fa' > './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' ");
	system ("./replace.sh ./$query/exonerate/final_proteins/${protein_fail}_${j}.fa");
	system("t_coffee $query.fa './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' > ./$query/tcoffee/${protein_fail}_${j}.tcoffee.fa ");
	}
	if ($transcript eq "2")
	{
	system ("fastatranslate -F 2 './$query/exonerate/cDNA/${protein_fail}_cDNA.fa' > './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' ");
	system ("./replace.sh ./$query/exonerate/final_proteins/${protein_fail}_${j}.fa");
	system("t_coffee $query.fa './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' > ./$query/tcoffee/${protein_fail}_${j}.tcoffee.fa ");
	}
	if ($transcript eq "3")
	{
	system ("fastatranslate -F 3 './$query/exonerate/cDNA/${protein_fail}_cDNA.fa' > './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' ");
	system ("./replace.sh ./$query/exonerate/final_proteins/${protein_fail}_${j}.fa");
	system("t_coffee $query.fa './$query/exonerate/final_proteins/${protein_fail}_${j}.fa' > ./$query/tcoffee/${protein_fail}_${j}.tcoffee.fa ");
	
	}
   $j = $j + 1;	
   print "\n";
   print "Tcoffee creado y guardado en su carpeta correspondiente!\n";
   print "\n";
  
   <>;
}
print "Finalizados todos los tcoffee! ENTER para realizar el genewise\n";
<>;

#genewise
opendir (SUBSEQDIR, "./$query/fastasubseq/");
while (my $arxiu = readdir(SUBSEQDIR))
  {
	next if ($arxiu =~ m/^\./);	
	chomp $arxiu;	
	my $protein_arxiu = $arxiu;
	$protein_arxiu =~ s/\_fastasubseq.fa//g;	
	chomp $protein_arxiu;
	system ("genewise -cdna -pep -pretty -gff -both './$query.fa' './$query/fastasubseq/$arxiu' > './$query/genewise/$protein_arxiu.genewise.fa' ");
  }	

print "Genewise creat!\n";
print "\n";
system ("mkdir -p ./$query/html");
system ("mkdir -p ./$query/aln");
system ("mkdir -p ./$query/$query.fa");
system ("mv $query.html $query/html/");
system ("mv $query.aln $query/aln/");
system ("mv $query.fa $query/$query.fa/");
print "Muchas gracias por su paciencia!!\n";
print "Adeu! Adios! Agur! T' hal'la frasec! Bye bye! ;)\n";
print "\n";
print "\n";

