#!/bin/bash

set -e

#INPUT VARIABLES DEFINITION
known_selenoproteins=$1 #SUPPOSED TO BE THE QUERY PROTEINS
genome_folder=$2 #SUPPOSED TO BE THE FOLDER WHERE THE GENOME IS STORED
outputDir=$3 #SUPPOSED TO BE THE FOLDER TO STORE ALL THE GENERATED RESULTS

[ -d $outputDir ] || mkdir $outputDir

#CHANGES U INTO X AND REMOVES SPECIAL CHARACTERS
cat $known_selenoproteins | tr "U" "X" | tr "@" " " | tr "#" " " | tr "%" " " > $outputDir/query.fa

#EXECUTES THE TBLASTN PROGRAM
tblastn -query $outputDir/query.fa -outfmt 6 -evalue 0.01 -db $genome_folder/genome.fa -out $outputDir/blast.tsv

#EXTRACTION OF THE QUERY ID AND THE SCAFFOLD FOR EACH DIFFERENT QUERY-SCAFFOLD HIT IN TBLASTN
query_contig=$(cut -f 1,2 $outputDir/blast.tsv | sort -u | awk '{ print $1"/"$2 }')

#COMMAND TO MAKE THE HEADING OF THE RESULTS TABLE
echo "QUERY ID\tSCAFFOLD\tSTART\tEND\tSTRAND\tNum_EXONS" > $outputDir/resultstable.csv

for hit in $query_contig
do
  query=$(echo $hit | cut -d "/" -f 1)
  contig=$(echo $hit | cut -d "/" -f 2)

  [ -d $outputDir/$query$contig ] || mkdir $outputDir/$query$contig
  cd $outputDir/$query$contig

#ESTABLISHMENT OF THE STRAND WHERE EACH HIT IS FOUND
  sense1=$(grep -E "$query.$contig" ../blast.tsv | cut -f 9 | head -n1)
  sense2=$(grep -E "$query.$contig" ../blast.tsv | cut -f 10 | head -n1)

  if [ $sense1 -gt $sense2 ]
  then
    sense="Reverse"
  else
    sense="Forward"
  fi

#DEFINITION OF THE BEGINNING AND THE END OF EACH GROUP OF HITS TO DETERMINE THE LENGTH AND THE GENOMIC COORDINATES OF EACH GENE
  min=$(grep -E "$query.$contig" ../blast.tsv | awk '{ inicio=$9>$10?$10:$9; final=$9>$10?$9:$10; print inicio, final }' | cut -d " " -f 1 | sort | head -n1)
  max=$(grep -E "$query.$contig" ../blast.tsv | awk '{ inicio=$9>$10?$10:$9; final=$9>$10?$9:$10; print inicio, final }' | cut -d " " -f 2 | sort | tail -n1)

#FASTAFETCH COMMAND TO EXTRACT THE SCAFFOLD FOR EACH HIT
  [ -f scaffold.fa ] || fastafetch $genome_folder/genome.fa $genome_folder/genome.index $contig > scaffold.fa

#DEFINITION OF THE START AND THE LENGTH OF THE FASTASUBSEQ FILE ACCORDING TO THE DIMENSIONS OF EACH SCAFFOLD
  scaffoldlength=$(grep -v ">" scaffold.fa | wc | awk '{print $3-$1}')

  if [ $min -gt 50000 ]
  then
    start=$(( $min - 50000 ))
  else
    start=1
  fi

  longlenght=$(($start + 150000))
  if [ $longlenght -gt $scaffoldlength ]
  then
    length=$(($scaffoldlength - $start))
  else
    length=150000
  fi

#EXECUTION OF THE FASTASUBSEQ COMMAND TO EXTRACT THE REGION OF INTEREST FOR EACH GROUP OF HITS
  [ -f fastasubseq.fa ] || fastasubseq scaffold.fa $start $length > fastasubseq.fa

#AWK PROGRAM TO EXTRACT THE SEQUENCE OF THE CURRENT QUERY PROTEIN IN THE WHOLE QUERY FILE
  awk -v id="$query" '
    $0~">" { recoja=0 }
    $0~id { recoja=1 }
    { if (recoja==1) print $0 }
  ' ../query.fa > proteintosearch.fa

#EXECUTION OF THE PREDICTION COMMANDS: EXONERATE AND EGREP
  [ -f exonerate.out ] || exonerate -m p2g --showtargetgff -n 1 -q proteintosearch.fa -t fastasubseq.fa > exonerate.out

  [ -f onlyexons.gff ] || egrep -w exon exonerate.out > onlyexons.gff | sleep .1
  numberofexons=$(wc -l onlyexons.gff | cut -d " " -f 1)

#IF COMMAND TO PERFORM T-COFFEE ANALYSIS ONLY WHEN EXONERATE HAS MADE A SIGNIFICANT PREDICTION
  if [ -s onlyexons.gff ]
    then
      [ -f predictedntseq.nuc ] || fastaseqfromGFF.pl fastasubseq.fa onlyexons.gff > predictedntseq.nuc
      [ -f predictedaaseq.fa ] || fastatranslate -f predictedntseq.nuc -F 1 | tr "*" "X" > predictedaaseq.fa
      [ -f tcoffee.fa ] || t_coffee proteintosearch.fa predictedaaseq.fa > tcoffee.fa

      if [ $sense = "Reverse" ]
      then
        [ -f genewise.fa ] || genewise -pep -pretty -cdna -gff -trev proteintosearch.fa fastasubseq.fa > genewise.fa
      else
        [ -f genewise.fa ] || genewise -pep -pretty -cdna -gff proteintosearch.fa fastasubseq.fa > genewise.fa
      fi

#ELSE COMMAND TO PERFORM ONLY GENEWISE PREDICTION WHEN EXONERATE HAS NOT MADE A SIGNIFICANT PREDICTION
  else

    if [ $sense = "Reverse" ]
    then
      [ -f genewise.fa ] || genewise -pep -pretty -cdna -gff -trev proteintosearch.fa fastasubseq.fa > genewise.fa
    else
      [ -f genewise.fa ] || genewise -pep -pretty -cdna -gff proteintosearch.fa fastasubseq.fa > genewise.fa
    fi

  fi

#COMMAND TO MAKE A NEW LINE IN THE RESULTS TABLE STORING THE VALUES OF THE MENTIONED VARIABLES
    echo "$query\t$contig\t$min\t$max\t$sense\t$numberofexons" >> $outputDir/resultstable.csv

  cd ..
  cd ..

done
