#!/bin/bash

export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH

cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH

export PATH=/cursos/BI/bin:$PATH

export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg

# Todos los archivos obtenidos en el fastasubseq los guardamos en la carpeta 'subseq'.

mkdir subseq

for fi in `grep 2012 /cursos/BI/genomes/protists/genomes_list_and_info.tab | cut -f 4`; do {

genome=`grep $fi /cursos/BI/genomes/protists/genomes_list_and_info.tab | cut -f 2`;

<

for seq in selprotT_human.fa selprotT_mouse.fa selR1_human.fa selR2_human.fa selR3_human.fa selR_cerevisiae.fa selR_droso.fa selR_elegans.fa; do {

# A continuación, definimos el principio y el final del hit de interés. De los archivos de la carpeta 'contigs' extraemos las columnas 9 y 10 que corresponden al principio y al final de la región respectivamente.

principio=`cat contigs/${seq}.CONTRA.${genome}.contig.fa | cut -f 9`

final=`cat contigs/${seq}.CONTRA.${genome}.contig.fa | cut -f 10`

# Para llevar a cabo el fastasubseq nos interesa el 'start' y el 'lenght' de la región de interés. En primer lugar, miramos si el 'principio' es mayor que el 'final' y si esto es así, definimos el principo com el final. En segundo lugar, para definir el 'start' y el 'lenght', extendemos la región de interés 1000 nucleótidos 'upstream' y extendemos 20000 nucleótidos desde el 'start'.

if [ "$principio" -gt "$final" ]; then

$principio="$final"

fi

start=$(( $principio - 10000 ))

length=$(( $start + 20000 ))

echo "$start" "$length"

# A continuación ejecutamos el fastasubseq para obtener la región genómica donde creemos que se encuentra el gen de interés.

# El siguiente programa nos sirve para solucionar los errores que puede haber al ejecutar el programa.

fastasubseq fetch/fastafetch.${seq}.CONTRA.${genome}.fa $start $length > subseq/subseq.${seq}.CONTRA.${genome}.fa 2> error_msg.temp

error=`cat error_msg.temp`

if [ "$error" != "" ]; then

error_largo=`cat error_msg.temp | egrep "before end"`

error_corto=`cat error_msg.temp | egrep "Unknown flag"`

# Si se pasa de largo, corregimos la longitud para que llegue al final.

if [ "$error_largo" != "" ]; then

longitud=$error_largo

longitud=${longitud##*(}

longitud=${longitud%)}

length=$(( $longitud - $start))

fastasubseq fetch/fastafetch.${seq}.CONTRA.${genome}.fa $start $length > subseq/subseq.${seq}.CONTRA.${genome}.fa

fi

# Si el principio está demasiado cerca, establecemos el 'start' en 0.

if [ "$error_corto" != "" ]; then

start="0"

fastasubseq fetch/fastafetch.${seq}.CONTRA.${genome}.fa $start $length > subseq/subseq.${seq}.CONTRA.${genome}.fa 2> error_msg_2.temp

error_2=`cat error_msg_2.temp`

if [ "$error_2" != "" ]; then

error_2_largo=`cat error_msg_2.temp | egrep "before end"`

if [ "$error_2_largo" != "" ]; then

longitud=$error_2_largo

longitud=${longitud##*(}

longitud=${longitud%)}

length=$(( $longitud - $start))

fastasubseq fetch/fastafetch.${seq}.CONTRA.${genome}.fa $start $length > subseq/subseq.${seq}.CONTRA.${genome}.fa 2> error_msg_3.temp

fi

fi

fi

fi

} done

} done