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

my $directory = "./humans_prot";
my $genome = "./genome_cool.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];
}
print "contigs_dir: @contigs_dir\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\n starts @starts\n lengths @lengths\n";
print "$scalar\n";
for(my $i=0; $i<$scalar; $i++){
print"loop i: $i\n";
system("fastafetch $genome /cursos/20428/BI/genomes/2017/Spermophilus_dauricus/genome.index  $contigs[$i] > $directory/$dir/contig_$contigs[$i]_$i.fa");
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];
}
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 exon $directory/$dir/raw_exonerate_$contigs[$i]_$i.fa > $directory/$dir/exonerate_final_$contigs[$i]_$i.gff");









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");	

}#hit

}#prot
closedir(DIR);