#!/bin/bash

# Script hecho por el grupo 4c de 4º de Biologia Humana del curso 2011-2012 de la Universitat Pompeu Fabra, para la asignatura de Bioinformatica.
## Antes de comenzar, se tiene que crear una carpeta "query" donde tengamos las secuencias que desamos analizar en los genomas.
## Los fallos del script más usuales son equivocarse al introducir los datos o no descargar los programas adjuntos: GFFtoFA.pl y fastaseqfromGFF.pl.

### Escogiendo la secuencia

echo -e "\033[1mBuenos días usuario! Hacemos unos blast? Te acordaste de crear la carpeta query?\033[0m"
echo
echo -e "\033[1mQue query quieres mirar?\033[0m"
ls ./query
echo
echo -en "\033[1mIntroduce el nombre exacto de la query:  \033[0m"
read query
echo

mkdir ./$query

### Escogiendo el organismo

echo -en "\033[1mLos organismos que quieres mirar, de que año son? \033[0m"
read year
echo
ls /cursos/BI/genomes/protists/$year

echo -en "\033[1mIntroduce el nombre exacto del organismo (te aconsejamos que lo copies):\033[0m"
read organismo
echo

mkdir ./$query/$organismo

### BLAST, una vez visualizado la sección de interés del BLAST se puede presionar "Q" para salir del more y continuar con el script.

echo -en "\033[1mHaciendo blast, procesando petición...\033[0m"

export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

mkdir ./$query/$organismo/blast

blastall -p tblastn -i ./query/$query -d /cursos/BI/genomes/protists/$year/$organismo/genome.fa -o ./$query/$organismo/blast/$query.$organismo.blast

echo -e "\033[1mAhora te mostraremos los resultados, con un more. Para salir del more y volver al programa pulsa Q\033[0m"
echo -en "\033[1mAsegurate de quedarte con la posición del hit de interés, si estas seguro dale al enter\033[0m"
read micasa

more ./$query/$organismo/blast/$query.$organismo.blast

### Fastafetch y Fastaindex, para realizar fastasubseq que es necesario para arreglar algunos genomas. 

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

rm ./$organismo.index

fastaindex /cursos/BI/genomes/protists/$year/$organismo/genome.fa ./$organismo.index

echo -en "\033[1mIntroduce donde se ecuentra tu hit:  \033[0m"
read region
echo
echo -e "\033[1mIdentifica la región de tu interés en la lista\033[0m"

grep $region ./$organismo.index > ./grep.$organismo.index
more ./grep.$organismo.index

echo -en "\033[1mInserta el nombre exacto de la region previamente identificada (copia solo la primera columna):  \033[0m"
read hit
echo

mkdir ./$query/$organismo/database

fastafetch /cursos/BI/genomes/protists/$year/$organismo/genome.fa $organismo.index $hit > ./$query/$organismo/database/$organismo.$query.$hit.db

rm grep.$organismo.index
rm $organismo.index

### Fastasubseq

echo -e "\033[1mRealizando fastasubseq ...\033[0m"
echo -en "\033[1mIntroduce el lugar de inicio (normalmente coje unas 10000 pares de bases por debajo de la posición de inicio del hit):  \033[0m"
read inicio
echo
echo -en "\033[1mIntroduce la longitud:  \033[0m"
read longitud
echo

mkdir ./$query/$organismo/fastasubseq

fastasubseq ./$query/$organismo/database/$organismo.$query.$hit.db $inicio $longitud > ./$query/$organismo/fastasubseq/$query.$organismo.genome.fa


echo -en "\033[1m Ahora verás tu secuencia, verifica que has hecho la seleccion que querias con fastasubseq. Recuerda que para salir del more y volver al programa has de pulsar Q. Por favor, dale ahora a enter para proseguir\033[0m"
read mitorre
more  ./$query/$organismo/fastasubseq/$query.$organismo.genome.fa
echo

echo -en "\033[1mQuieres hacer el exonerate, pon SI o NO (recuerda que va en mayusculas)   \033[0m" 
read gafas
if [ "$gafas" = "SI" ]; then


###para el exonerate modificar la U

echo -e "\033[1m*******************************************************************************************************************\033[0m"
echo -e "\033[1mPara realizar el exonerate cambia la U de tu query por una X, y elimina simbolos como # o @ que no se identificaran\033[0m" 
echo -e "\033[1mAprovecha para descargar el programa fastaseqfromGFF.pl,si aun no lo tienes\033[0m"
echo -e "\033[1m*******************************************************************************************************************\033[0m"
echo -e "\033[1mSi no lo haces exonerate no se ejecutara\033[0m"
echo -en "\033[1mRecuerda que es tu query es" $query ". ¿Hecho? Apreta Enter\033[0m"
read micoche

###haciendo el exonerate, si no sale seguir dando al enter, hasta que el programa realice el Genewise.

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

mkdir ./$query/$organismo/exonerate

echo -e "\033[1mExonerate\033[0m"
echo

exonerate -m p2g --exhaustive yes --showtargetgff -q ./query/$query -t ./$query/$organismo/fastasubseq/$query.$organismo.genome.fa > ./$query/$organismo/exonerate/$query.$organismo.exonerate.gff

more  ./$query/$organismo/exonerate/$query.$organismo.exonerate.gff

echo -e "\033[1m*******************************************************************************************************************\033[0m"
echo -e "\033[1mSi te interesa modificar, donde empiezan los exones manualmente para forzar que exonerate te prediga una region mas grande donde quizas se alinee tu selenocisteina, ¡Es ahora o nunca! Recuerda que el archivo lo encontraras en"./$query/$organismo/exonerate/$query.$organismo.exonerate.gff"\033[0m" 
echo -e "\033[1m*******************************************************************************************************************\033[0m"
echo -en "\033[1m ¿Hecho o no te interesa? Pues apreta Enter\033[0m"
read misllaves

egrep -w exon ./$query/$organismo/exonerate/$query.$organismo.exonerate.gff  > $query.$organismo.exo.gff

echo
echo -e "\033[1mObteniendo el cDNA\033[0m"
echo

###haciendo la extraccion del cDNA, con fastaseqfromGFF.pl

mkdir ./$query/$organismo/cDNA

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

chmod u+x ./fastaseqfromGFF.pl

./fastaseqfromGFF.pl ./$query/$organismo/fastasubseq/$query.$organismo.genome.fa $query.$organismo.exo.gff

./fastaseqfromGFF.pl ./$query/$organismo/fastasubseq/$query.$organismo.genome.fa $query.$organismo.exo.gff > ./$query/$organismo/cDNA/$query.$organismo.cDNA

rm $query.$organismo.exo.gff

###haciendo la conversion del cDNA a proteina

echo
echo -e "\033[1mEjecutando fastatranslate\033[0m"
echo

mkdir ./$query/$organismo/proteina

fastatranslate ./$query/$organismo/cDNA/$query.$organismo.cDNA > ./$query/$organismo/proteina/$query.$organismo.fa

echo -e "\033[1m*******************************************************************************************************************\033[0m"
echo -e "\033[1mSe ha creado el archivo con las proteinas en:"./proteina/$query.$organismo.fa" Visita el archivo, y modificalo de forma que te quedes con 1 único patron de lectura. Aprovecha para cambiar el * por una X para visualizar mejor tus resultados\033[0m"
echo -e "\033[1mSi tienes dudas sobre que patron escoger, mira el .gff del exonerate ("./$query/$organismo/exonerate/$query.$organismo.exonerate.gff") ya que te indica con que aminoacidos empieza y acaba tu hit.\033[0m" 
echo -e "\033[1mPor otro lado, aprovecha para revertir la X por la U en tu " $query" al reves que hiciste antes\033[0m"
echo -e "\033[1m*******************************************************************************************************************\033[0m"
echo -en "\033[1m¿Hecho? Para proseguir Apreta Enter\033[0m"
read yes

###haciendo el t_coffee

echo
echo -e "\033[1mHaciendo el t_coffee (si te pone mail, inventatelo en plan: a@a)\033[0m"
echo
mkdir ./$query/$organismo/t_coffee

t_coffee ./query/$query ./$query/$organismo/proteina/$query.$organismo.fa

mv *html ./$query/$organismo/t_coffee 
mv *aln ./$query/$organismo/t_coffee  

echo -e "\033[1mAhora te proponemos que compares tu hit con la base de datos no redundante de ncbi. Esto puede serte muy util, sin embargo tarda un tiempo y si lo deseas puedes saltarte la busqueda y pasar al Genewise.\033[0m"
echo -en "\033[1m Quieres hacer el blast NCBI? Escribe exactamente y en mayusculas: SI o NO   \033[0m"
read respuesta

if [ "$respuesta" = "SI" ]; then

###buscando en NCBI

echo -e "\033[1m Buscando en base de datos NCBI\033[0m"

mkdir ./$query/$organismo/ncbi

export PATH=/cursos/BI/bin/netblast/bin:$PATH
blastcl3 -p blastp -i  ./$query/$organismo/proteina/$query.$organismo.fa -d nr > ./$query/$organismo/ncbi/$query.$organismo.ncbi

fi

fi
###haciendo el genewise

echo -e "\033[1mQuieres hacer el genewise, pon SI o NO (recuerda que va en mayusculas)   \033[0m" 
read pasta
if [ "$pasta" = "SI" ]; then

echo
echo -e "\033[1mHaciendo el Genewise\033[0m"
echo

mkdir ./$query/$organismo/genewise

export PATH=/cursos/BI/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg

echo -e "\033[1mRealizando Genewise ... si tu secuencia es invertida,escribe: -trev. Si no dale al enter y ya esta   \033[0m"
read inverso

genewise -pep -pretty -cdna -gff ./query/$query ./$query/$organismo/fastasubseq/$query.$organismo.genome.fa $inverso > ./$query/$organismo/genewise/$query.$organismo.genewise.gff

echo
echo -e "\033[1mGenerando secuencia peptídica, asegurese de tener el programa GFFtoFA.pl y dar al enter\033[0m"
read si

chmod u+x ./GFFtoFA.pl

mkdir ./$query/$organismo/genewise_prot/

./GFFtoFA.pl < ./$query/$organismo/genewise/$query.$organismo.genewise.gff > ./$query/$organismo/genewise_prot/$query.$organismo.genewise.fa

echo -en "\033[1mA veces nuestro programa no extrae correctamente la secuencia proteica del .gff predicho por Genewise. Por favor, verifica que la secuencia se muestra correctamente, sin lineas raras en las primeras filas hasta el '>'. Encontrarás este archivo en "./$query/$organismo/genewise_prot/$query.$organismo.genewise.fa", modificalo si es necesario y dale al enter para continuar\033[0m"
read miperro

###haciendo el t_coffee del genewise

echo
echo -e "\033[1mHaciendo el t_coffee\033[0m"
echo
mkdir ./$query/$organismo/genewise_t_coffee

t_coffee ./query/$query ./$query/$organismo/genewise_prot/$query.$organismo.genewise.fa

mv *html ./$query/$organismo/genewise_t_coffee
mv *aln ./$query/$organismo/genewise_t_coffee

echo -e "\033[1mAhora te proponemos que compares tu hit con la base de datos no redundante de ncbi. Esto puede serte muy util, sin embargo tarda un tiempo y si lo deseas puedes saltarte la busqueda y finalizar el programa.\033[0m"
echo -en "\033[1m Quieres hacer el blast NCBI? Escribe exactamente y en mayusculas: SI o NO   \033[0m"
read resposta

if [ "$resposta" = "SI" ]; then

###buscando en NCBI_genewise

echo
echo -e "\033[1mBuscando la proteina nueva en el ncbi\033[0m"
echo

mkdir ./$query/$organismo/genewise_ncbi

export PATH=/cursos/BI/bin/netblast/bin:$PATH
blastcl3 -p blastp -i  ./$query/$organismo/genewise_prot/$query.$organismo.genewise.fa -d nr > ./$query/$organismo/genewise_ncbi/$query.$organismo.ncbi

fi
fi
echo
echo -e "\033[1mLa anotacion de tu "$query" en el genoma de "$organismo" ha finalizado com exito\033[0m"
echo




