#!/usr/bin/perl -w
use strict;
use Data::Dumper qw(Dumper); 

### This program belongs to the Bioinformatics final project. Lovingly coded by Laia López, María Ortuño, Mònica Biarnés and Irene García.  

### ALIAS FOR PROGRAMS AND REMOTE DOCUMENTS  

	my $directory = $ARGV[0]; #getting the path to the folder with the proteins, it is writen on the terminal
	my $tblastn="/mnt/NFS_UPF/soft/ncbi-blast-2.7.1+/bin/tblastn"; # path to Tblastn
	my $exonerate="/mnt/NFS_UPF/soft/exonerate/exonerate-2.2.0-x86_64/bin/exonerate"; #path to Exonerate
	my $genoma_mungo="/mnt/NFS_UPF/soft/genomes/2019/Mungos_mungo/genome.fa"; #path to the Mungos mungo genome
	my $indicegenoma="/mnt/NFS_UPF/soft/genomes/2019/Mungos_mungo/genome.index"; #path to the indexed Mungos mungo genome
	my $lengthgenoma="/mnt/NFS_UPF/soft/genomes/2019/Mungos_mungo/genome.lengths"; #path to the document with scaffold lengths 


### AUTOMATION LOOP 

		print "YOU ARE NOW USING A PROGRAM DESIGNED TO ANNOTATE THE MUNGOS MUNGO GENOME.\n\n";

my @fasta=glob($directory.'/*.fa'); # We create a vector that beholds the documents ended in .fa found in the directory

my $num_docs=scalar(@fasta); # We define the number of proteins we will iterate
		
for (my $p='0'; $p < $num_docs; $p++) { 

	my @query = split(/[\/\.]/,$fasta[$p]); # We split the directory path by "/" and by "." 

	my $array_length=scalar(@query); # We get the length of the array

	my $folder=$query[$array_length-'2']; # We got the file name from the path, the one we will use to make a new folder and name everything

			print "THE SELENOPROTEIN YOU ARE STUDYING IS: $folder \n\n";

	my $directory_out=$directory."/".$folder; #We do create a variable that beholds the directory path

	mkdir $directory_out; #We create the new folder using the folder path where everything is stored and adding the new name

	my $protein_file = "$folder".".fa"; #We define the name of a new document that will contain a modified version of each protein

	system ("sed 's/U/X/g' $fasta[$p] > '$directory_out/$protein_file'"); #We substitute U's for X's in our protein sequences and add this document to the new folder

			print "NOW YOU WILL BE REQUESTED TO ENTER A SERIES OF VALUES THAT SHALL BE PRINTED ON-SCREEN. THE FIRST COLUMN CORRESPONDS TO YOUR HITS' SCAFFOLD. THE STARTING POSITION OF YOUR HIT IS THE LOWEST NUMBER BETWEEN THE SECOND AND THIRD COLUMNS. THE FOURTH COLUMN INDICATES THE E-VALUE OF EACH HIT.\n\n";

	###TBLASTN

	my $tblastn_out="$folder"."_tblastn.txt";

	system("$tblastn -query $directory_out/$protein_file -db $genoma_mungo -outfmt 6 -evalue 0.001 -out '$directory_out/$tblastn_out'"); # We only select hits with an e-value lower than 0.001


	###FASTAFETCH 

	my $cut_out="$folder"."_cut.txt";

	system("cut -f '2','9','10','11' $directory_out/$tblastn_out > '$directory_out/$cut_out'"); #We select interesting TBLASTN output columns 

	my $decision='0'; #We create a loop that lets us choose whether we want to analyse more than one hit for each protein or not 

	my $iteracio='1'; #We use this variable in order to name our outputs so the outputs from hits from the same scaffold don't get overwritted  

	while ($decision=='0'){
	 
			#We define our scaffold:

					print "Your hits are: \n\n";
					
					system ("more '$directory_out/$cut_out'");

					print " \n";

					print "Enter your hit's scaffold: \n";

					my $scaffold=<STDIN>;
					
					chomp($scaffold);

			my $fastafetch_out="$folder"."_"."$scaffold"."_"."$iteracio"."_fastafetch.fa"; 

		system("fastafetch $genoma_mungo $indicegenoma '$scaffold' > '$directory_out/$fastafetch_out'"); 


		###FASTASUBSEQ 

			#We will obtain two sets of data for two sets of documents: the normal ones we will use to further analyze with programs such as T-Coffee, and a shorter version for its use in Seblastian, an algorism that does not permit sequences greater than 100.000.  

			#We define the start location of our hit
				
					print " \n";

					print "Enter the start location of your hit: \n";

					my $inici_hit=<STDIN>;

					chomp($inici_hit);


			my $inici=$inici_hit-'100000';

				if ($inici>='0') {
					$inici=$inici;
						}
				else {$inici='0'}


			my $inici_seblastian=$inici_hit-'50000';

				if ($inici_seblastian>='0') {
					$inici_seblastian=$inici_seblastian;
					}
				else {$inici_seblastian='0'}


			#We define our length

			my $length_scaffold = `grep '$scaffold' '$lengthgenoma' | tr -s ' ' | cut -f 1 -d ' '`;   # This directly saves the length of the scaffold into the variable length_scaffold


			#We check if the whole length we are selecting fits inside of our scaffold

				my $length;
						
				if ($length_scaffold > ($inici+'200000')){
					
					$length='200000';

					}

				else {$length=($length_scaffold-$inici);
							
					}


				my $length_seblastian;
						
				if ($length_scaffold > ($inici_seblastian+'100000')){

					$length_seblastian='100000';

					}

				else {$length_seblastian=($length_scaffold-$inici_seblastian);

					}


		my $subseq_out="$folder"."_"."$scaffold"."_"."$iteracio"."_fsubseq.fa"; 

		my $subseq_out_seblastian="$folder"."_"."$scaffold"."_"."$iteracio"."_fsubseq_seblastian.fa";  

				
		system("fastasubseq '$directory_out/$fastafetch_out' '$inici' '$length' > '$directory_out/$subseq_out'"); 

		system("fastasubseq '$directory_out/$fastafetch_out' '$inici_seblastian' '$length_seblastian' > '$directory_out/$subseq_out_seblastian'");


		###EXONERATE

		my $exonerate_out="$folder"."_"."$scaffold"."_"."$iteracio"."_exonerate.gff";

		system ("exonerate -m 'p2g' --showtargetgff -q $directory_out/$protein_file -t '$directory_out/$subseq_out' | egrep -w 'exon' > '$directory_out/$exonerate_out'");


		###FASTASEQFROMGFF

		my $fseqfromgff_out="$folder"."_"."$scaffold"."_"."$iteracio"."_fseqfromgff.fa";

		system ("fastaseqfromGFF.pl '$directory_out/$subseq_out' '$directory_out/$exonerate_out' > '$directory_out/$fseqfromgff_out'"); 


		###FASTATRANSLATE

		my $translate_out="$folder"."_"."$scaffold"."_"."$iteracio"."_translate.fa";

		system ("fastatranslate -f '$directory_out/$fseqfromgff_out' -F 1 | sed 's/\*/X/g' > '$directory_out/$translate_out'"); 


		###TCOFFEE

		my $tcoffee_out="$folder"."_"."$scaffold"."_"."$iteracio"."_tcoffee.txt"; 

		system("t_coffee $directory_out/$protein_file '$directory_out/$translate_out' > '$directory_out/$tcoffee_out' 2>/dev/null"); 

				print "\n";

		
		#We decide whether we want to keep analyzing the same protein or if we want to change

			print "DO YOU WANT TO KEEP ANALYZING THIS PROTEIN? (YES/NO)\n";

			my $resposta=<STDIN>;
			
			chomp($resposta);

			if ($resposta eq 'YES'){

				$decision='0';

				$iteracio=$iteracio+1;

				print " \n\n\n";

				print "Your protein was: $folder\n\n";

				}

			elsif ($resposta eq 'NO') {

				$decision='1';

					print " \n";

					print "HERE FINISHES THE PROCESS FOR THE SELENOPROTEIN $folder\n";
					
					print " \n";
				}

	}

}

print "YOU HAVE ANALYZED ALL OF THE PROTEINS\n";

print " \n";
