L'objectiu d'aquest treball ha estat identificar i localitzar les selenoproteïnes existents en el genoma d'Aotus nancymaae, així com la maquinària relacionada amb la seva síntesi. Per tal d'elaborar el projecte, s'han utilitzat com a queries les selenoproteïnes anotades en el genoma d'Homo sapiens, l'espècie més propera a Aotus nancymaae filogenèticament i les selenoproteïnes de la qual estan més ben descrites. S'han utilitzat diferents programes a través del sistema operatiu Ubuntu.
A continuació, presentarem una explicació detallada sobre les comandes i el procediment seguit per tal d'assolir una correcta anotació de selenoproteïnes en el genoma d'Aotus nancymaae.
1. Obtenció de les queries
Totes les proteïnes incloent selenoproteïnes, homòlegs amb cisteïna o altres aminoàcids i maquinària de biosíntesis han estat extretes de la base de dades SelenoDB d'Homo sapiens.
Per tal de realitzar l'anàlisi, s'han hagut de modificar aquelles queries que contenien una selenocisteïna, representada amb una U, ja que alguns dels programes utilitzats no eren capaços d'interpretar-la. Per tant, la U ha estat substituïda per una X. També s'han hagut d'eliminar els senyals de finalització de la seqüència protèica #,@ i %.
2. Obtenció del genoma d'Aotus nancymaae
El genoma de l'espècie a analitzar, ésa dir, de l'Aotus nancymaae, ha estat proporcionat pels professors de l'assignatura. Aquest es troba en el següent directori:
| /cursos/BI/genomes/2015/Aotus_nancymaae/genome.fa |
3. BLAST
Per tal de localitzar les selenoproteïnes en el genoma d'Aotus nancymaae, el primer que s'ha de fer és localitzar regions en aquest genoma que siguin susceptibles de contenir la seqüència codificant per la query. Amb aquest objectiu, s'utilitza el programa BLAST (Basic Local Alignment Search Tool).
Aquest és un programa informàtic que alinea localment una seqüència problema (query) amb seqüències provinents d'una base de dades (en aquest cas, seqüències provinents del genoma d'Aotus nancymaae). Per tal de realitzar l'alineament, el programa segueix un algoritme heurístic que selecciona ràpidament les seqüències (hits) amb més homologia respecte la seqüència problema. S'ha de tenir en compte, que el fet d'utilitzar un algoritme heurístic pot comportar la pèrdua de regions amb homologia real que no presenten una similitud molt elevada amb la query introduïda.
Hi ha diferents tipus de BLAST que permeten alinear diferents tipus de seqüències (DNA, RNA o proteïna) amb una base de dades, d'entre els quals l'utilitzat ha estat tblastn, que permet comparar una query formada per una seqüència proteica amb una base de dades nucleotídica.
L'output proporcionat per aquest programa consisteix en una llista de hits amb el seu e-value i el seu score, valors que s'hauran de tenir en compte per tal de seleccionar els hits òptims. Com menor sigui l'e-value i major el score, més probabilitats hi ha de que el hit trobat contingui la seqüència codificant per la query problema i, per tant, menys probabilitats de que el hit obtingut hagi estat donat per atzar.
A l'hora d'executar aquest programa, es pot definir un valor d'e-value mínim a partir del qual mostri els hits. En aquest projecte, el mínim ha estat definit com a 10-3. Tot i això, a l'hora d'interpretar els resultats, els hits que preferentment s'han tingut en compte són aquells que presentaven un e-value menor de 10-10.
En el cas de tenir múltiples hits amb bons e-values, van ser analitzats tots ells per tal de després poder destriar el més adequat a partir de la comparació dels resultats amb totes les queries. Concretament, en el cas de les famílies de selenoproteïnes analitzades, sovint trobàvem hits en les mateixes regions per a diferents proteïnes; això s'explica per l'alta homologia entre proteines de la mateixa família degut a un orígen filogenètic comú. En aquests casos, ha estat important comparar els resultats del BLAST de tots els membres de la família, juntament amb la interpretació dels següents passos de la cerca, per tal d'assignar la regió més semblant a la query.
Primer de tot, s'havia de definir un "camí" al shell per tal de que trobés la carpeta en la qual es localitzava el programa, i es va fer a través de les següents comandes:
| $ export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH |
| $ cp /cursos/BI/bin/ncbiblast/.ncbirc ~/ |
A partir d'aquí, l'ordre donada al shell per tal d'executar el programa és la següent:
| $ blastall -p tblastn -i query.fa -d/cursos/BI/genomes/2015/Aotus_nancymaae/genome.fa -o output.fa |
| -p: tipus de blast
-i: ubicació de la query -d: ubicació del genoma d'Aotus nancymaae -p: ubicació que volem donar a l'output (resultats del BLAST) output.fa: fitxer amb els resultats del BLAST |
3.1. Extracció de la regió genòmica
Primer de tot, s'ha d'indexar el genoma per tal de poder extreure les regions genòmiques d'Aotus nancymaae on l'alineament proporcionat pel tBLASTn indica que hi ha possibilitat de trobar una selenoproteïna. Per fer-ho, s'utilitza la comanda fastaindex. Aquest programa organitza el fitxer del genoma en diferents regions, és a dir, organitza els fitxers multifasta a una forma indexada per tal de poder extreure les parts a analitzar.
Per tal d'executar aquest programa, l'ordre utilitzada és la següent:
| $ fastaindex /cursos/BI/genomes/2015/Aotus_nancymaae/genome.fa output.index |
| output.index: ubicació que volem donar a l'index. |
Per a la realització d'aquest projecte, els professors de l'assignatura ja ens van donar aquest index fet. A partir de l'index, vam haver d'extreure les regions genòmiques d'interès per a cada query, és a dir, aquelles regions on el resultat del BLAST mostrava com a mínim un hit per a la query. Les regions d'interès extretes, anomenades scaffolds, van ser guardades en un fitxer amb format fasta amb el programa fastafetch. L'ordre utilitzada va ser la següent:
| $ fastafetch /cursos/BI/genomes/2015/Aotus_nancymaae/genome.fa output.index scaffold > scaffold.fa |
| output.index: ubicació de l'index
scaffold: regió que volem extreure scaffold.fa: fitxer fasta que conté la regió extreta |
Aquesta regió extreta pel programa fastafetch encara era massa gran per tal de poder treballar amb comoditat. Per aquest motiu, es va utilitzar el programa fastasubseq per extreure una zona encara més limitada dins d'aquesta regió. Per assegurar-nos de què ens quedàvem amb la regió potencialment codificant per a la proteïna cercada, en aquest programa consideravem la posició més baixa i la més alta d'entre tots els hits trobats en cada scaffold i, a més a més, hi sumàvem un marge de 50.000 nuceòtids tant upstream com downstream. L'ordre utilitzada va ser la següent:
| $ fastasubseq scaffold.fa start length > subseq.fa |
| scaffold.fa: resultat del fastafetch
start: nucleòtid inicial – 50.000 lenght: 100.000 nucleòtids + llargada del hit subseq.fa: fitxer fasta que conté la zona de la regió extreta pel fastafetch que conté el gen |
4. Anotació genòmica
4.1. Exonerate
El següent pas consisteix en assegurar que el hit extret es troba en una regió exònica i que, per tant, codifica per una proteïna. Per aquest motiu, es realitzen les anotacions dels gens d'interès amb el programa Exonerate, alineant el fragment de DNA extret amb fastasubseq amb la query inicial. S'utilitza aquest programa ja que fa alineament de seqüències d’una manera més exacta que el BLAST, i també més informativa ja que indica més característiques del gen (introns, exons, zones d’splicing, etc).
Abans de poder executar el programa, s'ha de definir un altre cop un "camí" al shell per tal que trobi la carpeta en la qual es localitza:
| $ export PATH=/cursos/BI/soft/exonerate/x86_64/bin:$PATH |
L'ordre utilitzada és la següent:
| $ exonerate -m p2g --showtargetgff --exhaustive -q query.fa -t subseq.fa > exonerate.gff |
| -m p2g: indica el model d'alineament (en aquest cas proteïna vs genoma)
--showtargetgff: indica el format amb el que es generarà el fitxer de sortida --exhaustive: comanda opcional (utilitzada en aquest anàlisi) per tal de tenir en compte també els alineaments subòptims -q: ubicació de la query inicial -t: resultat del fastasubseq -exonerate.gff: alineament aconseguit en format GFF |
4.1.1.FastaseqfromGFF
A continuació s'ha d'extreure la seqüència d'exons obtinguda a partir del Exonerate en format fasta, de manera que s'utilitza el programa FastaseqfromGFF, però abans és necessari saber les posicions dels exons per tal de poder-los extreure. Això es duu a terme amb la funció egrep.
L'ordre utilitzada és la següent:
| $ egrep -w exon exonerate.gff > cDNA.gff |
|
egrep: selecciona les línies on aparegui el patró definit
-w exon: patró que es vol buscar (exons) exonerate.gff: fitxer proporcionat per l'exonerate amb la seqüència d'exons cDNA.gff: fitxer de sortida que indica la posició dels exons |
Un cop se saben les posicions dels exons, es pot procedir a extreure'ls amb fastaseqfromGFF. Abans d'executar el programa, però, caldrà donar al shell la seva localització:
| $ export PATH=/cursos/BI/bin:$PATH |
Després s'executa amb la comanda:
| $ fastaseqfromGFF subseq.fa cDNA.gff > cDNA.fa |
|
subseq.fa: fitxer extret pel fastasubseq
cDNA.gff: fitxer de sortida que indica la posició dels exons cDNA.fa: fitxer de sortida en format fasta |
4.1.2. Fastatranslate
Per tal de poder realitzar un tercer alineament amb el programa T-coffee i Genewise caldrà traduir la seqüència extreta amb els passos anteriors a proteïna. Per fer-ho s'utilitza el programa fastatranslate, el qual s'executa amb la següent comanda:
| $ fastatranslate cDNA.fa > translate.mfa |
L'arxiu de sortida (translate.mfa) conté tots els marcs de lectura possibles per la seqüència en format multifasta. L'arxiu conté sis seqüències fasta que corresponen a les sis pautes de lectura possibles (tres forward i tres reverse) de les quals s'ha de seleccionar el marc de lectura que correspongui. En la realització d'aquest treball sempre s'ha tingut en compte només el primer marc de lectura, cosa que s'aconsegueix mitjançant l'ordre -F.
| $ fastatranslate cDNA.fa -F X > translate.fa |
4.2. T-COFFEE
Un cop la seqüència ha estat traduïda a proteïna, és a dir, no té introns i es troba en forma de seqüència d'aminoàcids, ja es pot fer un alineament d'aquesta amb la query i, d'aquesta manera, determinar la presència o absència d'homologia entre les dues seqüències. Per poder dur a terme aquest alineament s'utilitza el programa T-coffee (Tree-based Consistency Objective Function for alignment Evaluation). Aquest programa permet fer alineaments múltiples de proteïna, DNA i RNA.
Abans d'executar el programa, caldrà donar al shell la seva localització:
| export PATH=/cursos/BI/soft/t_coffee/x86_64/bin:$PATH |
S'executa amb la comanda:
| $ t_coffee query.fa translate.fa > t_coffee.fa |
| t_coffee.fa: fitxer de sortida amb l'alineament entre les dues seqüències |
4.3. Genewise
Per tal de contrastar els resultats obtinguts amb l'Exonerate i el T-coffee, es pot utilitzar un programa alternatiu com ara el Genewise, que genera una nova anotació del gen utilitzant un algorisme diferent al d'Exonerate. No obstant, nosaltres no hem utilitzat aquest programa ja que hem considerat suficient la informació obtinguda amb l'execució dels programes anteriors
5. Automatització
Per tal d'agilitzar el procés d'anotació de totes les selenoproteïnes i de la seva maquinària, s'ha desenvolupat un programa que automatitza l'execució de totes les comandes anteriorment descrites per a una sola query, especificada des del shell a l'hora de cridar el programa. El programa, escrit en llenguatge perl, es pot trobar en el següent link.
Al llarg del fitxer s'inclouen explicacions i comentaris sobre el codi desenvolupat per a cada pas del procés. Abans de cridar el programa cal executar des del shell les comandes (exports) especificades anteriorment per tal d'instal·lar tots els programes inclosos en l'automatització.
Degut a l'organització de directoris i fitxers que vam utilitzar des del començament del desenvolupament del projecte, cal córrer el programa des d'una carpeta separada dels directoris on emmagatzarem els inputs i els outputs. Aquesta carpeta s'ha de trobar dins d'un directori que contingui també dues altres carpetes: inputs, on emmagatzemarem un fitxer en format fasta per a cada query amb la seqüència d'aminoàcids de la proteïna humana, i outputs, on es direccionaran tots els fitxers generats pels diferents programes utilitzats. Dins de la carpeta outputs, cal crear una carpeta per a cadascun dels programes: Blast, Fastafetch, Fastasubseq, Exonerate i T-COFFEE. Per facilitar l'organització dels fitxers de cara a la interpretació dels resultats, el programa crearà una carpeta per a cada query dins de cadascuna d'aquestes, on es guardaran tots els fitxers generats per a cada query.
El programa desenvolupat presenta una sèrie de limitacions. En primer lloc, a l'hora de determinar la regió d'interès del genoma, es té en compte la posició més baixa i la més alta d'entre tots els hits trobats en un mateix scaffold, tenint en compte si aquests es troben en sentit forward o reverse. Però en el cas (poc freqüent) que en un mateix scaffold hi hagi hits en les dues direccions, a l'hora de predir la proteïna codificada, només es tindrà en compte un dels dos sentits i, per tant, podem perdre informació rellevant a l'hora de realitzar l'alineament. En els casos en què ens hem trobat amb aquest problema, hem realitzat la cerca de forma manual.
D'altra banda, el fet d'incloure el T-COFFEE en l'automatització presenta certes limitacions a l'hora d'avaluar l'alineament. La informació que ens proporciona aquest programa no sempre és suficient, i cal tenir en compte també els resultats dels passos previs. Per exemple, tant els codons STOP com els codons codificants per SelenocisteÏna són expressats amb un asterisc (*) al fitxer generat pel fastatranslate, i aquest símbol pot generar confusions a l'hora d'alinear la seqüència predita amb la query humana. Per això, de vegades ha estat necessari realitzar el T-COFFEE manualment, havent substituït anteriorment els asteriscs susceptibles de ser una Selenocisteina per U.
6. Cerca d'elements SECIs
Per tal d'assegurar que la predicció de selenoproteïnes és correcta, es pot fer una cerca d'elements SECI en la seqüència, ja que aquests són estructures essencials per a la síntesi de selenoproteïnes i es troben en la regió 3'-UTR dels gens que les codifiquen.
Per tal de predir els elements SECIs, s'han utilitzat dos programes online: SECISearch3 i Seblastian, per poder contrastar els resultats d'un amb els de l'altre. El fitxer que cal donar com a entrada a aquests dos programes és el el document generat pel fastasubseq, el qual conté la seqüència nucleotídica del genoma d'Aotus nancymaae potencialment codificant per la proteïna que busquem.
7. Arbres filogenètics
En el cas de les famílies de proteïnes, existeix una gran homologia entre els seus membres. Això fa que l’única manera de poder assignar la identitat de les proteïnes candidates a ser un membre determinat de la família sigui mitjançant l’elaboració d’un arbre filogenètic. Per tal d’elaborar els arbres de les diferents famílies, hem emprat l’opció One Click del recurs online d’anàlisis filogenètic Phylogeny.fr. Per tal de poder utilitzar aquest recurs, cal crear un arxiu multifasta amb les seqüències d’aminoàcids de totes les proteïnes de la família de l’espècie de referència (Homo sapiens) i de l’espècie a estudiar (Aotus nancymaae). Abans d’introduir l’arxiu a la web cal canviar totes les “U” presents en el fitxer per “X” i eliminar tots els símbols presents al final de les seqüències (%, # i @).