Selenoproteïnes

Automatització

A continuació s'exposen els diferents scripts creats, dins els quals s'ha explicat, pas per pas, quins programes es van executant i quines carpetes i fitxers de sortida es van creant.

Primer de tot, però, cal donar els següents permisos per tal que es puguin executar els programes NCBI Blast, exonerate, GeneWise, fastaseqfromGFF.pl i el TCoffee des del nostre ordinador local:

$ 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, el fastaseqfromGFF.pl i el t_coffee
$ export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg	# pel GeneWise

Amb el script següent hem automatitzat tots els passos per a realitzar tot el procés per a la família MsrA, cal dir però, que hem fet servir el mateix script per a SelW canviant simplement els noms dels fitxers.

automatitzacio_msra.sh


#!/bin/bash

## Creem totes les carpetes necessàries 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 conté 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 és $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 dóna una informació resumida (només les dades de cada hit) i el segon ens dóna la informació 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 (únicament 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. Gràcies 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 únicament 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 regió 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'última posició 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 posició 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 exó
				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 anàlisi 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 dóna únicament 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üència 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

Tornar a dalt »

Com hem comentat al script anterior, s'han creat dos programes Perl per tal de realitzar les accions repetitives comentades anteriorment. Els programes són els següents:

linies.pl

   #!/usr/bin/perl -w

   use strict;
   my $contador;

   open F0, $ARGV[0];
   $contador=1;

   while (<FO>)
   {
   	open FICHERO, ">hits/hit$contador";
   	print FICHERO "$_\n";
   	$contador++;
   }
   close FICHERO;
   close F0;
   $contador=$contador-1;

   print "hi ha $contador hits\n";

Tornar a dalt »

pep.pl

   #!/usr/bin/perl -w

   use strict;

   my @v;
   my $inici=0;
   my $long=0;
   my $i=0;
   my $cont1=0;
   my $cont2=0;
   my $seq1="";
   my $seq2="";
   my $l;

   open F0, $ARGV[0];

   while (<FO>){$l.=$_;}
   close F0;

   @v=split(//,$l);
   $long=scalar(@v);
   $i=0;

   while ($i < $long)
   {
   	if ($v[$i] eq ">")
   	{
   		$inici++;
   	}
		
   	if ($inici==2)
   	{
   		$seq2.=$v[$i];
   		$cont2++;		
   	}
   	else
   	{
   		$seq1.=$v[$i];
   		$cont1++;
   	}
   	$i++;
		
   }

   if ($cont1 > $cont2 || $cont1 == $cont2)
   {
   	print $seq1;
   }
   else
   {
   	print $seq2;
   }

La cerca dels elements SECIS la hem realitzat manualment a partir de la pàgina web SECISearch , com ja ha estat comentat anteriorment (tornar).

Tornar a dalt »

Automatització

Recerca de selenoproteïnes en protistes:
SelW,MrsA

Carlevaro J.; Casas A.; Eraso A.; Farrés A.; Gamallo B.

Recerca de selenoproteïnes en protistes: SelW,MrsA

Carlevaro J.; Casas A.; Eraso A.; Farrés A.; Gamallo B.

Automatització

Recerca de selenoproteïnes en protistes:
SelW,MrsA