#!/bin/bash

## Creem totes les carpetes necessaries per a poder executar el programa
mkdir blast_msra
mkdir blastdetallat_msra
mkdir fastafetch_msra
mkdir fastasubseq_msra
mkdir exonerate_msra
mkdir cDNA_msra
mkdir fastatranslate_msra
mkdir fastatranslate_msra_simple
mkdir genewise_msra
mkdir genewisedetallat_msra
mkdir SECIS_msra
mkdir genewise_msra_PEP


### Comencem creant dos bucles per tal de fer totes les combinacions query-genoma

## Assignem a la variable query tot el contingut de la carpeta msra (que es troba dins el directori on ens trobem), carpeta que conte les queries seleccionades de msra

for query in msra/* ; 
do {
	query=`basename $query`
	## Treiem el ".fa" de tots els noms
	query=${query%.fa}
	echo la query que estic fent servir es $query	
	mkdir blast_msra/$query.carpeta
	mkdir blastdetallat_msra/$query.carpeta

## Assignem a la variable genome cadascun dels genomes dels protistes que se'ns han assignat aquest any
	for genome in A.laibachii_Nc14 A.rara C.fasciculata D.discoideum_AX4 D.fasciculatum F.cylindrus G.niphandrodes L.donovani_BPK282A1 L.tarentolae P.capsici P.polycephalum I.multifiliis_strain_G5 S.arctica T.congolense ; 
	do {
		echo el genome que estic fent servir es $genome	

## Per a cada cas (query-genome) realitzem dos blast: el primer ens dona una informacio resumida (només les dades de cada hit) i el segon ens dona la informacio més detallada	
		blastall -p tblastn -i msra/$query.fa -d /cursos/BI/genomes/protists/2012/$genome/genome.fa -o blast_msra/$query.carpeta/$query.$genome.tblastn -m 9 -F F
		blastall -p tblastn -i msra/$query.fa -d /cursos/BI/genomes/protists/2012/$genome/genome.fa -o blastdetallat_msra/$query.carpeta/$query.$genome.tblastn -F F


## Creem un document evalues.fa on recollim en tres files: el nom del genome, el nom de la query i a la tercera fila tenim tres columnes: l'identificador del genoma, el punt d'inici (start) i el e-value, de cada hit
		echo $query >> evalues_msra.fa
		echo $genome >> evalues_msra.fa
## Escollim els evalues que continguin "e-" per seleccionar els hits significatius
		egrep -v '#' blast_msra/$query.carpeta/$query.$genome.tblastn | egrep 'e-' | awk '{print $2 "\t" $9"\t" $10"\t" $11}' >> evalues_msra.fa
		echo  >> evalues_msra.fa
		echo  >> evalues_msra.fa


## Creem un segon document temporal (blast_rotatiu.fa) on emmagatzemem cada vegada el resultat que surt del blast (unicament seleccionem el corresponent a la tercera fila del document evalues.fa) 
		egrep -v '#' blast_msra/$query.carpeta/$query.$genome.tblastn | egrep 'e-' | awk '{print $2 "\t" $9"\t" $10}' > blast_rotatiu_msra.fa
## Creem una carpeta temporal (hits) que més endavant borrarem. Gracies a un programa perl (linies.pl), posem cada linia que entra al document blast_rotatiu_msra.fa a un document (hit$contador) diferent dins la carpeta hits. Amb el if seleccionem unicament els arxius blast_rotatiu_msra.fa que contenen hits. 
		mkdir hits
		x=`egrep -v '#' blast_msra/$query.carpeta/$query.$genome.tblastn | egrep 'e-' | awk '{print $2 "\t" $9"\t" $10}'`
		if [ "$x" != "" ]; then
			echo La query que estic utilitzant es $query i el genoma es $genome
			./linies.pl blast_rotatiu_msra.fa
			for hit in hits/* ;
			do {
				cat $hit
			} done
		
			for hit in hits/* ;
			do {
				hit=`basename $hit`
				echo el nom del hit es $hit
		
## Creem les variables contig, start i end de cada hit, on cada una conté una de les columnes corresponents de cada document hit$contador (identificador, punt d'inici i punt de final)
				contig=`cut -f 1 $hit`
				start=`cut -f 2 $hit`
				end=`cut -f 3 $hit`
				echo $contig
				echo $start
				echo $end

				if [ $start -gt $end ]; then
					start=$end
				fi
				start=$(( $start - 10000 ))
				echo $start
				length="30000"

## Amb el fastafech seleccionem el contig que volem
				mkdir fastafetch_msra/$hit.$query.$genome
				fastafetch /cursos/BI/genomes/protists/2012/$genome/genome.fa index/$genome.index "$contig" > fastafetch_msra/$hit.$query.$genome/$hit

## Fem el fastasubseq per a retallar la regio del voltant del gen d'interés, derivant els FATAL ERROR a un document error_mrsa.txt
				mkdir fastasubseq_msra/$hit.$query.$genome
				fastasubseq fastafetch_msra/$hit.$query.$genome/$hit $start $length > fastasubseq_msra/$hit.$query.$genome/$hit 2> error_msra.txt

## En els casos en que els llocs d'inici o final definits es surten del contig corresponent, es tornen a definir a partir del document error_msra.txt tot depenent del tipus d'error:	
				error=`cat error_msra.txt`
				if [ "$error" != "" ]; then
					error_final=`cat error_msra.txt | egrep "end before end"`
					error_inici=`cat error_msra.txt | egrep "Unknown flag"`
						
				## Si el punt final es troba fora del contig, establim com a final l'ultima posicio del contig
					if [ "$error_final" != "" ]; then
						llargada=$error_final						
						llargada=${llargada##*(}
						llargada=${llargada%)}
						length=$(( $llargada - $start))
						fastasubseq fastafetch_msra/$hit.$query.$genome/$hit $start $length > fastasubseq_msra/$hit.$query.$genome/$hit
					fi
						
				## Si el punt d'inici es troba fora del contig, establim com a inici la primera posicio del contig
					if [ "$error_inici" != "" ]; then
						start="0"
						fastasubseq fastafetch_msra/$hit.$query.$genome/$hit $start $length > fastasubseq_msra/$hit.$query.$genome/$hit 2> error_msra_2.temp
				## En el cas que es trobin els dos errors, cal tornar a derivar l'error a un altre document
						error_2=`cat error_msra_2.temp`
							if [ "$error_2" != "" ]; then
								error_2_llarg=`cat error_msg_2.temp | egrep "before end"`
								if [ "$error_2_llarg" != "" ]; then
									llargada=$error_2_llarg						
									llargada=${llargada##*(}
									llargada=${llargada%)}
									length=$(( $llargada - $start))
									fastasubseq fastafetch_msra/$hit.$query.$genome/$hit $start $length > fastasubseq_msra/$hit.$query.$genome/$hit 
								fi
							fi 
					fi
		
				fi

## Realitzem l'exonerate i emmagatzemem les dades a un document (resultats_exonerate1.fa) per a saber quins hits tenen resultats i quins no 
				mkdir exonerate_msra/$hit.$query.$genome
				exonerate -m p2g --showtargetgff -q msra/$query.fa -t fastasubseq_msra/$hit.$query.$genome/$hit > exonerate_msra/$hit.$query.$genome/$hit
				echo $query >> resultats_exonerate1.fa
				echo $genome >> resultats_exonerate1.fa
				egrep -w exon exonerate_msra/$hit.$query.$genome/$hit >> resultats_exonerate1.fa
				echo  >> resultats_exonerate1.fa

## Generem un document rotatiu amb el resultat de l'exonerate que s'obté cada vegada per cada hit
				egrep -w exon exonerate_msra/$hit.$query.$genome/$hit > exonerate_rotatiu1.fa
## Generem la sequencia FASTA només dels hits que tenen algun exo
				y=`cat exonerate_rotatiu1.fa`
				if [ "$y" != "" ]; then
## Obtenim el cDNA i amb el fastatranslate traduïm a les 6 pautes de lectura la possible proteïna
				    fastaseqfromGFF.pl fastasubseq_msra/$hit.$query.$genome/$hit exonerate_rotatiu1.fa > cDNA_msra/$hit.$query.$genome
				    ## Amb la següent comanda realitzem un altre fastatranslate però que ens dóna directament la que considera millor pauta de lectura	
				fastatranslate cDNA_msra/$hit.$query.$genome -F 1 > fastatranslate_msra_simple/$hit.$query.$genome
				    fastatranslate cDNA_msra/$hit.$query.$genome > fastatranslate_msra/$hit.$query.$genome
				fi


## Analitzem amb genewise. Primer fem un analisi detallat: 
				mkdir genewisedetallat_msra/$hit.$query.$genome
				genewise -cdna -pretty -both -gff msra/$query.fa fastasubseq_msra/$hit.$query.$genome/$hit > genewisedetallat_msra/$hit.$query.$genome/$hit
## El segon ens dona unicament la proteïna resultant (-pep) tant el forward com en reverse (-both), i és el que fem servir per seguir amb el programa
				mkdir genewise_msra/$hit.$query.$genome
				genewise -pep -both msra/$query.fa fastasubseq_msra/$hit.$query.$genome/$hit > genewise_msra/$hit.$query.$genome/$hit


## Seleccionem la seqüencia proteica més llarga del genewise mitjançant un programa perl (pep.pl)
				mkdir genewise_msra_PEP/$hit.$query.$genome
				prot=`cat genewise_msra/$hit.$query.$genome/$hit`
				if [ "$prot" != "" ]; then
					./pep.pl genewise_msra/$hit.$query.$genome/$hit > genewise_msra_PEP/$hit.$query.$genome/$hit
					echo $hit $query $genome fet
				fi


## Busquem els elements SECIS a partir de la subseq utilitzada per predir la proteïna
				mkdir SECIS_msra/$hit.$query.$genome
				./SECISearch.pl < fastasubseq_msra/$hit.$query.$genome/$hit > SECIS_msra/$hit.$query.$genome/$hit

			} done
		
		fi
        

## Eliminem la carpeta temporal hits
		rm -r hits

	} done

} done

## Fem el tCoffee entre les queries i les proteïnes predites			
mkdir tcoffee_msra_exonerate
mkdir tcoffee_msra_genewise

for query in msra/* ; 
do {
	query=`basename $query`
	query=${query%.fa}

	for arxiu in fastatranslate_msra_simple/* ; do {
		arxiu=`basename $arxiu`
		echo Fem el tcoffee_exonerate de $arxiu
		mkdir tcoffee_msra_exonerate/$arxiu
		t_coffee msra/$query.fa fastatranslate_msra_simple/$arxiu -outfile=tcoffee_msra_exonerate/$arxiu/$query
			
	} done 
} done

for carpeta in genewise_msra_PEP/* ; do {
	carpeta=`basename $carpeta`

	for hit in genewise_msra_PEP/$carpeta/* ; do {
		hit=`basename $hit`
		
		for query in msra/* ; do {
			query=`basename $query`
			query=${query%.fa}
			
			mkdir tcoffee_msra_genewise/$carpeta
			echo la query que estic fent servir es $query

			t_coffee msra/$query.fa genewise_msra_PEP/$carpeta/$hit -outfile=tcoffee_msra_genewise/$carpeta/tcoffee_$query
		} done
	} done
} done