#!/bin/bash

# Posada en funcionament el 11/11/2013

# Genoma de la Foca de Weddell: "/cursos/BI/genomes/project_2014/Leptonychotes_weddellii/genome.fa"

./cambiarU.pl  #Posem en marxa el programa que cambiarà les "U" per "X". 


# Comencem afegint els permisos amb els que treballarem
export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH            # pel NCBI Blast
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/                    # pel NCBI Blast
export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH      # per l'exonerate


for query in ./*.fa; # La variable query contindrà només el nom de la selenoproteina amb la que estem corrent el programa. 
do {
   query=`basename $query`
   query=${query%.fa}     # Eliminem el ".fa" del nom de l'arxiu de la selenoproteina
   echo Selenoproteina es $query .
   mkdir $query.resultat  # Creem una carpeta amb el seu nom.resultat on enmagatzemem a els resultats obtinguts.

   # Trobem els hits, és a dir a on la selenoproteina pot trobar-se en el genoma amb una certa probabilitat
   blastall -p tblastn -m9 -i ./$query.fa -d /cursos/BI/genomes/project_2014/Leptonychotes_weddellii/genome.fa -o $query.resultat/blast

   # Seleccionem el hit, el més elevat, a on la selenoproteina es pot trobar amb major probabilitat
   # Fins ara sempre el 1r scaffold de l'arxiu blast és el de major probabilitat però ho verifiquem.
   
   # El 1r egrep ignora les línies que comencen amb #
   # El 2n egrep selecciona la 1a línia amb el text 'e-' ja que els evalues que continguin 'e-' són els hits significatius.
   # L'awk selecciona el 2n valor de la taula que és el SubjectID
   # Guardem en un nou arxiu temporal de nom blast.evalues només els SubjectID i els evalue
   egrep -v '#' $query.resultat/blast | egrep 'e-' | awk '{print $2 "\t" $11}' > $query.resultat/blast.evalues

   # buscarem el mínim evalue de l'arxiu utilitzant un programa en Perl
   ./evalue.pl $query.resultat/blast.evalues
   
   # La variable SubjectID contindrà l'identificador de l'scaffold de l'arxiu blast   
   SubjectID=`more evaluebuscat.tmp`
       
# Obtenim l'scaffold trobat en el pas anterior del genoma (DNA)
   fastafetch /cursos/BI/genomes/project_2014/Leptonychotes_weddellii/genome.fa ~/luke/genome.index "$SubjectID" > $query.resultat/genomic.fa

# Tornem a realitzar els permisos per si de cas.   
export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH

   # Realitzem l'exonerate
   exonerate -m p2g --showtargetgff -q ./$query.fa -t $query.resultat/genomic.fa > $query.resultat/exonerate.gff
   
   # Extreiem els exons de l'exonerate anterior
   egrep -w exon $query.resultat/exonerate.gff  > $query.resultat/exonerateexons.gff

# Novament, tornem a executar els permisos.    
export PATH=/cursos/BI/bin:$PATH    
    
 # Passem de fitxer .gff a .fa (DNA)
   fastaseqfromGFF.pl $query.resultat/genomic.fa $query.resultat/exonerateexons.gff > $query.resultat/fastaseqfromGFF.fa

   # Predim la seqüència proteica de la regió trobada.

   fastatranslate $query.resultat/fastaseqfromGFF.fa -F 1 > $query.resultat/fastatranslateF1.fa
   # Comparem la proteina predida amb la nostra proteina.

   t_coffee $query.resultat/fastatranslateF1.fa ./$query.fa > $query.resultat/tcoffee.txt
   
} done

#Al final obtenim una carpeta amb el nom de la proteina que hem mirat en cada cas a on hi trobarem tots els arxius d'interès.