#!/usr/bin/perl -w

use strict;


############FEM UN PROGRAMA QUE OBRI ELS BLASTS I FILTRI PER VEURE SI HI HA DUPLICACIONS#############


opendir(DIR, "./selenoproteines") or die "Could not open la carpeta\n";
while (my $protein = readdir(DIR)){
    next if $protein =~ /^\./;

##el programa va llegint les carpetes i va obrint els documents .blast, et mostra per pantalla el document output del blast i et pregunta si el vols analitzar. Si esculls que no, passa al ##següent pero si detectes que hi pot haver més d’un contig que val la pena analitzar i esculls que si, el programa et preguntarà quin contig vols analitzar i l’assignarà a la variable 
##$contig, realizant tots els passos necessaris per a la predicció de gens

open (BLAST, "./selenoproteines/$protein/$protein.canviat.blast") or die "no es pot obrir el blast";

  print "aquest es el blast: $protein\n"; 
    system ("more ./selenoproteines/$protein/$protein.canviat.blast");

    print "vols analitzar el blast $protein ? (y or n)\n";

    my $question = <STDIN>;
    chomp $question;
    if ($question eq "y"){
	print "ara anem a analitzar el blast\n";

	print "quin contig es? \n";
	my $contig = <STDIN>;
	chomp $contig;

	system ("mkdir ./selenoproteines/$protein/$contig");

	system ("egrep $contig ./selenoproteines/$protein/$protein.canviat.blast > ./selenoproteines/$protein/$contig/$contig.retallat.txt");

#################FASTAFETCH#############################

 system ("fastafetch /cursos/20428/BI/genomes/2017/Nannopterum_harrisi/genome.fa  /cursos/20428/BI/genomes/2017/Nannopterum_harrisi/genome.index '$contig' > ./selenoproteines/$protein/$contig/$contig.fa ");

	print "fetch fet\n";


#########FASTA SUBSEQ##########

system ("cut ./selenoproteines/$protein/$contig/$contig.retallat.txt -f 9,10 | sort -n >./selenoproteines/$protein/$contig/$contig.posicions.txt"); #tallem les dos columnes on hi ha les posicions i les posem en un nou arxiu
system ("xargs -n1 < ./selenoproteines/$protein/$contig/$contig.posicions.txt | sort -n > ./selenoproteines/$protein/$contig/$contig.minmax.txt"); #ajuntem les dos columnes en una sola columna, ho ordenem numericament

#	open (MINMAX, "./selenoproteines/$protein/$protein.minmax.txt") or die "no es pot obrir el document";
	my $minim = `head -1 ./selenoproteines/$protein/$contig/$contig.minmax.txt` - 50000;# `aixo fa que actuis en el system ia  la vegada puguis capturar el resultat, aleshores agafem el inici de la sequencia i li restem les 50000 bases de marge#
	my $maxim = `tail -1 ./selenoproteines/$protein/$contig/$contig.minmax.txt` + 50000;

	if ($minim < 0){
	    $minim = 0;
	}
	my $length = $maxim - $minim;
	    print "el minim es: $minim\n";
	    print "el maxim es: $maxim\n";
	print "la llargada es: $length\n";

	system ("fastasubseq './selenoproteines/$protein/$contig/$contig.fa' $minim $length > './selenoproteines/$protein/$contig/$contig.fastasubseq.fa' ");
	print "fastasubseq terminao\n";# de l'arxiu fasta del contig del nostre pajaro extreiem la zona que ens interesa comparar #



###################EXONERATE##########

	print "comencem amb el exonerate\n";

	system ("exonerate -m p2g --exhaustive yes --showtargetgff -q './selenoproteines/$protein/$protein.canviat.fa' -t './selenoproteines/$protein/$contig/$contig.fastasubseq.fa' > './selenoproteines/$protein/$contig/$contig.exonerate.gff' ");#predim la proteina comparant la proteina que nosaltres hem tret d'internet amb el dna del nostre pajaro#

	print "exonerate done\n";
	system ("egrep -w exon './selenoproteines/$protein/$contig/$contig.exonerate.gff' > './selenoproteines/$protein/$contig/$contig.exon.gff' ");#obtenim nomes els exons que codifiquen per la nostra proteina predita#


######################FASTASEQFROM GFF AND FASTATRANSLATE##################

	system ("fastaseqfromGFF.pl './selenoproteines/$protein/$contig/$contig.fastasubseq.fa' './selenoproteines/$protein/&contig/$contig.exon.gff' > './selenoproteines/$protein/$contig/$contig.prednt.fa' ");#traduim la proteina predita a partir del exons a nt#

	print "hem realitzat el fastaseq\n";
	system ("fastatranslate -f './selenoproteines/$protein/$contig/$contig.prednt.fa' -F 1 > './selenoproteines/$protein/$contig/$contig.predaa.fa' ");
	print "ara tenim la sequencia en aa\n"; #pasem la proteina predita en nt a aa#

	system ("sed 's/*/X/g' './selenoproteines/$protein/$contig/$contig.predaa.fa' > './selenoproteines/$protein/$contig/$contig.Xpredaa.fa' ");
	print "hem canviat els * a U\n";#canviem els llocs on hi havia * per U#


###########################TCOFFEE#################################

	system ("t_coffee './selenoproteines/$protein/$protein.canviat.fa'  './selenoproteines/$protein/$contig/$contig.Xpredaa.fa' > './selenoproteines/$protein/$contig/$contig.tcoffee.fa' "); #comparem al proteina baixada del selenoDB sense U amb la proteina predita per nosaltres#
##!!!!!!!!!!!!! larxiu ha de ser .fa o .aln???###
##ens surt a pantalla el tcoffee i no volem que pasi##

	print "tcoffee acabao\n";

##################GENEWISE##################

	system ("egrep -w gene './selenoproteines/$protein/$contig/$contig.exonerate.gff' > './selenoproteines/$protein/$contig/$contig.gene.gff' ");

	my $direction = `cut ./selenoproteines/$protein/$contig/$protein.gene.gff -f 7`;
	    chomp $direction;
	print " la direccio és: $direction\n";


	if ($direction eq "+"){
	   print "el gen es troba a la cadena fwd\n";
 system("genewise -pep -pretty -cdna -gff ./selenoproteines/$protein/$protein.canviat.fa ./selenoproteines/$protein/$contig/$protein.fastasubseq.fa > ./selenoproteines/$protein/$protein.genewise.fa ");
	}
	elsif ($direction eq "-"){
 print "el gen es troba a la caden rev\n";
 system ("genewise -pep -pretty -trev -cdna -gff ./selenoproteines/$protein/$protein.canviat.fa ./selenoproteines/$protein/$contig/$protein.fastasubseq.fa > ./selenoproteines/$protein/$contig/$protein.genewise.fa ");

	}
	else {
	    print "Alguna cosa va malament\n";
	}



	

 }




}