#!/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/green_proteins/allseqgreen	#Anole lizard Selenoproteins from SelenoDB
PATH_hits="~/Escriptori/vipera/green_proteins/allhitsgreen/"	#Hits obtained with blast software

#To change directory to the hit's path and to define the hits variable
cd $PATH_hits
hits="eEFsec(1)(hit1) GPx8(hit2) MsrA(hit5) SecS(hit1) SelR3(hit1) eEFsec(2)(hit1) GPx(hit1) Sel0(1)(hit1) Sel0(2)(hit1) eEFsec(3)(hit1) MsrA(hit1) PSTK(hit1) SelU2(hit1) GPx4(hit1) MsrA(hit2) SBP2(1)(hit1) SPS1(hit1) MsrA(hit3) SBP2(2)(hit1) SelI(hit1) SPS1(hit2)
GPx8(hit1) MsrA(hit4) SECp43(hit1) SelR2(hit1)"

#To declare a dictonary with each hit of each protein
declare -A hits_proteins
hits_proteins=(["eEFsec(1)(hit1)"]="eEFsec(1)" ["GPx8(hit2)"]="GPx8" ["MsrA(hit5)"]="MsrA" ["SecS(hit1)"]="SecS" ["SelR3(hit1)"]="SelR3" ["eEFsec(2)(hit1)"]="eEFsec(2)" ["GPx(hit1)"]="GPx" ["Sel0(1)(hit1)"]="Sel0(1)" ["eEFsec(3)(hit1)"]="eEFsec(3)" ["MsrA(hit1)"]="MsrA" ["PSTK(hit1)"]="PSTK" ["SelU2(hit1)"]="SelU2" ["GPx4(hit1)"]="GPx4" ["MsrA(hit2)"]="MsrA" ["SBP2(1)(hit1)"]="SBP2(1)" ["Sel0(2)(hit1)"]="Sel0(2)" ["SPS1(hit1)"]="SPS1" ["MsrA(hit3)"]="MsrA" ["SBP2(2)(hit1)"]="SBP2(2)" ["SelI(hit1)"]="SelI" ["SPS1(hit2)"]="SPS1"
["GPx8(hit1)"]="GPx8" ["MsrA(hit4)"]="MsrA" ["SECp43(hit1)"]="SECp43" ["SelR2(hit1)"]="SelR2")


#To enter a loop for each hit
for hit in $hits
do
	#To make a directory for each hit
	mkdir -p $PATH_vipera/green_proteins/resultsgreen/$hit/
	cd $PATH_vipera/green_proteins/resultsgreen/$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=1
	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
