#!/usr/bin/perl
use strict;
use warnings;

my $directory = "./Proteins_musmusculus";
my $genome = "/cursos/20428/BI/genomes/2016/Mus_spretus/genome.fa";


opendir (DIR, $directory) or die "Not possible to open $directory directory\n"; #obre la carpeta on tenim les subcarpetes amb les proteines en .fa
while (my $dir = readdir(DIR)) {
	next if $dir eq '.' or $dir eq '..'; #fa que no entri al directori actual (.) o a l'anterior (..). Next fa que s'ometi el codi de baix i salta a la seguent iteracio

	my $fasta_file = "$directory/$dir/$dir.fa"; #posa dins la variable $fasta_file el fitxer fasta amb la proteina
	
	for ($fasta_file) {
	    system ("sed 's/U/X/g' $fasta_file > $directory/$dir/prot_nclean.fa"); #s=substituir U per X en g=global (tota la sequencia)
	    system ("sed 's/[#@%]//g' $directory/$dir/prot_nclean.fa > $directory/$dir/prot_clean.fa"); #treure simbols raros
	}

	my $fasta_file_X = "$directory/$dir/prot_clean.fa"; #$fasta_file_X sera la proteina on s'ha canviat U per X
	print ("Performing tblastn for $fasta_file_X \n");

	system ("tblastn -query $fasta_file_X -db $genome -out $directory/$dir/resultat_blast.fa -outfmt 6");
	system ("cut -f 2,9,10,11 $directory/$dir/resultat_blast.fa | sort -n > $directory/$dir/filtrat_blast.fa");



##########PROCESSING BLAST##################

open(BLAST, '<', "$directory/$dir/filtrat_blast.fa") || die("can't open filtrat_blast.fa");

#el '<' es perq estas obrint un fitxer per llegirlo

#el nom del fitxer que vol obrir es posa entre cometes (en aquest cas es una taula amb nom de contig ordenat, coordenades i E-value)

#die et mostra si hi ha algun problema i no pot obrir el fitxer multifasta


my $start;
my $length;
my @coordbo = ("empty");
my @starts = ("0");
my @lengths = ("0");
my $n = 0;			#nombre de contigs (començant per zero)
my @contigs = ("empty");			#nom dels diferents contigs
my @direction = ("empty");	#direccio cadena, F es sense, R antisense
my @all_contigs;		#vector amb tots els contigs
my $num_hit = 0;		#comptador dels hits bons (e value < 0.01)
my @all_starts;
my @all_ends;
my @contigs_dir;
my $contig;

while(<BLAST>){

	my @hit = split(/\t/,$_);		#creem un vector que te tot el contingut d'una fila de la taula del blast, a cada posicio del vector hi ha el que esta separat per tabulador: \t)
	
	if($hit[3] < 0.01){			#escollim nomes hits amb E-value inferior a 0.01
		if ($hit[1] < $hit[2]){		#si esta en direccio sense
			$direction[$num_hit] = "F";
			$all_starts[$num_hit] = $hit[1];
			$all_ends[$num_hit] = $hit[2];
		}
		elsif ($hit[2] < $hit[1]){	#si esta en direccio reverse
			$direction[$num_hit] = "R";
			$all_starts[$num_hit] = $hit[2];
			$all_ends[$num_hit] = $hit[1];
		}			
		$all_contigs[$num_hit] = $hit[0];
		$num_hit++;
	}
}

close(BLAST);		
#print ("@all_contigs, starts: @all_starts, ends: @all_ends\n");
#print "$all_contigs[$n]\n";

for ($n = 0; $n < $num_hit; $n++){
	$contigs_dir[$n] = $all_contigs[$n].$direction[$n];
}

$n = 0;
for (my $comptador = 0; $comptador < $num_hit; $comptador++){
	if ($contigs[$n] eq "empty" || $contigs_dir[$comptador] ne $contigs_dir[$comptador-1]){ #escriure un contig nou
		if ($coordbo[0] ne "empty"){
			$start = $coordbo[0] - 100000;	#per si abans haviem detectat un hit que ja s'ha acabat, fem print del start i lentgh del hit
			$length = ($coordbo[1] + 100000) - $start;
			if($start < 0){
				$start = 1;		
			}			
		$starts[$n] = $start;
		$lengths[$n] = $length;
		$n++;
		}
		$contigs[$n] = $all_contigs[$comptador];
	
		$coordbo[0] = $all_starts[$comptador];
		$coordbo[1] = $all_ends[$comptador];
	}

	elsif ($contigs_dir[$comptador] eq $contigs_dir[$comptador-1]){ #afegir dades al contig

		if ($coordbo[0] < $all_starts[$comptador] && $all_ends[$comptador] > $coordbo[1]){	#si el hit esta despres del que ja hem definit
			$coordbo[1] = $all_ends[$comptador];						#allarguem per la dreta
		}
		elsif ($all_starts[$comptador] < $coordbo[0] && $coordbo[1] > $all_ends[$comptador]){	#si el hit esta abans del que ja hem definit
			$coordbo[0] = $all_starts[$comptador];						#allarguem per l'esquerra
		}	
	}
}

$start = $coordbo[0] - 100000;
$length = ($coordbo[1] + 100000) - $start;
	if ($start < 0){
		$start = 1;
	}
$starts[$n] = $start;
$lengths[$n] = $length;

my $scalar = scalar(@contigs);

#print"contigs @contigs starts @starts lengths @lengths\n";

close (BLAST);

#print"ha acabat MarinaMarta\n";


######################FASTASUBSEQ ETC ETC TCOFFEE########################




#loop per a cada contig: farem subseq, exonerate, genewise, Tcoffee
#tot va aqui dins!!!!!

#el $i son el nombre de contigs començant per zero

for(my $i=0; $i<$scalar; $i++){

print"loop i: $i\n";
system("fastafetch $genome /cursos/20428/BI/genomes/2016/Mus_spretus/genome.index $contigs[$i] > $directory/$dir/contig_$contigs[$i]_$i.fa");


#contar nombre de bases del contig




##########nou: comptar llargada contig a partir del fastalength

system("fastalength $directory/$dir/contig_$contigs[$i]_$i.fa > $directory/$dir/length_$contigs[$i]_$i.fa"); #fastalength et posa llargada i nom contig separats per espai

my $cont_length;

open (CONTIG, '<', "$directory/$dir/length_$contigs[$i]_$i.fa");
while(<CONTIG>){
	my @prova = split(/ /,$_); #definir vector, posicio 0 es llargada, posicio 1 nom contig
	$cont_length = $prova[0]; #definir variable amb la largada contig

}
close (CONTIG);



#definir nou lentgh si sobrepassa el propi length del contig

if ($starts[$i] + $lengths[$i] > $cont_length){
	$lengths[$i] = $cont_length - $starts[$i];
}



#ara farem fastasubseq etc etc hehe
system("fastasubseq $directory/$dir/contig_$contigs[$i]_$i.fa $starts[$i] $lengths[$i] > $directory/$dir/subseq_$contigs[$i]_$i.fa");


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

print ("Performing exonerate\n");

system("exonerate -m p2g --showtargetgff -q $directory/$dir/prot_clean.fa -t $directory/$dir/subseq_$contigs[$i]_$i.fa > $directory/$dir/raw_exonerate_$contigs[$i]_$i.fa ");

#system("exonerate -m p2g --showtargetgff --exhaustive yes -q $directory/$dir/prot_clean.fa -t $directory/$dir/subseq_$contigs[$i].fa > $directory/$dir/raw_exonerate_$contigs[$i].fa ");

system("egrep -w gene $directory/$dir/raw_exonerate_$contigs[$i]_$i.fa > $directory/$dir/genes_exonerate_$contigs[$i]_$i.gff");
system("egrep -w exon $directory/$dir/raw_exonerate_$contigs[$i]_$i.fa > $directory/$dir/all_exonerate_$contigs[$i]_$i.gff");

##### per cadascun dels gens predits per exonerate
#Percal Miki

system("cut -f 4 $directory/$dir/genes_exonerate_$contigs[$i]_$i.gff > $directory/$dir/exonerate_starts_$contigs[$i]_$i.txt");
system("cut -f 5 $directory/$dir/genes_exonerate_$contigs[$i]_$i.gff > $directory/$dir/exonerate_ends_$contigs[$i]_$i.txt");
# crea un vector @gene_starts que conté el principi de cada exó de cada gen 
my @gene_starts = ("0");
open(EXONERATE_GENES, '<', "$directory/$dir/exonerate_starts_$contigs[$i]_$i.txt") || die("can't open exonerate_starts_$contigs[$i]_$i.txt");
my $c = 0;
#print"c = $c\n";
while(<EXONERATE_GENES>){
	my $c_start = $_;
	chomp($c_start);
	$gene_starts[$c] = $c_start;	
	$c = $c + 1;
}
my $n_genes = scalar(@gene_starts);
close(EXONERATE_GENES);
#print"gene_starts @gene_starts\n";
#print"n_genes $n_genes\n";
# crea un vector @gene_ends que conté el principi de cada exó de cada gen 
my @gene_ends = ("0");
open(EXONERATE_GENES, '<', "$directory/$dir/exonerate_ends_$contigs[$i]_$i.txt") || die("can't open exonerate_ends_$contigs[$i]_$i.txt");
$c = 0;
#print"c = $c\n";
while(<EXONERATE_GENES>){
	my $c_end = $_;
	chomp($c_end);
	$gene_ends[$c] = $c_end;	
	$c = $c + 1;
}
close(EXONERATE_GENES);
#print"gene_ends @gene_ends\n";
#############
my $n_g = 0; #posicio en @gene_starts
my @start_ids = (0); #contindra un vector [1,1,1,1,2,2,2,2,2,n,n,n] agrupant per cada contig
for ( $n_g = 0; $n_g < $n_genes; $n_g++){ #per cada gene trobat
open(EXONERATE_RESULTS, '<', "$directory/$dir/all_exonerate_$contigs[$i]_$i.gff") || die("can't open all_exonerate_$contigs[$i]_$i.gff");
my $found = 0; #bandera de que ja has trobat aquell exó
my $e_pos = 0; #comptador d'exons trobats per l'exonerate
while(<EXONERATE_RESULTS>){
	my $exon = $_;
	#print"e_pos $e_pos\n";
	chomp($exon);
	#print"found $found\n";
	my @exon_vars = split(/\t/,$exon);
	my $exon_start = $exon_vars[3];
	my $exon_end = $exon_vars[4];
	#print "exon_start $exon_start\n";
	#print "exon_end $exon_end\n";
	#print "gene_starts $gene_starts[$n_g]\n";
	if ($exon_start == $gene_starts[$n_g] or $exon_end == $gene_ends[$n_g]){	#si el start o l'end coincideixen amb els donats per $n_g			
		$found = 1;		
	}
	if ($found == 1){
		$start_ids[$e_pos] = $n_g + 1;	
	}
	$e_pos++;						
}
close(EXONERATE_RESULTS);
#print "gene_starts @gene_starts\n";
}
#print "start_ids @start_ids\n";
# guarda un fitxer per cada gen



for ($n_g = 0; $n_g < $n_genes; $n_g++){
	my $id = $i.$n_g;
	print "id $id n_g $n_g\n";
	open(GENE, '>', "$directory/$dir/exonerate_final_$contigs[$i]_$id.gff");
	open(RESULTS, '<', "$directory/$dir/all_exonerate_$contigs[$i]_$i.gff") || die("can't open all_exonerate_$contigs[$i]_$i.gff");
	my $e_pos = 0;
	while(<RESULTS>){
		if($start_ids[$e_pos] == $n_g+1){	
			print GENE "$_\n";
		}		
		$e_pos++;
	}		
	close(RESULTS);
	close(GENE);	

	# en cada loop hi ha un gen
	##################

	system("fastaseqfromGFF.pl $directory/$dir/subseq_$contigs[$i]_$i.fa $directory/$dir/exonerate_final_$contigs[$i]_$id.gff > $directory/$dir/prednuc_$contigs[$i]_$id.fa");
	
	system("fastatranslate -f $directory/$dir/prednuc_$contigs[$i]_$id.fa -F 1 > $directory/$dir/predaa_$contigs[$i]_$id.fa");
	
	system ("sed 's/*/X/g' $directory/$dir/predaa_$contigs[$i]_$id.fa > $directory/$dir/clean_predaa_$contigs[$i]_$id.fa"); #substituir * per X

	system("t_coffee $directory/$dir/prot_clean.fa $directory/$dir/predaa_$contigs[$i]_$id.fa > $directory/$dir/tcoffee_exonerate_$contigs[$i]_$id.aln");
	
	##################
	
}












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

print ("Performing genewise\n");

if ($direction[$i] eq "F"){
	system("genewise -pep -pretty -cdna -gff $directory/$dir/prot_clean.fa $directory/$dir/subseq_$contigs[$i]_$i.fa > $directory/$dir/genewise_$contigs[$i]_$i.fa");

}

else{
	system ("genewise -pep -pretty -cdna -gff -trev $directory/$dir/prot_clean.fa $directory/$dir/subseq_$contigs[$i]_$i.fa > $directory/$dir/genewise_$contigs[$i]_$i.fa");

}


#######T-COFFE GENEWISE########

my $genewise_file = "$directory/$dir/genewise_$contigs[$i]_$i.fa";

print ("Creating genewise prediction\n");

open(GENEWISE, '<', $genewise_file) || die("can't open $genewise_file");

my $counter = 0;
while (<GENEWISE>) {
	my @genewise_vector = split(' ', $_);				#separem l'output de genewise per espais
	next if $genewise_vector[0] eq "";
	if ($genewise_vector[0] eq '//' && $counter == 1) {		#quan trobem // (esta al final de la proteina)
		close(OUT);						#tanca el document
		$counter = $counter + 1;				#augmentar $count perque no torni a entrar a aquest if
	}
	elsif ($genewise_vector[0] eq ">$contigs[$i]:subseq($starts[$i],$lengths[$i]).pep") {
		$counter = $counter + 1;
		open(OUT, '>', "$directory/$dir/genewise_pred_$contigs[$i]_$i.fa");
	}
	if ($counter == 1) {
		print OUT $_;
	}	
}
close(OUT);

print ("Performing t-coffee\n");

system ("t_coffee $directory/$dir/prot_clean.fa $directory/$dir/genewise_pred_$contigs[$i]_$i.fa > $directory/$dir/t_coffee_$contigs[$i]_$i.aln"); #fer t-coffee











#no borrar el } d'aqui baix!!!!!!!
}









} #no borrar (tancament del primer while del mon mundial)


closedir(DIR);

