#PEDRETTI_SCRIPT (hasta tcoffee)
#Programa creado por grupo "Ey" (By David Pedreño, Sara Vicente, Alberto Bofill, Laura Vidaurre y Jorge Sanchez). 16mar-2011, Bioinformatica, UPF.

#Requisitos programa: tener creada carpeta "Semillas" con las querys que queramos usar para correr el programa. Tener en la misma carpeta los programas perl Evalues.pl y fastaseqfromGFF.pl Hemos diseñado el programa de manera que procesará todo el análisis hasta tcoffee en las querys de las selenoproteinas que añadamos a la carpeta Semillas. 




#!/bin/bash

export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH             # pel NCBI Blast
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/                     # pel NCBI Blast
export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH	   # per l'exonerate
export PATH=/cursos/BI/bin:$PATH			   # pel genewise, fastaseqfromGFF.pl i t_coffee
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg    # pel genewise	




#CREACION ARCHIVOS FASTAINDEX
#Obtendremos un fasta index para cada organismo, separados entre ellos por carpetas con sus nombres de especie, que se crearan gracias al programa)

directori="/cursos/BI/genomes/protists/2011"      																					# Asignamos la ruta hasta nuestros genomas a nuestra nueva variable declarada $directori. Para cada año que corramos, deberemos cambiar el año final de carpeta a elegir.

mkdir Fastaindex        		# Creamos carpeta Fastaindex   

#Declaramos la variable genome, donde contendrá cada uno de los nombres de las especies a estudiar. Esta variable nos servira para encontrar las carpetas de cada especie en la ruta designada por $directori. De esta manera, podremos llegar al archivo genome.fa, i redirigirlo finalmente a carpetas con los nombres de cada especie en Fastaindex.

for genome in A.taiwanensis B.hominis C.merolae C.muris C.owczarzaki C.parvum E.siliculosus H.arabidopsidis N.gruberi P.pallidum P.tricornutum P.ultimum S.parasitica T.trahens ; do {

 mkdir Fastaindex/$genome 		
 fastaindex  $directori/$genome/genome.fa Fastaindex/$genome/genome.index  
} done







#CREACION NOMBRES DIRECTORIOS PARA TBLASTN. 
#Usaremos los nombres de los archivos ".fa" de cada query dentro de la carpeta Semillas para declarar la variable semilla_fa, la cual se obtendrá por el uso del comando basename, que nos proporciona como variable todo aquello que esté posterior al ultimo simbolo "/". De tal manera que nos copiara el nombre entero de "archivo.fa". Posteriormente, obtendremos la variable semilla, en la cual hemos hecho desaparecer la extensión terminal del archivo ".fa", y las hemos acumulado en la variable semillas, conteniendo esta ultima todas las querys procedentes de la variable semilla separadas entre ellas por un espacio.

directori="/cursos/BI/genomes/protists/2011"

for carpetes in ./Semillas/*; do {
 semilla_fa=`basename $carpetes`
 semilla=${semilla_fa%.fa}
 semillas=$semillas" $semilla"
} done


mkdir Tblastn 
mkdir Evalues


# OBTENCION TBLASTN

# Creamos variable semilla con listas con las rutas de las semillas, que es la variable anteriormente creada ($semillas)
# Creamos variable genome para numerar genomas, como anteriormente hemos hecho para Fastaindex
# Para el analisis de los hits creados por Tblastn usaremos la variable query para agrupar los hits por query, ordenados en el archivo de texto blastresults.txt, creado para poder analizar mas rapidamente el listado de hits por querys y especies.
# Corremos Tblastn con filtro de -e 0.001, para evitar alineamientos que sean debidos por azar. Posteriormente creamos variable $filtreevalue la cual nos permitirá preparar los archivos para correr el programa Evalues.pl, el cual renombrará los archivos temporales segun el score de su Evalue.
# El programa Evalues.pl nos permite dosificar los evalues archivos individuales, numerados por Evalues (1,2,3...) , por el genoma del cual procede, y por la query de la cual procede también. ARGV 0 i 1 son las querys y genomas, respectivamente, y hacemos que los datos almacenados temporalmente en el archivo temporal (el Evalue aceptado) se introduzca en el archivo. Todos estos archivos se estructuraran en carpetas asignadas por querys i genomas, dentro de los cuales numerados por Evalues.

rm blastresults.txt

for semilla in $semillas ; do {
 for genome in A.taiwanensis B.hominis C.merolae C.muris C.owczarzaki C.parvum E.siliculosus H.arabidopsidis N.gruberi P.pallidum P.tricornutum P.ultimum S.parasitica T.trahens ; do {
  
 queryinicial=`cat Semillas/$semilla.fa | egrep ^">"`
 queryprocessada=${queryinicial#*>}
 echo $queryprocessada > temporal
 query=`cut -f 1 -d " " temporal`
 rm temporal

 blastall -p tblastn -i ./Semillas/$semilla.fa -d $directori/$genome/genome.fa -o  Tblastn/$semilla.$genome.txt -F F -m9 -e 0.001

       	  
 echo ">${semilla}_$genome:" >> blastresults.txt
 egrep $query Tblastn/$semilla.$genome.txt | egrep e- >> blastresults.txt
                       
 filtreevalue=`egrep $query Tblastn/$semilla.$genome.txt | egrep e-`
 if [ "$filtreevalue" != "" ]; then
  egrep $query Tblastn/$semilla.$genome.txt | egrep e- > Evalueacceptat.temporal
  mkdir Evalues/$semilla 
  mkdir Evalues/$semilla/$genome
  ./Evalues.pl $semilla $genome < Evalueacceptat.temporal
 fi
 echo ${semilla}_$genome "fet"

 } done
} done

rm ./*.temporal



#PROGRAMA EXONERATE_CDNA_TCOFFEE


rm resultats_exonerate.txt


mkdir Fastafetch
mkdir Fastasubseq
mkdir Exonerate
mkdir Exonerate/Seleccionados
mkdir cDNA
mkdir proteina
mkdir proteina/completa
mkdir tcoffee	

directori="/cursos/BI/genomes/protists/2011"

# Creamos asignacion que nos permitirá correr las comandas de fastafetch, fastasubseq, exonerate, fastaseqfromGFF, fastatranslate i tcoffee para cada query, genoma, y Evalue que tengamos dentro de la carpeta Evalues creada anteriormente, siempre yendo recursiva y jerarquicamente descendiendo (de la query Dsba_ciona_intestinalis, se realizaran todos los genomas, i dentro de cada genoma, todos sus evalues, para pasar posteriormente a la siguiente query (Dsba_psicromonas_ingrahamii), repitiendo el mismo proceso recursivo.

for semillas in ./Evalues/* ; do {
 for genome in $semillas/* ; do {
  for Evalue in $genome/* ; do {

# Assignamos a las variables los nombres sin rutas ni extensiones (sin "/" ni ".txt") mediante las opciones basename i %.txt (solo para la variable Evalue, puesto que contienen extension .txt), respectivamente.
               
   semillas=`basename $semillas`
   genome=`basename $genome`
   Evalue=`basename $Evalue`
   Evalue=${Evalue%.txt}
   echo ${semillas}_${genome}_${Evalue}:

# Detección query (aparece a primera linea del arxivo Evalue)

   queryinicial=`cat Evalues/$semillas/$genome/$Evalue.txt | egrep ^">"`
   queryprocessada=${queryinicial#*>}
   echo $queryprocessada > temporal
   query=`cut -f 1 -d " " temporal`
   rm temporal
 

# FASTAFETCH
# Creamos variable $identificador con la segunda columna (-f 2) de cada archivo Evalue.txt de las carpetas, conteniendo justamente el identificador de la region genomica.
   mkdir Fastafetch/$semillas
   mkdir Fastafetch/$semillas/$genome
   identificador=`cat $query Evalues/$semillas/$genome/$Evalue.txt | cut -f 2`
   fastafetch $directori/$genome/genome.fa Fastaindex/$genome/genome.index $identificador > Fastafetch/$semillas/$genome/$Evalue.fa

# FASTASUBSEQ
# Definiremos el inicio i final del hit en la region identificadora mediante variables inici i final respectivamente. En el caso de secuencias con sentido inverso, solucionaremos el problema haciendo que la variable inicio empieze con el valor de la variable final. Posteriormente, restaremos 20000 bases a la posicion inicial de partida y realizaremos un fastasubseq con una longitud de 50000 declarada en la variable $lenght.

   inici=`cat $query Evalues/$semillas/$genome/$Evalue.txt | cut -f 9`
   final=`cat $query Evalues/$semillas/$genome/$Evalue.txt | cut -f 10`

   if [ "$inici" -gt "$final" ]; then
     inici="$final"
   fi             
   start=$(( $inici - 20000 ))
   length="50000"
   fastasubseq Fastafetch/$semillas/$genome/$Evalue.fa $start $length > Fastasubseq/${semillas}_${genome}_${Evalue}.fastasubseq.fa 


# EXONERATE
# Analizaremos el exonerate buscando regiones exonicas, creando un filtro con la variable $cds (encargada de solo escoger los exonerates aptos y con regiones codificantes). De tal manera que en la carpeta Seleccionados obtendremos la seleccion de exonerates aptos para el procesamiento, donde solo agruparemos los exones y sus regiones.

   echo "Iniciant exonerate de" ${semillas}_${genome}_${Evalue}

   exonerate -m p2g --showtargetgff --exhaustive yes -q Semillas/$semillas.fa -t Fastasubseq/${semillas}_${genome}_${Evalue}.fastasubseq.fa  > Exonerate/EXONERATE_completo_${semillas}_${genome}_${Evalue}.gff 
		
   cds=`grep cds Exonerate/EXONERATE_completo_${semillas}_${genome}_${Evalue}.gff`
   if [ "$cds" != "" ]; then
     echo ${semillas}_${genome}_${Evalue} : >> resultats_exonerate.txt
     exonerate -m p2g --showtargetgff --exhaustive yes -q Semillas/$semillas.fa -t Fastasubseq/${semillas}_${genome}_${Evalue}.fastasubseq.fa | egrep -w exon > Exonerate/Seleccionados/EXONERATE_${semillas}_${genome}_${Evalue}.gff 
   fi
   echo "ACABAT EXONERATE de" ${semillas}_${genome}_${Evalue} 
 
# cDNA

# Usaremos un filtro para seleccionar aquellos archivos exonerate que contengan regiones exonicas, de manera que eliminaremos posibles malos candidatos para la extraccion de cDNA para traducirlo a proteinas mediante los comandos fastaseqfromGFF (mediante el uso del programa fastaseqfromGFF.pl) y fastatranslate. De este ultimo comando realizamos dos variantes (la normal i una con comando -F 1, capaç de darnos el mejor transcrito) Para comprobar de que la opcion -F 1 sea fiable, realizamos ambas acciones y comparamos en un numero significativo de archivos si la proteina que escogía -F 1 era la mejor siempre, y realmente lo era, asi que los usamos en el posterior procesamiento de datos para tcoffee como mejor translate.                               
               
   echo "Haciendo t_coffee:"
                               
   exon=`grep -w exon Exonerate/Seleccionados/EXONERATE_${semillas}_${genome}_${Evalue}.gff`
   grep -w exon Exonerate/Seleccionados/EXONERATE_${semillas}_${genome}_${Evalue}.gff > goodones.temporal
   if [ "$exon" != "" ]; then
                       

       fastaseqfromGFF.pl Fastasubseq/${semillas}_${genome}_${Evalue}.fastasubseq.fa goodones.temporal > cDNA/cDNA_${semillas}_${genome}_${Evalue}.fa
       fastatranslate -F 1 cDNA/cDNA_${semillas}_${genome}_${Evalue}.fa > proteina/proteina_${semillas}_${genome}_${Evalue}.fa
       fastatranslate  cDNA/cDNA_${semillas}_${genome}_${Evalue}.fa > proteina/completa/proteina_completa_${semillas}_${genome}_${Evalue}.fa

                

# T_COFFEE              		
# Creamos carpeta tcoffe i los almacenamos por carpetas organizadas por querys y genomas.
			
       mkdir tcoffee/$semillas
       mkdir tcoffee/$semillas/$genome
       t_coffee Semillas/$semillas.fa proteina/proteina_${semillas}_${genome}_${Evalue}.fa > tcoffee/$semillas/$genome/tcoffee_${semillas}_${genome}_${Evalue}.txt
              	                                    
		
   fi

  } done 
 } done 
} done 

rm ./*.temporal

