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

my $directory = "./proteins_human";
my $genome = "/cursos/20428/BI/genomes/2017/Spermophilus_dauricus/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

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



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

my $scalar1 = scalar (@direction);

for ($n = 0; $n < $num_hit; $n++){
	$contigs_dir[$n] = $all_contigs[$n].$direction[$n];
}
#print "contigs_dir: @contigs_dir\n";
	$n = 0;

my @final_directions = ("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];
		$final_directions[$n] = $direction[$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);

for(my $i=0; $i<$scalar; $i++){
print"loop i: $i\n";
print "$dir\n";

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

print ("Performing genewise\n");

if ($final_directions[$i] eq "F"){
	system("genewise -pep -pretty -cdna -gff $directory/$dir/prot_clean.fa $directory/$dir/subseq_$contigs[$i]_$i.fa > $directory/$dir/genewise_bo_$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_bo_$contigs[$i]_$i.fa");

}


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

my $genewise_file = "$directory/$dir/genewise_bo_$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_bo_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_bo_pred_$contigs[$i]_$i.fa > $directory/$dir/t_coffee_bo_$contigs[$i]_$i.aln"); #fer t-coffee
}




	}
close (DIR);