#!/usr/bin/perl -w
use strict;
#Programa optimitzat per Ariadna Cobos, Marta Dacosta, Júlia Garcia i Gemma Nuñez, a partir dels programes de cursos anteriors.

print "PROGRAMA PROVA. Per continuar, prem enter\n";
<>;

#creem una carpeta selzebrafish on hi tenim tantes carpetes com proteïnes hi hagi a la selenodb

######################### INTRODUCCIÓ DE LA SEQÜÈNCIA ##########################
my $seq;
print "Introduïm el nom de la seqüència (sense el .fa):";
$seq = <STDIN>;
chomp $seq;

######################## Canviar U per X ###################################

my $x;
my @v;
my $l;
my $input;
my $output;

open ($input, "<$seq/$seq.fa");       #fitxer d'entrada serà la query que estem analitzant
open ($output, ">$seq/${seq}_x.fa");   #fitxer sortida serà la query que haurem modificat amb el miniprograma creat

while (<$input>){                         
	@v = split ("",$_);             
	$x = 0;
	$l = scalar(@v);                #la longitud serà la del vector creat, és a dir, la de la nostra query
	while ($l > 0){
		if ($v[$x] eq 'U'){     #si troba una U la canviarà per una X. Si troba un #, @ o % ho eliminarà.
			$v[$x] = 'X';			
		}
		if ($v[$x] eq '#') {
		    $v[$x] = '';
		}
		if ($v[$x] eq '@') {
		    $v[$x] = '';
		}
		if ($v[$x] eq '%') {
		    $v[$x] = '';
		}
		$x = $x + 1;
		$l = $l - 1;
	}
	print $output @v;

}

close ($input);
close ($output);

print "Ja tens la query modificada amb X enlloc de U\n";
print "Prem enter per començar a realitzar el blast\n";
<>;

########################## REALITZACIÓ DEL BLAST #################################################

if (-e "$seq/BLAST/$seq.blast") {
	print "Ja hi ha un blast per aquesta proteïna\n";
} 
else {

	system ("mkdir $seq/BLAST");      #creem un directori nou anomenat BLAST, on hi posarem tot el que creem.

	print "S'està creant el blast\n";   

	system ("tblastn -query $seq/${seq}_x.fa -db /cursos/20428/BI/genomes/2016/Ictalurus_punctatus/genome.fa -outfmt 7 -out $seq/BLAST/$seq.blast -evalue 0.0001");      #s'esta realitzant el blast a partir de la query modificada. Tambe ho tabulem i ens quedem nomes amb les queries amb un evalue maxim de 0.001

open ('BLAST', "<$seq/BLAST/$seq.blast");
print "Aquesta és la query que ens surt amb el Blast: $seq.blast\n";
system ("cat $seq.blast $seq/BLAST/$seq.blast");
close ('BLAST');

	print "Prem enter per continuar amb el FASTAFETCH\n";
	<>;
	
}

######################## REALITZACIÓ FASTAFETCH ####################################################

if (-e "$seq/FETCH"){
	print "Ja hi ha un fetch per aquesta query\n";
} 
else {
	system ("mkdir $seq/FETCH"); #creem un directori anomenat FETCH on hi guardarem el resultat del fastafetch.
    	system ("cut $seq/BLAST/$seq.blast -s -f 2 | sort | uniq > $seq/FETCH/regions.txt");   #tallem la segona columna del fitxer resultant del blast, que conté els noms dels cromosomes dels hits per a aquella proteina, i els posem en un fitxer de text.
    	open ('regions', "<$seq/FETCH/regions.txt");    #obrim el fitxer de text regions que conte el nom dels cromosomes
	my @v = <regions>;     #passem els noms dels cromosomes del fitxer de text a un array (vector).
	close ('regions');     #tanquem el fitxer de text amb els noms.
	my $f = 0;
	for (@v) {             #per cada posició del vector, és a dir, per cada cromosoma, farem:
		chomp;
       		system ("fastafetch /cursos/20428/BI/genomes/2016/Ictalurus_punctatus/genome.fa /cursos/20428/BI/genomes/2016/Ictalurus_punctatus/genome.index '$_' > '$seq/FETCH/$_.txt'");     #fem el fastafetch. El fitxer resultant es guardarà dins la carpeta FETCH amb el nom del cromosoma.
	$f = $f + 1;		
	}
	
	print "Ja s'ha realitzat el Fetch. Nombre de hits: $f\n";

}
print "Prem enter per fer els fastasubseq\n";
<>;


######################### REALITZACIÓ FASTASUBSEQ ################################################
my $pos_start;
my $pos_stop;
my $start;
my $length;
my $numcaracter;
my $numregiocar;
my $longdef;
my $i;
my %hits;
$i = 0;

if ( -e "$seq/SUBSEQ") {
	print "Ja tens els teus hits YO!\n";
}

else {

	system ("mkdir $seq/SUBSEQ");
	open ('BLAST', "<$seq/BLAST/$seq.blast");
	while (<BLAST>){

		#Tall de cada línia del blast

		chomp; 
		next if ($_ =~ m/^#/);
		my @blast = split ("\t", $_); 

		# Establiment de les posicions d'inici i final. Es té en compte si el sentit és positiu o negatiu.
  		if ($blast[8] < $blast[9]) {         #si el valor de la columna 8 del blast es mes petita que el valor de la columna 9:
			$pos_start = $blast[8];      #fes que la posicio inicial sigui el de la columna 8
			$pos_stop = $blast[9];       #fes que la posicio de la columna 9 sigui el final
    		}
    		else {                               #sino es aixi
			$pos_start = $blast[9];      #fem que el valor de la columna 9 sigui l'inici
			$pos_stop = $blast[8];       #i el valor de la columna 8 sigui el final
    		}

		if (exists($hits{$blast[1]})){       #si ja existex un array (anomenat hits) amb el nom del scaffold (columna 1). 
   
	                                             #en el moment que ja existeixi l'array: 
   		if ($hits{$blast[1]}{"START"} > $pos_start){     #si el start del nom del scaffold de la columna 1 es mes gran que la posicio inicial que ja teniem guardada
		    $hits{$blast[1]}{"START"} = $pos_start;      #deixa la posicio inicial tal i com estava
			}
    		if ($hits{$blast[1]}{"STOP"} < $pos_stop){       #si el final del nom del scaffold de la columna 1 es mes gran que la posicio final que ja teniem guardada
		    $hits{$blast[1]}{"STOP"} = $pos_stop;        #deixa la posicio final tal i com estava 
			}
  		}
  		else{                                            #en el cas que no sigui aixi
   		
	                                                         #Si no existeix un array amb el nom del scaffold (columna 1)
   		$hits{$blast[1]}{"START"} = $pos_start;  
    		$hits{$blast[1]}{"STOP"} = $pos_stop;
  		}

	                                                 	# Guardem el nombre de caràcters total de cada fitxer fetch a l'arxiu caracterseq.txt
		system ("wc '$seq/FETCH/${blast[1]}.txt' -m > '$seq/SUBSEQ/caracterseq.txt' ");      #compte el nombre de caracters total del fitxer amb la comanda wc -m. Ho guardem en un fitxer anomenat caracterseq
		system ("cut -d ' ' -f1 '$seq/SUBSEQ/caracterseq.txt' > '$seq/SUBSEQ/numcaracter.txt' ");   #Ens quedem unicament amb la columna 1, ja que es on hi ha el nombre de caracters total del fitxer

	                                                   
		system ("grep '>' './$seq/FETCH/${blast[1]}.txt' > '$seq/SUBSEQ/nomregio.txt' ");  #ens quedem amb la línia que comença amb el simbol '>', que correspon al nom de capçalera del fitxer (nom scaffold...). Ho guardem en un fitxer nou: nomregio
		system ("wc '$seq/SUBSEQ/nomregio.txt' -m > '$seq/SUBSEQ/caracternom.txt' ");      #compte el nombre de caracters del fitxer nom regio. Es a dir, compte els caracters de la capçalera del fitxer.
		system ("cut -d ' ' -f1 '$seq/SUBSEQ/caracternom.txt' > '$seq/SUBSEQ/numregiocar.txt' ");   #talla el primer valor de la columna, que correspon al nombre de caracters total de la capçalera del fitxer.

	                                                            

		open ('caracter', "<$seq/SUBSEQ/numcaracter.txt");   #obra el fitxer que conte el nombre de caracters total del scaffold
		my $numcaracter = <caracter>;                        #Posem la longitud del scaffold amb el nom a la variable $numcaracter
		close ('caracter');
		$numcaracter = $numcaracter - 1;                     #li restem 1 ja que comença per la posicio 0.

	                                                             	
		open ('regcaracter', "<$seq/SUBSEQ/numregiocar.txt");
		my $numregiocar = <regcaracter>;                    #Posem la longitud del nom de la regió a la variable $numregiocar

		close ('regcaracter');
		$numregiocar = $numregiocar - 1;

	                                                    	

		$longdef = $numcaracter - $numregiocar;         #Posem la longitud definitiva del scaffold a la variable $longdef. Hem de restar a la longitud total del fitxer el numero de caracters de la capçalera del fitxer

		# Esborrem tots els arxius intermitjos generats

		system ("rm '$seq/SUBSEQ/caracterseq.txt' "); 
		system ("rm '$seq/SUBSEQ/numcaracter.txt' ");
		system ("rm '$seq/SUBSEQ/nomregio.txt' "); 
		system ("rm '$seq/SUBSEQ/caracternom.txt' "); 
		system ("rm '$seq/SUBSEQ/numregiocar.txt' ");

		#Establim el marge de 25000 i Es compara $start + $length amb la longitud del scaffold sense el nom, i si és més gran farem que $length adquireixi el valor de (longitud de l'scaffold menys $stop).

		$start = $pos_start - 25000;                    #La posicio incial sera la calculada i anem 25000  posicions abans
		if ($start < 0) {                               #No ens interessen les posicons negatives, per tant, en el cas que anant 25000 posicions abans arribem a numeros negatius, farem que la posicio inicial sigui 0.
			$start = 0;
		}
		$length = $pos_stop - $pos_start + 25000;
		if ($length > $longdef) {
			$length = $longdef - $pos_start; 

		}
		
		# Crida de la comanda Fastasubseq
		system ("fastasubseq '$seq/FETCH/${blast[1]}.txt' $start $length > '$seq/SUBSEQ/${blast[1]}_${i}_subseq.txt' ");

		$i = $i + 1;
	}

close (BLAST);
}

print "els hits estan a la  carpeta SUBSEQ\n";

print "Prem enter per veure si els teus hits codifiquen per una proteïna\n";
<>;

##############################REALITZACIÓ DE L'EXONERATE########################

my $file;
my $subseq_file;

if (-e "$seq/EXONERATE"){
		print "Ja s'ha fet exonerate per aquesta query\n";
}
else {
                                             #creem totes les carpetes que necessitarem 
	system ("mkdir $seq/EXONERATE");
	system ("mkdir $seq/EXONERATE/GFF");
	system ("mkdir $seq/EXONERATE/cDNAGFF");
	system ("mkdir $seq/EXONERATE/fastaseqfromGFF");
	system ("mkdir $seq/EXONERATE/fastatranslate");

	opendir (SUBSEQ, "$seq/SUBSEQ/");     #obra el directori SUBSEQ (on es troben totes les sequencies de les proteines acotades)
	while ($file = readdir (SUBSEQ)){     #assigna un fitxer qualssevol del directori SUBSEQ a la variable $file. Aixo ho faras per tots els fitxers del directori.Mentre estiguis treballant amb aquest fitxer....
		next if ($file =~ m/^\./);    #descartem aquells fitxers que no ens interessen (estranys, ocults...). 
		chomp $file;                
		$subseq_file = $file;         #assignem a la variable subseq_file el nom del fitxer (file). La variable subseq nomes servira per donar nom al fitxer de sortida.
		$subseq_file =~ s/\_subseq.txt//g;  #eliminem l'extensio (_subseq.txt) del fitxer d'origen. D'aquesta manera, a la variable subseq_file tindrem el nom del scaffold com a tal.
		chomp $subseq_file;                 #guardem en memoria el nom del scaffold dins de la variable subseq_file
		system ("exonerate -m p2g --showtargetgff -q '$seq/${seq}_x.fa' -t '$seq/SUBSEQ/$file' --exhaustive > '$seq/EXONERATE/GFF/$subseq_file.gff'");                         #Primer pas del exonerate. Partim de la query modificada amb la X
		system ("egrep -w exon '$seq/EXONERATE/GFF/$subseq_file.gff' > '$seq/EXONERATE/cDNAGFF/${subseq_file}_cDNA.gff'");  #ens quedem amb els exons
		system (" fastaseqfromGFF.pl '$seq/SUBSEQ/$file' '$seq/EXONERATE/cDNAGFF/${subseq_file}_cDNA.gff' > '$seq/EXONERATE/fastaseqfromGFF/${subseq_file}_cDNA.fa' ");                #ho tenim tot a nucleotids
		system (" fastatranslate -F 1 '$seq/EXONERATE/fastaseqfromGFF/${subseq_file}_cDNA.fa' > '$seq/EXONERATE/fastatranslate/${subseq_file}_aa.fa'");                                     #fem el fastatranslate per passar-ho a proteines

	}
}

print "Ja tenim fets els cDNA's i la seva traducció\n";

print "Prem enter per continuar canviant les * per les X\n";
<>;


############################### CANVI * per X #######################

my $fastatranslate="./$seq/EXONERATE/fastatranslate"; #assignem a la variable fastatranslate els fitxers que es troben dins el directori fastatranslate
				
opendir my $ft, $fastatranslate or die "Could not open '$fastatranslate' for reading: $!\n";		#obrim el directori.  Si no es pot obrir, el programa mor.
my @filtered_doc;									 

    while (my $fitxer_fastatranslate = readdir $ft) {          #Assignes un fitxer del directori a la variable fitxer_fastatranslate per treballar sobre aquest. Aixo ho faras per tots els fitxers del directori
	next if ($fitxer_fastatranslate =~ m/^\./); #evites els fitxers extranys, ocults...
	system("mkdir -p ./$seq/EXONERATE/fastatranslate_X"); #creem un directori on aniran els fitxers que ja tinguin el * canviat per una X.
open (TRANS, "./$seq/EXONERATE/fastatranslate/$fitxer_fastatranslate") or die "can't open '$fitxer_fastatranslate'\n";
my $nom_fitxersortida = "./$seq/EXONERATE/fastatranslate_X/$fitxer_fastatranslate.X";		# Nom amb el que anomenare el fitxer de sortida
open(OUTPUT, ">", $nom_fitxersortida);						# Obro l'output per escriure-hi
while(<TRANS>) {								# Llegeixo linia per linia
										
										# l'identificador no el canvio pero tambe l'imprimeixo al fitxer de sortida
	
	if (/^>/) {                                                             #si el fitxer comença amb el simbol ">"
		print OUTPUT $_;                                                
	}

	elsif ($_ !~/^>/){							# Si no comença amb el simbol ">"
		s/\*/X/g;							# La variable "$_" es cada linia per defecte. El que faig es mirar si hi ha * i substuir per X	
		print OUTPUT $_;	
	}
}
close(OUTPUT);
}

print "Ja tenim canviades els * per les X\n";

print "Prem enter per continuar amb T-coffee\n";
<>;


######################## REALITZACIO DEL T-COFFEE ###########################
system ("mkdir $seq/TCOFFEE");
my $fastatranslate_mod="./$seq/EXONERATE/fastatranslate_X"; #assignem a la variable fastatranslate_mod els fitxers que es troben dins el directori fastatranslate_X (aquells modificats amb la X enlloc del *)
				
opendir my $fm, $fastatranslate_mod or die "Could not open '$fastatranslate_mod' for reading: $!\n";		#obrim el directori. Si no es pot obrir, el programa mor.
my @filtered_arxius;									 

    while (my $fitxer_fastatranslate_mod = readdir $fm) { #Assignes un fitxer del directori a la variable file per treballa sobre aquest. Aixo ho faras per tots els fitxers del directori
	next if ($fitxer_fastatranslate_mod =~ m/^\./);    #eliminem els fitxers extranys, ocults...

system ("t_coffee $seq/${seq}_x.fa './$seq/EXONERATE/fastatranslate_X/$fitxer_fastatranslate_mod' > './$seq/TCOFFEE/$fitxer_fastatranslate_mod.tcoffee' ");
}
print "t-coffee fet\n"; #Un cop acabat, printem el missatge informatiu.
