Preparació del genoma
El genoma de Mustela putorius furo el vam obtenir d'una base dades creada pels professors de l'assignatura. El nostre genoma el trobem a:
/cursos/BI/genomes/project_2013/Mustela_putorius_furo/genome.fa
Tot i que la carpeta ja contenia aquest arxiu un pas molt important és la indexació del genoma (arxiu anomenat genoma.index). Aquesta indexació separarà el genoma en Scaffolds, que són segments del genoma separats i ordenats. Tornar a l'inici de la pàginaObtenció i preparació de les QUERYs
L'objectiu del treball és obtenir totes les selenoproteïnes i proteïnes relacionades del genoma de la fura. Per això, en primer lloc buscarem totes les selenoproteïnes anotades en espècies de referència i les enfrontarem contra tot el genoma de la fura. Farem la cerca a la base de dades SelenoDB
D’aquesta base de dades, hem escollit el mamífer més proper, el Mus musculus. Més endavant, i degut a que el genoma humà és un dels millors estudiats, hem buscat la resta de selenoproteïnes que trobàvem a l’humà.
Les seqüències d'aminoàcids de les proteïnes seleccionades seran les QUERY i van ser guardades amb fitxers FASTA anomenats “query.fa”
Tornar a l'inici de la pàginaAlineament amb la realització i anàlisi del BLAST
El Basic Local Alignment Search Tool o BLAST és una eina que ens permet comparar seqüències biològiques, com ara aminoàcids o DNA. Aquest algorisme troba similituds entre dues seqüències, realitzant de forma heurística alineaments locals entre les seqüències donades.
D’entre tots els tipus diferents de BLAST utilitzarem el tBLASTn, el qual permet comparar la seqüència d'aminoàcids (que en el nostre cas serà l’anomenada QUERY) d'una proteïna contra una seqüència de nucleòtids o contra una base de dades de DNA traduïda en tots els possibles marcs de lectura (que en el nostre cas serà el genoma de la fura).
Així, el tBLASTn ens proporcionarà els millors alineaments possibles entre cada QUERY i el genoma. Cada alineament possible ens donarà un hit, amb un E-value i un Score determinats. Quan més petit sigui l’E-value i més gran sigui l’Score menys possibilitats hi haurà que l'alineament s'hagi produit per atzar.
Per poder utilitzar el BLAST haurem d'extreure el software necessari. Això és possible introduint les dues comandes següents al Shell:
$ export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
$ cp /cursos/BI/bin/ncbiblast/.ncbirc ~/
*Com que en el nostre cas ja està fet, no caldrà fer una base de dades de BLAST (amb la comanda formatdb) per al genoma de la fura.Un cop introduïdes aquestes comandes al Shell ja és possible realitzar el BLAST amb la comanda següent:
$ blastall -p tblastn -i query.fa -d /cursos/BI/genomes/project_2013/Mustela_putorius_furo/genome.fa -o fitxerdesortida
on el tipus de blast és el tBLASTn (indicat amb l'ordre “-p”), “query.fa” fa referència al fitxer on es troba la nostra QUERY seleccionada en cada cas (indicat amb l'ordre “-i”), “genome.fa” es refereix al fitxer que conté el genoma amb el qual volem fer la cerca (indicat amb “-d”) i “fitxerdesortida” és el nom que volem donar al fitxer en el qual s'emmagatzemarà la informació del blast (indicat amb l'ordre “-o”).Una vegada obtingut aquest fitxer, serà analitzat escollint els millors E-value i Score. El llindar d’elecció serà aquell imposat per nosaltres. El llindar que nosaltres hem establert per l'E-value és de valors inferiors a 0,0001.
Un cop seleccionats els hits, apuntarem l’identificador de la regió del genoma on es troba aquest i la posició del genoma que ocupa.
Tornar a l'inici de la pàginaSelecció del Scaffold
Una vegada seleccionat el hit, el que volem és passar a tenir només l’ Scaffold on es troba cada alineament en lloc de treballar amb tot el genoma que representa més exigències computacionals. Per això executarem el programa fastafetch amb la següent comanda:
$ fastafetch /cursos/BI/genomes/project_2013/Mustela_putorius_furo/genome.fa genome.index "Scaffold-cromosoma" > cromosoma_query.fa
On primer s’indica la ubicació del genoma i a continuació l’arxiu on es troba la separació en Scaffolds i l’ Scaffold escollit entre cometes. Aquest Scaffold seran les quatre primeres columnes de l’arxiu que hauran estat seleccionades prèviament. Tornar a l'inici de la pàginaSelecció de la regió del hit
Finalment, amb el programa fastasubseq acabem de acotar la regió del hit, a partir de l’arxiu on ja tenim només l’ Scaffold on es troba l’alineament. Per fer això,utilitzarem el paràmetre del blast que ens informa de la possició de la seqüència alineada. Per a no perdre informació i assegurar-nos la presència del gen upstream i downstream, augmentaren el nostre hit en 30.000 nucleòtids per cada extrem.
En la comanda del fastasubseq indiquem primer la ubicació de l’arxiu on es troba l’ Scaffold seleccionat, “start” serà la posició del primer nucleòtid que volem contemplar, “llargada” el nombre de nucleòtids que agafem en total i “regio_query.fa” el nom del fitxer de sortida.
$ fastasubseq cromosoma_query.fa start length > regio_query.fa
Tornar a l'inici de la pàginaAnotació amb el programa Exonerate
El programa exonerate ens permet fer una anotació del gen a partir de l'alineament entre la QUERY i la regió genòmica, és a dir, prediu tots els elements funcionals continguts en una seqüència genòmica.
Per poder utilitzar l'Exonerate necessitem instal•lar uns programes determinats que s’introduiran al Shell la següent manera:
$ export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH
Ara ja podem executar l'Exonerate amb la comanda:$ exonerate -m p2g --showtargetgff -q query.fa -t regio_query.fa > query_info.exonerate.gff
on “p2g” és el model de l'alineament (indicat amb l'ordre “-m”), en aquest cas proteïna versus genoma, “--showtargetgff” inclou el resultat en format GFF al fitxer de sortida, “query.fa” és la nostra query (indicada amb l'ordre “-q”), “regio_query.fa” es refereix al fitxer de sortida del fastasubseq (indicat amb “-t”) i “query_info.exonerate.gff” és el fitxer de sortida. *En el cas que la nostre query contingui una U l’haurem de canviar al fitxer query.fa o ens donarà error.Tornarem a executar l’Exonerate, però aquest cop utilitzant la comanda "egrep" per a seleccionar només els exons, cosa que ens facilitarà el posterior anàlisi. $ exonerate -m p2g --showtargetgff -q query.fa -t regio_query.fa | egrep -w exon > query_exon.exonerate.gff
Amb un programa Perl, fastaseqfromGFF.pl, extraurem la seqüència de cDNA corresponent al fitxer GFF que ens genera l'Exonerate. Per executar aquest programa hem d'introduir primer al Shell la comanda següent:
$ export PATH=/cursos/BI/bin:$PATH
Ara ja podem cridar el Perl, on “query.cdna” serà el fitxer on hi serà el cDNA:$ fastaseqfromGFF.pl regio_query.fa query_exon.exonerate.gff > query.cdna
A continuació, executarem el programa fastatranslate, que el que fa es traduir en els sis possibles marcs de lectura (ORFs) una seqüència de DNA donada a proteïna. Per executar el programa cal la comanda:$ fastatranslate cDNA.fa > query_pautes.fa
Si a més afegim l'ordre “-F 1” a la comanda, el programa ens dóna només la proteïna traduïda en la que es considera la millor pauta de lectura. La comanda quedaria així:$ fastatranslate query.cdna -F 1 > query_regio_seqaa.fa
En general nosaltres ens fixarem en el resultat que ens dóna aquesta segona opció, però en casos en què ens sembli necessari revisarem les altres pautes de lectura. Tornar a l'inici de la pàginaGeneració d'una anotació amb el programa GeneWise
El GeneWise és un programa alternatiu a l'Exonerate, que té com a objectiu fer un alineament més precís que el BLAST entre la QUERY i la regió genòmica obtinguda, per predir l'estructura exònica d'una seqüència nucleotídica donada. En el nostre cas no sempre ha estat utilitzat. Només quan l’anotació amb l’Exonerate no era del tot bona. Per poder executar-lo també necessitem introduir primer comandes al Shell:
$ export PATH=/cursos/BI/bin:$PATH
$ export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg
Per utilitzar-lo el cridarem de la següent forma:$ genewise -pep -pretty –cdna -both -gff query.fa regio_query.fa > query.genewise
on “-pep” ens donarà la seqüència peptídica predita, “-pretty” és per a què mostri l’alineament, “-cdna” per a què mostri la seqüència genòmica alineada, “-both” és per a què analitzi en els dos sentits, forward i reverse, “-gff” per a què mostri la informació en format GFF, “query.fa” és la ubicació de la nostra query, “regio_query.fa” és la subseqüència i “query.genewise.fa” és el nom del fitxer de sortida. Tornar a l'inici de la pàginaAnàlisi de les proteïnes obtingudes: TCoffee
Per últim, utilitzem el programa TCoffee per a fer un alineament global entre la nostra query inicial i la proteïna predita al final del procés. Així, aquest programa ens donarà una idea de com de similars són aquestes dues proteïnes i amb què s'està alineant la selenocisteïna de la query en el genoma que estudiem. A més, ens facilita una puntuació que ens dóna idea de la fiabilitat del resultat.
Per executar TCoffee s’ha d’escriure la següent ordre al Shell:
$ t_coffee query.fa fitxerFASTAfinal
on "query.fa" és la proteïna query i "fitxerFASTAfinal" és la proteïna obtinguda amb l'Exonerate o el GeneWise segons cada cas. Tornar a l'inici de la pàginaCorrecció i revisió
És molt important destacar que tot el procés de recerca ha estat fet manualment d’una manera acurada. La verificació en l’obtenció correcta dels diferents arxius i fitxers ha estat feta en els diferents passos, cosa que ha permès una millor anàlisi posterior i una correcció de la metodologia i de la recerca en cas necessari.
Si no s’han obtingut els resultats esperats com obtenir una seqüència homòloga, on la proteïna trobada comencés per metionina (M) o no s’ha trobat una regió de la query amb el nostre genoma, s’han dut a terme diferents correccions com:
- Allargar la seqüència si eren pocs aminoàcids
- Valorar si la regió seleccionada podria ser limitant a l’hora d’anotar la proteïna
- Fer un alineament de nou de la regió de la query no trobada amb BLAST de tot el genoma per veure si una part de la proteïna es trobava en un Scaffold diferent
- Executar un tblastn des del NCBI contra els genomes de tots els vertebrats
Cerca dels elements SECIS
Per a què una selenoproteïna pugui incorporar una selenocisteïna en lloc d'un codó stop, cal la presència dels elements SECIS, ubicats a poques bases de distància de la proteïna i en l'extrem 3' de la cadena. Per això, una altra eina útil a l'hora de valorar la presència o no de selenoproteïnes són aquests elements.
Una pàgina web on podem buscar aquests elements és la pàgina SECISearch. Introduïm la regió del genoma proporcionat pel fastasubseq on es troba el hit i el programa ens prediu la presència o no d'elements SECIS.
Si la regió és més gran de 50000 nucleòtids el programa ens donarà error. És per això, que nosaltres hem seleccionat a partir de la proteïna predita la regió on ens interessava trobar-lo, mitjançant el programa fastasubseq. El que hem fet ha sigut seleccionar 15000 nucleòtids des del final de la nostra proteïna, utilitzant la informació del programa Exonerate. D'aquesta manera, seleccionem de forma més acurada la regió on potencialment podem trobar un element SECIS.
$ fasatsubseq regio_query.fa start length > SECIS.query.fa
Aquest programa també ens dóna una puntuació o Score en funció de la fiabilitat per a cada element SECIS trobat. En principi, el llindar a partir del qual ens podríem fiar de què una seqüència realment sigui un SECIS és una puntuació de 15.
Tornar a l'inici de la pàgina