#!/bin/bash 

#Carreguem els mòduls 

module load modulepath/goolf-1.7.20
module load BLAST+/2.2.30-goolf-1.7.20
module load Exonerate/2.2.0-goolf-1.7.20
module load T-Coffee/11.00.8cbe486-goolf-1.7.20
export PATH=/cursos/20428/BI/bin:$PATH
export PATH=/cursos/20428/BI/soft/genewise/x86_64/bin:$PATH
export WISECONFIGDIR=/cursos/20428/BI/soft/genewise/x86_64/wise2.2.0/wisecfg/

#Escrivim els paths dels diferents fitxers que utilitzarem al llarg del programa, els escrivim en forma de variable per evitar escriure tot el path cada cop que el necessitem

path_db="/cursos/20428/BI/genomes/2017/Taeniopygia_guttata/"
path_BI="/export/home/u114504/"
path_carpeta="${path_BI}/treball/"

#Iniciem el primer loop, ens llegirà les diferents proteïnes
while read proteina; do

path_proteines="${path_carpeta}/${proteina}/${proteina}"

echo "Comencem amb:"${proteina}

#Fem el blast i seleccionem, amb la opció de blast, els e-values significatius, això ho fem posant el filtre de 0.001

tblastn -query ${path_proteines}.fa -db ${path_db}/genome.fa -evalue 0.001 -out ${path_proteines}.blast6 -outfmt 6

echo "Blast fet"

#Extraiem l'e-value, l'scaffold, l'inici i el final en un nou fitxer, per tal de poder analitzar TOTS els hits significatius

cut -f2,9,10,11 ${path_proteines}.blast6 > ${path_proteines}.info.scaffolds

#Treiem el .1 del final dels scaffolds i el MUGN0100 perquè creiem que així és més fàcil de comparar

sed 's/\.1//g' ${path_proteines}.info.scaffolds > ${path_proteines}.nopunt.scaffolds
sed 's/\MUGN0100//g' ${path_proteines}.nopunt.scaffolds > ${path_proteines}.dades.scaffolds

#Definim la variable old scaffold perquè volem analitzar cada scaffold i no cada exó que ha trobat el blast 

oldscaffold="0"

#Definim el contador que ens permetrà emmagatzemar els diferents resultats sense sobreescriure fitxers

contador="1"

#Iniciem el 2n loop, que ens farà tot el procés amb cada scaffold de cada proteïna 

while read hits; do

echo "${contador}"

#Extreiem l'scaffold del fitxer dels hits amb el que realitzarem el 2n loop 

scaffold="$(awk NR==${contador} ${path_proteines}.dades.scaffolds | cut -f1)"

if (( ${oldscaffold} == ${scaffold} )); then 
echo "L'scaffold ja l'hem estudiat"
else echo "Scaffold:"${scaffold}

#Extreiem l'e-value del hit estudiat

evalue="$(awk NR==${contador} ${path_proteines}.dades.scaffolds | cut -f4)" 
echo "e-value:"${evalue}

#Busquem el punt d'inici i de final del hit 

i="$(awk NR==${contador} ${path_proteines}.dades.scaffolds | cut -f2 )" 
f="$(awk NR==${contador} ${path_proteines}.dades.scaffolds | cut -f3 )" 

#Busquem saber si la seqüència on es troben els hits és la cadena reversa o no.

if ((${i}>${f})); then
	inici1="${f}" 
	final1="${i}"
	echo "La seqüència és reversa" 
else 
	inici1="${i}" 
	final1="${f}"
	echo "La seqüència no és reversa"
fi

inici=$((${inici}-50000))
final=$((${final}+50000))
longitud=$((${final}-${inici}))


#Per assegurar-nos que l'inici i el final estan dins de la nostra seqüència mirem que l'inici no sigui negatiu i que el final no surti de la longitud de l'scaffold

if ((${inici}<0));then
	inici2="${inici}"
else
	inici2="1"
fi

#Hit es la variable del nom complet de l'scaffold, que es necessaria per poder fer el Fastafetch 

hit=$((awk NR==${contador} ${path_proteines}.info.scaffolds | cut -f1))

#Longscaffold es la variable que farem servir per trobar la longitud de l'scaffold, que tenim emmagatzemada a un fitxer amb el nom de l-scaffold meny el .1 del final 

longscaffold=$((awk NR==${contador} ${path_proteines}.nopunt.scaffolds | cut -f1))
longitudseq=$((egrep ${longscaffold} ${path_carpeta}/llargada.scaffolds | cut -f1 -d' ')) 

if ((${final}<${longitudseq})); then 
	final2="${final}"
else
	final2="${longitudseq}-${inici2}" 	
fi

#Recalculem la longitud del hit

longitud="${final2}-${inici2}"

echo "inici2:"${inici2}
echo "final2:"${final2}
echo "longitud:"${longitud}

#Fem el fastafetch per extreure l'scaffold que ens interessa del genoma de referència

fastafetch ${path_db}/genome.fa ${path_db}/genome.index ${hit} > ${path_proteines}.${contador}.fastafetch.fa

#Extraiem el nostre hit de l'scaffold, seleccionat al pas anterior, fent el fastasubseq

fastasubseq ${path_proteines}.${contador}.fastafetch.fa ${inici2} ${longitud} > ${path_proteines}.${contador}.subseq.fa

#Exonerate, per obtenir només les seqüències exòniques fem servir la funció egrep

exonerate -m p2g --showtargetgff -q ${path_proteines}.fa -t ${path_proteines}.${contador}.subseq.fa | egrep -w exon > ${path_proteines}.${contador}.exonerate.gff

#Fem el fastaseqfromgff per passar el fitxer de l'exonerate de format .gff a .fasta (el fitxer resultant és de nucleòtids)

fastaseqfromGFF.pl ${path_proteines}.${contador}.subseq.fa ${path_proteines}.${contador}.exonerate.gff > ${path_proteines}.${contador}.nuc.fa

#Per passar de nucleòtids a aminoàcids utilitzem el fastatranslate

fastatranslate -F 1 ${path_proteines}.${contador}.nuc.fa > ${path_proteines}.${contador}.pred.aa.fa

#Per comparar la seqüència de la selenoproteïna coneguda amb la que hem predit del genoma de la nostra espècie fem servir el t-coffee

t_coffee ${path_proteines}.${contador}.pred.aa.fa ${path_proteines}.fa > ${path_proteines}.${contador}.tcoffee

#L'anàlisi de resultats el farem manualment amb els fitxers resultants del t-coffee i dels altres passos si cal 

echo "Ja has acabat amb el scaffold" ${scaffold}

oldscaffold="${scaffold}" #Emmagatzemem  el nom de l'últim scaffold estudiat

contador=$(($contador+1)) #Emmagatzemem el contador i aixi podem saber a quina linia estem 
fi

done < ${path_proteines}.info.scaffolds #Acaba el 2n loop

echo "Ja has acabat amb la proteina" ${proteina} #Acaba el 1r loop

done < llista 
#Llista és on tenim el nom de totes les proteïnes que volem analitzar i on el programa anirà agafant una per una i farà tot el loop agafant el nom com la variable proteina  

