pseudoalineamientos y conteos de tránscritos

Hola,
durante el análisis de unos experimentos de RNAseq mi colega Carlos Cantalapiedra y yo nos desesperábamos al calcular los valores de expresión de cada transcrito en una serie de condiciones, dado que para cada una de ellas era necesario alinear las lecturas/reads originales (en formato FASTQ) contra los transcritos ensamblados. Al menos éste era el protocolo habitual antes de aparecer el algoritmo genérico de pseudoalineamiento, que se describe en la siguiente figura:

Pseudoalineamiento de reads a partir de su composición en k-meros (k=3 en el ejemplo) y un  vector de sufijos de un transcriptoma. Figura tomada de http://bioinformatics.oxfordjournals.org/content/32/12/i192.full.  

Este tipo de algoritmos, implementados en software como kalllisto o RapMap (éste último con licencia GPL), permiten estimar de manera muy eficiente y precisa con qué transcritos son compatibles las lecturas de un archivo FASTQ, es decir, a qué secuencias de un transcriptoma mapean, sin necesidad de alinear base a base los reads. Los alineamientos se pueden calcular, si fuera preciso, después del mapeo. 

Hasta la fecha se han propuesto al menos dos implementaciones del algoritmo genérico de pseudoalineamiento, que se diferencian fundamentalmente en las estructuras de datos utilizadas, como se explica en estos blogs (1 , 2) [en el segundo se hace un repaso a la evolución de este tipo de algoritmos y a su nomenclatura].

En los próximos párrafos muestro dos aplicaciones usando ejecutables de RapMap y como datos un fichero (tr.fna) con 67K transcritos de Arabidopsis thaliana, ensamblados de novo con trinity, y otro (1.f1) con 89M de lecturas Illumina SE en formato FASTQ.

1. Mapeo de reads (20 cores, fichero de salida SAM)

          operación   tiempo   comando      

BWAmem       índice    1m28s   bwa-0.7.12/bwa index tr.fna    

BWAmem        mapeo    9m57s   bwa-0.7.12/bwa mem -t 20 tr.fna 1.fq > 1.sam 

RapMap       índice      58s   RapMap-0.2.2/bin/rapmap quasiindex -t tr.fna -i ./ 

RapMap        mapeo    5m56s   RapMap-0.2.2/bin/rapmap quasimap -t 20 -i ./ -r 1.fq \

                                 -o 1.sam 

Por comparación con BWA queda claro que incluso generando alineamientos SAM estos algoritmos son mucho más rápidos. 


2. Conteo de transcritos (20 cores, implementación Sailfish)

          operación   tiempo   comando      

Sailfish     índice    1m06s   SailfishBeta-0.10.0/bin/sailfish index -t tr.fna -o qindex \

                                  -p 20   

Sailfish      mapeo    1m09s   SailfishBeta-0.10.0/bin/sailfish quant -i qindex/ -r 1.fq \

                                  -p 20 -o 1.bur-0.sf --auxDir tmpsf --dumpEq --libType SF
  

Como se muestra en la tabla, la estimación de valores de expresión de transcritos por pseudoalineamiento y conteo de reads es muy eficiente. En este ejemplo se genera un archivo de salida (quant.sf) con este contenido:

Name      Length EffectiveLength TPM NumReads
TR1|c0_g1_i1 517 316.623      5.3984 214.782
TR2|c0_g1_i1 390 191.501     3.36608     81
TR3|c0_g1_i1 699 498.611      10.502 658
...

Entre las opciones del program está la posibilidad de calcular los conteos con bootstraping, como hace kallisto, aunque en la versión que yo he probado es experimental.

En la siguiente entrada veremos más aplicaciones,
un saludo,
Bruno

Correlación entre transcriptoma y proteoma

Hola,
una cuestión que ya dio qué hablar en los tiempos de los chips de RNA (microarrays) y que ha vuelto a resurgir recientemente es hasta qué punto podemos hacernos una idea de la actividad celular si sólo medimos la expresión de los genes en forma de mRNAs, cuando en realidad muchas de las funciones las desempeñan proteínas.
Por ejemplo, Tian et al.(2004) publicaron correlaciones del 40% en líneas celulares hematopoyéticas de mamíferos.

Sin embargo, ahora en muchos ámbitos los microarrays están siendo sustituidos por experimentos de RNAseq, y por tanto es pertinente volver a hacerse la pregunta, dado que se espera que esta tecnología sea superior. Dos artículos muy recientes nos ayudan precisamente a valorar esto.

El trabajo de Jiang et al.(2011), que propone el uso de controles para corregir los niveles de expresión medidos por RNAseq, incluye una figura que resume la distribución típica de valores de expresión génica que se obtienen por RNAseq:

Figura original de Jiang et al. (2011). En el A panel se muestra la equivalencia entre 'copias por célula' y 'fragmentos por kilobase mapeada' (FPKM). En el B se muestra el histograma de genes expresados, que se solapa a la izquierda (un hombro) con los valores de genes que se expresan menos de una copia por célula. Los datos se obtuvieron a partir de material extraído de una línea celular de Drosophila melanogaster.

El trabajo de Nagaraj et al. (2011), que se basa en material de la archiconocida línea celular humana HeLa, encuentra una distribución similar de mRNAs, lo que sugiere que debe ser una distribución estándar, pero además la correlaciona con valores de abundancia de péptidos medidos por espectrometría de masas de alta resolución. La figura siguiente resume un montón de datos experimentales:

Figura original de Nagaraj et al. (2011).

En el panel A efectivamente se reproduce el hombro observado en los datos de RNAseq de D.melanogaster y en el B se muestran los datos equivalentes de proteómica. El diagrama de Venn C es interesante, porque muestra que hay muchos mensajeros que no se observan como proteína, y algunas proteínas cuyo mensajero está ausente. El panel D muestra la distribución funcional de los genes expresados y, finalmente, el diagrama de dispersión E muestra la correlación observada entre ambos tipos de mediciones, con un coeficiente de correlación de Spearman de 0.6,
un saludo,
Bruno

Actualizaciones posteriores:
[1] En este trabajo se estudia cómo se heredan los niveles de expresión en mosca
[2] En este artículo se comenta que los genes ortólogos entre distintos organismos tienen mayores correlaciones en sus [proteína] que en [mRNA]