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

# Selenocuqui5.pl 
# Programa editado por los alumnos de la UPF de la asignatura de Bioinformática 2016-2017: Pedro Fortes, Sergi Garcia, Nuria Güil, Alba Ponce.
# Permite hacer predicciones de selenoproteínas del genoma de Cebus capucinus imitator por homología con otras especies. 
# Para cualquier duda, contactar con nuria.guil01@estudiant.upf.edu

# Introducción de las variables ###############################################################################################################
 
my $index = "/cursos/20428/BI/genomes/2016/Cebus_capucinus_imitator/genome.index"; #Directorio del genoma indexado
my $genome = "/cursos/20428/BI/genomes/2016/Cebus_capucinus_imitator/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 = ; #Introducimos por pantalla el nombre de la Query.
chomp $query_aa;

# Tblastn ####################################################################################################################################
# Se genera un BLAST simplificado tabulado dentro de una carpeta con el nombre de la Query, que escogerá solo los Scaffolds con un e-value menor a 10-2, de la posición de inicio y de la posición final. 


system ("mkdir ./$query_aa");
system ("mkdir ./$query_aa/BLAST");
system ("tblastn -query $query_aa.fa -db $genome -out ./$query_aa/BLAST/${query_aa}_blast.fa -evalue 0.01 -outfmt 6");

print "Cuqui, el archivo ${query_aa}_blast.fa se ha creado y almacenado correctamente en la carpeta $query_aa, subcarpeta BLAST\n";


# Fastafetch ####################################################################################################################################
# Seleccionamos el Scaffold en el que se encuentra el Hit. El archivo nombrehits.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.

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 =  ; #Creamos un vector en que cada Hit está en una posición del vector.
close (NOMBLAST);

system ("mkdir ./$query_aa/fastafetch");

for (@v)
{
  chomp;
	print '$_';
	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 cogeremos 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 %dic; 

#my open (scaffold_length, "/cursos/20428/BI/genomes/2016/Cebus_capucinus_imitator/genome.length"); 

system ("mkdir ./$query_aa/fastasubseq"); #Subcarpeta donde guardamos el Fastasubseq.
open (TBLASTN, "<./$query_aa/BLAST/${query_aa}_blast.fa");#Abrimos el Blast y lo metemos en el objeto TBLASTN.
while() 
{
	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;
	#Ampliamos el tamaño de la longitud + 100000 nucleótidos. 

#En esta condición del "if" estamos indicando que si la longitud obtenida es menor a la longitud del Scaffold, que ésta sea igual a la del Scaffold. 

if ($length < ($final_blast - $start_blast)) {

	$length = $final_blast - $start_blast; 
	$dic{$_} = $length;
	system ("fastasubseq './$query_aa/fastafetch/${blast[1]}.fa' $start $length > './$query_aa/fastasubseq/${blast[1]}_fastasubseq.fa' ");
}

else { 
$dic{$_} = $length;
	system ("fastasubseq './$query_aa/fastafetch/${blast[1]}.fa' $start $length > './$query_aa/fastasubseq/${blast[1]}_fastasubseq.fa' ");

}


close (TBLASTN);
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 Cebus Capucinus.

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 = ;
        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"; 
}