Selenoproteïnes En genomes de protistes

La base de dades SelenoDB ens proporciona les seqüències proteiques de GPx i SelN de diferents organismes.
Tot i que el criteri teòric sigui usar com a queries les seqüències proteiques filogenèticament més properes als protistes investigats, en el nostre treball hem usat seqüències d'H.sapiens per ser l'organisme en què s'han descrit més selenoproteïnes en la família Glutathione Peroxidase. En el cas de la SelN hem utilitzat, a més de les seqüències que ens proporciona SelenoDB d'humà i ratolí, dues seqüències més extretes de la base de dades de l'NCBI, una de Ciona intestinalis i una altra de Dario rerio.

Tornar

Els genomes a investigar van ser facilitats pels professors de l'assignatura de bioinformàtica i pertanyen a les següents espècies: A.taiwanensis, C.muris, E.siliculosus, N.gruberi, P.ultimum, C.owczarzaki, P.pallidum, S.parasitica, C.merolae, C.parvum, H.arabidopsidis, P.tricornutum i T.trahens.

Tots aquest genomes es poden consultar indexats a la carpeta Genomes.

Tornar

El BLAST (Basic Local Alignment Search Tool) és una eina bioinformàtica que ens permet trobar similituds locals entre seqüències usant un principi heurístic i programació dinàmica. Dins de les diferents modalitats d'aquest recurs, per a realitzar aquest treball hem utilitzat el tBLASTN, el qual compara seqüències proteiques com a queries contra bases de dades nucleotídiques. Per fer-ho, tradueix a proteïna les seqüències genòmiques de la base de dades en els 6 possibles marcs de lectura (ORFs) i alinea els resultats amb la nostra query. Finalment, si algun dels alineaments compleix els valors estadístics establerts, els reporta com a hits.
En el nostre treball usem com a bases de dades genomes seqüenciats recentment i com a queires les seqüències aminoacídiques de les proteïnes assignades. Aquestes seqüències estan emmagatzemades en format FASTA, un format que no podem usar per fer el BLAST, de manera que el primer que farem serà formatejar la base de dades amb la següent comanda:

Instrucció	Detalls
$ formatdb -i /cursos/BI/genomes/protists/any/nom_protista/genome.fa -p F -n genoma_protista.fa	on, l'argument -i indica la ruta d'arxius per accedir al genoma, el paràmetre -p F (false) ens informa que la base de dades no és un arxiu de proteïna i on l'argument -n ens permet donar nom a l'arxiu de sortida de la base de dades.

Ara, per poder executar el BLAST cal posar les següents ordres al shell:

Instrucció

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

$ cp /cursos/BI/bin/ncbiblast/.ncbirc ~/ da:

A continuació, ja podem executar el BLAST amb la següent ordre:

Instrucció	Detalls
$ blastall -p tblastn -i query.fa -d genoma_protista.fa -o tblastn_queryVSgenoma.fa	on l'argument -p especifica la modalitat de BLAST a usar, el paràmetre -i especifica la query, el paràmetre -d indica la base de dades i -o especifíca el nom de l'arxiu de sortida del BLAST

Tornar

Un cop executat el tBLASTn i reportats els hits, cal comprovar que siguin significatius estadística i biològicament. Això es fa en funció de dos criteris:
- E-value menor a 10^-4 ja que aquest paràmetre indicia el nombre de hits igual o més bons que l'obtingut que esperem trobar per atzar en la nostra base de dades.
- Cal assegurar-nos que la selenocisteïna de la query alineï amb una cisteïna o una selenocisteïna del genoma analitzat.

Això pot fer-se manualment obrint els arxius de resultats i revisant cada hit o pot automatitzar-se en diferents graus. En el nostre cas hem creat un programa que analitza els E-values i arxiva per separat cada hit significatiu. Posteriorment hem revisat manualment els alineaments de cada hit.

Tornar

Els resultats significatius del tBLASTn, guardats en arxius independents, contenen informació sobre el punt d'inici i l'extensió de cada regió alineada. Aquestes dades són les que usem per extreure fragments d'interès dels genomes estudiats (que potencialment contenen la proteïna buscada) i guardar-los indexats per facilitar-ne l'anàlisi posterior.
Per fer això, primer indexem els genomes amb fastaindex usant la comanda, de manera que la comanda d'extracció pugui treballar amb menys exigències computacionals:

Instrucció

$ fastaindex /cursos/BI/ genomes/protists/any/nom_protista/genome.fa genoma_protista.index

A continuació, delimitem la regió on se situen els hits amb fastafetch:

Instrucció

$ fastafetch /cursos/BI/ genomes/protists/any/nom_protista/genome.fa genoma_protista.index nomseq > nomseq.fa

Finalment, amb el fastasubseq ja podem extreure la regió genòmica d'interés i seqüències flanquejants upstream i downstream:

Instrucció	Detalls
$ fastasubseq /cursos/BI/genomes/2011/genome.fa start length > genomic.fa	$ on start fa referència a la posició de començament del hit que estem analitzant i length fa referència a la quantitat de nucleòtids que volem extreure a partir d'aquella posició

Tornar

Un cop tenim en un arxiu específic la regió genòmica on es localitza el hit que volem comprovar, executem Exonerate o Genewise, programes d'anàlisi de seqüències que, entre d'altres, permeten fer una anotació¹ a partir d'un alineament precís de la proteïna query i la seqüència genòmica extreta.
¹Entenem per anotació la identificació de tots els elements funcionals continguts en una seqüència genòmica.

Tornar

El software TCofee fa alineaments globals múltiples usant un mètode progressiu que aparella seqüències similars i usa el resultat com a punt de partida pel següent aparellament. L'alineament final té una puntuació que ens dóna idea de la fiabilitat del resultat.

Tornar

En aquest treball executem repetidament un joc de seqüències de comandes sobre diferents arxius. En el nostre cas a més, havíem de realitzar el mateix procés amb moltes queries diferents. Per aquest motiu hem cregut convenient l'automatització dels procediments amb scripts del shell d'UNIX. Per fer-ho, primer cal analitzar l'estructura d'arxius que usem per tal de crear cicles que rastregin les carpetes indicades i executin les comandes sobre els arxius que trobin. Redirigim els resultats a arxius dins de carpetes específiques, de manera que sempre ens sigui possible accedir a aquestes dades per tal d'usar-les per altres tasques o per extreure conclusions. Les taules de resultats també es generen automàticament a partir dels arxius obtinguts.

A continuació es detallen els diferents scripts creats, quina part del procés fan i quins arxius generen. Algunes característiques no s'exposen a la descripció i estan explicades com a comentaris a l'interior de l'script

rutes.sh

Aquest programa carrega de cop totes les rutes necessàries per executar els programes que no es troben en el nostre ordinador local sinó que estan carregats al servidor.

#!/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

blast.sh

En aquest script automatitzem tots els BLAST que hem de realitzar. Per a totes les proteïnes en la carpeta Llibreria fa un blast contra cada un dels genomes (la ruta es defineix com a variable).
Els resultats s'emmagatzemen de dues maneres. Per un banda, es crea l'estructura de directoris següent: dins de la carpeta hits, cada proteïna té una subcarpeta, i dins de cada una d'elles hi ha una altra carpeta per cada genoma que hagi produït hits significatius. Cada un d'aquests hits té un arxiu propi (hitX.txt) amb una sola línia, que inclou tots els detalls del hit. Per exemple, un hit de GPx1 en H.arabidopsis es guardaria així: hits/GPx1/H.arabidopsis/hit0.txt. També generem un arxiu recopilatori de tots els hits obtinguts.
Per altra banda, es genera un arxiu PROTEINA_output_blast_NOMdelGENOMA.txt amb tots els hits de cada genoma per a cada proteïna detallats, de manera que es poden observar els alineaments. Aquests resultats es guarden dins del directori output_gran_blast.

#!/bin/bash

#Definim la ruta on es troben els genomes

adreça="/cursos/BI/genomes/protists/2011"

#Detectem quines proteïnes hi ha a la carpeta Llibreria

for proteines_fa in ./Llibreria/*; do {
basename_fa=`basename $proteines_fa`
proteina=${basename_fa%.fa}
proteines=$proteines" $proteina"
} done

#INICI: Dos bucles per a fer les combinaciones entre proteïnes i genomes

rm resultats_blast.txt 2> error.txt

       
       for protein in $proteines ; do {
               for genomes 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 {

                       patroalfa=`cat ./Llibreria/$protein.fa | egrep ^">"`
                       patrobeta=${patroalfa#*>}
                       echo $patrobeta > temppatro
                       patro=`cut -f 1 -d " " temppatro`
                       rm temppatro
			echo patro es $patro
                       
			#Per a cada cas realitzem el blast de dues maneres, una perquè ens doni els resultats resumits (només les dades de cada hit) i una segona per poder observar els alineaments detalladament.
                       blastall -p tblastn -i ./Llibreria/$protein.fa -d /cursos/BI/genomes/protists/2011/$genomes/genome.fa -o output_blast/$protein.$genomes.txt -m 9 -F F

			blastall -p tblastn -i ./Llibreria/$protein.fa -d /cursos/BI/genomes/protists/2011/$genomes/genome.fa -o output_gran_blast/$protein.$genomes.txt -F F
			
			#Generem un arxiu (hit) per a cada hit creat. Generem també l'estructura de directoris
                       echo ${protein}_$genomes: >> resultats_blast.txt
                       egrep "$patro" output_blast/$protein.$genomes.txt | egrep e- >> resultats_blast.txt
			egrep "$patro" output_blast/$protein.$genomes.txt | egrep e- > temp
			x=`egrep "$patro" output_blast/$protein.$genomes.txt | egrep e-`			
			if [ "$x" != "" ]; then			
			mkdir hits 2> error.txt
			mkdir hits/$protein 2> error.txt
			mkdir hits/$protein/$genomes 2> error.txt
			./uni.pl $protein $genomes < temp 
			fi
                       echo ${protein}_$genomes fet
               } done
       } done

exonerate.sh

Aquest és el programa més llarg i que automatitza més accions. La funció general del programa és que per a cada hit obtingut, EXONERATE busqui possibles proteïnes en una regió de 30000bp que inclogui el hit. Ja hem comentat abans els passos que cal realitzar per fer-ho i per tant l'únic que fem en aquest script és encadenar-les totes perquè ho faci per a cada hit obtingut.

#!/bin/bash

#Escollim si volem usar o no un mètode exhaustiu en l'exonerate

if [ "$1" = "exhaustive" ]; then
exhaustive="--exhaustive yes"
fi

#Borra l'arxiu previ de resultats
rm resultats_exonerate.txt 2> errors.txt

#Detecció de cada proteïna
for proteins in hits/* ; do {
	#Detecció de cada genoma
	for genome in $proteins/* ; do {
		#Detecció de cada hit
		for hit in $genome/*; do {
		
		#Obtenim els noms senzills de cada proteïna, genoma i hit (sense la ruta)
		contig=`egrep SPP $hit| cut -f2`
		proteins=`basename $proteins`
		genome=`basename $genome`
		hit=`basename $hit`
		
		#Creem les regions per a cada hit, els hits del mateix contig es sobrescriuen			
		fastafetch /cursos/BI/genomes/protists/2011/$genome/genome.fa /cursos/BI/genomes/protists/2011/$genome/genome.index "$contig" > regions/DNA_${contig}_$genome.fa
		echo "region $genome done"
		
		#Identifiquem els punts d'inici i final, a partir de les dades del hit

		inici=`cat hits/$proteins/$genome/$hit | cut -f 9`
		
		final=`cat hits/$proteins/$genome/$hit | cut -f 10`
		
		if [ $inici -gt $final ]; then
		inici=$final
		fi		
		start=$(( $inici - 20000 ))
		length="30000"
		fastasubseq regions/DNA_${contig}_$genome.fa $start $length > subseq/${proteins}_${genome}_${hit%.txt}_genomic.fa 2> error_msg.txt	
		echo "subregion $genome done"	

	#En cas que els llocs d'inici o final definits es surtin del contig corresponent, es generen de nou adaptan-se a la posició on es troben
					
	error=`cat error_msg.txt`
	if [ "$error" != "" ]; then
			error_llarg=`cat error_msg.txt | egrep "before end"`
			error_curt=`cat error_msg.txt | egrep "Unknown flag"`
						
		#Si és massa llarg, establim que el punt final és el final del contig
		if [ "$error_llarg" != "" ]; then
			llargada=$error_llarg						
			llargada=${llargada##*(}
			llargada=${llargada%)}
			length=$(( $llargada - $start))
			fastasubseq regions/DNA_${contig}_$genome.fa $start $length > subseq/${proteins}_${genome}_${hit%.txt}_genomic.fa
		fi
						
		#Si l'inici està fora del contig, establim que l'inici és la posició zero del contig
		if [ "$error_curt" != "" ]; then
			start="0"
			fastasubseq regions/DNA_${contig}_$genome.fa $start $length > subseq/${proteins}_${genome}_${hit%.txt}_genomic.fa
		fi
	fi


	#EXONERATE

	exonerate -m p2g --showtargetgff -q Llibreria/$proteins.fa -t subseq/${proteins}_${genome}_${hit%.txt}_genomic.fa --exhaustive yes > exonerate_output/exonerate_${proteins}_${genome}_${hit%.txt}.gff 2> error_de_substitucio.txt
	
	echo "exonerate $genome done"
	echo
	#Substituïm les U de la seqüència query per X, en cas de que hi siguin. (Exonerate no entén la U).				

	error_substitucio=`cat error_de_substitucio.txt | egrep "Unknown amino acid"`
	if [ "$error_substitucio" = "** FATAL ERROR **: Unknown amino acid [U]" ]; then
			if [ "$int_UX" = "0" ]; then
			./substitucio.pl < Llibreria/$proteins.fa > ${proteins}_temporal.txt
			int_UX="1"
			fi
	exonerate -m p2g --showtargetgff -q ${proteins}_temporal.txt -t subseq/${proteins}_${genome}_${hit%.txt}_genomic.fa --exhaustive yes > ./exonerate_output/exonerate_${proteins}_${genome}_${hit%.txt}.gff 2> error_de_substitucio2.txt		
	echo "exonerate subs $genome done"		
			fi

	#Resultats obtinguts amb Exonerate
	gene=`grep -w gene exonerate_output/exonerate_${proteins}_${genome}_${hit%.txt}.gff`
               if [ "$gene" != "" ]; then
                       echo ${proteins}_${genome}_${hit%.txt} : >> resultats_exonerate.txt
               fi

	#Generem un arxiu temporal amb els resultats amb exon de cada output
	grep -w exon exonerate_output/exonerate_${proteins}_${genome}_${hit%.txt}.gff > temp.exonerate.gff
	reso=`cat temp.exonerate.gff`
	
	#Si hi ha algun exó, generem la sequencia FASTA d'aquella regió
	if [ "$reso" != "" ]; then
	#Obtenim les dades de la regió que conté l'exó
	./fastaseqfromGFF.pl subseq/${proteins}_${genome}_${hit%.txt}_genomic.fa temp.exonerate.gff > cDNA/${proteins}_${genome}_${hit%.txt}_exons.fa
	#Traduïm en les 6 pautes de lectura la proteina possibles d'aquella regió	
	fastatranslate cDNA/${proteins}_${genome}_${hit%.txt}_exons.fa > proteins/${proteins}_${genome}_${hit%.txt}_proteins.fa
	fi


} done
echo $genome done
echo
} done
echo $proteins done
echo
} done

genewise.sh

Aquest script fa un procés igual que l'anterior, però enlloc de fer-ho amb l'Exonerate, ho fa amb el Genewise. De tota manera utilitzem les subseqüències generades per l'anterior script, per tant sempre s'han d'executar en aquest ordre.
El Genewise fa una anàlisi una mica diferent que l'Exonerate, per això sovint les proteïnes predites són diferents. A més el Genewise quan troba regions exoǹiques, ens produeix dues proteïnes, una en cada pauta de lectura. Per això afegim un programa perl en aquest script que ens selecciona la proteïna correcta (agafa sempre la més llarga ja que en la major part de casos la segona proteïna produida és molt curta ja que té un codó stop molt proper a l'inici). Com que vàrem veure que en la majoria de casos el resultat del Genewise era millor que el del Exonerate utilitzem la proteïna predita pel Genewise per fer un alineament automàtic utilitzant TCoffee entre la query i la proteïna predita. Per últim busquem elements SECIS mitjançant un programa perl que ens vem descarregar de la web de bioinformàtica i que hem afegit també a la nostra carpeta arrel.

#!/bin/bash

#rm genewise/*.fa
#rm genewise/*.txt

#Bucles que escanejen tots els hits per a cada organisme i proteïna (igual que Exonerate)
for proteins in hits/* ; do {
       for genome in $proteins/* ; do {

               for hit in $genome/*; do {
		
		proteins=`basename $proteins`
        	genome=`basename $genome`
        	hit=`basename $hit`
		hit=${hit%.txt}
		
		echo $proteins
		echo $genome
		echo $hit
		
		#Si tenim subseq per aquell hit
		if [ -f subseq/${proteins}_${genome}_${hit}_genomic.fa ]; then
			#Realitzem anàlisi amb genewise: -pep mostra únicament com a resultat la proteïna i -both ens genera la proteïna en forward i reversa
			genewise -pep -both Llibreria/$proteins.fa subseq/${proteins}_${genome}_${hit}_genomic.fa > genewise/genewise_${proteins}_${genome}_${hit}.txt	
	
			#Seleccionem la proteïna més llarga
			prot=`cat genewise/genewise_${proteins}_${genome}_${hit}.txt`
		
			if [ "$prot" != "" ]; then
		
				./llarg.pl < genewise/genewise_${proteins}_${genome}_${hit}.txt > genewise/PEP_${proteins}_${genome}_${hit}.fa
			fi


			echo ${proteins}_${genome}_${hit} fet

#t_coffee 1aab_ref1.pep -run_name / >/dev/null 2>/dev/null
			#Realitzem el T-coffee entre la query i la proteïna predita
			t_coffee Llibreria/$proteins.fa genewise/PEP_${proteins}_${genome}_${hit}.fa > t_coffee/restcoffee_${proteins}_${genome}_${hit}.txt
		
		fi

		#SECIS search-Busca elements SECIS en la subseq que hem utilitzat per predir la proteïna
		
		./SECISearch.pl < subseq/${proteins}_${genome}_${hit}_genomic.fa > SECIS/${proteins}_${genome}_${hit}.txt



		} done

		echo all hits of $genome done
		echo
	} done
} done

Programes Perl

Com ja hem comentat en alguns dels scripts anteriors, hem creat petits programes perl per automatitzar algunes accions repetitives. A continuació es mostra el codi de cada un d'ells, la funció de cada un d'ells ja ha estat explicada abans.

uni.pl

#!/usr/bin/perl

use strict;
my $contador=0;

while ()
{
	open (openfile,">hits/$ARGV[0]/$ARGV[1]/hit$contador.txt");
	print openfile $_;
	close openfile;
	$contador++;
}

llarg.pl

#!/usr/bin/perl -w

#PROGRAMA QUE MESURA LES DUES SEQÜÈNCIES PROTEIQUES I DÓNA LA MÉS LLARGA
#CAL QUE L'INPUT SIGUI EN FORMAT GENEWISE

use strict;

my @v;
my $inici=0;
my $longv=0;
my $c1=0;
my $par1=0;
my $par2=0;
my $temp1="";
my $temp2="";


while ()
{
	@v=split(//,$_);
	$longv=scalar(@v);
	$c1=0;
	while ($c1 < $longv)
	{
		if ($v[$c1] eq ">")
		{
			$inici++;
		}
		
		if ($inici==2)
		{
			$temp2.=$v[$c1];
			$par2++;		
		}
		else
		{
			$temp1.=$v[$c1];
			$par1++;
		}
		$c1++;
		
	}

}

if ($par1 > $par2)
{
	print $temp1;
}
else
{
	print $temp2;
}

substitucio.pl

#!/usr/bin/perl

use strict;

my @v;
my $string;
my $counter = 0;

while ()
{
	if ($counter == 0)
	{
		print $_;
	}
	if ($counter > 0)
	{

		@v=split(//,$_);
	
		while (<@v>)
		{
			if ($_ eq "U")
			{		
				$string.= "X";
			}		
			else
			{
				$string.=$_;
			}
		}
	}
	$counter++;
	$string.="\n";
}

print $string;

Tornar