La cerca de selenoproteïnes en un genoma és bastant complexa ja que, malgrat disposem de programes de predicció de gens, no sempre s'obtenen bons resultats. Per aquest motiu, doncs, cal una visió a l'hora d'abordar els problemes que van sorgint al llarg del procés analític.

A continuació es detalla el protocol seguit per a obtenir les selenoproteïnes d'Aureococcus anophagefferens.

1. Obtenció del genoma d'Aureococcus anophagefferens
2. Obtenció de la seqüència de selenoproteïnes
3. Cerca de similaritat: BLAST
4. Extracció de la regió genòmica que potencialment conté el gen que busquem
5. Anotació del genoma mitjançant Exonerate
6. Anotació del genoma mitjançant Genewise
7. Predicció de gens
8. Cerca de regions conservades
9. Cerca de dominis
10. Cerca d'elements SECIS
11. Anàlisi del tRNA
12. Cerca de noves selenoproteïnes

El genoma d'aquest organisme va ser aportat pels professors de l'assignatura de Bioinformàtica, ja que no es troba disponible en les bases de dades més habituals com el NCBI o GeneDB.

Per a obtenir la seqüència de les selenoproteïnes conegudes vam recórrer a la base de dades SelenoDB i a diverses referències bibliogràfiques sobre selenoproteïnes a prostistes. Juntament amb les selenoproteïnes, també es va analitzar la seva maquinària de síntesi.

Donat que Aureococcus anophagefferens es troba filogenèticament bastant allunyat de la resta d'organismes, vam decidir partir de l'anàlisi de totes les selenoproteïnes conegudes que estaven disponibles. Tots els arxius van ser descarregats individualment en format FASTA (.fa).

Per tal de cercar els possibles alineaments entre el nostre genoma i les selenoproteïnes conegudes, vam utilitzar el programa tBLASTn, que ens compara DNA (seqüència nucleotídica) contra proteïnes (seqüència aminoacídica).

Vam haver de canviar les U, corresponents a la selenocisteïnes, per X, ja que el programa no la reconeix i dóna errors, aquest canvi també ens serà útil per treballar amb els següents programes.

Per treballar amb aquest programa, vam escriure les següents comandes a la terminal:

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

Després es va procedir a fer una base de dades de BLAST del genoma d'Aureococcus anophagefferens a partir del fitxer que conté aquest genoma amb la següent ordre:

$ formatdb -i /disc8/genomes/A.anophagefferens/genome.fa -p F -n genome.fa

Per fer anar el programa amb la seqüència aminoacídica de cada selenoproteïna contra el genoma d'Aureococcus anophagefferens:

$ blastall -p tblastn -i SelenoproteinaX.fa -d genome.fa -o tblastn_SelenoproteinaX.fa

En els casos en que l'alineament no mostrava cap hit significatiu, per tal de confirmar realment que no hi havia homologia, es va procedir a eliminar un dels filtres que es troba per defecte en el programa tBLASTn: el query filter. Aquest filtre exclou les seqüències amb regions de baixa complexitat composicional per tal d'eliminar els hits potencialment confusos:

$ blastall -p tBLASTn -d SelenoproteinaX.fa -i genome.fa -F F

A més a més, ens calia saber les posicions dels alineaments trobats. Aquestes posicions serveixen per a estimar la regió cromosòmica que ocupa el gen que codifica per a cada selenoproteïna, és a dir, per a saber on comença i on acaba. Per això, es va utilitzar la següent comanda, que mostra els alineaments ordenats i tabulats, mostrant informació de cadascun:

$ blastall -p tblastn -i SelenoproteinaX.fa -d genome.fa -m 9 > SelenoproteinaX_posicions.fa

A partir de les dades de localització cromosòmica obtingudes amb comanda anterior, i tenint en compte només els hits significatius, es va extreure la regió genòmica on l'alineament amb el tBLASTn ens indicava que es trobava la selenoproteïna:

$ fastafetch A.anophagefferens/genome.fa A.anophagefferens/genome.index "scaffold_X" > scaffoldX.fa

Un cop ja hem delimitat el scaffold on tenim la proteïna, dins d'aquest cal acotar millor la regió on es troba el gen d'interès. Amb una seqüència més curta podrem treballar amb els gens amb més facilitat. Per assegurar-nos que agafem el gen sencer, tallarem 3000 nucleòtids més abans del inici del alineament que hem obtingut amb tblastn i 3000 nucleòtids més al final. Fem aquesta suposició perquè, en un organisme com el nostre no s'esperen introns llargs. Vam utilitzar el Fastasubseq per tallar la regió, utilitzant la següent comanda:

$ fastasubseq scaffold_X.fa inici longitud > nomdelaproteina_scaffoldX_seqDNA.fa

Com a nucleòtid d'inici vam agafar el més petit independentment de si l'alineament es produïa en sentit forward o reverse.

Tenint ja la seqüència genòmica on es troba el gen d'interès, el que fem amb aquest programa és tornar a fer un alineament més precís de les seqüències proteiques dels diferents organismes amb la seqüència proteica resultant de traduir el nostre DNA en els diferents frames.

Abans de fer res caldrà revisar les seqüències de les proteïnes que volem comparar, ja que hem vist que Selenodb posa símbols al final que provoquen un error en el programa. Per a realitzar la comanda haurem d'escriure el següent:

$ exonerate -m p2g --showtargetgff -q proteïna.fa -t nomdelaproteïna_scaffoldX_seqDNA.fa > exonerate_nomdelaproteïna_scaffoldX.fa

Genewise també alinea les proteïnes dels diferents organismes en la regió genòmica acotada. La diferència amb el mètode anterior és que Genewise només compara una seqüència d'un organisme en cada comanda. Caldrà doncs separar els arxius en els que tenim les seqüències d'una proteïna dels diferents organismes en tants arxius com organismes tinguem. Abans d'executar el programa caldrà posar la següent comanda:

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

Es important remarcar que també utilitzem aquesta anotació perquè en alguns casos no s'obtenen resultats mitjançant el software Exonerate. Cal tenir en compte que Genewise, a diferència de l'Exonerate, sí que és sensible a la direccionalitat del possible alineament, així doncs caldrà fer una comanda segons tinguem strand forward o reverse.

Forward:

$ genewise -pep -pretty -cdna -gff proteïna_org.fa nomdelaproteïna_scaffoldX_seqDNA.fa > genewise_nomdelaproteïna_scaffoldX_organisme.fa

Reverse:

$ genewise -pep -pretty -cdna -gff -trev proteïna_org.fa nomdelaproteïna_scaffoldX_seqDNA.fa > genewise_nomdelaproteïna_scaffoldX_organisme.fa

Arribats a aquest punt disposem dels resultats obtinguts amb el tblastn, l'Exonerate i el Genewise, els observem per cercar selenoproteïnes o homòlegs en cisteïna.

Dins l'exonerate sabrem si hi ha selenoproteïnes quan trobem en la seqüència subject, la d'Aerococcus anophagefferens, tres asteriscs. En les seqüències query veurem Ukn enlloc de la selenocisteïna. En Genewise veurem una X allà on tinguem una selenocisteïna.

Els homòlegs amb cisteïna es caracteritzaran per tenir una cisteïna en el lloc on altres seqüències tenen una selenocisteïna. També podem trobar proteïnes homòlogues amb altres aminoàcids. En el cas de trobar una selenoproteïna o bé una proteïna homòloga caldrà extreure'n la seqüència més llarga possible, utilitzant per això el resultat obtingut amb els tres programes, complementant les seqüències. Aquests arxius els desarem amb el següent nom:

Nomdelaproteïna_scaffoldX.fa

En el cas que la predicció no sigui gaire bona, perquè observem que és molt curta comparada amb les que tenim dels altres organismes, caldrà que fem un blastp (Protein Blast) contra la base de dades del NCBI. D'aquesta manera, obtindrem una seqüència de la proteïna que busquem més semblant al fragment que nosaltres hem predit anteriorment.

Utilitzarem la seqüència proteica del organisme que no havíem utilitzat per tornar a realitzar els passos 5 i 6. Si la seqüència més propera que ens proposa NCBI és una que ja havíem utilitzat significa que el fragment obtingut és la millor predicció a la que podem arribar amb aquests programes.

Un cop tenim la seqüència allargada fem un alineament múltiple entre la nostra selenoproteïna i totes les selenoproteïnes conegudes utilitzant el programa Tcoffee. El resultat d’aquest alineament ens dirà si la seqüència de la nostra selenoproteïna ha divergit molt o poc respecte les altres ja que ens mostrarà les regions que presenten homologia en vermell, és a dir, ens mostrarà les regions conservades entre les diferents espècies.

Per comprovar que la selenoproteïna que hem trobat en el nostre genoma realment és aquella selenoproteïna coneguda que estavem buscant farem ús del programa Interpro. En aquest programa introduirem la seqüència d'aminoàcids que tenim i farem una cerca de a quins dominis correspon. El tipus de domini que aparegui als resultats hauria de concordar amb la selenoproteïna que estavem buscant en el nostre organisme.

Quan ja tenim la selenoproteïna hem de confirmar que el codó TGA s'està traduint realment com a selenoproteïna i no com un codó d'STOP dins la nostra seqüència. Per esbrinar-ho busquem la presència d'elements SECIS, ja que aquests són els elements responsables d'aquest canvi en la traducció, dins l'scaffold que conté la selenoproteïna utilitzant el programa SECISearch.

En el cas d'obtenir més d'un element SECIS dins de l'scaffold a on es troba la nostra seqüència hem de determinar quin és el que correspon a la nostra selenoproteïna o almenys establir quins són els millors candidats. Per fer-ho ens fixem en quina posició es troba, és a dir, mirar quin és el SECIS més proper a la selenoproteïna en sentit “dowstream” (en el cas de que la nostra selenoproteïna sigui “reverse” buscarem el SECIS en sentit “upstream”). Tant si mirem en un sentit com un altre els SECIS que poden codificar per la nostra selenoproteïna es trobaran a una distància màxima d'uns 1000 nuclèotids. Un cop seleccionem aquells possibles candidats per proximitat podem determinar quins són els millors a partir del valor de la seva energia (interesa un valor elevat) i de la seva conformació (que presenti els dos loops i les bases característiques conservades).

Per tal de trobar el tRNASec en el genoma d'A.anophagefferens, es va obtenir la seqüència genòmica que codifica pel tRNASec en diferents organismes en el tRNADB i es va procedir a realitzar un alineament (DNA-DNA) amb BLASTn de la següent manera:

$ blastall -p blastn -i tRNA.fa -d genome.fa -o blastn_tRNA.fa

Un cop identificat el millor alineament i tallada la seqüència del scaffold corresponent tal i com es va realitzar per cercar les selenoproteïnes conegudes, es va realitzar un segon alineament amb el programa Exonerate utilitzant la comanda (per comparar DNA-DNA):

$ bash-3.2$ export PATH=$PATH:/disc8/bin/exonerate/bin/

$ exonerate --model est2genome --bestn 10 tRNA.fa seqDNA_tRNA_scaffold14.fa

Primer de tot, s'han de buscar en el nostre genoma tots els possibles elements SECIS que pot contenir. Això es realitza amb el programa SECISearch.pl, i la comanda seria la següent:

$ SECISearch.pl -p 's' genome.fa

El símbol 's' ens indica que realitzarà una cerca estàndard. Després d'obtenir tots els elements SECIS, es procedeix a seleccionar els òptims segons la seva energia lliure, i en el nostre cas, s'han seleccionat aquells que tenen una energia lliure inferior a -20 amb aquesta comanda:

$ SECISearch.pl -p s -e -20 -I genome.fa > secis20.fa

A partir d'aquí, es guarda el nom de cada scaffold que conté SECIS, tant de l'strand positiu com en el negatiu. Simplement, es busquen les línies amb “>” que són les que contenen el nom de l'scaffold. Per SECIS en l'strand positiu es porta a terme aquesta comanda:

$ egrep ">" genome.fa.std.secis | grep -v complemen | gawk -F: '{print $1}' | sed 's/>//' | sort | uniq > nombres_mas

on “nombre_mas” correspon a la carpeta on es redireccionen els scaffolds amb aquests SECIS. Amb SECIS en l'strand negatiu es fa servir la següent:

$ egrep ">" genome.fa.std.secis | grep complemen | gawk -F: '{print $1}'| sed 's/>//' | sort | uniq > nombres_menos

on es redirrecionen els scaffolds amb aquests SECIS a la carpeta “nombres_menos”. Ara ja ho tenim tot distribuït en les carpetes corresponents, i s'ha de procedir a extreure la seqüència dels scaffolds que ens interessen. Això es realitza mitjançant el programa retrieveseqs.pl amb la següent comanda:

$ perl retrieveseqs.pl -vf /disc8/genomes/A.anophagefferens/genome.fa nombres_mas > contigs_mas.fa

$ perl retrieveseqs.pl -vf /disc8/genomes/A.anophagefferens/genome.fa nombres _menos > contigs_menos.fa

cada una per la respectiva carpeta. Així s'obtenen els scaffolds d'interès i s'emmagatzemen en els arxius de “contigs_mas/menos.fa” en format FASTA. Per simplificar, es canvien els scaffolds amb SECIS en l'strand negatiu al positiu seguint aquesta comanda:

$ /disc8/bin/exonerate/bin/fastarevcomp contigs_menos.fa > contigs_menos.revcomp.fa

I per concatenar-los tots en un únic fitxer anomenat “contigs_todos.fa” es fa:

cat contigs_mas.fa contigs_menos.fa > contigs_todos.fa

Després de realitzar tot aquest procediment, s'ha de crear una carpeta anomenada “salidas” perquè es puguin redireccionar tots els fitxers obtinguts després d'haver extret 500 nucleòtids de cada SECIS predit. Per els que tenim en l'strand positiu es fa la següent comanda:

grep ">" secis20.fa | grep -v comple | perl -ne '/>(.+?):\[(\d+)/; my $a=$2-500; my $b=$2; my $name=$1; $name=~/(\d+)/; my $outname=$1; if ($b>=500){ system("perl retrieveseqs.pl -vfn contigs_todos.fa \"$name\" > contig_temp.fa; /disc8/bin/exonerate/bin/fastasubseq -s $a -l 500 contig_temp.fa> salidas/$outname.$b.subseq.fa") }'

i per el negatiu:

grep ">" secis20.fa | grep comple | perl -ne '/>(.+?):\[\d+,(\d+)/; my $a=$2-500; my $b=$2; my $name=$1; $name=~/(\d+)/; my $outname=$1; if ($b>=500){ system("perl retrieveseqs.pl -vfn contigs_todos.fa \"$name\" > contig_temp.fa; /disc8/bin/exonerate/bin/fastasubseq -s $a -l 500 contig_temp.fa> salidas/$outname.$b.menos.subseq.fa") }'

Si tot ha sortit com hauria, en la carpeta “salidas” ha d'haver molts fitxers per cada element SECIS anomenats “scaffolds_12345.678.subseq.fa” on 12345 és el número de l'scaffold i 678 el punt de començament del SECIS. També tenim fitxers *.menos que es refereixen als SECIS en l'strand negatiu, però que ja s'han canviat anteriorment.
A partir d'ara, s'assumeix que en el cas que el SECIS predit és real, aquest tindrà algun exó codificant i que conté selenocisteïna en els primers 500 nucleòtids al SECIS. Per això, s'han de traduir aquestes seqüències en les tres pautes de lectura i fer un BLAST contra la base de dades de NCBI per poder buscar seqüències conservades. Només s'han de mirar aquelles que no corresponguin a la nostra espècie a estudiar, ja que correspondrà a les selenoproteïnes que s'han trobat. Es volen trobar hits contra altres espècies que tinguin una U de la query alineada amb un * o una C en el subject. Això serà indicatiu que aquesta seqüència està conservada, i per tant, pot ser codificant.
Per poder traduir-les en les tres pautes de lectura es fa servir el programa trans.pl que tradueix per defecte en les tres pautes (no en sis, perquè no mira l'strand complementari) i po sa el codó TGA com una U i els STOP amb *. La comanda per traduir és la següent:

for n in $(/bin/ls salidas/*fa); do echo "traduciendo $n..."; perl trans.pl $n > $n.pep; done

i per concatenar tot el fitxer:

cat salidas/*pep > salidas/contigs_todos.pep

Ara ja sí, es porta a terme el BLAST contra NCBI, que és el blastcl3, amb aquesta comanda:

blastcl3 -p tblastn -i salidas/contigs_todos.pep -d nr > contigs-nr.out

Una vegada obtingut el fitxer del blast, es procedeix a revisar tots els alineaments i buscar en tots els organismes que no correponguin al nostre una alineació d'U amb * o C. Això es realitza amb un programa que selecciona aquest tipus d'alineaments que és el tblastn_hits.pl . Finalment, s'hauria de realitzar el genewise amb el millor hit obtingut del Blastcl3 contra l'scaffold d'on s'ha extret. Això ens donarà la predicció d'una proteïna en el genoma d'A. anophagefferens i després, s'ha de fer un BlastP d'aquesta contra la base de dades de NCBI per veure si troba proteïnes conegudes i descobrir alguna selenoproteïna nova.