
echo "Enter the fasta file of protein and hit enter"
read protein #bash accept the input

prefix="$(echo ${protein}|cut -f1 -d "_")"

#PROTEIN'S DIRECTORY 
mkdir -p ${prefix}
cp ${protein} ${prefix}/
cd ${prefix}

#TBLASTN
echo "Performing tblastn"
echo "tblastn -outfmt 6 -query ${protein} -db /mnt/NFS_UPF/soft/genomes/2019/Varanus_komodoensis/genome.fa -out ${prefix}_blast.txt"
tblastn -outfmt 6 -query ${protein} -db /mnt/NFS_UPF/soft/genomes/2019/Varanus_komodoensis/genome.fa -out ${prefix}_blast.txt
echo "done"

echo "Filtering blast output"
awk '{ if(($3 > 50) && ($11 < 0.00001)) { print }}' ${prefix}_blast.txt > ${prefix}_blast_filtered.txt
echo "Done"


echo "Below is the filtered blast output"
cat ${prefix}_blast_filtered.txt


echo "Enter the 'scaffold_name' 'start' 'end' of the region you want to analyze. Otherwise type exit to quit" #write spaces between
read region_info


while [ "${region_info}" != "exit" ];
do
SCAFFOLD="$(echo ${region_info} | cut -f 1 -d " ")"
START="$(echo ${region_info} | cut -f 2 -d " ")"
END="$(echo ${region_info} | cut -f 3 -d " ")"

#FASTAFETCH
echo "Running fastafetch"
echo "fastafetch /mnt/NFS_UPF/soft/genomes/2019/Varanus_komodoensis/genome.fa /mnt/NFS_UPF/soft/genomes/2019/Varanus_komodoensis/genome.index ${SCAFFOLD} > ${SCAFFOLD}.fa"
fastafetch /mnt/NFS_UPF/soft/genomes/2019/Varanus_komodoensis/genome.fa /mnt/NFS_UPF/soft/genomes/2019/Varanus_komodoensis/genome.index ${SCAFFOLD} > ${SCAFFOLD}.fa
echo "Done"

SCAFFOLD_LEN="$(fastalength ${SCAFFOLD}.fa | cut -f 1 -d " ")"

#NEW LIMITS
#FORWARD
if (( $START < $END )); then
	echo "${START} ${END} Region is forward"
					
	START_COORD=$(($START-50000))
	if (($START_COORD < 0)); then
		START_COORD=1 
	fi
		
	END_COORD=$(($END+50000))		
	if (($END_COORD > $SCAFFOLD_LEN)); then
		END_COORD=$SCAFFOLD_LEN
	fi
#REVERSE
else
	echo "${START} ${END} Region is reverse"
	END_COORD=$(($START + 50000))
	if (($END_COORD > $SCAFFOLD_LEN)); then
		END_COORD=$SCAFFOLD_LEN
	fi
	START_COORD=$(($END - 50000))
	if (($START_COORD < 0)); then
		START_COORD=1
	fi

fi

LENGTH=$(($END_COORD-$START_COORD))

#FASTASUBSEQ: extraction
echo "Running fastasubseq"
echo "fastasubseq ${SCAFFOLD}.fa ${START_COORD} ${LENGTH} > ${SCAFFOLD}_${START_COORD}_${END_COORD}.fa"
fastasubseq ${SCAFFOLD}.fa ${START_COORD} ${LENGTH} > ${SCAFFOLD}_${START_COORD}_${END_COORD}.fa 
echo "done"	

#EXONERATE: annotation
echo "Running exonerate"
echo "exonerate -n 1 -m p2g --showtargetgff -q ${protein} -t ${SCAFFOLD}_${START_COORD}_${END_COORD}.fa | egrep -w exon > ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons.gff"
exonerate -n 1 -m p2g --showtargetgff -q ${protein} -t ${SCAFFOLD}_${START_COORD}_${END_COORD}.fa | egrep -w exon > ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons.gff
echo "done"

#FASTASEQ:cDNA + FASTATRANSLATE:aa
if [ ! -s ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons.gff ]; then #! -s is empty
	echo "Exonerate output is empty! Removing the region file and gff file!"
	rm ${SCAFFOLD}_${START_COORD}_${END_COORD}.fa
	rm ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons.gff #stop?
else

	echo "Converting GFF to Fasta"
	echo "fastaseqfromGFF.pl ${SCAFFOLD}_${START_COORD}_${END_COORD}.fa ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons.gff > ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons_nucl.fa"
	fastaseqfromGFF.pl ${SCAFFOLD}_${START_COORD}_${END_COORD}.fa ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons.gff > ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons_nucl.fa
	echo "done"
	
	
	echo "Translating nucleotides to protein and substituting * by X"
	echo "fastatranslate -f ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons_nucl.fa -F 1 > ${SCAFFOLD}_${START_COORD}_${END_COORD}_aa.fa"
	fastatranslate -f ${SCAFFOLD}_${START_COORD}_${END_COORD}_exons_nucl.fa -F 1 > ${SCAFFOLD}_${START_COORD}_${END_COORD}_aa.fa
	#X substitution
	sed 's/*/X/g' ${SCAFFOLD}_${START_COORD}_${END_COORD}_aa.fa > ${SCAFFOLD}_${START_COORD}_${END_COORD}_aa_U_repl.fa
	echo "done"
	
	echo "Running T-COFFEE"
	echo "t_coffee ${protein} ${SCAFFOLD}_${START_COORD}_${END_COORD}_aa_U_repl.fa > ${SCAFFOLD}_${START_COORD}_${END_COORD}_tcoffee"
	t_coffee ${protein} ${SCAFFOLD}_${START_COORD}_${END_COORD}_aa_U_repl.fa > ${SCAFFOLD}_${START_COORD}_${END_COORD}_tcoffee
	echo "done"

fi

#ENDING
echo "Analysis for ${region_info} is finished"
echo "The filtered blast output was "

#REPEAT FOR OTHER SCAFFOLDS
cat ${prefix}_blast_filtered.txt

echo "Type another region information or type exit to quit"

read region_info
done

cd ..
