#!/bin/bash

#To export all software required (exonerate, fastaseqfromGFF.pl and t_coffee)
export PATH=/cursos/BI/soft/exonerate/x86_64/bin:$PATH	#exonerate
export PATH=/cursos/BI/bin:$PATH	#fastaseqfromGFF.pl/cursos/BI/soft/t_coffee/x86_64/bin:$PATH.
export PATH=/cursos/BI/soft/t_coffee/x86_64/bin:$PATH	#t_coffee

#To define the paths of the inputs and outputs
PATH_vipera=~/Escriptori/vipera	#Directory that will contained all the outputs directories and files
PATH_files=/cursos/BI/genomes/2015/Vipera_berus_berus	#Vipera berus genome, genome.index and genome.lenghts
PATH_proteins=~/Escriptori/vipera/red_proteins/allseqred	#Anole lizard Selenoproteins from SelenoDB
PATH_hits="~/Escriptori/vipera/red_proteins/allhitsred/"	#Hits obtained with blast software

#To change directory to the hit's path and to define the hits variable
cd $PATH_hits
hits="Sel15(hit1) SelP1(hit1) SPS(hit1) DI2(hit1) GPx(hit1) SelH(hit1) SelP2(hit1) DI1(hit1) DI2(hit2) GPx(hit2) SelI(hit1) SelR1(hit1) TR3(hit1) DI1(hit2) GPx2(hit1) GPx(hit3) SelN(hit1) SelT(hit1) TR(hit1) DI1(hit3) GPx3(hit1) SelU1(hit1)"


#To declare a dictonary with each hit of each protein
declare -A hits_proteins
hits_proteins=(["DI1(hit2)"]="DI1" ["DI1(hit3)"]="DI1" ["DI1(hit1)"]="DI1" ["DI2(hit2)"]="DI2" ["DI2(hit1)"]="DI2" ["GPx(hit2)"]="GPx" ["GPx(hit3)"]="GPx" ["GPx2(hit1)"]="GPx2" ["GPx3(hit1)"]="GPx3" ["GPx(hit1)"]="GPx" ["SelH(hit1)"]="SelH" ["SelI(hit1)"]="SelI" ["SelN(hit1)"]="SelN" ["SelP1(hit1)"]="SelP1" ["SelP2(hit1)"]="SelP2" ["SelR1(hit1)"]="SelR1" ["SelT(hit1)"]="SelT" ["SelU1(hit1)"]="SelU1" ["SPS(hit1)"]="SPS" ["TR3(hit1)"]="TR3" ["TR(hit1)"]="TR" ["Sel15(hit1)"]="Sel15")
	
#To enter a loop for each hit
for hit in $hits
do
	#To make a directory for each hit
	mkdir -p $PATH_vipera/red_proteins/resultsred/$hit/
	cd "$PATH_vipera/red_proteins/resultsred/$hit/" 

	#To define the sequence ($seq variable) from where the hit comes from
	seq=$(cat $PATH_hits/$hit | cut -f1)

	#To run the fastafetch in order to obtain the sequence of the hit
	fastafetch $PATH_files/genome.fa $PATH_files/genome.index $seq > $hit.seq.fetch	
	
	##To define the start, the end and the lenght variables for the fastasubseq software:
	#1. To cut the variables $start and $end from the hits file and $lenghtseq from the genome.lenghts file
	start1=$(cat $PATH_hits/$hit | cut -f2) 
	end1=$(cat $PATH_hits/$hit | cut -f3)	
	lengthseq=$(cat $PATH_files/genome.lengths | egrep $seq | cut -d ' ' -f1)
	
	#2. To take 50000 nucleotides more at the start and at the end in order to ensure the region kept contains all the selenoprotein gene 
	start2=$(($start1-50000))
	end=$(($end1+50000))
	length=$(($end-$start2))

	#3. To make sure both the start and the end take possible positions of the sequence, they're not out of range
	if (( $start2>0)); then
		start=$start2
	else
		start=0
	fi

	if [ "$end">"$lengthseq" ]; then
		lengfinal=$(($lengthseq-$start))
	fi

	#To run the fastasubseq in order to obtain the scaffold?? (genomic.fa) of each hit
	fastasubseq  $hit.seq.fetch $start $lengfinal > $hit.genomic.fa

	#To run the exonerate software in order to predict the gene structure
	exonerate -m p2g --showtargetgff -q $PATH_proteins/${hits_proteins["$hit"]} -t $hit.genomic.fa > $hit.exonerate.gff	

	#To keep only the exon region from the exonerate prediction
	egrep -w exon $hit.exonerate.gff > $hit.exon.gff
	
	#To change the file format from gff to fasta	
	fastaseqfromGFF.pl $hit.genomic.fa $hit.exon.gff > $hit.pred.nuc.fa #gene prediction in fasta format (only exons)

	fastaseqfromGFF.pl $hit.genomic.fa $hit.exonerate.gff > $hit.SECIS.nuc.fa #gene prediction in fasta format (introns and exons)
	

	#To run the fastatranslate software to translate the nucleotide sequence obtained with exonerate to the aminoacids sequence
	fastatranslate -F 1 $hit.pred.nuc.fa > $hit.pred.aa.fa 

	#To run t_coffee to compare the aminoacid sequence predicted from the Vipera berus genome with the Anole lizard one
	t_coffee $hit.pred.aa.fa $PATH_proteins/${hits_proteins["$hit"]} > $hit.tcoffee 

done
