Estructura secundaria de los RNA

                Predicción de estructuras secundarias de los RNA

Por Martí Nolla y Ferran Pons


Material y Métodos

Resultados

Discusión

Referencias

Agradecimientos



Material y Métodos

  1. IREs humanos
  2. tBLASTn IREs humanos-cDNAs ratón
  3. RNAfold



  1. Basándonos en la base de datos NCBI hemos encontrado un total de 61 secuencias proteicas humanas que contienen IREs. La búsqueda la realizamos mediante las palabras clave "IRE", "transferrin" y "ferritin". Estas dos últimas proteínas sabíamos previamente que contenían IREs. Las 61 secuencias están recopiladas en el fichero secuencias


  2. El siguiente paso consistió en hacer un BLAST tirando las 61 proteínas encontradas que contenían IREs contra los 60.000 cDNAs de ratón. El programa utilizado fue el tBLASTn, el cual nos permite comparar secuencias proteicas con secuencias nucleotídicas. El objetivo de este procedimiento fue encontrar proteínas de ratón que tuvieran el mismo patrón que las proteínas con IREs humanas.

    Partiendo de la base que las proteínas de ratón que contuvieran IREs tenían que tener una gran similaridad de secuencia intentamos hacer varios BLASTs con E values distintos. Finalmente nos quedamos con un E value = 0.00001 por la razón expuesta, puesto que de esta manera conseguíamos reducir el número de apareamientos al azar. 

    blastall -p tblastn -i ires.fa -d db/M.musculus/cDNA/fantom2.00.seq.ri -e 0.00001 > prova3.tblastn 
    blastall -p tblastn -i ires.fa -d db/M.musculus/cDNA/fantom2.00.seq.ri -e 0.001 > prova23.tblastn

    Para simplificar los resultados hemos agrupado los querys (IREs humanos) sólo con el número del identificador, y los hemos relacionado con el identificador del subject (cDNAs). 

    grep "Query=|>ri_" prova3.tblastn
    Ejemplo 1 de query y sus subjects relacionados.

    A continuación hemos ordenado los subjects y contado el número total que nos ha dado el BLAST. Obteniendo 424 cDNAs. 

    egrep ">ri_" prova3.tblastn | sort | uniq | wc -l (424)
    Una vez llegados a este punto, lo que hemos hecho es ver cuántas pautas de lectura encontrábamos para todos los cDNAs del resultado del BLAST.


    check frame: 

    grep "Frame" prova3.tblastn | wc -l (el resultado fue de 1045)

    comprobar cuantos están en forward 

    grep "Frame" prova3.tblastn | grep "+" | wc -l (el resultado fue de 1016)

    comprobar cuantos están en reverse 

    grep "Frame" prova3.tblastn | grep "-" | wc -l (el resultado fue de 29)


    Encontrando un total de 29 cDNAs que estaban en reverse hemos procedido a identificarlos 

    gawk '{if (substr($1,1,1)==">") printf "\n%s ",$1; if ($1=="Frame") printf "%s ", $3}' prova3.tblastn | grep "-" | more

    Hemos elegido uno de estos cDNAs en reverse para estudiarlo un poco más, ya que a priori no esperaríamos encontrar secuencias IREs en reverse, pero por sorpresa nuestra éste presentaba muy buen score. 

    cDNA in reverse:

>ri_A230097D22_PX00130O17_1391
          Length = 1391

 Score =  100 bits (249), Expect = 1e-20
 Identities = 46/46 (100%), Positives = 46/46 (100%)
 Frame = -3

Query: 142  QFWQYGEWVEVVVDDRLPTKDGELLFVHSAEGSEFWSALLEKAYAK 187
            QFWQYGEWVEVVVDDRLPTKDGELLFVHSAEGSEFWSALLEKAYAK
Sbjct: 1116 QFWQYGEWVEVVVDDRLPTKDGELLFVHSAEGSEFWSALLEKAYAK 979



 Score = 71.6 bits (174), Expect = 6e-12
 Identities = 33/34 (97%), Positives = 33/34 (97%)
 Frame = -3

Query: 406 TFLVGLIQKHRRRQRKMGEDMHTIGFGIYEVPEE 439
           TFLVGLIQKHRRRQRKMGEDMHTIGFGIYEVP E
Sbjct: 114 TFLVGLIQKHRRRQRKMGEDMHTIGFGIYEVPAE 13


  3. Utilizando el programa RNAfold calculamos la energía libre de Gibbs para una serie de cDNAs que pueden contener IREs. Estos cDNAs fueron obtenidos mediante un programa que selecciona aquéllos que presentan un patrón que presumiblemente está en los IREs, el patrón laxo . Sabemos que los IREs son muy estables, y por lo tanto tienen un valor de energía libre de Gibbs negativo, ya que en teoría cuanto más bajo es este valor, más terodinámicamente estable son.

    El fichero laxo contiene aquellos cDNAs en forward filtrados por el patrón laxo. Se han seleccionado 601 cDNAs de los 60.000 iniciales. El resultado del RNAfold lo guardamos en sortida3.RNAfold 

     RNAfold < IRESpatrolaxeforward.fasta > sortida3.RNAfold

    Para facilitar la manipulación de los datos hemos quitado los puntos, paréntesis y la secuencia en FASTA del resultado obtenido por el RNAfold, quedándonos sólo con el identificador del cDNA y el valor de la energía libre de Gibbs en el fichero sortida3.RNAfold.tbl 

     sed -e 's/[)(]//g' -e /^[UGCA]/d sortida3.RNAfold | gawk '{if (substr($1,1,1)==">"); printf "\n%s "$1;if (substr($1,1,1)!=">") {$1
=""; print $0}}' | more

    Con el fichero limpio, lo que realizamos a continuación es un filtrado del resultado dado por el RNAfold. Leyendo los resultados obtenidos en el trabajo sobre IREs 2001-2002 nos damos cuenta que no encontraron ningún IREs con un valor de enegía libre de Gibbs superior a -1.80, por lo cual decidimos realizar el primer filtraje con un cut off de -2.00, obteniendo un total de 530 secuencias sobre el total de 601 de que partíamos. Introducimos las secuencias en sortida5.RNAfold . 

    gawk '{if($2 <= -2) print $0}' sortida3.RNAfold.tbl >sortida5.RNAfold  (530)

    Repetimos este paso pero cogiendo un cut off de -3.00 para evitar los falsos positivos que se encontrarían en el filtraje anterior. Introducimos los resultados en sortida6.RNAfold. Encontramos un total de 428 secuencias. 

    gawk '{if($2 <= -3) print $0}' sortida3.RNAfold.tbl >sortida6.RNAfold  (428)

    Repetición del filtraje para un cut off de -6.00. Colocamos los resultados en sortida7.RNAfold. Obtenemos un total de 166 secuencias. 

    gawk '{if($2 <= -6) print $0}' sortida3.RNAfold.tbl >sortida7.RNAfold  (166)


    Para elegir el cut off adecuado nos hemos basado en el porcentaje de IREs reales totales contenidos en cada uno de los cut off. Estos valores fueron obtenidos del trabajo sobre IREs 2001-2002 

     deltaG<=-2 (90.6%)

 deltaG<=-3 (53.12%)

 deltaG<=-4 (46.88%)

 deltaG<=-5 (43.75%)

 deltaG<=-6 (25%)

 deltaG<=-7 (15.63%)

    En un principio esperamos coger el cut off menos negativo, puesto que de esta manera nos aseguraríamos de no perder ninguna secuencia con IREs. Esta elección tendría el problema que cogeríamos también muchos falsos positivos que tendrían que ser filtrados posteriormente. Para realizar esto es necesario hacer un filtraje mediante otro método, el cual necesita analizar secuencias mínimamente largas para poder campararlas, lo cual no sucede con las secuencias IREs. Por tanto, de esta manera obtendrímos una validación de los resultados poco sensible, recogiendo muchos falsos positivos. Así pues, hemos decidido basarnos en un cut off de -3 presumiendo que vamos a perder algunas de las secuencias IREs pero que, a su vez, vamos a reducir mucho el número de falsos positivos.

    Resultados

    Una vez tenemos el resultado a partir de nuestro tBLASTn y del patrón laxo corrido por el programa RNAfold y filtrado procedemos a comparar las secuencias IREs obtenidas por uno y otro método. Esto lo podemos realizar con los distintos cut off que hemos utilizado, de forma que comparamos de mayor o menor número de secuencias.

    Para poder hacer la comparación anterior correctamente de las secuencias tenemos que limpiar nuestros ficheros mediante la orden: 

    sort sortida5.RNAfold | gawk -F"[ :]" '{print $1,$4}' > sortida5.RNAfold.sort
sort sortida6.RNAfold | gawk -F"[ :]" '{print $1,$4}' > sortida6.RNAfold.sort  
sort sortida7.RNAfold | gawk -F"[ :]" '{print $1,$4}' > sortida7.RNAfold.sort

    Ahora ya podemos comparar los resultados de los dos métodos y obtenemos: Por un cut off de -2 obtenemos 14 secuencias con IREs 

    join prova3.tblastn.cDNAs sortida7.RNAfold.sort | more 
ri_0610006G13_R000001C10_918 -2.50
ri_0610006K08_R000001M11_922 -2.50
ri_0610011I07_R000002J21_403 -2.50
ri_1110001P10_R000013E13_922 -2.50
ri_1110070P14_ZA00011E19_866 -3.20
ri_1300002E02_R000010J05_2911 -6.80
ri_1300011K05_R000011G10_2966 -6.80
ri_2510027K02_ZX00047P10_923 -2.50
ri_2510030C10_ZX00082I24_921 -2.50
ri_2600017I12_ZX00048A04_924 -2.50
ri_E430019P14_PX00099O06_3004 -6.50
ri_E430022I20_PX00099B18_927 -2.50
ri_E430032I15_PX00101O22_2071 -3.20
ri_F630105G17_PL00015G24_1979 -3.20

    Por un cut off de -3 obtenemos 6 secuencias con IREs 

    join prova3.tblastn.cDNAs sortida6.RNAfold.sort | more
ri_1110070P14_ZA00011E19_866 -3.20
ri_1300002E02_R000010J05_2911 -6.80
ri_1300011K05_R000011G10_2966 -6.80
ri_E430019P14_PX00099O06_3004 -6.50
ri_E430032I15_PX00101O22_2071 -3.20
ri_F630105G17_PL00015G24_1979 -3.20

    Por un cut off de -6 obtenemos 3 secuencias con IREs 

    join prova3.tblastn.cDNAs sortida7.RNAfold.sort | more
ri_1300002E02_R000010J05_2911 -6.80
ri_1300011K05_R000011G10_2966 -6.80
ri_E430019P14_PX00099O06_3004 -6.50

    El siguiente paso a realizar, es una validación comparativa human-mouse, de los datos del RNAfold una vez filtrados, y un blast de los 601 IREs contra el genoma humano. Pero por problemas técnicos decidimos hacer el blast contra ESTs para mejorar la especificidad, pero los problemas continuan debido a la corta longitud de los IREs. validacion

    Por esta razón realizamos una validación termodinámica de aquellas 51 posibles proteínas con un IRE en el mRNA, que contengan IREs lo suficientemente estables. Una vez lo hemos hecho por las distintos valores de delta de Gibbs, nos quedamos con el deltaG<=-2. 

    join sortida3.RNAfold.tbl_2.id IRESfantom.geneid_ORF_IRE.id | wc -l (45 cDNAs con deltaG<=-2)  
join sortida3.RNAfold.tbl_3.id IRESfantom.geneid_ORF_IRE.id | wc -l (36 cDNAs con deltaG<=-3)
join sortida3.RNAfold.tbl_4.id IRESfantom.geneid_ORF_IRE.id | wc -l (30 cDNAs con deltaG<=-4)      
join sortida3.RNAfold.tbl_5.id IRESfantom.geneid_ORF_IRE.id | wc -l (25 cDNAs con deltaG<=-5)
join sortida3.RNAfold.tbl_6.id IRESfantom.geneid_ORF_IRE.id | wc -l (21 cDNAs con deltaG<=-6)
join sortida3.RNAfold.tbl_7.id IRESfantom.geneid_ORF_IRE.id | wc -l (16 cDNAs con deltaG<=-7)

    Como se puede comprobar el número total de cDNAs que contienen IREs no es demaiado elevado para ninguno de los delta de Gibbs, de manera que decidimos coger el valor de -2, ya que de esta manera nos aseguramos una sensibilidad del 90%, es decir, sabemos con casi total seguridad que todos los IREs predichos lo son realmente y a la vez nos quedamos con muy pocas secuencias (45) respecto del total que partíamos (51). Por tanto, en un principio, tambiéen alcanzamos una especificidad elevada, ya que con mucha sensibilidad obtenemos muy pocos candidatos menos.

    Ahora necesitamos saber la descripción de estos cDNAs en un base de datos proteica, de manera que, una vez transformados estos cDNAs a la proteína predicha por el programa geneid, realizamos un BLASTP. 

    sed -n '/>/,/^$/p' IRESfantom.geneid | gawk -F"|" '{print $1}' | sed 's/_[0-9]$//' | FastaToTbl | gawk 'length($2)>=15' | sort |
join -IRESfantom.geneid_ORF_IRE.id_RNAfold.2 | TblToFasta >IRESfantom.geneid_ORF_IRE.id_RNAfold.2.prot.fa

    blastall -p blastp -i IRESfantom.geneid_ORF_IRE.id_RNAfold.2.prot.fa -d IPI/ipi.HUMAN.fasta -e 0.00001 > IRESfantom.geneid_ORF_
IRE.id_RNAfold.2.prot.ipi.blastp

    Con el BLASTP podemos dividir los resultados obtenidos en tres grupos:




    Discusión



    Una vez analizadas las 7 secuencias expuestas anteriormente, es obvio que la validación termodinámica no es todo lo óptima que desearíamos. Por lo tanto, proponemos seguidamente una serie de opciones que permitirían mejorar la predicción final de estructuras IREs:

    - Como conocemos perfectamente la estructura que adoptan los IREs podríamos restringir las condiciones de análisis del programa Patscan forzando a que éste solamente haga las predicciones junto con su estabilidad termodinámica de aquellas secuencias que puedan adoptar la estructura que nosotros previamente determinamos. Es decir, apareamiento, bucle con una C desapareada, apareamiento y bucle final.

    - Una vez obtenidos los resultados de nuestra validación termodinámica, podrímos lanzar cada una de nuestras 7 secuencias que presumiblemente contienen IREs desconocidos contra una base de datos de dominios proteicos, como por ejemplo Interpro. De esta manera, si alguna de estas secuencias tuviera un dominio de unión a hierro podríamos confirmarlo con la base de datos y nos permitiría reafirmar la validación termodinámica predicha anteriormente.

    Referencias



    Agradecimientos


    Queremos agradecer:

    por su paciencia y soporte técnico a Castellano, S.

    por sus magníficas prácticas a Castelo, R.

    por sus clases teóricas a Guigó, R.