# !/usr/bin/perl -w
# El programa en perl és elaborat per: Eloi Mallofré, Isabel Pérez, Aida Ripoll i Carla San Martín.
# Aquest script analitza si una query donada es troba dins del genoma de Nanorana parkeri.
# Per algun altre dubte, poseu-vos en contacte amb algun dels membres del grup.

 
my $master;
my $mfile;
my $outdir;

print "Programa en perl d'anotació del genoma de Nanorana parkeri\n";


### Pel funcionament correcte del PROGRAMA necessitem haver creat una carpeta "master" que contingui:
# query creada: amb les U canviades per X i sense cap caràcter final.
# programa.sh
# programa.pl


### Creem un loop per analitzar totes les queries de la carpeta "master":

print "Nom de la carpeta $master on hi ha les queries\n:";
$master = <STDIN>;
$master = "./master";
chomp $master;

opendir(DIR, $master) or die $!;
while (my $mfile = readdir(DIR)){
	if($mfile=~/.fa$/){
		print "$mfile\n"; #arxiu fasta ("nomquery".fa)

		my ($outdir) = $mfile=~/(.+).fa/;
		print "$outdir\n"; 

		system("mkdir $outdir"); #directori ("nomquery) 

		### BLASTALL: obtindrem un blast simplificat (només donarà el nom dels scaffolds amb un e-value <0.05) que es guardarà en una carpeta amb el nom de la query.
		# Si el blast ja està creat, no el tornarà a fer. 
		if (-e "./$outdir/$mfile.blast.fa"){
			print "el blast ja existeix\n";
		} else {
			system ("blastall -p tblastn -i $master/$mfile -d /cursos/BI/genomes/2015/Nanorana_parkeri/genome.fa -o ./$outdir/$mfile.blast.fa -e 0.05 -m 8"); #("nomquery".blast.fa)
			print "fitxer ./$outdir/$mfile.blast.fa creat\n"; 
		}

		### FASTAINDEX: Per agilitzar el procés es crea un índex del genoma.
		# Ens assegurem que l'índex encara no existeixi, ja que només cal crear-lo una vegada.
		if (-e "./Nanorana_parkeri.index") {
		    print "l'índex ja existeix\n";
		} else {
		    system ("fastaindex /cursos/BI/genomes/2015/Nanorana_parkeri/genome.fa Nanorana_parkeri.index"); #("nomgenoma".index)
		    print "fitxer Nanorana_parkeri.index creat\n";
		}
		
		### FASTAFETCH: A partir de l'índex del genoma creat s'extreu la seqüència d'un scaffold concret.
		my $i;
		my @v;

		system ("cut ./$outdir/$mfile.blast.fa -s -f 2 > ./$outdir/$mfile.hits.txt"); # Es crea un document amb el nom de cada hit (hits.txt) 

		my $i = 0;
		open (NOMSBLAST,"<./$outdir/$mfile.hits.txt");
		my @v = <NOMSBLAST> ;
		close (NOMSBLAST);

		for (@v) {
			chomp;
			unless (-e "./$outdir/${_}.fetch.fa") {
				system ("fastafetch /cursos/BI/genomes/2015/Nanorana_parkeri/genome.fa Nanorana_parkeri.index '$_' > './$outdir/${_}.fetch.fa' "); #("nomcromosoma".fetch.fa)
				print "fitxer ./$outdir/${_}.fetch.fa creat a la carpeta\n";
				$i = $i + 1;
			}
		}


		### FASTASUBSEQ: A partir del fitxer blast generat prèviament, agafem els valors de les columnes 8 (inici hit) i 9 (final hit).
		my %hits;
		my $start_hit;
		my $stop_hit;
		my $start;
		my $length;

		system ("mkdir -p ./$outdir/subseq");
		print "carpeta subseq creada\n";
	
		my %lengths;
		open(LENGTHS,"</cursos/BI/genomes/2015/Nanorana_parkeri/genome.lengths");
		while (my $row = <LENGTHS>) { 
			chomp;
			my @line = split(" ",$row);
			$lengths{$line[1]} = $line[0];
		}

		close (LENGTHS);

		open (BLAST, "<./$outdir/$mfile.blast.fa");
		while (my $row = <BLAST>) {
			chomp;
			my @blast = split ("\t",$row);
			# assignem la posició inici i final segons el sentit de la seqüència sigui positiu o negatiu.
			if ($blast[8]<$blast[9]) {
				$start_pos = $blast[8];
				$stop_pos = $blast[9];
			} else {
				$start_pos = $blast[9];
				$stop_pos = $blast[8];
			}

			# actualitzar o iniciar les posicions d'inici i final.
			if (-e ($hits{$blast[1]})) {
				if ($hits{$blast[1]}{"START"} > $start_pos) {
					$hits{$blast[1]}{"START"} = $start_pos;
				}
				if ($hits{$blast[1]}{"STOP"} < $stop_pos) {
					$hits{$blast[1]}{"STOP"} = $stop_pos;
				}
			} else {
				$hits{$blast[1]}{"START"} = $start_pos;
				$hits{$blast[1]}{"STOP"} = $stop_pos;
			}
			print "$row\n";

			# per cada hit s'assigna un marge de 50000 nucleòtids als extrems. 

			# condicions a seguir en cas de que: la variable $start sigui inferior a 0 i/o la variable $stop sigui superior a la llargada de l'scaffold.
			$start = $start_pos - 50000;
			$length = ($stop_pos - $start_pos) + 100000;
			$stop = $start + $length;
			if ($start < 0) {
				$start = 0;
				$length = ($stop_pos-$start_pos) +50000;
			}
			if ($stop > $lengths{$blast[1]}) {
				$stop = $lengths{$blast[1]};
				$length = $stop - $start;
			}
	
			system ("fastasubseq './$outdir/${blast[1]}.fetch.fa' $start $length > './$outdir/subseq/${blast[1]}_subseq.fa' ");
	
#("nomcromosoma".subseq.fa) ./$outdir/${_}.fetch.fa vs './$outdir/${blast[1]}.fetch.fa'

			print "fitxers subseq creats a la carpeta $outdir/subseq\n";

		} #tanca while (BLAST)
		close (BLAST);

		# crear carpetes pels seguents passos: exonerate, genewise i tcoffee
		system ("mkdir -p ./$outdir/exonerate");
		system ("mkdir -p ./$outdir/exonerate/gff");
		system ("mkdir -p ./$outdir/exonerate/cDNAgff");
		system ("mkdir -p ./$outdir/exonerate/genespredict");
		system ("mkdir -p ./$outdir/exonerate/protspredict");
		system ("mkdir -p ./$outdir/exonerate/protsdef");
		system ("mkdir -p ./$outdir/genewise");
		system ("mkdir -p ./$outdir/tcoffee");


		opendir (SUBSEQ, "$outdir/subseq/");
		while (my $file = readdir (SUBSEQ)) {
			next if ($file =~ m/^\./);
			chomp $file;
			my $protein_file = $file;
			$protein_file =~ s/\_subseq.fa//g;
			chomp $protein_file;

			### EXONERATE: Eïna genèrica per fer l'alineament de seqüències. 
			system ("exonerate -m p2g --showtargetgff -q '$master/$mfile' -t './$outdir/subseq/$file' > './$outdir/exonerate/gff/$protein_file.gff' ");

			system ("egrep -w exon './$outdir/exonerate/gff/$protein_file.gff' > './$outdir/exonerate/cDNAgff/$protein_file.cDNA.gff' ");

			print "fitxer exonerate creat a la carpeta $outdir/exonerate\n";

			### FASTASEQfromGFF:
			system ("fastaseqfromGFF.pl './$outdir/subseq/$file' './$outdir/exonerate/cDNAgff/$protein_file.cDNA.gff' > './$outdir/exonerate/genespredict/$protein_file.genespredict.fa' "); 

			print "fitxer fastaseqfromgff creat a la carpeta $outdir/exonerate\n";

			### FASTATRANSLATE:
			# s'escull unicament el "reading frame" 1 (-F1) ja que s'ha comprovat empiricament que l'exonerate mostra com a frame correcte l'1.
			system ("fastatranslate -f './$outdir/exonerate/genespredict/$protein_file.genespredict.fa' -F 1 > './$outdir/exonerate/protspredict/$protein_file.protspredict.fa' ");

			print "fitxer fastatranslate creat a la carpeta $outdir/exonerate\n";
			
			# obrim la carpeta on hi ha els arxius translate i canviem les X de la seqüència d'aminoàcids de la proteïna predita per U. 
			open (TRANS, "$outdir/exonerate/protspredict/$protein_file.protspredict.fa");
			my $name = "./$outdir/exonerate/protsdef/$protein_file.protsdef.fa";

			open(OUTPUT, ">", $name);
			while(<TRANS>) {
				if (/^>/){
					print OUTPUT $_;
				} elsif ($_ !~/^>/){
					s/\*/U/g;
					print OUTPUT $_;
				}
			}
			close (OUTPUT);

			### TCOFFEE:
			system("t_coffee './master/$mfile' './$outdir/exonerate/protsdef/$protein_file.protsdef.fa' > './$outdir/tcoffee/$protein_file.tcoffee.fa' "); 
			print "fitxer t-coffee creat a la carpeta $outdir/tcoffee\n";

			### GENEWISE:
			system ("genewise -cdna -pep -pretty -gff -both './master/$mfile' './$outdir/subseq/$file' > './$outdir/genewise/$protein_file.genewise.fa' ");
			print "fitxer genewise creat a la carpeta $outdir/genewise\n";


		} #tanca el while del readdir SUBSEQ (abans exonerate)

	} #tanca el primer if del programa
} #tanca el while (loop grandullon)

