#!/bin/bash

set -e


known_selenoproteins=$1
genome_folder=$2
outputDir=$3

[ -d $outputDir ] || mkdir $outputDir

cat $known_selenoproteins | tr "U" "X" | tr "@" " " | tr "#" " " | tr "%" " " > $outputDir/query.fa

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

query_contig=$(cut -f 1,2 $outputDir/blast.tsv | sort -u | awk '{ print $1"/"$2 }')

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

  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

  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)

  [ -f scaffold.fa ] || fastafetch $genome_folder/genome.fa $genome_folder/genome.index $contig > scaffold.fa

  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

  [ -f fastasubseq.fa ] || fastasubseq scaffold.fa $start $length > fastasubseq.fa

  awk -v id="$query" '
    $0~">" { recoja=0 }
    $0~id { recoja=1 }
    { if (recoja==1) print $0 }
  ' ../query.fa > proteintosearch.fa

  [ -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 [ -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

    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

  cd ..
  cd ..

  echo "$query\t$contig\t$min\t$max\t$sense\t$numberofexons" >> $outputDir/resultstable.csv

done
