1. Búsqueda de querys para anotar los genomas de protistas proporcionados.
Accedimos a la base de datos de selenoproteínas SelenoDB para adquirir la secuencia proteica de las proteínas SelR y MrsA que se habían identificado hasta ahora en diferentes organismos, para poder utilizarlas como query.
En primer lugar, encontramos la secuencia de SelR en el genoma de H. sapiens, S. cerevisiae, D. melanogaster y C. elegans pero la única que era selenoproteína era una de las tres halladas en humano (SelR1) ya que en los demás organismos eran homólogas en cisteína. Como la única proteína con selenocisteína que teníamos, SelR1, era de humano y estaba demasiado alejada filogenéticamente a los protistas, al empezar la anotación del genoma la secuencia no era reconocida por los programas utilizados, y decidimos buscar en la base de datos TARBALL otra query más cercana a los protistas. No obstante, en TARBALL no encontramos la secuencia proteica de SelR, de modo que optamos por buscar secuencias de SelR en los proyectos de los años anteriores ya que los estudiantes la habían identificado en los genomas de algunos protistas.
Así pues, vimos que identificaron diversos homólogs con Cys de SelR en varios genomas, pero sólo la secuencia predicha en el genoma de Dictyostelium purpureum contemplaba el inicio de la secuencia de la proteína, porque empezaba por metionina.
Utilizamos esta secuencia como query para buscar en los genomas propuestos la presencia o ausencia de la familia SelR, pero para asegurarnos que al no utilizar el resto de querys predichas en los otros genomas no estábamos obviando posibles selenoproteínas en los genomas problema, decidimos hacer una predicción en un mismo genoma a partir de dos querys diferentes, comprobando de esta manera que predecíamos lo mismo. Concretamente, utilizamos el homólogo identificado en D.purpureum y otro homólogo identificado en Aureococcus anophagefferens como querys, y predijimos en el genoma de Trypanosoma cruzi la presencia o ausencia de SelR. El resultado obtenido resultó ser el mismo, es decir, predijimos la misma selenoproteína a partir de las dos querys, y además la predicción mediante D.purpureum fue más completa. El resultado del alineamiento global realizado puede ser visualizado haciendo clic aquí.
La secuencia de MrsA la encontramos en la base de datos TARBALL, en el genoma de Chlamydomonas reinhardtii. Dado que se trataba de un protista y empezaba por metionina la utilizamos como query.
2. Alineamiento de secuencias: TBLASTN
Con el programa BLAST (Basic Local Alignment Search Tool) buscamos secuencias homólogas a nuestra query en los genomas de los organismos que nos habían proporcionado. El BLAST es un algoritmo para comparar secuencias biológicas de aminoácidos de varias proteínas, o secuencias de DNA. En concreto el TBLASTN hace alineamientos locales a partir de una query proteica y de una secuencia nucleotídica de un genoma que traduce a proteína.
Conseguimos automatizar el proceso e hicimos el tblastn para todos los genomas del 2009 y del 2010 comparándolos con nuestra query de D. purpureum y C. reinhardtii. Después hicimos una primera selección de los posibles hits basándonos en su e-value. El e-value es el número esperado de alineamientos que podemos obtener con un score igual o superior por azar en un alineamiento múltiple. Cómo mayor sea el e-value, menos significativo será el alineamiento. Para considerar un apareamiento como hit, pusimos como criterio de selección que tuvieran un e-value inferior a e-04. Como segunda condición, en el caso de MrsA miramos que en la posición del genoma problema correspondiente a la selenocisteína de la query hubiera o una U, un codón stop (posible predicción de selenoproteína) o una C (posible homólogo con cisteína). En el caso de SelR, miramos que en la posición del genoma problema correspondiente a la cisteína homóloga a la selenocisteína de humano (SelR1), hubiera una U, un codón stop (posible predicción de selenoproteína) o una C (posible homólogo con cisteína).
3. Obtención de la secuencia genómica de interés.
Una vez hecho el TBLASTN, tenemos que extraer la secuencia genómica de la región encontrada donde se encuentra el alineamiento o alineamientos en un fichero FASTA mediante el programa fastasubseq. Fastasubseq solo funciona sobre este tipo de fichero cuando contiene una única secuencia, por tanto para poder ejecutarlo es necesario primero extraer la secuencia que contiene la región genómica con fastafetch, utilizando las siguientes comandas:
$ fastaindex /cursos/BI/genomes/protists/2010/G.intestinalis/genome.fa giardia.genome.index
$ fastafetch /cursos/BI/genomes/protists/2010/G.intestinalis/genome.fa giardia.genome.index subjectID > fastafetch.fa
Una vez tenemos el fichero en formato FASTA, utilizamos la comanda del fastasubseq en la cual tenemos que establecer los límites upstream y dowstream de la región genómica que queremos extraer. Para saber cuántos nucleótidos de más teníamos que extraer para asegurarnos de coger todo el gen, buscamos en la literatura el tamaño medio de un gen protista, y era de unas 1000 bases. Por tanto, decidimos coger 2000 bases de más por cada extremo.
$ fastasubseq fastafetch.fa start length > genomic.fa
En el caso de que el transcrito del gen fuera reverso, es decir, de 3' a 5', después de hacer el fastasubseq se tiene que ejecutar la comanda "fastarevcomp" de la siguiente manera:
$ fastarevcomp genomic.fa > fastarevcomp.fa
En la mayoría de los genomas, las secuencias que obteníamos del TBLASTN eran contigs cortos y decidimos omitir este paso y correr directamente el Exonerate o el Genewise con todos los contigs donde se encontraban los alineamientos. En los casos que teníamos varios scores dentro de un mismo contig, tampoco realizábamos el fastasubseq ya que de esta forma considerábamos el gen entero y el mejor alineamiento dentro de todos los scores posibles. De todas formas, en algunos casos sí que se tuvo que hacer fastasubseq ya que aunque fuesen diferentes scores dentro de un mismo contig las posiciones de estos eran muy distantes y nos aconsejaron analizarlos por separado.
4. Predicción de la estructura exónica: Exonerate o GeneWise.
A partir de la query de D. purpureum o de la C. Reinhardtii y de la región genómica que hemos extraído anteriormente, generaremos la anotación del gen que da lugar a la proteína mediante el programa exonerate. Exonerate es un programa que alinea secuencias prediciendo cuál podría ser la estructura exónica de la secuencia problema. Utilizamos las siguientes comandas:
$ exonerate -m p2g --showtargetgff -q sel15human.aa.fa -t genomic.fa
$ exonerate -m p2g --showtargetgff -q proteinaquery.fa -t genomic.fa | egrep -w exon > exonerateoutput.gff
La versión actual de este programa tiene en cuenta las diferentes pautas de lectura y además proporciona las secuencias de los intrones, por tanto para extraer el cDNA a partir del fichero GFF obtenido, utilizamos otro programa en Perl denominado fastaseqfromGFF.pl y para ello escribimos la siguiente orden:
$ fastaseqfromGFF.pl genomic.fa exonerateoutput.gff
El Genewise lo utilizamos sólo con los primeros hits. Aunque ambos programas hacen lo mismo, nos recomendaron que tuviéramos preferencia por el Exonerate ya que a pesar de tener la misma función usan criterios diferentes, y por la experiencia de los tutores nos contaron que siempre predecían más genes con el exonerate. Por tanto, solo usamos el Genwise en los primeros momentos para comprobar que nos daba el mismo resultado que el exonerate pero luego no lo seguimos utilizando (salvo que con el exonerate no obtuviéramos ningún resultado).
Algunas veces el programa Exonerate nos mostraba --completed exonerate analysis y no aparecía ningún resultado en el shell. Es en estas situaciones cuando usamos el --exhaustive yes con el que ordenábamos al programa realizar una búsqueda más exhaustiva de los exones. Sin embargo, no es bueno realizar este método por defecto ya que a veces el Exonerate encuentra exones incorrectos o crea intrones imposibles.
5.Traducción de la proteína.
Una vez habíamos obtenido la secuencia de cDNA, lo traducíamos mediante la orden:
$ fastatranslate -F 1 cDNA.fa > proteína.fa
La orden -F 1 se pone porque los tres primeros nucleótidos constituyen un codón, de manera que ya se puede empezar la traducción a partir del primer nucleótido.
6.Alineamiento de las proteínas resultantes: T-Coffee.
Posteriormente hicimos el alineamiento global de la proteína predicha con la query original (D. Purpureum o C. reinhardtii) mediante el programa TCOFFEE que se ordena de la siguiente forma:
$ t_coffee ficheroFASTAsecuencia1 ficheroFASTAsecuencia2
En esta web os mostraremos el T-Coffee conjunto de todos las proteínas predichas en el apartado Resultados.
Cuando al hacer el alineamiento global la U o la C homóloga a la selenocisteína no se alineaba en el genoma problema, hacíamos un híbrido y lo utilizábamos como query para facilitar que el exonerate extendiera la predicción al encontrar un trozo de secuencia idéntica (cogemos el trozo que se ha predicho en el genoma estudiado y el resto de secuencia de la query utilizada).
Otro método que realizamos para extender la predicción fue buscar en el NCBI qué proteína contenía la secuencia predicha en el genoma problema a partir de la query. Si esta parte de la proteína correspondía con SelR, MrsA u otro miembro de la misma familia, la tomábamos como query; mientras que si correspondía a otra familia diferente descartábamos esta posibilidad.
7. Búsqueda de elementos SECIS en posibles selenoproteínas.
Una vez obtenidas las posibles selenoproteínas utilizamos el programa SECISearch para comprobar si contenían el elemento SECIS en 3'. SECISearch localiza posibles dominios de elementos SECIS y valora su estructura tridimensional. Usamos criterios poco restrictivos en la búsqueda para poder identificar dichos elementos.
8. Identificación del codón TGA enfretado con la U o C homóloga.
El programa ExPASy traduce la secuencia de cDNA a secuencia aminoacídica en los 6 frames posibles. Lo utilizamos cuando teníamos una C frente a un codón stop y queríamos averiguar si se trataba de un TGA ya que el codón de la selenocisteína es el TGA.