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

# Selenodeath.pl 
# Programa creado por los alumnos de la UPF de la asignatura de Bioinformática 2015-2016: Carla Cuní, Hugo Ramos, Miriam Rivera, Paola Santisteban, Anna Sochilina.
# Permite hacer predicciones de selenoproteínas del genoma de Macaca nemestrina por homología con otras especies. 
# Para culaquier duda, contactar con paola.santisteban01@estudiant.upf.edu

# Introducción de las variables ###############################################################################################################
 
my $index = "/cursos/BI/genomes/2015/Macaca_nemestrina/genome.index"; #Directorio del genoma indexado
my $genome = "/cursos/BI/genomes/2015/Macaca_nemestrina/genome.fa"; #Directorio del genoma
my $query_aa; #Variable donde introducimos la Query
print "Introduzca el nombre de la query sin extension .fa:";
$query_aa = <STDIN>; #Introducimos por pantalla el nombre de la Query.
chomp $query_aa;

# Blastall ####################################################################################################################################
# Se genera un BLAST simplificado tabulado (indicado en la mediante Blastall mediante -m 8), dentro de una carpeta con el nombre de la Query, que escogerá solo los Scaffolds con un e-value menor a 10-4, de la posición de inicio y de la posición final. Si el Blastall ya se había generado en otra ocasión, no volverá a procesarlo.

if (-e "$query_aa.fa") #Este "if" mira si existe el archivo con el nombre de la query introducido.
{
   print "Se ha encontrado la Query. Vamos a ver si el Blast ya se había generado anteriormente.\n";
   if (-e "./$query_aa/BLAST/${query_aa}_blast.fa") #Comprovamos con este "if" que el Blast existe.
   {
        print "El Blast ya se había generado anteriormente.\n";
   }
   else 
   {
	system ("mkdir -p ./$query_aa");
        system ("mkdir -p ./$query_aa/BLAST");
	print "No existia el archivo BLAST. Procedemos a ejecutarlo. Espere un momento.\n";
	system ("blastall -p tblastn -i $query_aa.fa -d $genome -o ./$query_aa/BLAST/${query_aa}_blast.fa -e 0.0001 -m 8");
	print "El archivo ${query_aa}_blast.fa se ha creado y almacenado correctamente en la carpeta $query_aa, subcarpeta BLAST\n";
   }
} 
else 
{   
   print "Error. No se ha encontrado la Query.\n"; #Si no habías creado la Query, el programa se detendrá para que la generes antes de volverlo a ejecutar.
   die; 
}

system ("more ./$query_aa/BLAST/${query_aa}_blast.fa"); #Enseñamos por el terminal el Blast.

# Fastafetch ####################################################################################################################################
# Seleccionamos el Scaffold en el que se encuentra el Hit. El archivo nombreshits.txt contiene los nombres de cada hit. Se leerá este archivo y se utilizará el nombre para crear el archivo fetch del Scaffold correspondiente. Si el archivo Fastafetch ya existe, no se creará uno nuevo.

system ("cut ./$query_aa/BLAST/${query_aa}_blast.fa -s -f 2 > ./$query_aa/nombrehits.txt"); #Cortamos los nombres del Scaffold y los pegamos en un archivo txt.
system ("more ./$query_aa/nombrehits.txt"); #Mostramos el txt creado.
my $hits = 0; #Variable para contar los hits.
open (NOMBLAST,"<./$query_aa/nombrehits.txt"); #Abrimos el txt y lo colocamos en el objeto NOMBLAST.
my @v = <NOMBLAST> ; #Creamos un vector en que cada Hit está en una posición del vector.
close (NOMBLAST);
for (@v)
{
  chomp;
  unless (-e "./$query_aa/fastafetch/${_}.fa") #Miramos si existe el Fastafetch. Si no, lo creamos.
     { 
        system ("mkdir -p ./$query_aa/fastafetch");
	system ("fastafetch $genome $index '$_' > './$query_aa/fastafetch/${_}.fa' ");
  	
	print "Archivo del Scaffold $_ Creado \n";
     }
  $hits = $hits + 1; #Contamos los Hits gracias al orden del vector.
}


print "$hits Hits encontrados en el Scaffold.\n";
print "A continuación realizaremos el Fastasubseq...\n";

# Fastasubseq #############################################################################################################################
# Leeremos el archivo generado por Blast y cogeremo los valores de inicio y final del Hit (columnas 8 y 9).

my %hits;
my $a;
my $start_blast;
my $final_blast;
my $start;
my $length;
my $leng;

system ("mkdir -p ./$query_aa/fastasubseq"); #Subcarpeta donde guardamos el Fastasubseq.
open (BLASTALL, "<./$query_aa/BLAST/${query_aa}_blast.fa");#Abrimos el Blast y lo metemos en el objeto BLASTALL.
while(<BLASTALL>) 
{
	chomp;
	my @blast = split("\t",$_);

        #Con este "if" determinamos cuales son las posiciones de inicio y final teniendo en cuenta la orientación de la cadena.
	if ($blast[8] < $blast[9]) {
		$start_blast = $blast[8];
		$final_blast = $blast[9];
    	}
    	else {
		$start_blast = $blast[9];
		$final_blast = $blast[8];
    	}	

        #Ampliamos el tamaño del scaffold +/- 50000 nucleótidos.
  	$start = $start_blast - 50000;
	if ($start < 0) {
		$start = 0;
	}
	$length = $final_blast - $start_blast + 100000;
	
	system ("fastasubseq './$query_aa/fastafetch/${blast[1]}.fa' $start $length > './$query_aa/fastasubseq/${blast[1]}_fastasubseq.fa' ");
}

close (BLASTALL);
print "Fastasubseq se ha creado correctamente y guardado en la carpeta correspondiente ;)\n";
print "Procediendo a realizar exonerate...\n";

#Generación de todas las carpetas ########################################################################################################

system ("mkdir -p ./$query_aa/exonerate");
system ("mkdir -p ./$query_aa/exonerate/gff");
system ("mkdir -p ./$query_aa/exonerate/cDNA_gff");
system ("mkdir -p ./$query_aa/exonerate/cDNApred");
system ("mkdir -p ./$query_aa/exonerate/proteins_gen");
system ("mkdir -p ./$query_aa/exonerate/transcripts");
system ("mkdir -p ./$query_aa/tcoffee");
system ("mkdir -p ./$query_aa/genewise");

# Exonerate ####################################################################################################################################
# Para buscar nuestra Query en el genoma de Macaca nemestrina.

print "Se está ejecuntando el Exonerate, espere un poco.\n";
opendir (SUBSEQDIR, "./$query_aa/fastasubseq/");
while (my $arch = readdir(SUBSEQDIR))
  	{
	next if ($arch =~ m/^\./);	
	chomp $arch;	
	my $arch_prot = $arch;
	$arch_prot =~ s/\_fastasubseq.fa//g;	
	chomp $arch_prot;
	system ("exonerate -m p2g --showtargetgff -q './$query_aa.fa' -t './$query_aa/fastasubseq/$arch' > './$query_aa/exonerate/gff/$arch_prot.gff' ");
	system ("egrep -w exon './$query_aa/exonerate/gff/${arch_prot}.gff' > './$query_aa/exonerate/cDNA_gff/${arch_prot}_cDNA.gff' ");

# FastaseqfromGFF #######################################################################################################################
# Para obtener la secuencia nucleotídica y Fastatranslate.

	system ("fastaseqfromGFF.pl './$query_aa/fastasubseq/$arch' './$query_aa/exonerate/cDNA_gff/${arch_prot}_cDNA.gff' > './$query_aa/exonerate/cDNApred/${arch_prot}_cDNApred.fa' ");
	system ("fastatranslate './$query_aa/exonerate/cDNApred/${arch_prot}_cDNApred.fa' > './$query_aa/exonerate/proteins_gen/${arch_prot}_proteins_gen.mfa' ");		
	}

print "Exonerate se ha creado y almacenado correctamente. En la carpeta Exonerate encontrarás el alineamiento, los exones, y la predicción.\n";
print "El siguiente paso es comparar la Query con la secuencia predicha. Un momento.\n";

# T-Coffee ##############################################################################################################################
# Alineará la predicción de Exonerate con la Query.

my $marclect = 1;
opendir (TRANSLATE, "./$query_aa/exonerate/proteins_gen/");
while (my $fail = readdir(TRANSLATE))
  {
	next if ($fail =~ m/^\./);	
	chomp $fail;
	my $protein_fail = $fail;
	$protein_fail =~ s/\_proteins_gen.mfa//g;
	chomp $protein_fail;
        system ("fastatranslate -F 1 './$query_aa/exonerate/cDNApred/${protein_fail}_cDNApred.fa' > './$query_aa/exonerate/transcripts/${query_aa}_${marclect}.fa' ");

        my $ux; #Cambio de * por X manual en el archivo Transcripts de la carperta Exonerate.
        print "Cuando hayas cambiado los * por X en el archivo Transcripts, pulsa Enter.";
        $ux = <STDIN>;
        chomp $ux;

        system("t_coffee $query_aa.fa './$query_aa/exonerate/transcripts/${query_aa}_${marclect}.fa' > './$query_aa/tcoffee/${query_aa}_${marclect}_tcoffee.fa' ");
	
        $marclect = $marclect + 1;	
        print "tcoffee creado\n";
   	print "\n";
}

print "T-coffee creado. Para hacer hacer el Genewise pulse Enter.\n"
;
<>;

# Genewise ############################################################################################################################
# Utiliza un algoritmo diferente a Exonerate para hacer una predicción de la secuencia, y el alineamiento con la Query.

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

print "Genewise generado.\n";
print "\n";

