Anotació de genomes (II)

Bioinformàtica - 2016/2017 - UPF

Exercici: trobar l’anotació d’una selenoproteïna humana a Drosophila Melanogaster

Intenteu trobar l’anotació de la selenoproteïna humana:

>SPS2_human
MAEASATGACGEAMAAAEGSSGPAGLTLGRSFSNYRPFEPQALGLSPSWRLTGFSGMKGU
GCKVPQEALLKLLAGLTRPDVRPPLGRGLVGGQEEASQEAGLPAGAGPSPTFPALGIGMD
SCVIPLRHGGLSLVQTTDFFYPLVEDPYMMGRIACANVLSDLYAMGITECDNMLMLLSVS
QSMSEEEREKVTPLMVKGFRDAAEEGGTAVTGGQTVVNPWIIIGGVATVVCQPNEFIMPD
SAVVGDVLVLTKPLGTQVAVNAHQWLDNPERWNKVKMVVSREEVELAYQEAMFNMATLNR
TAAGLMHTFNAHAATDITGFGILGHSQNLAKQQRNEVSFVIHNLPIIAKMAAVSKASGRF
GLLQGTSAETSGGLLICLPREQAARFCSEIKSSKYGEGHQAWIVGIVEKGNRTARIIDKP
RVIEVLPRGATAAVLAPDSSNASSEPSS

en el genoma de Drosophila Melanogaster (mosca). Les seqüències del genoma d’aquest insecte les trobarem en un únic fitxer FASTA en el següent lloc del sistema de fitxers:

/cursos/20428/BI/genomes/D.melanogaster/flybase/dmel-all-chromosome-r6.02.fasta

Seguiu les següents passes:

  1. necessitareu fer anar el software del NCBI Blast, el software exonerate, el software GeneWise, el script en Perl fastaseqfromGFF.pl i el programa d’alineament global de seqüències T-COFFEE, la instal.lació dels quals a les aules del campus del mar requereix que hi introduiu primer les següents ordres cada cop que els aneu a utilitzar en una nova finestra terminal:
module load modulepath/goolf-1.7.20
module load BLAST+/2.2.30-goolf-1.7.20
module load Exonerate/2.2.0-goolf-1.7.20
module load T-Coffee/11.00.8cbe486-goolf-1.7.20
export PATH=/cursos/20428/BI/bin:$PATH
export PATH=/cursos/20428/BI/soft/genewise/x86_64/bin:$PATH
export WISECONFIGDIR=/cursos/20428/BI/soft/genewise/x86_64/wise2.2.0/wisecfg/

tal i com feiem a la pràctica anterior.

  1. feu una base de dades de BLAST del genoma de mosca a partir del fitxer donat amb la comanda:
makeblastdb -in /cursos/20428/BI/genomes/D.melanogaster/flybase/dmel-all-chromosome-r6.02.fasta -dbtype nucl -out dm2.fa
  1. localitzeu la regió genòmica on es troba la selenoproteïna humana mitjançant una cerca TBLASTN. Hi ha algun HSP on poguem observar la selenocisteina alineada amb un altre aminoàcid? Observant aquest aminoàcid, quin HSP és més versemblant?

  2. extreieu la seqüència genòmica de la regió trobada. Aquest cop necessitareu l’ajuda dels següents programes que acompanyen a exonerate, per tal de poder extreure primer la seqüència que conté la regió, i després extreure la regió a partir de la seqüència (fastasubseq només funciona sobre un fitxer FASTA que contingui una única seqüència):

fastaindex /cursos/20428/BI/genomes/D.melanogaster/flybase/dmel-all-chromosome-r6.02.fasta dm2.index
fastafetch /cursos/20428/BI/genomes/D.melanogaster/flybase/dmel-all-chromosome-r6.02.fasta dm2.index "nomseq" > nomseq.fa

on nomseq es refereix al nom de la seqüència en la qual estem interessats. A partir d’aquí ja podem extreure la regió genòmica amb fastasubseq com haviem vist abans.

  1. genereu una anotació amb exonerate, emmagatzemant l’anotació de les coordenades dels exons en un fitxer en format GFF. Si creieu que haurieu de trobar el gen en la regió extreta pero l’anotació resultant no inclou tots els exons que esperaveu, proveu d’utilitzar l’opció --exhaustive yes (si esteu interessats en saber que fa aquesta opció consulteu la següent plana d’ajuda d’aquest software que trobareu aquí.

  2. genereu una anotació amb GeneWise, emmagatzemant l’anotació de les coordenades dels exons en un fitxer en format GFF.

  3. comproveu si tots dos programes ens donen la mateixa anotació comparant els fitxers GFF.

  4. decidiu amb quina anotació us quedeu o si construiu un fitxer GFF on l’anotació final la fessim a partir de les anotacions d’exons que trobem més verosímils d’entre tots dos programes.

  5. fent anar el programa fastaseqfromGFF.pl emmgatzemeu en un fitxer FASTA el cDNA corresponent i, mitjançant el programa fastatranslate, en un altre fitxer FASTA emmagatzemeu la proteïna corresponent.

  6. feu un alineament global de la selenoproteïna humana amb la proteïna predita a partir de la vostra anotació mitjançant el programa T-COFFEE. Recordeu que aquest programa el podeu fer anar tant desde la seva web com en linia de comandes fent:

t_coffee <fitxerFASTAsequencia1> <fitxerFASTAsequencia2>
  1. comproveu si la proteïna que heu obtingut te homòlegs en altres espècies (concretament hauriem de poder recuperar la selenoproteïna humana) mitjançant una cerca BLASTP contra el conjunt no-redundant de totes les proteines disponibles a NCBI utilitzant el software netblast de la següent forma:
blastp -query <fitxerFASTAproteïna> -db /cursos/20428/BI/db/blast/nr.fa

hi ha algun hit a una versió sencera d’aquesta mateixa proteïna en Drosophila Melanogaster? si és així recupereu la seqüència d’aquesta proteïna a traves de GenBank i alinieu-la amb la regió genòmica en la que hem estat treballant: alinea sencera? us permet generar una millor anotació?

  1. aneu a la bases de dades de selenoproteines SelenoDB, busqueu l’anotació de la selenoproteïna Selenophosphate synthetase (SEPHS2) en Drosophila Melanogaster i compareu-la amb la que heu obtingut vosaltres mateixos.

Altres eines i consells per la cerca de selenoproteines

  • altres programes que acompanyen l’exonerate i que us poden ser d’utilitat són el fastatranslate per traduir d’ADN a proteïna i el fastarevcomp per calcular la seqüència complementaria inversa.
  • cerca d’elements SECIS, ho podeu fer a través de la plana del software SECISearch3
  • l’aminoàcid de la selenocisteina es representa amb U, pero alguns programes el poden representar amb un asterisc (*).
  • dels programes que llegeixen seqüències d’aminoacids, alguns no poden treballar amb una U i per tant l’haurem de reemplaçar per un altre símbol com ara, un asterisc (*) o una X. Alguns altres no poden treballar tampoc amb un asterisc.
  • aquelles selenoproteines que continguin les anomenades caixes redox [UC]XX[UC] on [UC] es refereix a o be un aminoàcid selenocisteina o be una cisteina i X es refereix a qualsevol aminoàcid, en general els seus homòlegs haurien de tenir també conservades aquestes caixes redox.

Automatització mitjançant scripts del shell

En un treball o projecte de l’àmbit de la bioinformàtica sovint ens trobarem que hem de repetir moltes vegades una mateixa comanda en el shell però canviant els arguments. Un exemple típic és el de fer una cerca BLAST d’una mateixa proteïna contra tot un seguit de differents bases de dades de BLAST. El fet de que poguem adreçar la qüestió que volem resoldre mitjançant línies de comanda en Unix permet automatizar aquest tipus de tasques, en contraposició a la utilizació d’eines a través de la web. La forma d’automatitzar tasques en Unix consisteix en escriure programes (scripts) en el llenguatge de programació del shell del Unix.

Podem escriure un script del shell en els següents dos passos:

  1. feu un fitxer amb un editor de text (com ara l’emacs) que tingui com a primera linia la següent:
#!/bin/bash

i a continuacio hi escriviu les ordres que voldrieu que executés el shell del Unix.

  1. doneu permissos d’execució (chmod u+x) i crideu l’script tal i com ho feu amb qualsevol altre programa.

A continuació hi trobareu una serie d’exemples que il.lustren la creació de scripts del shell a partir dels quals heu de poder escriure els que pogueu necessitar al llarg d’aquesta assignatura. Tots els exemples fan anar el següent fitxer:

/cursos/20428/BI/genomes/vertebrates/genomes_list.tab

el qual conte informació sobre els genomes de protistes emmagatzemats al sistema de fitxers del campus del mar. Haureu d’examinar primer el seu contingut per tal d’entendre els exemples. Per tal de provar els exemples, els heu de copiar en un fitxer de texte seguint els dos passos esmentats anteriorment per poder-los executar:

Exemple 1:

for any in 2012 2013 ; do {
  echo Els genomes de vertebrats pels treballs de l\'any $any son
  grep $any /cursos/20428/BI/genomes/vertebrates/genomes_list.tab | cut -f 2 ;
} done

Exemple 2:

for genome in Ailuropoda_melanoleuca Cricetulus_griseus ; do {
  blastdb=`grep $genome /cursos/20428/BI/genomes/vertebrates/genomes_list.tab | cut -f3`
  echo la bb.dd. de blast del genoma de $genome esta a $blastdb
} done

Exemple 3:

for genome in Ailuropoda_melanoleuca Cricetulus_griseus ; do {
  blastdb=`grep $genome /cursos/20428/BI/genomes_list.tab | cut -f 3`
  tblastn -query sel15human.aa.fa -db $blastdb -out sel15humanCONTRA$genome.tblastn
} done