#!/usr/bin/perl -w 

use strict;
use warnings;

#write files
#abans d'executar este programa el tblastn ja s'ha realitzat a mà des del shell.



#################################### AIXÒ ÉS INTRODUCCIÓ DE DADES DEL TBLASTN###################################



my $prot; #proteïnes que anem a analitzar; s'introduirà cada vegada que comencem una nova proteïna (38 en total). 
my $scaff; #scaffold que volem analitzar. 
my $pi; #posició inicial de cada scaffold.
my $pf; #posició final de cada scaffold.
my $llarg; #llargada de l'scaffold. 
my @scafflength; #vector de les llargades de tots els scaffolds
my $linia; #que es el nostre row en aquest cas 
my $sl; #llargada de l'escaffold en qüestió
my $length; #llargada del vector scafflength

#print "Which is the protein you want to analyse?\n ";
#$prot = <STDIN>; 
$prot = $ARGV[0];
chomp($prot);


#print "Which is the specific scaffold of your protein you want to analyse?\n";
#$scaff = <STDIN>; 
$scaff = $ARGV[1];
chomp($scaff);

#system ("mkdir $prot.$scaff");

# Determinació de la posició inicial del scaffold, la posició final i la llargada. Li sumem 100.000 per darrere i li'n restem 50.000 per davant. 

#print "Which is the start position of the scaffold?\n";
#$pi = <STDIN>;
$pi = $ARGV[2];
$pi = $pi - 50000;
chomp($pi);

#print "Which is the end position of the scaffold?\n";
#$pf = <STDIN>;
$pf = $ARGV[3];
$pf = $pf + 50000;
chomp($pf);

if ($pi < 0 ){
	$pi = 1;
}


open(IN, "<", "scaffoldlength.fa"); #variable in de fitxer, es una direcció a este fitxer del teu path per al length

while ($linia = <IN>) {		

chomp $linia; # eliminem el \n de $row

@scafflength = split(/\ /, $linia);

#$length = scalar (@scafflength);

		#if ($length >= 1){		    
			if ($scafflength[1] eq "$scaff") {

    			$sl = $scafflength[0];
			}
		#}
}
close(IN);

print "This is the length of the scaff $sl \n"; 

if ($pf > $sl) {
	$pf = $sl;
}                            
print "This is the last position $pf\n";

$llarg = ($pf)-($pi);

################################ AIXÒ ÉS FASTAFETCH I FASTASUBSEQ ###########################################

system ("fastafetch /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Urocitellus_parryii/genome.fa /mnt/NFS_UPF/bioinfo/BI/genomes/2018/Urocitellus_parryii/genome.index $scaff > $scaff.scaffold.fa");

print " Congrats! Fastafetch has been properly done\n";

system ("fastasubseq $scaff.scaffold.fa $pi $llarg > $scaff.subseq.fa");  

print " Congrats! Fastasubseq has been properly done\n";



############################# AIXÒ ÉS EXONERATE AUTOMATITZAT ejejejjejje ####################################


system ("exonerate -m p2g --showtargetgff -q ./proteins_human/$prot/$prot.fa -t $scaff.subseq.fa > EXONERATE.raw"); #ojo en el path que la u anira variant


my $i = 0; #i es el número de genes
my $row; #la linea en questio
my @row_array; #vector de cada linia
my $filename = "EXONERATE.raw"; 
my $l;

# fem un grep de les línies que ens interessen content: les que contenen 'exonerate:protein2genome:local'

system("grep 'exonerate:protein2genome:local' $filename > corrected_exonerate.gff");

# obrir el fitxer
open(IN, "<", "corrected_exonerate.gff"); # or die "Could not open file corrected_exonerate.gff!"; 

while ($row = <IN>) {		

	chomp $row; # eliminem el \n de $row

	@row_array = split(/\t/, $row);

	$l = scalar (@row_array);

		if ($l >= 2){
		    
			if ($row_array[2] eq "gene") {

    				if ($i != 0){ close(OUT); } # només a partir del 2n gen que trobem hem de tancar el OUT anterior
		
				# obrir un fitxer per a un gen nou

				$i++; # sumem 1 a la $i 
			
				open(OUT,">" , "$prot.$scaff.gene_$i.gff");
	
			}
	}
	print OUT $row."\n"; # escrivim totes les línies en algun fitxer. Recordem que al fer chomp hem eliminat el \n, així que aquí l'afegim

}
close(OUT); # tanquem l'últim fitxer
close(IN); # hem de tancar el que estavem llegint

# anem per tots els fitxer gene_1, gene_2 ... I fem les annotacions
my $gID = 1; # comptador de gens predts
my $gene_filename; # el nom del fitxer

while ($gID <= $i){

	$gene_filename = "$prot.$scaff.gene_$gID.gff"; # el fitxer que volem treballar

	system ("egrep -w exon $gene_filename > $scaff.$gID.exon.gff");

	print "Analyzing $gene_filename\n";

	system ("fastaseqfromGFF.pl $scaff.subseq.fa $scaff.$gID.exon.gff > $prot.$scaff.$gID.nucl.fa"); # el que fem es traduir de aa a nucleotids; de format gff a fasta (.fa)
	

	print "Congrats! FastaseqfromGFF has been properly done\n";

	system ("fastatranslate $prot.$scaff.$gID.nucl.fa -F 1 > $scaff.$gID.pep.fa"); # s'ha traduit nucleotids a aa

	print "Congrats! Fastatranslate has been properly done\n ";

# Canvi de * per X

	system ("sed 's/*/X/g' $scaff.$gID.pep.fa > $scaff.$gID.pepX.fa");
	
# system #genewise

	system ("t_coffee $scaff.$gID.pepX.fa ./proteins_human/$prot/$prot.fa > ./proteins_human/$prot/$scaff.$gID.tcoffee.txt 2> ./proteins_human/$prot/t_coffee_$gID.stderr.txt"); # s'executa el tcoffee	
	
	system ("cp $scaff.$gID.exon.gff ./proteins_human/$prot/");
	system ("cp $scaff.subseq.fa ./proteins_human/$prot/");
	system ("cp $scaff.$gID.pep.fa ./proteins_human/$prot/");
	system ("cp $scaff.$gID.pepX.fa ./proteins_human/$prot/");
   	
	$gID++; # sumem 1 al $gID

}