Sin embargo, en algunos casos es posible quedarse digamos en el tejido conectivo, una mezcla de información proveniente de documentación, listas de correo, foros, logs del programa, diversos tests y publicaciones que explican ciertos entresijos del programa sin llegar al detalle del código. Así avanzaremos un poco en la compresión de su funcionamiento y, lo que es más importante, de los resultados obtenidos. Sin embargo, puede llevar, como veremos, a nuevas preguntas sin resolver, la identificación de posibles bugs y a ser muy cauto cuando uno tiene que interpretar los resultados de un programa que no conoce en detalle. Sin embargo, lo positivo es que podremos estar al tanto de dichas peculiaridades, e incluso podremos colaborar a mejorar los programas indirectamente informando a los meritorios desarrolladores de los resultados de nuestras pruebas. Podremos también ser más críticos ante los resultados generados con uno de estos programas.
En el caso de esta entrada, vamos a hablar un poco de TopHat2.0.9, que es el programa más utilizado para mapear reads cortos de secuenciación masiva de cDNA a una referencia, normalmente un genoma ensamblado. Es una entrada bastante prescindible para quien no vaya a trabajar con TopHat2 o programas de alineamiento similares, aviso.
Digamos que TopHat2 tiene 3 formas básicas de operación, según vayamos añadiendo una nueva capa de filtrado.
La fase común a toda ejecución de TopHat2 consiste en un mapeo contra el genoma. Todos aquellos reads no mapeados serán fragmentados, y dichos segmentos volverán a ser alineados, con el objetivo de encontrar posibles nuevas uniones entre exones, de forma que finalmente podamos mapear dichos reads en toda su longitud.
Sin embargo, incluyendo la anotación de un transcriptoma, en formato GFF o GTF, TopHat2 mapeará antes de nada nuestros reads sobre las regiones especificadas, de forma que los reads que alineen en el transcriptoma los va a considerar directamente como mapeos buenos y dejará de buscar alineamientos para esos reads ("Only the reads that do not fully map to the transcriptome will then be mapped on the genome", http://tophat.cbcb.umd.edu/manual.shtml). Posteriormente, todo lo que no alineó contra el transcriptoma pasa por la fase descrita en el párrafo anterior: mapeo contra el genoma, realineamiento de segmentos de los reads sin hits, definición de nuevas uniones entre exones, alineamiento final de reads.
Por último, para ciertas aplicaciones, como variant calling, ensamblaje o caracterización estructural del transcriptoma; también según con la naturaleza de los datos y referencias con los que se trabaje; puede convenir filtrar aquellos reads con mapeos múltiples. Especialmente cuando se trabaja con referencias ricas en repeticiones, como es el caso de los genomas de algunas plantas, por ejemplo. Con TopHat2 podemos utilizar la opción -M (o --prefilter-multihits) para realizar un prefiltrado mapeando reads contra el genoma. En esta fase, todos aquellos reads que superen el número de alineamientos indicado por el parámetro -g serán descartados. Así, las fases de alineamiento al transcriptoma y subsiguientes, se realizarán únicamente con los reads con el número máximo de alineamientos que deseemos. Obviamente, si queremos reads que alineen una única vez al genoma deberíamos utilizar la opción "-g 1".
Éste es un pequeño pero importante detalle del funcionamiento de TopHat2 cuando se quiere utilizar la opción -M. (Ojo que vienen curvas a partir de aquí :) ). Hay que ser conscientes de que el parámetro -g (o --max-multihits) se comporta como un threshold en la fase de prefiltrado, de forma que los reads que superen dicho threshold serán descartados ("[...] directs TopHat to first align the reads to the whole genome in order to determine and exclude such multi-mapped reads (according to the value of the -g option)", http://tophat.cbcb.umd.edu/manual.shtml). Sin embargo, el parámetro -g en las fases subsiguientes sirve simplemente para indicar a TopHat2 cuántos alineamientos como máximo queremos que aparezcan en el BAM para cada read. Si hay más alineamientos que el valor de -g, simplemente se elige cuáles incluir en la salida aleatoriamente ("[...] will report the alignments with the best alignment score. If there are more alignments with the same score than this number, TopHat will randomly report only this many alignments", http://tophat.cbcb.umd.edu/manual.shtml).
Esto trae efectos secundarios, quizás inesperados. Sólo podemos definir -g una vez, y será utilizado en todas las fases. En el prefiltrado, utilizando "-g 1", eliminaremos todos los reads con alineamientos múltiples, al genoma. Pero no mapeos múltiples generados en fases posteriores*. Si queremos detectar dichos mapeos múltiples no podremos en principio, porque con la opción "-g 1" se nos informa de un sólo alineamiento por read, y además los valores MAPQ y NH:i:1 del BAM indican todos y cada uno de los reads como de alineamiento simple (!). Esto implicaría que no podríamos diferenciar dichos reads tampoco en programas que utilicen el BAM como entrada (como samtools, por ejemplo).
Hablando de NH:i, esta es otra peculiaridad de TopHat2 que hay que conocer. A pesar de que utiliza Bowtie2 como algoritmo de alineamiento, de alguna forma se ha descartado el uso de la tag XS:i (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#bowtie2-build-opt-fields-xs) y se ha optado por utilizar NH:i. Sin entrar en las tripas del programa no sabemos exactamente en qué momento se realiza dicha sustitución. Ni por supuesto el porqué o para qué. Esto no está recogido en la documentación de TopHat2.
Volviendo al prefiltrado: ¿por qué no obviar el uso de la opción -M y directamente filtrar los mapeos múltiples marcados tras el mapeo al transcriptoma y al genoma? Utilizando la opción del prefiltrado obtendremos muchos más mapeos múltiples, ya que no se da por sentado que el mapeo al transcriptoma es el correcto para un read en particular (como se ha explicado más arriba, al comentar la fase de mapeo al transcriptoma). Véase la comparación de un mapeo con y sin -M:
| contigs mapped | Reads mapped | ||
| tophat_s_04 | 508,193 | 28,660,380 | 87.95% |
| tophat_s_04_M_g2 | 41,321 | 16,375,099 | 50.25% |
La diferencia de las dos filas anteriores se debe únicamente al uso de "-M -g 2" (activos en la segunda fila). Como puede verse, el resultado final es muy diferente. Esta diferencia podría ser mayor incluso con la opción "-g 1". La disminución en el número de reads mapeados es importante, pero cómo se refleja esto en la diversidad de hits sobre la referencia es abismal. En el segundo caso estamos mapeando con una confianza mucho mayor a una región más reducida del genoma, quizás indicando que está compuesta por zonas de alta mapabilidad.
Ahí ya depende del criterio de cada uno, del objetivo del análisis, de la mapabilidad del genoma y la completitud de su anotación. Quizás para genomas con un transcriptoma muy bien anotado, en el que todo lo que quede fuera puedan ser pseudogenes o repeticiones que se sabe no son codificantes y si se están análizando únicamente genes que codifiquen para proteínas, por poner un ejemplo, sea más conveniente no usar el prefiltrado. Quizás para genomas sin acabar, con regiones de transcriptoma no del todo delimitadas sea más seguro utilizar la opción -M, especialmente si se están tratando de detectar loci que puedan ser regiones transcritas no identificadas previamente, como posibles parálogos.
Otra cosa interesante que me ha quedado por resolver respecto a los mapeos múltiples es en referencia a la nueva salida que genera TopHat2.0.9: align_summary.txt. Utilizando "-g > 1" los reads con alineamientos múltiples aparecen con MAPQ < 50 y NH:i>1 en el BAM. Sin embargo, no he conseguido ningún alineamiento en el que el número de reads así indicado coincida con el número de reads al que se refiere en el fichero align_summary.txt.
Respecto al pre-filtrado y su compatibilidad con la detección de mapeos múltiples en fases posteriores (debido a lo comentado antes), podríamos de nuevo cuestionarnos ¿Por qué no realizar un mapeo inicial mediante Bowtie2 para detectar dichos mapeos y luego ya utilizar TopHat2 de forma habitual? Esto puede realizarse, sin duda. Pero una de las características más interesantes de TopHat2 es el control que tenemos sobre el número de mismatches y bases en gaps, así como la edit distance final en el alineamiento (como vamos a ir viendo). Estos parámetros además son comunes para todas las fases de alineamiento, lo que nos permite ser consistentes en cierta medida en cuanto a qué consideramos un alineamiento válido para nuestras muestras. Veamos un poco más en detalle qué pasos realiza TopHat2 para los alineamientos, gracias a la información presente en logs/run.log.
Para cada alineamiento, TopHat2 ejecuta Bowtie2 mediante un wrapper (bowtie2_align). En principio, TopHat2 realiza los alineamientos con los reads por separado, sean reads pareados o no ("For paired-end reads, TopHat2 processes the two reads separately through the same mapping stages described above. In the final stage, the independently aligned reads are analyzed together to produce paired alignments, taking into consideration additional factors including fragment length and orientation", http://genomebiology.com/2013/14/4/R36). Para la ejecución de Bowtie2 podemos pasarle ciertos parámetros que serán utilizados siempre por el wrapper, como los presets de sensibilidad (que incluyen -D, -R, -N, -L, -i; todos parámetros relacionados con las seeds del alineamiento de Bowtie2) o los de scoring (mismatches, gaps, etc.). También podemos incluir --b2-score-min, que hará las veces de filtro de alineamientos que no superen dicho score.
Es importante diferenciar el uso que hace el wrapper de -k en las distintas fases de alineamiento, porque éste es un parámetro clave de la configuración de Bowtie2 que no podemos cambiar en TopHat2, al menos no de forma directa. Parece, junto a la actividad de fix-map-ordering (que vemos un poco más abajo), uno de los pilares fundamentales del programa, de lo que lo hace un mapeador especializado en datos de RNAseq.
Para la fase de alineamiento al transcriptoma se utiliza -k 60, mientras que para el alineamiento posterior al genoma se utiliza -k 20. Para los segmentos se pasa a -k 41, al igual que para el mapeo de los reads a las nuevas uniones de exones. Curiosamente, en el paso de prefiltrado (-M) se utiliza -k 2. ¿Qué demonios es -k? Pues el -k de Bowtie2 ("[...] it searches for at most
En principio tengo curiosidad por ver qué pasaría si al prefiltrado le indicáramos "-g 3". Pero no lo he probado, ni tiene pinta de que lo haré. Por otro lado, digamos que la búsqueda más amplia que va a realizar TopHat2 para encontrar el mejor alineamiento posible va a ser sobre el transcriptoma, aprovechando además que es un conjunto de secuencias más reducido que el genoma. Unos cuantos alineamientos menos realizará TopHat2 con los segmentos y los realimientos a uniones de exones. Respecto al alineamiento al genoma, si vamos a jugar con un genoma altamente repetitivo, qué garantías nos da -k 20? ¿Habría que modificarlo al trabajar con un genoma de bacteria, o con el genoma del pino? Preguntas sin resolver ¿tú qué opinas?
Siguiendo con los parámetros que le pasamos al wrapper, es importante no ser demasiado estricto en cuanto al score mínimo permitido (--b2-score-min), o incluso el scoring de los mismatches y gaps, si se plantea uno el alinear reads con 3 mismatches (http://genomebiology.com/2013/14/4/R36, Table S2).
El siguiente paso que TopHat2 realiza sobre cada alineamiento es ejecutar el script fix_map_ordering, al que le pasa los parámetros --read-mismatches, --read-gap-length, --read-edit-dist, --read-realign-edit-dist, --bowtie2-min-score. Éste último en principio no lo podemos cambiar. Con lo cual parece un threshold fijo, aunque no me ha quedado claro aún que relación tiene con el score-min de Bowtie2 que podemos cambiar para el wrapper bowtie2_align. Este --bowtie2-min-score cambia en las distintas fases, siendo de 15 para el transcriptoma y el genoma, de 10 para los segmentos y realineamientos a junctions. Su valor es 35 para el alineamiento al genoma de la opción -M. El resto de parámetros sí los podemos cambiar, y esto puede ser realmente útil si queremos controlar el número de mismatches y bases en gaps que queremos aceptar. Cuando cambiemos estos parámetros, serán utilizados en todas las fases de alineamiento, eso sí.
Reads mapped
tophat_s_a 12,806,680 39.30% --read-mismatches 0 --read-gap-length 0 --read-edit-dist 0
tophat_s_b 12,789,305 39.25% --read-mismatches 0 --read-gap-length 0 --read-edit-dist 6
tophat_s_c 15,038,778 46.15% --read-mismatches 1 --read-gap-length 0 --read-edit-dist 6
tophat_s_d 15,056,015 46.20% --read-mismatches 1 --read-gap-length 0 --read-edit-dist 1
tophat_s_e 15,539,057 47.69% --read-mismatches 1 --read-gap-length 3 --read-edit-dist 6
tophat_s_f 16,388,235 50.29% --read-mismatches 2 --read-gap-length 3 --read-edit-dist 6
tophat_s_g 16,478,122 50.57% --read-mismatches 2 --read-gap-length 6 --read-edit-dist 6
tophat_s_h 17,038,705 52.29% --read-mismatches 4 --read-gap-length 6 --read-edit-dist 6
tophat_s_ia 17,046,981 52.31% --read-mismatches 4 --read-gap-length 6 --read-edit-dist 12
tophat_s_ib 17,053,281 52.33% --read-mismatches 4 --read-gap-length 12 --read-edit-dist 12
tophat_s_l 17,188,433 52.75% --read-mismatches 8 --read-gap-length 12 --read-edit-dist 18
--b2-very-sensitive --b2-score-min C,-28,0 -M -g 1 (Parámetros usados en todos los casos)
Esta es una tabla con alineamientos cambiando los parámetros de fix-map-ordering. En éste caso el cambio en el número de mismatches aceptados tiene un impacto mayor en el resultado que el aumento del número de bases en gaps. Especialmente aumenta el número de reads mapeados en el cambio a 1 mismatch, 2 mismatches, 4 mismatches... Ya no aumenta tanto al utilizar 8 mismatches, con 12 bases en gaps, y apenas unos 20,000 alineamientos de los más de 100,000 reads mapeados nuevos, eran de mapeo único (datos no mostrados).
Para el caso de tophat_s_h, se generaron alineamientos modificando --b2-score-min. Como se observa, el número de reads mapeados contra el genoma parece saturarse ya con score entre 8 y 16.
mapped
0 13,848,299 42.50%
4 14,205,101 43.59%
8 16,350,857 50.18%
16 17,118,163 52.53%
32 17,486,162 53.66%
64 17,351,895 53.25%
¿Es esta una saturación real, debido a que los reads mapeables independientemente son cerca de los 17 M en esta muestra? O se debe quizás a que está actuando el score mínimo de fix_map_ordering (--bowtie2-min-score 15). En la tabla anterior vemos los alineamientos del prefiltrado al genoma. Sin embargo, a medida que se reduce el número de reads sin mapear, aumenta el número de reads con mapeos múltiples (datos no mostrados).
De esta forma, el número reportado con TopHat2 es bastante robusto frente a --b2-score-min usando la opción -M (con el genoma de éste ejemplo, etcétera.), como se ve en la tabla siguiente (--b2-score-min de 0 a 64).
De esta forma, el número reportado con TopHat2 es bastante robusto frente a --b2-score-min usando la opción -M (con el genoma de éste ejemplo, etcétera.), como se ve en la tabla siguiente (--b2-score-min de 0 a 64).
| 0 | 16,796,134 | 51.54% |
| 4 | 16,756,691 | 51.42% |
| 8 | 16,573,331 | 50.86% |
| 16 | 16,649,626 | 51.09% |
| 32 | 17,091,930 | 52.45% |
| 64 | 17,125,504 | 52.55% |
Otras opciones importantes, en éste caso para el uso con reads pareados, son --no-discordant y --no-mixed (http://tophat.cbcb.umd.edu/manual.shtml). Curiosamente, estas opciones no aparecen en logs/run.log, por lo que no queda claro cómo y en qué momento las usa TopHat2. En principio parecería que lo lógico es que sean las opciones de Bowtie2 (del mismo nombre), pero por otro lado esto parece poco probable si realmente TopHat2 mapea cada pareja por su lado y finalmente acaba realizando un chequeo de las parejas (momento en el cual podría ser que se apliquen estos parámetros). En cualquier caso, su reflejo en los resultados es importante:
contigs Reads unmapped filtered multiple Reads mapped
M_1_b 54,537 17,447,158 19.62% 36,442,332 40.99% 35,017,524 39.39%
M_1_c 31,550 24,460,620 27.51% 39,179,540 44.07% 25,266,854 28.42% --no-discordant --no-mixed
M_1_b 54,537 17,447,158 19.62% 36,442,332 40.99% 35,017,524 39.39%
M_1_c 31,550 24,460,620 27.51% 39,179,540 44.07% 25,266,854 28.42% --no-discordant --no-mixed
Hasta aquí un breve vistazo a partes del tejido conjuntivo de TopHat2.0.9. Hemos visto que podemos utilizar básicamente 3 modos de operación, y que su uso tiene un reflejo importante en los resultados obtenidos. También hemos resuelto algunas dudas respecto al funcionamiento de parámetros como -g y -M, o al uso de las tags NH:i con "-g 1" (probablemente erróneo). Hemos averiguado también que hay ciertos parámetros que no podemos modificar directamente (como -k o --bowtie2-min-score), y parecen ser características intrínsecas de la naturaleza de TopHat2 en unos casos, auténticas incógnitas en otros. Finalmente, se han generado algunas dudas sin resolver en la interpretación de los resultados, como por qué aumenta el número de reads con alineamiento múltiple con la opción --no-mixed, o qué demonios intenta decirnos el fichero align_summary.txt acerca de nuestros reads y su singularidad.
En fin, si alguien se anima a meterse en las tripas y hacer una segunda parte a éste post, sería fantástico. Por el momento, a ver si en la próxima entrada hago algo un poco más light y os cuento un poco sobre www.circos.ca.
Hasta la vista!
* No está muy claro de dónde proceden dichos mapeos múltiples, pero aparecen. Quizás sea el resultado de varios efectos. Por ejemplo, el mapeo a las nuevas junctions generadas durante la fase de re-alineamiento de segmentos; o que no todo lo que mapeaba de forma múltiple fue encontrado en la fase de prefiltrado. La significación cuantitativa tampoco queda del todo clara. Un ejemplo de los que tengo, de ~32,5 M de reads se prefiltraron ~ 9.7 M. Además, ~ 5.8 M no alinearon al genoma en esta primera fase por lo que fueron descartados también. El resto de reads se volvieron a pasar por TopHat2, ahora sin -M y con -g 2. De unos 17 M de reads, 0,8 M (~ 4,7% en el peor de los casos) tenían mapeos múltiples según align_summary.txt; unos 31 K (~ 0,18% en el mejor) según NH:i y MAPQ.
